R script to get a digital elevation model with sf and elevatr, create a hillshaded raster with rayshader and plot it with ggplot2 and ggspatial

elevatr_rayshader_ggplot2.R

library(magrittr)
library(ggplot2)
library(rayshader)

#### get some point data of c14 dates from Switzerland (with c14bazAAR) ####
c14dates <- c14bazAAR::get_RADON() %>%
  dplyr::filter(
    country == "Switzerland",
    culture == "Cortaillod"
  ) %>% 
  dplyr::group_by(site) %>%
  dplyr::arrange(desc(c14age)) %>%
  dplyr::filter(dplyr::row_number() == 1) %>%
  dplyr::ungroup() %>%
  dplyr::select(
    labnr, c14age, c14std, lat, lon, site
  ) %>%
  c14bazAAR::as.c14_date_list() %>%
  c14bazAAR::as.sf() %>%
  sf::st_transform(crs = 32632)
  
#### get Switzerland country border (with rnaturalearth) ####
switzerland <- rnaturalearth::ne_countries(country = 'switzerland', scale = 'large', returnclass = "sf") %>%
  sf::st_transform(32632)

#### determine area of interest (with sf) ####
area_sf <- sf::st_bbox(c14dates) %>% sf::st_as_sfc() %>% sf::st_buffer(1000)
area_sp <- area_sf %>% sf::as_Spatial()

#### get digital elevation data (with elevatr) ####
dem <- elevatr::get_elev_raster(locations = area_sp, prj = sf::st_crs(32632), z = 7, clip = "bbox")

#### prepare hillshaded raster (with rayshader) ####
elmat <- matrix(
  raster::extract(dem, raster::extent(dem)),
  nrow = ncol(dem), 
  ncol = nrow(dem)
)

raymat <- ray_shade(elmat)
ambmat <- ambient_shade(elmat)

hillshade <- elmat %>%
  sphere_shade(texture = "imhof1") %>%
  add_shadow(raymat) %>%
  add_shadow(ambmat)

#### tranform hillshade to georeferenced raster (with raster) ####
e <- raster::extent(dem)
hillshade_raster <- raster::brick(hillshade, xmn = e[1], xmx = e[2], ymn = e[3], ymx = e[4], crs = raster::crs(dem))

#### plot raster and points (with ggplot2) ####
ggplot() +
  ggspatial::layer_spatial(hillshade_raster) +
  geom_sf(
    data = switzerland,
    colour = "red",
    fill = NA
  ) +
  geom_sf(
    data = c14dates,
    mapping = aes(colour = data.c14age, size = data.c14age),
    show.legend = 'point'
  ) +
  geom_sf_label(
    data = c14dates[1,],
    mapping = aes(label = data.site),
    size = 3,
    alpha = 0.3,
    nudge_y = -5000
  ) +
  theme_bw() +
  theme(
    axis.title = element_blank()
  )