This vignette demonstrates how to compute building shadow footprints,
shadow height, and surface radiation with sf and
terra objects. The examples use a small synthetic building
layer (sampled in Detroit) so the vignette can run without downloading
external data.
sf building polygon layer with a numeric
Height column.solar_time and time_zone so the
functions estimate solar position from the building-layer centroid.get_shadow_footprint() for shadow polygons.get_shadow_height() with a
terra::SpatRaster.get_radiation() to sample roof and facade radiation
and visualize the returned sf layer.library(gloBFPr)
library(sf)
library(terra)The shadow and radiation functions expect an sf polygon
layer with a numeric height column. Building data returned by
search_3dglobdf() already uses the Height
field, so the same default works here.
The package includes a small sf example dataset with
building footprints, unique IDs, and building heights.
data(globfp_example)
data(globfp_example_dem)
data(globfp_example_canopy_height)
buildings <- globfp_example
# buildings <- buildings[seq_len(min(10, nrow(buildings))), ]
dem <- rast(globfp_example_dem)
canopy_height <- rast(globfp_example_canopy_height)For your own area of interest, first retrieve building footprints
with search_3dglobdf(), then pass the returned polygon
layer to the metric functions.
Start with solar_time and time_zone. The
functions estimate sun azimuth and elevation from the time and the
centroid of the building layer.
If you already have solar geometry from another source, you can use
manual azimuth and elevation values instead.
Azimuth is measured clockwise from north, and elevation is measured
above the horizon.
solar_time <- "2026-06-21 15:00:00"
time_zone <- "America/Detroit"
azimuth <- 135
elevation <- 35get_shadow_footprint() returns an sf
polygon layer containing one shadow footprint for each building.
shadow_footprints <- get_shadow_footprint(
buildings,
solar_time = solar_time,
time_zone = time_zone,
plot = TRUE,
quiet = TRUE
)To compare shadows from multiple times, pass a vector of
solar_time values. When
plot_overlap_gradient = TRUE, the plot uses transparent
gray shadows so overlapping shadow areas appear darker.
solar_times <- data.frame(
sun_id = c("morning", "midday", "afternoon"),
solar_time = c(
"2026-06-21 09:00:00",
"2026-06-21 12:00:00",
"2026-06-21 16:00:00"
)
)
multi_shadow_footprints <- get_shadow_footprint(
buildings,
solar_time = solar_times$solar_time,
time_zone = time_zone,
plot = TRUE,
plot_overlap_gradient = TRUE,
quiet = TRUE
)
unique(multi_shadow_footprints$sun_id)Set overlap_shadow = TRUE when you want a dissolved
cumulative shadow footprint. This is useful for seeing the total ground
area affected by all supplied solar positions.
combined_shadow_footprints <- get_shadow_footprint(
buildings,
solar_time = solar_times$solar_time,
time_zone = time_zone,
overlap_shadow = TRUE,
plot = TRUE,
quiet = TRUE
)Pass a terra::SpatRaster template to compute a gridded
shadow-height surface. If shadow_locations = NULL,
get_shadow_height() creates a template automatically.
template <- rast(
xmin = st_bbox(shadow_footprints)[["xmin"]],
xmax = st_bbox(shadow_footprints)[["xmax"]],
ymin = st_bbox(shadow_footprints)[["ymin"]],
ymax = st_bbox(shadow_footprints)[["ymax"]],
resolution = 4,
crs = st_crs(buildings)$wkt
)
shadow_height_surface <- get_shadow_height(
buildings,
shadow_locations = template,
solar_time = solar_time,
time_zone = time_zone,
quiet = TRUE
)
plot(shadow_height_surface, main = "Shadow height")
plot(st_geometry(buildings), col = NA, border = "black", add = TRUE)Tree canopy can be included as an additional shadow obstacle by passing a canopy height map. If you also provide a DEM, the function compares canopy and building shadows in absolute elevation and returns shadow heights above local ground.
If you do not already have these rasters, the shadow and radiation
functions can retrieve them internally. Canopy height currently supports
datasource_canopy_height = "metachm" or
"ethCHM", and DEM retrieval requires an OpenTopography API
key.
shadow_height_with_downloaded_trees <- get_shadow_height(
buildings,
shadow_locations = template,
solar_time = solar_time,
time_zone = time_zone,
datasource_canopy_height = "metachm",
key = "YOUR_OPENTOPOGRAPHY_API_KEY",
min_tree_height = 2,
quiet = TRUE
)
radiation_with_downloaded_trees <- get_radiation(
buildings,
solar_time = solar_time,
time_zone = time_zone,
solar_normal = 850,
solar_diffuse = 120,
datasource_canopy_height = "metachm",
key = "YOUR_OPENTOPOGRAPHY_API_KEY",
min_tree_height = 2,
canopy_transmissivity = 0.15,
quiet = TRUE
)plot(dem, main = "Sample DEM")
plot(st_geometry(buildings), col = NA, border = "black", add = TRUE)
plot(canopy_height, main = "Canopy height")
plot(st_geometry(buildings), col = NA, border = "black", add = TRUE)The same get_shadow_height() call can now include tree
canopy.
shadow_footprints_with_trees <- get_shadow_footprint(
buildings,
solar_time = solar_time,
time_zone = time_zone,
canopy_height = canopy_height,
dem = dem,
min_tree_height = 1.5,
quiet = TRUE,
plot = TRUE
)shadow_height_with_trees <- get_shadow_height(
buildings,
shadow_locations = template,
solar_time = solar_time,
time_zone = time_zone,
canopy_height = canopy_height,
dem = dem,
min_tree_height = 2,
quiet = TRUE
)
plot(shadow_height_with_trees, main = "Building and canopy shadow height")
plot(st_geometry(buildings), col = NA, border = "black", add = TRUE)get_radiation() estimates direct, diffuse, and total
radiation on building roofs and facades. When no custom grid is
supplied, the function automatically samples the full 3D building
surface — both horizontal roof points and vertical facade strips.
solar_normal is direct normal irradiance and
solar_diffuse is diffuse horizontal irradiance, each as a
numeric vector with one value per solar_time.
Use plot = TRUE for the default 2D radiation map. Set
plot_3d = TRUE only when you also want the isometric
direct, diffuse, and total radiation view.
grid_res controls the horizontal and vertical sampling
density of the surface grid. Smaller values create more strips per wall
and a smoother color gradient:
grid_res |
Approx. strips per 10 m wall |
|---|---|
| 8 (default) | 1–2 |
| 4 | 2–3 |
| 2 | ~5 |
radiation <- get_radiation(
buildings,
solar_time = solar_time,
time_zone = time_zone,
solar_normal = 850,
solar_diffuse = 120,
grid_res = 4,
plot_3d = TRUE
)
head(st_drop_geometry(radiation))Single vs. multiple solar times. At a single time step, direct radiation on horizontal rooftops is binary — each roof is either fully lit (
solar_normal × sin(elevation)) or in shadow (0). Spatial variation in the direct panel therefore reflects the shadow pattern at that moment. For a continuous gradient across the scene, accumulate several hours:solar_day <- format( seq(as.POSIXct("2026-06-21 07:00", tz = time_zone), as.POSIXct("2026-06-21 19:00", tz = time_zone), by = "hour"), "%Y-%m-%d %H:%M:%S" ) radiation_day <- get_radiation( buildings, solar_time = solar_day, time_zone = time_zone, solar_normal = rep(850, length(solar_day)), solar_diffuse = rep(120, length(solar_day)), grid_res = 4, plot_3d = TRUE )
To include tree shade, pass a canopy height map and optionally a DEM.
canopy_transmissivity controls how much direct radiation
penetrates the canopy; 0 is fully opaque and 1
has no effect.
radiation_with_trees <- get_radiation(
buildings,
solar_time = solar_time,
time_zone = time_zone,
solar_normal = 850,
solar_diffuse = 120,
canopy_height = canopy_height,
dem = dem,
min_tree_height = 2,
canopy_transmissivity = 0.15,
grid_res = 4,
plot_3d = TRUE,
quiet = TRUE
)Compare mean roof radiation with and without tree canopy:
roof_building_only <- aggregate(
total ~ building_id,
data = st_drop_geometry(radiation[radiation$surface == "roof", ]),
FUN = mean
)
roof_with_trees <- aggregate(
total ~ building_id,
data = st_drop_geometry(radiation_with_trees[radiation_with_trees$surface == "roof", ]),
FUN = mean
)
roof_compare <- merge(
roof_building_only,
roof_with_trees,
by = "building_id",
suffixes = c("_building_only", "_with_trees")
)
head(roof_compare)Because both radiation objects share the same surface grid (same
buildings, same grid_res, same solar_time),
rows correspond one-to-one and the difference can be computed directly.
Negative values mean tree shade reduced radiation at that sample
point.
radiation_diff <- radiation
radiation_diff$direct_diff <- radiation_with_trees$direct - radiation$direct
radiation_diff$total_diff <- radiation_with_trees$total - radiation$total
# Summary by surface type
aggregate(
cbind(direct_diff, total_diff) ~ surface,
data = st_drop_geometry(radiation_diff),
FUN = function(x) round(mean(x), 1)
)Map the total radiation reduction. The palette is centered at 0 using symmetric breaks so blue always means reduction and red means increase, regardless of the asymmetric range.
diff_pal <- hcl.colors(100, "Blue-Red 3")
max_abs <- max(abs(radiation_diff$total_diff), na.rm = TRUE)
diff_breaks <- cut(radiation_diff$total_diff,
breaks = seq(-max_abs, max_abs, length.out = 101),
include.lowest = TRUE, labels = FALSE)
plot(st_geometry(buildings), col = "grey95", border = "grey45",
main = "Total radiation change due to tree canopy (W/m²)")
plot(st_geometry(radiation_diff), pch = 16, cex = 0.45,
col = diff_pal[diff_breaks], add = TRUE)
legend("topright",
legend = c(paste0("≤ ", -round(max_abs, 0)), "0",
paste0("≥ +", round(max_abs, 0))),
pch = 16,
col = c(diff_pal[1], diff_pal[50], diff_pal[100]),
bty = "n")The returned sf object is a point layer, so standard
sf map plots work directly. Roof and facade points are most
readable when shown separately.
roof_radiation <- radiation[radiation$surface == "roof", ]
facade_radiation <- radiation[radiation$surface == "facade", ]
plot(roof_radiation["total"], pch = 16, cex = 0.8, key.pos = 4)
plot(st_geometry(buildings), col = NA, border = "grey30", add = TRUE)
plot(facade_radiation["total"], pch = 16, cex = 0.45, key.pos = 4)
plot(st_geometry(buildings), col = NA, border = "grey30", add = TRUE)Passing multiple solar times accumulates radiation across the day and produces more spatial variation in the direct component — surfaces that spend more time in sunlight receive proportionally higher totals.
solar_day <- c(
"2026-06-21 08:00:00",
"2026-06-21 11:00:00",
"2026-06-21 14:00:00",
"2026-06-21 17:00:00"
)
radiation_day <- get_radiation(
buildings,
solar_time = solar_day,
time_zone = time_zone,
solar_normal = c(500, 850, 900, 650),
solar_diffuse = c(180, 120, 110, 150),
grid_res = 4,
plot_3d = TRUE,
quiet = TRUE
)roof_day <- radiation_day[radiation_day$surface == "roof", ]
facade_day <- radiation_day[radiation_day$surface == "facade", ]
roof_day_breaks <- cut(roof_day$total, breaks = 100,
include.lowest = TRUE, labels = FALSE)
plot(st_geometry(buildings), col = "grey95", border = "grey45")
plot(st_geometry(facade_day), pch = 16, cex = 0.25,
col = "#2563EB55", add = TRUE)
plot(st_geometry(roof_day), pch = 16, cex = 0.75,
col = roof_cols[roof_day_breaks], add = TRUE)
legend("topright",
legend = c("lower roof total", "higher roof total", "facade samples"),
pch = 16,
col = c(roof_cols[1], roof_cols[100], "#2563EB55"),
bty = "n")Setting ground = TRUE adds a regular grid of sample
points at street/terrain level to the radiation object. Ground points
lie outside all building footprints and have an upward-facing normal
(nz = 1), so they receive:
solar_normal × sin(elevation) when not in a building or
canopy shadow, and 0 (or attenuated by
canopy_transmissivity) when shaded.solar_diffuse × SVF, where SVF is the Sky View Factor at
ground level. Adjacent buildings reduce SVF, lowering diffuse radiation
near walls and in narrow alleys.ground_res sets the ground grid spacing independently
from grid_res (the building surface spacing). A finer
ground_res produces a denser radiation map but takes longer
to compute.
radiation_ground <- get_radiation(
buildings,
solar_time = c("2026-06-21 08:00:00", "2026-06-21 11:00:00",
"2026-06-21 14:00:00", "2026-06-21 17:00:00"),
time_zone = time_zone,
solar_normal = c(500, 850, 900, 650),
solar_diffuse = c(180, 120, 110, 150),
grid_res = 5,
ground = TRUE,
ground_res = 5,
plot = TRUE
)With ground = TRUE, plot = TRUE draws three
2D maps: ground surface, facade, and roof total radiation. The maps
share one W/m² color bar.
Separate the ground layer for further analysis:
ground_rad <- radiation_ground[radiation_ground$surface == "ground", ]
summary(st_drop_geometry(ground_rad[, c("svf", "direct", "diffuse", "total")]))Passing a canopy height map attenuates shaded samples using
canopy_transmissivity. With plot = TRUE,
get_radiation() also plots the canopy impact as
canopy - no_canopy total radiation.
radiation_ground_trees <- get_radiation(
buildings,
solar_time = c("2026-06-21 08:00:00", "2026-06-21 11:00:00",
"2026-06-21 14:00:00", "2026-06-21 17:00:00"),
time_zone = time_zone,
solar_normal = c(500, 850, 900, 650),
solar_diffuse = c(180, 120, 110, 150),
canopy_height = canopy_height,
dem = dem,
min_tree_height = 2,
canopy_transmissivity = 0.15,
grid_res = 5,
ground = TRUE,
ground_res = 5,
quiet = TRUE,
plot = TRUE
)ground_trees <- radiation_ground_trees[
radiation_ground_trees$surface == "ground", ]
# Mean radiation reduction at ground level due to tree canopy
mean(ground_trees$total - ground_rad$total, na.rm = TRUE)Combining ground and roof radiation in one call lets you compare solar exposure across all horizontal surfaces at once:
roof_ground <- radiation_ground[radiation_ground$surface %in% c("roof", "ground"), ]
aggregate(
cbind(direct, diffuse, total) ~ surface,
data = st_drop_geometry(roof_ground),
FUN = function(x) round(mean(x, na.rm = TRUE), 1)
)