Quick elevation maps with R

relief_maps_rspatial.R

library(sf)
library(terra)

## Define area
coords <- data.frame(x = c(-5.5, -5.5, -5.3, -5.3),
                     y = c(36.7, 36.8, 36.7, 36.8))
coords.sf <- st_as_sf(coords, coords = c("x", "y"), crs = 4326)

## Download elevation data
elev.ras <- elevatr::get_elev_raster(coords.sf, z = 13)
elev <- rast(elev.ras)


## Calculate hillshade
slopes <- terrain(elev, "slope", unit = "radians")
aspect <- terrain(elev, "aspect", unit = "radians")
hs <- shade(slopes, aspect)

## Plot hillshading as basemap
# (here using terra::plot, but could use tmap)
plot(hs, col = gray(0:100 / 100), legend = FALSE, axes = FALSE,
     xlim = c(-5.50, -5.30), ylim = c(36.69, 36.82))
# overlay with elevation
plot(elev, col = terrain.colors(25), alpha = 0.5, legend = FALSE,
     axes = FALSE, add = TRUE)
# add contour lines
contour(elev, col = "grey40", add = TRUE)




#### Overlaying orthoimage ####

## Could use own image, but here downloading w/ maptiles

library(maptiles)
ortho <- get_tiles(ext(-5.50, -5.30, 36.7, 36.8),
                   provider = "Esri.WorldImagery", zoom = 13)

## Plot
plot(hs, col = gray(0:100 / 100), legend = FALSE, axes = FALSE,
     xlim = c(-5.50, -5.30), ylim = c(36.69, 36.82))
# overlay with elevation
plot(ortho, alpha = 150, add = TRUE)




#### 3-D map with rayshader and rayvista ####

library(rayshader)
library(rayvista)

graz.3D <- plot_3d_vista(dem = elev)
render_snapshot(filename = "rayvista.png")

An example showing ocean bathymetry:

library(sf)
library(terra)

## Define area
coords <- data.frame(x = c(-6, -6, -5.2, -5.2),
                     y = c(35.85, 36.1, 35.85, 36.1))
coords.sf <- st_as_sf(coords, coords = c("x", "y"), crs = 4326)

## Download elevation data
elev.ras <- elevatr::get_elev_raster(coords.sf, z = 11)
elev <- rast(elev.ras)
elev[elev > 0] <- NA

## Calculate hillshade
slopes <- terrain(elev, "slope", unit = "radians")
aspect <- terrain(elev, "aspect", unit = "radians")
hs <- shade(slopes, aspect)


## Plot hillshading as basemap
# (here using terra::plot, but could use tmap)
plot(hs, col = gray(0:100 / 100), legend = FALSE, axes = FALSE, colNA = "grey")
# overlay with elevation
plot(elev, col = rev(RColorBrewer::brewer.pal(9, "Blues")), alpha = 0.5, legend = FALSE,
     axes = FALSE, add = TRUE, colNA = "grey")
# add contour lines
contour(elev, col = "grey40", add = TRUE)

Hi, I just updated all packages and seem to be getting a CRS error - notice anything that looks off below?

library(sf)
Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1
library(terra)
terra version 1.4.22

Define area

coords <- data.frame(x = c(-5.5, -5.5, -5.3, -5.3),

•                  y = c(36.7, 36.8, 36.7, 36.8))
  

coords.sf <- st_as_sf(coords, coords = c("x", "y"), crs = 4326)

Download elevation data

elev.ras <- elevatr::get_elev_raster(coords.sf, z = 13)
Error in h(simpleError(msg, call)) :
error in evaluating the argument 'CRSobj' in selecting a method for function 'spTransform': NA
sessionInfo()
R version 4.1.2 (2021-11-01)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 22000)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.1252
[2] LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C
[5] LC_TIME=English_United States.1252

attached base packages:
[1] stats graphics grDevices utils
[5] datasets methods base

other attached packages:
[1] terra_1.4-22 sf_1.0-4

