Elevation data processing and analysis in R

1. Introduction
2. Setup
3. Introduction to elevatr
4. Get elevation data for your department
5. Visualize the elevation data
6. Bibliography
- Reproducibility

1. Introduction

This notebook illustrates several functionalities to obtain, process and visualize digital elevation models (DEMs) in R. It aims to help Geomatica Basica students at Universidad Nacional de Colombia to get acquainted with DEMs and geomorphometric variables.

For geospatial data management, there are two main packages in R:

the simple features (sf) library that provides functionalities for vector data
the terra library that provides functionalities for raster data.

In addition, for interactive mapping, a good option is the leaflet library.

Before proceeding, a few tips for writing your own notebook follow:

Write and run every chunk step-by-step
In case of errors, read the message, interpret it and try to fix it
In case the error persists, ask the lecturer before the deadline
Once your html file is produced, it is time to publish it at RPubs
After publishing the notebook, review it, fix poor quality wording, and republish it using the update option
In case of trouble, please review the R documentation.

When creating your notebook, make sure to add as many text as needed in order to: (i) explain each code chunk; and (ii) describe the corresponding output. These two criteria will be comsidered to assessing the quality of your notebook.

2. Setup

The required packages should be installed in advance from the Console.

## SETUP
# INSTALL THESE PACKAGES FROM THE CONSOLE, NOT FROM THIS CHUNK
# paquetes = c("terra","elevatr","sf", "leaflet")
# install.packages(paquetes)

It is advisable to start cleaning the memory:

rm(list=ls())

Now, let’s load the libraries:

library(elevatr)
library(sf)
library(leaflet)
library(terra)

When running a chunk, pay attention to any message indicating error or warning eventual conflicts.

Conflicts may be caused by functions, from differents packages, which have the same names (e.g. terra::extract and tidyr::extract). In order to avoid such problems, we need to call several functions using the long way, i.e. package_name::function

3. Introduction to elevatr

Elevation data is used for a wide array of applications, including, for example, visualization, hydrology, and ecological modelling. Gaining access to these data in R has not had a single interface, is made available through functions across many packages, or requires local access to the data. This is no longer required as a variety of APIs now exist that provide programmatic access to elevation data. The elevatr package was written to standarize access to elevation data from web APIs.

To access raster elevation data (e.g. a DEM) the elevatr package uses the Amazon Web Services Terrain Tiles. We will explore this functionality in this notebook.

There are several sources for digital elevation models such as the Shuttle Radar Topography Mission (SRTM), the USGS National Elevation Dataset (NED), Global DEM (GDEM), and others. Each of these DEMs has pros and cons for their use. Prior to its closure in January of 2018, Mapzen combined several of these sources to create a synthesis elevation product that utilizes the best available elevation data for a given region at given zoom level. Although closed, these data compiled by Mapzen are currently accessible through the Terrain Tiles on Amazon Web Services (AWS).

The main function is get_elev_raster() which needs a grpspatial object such a data.frame with x (longitude) and y (latitude) locations as the first two columns, any simple feature (sf) object, or any raster object and it returns a RasterLayer of the tiles that overlap the bounding box of the input. If multiple tiles are retrieved, the resultant output is a merged Raster Layer.

Using get_elev_raster() to access the Terrain Tiles on AWS.

Note that the geospatial object requires a CRS (i.e. a coordinate reference system). The zoom level (z) defaults to 9 (a trade off between resolution and time for download). You could start trying a higher zoom level (e.g. 10).

4. Get elevation data for your department

First, we will load a shapefile or a geopackage representing municipalities of our department. In this notebook, I will use a geopackage with the official boundaries of Cordoba.

Let’s read the geospatial file using a function provided by the sf package:

4.1 Import department boundaries

As this notebook lives together to a folder named cordoba, it is simple to get a list of the available data there.

list.files("./cordoba")

## [1] "cordoba_cities.gpkg"      "cordoba_roads.gpkg"      
## [3] "cordoba_water_areas.gpkg" "cordoba_water_lines.gpkg"
## [5] "cordoba.gpkg"             "cordoba.qgz"             
## [7] "elev_cordoba_z10.tif"

