Mountains

# Config ------------------------------------------------------------------

library(tidyverse)
library(elevatr)
library(sf)
# remotes::install_github("clauswilke/ggisoband")
library(ggisoband)
library(terra)
library(tidyterra)
library(glue)
library(janitor)
library(memoise)


# Parameters --------------------------------------------------------------

# Opentopo API
# https://portal.opentopography.org/myopentopo
# elevatr::set_opentopo_key("xxx")


# Data --------------------------------------------------------------------

# Locations (decimal degrees, WGS84)
locations <- tribble(
  ~name,                     ~lon,       ~lat, ~alti, ~radius, ~iso_interval, ~source,
  "Aconcagua",         -70.017650, -32.658564,  6961,   10000,           150, "aws",
  "Эльбрус",            42.437876,  43.352404,  5643,   10000,           150, "aws",
  "Macizo Vinson",     -85.617388, -78.525399,  4892,   20000,           150, "aws",
  "Puncak Jaya",       137.158613,  -4.078606,  4884,   10000,           150, "aws",
  "Mount Kosciuszko",  148.263510, -36.455832,  2228,   15000,           100, "aws",
  "K2",                 76.513739,  35.881869,  8611,    5000,           180, "alos")
 

# Area of interest --------------------------------------------------------

build_target <- function(lon, lat, radius) {
    tibble(x = lon, 
           y = lat) |> 
    st_as_sf(coords = c("x", "y"), 
             crs = "EPSG:4326") |> 
    st_buffer(radius) |>
    st_convex_hull() 
}

#' Get elevation data
#' 
#' This function is memoised, so that we won't request the data if it has 
#' already been downloaded
#'
#' @param target (sf) : area of interest
#' @param zoom (int) : zoom level ; z = 10 is good, z = 14 max resolution
#' @param source (char) : 
#'
#' @return (df) : with x, y, z, coords
get_elevation <- memoise(
  function(target, zoom = 10, source = "aws") {
    get_elev_raster(target, z = zoom, src = source) |> 
      rast() |> 
      fortify() |> 
      rename(z = 3)
})

#' convert decimal degrees to DMS
#'
#' @param coord (vector of num, length=2)
#'
#' @return (char) : DD°MM'SS.SSH
#'
#' @examples dec_to_sex(c(1.23, 2.55))
dec_to_sex <- function(coord) {
  conv_dec_to_sex <- function(deg_dec, type) {
    deg <- abs(trunc(deg_dec))
    min_dec <- (abs(deg_dec) - deg) * 60
    min <- trunc(min_dec)
    sec <- round(((min_dec - min) * 60), 2)
    h <- case_when(deg_dec < 0  & type == "lon" ~ "W",
                   deg_dec >= 0 & type == "lon" ~ "E",
                   deg_dec < 0  & type == "lat" ~ "S",
                   deg_dec >= 0 & type == "lat" ~ "N")
    sprintf("%02i°%02i′%.2f″%s", deg, min, sec, h)
  }
  glue("{conv_dec_to_sex(coord[2], type = 'lat')} {conv_dec_to_sex(coord[1], type = 'lon')}")
} 


# Final plot --------------------------------------------------------------

plot_isolines <- function(name, lon, lat, alti, radius, iso_interval, 
                          zoom = 10, source = "aws") {
  target <- build_target(lon, lat, radius)
  target_bbox <- st_bbox(target)
  
  p <- get_elevation(target, zoom, source) |> 
    ggplot() +
    geom_sf(data = target, fill = NA, color = NA) +
    geom_isobands(aes(x, y, z = z, fill = z),
                  binwidth = iso_interval, fill = NA) +
    coord_sf(xlim = target_bbox[c(1, 3)],
             ylim = target_bbox[c(2, 4)]) +
    labs(title = glue("{name} — {alti} m"),
         subtitle = dec_to_sex(c(lon, lat)),
         caption = "https://shirt.iresmi.net/") +
    theme_void() +
    theme(plot.title = element_text(family = "Arial narrow", 
                                    face = "bold", size = 22),
          plot.title.position = "plot",
          plot.subtitle = element_text(family = "Arial narrow",
                                       face = "bold", size = 12),
          plot.caption = element_text(family = "Arial", size = 9))
  
    ggsave(glue("{janitor::make_clean_names(name)}_{alti}_m.png"), p,
           dpi = 300, width = 20, height = 20, units = "cm")
}

locations |> 
  pwalk(plot_isolines, .progress = TRUE)