Mountains

Isohypses and T-shirts
R
spatial
datavisualization
Author

Michaël

Published

2023-09-21

Modified

2024-04-12

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)
library(magick)


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

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


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

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

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

#' Create a {sf} layer from a point and radius
#'
#' @param lon (num) longitude (decimal degrees WGS84)
#' @param lat (num) latitude (decimal degrees WGS84)
#' @param radius (num) radius (meters)
#'
#' @return (sf) polygon
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
#'
#' @param target (sf) : area of interest
#' @param zoom (int) : zoom level ; z = 10 is good, z = 14 max resolution
#' @param source (char) : DEM source see {elevatr}. default "aws". "alos" is 
#'   sometimes more detailed and with less artefacts but need an API key
#'
#' @return (data.frame) : with x, y, z coordinates
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
#' 
#' For displaying.
#' Note that we generally use (lon, lat) in computers, like (abscisse, 
#' ordinates) or (x, y) but coordinates are generally displayed as "lat lon"
#'
#' @param coord (vector of num, length=2) : coordinates in decimal degrees : 
#'   c(longitude, latitude)
#'
#' @return (char) : lat+lon coordinates : DD°MM′SS.SS″[NS] DDD°MM′SS.SS″[EW]
#'
#' @examples dec_to_sex(c(-155.468071, 19.820665))
dec_to_sex <- function(coord) {
  conv_dec_to_sex <- function(deg_dec, type) {
    deg <- abs(trunc(deg_dec))
    if (type == "lat" & deg > 90) {
      stop(glue("Latitude ({deg_dec}) shouldn't be over ±90°.
                   Did you mix latitude and longitude ?"))
    }
    if (type == "lon" & deg > 180) {
      stop(glue("Longitude ({deg_dec}) shouldn't be over ±180°.
                   Check your coordinates."))
    }
    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 --------------------------------------------------------------

#' Create a PNG (and its negative) of isolines around a place
#'
#' @param name (char) : place name
#' @param lon (num) : longitude (decimal degrees WGS84)
#' @param lat (num) : latitude (decimal degrees WGS84)
#' @param alti (num) : altitude (metres)
#' @param radius (num) : radius display range (metres)
#' @param iso_interval : elevation between isolines (metres)
#' @param zoom (int) : zoom level (1-14). default 10
#' @param source (char) : DEM source see {elevatr}. default "aws". "alos" is 
#'   sometimes more detailed and with less artefacts but need an API key
#'
#' @return NULL (write on disk in the working directory)
#' @export
#'
#' @examples plot_isolines("Matterhorn", 7.658573, 45.976440, 4478, 5000, 100)
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)
  filename_image <- glue("{make_clean_names(name)}_{alti}_m")
  
  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("{filename_image}.png"), p,
         dpi = 300, width = 20, height = 20, units = "cm")
  
  # create negative image
  image_read(glue("{filename_image}.png")) |> 
    image_negate() |> 
    image_write(glue("{filename_image}_neg.png"))
}

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