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

Hi, thank you for sharing the code. While trying to reproduce your code as-is, I find that the osm data overlay does not match the base map. There is a small shift. Could you please help me on this, any changes on packages so that your code does not work.
Thank you very much,
Luis

For a better understanding of what is happenin

g I send you an image of whats happening.

OK, it looks like something might be off with the projection one or other of the layers is using.

I haven't looked at this code for a long time, but I'll try to take a look at the weekend.

Thank you so much. I have been trying for a long time to figure out what is the problem.

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.