Now, we will read the needed geopackage using the sf library:

(munic <-  sf::st_read("cordoba/cordoba.gpkg"))

## Reading layer `mgn_adm_mpio_grafico' from data source 
##   `/Users/ials/Documents/unal/GB2024/data/cordoba/cordoba.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 30 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -76.51453 ymin: 7.347087 xmax: -74.78095 ymax: 9.447748
## Geodetic CRS:  MAGNA-SIRGAS

ABCDEFGHIJ0123456789

dpto_ccdgo <chr>	mpio_ccdgo <chr>	mpio_cdpmp <chr>	mpio_cnmbr <chr>	geom <s_MULTIP>
23	001	23001	MONTERÍA	<s_MULTIP>
23	068	23068	AYAPEL	<s_MULTIP>
23	079	23079	BUENAVISTA	<s_MULTIP>
23	090	23090	CANALETE	<s_MULTIP>
23	162	23162	CERETÉ	<s_MULTIP>
23	168	23168	CHIMÁ	<s_MULTIP>
23	182	23182	CHINÚ	<s_MULTIP>
23	189	23189	CIÉNAGA DE ORO	<s_MULTIP>
23	300	23300	COTORRA	<s_MULTIP>
23	350	23350	LA APARTADA	<s_MULTIP>

Note the CRS of this dataset: MAGNA-SIRGAS with no details. As we want more info about CRS we will ask:

(cordoba_crs = st_crs(st_geometry(munic)))

## Coordinate Reference System:
##   User input: MAGNA-SIRGAS 
##   wkt:
## GEOGCRS["MAGNA-SIRGAS",
##     DATUM["Marco Geocentrico Nacional de Referencia",
##         ELLIPSOID["GRS 1980",6378137,298.257222101,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["Horizontal component of 3D system."],
##         AREA["Colombia - onshore and offshore. Includes San Andres y Providencia, Malpelo Islands, Roncador Bank, Serrana Bank and Serranilla Bank."],
##         BBOX[-4.23,-84.77,15.51,-66.87]],
##     ID["EPSG",4686]]

Take note of the attribute what stores municipalities’ names. In my case, it is mpio_cnmbr. We will need it later for plotting.

Note also that, in my case, the municipalities do not have any attribute with the corresponding area. Thus, it is convenient to compute it:

# This code calculates area in square meters
# Can you tell us why?
munic$area_m2 = sf::st_area(munic)

Now, the area expressed in square kilometers:

# This code calculates area in square kilometers
# Can you tell us why?
munic$area = munic$area_m2/1000000
# Now, let's round the area to two digits
munic$area = signif(munic$area, digits = 6)

Let’s call again our municipalities object:

munic

ABCDEFGHIJ0123456789

	dpto_ccdgo <chr>	mpio_ccdgo <chr>	mpio_cdpmp <chr>	mpio_cnmbr <chr>	geom <s_MULTIP>
1	23	001	23001	MONTERÍA	<s_MULTIP>
2	23	068	23068	AYAPEL	<s_MULTIP>
3	23	079	23079	BUENAVISTA	<s_MULTIP>
4	23	090	23090	CANALETE	<s_MULTIP>
5	23	162	23162	CERETÉ	<s_MULTIP>
6	23	168	23168	CHIMÁ	<s_MULTIP>
7	23	182	23182	CHINÚ	<s_MULTIP>
8	23	189	23189	CIÉNAGA DE ORO	<s_MULTIP>
9	23	300	23300	COTORRA	<s_MULTIP>
10	23	350	23350	LA APARTADA	<s_MULTIP>

4.2 Get centroids of municipalities

Then, let’s create a new object with centroids of municipalities:

(centers <- st_centroid(munic))

ABCDEFGHIJ0123456789

	dpto_ccdgo <chr>	mpio_ccdgo <chr>	mpio_cdpmp <chr>	mpio_cnmbr <chr>	geom <sf_POINT>
