Downloading elevation data in R

Installing the whatarelief package

We will be using a relatively new R package called whatarelief. This package is hosted on Github, but is not currently on CRAN so we have to use the devtools package to install (you may need to install this package first).

The primary function in this package is called elevation(), which downloads a elevation raster directly into our R session. This raster can be automatically formatted for specific locations using a reference extent, shapefile, or other raster file.

The package documentation provide all the information you’ll need to get started, including how to select difference elevation sources. It’s worth looking through the documentation, as there are some quarks.

library(whatarelief)
image(im <- elevation())

[1] "/vsicurl/https://public.services.aad.gov.au/datasets/science/GEBCO_2019_GEOTIFF/GEBCO_2019.tif"

image(t(im[nrow(im):1, ]))

The next section will demonstrate a brief example for getting elevation data from Zambia

Example: Elevation raster for Zambia

To get started, let’s load some useful packages for working with rasters and shapefiles, then we will download a shapefile for Zambia.

library(raster)   # Raster package
library(sf)       # Shapefile package
library(PATHtoolsZambia)  # PATH data for Zambia

# Load reference shapefile
shp <- retrieve("province-shp")

# Load reference raster
rst <- retrieve("grid3-pop-rescaled")

# Visualize
plot(rst, col = viridis::plasma(100))
plot(st_geometry(shp), add = T)

To use a reference raster to select an area for interest, all we need to do is include the raster object in the elevation() functions.

zm_elv1 <- elevation(rst)

[1] "/vsicurl/https://public.services.aad.gov.au/datasets/science/GEBCO_2019_GEOTIFF/GEBCO_2019.tif"

plot(zm_elv1)
plot(st_geometry(shp), add = T)

Conveniently, this elevation raster is at the same resolution and alignment as our reference.

compareRaster(rst, zm_elv1)

[1] TRUE

stack(rst, zm_elv1)

class      : RasterStack 
dimensions : 1177, 1406, 1654862, 2  (nrow, ncol, ncell, nlayers)
resolution : 0.008333333, 0.008333333  (x, y)
extent     : 21.99542, 33.71208, -18.07958, -8.27125  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
names      : grid3_1km_rescaled.1, grid3_1km_rescaled.2 
min values :             1.042396,           118.241821 
max values :            30135.733,             2726.069 

To test this functionality, let’s aggregate our reference raster and get a new elevation raster.

# Aggregate by a factor of 5
rst_agg <- aggregate(rst, 5, fun = "sum")
plot(rst_agg, col = viridis::plasma(100))
plot(st_geometry(shp), add = T)

# Pull aggregated version of elevation 
zm_elv2 <- elevation(rst_agg)

[1] "/vsicurl/https://public.services.aad.gov.au/datasets/science/GEBCO_2019_GEOTIFF/GEBCO_2019.tif"

plot(zm_elv2)
plot(st_geometry(shp), add = T)

# Compare
compareRaster(rst_agg, zm_elv2)

[1] TRUE

stack(rst_agg, zm_elv2)

class      : RasterStack 
dimensions : 236, 282, 66552, 2  (nrow, ncol, ncell, nlayers)
resolution : 0.04166667, 0.04166667  (x, y)
extent     : 21.99542, 33.74542, -18.10458, -8.27125  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
names      : grid3_1km_rescaled.1, grid3_1km_rescaled.2 
min values :             1.065404,           131.133560 
max values :           226794.442,             2450.497 

try(compareRaster(zm_elv1, zm_elv2))

Error in compareRaster(zm_elv1, zm_elv2) : different extent

It is also possible to download rasters using just the extent as a numeric vector, which means it is also possible to use a shapefile. However, you will also have to provide the projection and potential other spatial information, which may be more difficult than just using a reference raster.

Citation

BibTeX citation:

@online{millar2022,
  author = {Justin Millar},
  title = {Downloading Elevation Data in {R}},
  date = {2022-07-13},
  url = {https://github.com/PATH-Global-Health/MNTD-tech-docs/posts/elevation-rasters},
  langid = {en}
}

For attribution, please cite this work as:

Justin Millar. 2022. “Downloading Elevation Data in R.” July 13, 2022. https://github.com/PATH-Global-Health/MNTD-tech-docs/posts/elevation-rasters.