Urban blocks — the land parcels enclosed by streets — are a natural unit for summarising building-level data into neighbourhood-scale indicators. This vignette walks through the two functions that support block-level analysis in gloBFPr:
generate_block() delineates block polygons from a road
network and assigns each building to its enclosing block.aggregate_block() rolls building-level metrics up to
the block level.library(gloBFPr)
library(sf)
library(dplyr)The package includes a small building footprint layer with heights. We use it throughout this vignette.
data(globfp_example)
buildings <- globfp_exampleFor your own area of interest, retrieve footprints first with
search_3dglobdf().
generate_block() uses a two-stage approach. In the first
stage it polygonizes the road network directly: each closed loop of
streets becomes a block polygon and buildings are assigned by
centroid-within spatial join. In the second stage, any buildings not
covered by the first stage (due to network gaps or dead-end pockets in
OSM data) are handled via a raster-based fallback: the network is
rasterized as a barrier, connected non-road patches are identified, and
the patches that contain at least one building are converted back to
polygons.
By default the road network is fetched from Overture Maps for the
bounding box of the input buildings (via a DuckDB parquet query, which
requires the duckdb and DBI packages). Set
network_source = "osm" to fetch from OpenStreetMap via the
Overpass API instead, which requires the osmdata
package.
block_result <- generate_block(buildings, quiet = FALSE)The return value is a named list. $blocks is an
sf polygon object with one row per block and a
block_id column. $buildings is the input layer
with block_id appended.
blocks <- block_result$blocks
buildings <- block_result$buildings
# How many blocks were found?
nrow(blocks)
# How many buildings were assigned?
sum(!is.na(buildings$block_id))Supply a pre-downloaded network to skip any remote fetch, which is useful for offline workflows or when you have a pre-processed network.
library(osmdata)
net <- opq(bbox = sf::st_bbox(buildings)) |>
add_osm_feature("highway",
value = c("motorway", "trunk", "primary", "secondary",
"tertiary", "residential", "unclassified")) |>
osmdata_sf()
net_lines <- net$osm_lines
block_result <- generate_block(buildings, network = net_lines, quiet = FALSE)
plot(block_result$blocks)Motorways and trunk roads in OSM are commonly mapped as dual
carriageways: two parallel lines for opposite directions of travel.
Without correction, polygonization treats the narrow strip between them
as a valid block, producing spurious slivers.
generate_block() removes these before polygonizing by
comparing each line in the target highway classes against its spatial
neighbours: if a shorter line overlaps the buffer of the longest
neighbour by more than dc_overlap_threshold, it is
dropped.
The default settings handle most urban road networks well. You can widen or narrow the affected road classes and adjust the overlap threshold if your data requires it.
# Include primary roads in dual-carriageway simplification
block_result <- generate_block(
buildings,
dc_highway_types = c("motorway", "trunk", "primary"),
dc_overlap_threshold = 0.7,
quiet = FALSE
)To disable dual carriageway simplification entirely, pass an empty character vector.
block_result <- generate_block(
buildings,
dc_highway_types = character(0),
quiet = FALSE
)min_block_area (in m²) controls the minimum size of a
valid block. In the polygonize stage, enclosures smaller than this
threshold are dropped before the building join. In the raster fallback
stage, small patches are merged into their nearest larger neighbour. The
default of 500 m² filters out residual slivers while keeping compact
urban blocks. Increase the value in areas with very fine-grained street
grids; decrease it for historic centres with small medieval blocks.
block_result <- generate_block(
buildings,
min_block_area = 1000,
quiet = FALSE
)Once buildings have been assigned to blocks,
aggregate_block() summarises any numeric or logical columns
in $buildings at the block level. It expects the list
returned by generate_block() directly.
Before calling aggregate_block() it is useful to compute
building-level metrics first so they are available for aggregation.
# Compute building metrics, then generate blocks
buildings_with_metrics <- buildings |>
get_morphology(quiet = TRUE)
block_result <- generate_block(buildings_with_metrics, quiet = FALSE)
# Aggregate to block level
block_metrics <- aggregate_block(block_result, quiet = FALSE)
block_metrics |>
st_drop_geometry() |>
select(block_id, n_buildings, coverage_ratio, g_area, vol) |>
head()aggregate_block() applies sensible defaults
automatically:
g_area, pmeter,
v_surf, t_surf, vol,
obb_vol (quantities that are additive across buildings in a
block).rec, fra,
cnv, and elongation ratios.n_buildings (count) and
coverage_ratio (total building footprint area divided by
block polygon area, computed in UTM for accuracy).Note that per-building residential classification columns
(res, res_pct, and related GHS sub-columns
from get_residential()) are excluded from the default block
aggregation. Use residential = TRUE instead to derive a
geometrically accurate residential proportion directly from GHS rasters
at block level (see below).
Override the default for any column by passing a named list to
.fns. The names must match column names in
$buildings.
# Use maximum volume instead of sum; use median height instead of mean
block_metrics_custom <- aggregate_block(
block_result,
.fns = list(vol = max, Height = median),
quiet = FALSE
)Set population = TRUE to fetch GHSL population counts
directly at block level. Population is extracted by overlaying each
block polygon against the GHS population raster and summing cell values
weighted by the overlap area, giving the total resident count per block.
Use population_year to select the GHSL epoch (1975–2030,
default 2025).
block_metrics <- aggregate_block(
block_result,
population = TRUE,
population_year = 2020,
quiet = FALSE
)
block_metrics |>
st_drop_geometry() |>
select(block_id, n_buildings, pop_total) |>
arrange(desc(pop_total)) |>
head()This approach queries the raster directly at block level, which avoids the error introduced by first allocating population to individual buildings and then summing back up.
Set residential = TRUE to compute res_prop
— the fraction of total built-up surface area within each block that is
classified as residential by the GHS built-up surface layer. Use
residential_year to select the GHS epoch (default
2020).
block_metrics <- aggregate_block(
block_result,
residential = TRUE,
residential_year = 2020,
quiet = FALSE
)
block_metrics |>
st_drop_geometry() |>
select(block_id, n_buildings, res_prop) |>
arrange(desc(res_prop)) |>
head()res_prop is computed as the area-weighted sum of the
residential built-up surface raster divided by the area-weighted sum of
the total built-up surface raster over each block polygon.
Because aggregate_block() returns an sf
polygon layer, block-level metrics can be mapped directly with
ggplot2::geom_sf(). Choropleth maps are useful for checking
whether high-density, high-volume, or high-population blocks cluster in
expected parts of the study area.
library(ggplot2)
ggplot(block_metrics) +
geom_sf(aes(fill = coverage_ratio), color = "white", linewidth = 0.15) +
scale_fill_viridis_c(
option = "magma",
labels = function(x) paste0(round(100 * x), "%"),
na.value = "grey90"
) +
coord_sf(datum = NA) +
labs(
fill = "Coverage",
title = "Building coverage ratio by block"
) +
theme_minimal()For count or volume metrics, a transformed scale usually makes compact urban blocks easier to compare because a few very large blocks can otherwise dominate the color range.
ggplot(block_metrics) +
geom_sf(aes(fill = n_buildings), color = "grey85", linewidth = 0.1) +
scale_fill_viridis_c(
option = "plasma",
trans = "sqrt",
na.value = "grey90"
) +
coord_sf(datum = NA) +
labs(
fill = "Buildings",
title = "Number of buildings per block"
) +
theme_minimal()
ggplot(block_metrics) +
geom_sf(aes(fill = vol), color = NA) +
scale_fill_viridis_c(
option = "cividis",
trans = "log10",
labels = function(x) format(x, scientific = TRUE),
na.value = "grey90"
) +
coord_sf(datum = NA) +
labs(
fill = "Volume",
title = "Total building volume by block"
) +
theme_minimal()When population or residential proportion has been requested, the same pattern can be used for demographic or land-use indicators:
ggplot(block_metrics) +
geom_sf(aes(fill = pop_total), color = "white", linewidth = 0.1) +
scale_fill_viridis_c(
option = "inferno",
trans = "sqrt",
na.value = "grey90"
) +
coord_sf(datum = NA) +
labs(
fill = "Population",
title = "Estimated population by block"
) +
theme_minimal()
ggplot(block_metrics) +
geom_sf(aes(fill = res_prop), color = "white", linewidth = 0.1) +
scale_fill_viridis_c(
option = "viridis",
limits = c(0, 1),
labels = function(x) paste0(round(100 * x), "%"),
na.value = "grey90"
) +
coord_sf(datum = NA) +
labs(
fill = "Residential",
title = "Residential built-up surface proportion by block"
) +
theme_minimal()To show how individual buildings relate to the block summary, join the block-level metric back to the building layer and draw buildings over the block polygons.
buildings_for_map <- block_result$buildings |>
left_join(
block_metrics |>
st_drop_geometry() |>
select(block_id, block_coverage_ratio = coverage_ratio),
by = "block_id"
)
ggplot() +
geom_sf(data = block_metrics, fill = "grey95", color = "white",
linewidth = 0.15) +
geom_sf(data = buildings_for_map,
aes(fill = block_coverage_ratio),
color = "grey25", linewidth = 0.05) +
scale_fill_viridis_c(
option = "magma",
labels = function(x) paste0(round(100 * x), "%"),
na.value = "grey90"
) +
coord_sf(datum = NA) +
labs(
fill = "Block coverage",
title = "Buildings colored by their block-level coverage ratio"
) +
theme_minimal()For exploratory analysis, it is often useful to scan several block-level indicators side by side before deciding which metric deserves a dedicated map. The pattern below reshapes selected metrics into a long table and draws a faceted choropleth with one shared layout. This keeps the block geometry, color scale, and map framing consistent across metrics.
metric_labels <- c(
coverage_ratio = "Coverage ratio",
n_buildings = "Building count",
vol = "Building volume"
)
metrics_long <- bind_rows(lapply(names(metric_labels), function(metric_col) {
out <- block_metrics[, c("block_id", metric_col, "geometry")]
names(out)[names(out) == metric_col] <- "value"
out$metric <- metric_labels[[metric_col]]
out
}))
ggplot(metrics_long) +
geom_sf(aes(fill = value), color = "white", linewidth = 0.1) +
scale_fill_viridis_c(
option = "mako",
trans = "sqrt",
na.value = "grey90"
) +
facet_wrap(vars(metric), nrow = 1) +
coord_sf(datum = NA) +
labs(fill = "Value", title = "Block-level metric comparison") +
theme_minimal() +
theme(
panel.grid = element_blank(),
strip.text = element_text(face = "bold")
)The typical workflow is to fetch buildings, compute building-level metrics, and then generate and summarize blocks in sequence.
bbox <- c(-83.065644, 42.333792, -83.045217, 42.346988)
# 1. Retrieve building footprints
buildings <- search_3dglobdf(bbox = bbox, out_type = "poly", quiet = TRUE)
# 2. Compute building-level metrics
buildings <- buildings |>
get_morphology(quiet = TRUE)
# 3. Generate blocks
block_result <- generate_block(buildings, quiet = FALSE)
# 4. Aggregate metrics to block level, including population and residential
block_metrics <- aggregate_block(
block_result,
population = TRUE,
population_year = 2025,
residential = TRUE,
residential_year = 2025,
quiet = FALSE
)
# 5. Inspect
block_metrics |>
st_drop_geometry() |>
select(block_id, n_buildings, coverage_ratio, g_area, vol,
pop_total, res_prop) |>
arrange(desc(n_buildings)) |>
head(10)