1	23	001	23001	MONTERÍA	<sf_POINT>
2	23	068	23068	AYAPEL	<sf_POINT>
3	23	079	23079	BUENAVISTA	<sf_POINT>
4	23	090	23090	CANALETE	<sf_POINT>
5	23	162	23162	CERETÉ	<sf_POINT>
6	23	168	23168	CHIMÁ	<sf_POINT>
7	23	182	23182	CHINÚ	<sf_POINT>
8	23	189	23189	CIÉNAGA DE ORO	<sf_POINT>
9	23	300	23300	COTORRA	<sf_POINT>
10	23	350	23350	LA APARTADA	<sf_POINT>

Note that the geom attribute encodes the point geometry. Thus, we need to do some work to split each into their x, y coordinates:

centers$x = st_coordinates(centers)[,1]

centers$y = st_coordinates(centers)[,2]

Let’s check the output:

centers

ABCDEFGHIJ0123456789

	dpto_ccdgo <chr>	mpio_ccdgo <chr>	mpio_cdpmp <chr>	mpio_cnmbr <chr>	geom <sf_POINT>
1	23	001	23001	MONTERÍA	<sf_POINT>
2	23	068	23068	AYAPEL	<sf_POINT>
3	23	079	23079	BUENAVISTA	<sf_POINT>
4	23	090	23090	CANALETE	<sf_POINT>
5	23	162	23162	CERETÉ	<sf_POINT>
6	23	168	23168	CHIMÁ	<sf_POINT>
7	23	182	23182	CHINÚ	<sf_POINT>
8	23	189	23189	CIÉNAGA DE ORO	<sf_POINT>
9	23	300	23300	COTORRA	<sf_POINT>
10	23	350	23350	LA APARTADA	<sf_POINT>

4.3 Get the elevation data

Now, it is time to get the DEM:

elevation <- get_elev_raster(munic, z = 10)

## Mosaicing & Projecting

## Note: Elevation units are in meters.

What is the elevation object?

elevation

## class      : RasterLayer 
## dimensions : 3568, 3090, 11025120  (nrow, ncol, ncell)
## resolution : 0.000682546, 0.000682546  (x, y)
## extent     : -76.64062, -74.53156, 7.013738, 9.449062  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs 
## source     : file15c314e519ec.tif 
## names      : file15c314e519ec

Take note of the characteristics of the elevation object. And understand each one. What is its spatial resolution (in meters)?

Notice a few things about this DEM:

dimensions: the “size” of the file in pixels
nrow, ncol: the number of rows and columns in the data (imagine a spreadsheet or a matrix).
ncells: the total number of pixels or cells that make up the raster.
resolution: the size of each pixel (in decimal degrees in this case). Let’s remind that 1 decimal degree represents about 111.11 km at the Equator.
extent: the spatial extent of the raster. This value will be in the same coordinate units as the coordinate reference system of the raster.
crs: the coordinate reference system string for the raster. This raster is in geographic coordinates with a datum of WGS 84.

In your notebook, write down the cell size of your raster (i.e. the spatial resolution), in meters.

Note that elevation is a temporary raster layer (i.e. it exists only in-memory). However, it is possible to write the DEM to your local directory using the writeRaster function provided by the rgdal library:

writeRaster(elevation, "cordoba/elev_cordoba_z10.tif", overwrite=TRUE)

In my case, elevation data at zoom 12 was about 340.3 MB. When zoom was 10, the file size lowered to about 28.5 MB.

Write down the size of your DEM.

Now, we need to convert the elevation to a terra’s raster object:

(elevation2 <- terra::rast(elevation))

## class       : SpatRaster 
## dimensions  : 3568, 3090, 1  (nrow, ncol, nlyr)
## resolution  : 0.000682546, 0.000682546  (x, y)
## extent      : -76.64062, -74.53156, 7.013738, 9.449062  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs 
## source      : file15c314e519ec.tif 
## name        : file15c314e519ec

5. Visualize the elevation data

First, we will need a color palette:

pal <- colorNumeric(c("cyan", "forestgreen","yellow","tan","orange", "brown"), values(elevation),
  na.color = "transparent")

Now, we will reduce the spatial resolution of the data:

(elevation3 <- terra::aggregate(elevation2, fact = 2))

## |---------|---------|---------|---------|=========================================

