Constructing Spatial Covariates with spatcovar

Introduction

In applied spatial research, researchers frequently start with polygon units (districts, municipalities, administrative boundaries, grid cells, or ecological zones) and need to attach spatial features derived from other geospatial layers:

Polygon geometric area
Distances to central places, coastlines, or infrastructure
Point counts (e.g., facilities, events, schools)
Linework lengths within each polygon (e.g., roads, rivers, borders)
Overlap area or coverage share with other polygon layers
Zonal summary statistics from raster grids (e.g., elevation, climate, population density)

While packages such as sf, terra, and exactextractr provide the low-level spatial primitives for these operations, combining them across multiple layers often requires substantial boilerplate code to manage coordinate reference systems, geometric validation, units, missing-value semantics, and row preservation.

spatcovar provides a consistent, pipeable interface for constructing these spatial covariates directly on sf polygon data frames.

Synthetic Example Data

spatcovar provides lightweight synthetic fixtures in projected coordinates (EPSG:32632, UTM Zone 32N) so you can explore all operations without external data dependencies:

library(spatcovar)

regions <- example_polygons()
sites   <- example_points()
routes  <- example_lines()
zones   <- example_grid()
rst     <- example_raster()

# Inspect the target polygon layer
regions
#> Simple feature collection with 6 features and 1 field
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 4e+05 ymin: 5400000 xmax: 430000 ymax: 5420000
#> Projected CRS: WGS 84 / UTM zone 32N
#>       name                       geometry
#> 1 region_a POLYGON ((4e+05 5400000, 41...
#> 2 region_b POLYGON ((410000 5400000, 4...
#> 3 region_c POLYGON ((420000 5400000, 4...
#> 4 region_d POLYGON ((4e+05 5410000, 41...
#> 5 region_e POLYGON ((412000 5410000, 4...
#> 6 region_f POLYGON ((420000 5410000, 4...

Core Covariate Pipeline

All spatcovar functions take the target sf polygons as the first argument, calculate the requested metric, and append the result as a new column while strictly preserving the original row count, row order, and geometry.

When name = NULL (the default), output column names are generated dynamically from the requested unit or measure (e.g., area_km2, dist_km, length_km, overlap_share).

covariates <- regions |>
  spat_area(unit = "km2") |>
  spat_count(sites, name = "n_sites") |>
  spat_length(routes, unit = "km") |>
  spat_distance(sites, method = "minimum", unit = "km") |>
  spat_overlap(zones, measure = "share") |>
  spat_raster(rst, stats = c("mean", "max"), name = "elevation")
#> Warning: Source polygons overlap each other with positive area. Overlap share
#> may exceed 1.0 for some target polygons.

# View the constructed covariate table
sf::st_drop_geometry(covariates)
#>       name area_km2 n_sites length_km dist_km overlap_share elevation_mean
#> 1 region_a      100       3  21.18034       0     0.2500000      221.60999
#> 2 region_b      100       2  20.00000       0     0.8000000      227.49998
#> 3 region_c      100       3  10.00000       0     0.6400000      233.39000
#> 4 region_d      110       2  11.18034       0     0.4218182       85.91727
#> 5 region_e       80       1  10.00000       0     0.9125000       94.07500
#> 6 region_f      100       3   0.00000       0     0.8000000      108.64603
#>   elevation_max
#> 1           288
#> 2           293
#> 3           299
#> 4           148
#> 5           153
#> 6           159

Methodological Details & Semantics

1. Polygon Area (`spat_area`)

spat_area() calculates the area of each polygon. When input polygons use a geographic coordinate reference system (longitude/latitude, e.g. WGS84), spat_area() calculates geodesic areas on the sphere/ellipsoid using sf (via s2). For projected coordinate systems, planar area is calculated.

Supported units include "m2", "km2", "ha", and "mi2".

regions |>
  spat_area(unit = "ha") |>
  subset(select = c(name, area_ha))
#> Simple feature collection with 6 features and 2 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 4e+05 ymin: 5400000 xmax: 430000 ymax: 5420000
#> Projected CRS: WGS 84 / UTM zone 32N
#>       name area_ha                       geometry
#> 1 region_a   10000 POLYGON ((4e+05 5400000, 41...
#> 2 region_b   10000 POLYGON ((410000 5400000, 4...
#> 3 region_c   10000 POLYGON ((420000 5400000, 4...
#> 4 region_d   11000 POLYGON ((4e+05 5410000, 41...
#> 5 region_e    8000 POLYGON ((412000 5410000, 4...
#> 6 region_f   10000 POLYGON ((420000 5410000, 4...

2. Distance to Reference Features (`spat_distance`)

spat_distance() calculates the distance from each target polygon to reference features. Three measurement methods are supported:

