A script to produce a 3D render of Fort Bourtange from LiDAR data

bourtange.R

## ----setup,include=FALSE,message=FALSE,warning=FALSE-----------------------------------------------------------------------------------------------------------------------------------------------------
library("here") # 1.0.1
library("sf") # 1.0-14
library("rgl") # 1.2.1
library("geoviz") # 0.2.2
library("raster") # 3.6-26
library("rayshader") # 0.35.7
library("osmdata") # 0.2.5
library("tidyverse") # 2.0.0

## ----load_elevation_data---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# set this path to wherever you've saved the elevation data - it can be downloaded from https://download.pdok.nl/rws/ahn3/v1_0/05m_dsm/R_13DN2.ZIP
elevation_data <- here::here("data/elevation_data/bourtange")

# this data can be downloaded from
elevation_file <- elevation_data %>%
  file.path("R_13DN2.TIF")

elevation_raster <- elevation_file %>%
  raster::raster()

raster::crs(elevation_raster) <- "EPSG:28992"

raster_crs <- elevation_raster %>%
  raster::crs()

## ----zoom_to_fortress------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
xmin <- 7.1869551422
xmax <- 7.1971968768
ymin <- 53.0042753099
ymax <- 53.0105807242

cropped_extent <- raster::extent(c(xmin, xmax, ymin, ymax)) %>%
  sf::st_bbox(crs = 4326) %>%
  sf::st_as_sfc() %>%
  sf::st_transform(crs = raster_crs) %>%
  sf::st_bbox() %>%
  raster::extent()

cropped_raster <- elevation_raster %>%
  raster::crop(y = cropped_extent)

plot(cropped_raster)

## ----scale_matrix----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
scaling_factor <- 1 # set > 1 to reduce size of matrix

raster_matrix <- cropped_raster %>%
  rayshader::raster_to_matrix(verbose = FALSE) %>%
  rayshader::resize_matrix(scale = 1 / scaling_factor)

# set missing values to -1 (assume this is water)
raster_matrix[is.na(raster_matrix)] <- -1

## ----calculate_zscale------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
zscale_scaling_factor <- 1

estimated_zscale <- cropped_raster %>%
  geoviz::raster_zscale()

# this can be scaled to exaggerate elevation
scaled_zscale <- estimated_zscale * zscale_scaling_factor

## ----build_base_map--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
base_map <- raster_matrix %>%
  rayshader::height_shade() %>%
  rayshader::add_overlay(rayshader::sphere_shade(raster_matrix, colorintensity = 5),
                         alphalayer = 0.5)  %>%
  rayshader::add_overlay(
    rayshader::sphere_shade(raster_matrix,
                            texture = "desert",
                            colorintensity = 5),
    alphalayer = 0.5
  ) %>%
  rayshader::add_shadow(rayshader::lamb_shade(raster_matrix), 0) %>%
  rayshader::add_shadow(rayshader::ambient_shade(raster_matrix), 0) %>%
  rayshader::add_shadow(
    rayshader::texture_shade(
      raster_matrix,
      detail = 8 / 10,
      contrast = 9,
      brightness = 11
    ),
    0.1
  )

## ----get_osm_overlays------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# yes, this is basically the same code copied and pasted - don't tell @hadley
bbox_long_lat <- sf::st_bbox(cropped_raster) %>%
  sf::st_as_sfc() %>%
  sf::st_transform(crs = 4326) %>%
  sf::st_bbox()

# get highways
bourtange_footways <- osmdata::opq(bbox_long_lat) %>%
  osmdata::add_osm_feature("highway") %>%
  osmdata::osmdata_sf()

footway_geom <- bourtange_footways$osm_lines %>%
  sf::st_transform(crs = raster_crs)

footway_overlay <- footway_geom %>%
  rayshader::generate_line_overlay(
    extent = cropped_extent,
    heightmap = raster_matrix,
    color = "gainsboro",
    linewidth = 2
  )

# get trees
bourtange_trees <- osmdata::opq(bbox_long_lat) %>%
  osmdata::add_osm_feature("landuse", "forest") %>%
  osmdata::osmdata_sf()

trees_geom <- bourtange_trees$osm_polygons %>%
  sf::st_transform(crs = raster_crs)

