Using {geodata} to get elevation raster maps

Exploring Global and Country-Specific Elevation Maps with {geodata}: high-resolution raster elevation data using functions like elevation_30s() and elevation_global(). Simplify geospatial analysis and visualization.

Background Map
{terra}
Raster Map
Author

Aditya Dahiya

Published

December 10, 2024

Package {geodata}

The geodata R package (Hijmans et al. 2024) provides functions to download elevation data for any country, primarily sourced from the Shuttle Radar Topography Mission (SRTM). The elevation_3s function retrieves high-resolution elevation data (approximately 90 meters, i.e., 3 arc seconds) for specified coordinates, while elevation_30s offers coarser resolution data (about 1 kilometre, i.e., 30 arc seconds) for entire countries. For global coverage, elevation_global allows users to obtain elevation data at resolutions ranging from 0.5 to 10 arc-minutes. These datasets are essential for various geospatial analyses, including topographic assessments and environmental modelling.

Code 
# Load spatial and environmental datasets.
library(geodata)

# Handle, analyze, and visualize raster and vector data.
library(terra)

# Tidy data workflows with 'terra' objects.
library(tidyterra)

# Data manipulation, visualization, and wrangling.
library(tidyverse)

Getting World Elevation Map

This code downloads global elevation data at a 10-degree resolution using the geodata package, calculates its memory size in R, and visualizes it using ggplot2. The geom_spatraster function plots the raster data, applying a colour scale with squished limits (0-6000 meters), and customizes the map’s title, subtitle, and legend placement.

Code 
world <- geodata::elevation_global(10, path = tempdir())

# With ggplot2
g <- ggplot() +
  geom_spatraster(data = world) +
  scale_fill_wiki_c(
    limits = c(0, 6000),
    oob = scales::squish
  ) +
  labs(
    title = "World Elevation Map",
    subtitle = "Resolution of 10 degrees",
    fill = "Elevation (metres)"
  ) +
  theme(
    legend.position = "bottom",
    legend.key.width = unit(50, "pt")
  )

ggsave(
  plot = g,
  filename = here::here("geocomputation",
                        "images",
                        "elevation_raster_maps_1.png"),
  height = 1400,
  width = 2000,
  units = "px"
)
Figure 1: A simple elevation map of the world using data from {geodata} and geom_spatraster() from {terra}

Countries in {geodata}

This code retrieves country codes from the geodata package, converts them into a tibble, and displays them in a styled, interactive table using the gt package. Column labels are formatted with snakecase and stringr, missing values are replaced with blank text, and a custom theme with a header title is applied for presentation.

Code 
geodata::country_codes() |> 
  as_tibble() |> 
  gt::gt() |>
  gt::cols_label_with(fn = snakecase::to_title_case) |> 
  gt::cols_label_with(fn = str_to_upper) |>
  gt::sub_missing(missing_text = "") |> 
  gt::opt_interactive() |> 
  gtExtras::gt_theme_538() |> 
  gt::tab_header(
    title = "List of countries available in {geodata}"
  )
Table 1
List of countries available in {geodata}

Country-specific raster elevation maps

This code downloads 30-arc-second resolution elevation data for Switzerland using the geodata package and visualizes it with ggplot2. The elevation is displayed as a raster map with a colour scale capped at 6000 meters, custom titles, subtitles, and a bottom-positioned legend with an adjusted width.

Code 
# Taking a smaller country to save data download time
switzerland_raster <- geodata::elevation_30s(
  country = "CHE", 
  path = tempdir()
  )

# With ggplot2
g <- ggplot() +
  geom_spatraster(data = switzerland_raster) +
  scale_fill_wiki_c(
    limits = c(0, 6000),
    oob = scales::squish
  ) +
  labs(
    title = "Elevation Map of Switzerland",
    subtitle = "Resolution of approx. 1 km (30 arc seconds)",
    fill = "Elevation (metres)"
  ) +
  theme(
    legend.position = "bottom",
    legend.key.width = unit(50, "pt")
  )

ggsave(
  plot = g,
  filename = here::here("geocomputation",
                        "images",
                        "elevation_raster_maps_2.png"),
  height = 1400,
  width = 2000,
  units = "px"
)
Figure 2: Base raster map of Switzerland from {geodata}

Getting outline map of Switzerland {sf}

This code creates two visualizations of Switzerland using ggplot2. The first plot is a vector map generated from the {rnaturalearth} package, displaying Switzerland’s boundaries at a large scale (1:10 million). The second plot combines this vector boundary with a high-resolution raster elevation map from switzerland_raster, highlighting elevation levels up to 6000 meters. The raster data is overlaid with a transparent outline of Switzerland and utilizes a squished colour scale for elevation visualization.

Code 
switzerland_vector <- rnaturalearth::ne_countries(
  country = "Switzerland",
  returnclass = "sf",
  scale = "large"
)

g <- ggplot() +
  geom_sf(data = switzerland_vector) +
  labs(
    title = "Vector Map of Switzerland from {rnaturalearth}",
    subtitle = "Scale of 1:10 million",
    fill = "Elevation (metres)"
  )

