Mountains

[This article was first published on r.iresmi.net, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

A map of the Aconcagua mountain consisting of isohypses

Aconcagua – 6991 m

If you like mountains, R and T-shirts, I got you covered. Here we use {elevatr} to get a DEM around some well known summits and make a map consisting only of isohypses. It produces (sometimes) very nice visualizations.

And from these designs we can make nice T-shirts. You can buy some on https://shirt.iresmi.net/ and if your favorite mountain is not there, ask me and I’ll add it (or make it yourself, of course).

I intended to automate the process but sadly the Spreadshirt API doesn’t have an upload feature (!?). If you know a better shop, let me know…

# 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)
To leave a comment for the author, please follow the link and comment on their blog: r.iresmi.net.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.