"minimum" (default): Shortest geometry-to-geometry distance between the polygon boundary/interior and the nearest feature in y.
"centroid": Euclidean or geodesic distance from the target polygon centroid to the nearest feature in y.
"point_on_surface": Distance from a point guaranteed to lie inside the polygon (useful for irregular or concave polygons where the centroid may fall outside).

For multi-feature reference layers, spat_distance() uses spatial indexing (sf::st_nearest_feature()) to locate nearest candidate features efficiently without calculating the full \(N \times M\) distance matrix.

regions |>
  spat_distance(sites, method = "minimum", unit = "km") |>
  spat_distance(sites, method = "centroid", unit = "km", name = "dist_cent_km") |>
  subset(select = c(name, dist_km, dist_cent_km))
#> Simple feature collection with 6 features and 3 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 4e+05 ymin: 5400000 xmax: 430000 ymax: 5420000
#> Projected CRS: WGS 84 / UTM zone 32N
#>       name dist_km dist_cent_km                       geometry
#> 1 region_a       0     0.000000 POLYGON ((4e+05 5400000, 41...
#> 2 region_b       0     1.000000 POLYGON ((410000 5400000, 4...
#> 3 region_c       0     0.000000 POLYGON ((420000 5400000, 4...
#> 4 region_d       0     2.519711 POLYGON ((4e+05 5410000, 41...
#> 5 region_e       0     1.000000 POLYGON ((412000 5410000, 4...
#> 6 region_f       0     0.000000 POLYGON ((420000 5410000, 4...

3. Intersecting Feature Counts (`spat_count`)

spat_count() counts the number of source features that spatially intersect each polygon.

Key semantics: - Features touching the polygon boundary are counted. - A MULTIPOINT, MULTILINESTRING, or MULTIPOLYGON record represents one source feature. To count individual constituent points, use sf::st_cast(sites, "POINT") before passing to spat_count(). - Polygons with no intersecting features receive an explicit count of 0L.

regions |>
  spat_count(sites, name = "site_count") |>
  subset(select = c(name, site_count))
#> Simple feature collection with 6 features and 2 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 4e+05 ymin: 5400000 xmax: 430000 ymax: 5420000
#> Projected CRS: WGS 84 / UTM zone 32N
#>       name site_count                       geometry
#> 1 region_a          3 POLYGON ((4e+05 5400000, 41...
#> 2 region_b          2 POLYGON ((410000 5400000, 4...
#> 3 region_c          3 POLYGON ((420000 5400000, 4...
#> 4 region_d          2 POLYGON ((4e+05 5410000, 41...
#> 5 region_e          1 POLYGON ((412000 5410000, 4...
#> 6 region_f          3 POLYGON ((420000 5410000, 4...

4. Linework Length inside Polygons (`spat_length`)

spat_length() clips line features to each polygon boundary and sums the total length of the intersecting segments. It explicitly accounts for the native linear units of the projected CRS (e.g. feet vs. metres).

CRS Requirement: spat_length() requires a projected (planar) coordinate reference system. If your data uses geographic coordinates (degrees), supply a projected CRS via crs = ...:

regions |>
  spat_length(routes, unit = "km") |>
  subset(select = c(name, length_km))
#> Simple feature collection with 6 features and 2 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 4e+05 ymin: 5400000 xmax: 430000 ymax: 5420000
#> Projected CRS: WGS 84 / UTM zone 32N
#>       name length_km                       geometry
#> 1 region_a  21.18034 POLYGON ((4e+05 5400000, 41...
#> 2 region_b  20.00000 POLYGON ((410000 5400000, 4...
#> 3 region_c  10.00000 POLYGON ((420000 5400000, 4...
#> 4 region_d  11.18034 POLYGON ((4e+05 5410000, 41...
#> 5 region_e  10.00000 POLYGON ((412000 5410000, 4...
#> 6 region_f   0.00000 POLYGON ((420000 5410000, 4...

5. Polygon Overlap Measures (`spat_overlap`)

spat_overlap() calculates relationships between target polygons and another polygon layer y:

measure = "area": Total area of intersection in the requested unit ("m2", "km2", "ha", "mi2").
measure = "share": Fraction of the target polygon’s area covered by source polygons. The share is typically between 0.0 and 1.0 when source polygons form a non-overlapping coverage, but may exceed 1.0 when source polygons overlap one another with positive area.
measure = "count": Number of source polygons intersecting the target polygon.

When source polygons in y have positive-area duplicate coverage among themselves, spat_overlap(..., measure = "share") issues an informative warning that the sum of overlapping parts may exceed \(1.0\). Source polygons that merely touch along edges or at vertices do not trigger the warning.

