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()
#
#

OK, it looks like manually specifying the projection used in the elevation data does the trick.

According to this paper, the data was uploaded using the EPSG:28992 projection.

So, adding the line:

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

after

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

should fix the alignment. I have no idea why it worked before, though I think some changes have happened over the last year or so with regards PROJ4 projections.

I'll update the gist itself when I get a chance. Hope that helps!

EDIT: I should add that I tested this using the most up-to-date versions of the packages listed at the top of the script. I'll update the version numbers there.

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