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")

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!