Code to take DSM data and imagery data to render a 3D map of Madison

madison_dsm_imagery.R

# Load packages ----
library(sf)
library(terra)
library(dplyr)
library(ggplot2)
library(rayshader)
library(rgl)
library(tidyverse)

# Data sources
# https://geodata.wisc.edu/catalog/67d0518f-d029-4d8f-a8ac-41d9138fb30d
# https://gis-countyofdane.opendata.arcgis.com/pages/imagery-download

# elevation data
dem_file <- tempfile()
download.file("https://web.s3.wisc.edu/wsco-wisconsinview/lidar/Dane/Madison_2022_City_Delivery/Raster_DSM_Tiles/Madison_0386.zip", dem_file)
unzip(dem_file)
mad <- rast("Madison_0386.tif")

bbox_mad <- st_bbox(mad)

msn_zip_names <- list.files("data/msn_img/", full.names = TRUE)
walk(msn_zip_names, ~ unzip(., exdir = "data/msn_img/"))
msn_tif_names <- list.files("data/msn_img/", pattern = "tif$", full.names = TRUE)

read_and_aggregate_tif <- function(tif_path) {
  rast(tif_path) |> aggregate(2)
}
msn_mosaic_list <- map(msn_tif_names, ~ read_and_aggregate_tif(.))
msn_mosaic <- exec("mosaic", !!!msn_mosaic_list)
msn_cropped <- crop(msn_mosaic, bbox_mad)

convert_to_array <- function(input_raster) {
  # rayshader needs a matrix/array as input
  # turn raster into array
  
  names(input_raster) = c("r", "g", "b", "i")
  
  matrix_r <- raster_to_matrix(input_raster$r)/256
  matrix_g <- raster_to_matrix(input_raster$g)/256
  matrix_b <- raster_to_matrix(input_raster$b)/256
  
  rgb_array <- array(0, dim = c(nrow(matrix_r), ncol(matrix_r), 3))
  rgb_array[, , 1] = matrix_r  #Red layer
  rgb_array[, , 2] = matrix_g  #Blue layer
  rgb_array[, , 3] = matrix_b  #Green layer
  
  aperm(rgb_array, c(2, 1, 3))
}
msn_array <- convert_to_array(msn_cropped)

msn_rast <- raster_to_matrix(mad)

plot_3d(msn_array, heightmap = msn_rast,
          zscale = 1,
          windowsize = c(1920,1080))
render_highquality("output/msn_with_imagery.png", samples = 256, light = F)