Skip to content

Instantly share code, notes, and snippets.

@norwegianblueparrot
Last active November 25, 2023 03:40
Show Gist options
  • Save norwegianblueparrot/b9d5d48f2d591d78a14320bf17459cc5 to your computer and use it in GitHub Desktop.
Save norwegianblueparrot/b9d5d48f2d591d78a14320bf17459cc5 to your computer and use it in GitHub Desktop.
A script to produce a 3D render of Fort Bourtange from LiDAR data
## ----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()
#
#
@lbarqueira
Copy link

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

@lbarqueira
Copy link

For a better understanding of what is happenin
image
g I send you an image of whats happening.

@norwegianblueparrot
Copy link
Author

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.

@lbarqueira
Copy link

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

@norwegianblueparrot
Copy link
Author

norwegianblueparrot commented Oct 22, 2023

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.

image

@lbarqueira
Copy link

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment