---
title: "Block-Level Analysis with gloBFPr"
author: "Xiaohao Yang"
vignette: >
  %\VignetteIndexEntry{Block-Level Analysis with gloBFPr}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

# 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.

```{r setup, message=FALSE}
library(gloBFPr)
library(sf)
library(dplyr)
```

The package includes a small building footprint layer with heights. We use it
throughout this vignette.

```{r}
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.

```{r eval=FALSE}
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.

```{r eval=FALSE}
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.

```{r eval=FALSE}
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.

```{r eval=FALSE}
# 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.

```{r eval=FALSE}
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.

```{r eval=FALSE}
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.

```{r eval=FALSE}
# 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`.

```{r eval=FALSE}
# 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`).

```{r eval=FALSE}
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`).

```{r eval=FALSE}
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.

```{r eval=FALSE}
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.

```{r eval=FALSE}
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:

```{r eval=FALSE}
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.

```{r eval=FALSE}
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.

```{r eval=FALSE}
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.

```{r eval=FALSE}
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)
```