Skip to content

Instantly share code, notes, and snippets.

@nevrome
Last active January 17, 2021 08:00
Show Gist options
  • Save nevrome/82b3107486985e59f6e7d6ea11c0c9c7 to your computer and use it in GitHub Desktop.
Save nevrome/82b3107486985e59f6e7d6ea11c0c9c7 to your computer and use it in GitHub Desktop.
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
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()
)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment