if(run_documentation()) {
# Load and visualize building footprints from Open Street Map
library(osmdata)
library(sf)
library(raster)
osm_bbox = c(-121.9472, 36.6019, -121.9179, 36.6385)
#Get buildings from OpenStreetMap
opq(osm_bbox) |>
add_osm_feature("building") |>
osmdata_sf() ->
osm_data
#Get roads from OpenStreetMap
opq(osm_bbox) |>
add_osm_feature("highway") |>
osmdata_sf() ->
osm_road
#Get extent
building_polys = osm_data$osm_polygons
osm_dem = elevatr::get_elev_raster(building_polys, z = 11, clip = "bbox")
e = extent(building_polys)
# Crop DEM, but note that the cropped DEM will have an extent slightly different than what's
# specified in `e`. Save that new extent to `new_e`.
osm_dem |>
crop(e) |>
extent() ->
new_e
osm_dem |>
crop(e) |>
raster_to_matrix() ->
osm_mat
#Visualize areas less than one meter as water (approximate tidal range)
osm_mat[osm_mat <= 1] = -2
osm_mat %>%
rayimage::render_resized(mag=4) |>
sphere_shade(texture = "desert") |>
add_overlay(generate_polygon_overlay(building_polys, extent = new_e,
heightmap = osm_mat,
linewidth = 6,
resolution_multiply = 50), rescale_original = TRUE) |>
add_overlay(generate_line_overlay(osm_road$osm_lines, extent = new_e,
heightmap = osm_mat,
linewidth = 6,
resolution_multiply = 50), rescale_original = TRUE) |>
plot_3d(osm_mat, water = TRUE, windowsize = 800, watercolor = "dodgerblue",
zscale = 10,
background = "pink")
#Render buildings
render_buildings(building_polys, flat_shading = TRUE,
angle = 30 , heightmap = osm_mat,
material = "white", roof_material = "white",
extent = new_e, roof_height = 3, base_height = 0,
zscale=10)
render_camera(theta=220, phi=22, zoom=0.45, fov=0)
render_snapshot()
}
if(run_documentation()) {
#Zoom in to show roof details and render with render_highquality()
render_camera(fov=110)
render_highquality(camera_location = c(18.22, 0.57, -50.83),
camera_lookat = c(20.88, -2.83, -38.87),
focal_distance = 13,
lightdirection = 45)
}
Run the code above in your browser using DataLab