regions |>
  spat_overlap(zones, measure = "area", unit = "km2") |>
  spat_overlap(zones, measure = "share") |>
  subset(select = c(name, overlap_km2, overlap_share))
#> Warning: Source polygons overlap each other with positive area. Overlap share
#> may exceed 1.0 for some target polygons.
#> Simple feature collection with 6 features and 3 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 4e+05 ymin: 5400000 xmax: 430000 ymax: 5420000
#> Projected CRS: WGS 84 / UTM zone 32N
#>       name overlap_km2 overlap_share                       geometry
#> 1 region_a        25.0     0.2500000 POLYGON ((4e+05 5400000, 41...
#> 2 region_b        80.0     0.8000000 POLYGON ((410000 5400000, 4...
#> 3 region_c        64.0     0.6400000 POLYGON ((420000 5400000, 4...
#> 4 region_d        46.4     0.4218182 POLYGON ((4e+05 5410000, 41...
#> 5 region_e        73.0     0.9125000 POLYGON ((412000 5410000, 4...
#> 6 region_f        80.0     0.8000000 POLYGON ((420000 5410000, 4...

6. Raster Zonal Statistics (`spat_raster`)

spat_raster() computes zonal statistics across single-layer raster grids using the exact coverage-fraction weighting engine provided by exactextractr.

Key properties: - Supported statistics: "mean", "median", "min", "max", "sum", "count", "stdev". - mean, sum, median, and stdev are weighted by the fraction of each cell covered by the polygon. - min and max are unweighted extrema over intersected cells. - count is the sum of coverage fractions of cells with non-NA values (effective covered cells) and can be fractional. - Valid zero values in the raster are preserved as valid observed data. - If a target polygon has no valid non-NA raster cells contributing to the extraction (outside raster extent or all covered cells are NA), all requested statistics evaluate to NA. - If the polygon CRS and raster CRS differ, polygon copies are automatically reprojected to the raster’s CRS during extraction.

regions |>
  spat_raster(rst, stats = c("mean", "median", "min", "max"), name = "topo") |>
  subset(select = c(name, topo_mean, topo_median, topo_min, topo_max))
#> Simple feature collection with 6 features and 5 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 4e+05 ymin: 5400000 xmax: 430000 ymax: 5420000
#> Projected CRS: WGS 84 / UTM zone 32N
#>       name topo_mean topo_median topo_min topo_max
#> 1 region_a 221.60999   223.66507      162      288
#> 2 region_b 227.49998   229.54295      168      293
#> 3 region_c 233.39000   235.41417      173      299
#> 4 region_d  85.91727    85.02845       22      148
#> 5 region_e  94.07500    91.84762       29      153
#> 6 region_f 108.64603   115.20147       33      159
#>                         geometry
#> 1 POLYGON ((4e+05 5400000, 41...
#> 2 POLYGON ((410000 5400000, 4...
#> 3 POLYGON ((420000 5400000, 4...
#> 4 POLYGON ((4e+05 5410000, 41...
#> 5 POLYGON ((412000 5410000, 4...
#> 6 POLYGON ((420000 5410000, 4...

Operation Diagnostics (`spat_diagnostics`)

Every spat_* function accepts diagnostics = TRUE. When enabled, lightweight metadata describing the operation is stored in the spat_diagnostics attribute:

result <- regions |>
  spat_distance(sites, unit = "km", diagnostics = TRUE)

spat_diagnostics(result)
#> --- spatcovar diagnostics ---
#> Operation: spat_distance
#> Target polygons: 6
#> Source features: 15
#> NA results: 0
#> Method: minimum
#> CRS used: EPSG:32632
#> Unit: km

Summary of Guarantees

Row Preservation: The output always has identical row count and row order as the input target polygons.
Original Geometry: Returned geometries remain identical to the input geometries even when internal geometry repair was necessary.
Collision Safety: Overwriting existing columns requires overwrite = TRUE.
Predictable Types: Counts return integers; metric quantities return standard base numerics in explicit units.
Safe Missing-CRS Policy: Operations that require knowing coordinate meaning or units error explicitly when CRS information is missing.
Standardised Missing-Data Semantics: Zero-support raster extractions evaluate to NA across all requested statistics.

Welcome to ClientVPS Mirrors

Constructing Spatial Covariates with spatcovar

Introduction

Synthetic Example Data

Core Covariate Pipeline

Methodological Details & Semantics

1. Polygon Area (spat_area)

2. Distance to Reference Features (spat_distance)

3. Intersecting Feature Counts (spat_count)

4. Linework Length inside Polygons (spat_length)

5. Polygon Overlap Measures (spat_overlap)

6. Raster Zonal Statistics (spat_raster)

Operation Diagnostics (spat_diagnostics)