Plotting 3d maps of the Netherlands with different water levels

water_NL.R

# Use the great rayshader package (www.rayshader.com)
library(rayshader)

# Download data
zip_loc <- tempfile() 
download.file("http://geodata.nationaalgeoregister.nl/ahn1/extract/ahn1_100m/ahn1_100.tif.zip", 
              zip_loc)
local_tif <- raster::raster(unzip(zip_loc, "ahn_100.tif"))
unlink(zip_loc)

# Create matrix of reasonable size
tif_small   <- raster::aggregate(local_tif, fact = 4, FUN = mean)
elmat_small <- matrix(raster::extract(tif_small, 
                                      raster::extent(tif_small), 
                                      buffer = 1000),
                         nrow = ncol(tif_small),
                         ncol = nrow(tif_small))

# Normalize / threshold
elmat_small[elmat_small < -2000000000] <- -10000
elmat_small[elmat_small < -1000]       <- -1000

# Compute shading layers
ambmat_small <- ambient_shade(elmat_small, zscale = 3)
raymat_small <- ray_shade(elmat_small, zscale = 3, maxsearch = 3000)

# Make images with three different water levels
for (i in c(-.5, -.2, 0)) {
  elmat_small %>%
    sphere_shade(texture = "imhof3", zscale = 3) %>%
    add_shadow(raymat_small, 0.5) %>%
    add_shadow(ambmat_small, 0.5) %>%
    plot_3d(elmat_small, 
            zscale     = 1000, 
            fov        = 0,
            theta      = -30,
            phi        = 45,
            windowsize = c(2560, 2560),
            zoom       = 1,
            water      = TRUE, 
            waterdepth = i)
  
  png(paste0("water_", i, ".png"), width = 2560, height = 2560)
  render_snapshot()
  dev.off()
  
  rgl::rgl.close()
}