ggsave(
  plot = g,
  filename = here::here("geocomputation",
                        "images",
                        "elevation_raster_maps_3.png"),
  height = 1400,
  width = 2000,
  units = "px"
)

g <- ggplot() +
  geom_spatraster(data = switzerland_raster) +
  geom_sf(
    data = switzerland_vector,
    linewidth = 1, 
    fill = NA,
    colour = "black",
    alpha = 0.8
  ) +
  coord_sf(
    crs = 4326
  ) +
  scale_fill_wiki_c(
    limits = c(0, 6000),
    oob = scales::squish
  ) +
  labs(
    title = "Elevation Map of Switzerland",
    subtitle = "Overlaid with a vector map with geom_sf()",
    fill = "Elevation (metres)"
  ) +
  theme(
    legend.position = "bottom",
    legend.key.width = unit(40, "pt")
  )


ggsave(
  plot = g,
  filename = here::here("geocomputation",
                        "images",
                        "elevation_raster_maps_4.png"),
  height = 1400,
  width = 2000,
  units = "px"
)
(a) Switzerland: Vector map using geom_sf{}
(b) Vector map of Switzerland, plotted using {sf}, overlaid with a raster map, plotted using {terra}
Figure 3

Raster Operations: Spatial Subsetting

This code calculates the mean elevation of Switzerland using a raster dataset and filters the raster to show only areas above 2000 meters. It then visualizes these high-altitude regions as a raster map overlaid with Switzerland’s vector boundaries using ggplot2. The map employs a custom color palette, minimal theme, and detailed labels, emphasizing Switzerland’s mountain ranges above 2000 meters elevation.

# Mean elevation above Sea Level in Switzerland
switzerland_raster |> 
  values() |> 
  as_tibble() |> 
  summarise(mean = mean(CHE_elv_msk, na.rm = T)) |> 
  pull(mean)

#> [1] 1289.186

The mean elevation of Switzerland is 1,289.19 metres.

Code 
# Plotting only those portions of Switzerland that are above 2000 metres above sea level
swit_mountains <- switzerland_raster > 2000
swit_mountains <- switzerland_raster[swit_mountains, drop = FALSE]

g <- ggplot() +
  geom_spatraster(data = swit_mountains) +
  geom_sf(
    data = switzerland_vector,
    linewidth = 1, 
    fill = NA,
    colour = "black",
    alpha = 0.8
  ) +
  coord_sf(
    crs = 4326
  ) +
  paletteer::scale_fill_paletteer_c(
    "ggthemes::Brown",
    na.value = "transparent"
    ) +
  labs(
    title = "Switzerland mountain ranges",
    subtitle = "Only showing areas above 2000 metres elevation",
    fill = "Elevation (metres)"
  ) +
  theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.key.width = unit(40, "pt")
  )


ggsave(
  plot = g,
  filename = here::here("geocomputation",
                        "images",
                        "elevation_raster_maps_5.png"),
  height = 1400,
  width = 2000,
  units = "px"
)
Figure 4: Subsetting rasters: Displaying areas in Switzerland with altitude over 2000 metres

Computing distance from highest point in Switzerland

This code identifies Switzerland’s highest elevation point at 4442 meters and filters for all three raster points above 4200 meters. It calculates distances from these points (and removes non-Switzerland areas: this is work-in-progress). Finally, it visualizes these extreme elevations on a map using ggplot2, overlaying the filtered raster data with Switzerland’s vector boundaries. The plot employs a custom brown color palette and highlights the topographic peaks within Switzerland exceeding 4200 meters.

Code 
# Highest Point in Switzerland
switzerland_raster |> 
  values() |> 
  max(na.rm = TRUE)

# It is 4442 metres above sea level.

# Let us select all points above 4200 metres
swit_4200 <- switzerland_raster == 4442
swit_4200 <- switzerland_raster[swit_4200, drop = FALSE]
swit_4200 <- terra::distance(swit_4200)
swit_4200 <- swit_4200 / 1000

# Remove non-Switzerland area
values(swit_4200) <- tibble(
  dist = values(swit_4200),
  actual = values(switzerland_raster)
) |> 
  mutate(
    masked = if_else(
      is.nan(actual),
      NA,
      dist
    )
  ) |> 
  pull(masked)

# Plotting the three points in Switzerland above 4200 m
g <- ggplot() +
  geom_spatraster(data = swit_4200) +
  coord_sf(
    crs = 4326
  ) +
  paletteer::scale_fill_paletteer_c(
    "grDevices::Inferno",
    direction = -1,
    na.value = "transparent"
  ) +
  labs(
    title = "Distance to Switzerland's highest points",
    subtitle = "Distance from the highest points (>4200 m elevation)",
    fill = "Distance (km)"
  ) +
  theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.key.width = unit(40, "pt")
  )


ggsave(
  plot = g,
  filename = here::here("geocomputation",
                        "images",
                        "elevation_raster_maps_6.png"),
  height = 1400,
  width = 2000,
  units = "px"
)
Figure 5: Distance from the highest point in Switzerland

References

Hijmans, Robert J., Márcia Barbosa, Aniruddha Ghosh, and Alex Mandel. 2024. “Geodata: Download Geographic Data.” https://CRAN.R-project.org/package=geodata.