loaded via a namespace (and not attached):
[1] Rcpp_1.0.7 magrittr_2.0.1
[3] units_0.7-2 tidyselect_1.1.1
[5] lattice_0.20-45 R6_2.5.1
[7] rlang_0.4.12 fansi_0.5.0
[9] dplyr_1.0.7 tools_4.1.2
[11] rgdal_1.5-27 grid_4.1.2
[13] KernSmooth_2.23-20 utf8_1.2.2
[15] e1071_1.7-9 DBI_1.1.1
[17] ellipsis_0.3.2 class_7.3-19
[19] digest_0.6.29 assertthat_0.2.1
[21] tibble_3.1.6 lifecycle_1.0.1
[23] crayon_1.4.2 progressr_0.9.0
[25] elevatr_0.4.1 purrr_0.3.4
[27] codetools_0.2-18 vctrs_0.3.8
[29] glue_1.5.1 sp_1.4-6
[31] proxy_0.4-26 compiler_4.1.2
[33] pillar_1.6.4 generics_0.1.1
[35] classInt_0.4-3 pkgconfig_2.0.3

Hi @BrentPease1

That looks again a problem with CRS handling in elevatr/sp.

Does it work if you convert the coords.sf object to Spatial before calling elevatr? That is:

library(sf)
library(terra)

## Define area
coords <- data.frame(x = c(-5.5, -5.5, -5.3, -5.3),
                     y = c(36.7, 36.8, 36.7, 36.8))
coords.sf <- st_as_sf(coords, coords = c("x", "y"), crs = 4326)
coords.sp <- as_Spatial(coords.sf)

## Download elevation data
elev.ras <- elevatr::get_elev_raster(coords.sp, z = 13)   # Note using coords.sp now
elev <- rast(elev.ras)

@Pakillo Sorry, I meant to mention that I also tried that yesterday with the same result:

Error in h(simpleError(msg, call)) :
error in evaluating the argument 'CRSobj' in selecting a method for function 'spTransform': NA

Alright. I can't reproduce the problem in my computer so I'm afraid I can't try proposed solutions.

But I'm guessing the below should work? Here using sp rather than sf, because {elevatr} still uses the former.

library(sp)
library(terra)

## Define area
coords <- data.frame(x = c(-5.5, -5.5, -5.3, -5.3),
                     y = c(36.7, 36.8, 36.7, 36.8))
coordinates(coords) <- c("x", "y")
proj4string(coords) <- CRS("EPSG:4326")

## Download elevation data
elev.ras <- elevatr::get_elev_raster(coords, z = 13)
elev <- rast(elev.ras)

An alternative would be to download elevation data using geodata::elevation_3s(), see https://github.com/rspatial/geodata

@Pakillo oddly, I still get an error. Not sure what is going on? Any advice for troubleshooting?

After running the proj4string() line, I get:

Error in CRS("EPSG:4326") : NA

The problem is with CRS handling in elevatr/sp (see jhollist/elevatr#48). I'm afraid I don't know how to fix it easily.

You can use {geodata} to get elevation data, as an alternative to {elevatr}:

library(sf)
library(terra)

#### Define area ####
coords <- data.frame(x = c(-5.5, -5.5, -5.3, -5.3),
                     y = c(36.7, 36.8, 36.7, 36.8))
coords.sf <- st_as_sf(coords, coords = c("x", "y"), crs = 4326)


#### Download elevation data ####

## Using elevatr
# elev.ras <- elevatr::get_elev_raster(coords, z = 10)
# elev <- rast(elev.ras)

## Using geodata
elev <- geodata::elevation_3s(lon = -5.4, lat = 36.75, path = getwd())
elev <- crop(elev, coords.sf)


#### Calculate hillshade ####
slopes <- terrain(elev, "slope", unit = "radians")
aspect <- terrain(elev, "aspect", unit = "radians")
hs <- shade(slopes, aspect)


#### Make map ####

## Plot hillshading as basemap
# (here using terra::plot, but could use tmap)
plot(hs, col = gray(0:100 / 100), legend = FALSE, axes = FALSE,
     xlim = c(-5.50, -5.30), ylim = c(36.69, 36.82))
# overlay with elevation
plot(elev, col = terrain.colors(25), alpha = 0.5, legend = FALSE,
     axes = FALSE, add = TRUE)
# add contour lines
contour(elev, col = "grey40", add = TRUE)

Hope it helps

I'm getting the same CRS/PROJ errors with elevatr as @BrentPease1 describes, but the geodata alternative worked.
Thanks for sharing this!