R code snippet: extract elevation matrix, clip it to a spatial polygon, render with rayshader

andes.R

library(sf)
library(elevatr)
library(raster)
library(rayshader)
library(dplyr)
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

# read in andes shapefile from https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5466557/
# idk exactly what i'm doing yet
# get extent we're interested in
andes_shp <- read_sf("RegionAndina/Region Andina.shp")

# get elevation for this area 
# z = low because otherwise a massive amount of data is downloaded
elevation_data <-get_elev_raster(andes_shp, z=3, src = "aws",
                                 clip='bbox')
# okay, this works
plot(elevation_data)

# amazing solution to only keep info in sf
# from https://gis.stackexchange.com/questions/293993/plotting-and-analyzing-extracted-elevation-data-in-r
library(fasterize)
library(rnaturalearth)
worldmap <- ne_countries(scale = "small", returnclass = "sf")
southam <- worldmap %>%
  filter(region_wb == "Latin America & Caribbean" &
           subregion == "South America" )

andes_shp$POLYID <- 1:nrow(andes_shp)
polymap <- fasterize(andes_shp, elevation_data, field = "POLYID")
## mask out elevation where there's no polygon: set to 1
elevation_data[is.na(values(polymap))] <- 1
plot(elevation_data)

# now repeat but for the outline of south america
southam$POLYID <- 1:nrow(southam)
polymap <- fasterize(southam, elevation_data, field = "POLYID")
elevation_data[is.na(values(polymap))] <- NA # set to NA
plot(elevation_data) # make sure this works

# generate elevation matrix
elmat <- raster_to_matrix(elevation_data) 

# render
elmat %>%
  height_shade(texture = (grDevices::colorRampPalette(c('#8a8a8a','#474656',"#4B7083", "#458293", "#768E80",'#FB7F14', "#ED580E",'#D54505')))(256)) %>%
  add_shadow(ambient_shade(elmat), 0.1) %>%
  plot_3d(elmat, zscale = 150, theta = 10, phi = 36, zoom = .47, windowsize = c(800, 800),
          background = '#aea49e', solid = T,  solidcolor = 'grey20')
# idk why the zscale has to be so intense, but otherwise the map looks a little wild

render_snapshot('andes_light',clear=F)