## class       : SpatRaster 
## dimensions  : 1784, 1545, 1  (nrow, ncol, nlyr)
## resolution  : 0.001365092, 0.001365092  (x, y)
## extent      : -76.64062, -74.53156, 7.013738, 9.449062  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs 
## source(s)   : memory
## name        : file15c314e519ec 
## min value   :          -1238.5 
## max value   :           3711.0

Furthermore, we will clip (or crop) the elevation3 object in order to fit our department’s boundaries:

elevation4 <- terra::crop(elevation3, munic, mask=TRUE)

Now, it is plotting time. Make sure you change several options in the code below in order to match the attributes’ values in your municipalities object:

leaflet(munic) %>% addTiles() %>% setView(-75.9, 8.7, 9) %>% 
    addPolygons(color = "white", weight = 1.0, smoothFactor = 0.5,
    opacity = 0.25, fillOpacity = 0.15,
    popup = paste("Municipio: ", munic$mpio_cnmbr, "<br>",
                          "Km2: ", munic$area, "<br>")) %>%
  addLabelOnlyMarkers(data = centers,
                    lng = ~x, lat = ~y, label = ~mpio_cnmbr,
                    labelOptions = labelOptions(noHide = TRUE, direction = 'top', textOnly = TRUE, textsize = "10px")) %>%
  addRasterImage(elevation4, colors = pal, opacity = 0.9)  %>%
  addLegend(pal = pal, values = values(elevation),
    title = "Elevation data for Cordoba (mts)")

Click on each municipality to get its name and its area [km2].

Your are free to experiment with other plotting options. You may benefit of visiting this site.

That’s all for now!

6. Bibliography

If you reuse this notebook’s code please cite this work as:

Lizarazo, I., 2025. Elevation data processing and analysis in R. Available at https://rpubs.com/ials2un/otro_dem

Reproducibility

sessionInfo()

## R version 4.3.2 (2023-10-31)
## Platform: x86_64-apple-darwin20 (64-bit)
## Running under: macOS Sonoma 14.7.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/Bogota
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] terra_1.7-78   leaflet_2.2.1  sf_1.0-17      elevatr_0.99.0
## 
## loaded via a namespace (and not attached):
##  [1] s2_1.1.7           sass_0.4.9         class_7.3-22       KernSmooth_2.23-22
##  [5] lattice_0.22-6     hms_1.1.3          digest_0.6.37      magrittr_2.0.3    
##  [9] evaluate_1.0.0     grid_4.3.2         fastmap_1.2.0      jsonlite_1.8.9    
## [13] progress_1.2.3     e1071_1.7-16       DBI_1.2.3          httr_1.4.7        
## [17] purrr_1.0.2        crosstalk_1.2.1    scales_1.3.0       slippymath_0.3.1  
## [21] codetools_0.2-19   jquerylib_0.1.4    cli_3.6.3          rlang_1.1.4       
## [25] crayon_1.5.3       units_0.8-5        munsell_0.5.1      cachem_1.1.0      
## [29] yaml_2.3.10        tools_4.3.2        raster_3.6-26      colorspace_2.1-1  
## [33] curl_5.2.1         mime_0.12          png_0.1-8          vctrs_0.6.5       
## [37] R6_2.5.1           proxy_0.4-27       lifecycle_1.0.4    classInt_0.4-10   
## [41] htmlwidgets_1.6.4  pkgconfig_2.0.3    progressr_0.14.0   bslib_0.8.0       
## [45] glue_1.8.0         Rcpp_1.0.13        xfun_0.47          rstudioapi_0.16.0 
## [49] knitr_1.48         farver_2.1.2       htmltools_0.5.8.1  rmarkdown_2.28    
## [53] wk_0.9.3           compiler_4.3.2     prettyunits_1.2.0  sp_2.1-4

Elevation data processing and analysis in R

Ivan Lizarazo

15.01.2025

1. Introduction

2. Setup

3. Introduction to elevatr

4. Get elevation data for your department

4.1 Import department boundaries

4.2 Get centroids of municipalities

4.3 Get the elevation data

5. Visualize the elevation data

6. Bibliography

Reproducibility