tree_overlay <- trees_geom %>%
  rayshader::generate_line_overlay(extent = cropped_extent,
                                   heightmap = raster_matrix,
                                   color = "darkolivegreen")

# get grassy areas
bourtange_grass <- osmdata::opq(bbox_long_lat) %>%
  osmdata::add_osm_feature("landuse", "grass") %>%
  osmdata::osmdata_sf()

grass_geom <- bourtange_grass$osm_polygons %>%
  sf::st_transform(crs = raster_crs)

grass_overlay <- grass_geom %>%
  rayshader::generate_line_overlay(extent = cropped_extent,
                                   heightmap = raster_matrix,
                                   color = "darkolivegreen")

# get buildings
bourtange_buildings <- osmdata::opq(bbox_long_lat) %>%
  osmdata::add_osm_feature("building") %>%
  osmdata::osmdata_sf()

buildings_geom <- bourtange_buildings$osm_polygons  %>%
  sf::st_transform(crs = raster_crs)

building_overlay <- buildings_geom %>%
  rayshader::generate_polygon_overlay(extent = cropped_extent,
                                      heightmap = raster_matrix,
                                      palette = "sienna1")

# get windmill
windmill <- osmdata::opq(bbox_long_lat) %>%
  osmdata::add_osm_feature("man_made", "windmill") %>%
  osmdata::osmdata_sf()

windmill_geom <- windmill$osm_polygons  %>%
  sf::st_transform(crs = raster_crs)

windmill_overlay <- windmill_geom %>%
  rayshader::generate_polygon_overlay(extent = cropped_extent,
                                      heightmap = raster_matrix,
                                      palette = "darkgrey")

## ----render_overhead_map----------------------------------------------------------------------------------------------------------------------------------------------------------------------
base_map %>%
  rayshader::add_overlay(overlay = footway_overlay, alphalayer = 0.95) %>%
  rayshader::add_overlay(overlay = building_overlay, alphalayer = 0.75) %>%
  rayshader::add_overlay(overlay = tree_overlay, alphalayer = 0.6) %>%
  rayshader::add_overlay(overlay = grass_overlay, alphalayer = 0.6) %>%
  rayshader::add_overlay(overlay = windmill_overlay, alphalayer = 0.85) %>%
  rayshader::plot_3d(
    raster_matrix,
    water = TRUE,
    watercolor = "dodgerblue",
    windowsize = c(1200, 1200),
    zscale = estimated_zscale
  )

rayshader::render_camera(theta = 0, phi = 89.999, zoom = 0.75) # setting phi = 90 seems to cause problems for rayshader?

rayshader::render_highquality(
  filename = "bourtange_overhead.png",
  lightaltitude = c(22.5, 75),
  lightintensity = c(1000, 200),
  lightdirection = c(250, 0),
  print_scene_info = TRUE,
  samples = 100,
  parallel = TRUE,
  min_variance = 0,
  verbose = TRUE
)
rgl::rgl.close()

## ----render_zoomed_map-------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# base_map %>%
#   rayshader::add_overlay(overlay = footway_overlay, alphalayer = 0.95) %>%
#   rayshader::add_overlay(overlay = building_overlay, alphalayer = 0.75) %>%
#   rayshader::add_overlay(overlay = tree_overlay, alphalayer = 0.6) %>%
#   rayshader::add_overlay(overlay = grass_overlay, alphalayer = 0.6) %>%
#   rayshader::add_overlay(overlay = windmill_overlay, alphalayer = 0.85) %>%
#   rayshader::plot_3d(
#     raster_matrix,
#     water = TRUE,
#     watercolor = "dodgerblue",
#     windowsize = c(1200, 1200),
#     zscale = estimated_zscale
#   )
#
# rayshader::render_camera(theta = 50, phi = 30, zoom = 0.2, fov = 135)
#
# rayshader::render_depth(focus = 0.9, focallength = 400, aberration = 0.4)
#
# rayshader::render_highquality(
#   filename = "bourtange_zoomed.png",
#   lightaltitude = c(22.5, 75),
#   lightintensity = c(1000, 200),
#   lightdirection = c(250, 0),
#   print_scene_info = TRUE,
#   samples = 500,
#   parallel = TRUE,
#   min_variance = 0,
#   verbose = TRUE
# )
# rgl::rgl.close()
#
#

Thank you very much for spending your time to help me on solving this issue.