Recreating Wapo Judiciary Square Helicopter Dataviz

helicopter.R

#Load all the libraries needed
library(sf)
library(dplyr)
library(rayrender)
library(elevatr)
library(rayshader)
library(geojsonsf)
library(raster)

#To recreate the wapo visualization, we need to visualize a few things in 3D: the buildings,
#the ground the buildings are standing on, and the helicopter path.

#Use legacy sf behavior
sf_use_s2(FALSE)

#Load the 3D DC building dataset (download from https://opendata.dc.gov/documents/buildings-in-3d) 
st_read("BldgPly_3D.shp") -> buildings
buildings

#This is what the raw data looks like
buildings$geometry[[1]]

#Set approximate location of helicopter descent from the story
heli_location = c(-77.018929,38.896128)

#Transform the coordinate system from WGS84 (lat/long) to the DC building coordinate system
heli_location |> 
  st_point() |> 
  st_sfc(crs = 4326) |> 
  st_transform(st_crs(buildings)) ->
helo_point

#Point is now in the same coordinate system as the DC building dataset
helo_point

#Need this to filter out only the MULTIPOLYGON entries
is_multipolygon = function(geometry) {
  return(unlist(lapply(geometry, inherits,"MULTIPOLYGON")))
}

#Calculate the centroids and only include buildings within 1000m of judiciary square
buildings |> 
  filter(is_multipolygon(geometry)) |> 
  mutate(centroid =  st_centroid(geometry)) |> 
  mutate(dist = st_distance(centroid, helo_point)) |> 
  filter(dist < units::as_units(1000,"m"))->
judiciary_square_buildings


#This function converts MULTIPOLYGON Z features into a 3D model
convert_multipolygonz_to_obj = function(sfobj, filename) {
  con = file(filename, "w")
  on.exit(close(con))
  num_polygons = nrow(sfobj)
  total_verts = 0
  cat_list = list()
  counter = 1
  for(j in seq_len(num_polygons)) {
    if(j %% 100 == 0) {
      print(j)
    }
    single_geom = sfobj$geometry[[j]]
    number = length(single_geom)
    for(i in seq_len(number)) {
      mat = single_geom[[i]][[1]][-1,]
      text_mat = matrix(sprintf("%0.4f", mat),ncol=ncol(mat), nrow=nrow(mat))
      indices = sprintf("%d", seq_len(nrow(mat))+total_verts)
      cat_list[[counter]] = sprintf("v %s", apply(text_mat,1,paste0, collapse=" "))
      counter = counter + 1
      cat_list[[counter]] = sprintf("f %s", paste0(indices,collapse = " "))
      counter = counter + 1
      total_verts = total_verts + nrow(mat)
    }
  }
  cat(unlist(cat_list),file=con, sep="\n")
}

convert_multipolygonz_to_obj(judiciary_square_buildings, "judiciary.obj")
center_point = st_coordinates(helo_point)

#Render the buildings in 3D
obj_model("judiciary.obj", load_material = F, material=diffuse(color="grey20")) |>  
  group_objects(translate = c(-center_point[1],-center_point[2],0)) |> 
  group_objects(angle=c(-90,0,0)) |> 
  render_scene(lookfrom=c(5000,5000,5000)*1.5, ambient_light = T, fov=10,
               width = 800, height= 800, sample_method="sobol_blue")

#No ground--let's fix that.

#Get the elevation data 
elevation_data = elevatr::get_elev_raster(locations = st_buffer(helo_point, 
                                                                dist=units::as_units(1000,"m")), 
                                          z = 12)

#Get the bounding box for the scene
scene_bbox = st_bbox(judiciary_square_buildings)

#Crop the elevation data to that bounding box
cropped_data = raster::crop(elevation_data, scene_bbox)

#Use rayshader to convert that raster data to a matrix
dc_elevation_matrix = raster_to_matrix(cropped_data)

#Turn the elevation data into a 3D surface and plot using rayshader
dc_elevation_matrix |> 
  sphere_shade(zscale=1/10)|> 
  plot_3d(dc_elevation_matrix,
          solid=T,
          shadow=F, soliddepth=0)

#Save the terrain as a 3D model
save_obj("judiciary_square_elevation.obj")


#Determine the amount we must scale the terrain to match the coordinate system of the map
scale_x = (scene_bbox[3]-scene_bbox[1])/nrow(dc_elevation_matrix)
scale_y = (scene_bbox[4]-scene_bbox[2])/ncol(dc_elevation_matrix)

#The coordinate we must translate the terrain model to (so it matches the building coordinate system)
center_coord = st_coordinates(helo_point)

#Load the helicopter flight path data
heli <- geojson_sf("~/Desktop/two_helis.geojson")
heli
  

#Transform to the coordinate system of the building data
heli_transformed = st_transform(heli, sf::st_crs(buildings))
xy_values = st_coordinates(heli_transformed$geometry)
heli_coords = data.frame(x=xy_values[,1],  y=xy_values[,2], 
                         z=heli$altitude_modeled_meters)

#Filter so we only get flight path data close to the point of interest
heli_coords |> 
  dplyr::filter((center_coord[1]-x)^2 + (center_coord[2]-y)^2 < 500^2) ->
heli_near_values


#Filter again so we only get flight path data before 
heli_near_values |> 
  dplyr::mutate(dist = (center_coord[1]-x)^2 + (center_coord[2]-y)^2) |> 
  dplyr::filter(row_number() < which.min(dist)) |>  
  dplyr::select(x,y,z) ->
heli_right_before

#Two trim to just points of interest
end_at_intersection = heli_right_before[1:21,]


#Render the scene in 3D
obj_model("judiciary.obj", load_material = F, material=diffuse(color="grey20")) |>
  add_object(extruded_path(end_at_intersection, width=2,
                           material=light(intensity=1,importance_sample = F,color="orange"))) |>
  group_objects(translate = c(-center_point[1],-center_point[2],0)) |>
  group_objects(angle=c(-90,0,0)) |>
  add_object(obj_model("judiciary_square_elevation.obj", load_material = FALSE,
                       scale = c(scale_x,1,scale_y), material=diffuse(color="grey20"))) |>
  render_scene(lookfrom=c(5000,5000,5000)*1.5, ambient_light = T, fov=10,
               width = 800, height= 800, sample_method="sobol_blue")

#Translate the camera path back to our rendering coordinate system
camera_path = end_at_intersection - data.frame(x=rep(center_point[1],nrow(end_at_intersection)),
                                               y=rep(center_point[2],nrow(end_at_intersection)),
                                               z=rep(0,nrow(end_at_intersection)))

rot_mat = rayrender:::generate_rotation_matrix(c(90,0,0),c(1,2,3))[1:3,1:3]
rotated_path = t(apply(camera_path,1,(function(x) x %*% rot_mat)))

#Generate the camera motion
heli_motion = generate_camera_motion(rotated_path, lookats=rotated_path, frames=360, 
                                     type="bezier", fovs=80,
                                     constant_step = T, offset_lookat=10)


#Generate an animation from the POV of the helicopter
obj_model("judiciary.obj", load_material = F, material=diffuse(color="grey50")) |> 
  group_objects(translate = c(-center_point[1],-center_point[2],0)) |> 
  group_objects(angle=c(-90,0,0)) |> 
  add_object(obj_model("judiciary_square_elevation.obj", load_material = FALSE,
                       scale = c(scale_x,1,scale_y), material=diffuse(color="grey20"))) |>
  render_animation( camera_motion = heli_motion,samples=1, min_variance = 1e-6,
                    ambient = TRUE, 
                    clamp_value=10, 
                    width=400,height=400,sample_method="sobol_blue", )

#Make it a rollercoaster!
obj_model("judiciary.obj", load_material = F, material=diffuse(color="grey50")) |> 
  add_object(extruded_path(end_at_intersection, width=2, z=-5,
                           material=light(intensity=1,importance_sample = F,color="orange"))) |>
  group_objects(translate = c(-center_point[1],-center_point[2],0)) |> 
  group_objects(angle=c(-90,0,0)) |> 
  add_object(obj_model("judiciary_square_elevation.obj", load_material = FALSE,
                       scale = c(scale_x,1,scale_y), material=diffuse(color="grey20"))) |>
  render_animation( camera_motion = heli_motion,samples=1, min_variance = 1e-6,
                    ambient = TRUE, 
                    clamp_value=10,
                    width=400,height=400,sample_method="sobol_blue", )