# Block-Level Analysis with gloBFPr

#### Xiaohao Yang

# Block-Level Analysis

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_example

For your own area of interest, retrieve footprints first with `search_3dglobdf()`.

## 1 Generating blocks

`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)

### Dual carriageway simplification

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
    )

### Controlling block size

`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
    )

## 2 Aggregating metrics to block level

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:

- Summed — `g_area`, `pmeter`, `v_surf`, `t_surf`, `vol`, `obb_vol` (quantities that are additive across buildings in a block).
- Averaged — all other numeric columns, including shape indices such as `rec`, `fra`, `cnv`, and elongation ratios.
- Always added — `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).

### Custom aggregation functions

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
    )

### Block-level population

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.

### Block-level residential proportion

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.

## 3 Visualizing block-level metrics

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()

### Comparing several block metrics in one figure

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")
      )

## 4 A complete block analysis pipeline

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)