---
title: "Streaming spatial operations"
author: "Gilles Colling"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Streaming spatial operations}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
has_sf <- requireNamespace("sf", quietly = TRUE) &&
  requireNamespace("s2", quietly = TRUE)
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval = has_sf,
  fig.width = 7,
  fig.height = 4.2,
  out.width = "100%",
  dpi = 96
)
```

This vignette is the hands-on companion to
[Out-of-core GIS with vectra](spatial.html). That article lays out the cost
model; this one runs the streaming vector verbs end to end on a real layer, so
every block below is code you can execute and check. It needs the optional
[sf](https://r-spatial.github.io/sf/) package, which the spatial verbs use to
take layers and coordinate reference systems as arguments, and which runs the
operations vectra has no native path for.

Installation is covered in [Get started](quickstart.html).

```{r libs, message = FALSE}
library(vectra)
library(sf)
```

## The problem

A desktop GIS holds one layer in memory, runs a tool, and writes the result.
That model breaks the moment the working layer is bigger than RAM: a national
occurrence set, a continental road network, every parcel in a country. The
geometry alone runs to tens of gigabytes before any operation touches it.

vectra keeps the same toolbox but changes what stays resident. A query is pulled
through the C engine in fixed-size batches; each spatial step works on the batch
in front of it and spills the transformed batch back to disk as a fresh lazy
node. Peak memory is one batch plus whatever small layer the step compares
against, so a billion-row layer flows past a fixed budget. The verbs in this
vignette are that streaming envelope: `spatial_map()`, `spatial_filter()`,
`spatial_clip()`, `spatial_join()`, `spatial_dissolve()`, `spatial_overlay()`,
and `rasterize()`.

## How geometry travels

vectra has no geometry type. A geometry rides through the engine as hex-encoded
WKB in an ordinary string column, and the coordinate reference system is carried
on the returned node rather than written into the `.vtr` file. The geometry
algorithms come from GEOS. For the operations vectra recognises (the common
predicates, clipping, buffering, measures and the like on planar data) it calls
the GEOS C API directly, through the
[libgeos](https://paleolimbot.github.io/libgeos/) package it links against, and
works straight off the WKB column in C. Everything else, such as an arbitrary
function passed to `spatial_map()` or geographic data with spherical geometry
switched on, runs through sf batch by batch.

On that sf path, one streamed step is a loop. Pull a batch, decode its WKB
column into an `sf` object, run the operation, encode the result back into a
WKB string column, and append it to a run-file. When the loop finishes, the run-files become a single
lazy node you can carry on querying. The memory a run holds is

```
peak  =  one batch  +  the resident comparison layer
```

independent of how many batches stream past. The batch size for the spill is
`flush_rows` (default 500,000 transformed rows); the resident layer is the small
`y` a join or filter tests against.

To make the examples concrete, load the North Carolina counties shipped with sf
and project them to the state plane in metres. Projecting matters here: in a
planar CRS distances, areas, and buffers are Euclidean, and the recognised
operations run natively on the GEOS C API straight off the WKB column instead of
decoding each batch back to `sf`.

```{r load-nc}
nc <- st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)
nc <- st_transform(nc, 32119)        # NAD83 / North Carolina, metres
crs_nc <- st_crs(nc)
nrow(nc)
```

Write the polygons to a `.vtr` file with their geometry as a hex-WKB column.
This is the shape every streamed layer takes: attribute columns plus one string
column of WKB. In a real out-of-core workflow this file would be written once
from a large source and reused; here it is 100 counties standing in for the
billion-row case.

```{r write-poly}
f_poly <- tempfile(fileext = ".vtr")
write_vtr(data.frame(
  NAME  = nc$NAME,
  BIR74 = nc$BIR74,
  SID74 = nc$SID74,
  geometry = st_as_binary(st_geometry(nc), hex = TRUE)
), f_poly)
```

`tbl(f_poly)` now opens a lazy node over that file. Nothing is read until a
verb pulls it. Two functions bring a streamed result back to R: `collect()`
returns the underlying data.frame, geometry still a WKB string; `collect()`
decodes that string and reattaches the node's CRS, giving an ordinary `sf`
object ready to plot or hand to any sf function.

## Per-feature transforms with spatial_map

The widest door is `spatial_map()`: apply any per-feature sf operation to a
streamed layer one batch at a time. Buffering, simplifying, computing centroids,
reprojecting, extracting coordinates, measuring area. The function you pass takes
one `sf` batch and returns an `sf` object, an `sfc`, or a plain data.frame. The
active geometry of the return becomes the streamed geometry.

Buffer every county centroid by 15 km. The first step builds a centroid layer
on disk, the second streams a buffer over it.

```{r map-buffer}
cent <- tempfile(fileext = ".vtr")
write_vtr(data.frame(
  NAME = nc$NAME,
  geometry = st_as_binary(st_centroid(st_geometry(nc)), hex = TRUE)
), cent)

buffered <- tbl(cent) |>
  spatial_map(~ st_buffer(.x, 15000), crs = crs_nc)
buffered
```

`buffered` is a lazy node, not a materialized layer: the buffer has not run yet.
`collect()` pulls it through and decodes the geometry.

```{r map-buffer-plot}
b_sf <- collect(buffered)
plot(st_geometry(nc), border = "grey80", col = NA,
     main = "County centroids buffered by 15 km")
plot(st_geometry(b_sf), border = "#3366cc", col = "#3366cc22", add = TRUE)
```

The formula syntax (`~ st_buffer(.x, 15000)`) is rlang's: `.x` is the batch. A
named function works identically and is clearer for anything longer than one
call. The return can also be a plain data.frame rather than a geometry, which
turns the spatial step into a streamed feature-summary with the geometry
dropped. That is the way to compute per-feature area over a layer too large to
hold whole.

```{r map-area}
areas <- tbl(f_poly) |>
  spatial_map(~ data.frame(NAME = .x$NAME,
                           area_km2 = as.numeric(st_area(.x)) / 1e6),
              crs = crs_nc)
head(collect(areas))
```

Because the result is a node, it chains. A buffer feeding a filter feeding a
join is three streamed steps, each spilling to its own run-files, with the CRS
carried forward so you state it once. `flush_rows` controls the spill batch:
raise it for fewer, larger temporary files when rows are small, lower it when a
single feature's geometry is heavy. The default of 500,000 suits point and
small-polygon work.

## Smooth geometry with spatial_smooth

`spatial_smooth()` rounds the corners of every line and polygon by Chaikin
corner-cutting, one batch at a time. Each pass replaces a vertex with two points
along its adjacent edges, so sharp angles become a smooth curve; more
`iterations` give a smoother result with more vertices. It is computed directly
on the coordinates with no GEOS call, and open lines keep their endpoints.

```{r smooth}
zig <- st_linestring(rbind(c(0, 0), c(1, 1), c(2, 0), c(3, 1), c(4, 0)))
f_zig <- tempfile(fileext = ".vtr")
write_vtr(data.frame(
  id = 1L, geometry = st_as_binary(st_sfc(zig), hex = TRUE)), f_zig)

tbl(f_zig) |> spatial_smooth(iterations = 3) |> collect()
unlink(f_zig)
```

## Select by location with spatial_filter

The most-used vector tool keeps the features standing in a spatial relation to
another layer: points inside a study region, parcels touching a road, patches
within a buffer of the coast. `spatial_filter()` streams the large layer and
tests each batch against a small resident layer with an sf predicate, keeping or
dropping whole features. It is the spatial counterpart to `semi_join()`: rows
are filtered, never duplicated, and no columns are added, so the output carries
the input schema unchanged.

Sample 500 points across the state and store them as plain x/y columns. Point
input as coordinate columns is the headline case and the fully sf-free path:
each point is built in C and matched against the resident layer's spatial index
with no per-batch round-trip through sf.

```{r sample-points}
set.seed(1)
pts <- st_coordinates(st_sample(st_union(nc), 500))
fp <- tempfile(fileext = ".vtr")
write_vtr(data.frame(id = seq_len(nrow(pts)), x = pts[, 1], y = pts[, 2]), fp)
```

Keep the points falling inside a five-county region in the northwest.

```{r filter-region}
region <- nc[nc$NAME %in% c("Ashe", "Alleghany", "Surry", "Wilkes", "Watauga"),
             "NAME"]

keep_xy <- tbl(fp) |>
  spatial_filter(region, coords = c("x", "y"), crs = crs_nc) |>
  collect()
nrow(keep_xy)
```

Twenty-five of the 500 points land in the region. Plot the kept points against
the dropped ones.

```{r filter-plot}
plot(st_geometry(nc), border = "grey85", col = NA, main = "Select by location")
plot(st_geometry(region), border = "#cc3344", col = "#cc334411", add = TRUE)
points(pts, pch = 16, cex = 0.5, col = "grey70")
points(keep_xy$x, keep_xy$y, pch = 16, cex = 0.6, col = "#cc3344")
```

`negate = TRUE` inverts the test, keeping the 475 points outside the region. A
different `predicate` changes the relation: `st_within`, `st_covered_by`,
`st_touches`, `st_crosses`, and the rest of the sf binary predicates are all
accepted. `st_is_within_distance` takes its radius through `dist`, in CRS units,
and keeps every feature within that distance of the layer.

```{r filter-distance}
near <- tbl(fp) |>
  spatial_filter(region, predicate = st_is_within_distance,
                 coords = c("x", "y"), crs = crs_nc, dist = 30000)
nrow(collect(near))
```

Sixty points sit within 30 km of the region, more than the 25 strictly inside.
For the recognised predicates on planar data the test runs in C off the WKB
column or the raw coordinates; an unrecognised predicate, geographic data with
spherical geometry switched on, or an extra predicate argument falls back to the
per-batch sf loop, which preserves sf's exact semantics at the cost of decoding
each batch.

## Cut geometry with spatial_clip

`spatial_filter()` keeps or drops whole features; `spatial_clip()` cuts their
geometry along a mask. Clipping intersects each batch with the mask and keeps the
part inside it, the streamed equivalent of `st_intersection()` against a fixed
boundary. `erase = TRUE` keeps the part outside instead.

```{r clip}
mask_region <- st_union(st_geometry(region))

clipped <- tbl(f_poly) |> spatial_clip(mask_region, crs = crs_nc)
c_sf <- collect(clipped)
nrow(c_sf)
```

Twelve counties have area inside the region's bounding shape, and each comes back
trimmed to the part that overlaps it. Plot the clipped slivers over the full
state.

```{r clip-plot}
plot(st_geometry(nc), border = "grey85", col = NA,
     main = "Counties clipped to the region")
plot(st_geometry(c_sf), border = "#2a9d5c", col = "#2a9d5c33", add = TRUE)
```

The mask is unioned once into a single resident geometry, then each batch is cut
against it in C. As with the join and filter, the cut runs natively on the WKB
column for planar data and through `st_intersection()` / `st_difference()` for
the geographic-with-s2 and coordinate-assembled cases. Erase reverses the
keep: `tbl(f_poly) |> spatial_clip(mask_region, erase = TRUE)` returns the 95
counties with any area outside the region, each trimmed to that outside part.

## Split geometry with spatial_split

`spatial_split()` cuts each feature against a small resident `blade` layer, the
QGIS "split with lines". A polygon is divided into the faces the blade carves
out, a line into the arcs between crossings, and each piece is emitted as its own
row with the source attributes copied. A feature the blade does not cross passes
through whole. With `extract = "points"` the verb returns the crossing points
instead.

```{r split}
square <- st_polygon(list(rbind(c(0, 0), c(4, 0), c(4, 4), c(0, 4), c(0, 0))))
blade  <- st_sfc(st_linestring(rbind(c(2, -1), c(2, 5))))

f_sq <- tempfile(fileext = ".vtr")
write_vtr(data.frame(
  id = 1L, geometry = st_as_binary(st_sfc(square), hex = TRUE)), f_sq)

tbl(f_sq) |> spatial_split(blade) |> collect()
unlink(f_sq)
```

## Tag features with spatial_join

To attach a layer's attributes rather than filter or cut, `spatial_join()`
streams the large left side and joins each batch against a small resident right
side. The dominant workload is tagging a huge point set with the polygon it
falls in: which county each occurrence sits in, which census tract each address
belongs to. The billion-row left stream never materializes while the polygon
layer stays in RAM. On the native path the join is a lazy node like any other
verb: the joined rows stream straight into whatever comes next, a `count()` or a
`write_vtr()`, without being spilled first.

```{r join-tag}
tagged <- tbl(fp) |>
  spatial_join(nc["NAME"], coords = c("x", "y"), crs = crs_nc)
tdf <- collect(tagged)
head(tdf[, c("id", "NAME")])
```

Every point now carries the `NAME` of its county. The default predicate is
`st_intersects` and the default is a left join, so an unmatched left row is kept
once with NA in the right columns; `left = FALSE` drops the unmatched rows for an
inner join. Colour the points by the county they were tagged with.

```{r join-plot}
tagged_sf <- st_as_sf(tdf, coords = c("x", "y"), crs = crs_nc)
plot(st_geometry(nc), border = "grey85", col = NA,
     main = "Points tagged with their county")
plot(tagged_sf["NAME"], pch = 16, cex = 0.5, add = TRUE)
```

When only the attributes matter, `keep_geom = FALSE` drops the left geometry
column. With `coords =` the point geometries are then never built at all, which
is the cheap way to tag points and aggregate:

```{r join-count}
tbl(fp) |>
  spatial_join(nc["NAME"], coords = c("x", "y"), crs = crs_nc,
               keep_geom = FALSE) |>
  count(NAME) |>
  collect() |>
  head()
```

Columns present on both sides are disambiguated with `suffix` (default
`c(".x", ".y")`), exactly as `st_join()` does. A different `join` predicate
changes the relation; `st_nearest_feature` finds the closest right feature to
each left one, which is how you snap points to the nearest road or station.

```{r join-nearest}
ncent <- st_sf(NAME = nc$NAME, geometry = st_centroid(st_geometry(nc)))
nearest <- tbl(fp) |>
  spatial_join(ncent, join = st_nearest_feature,
               coords = c("x", "y"), crs = crs_nc)
nrow(collect(nearest))
```

### When both sides are larger than RAM

The resident-`y` path assumes the right side fits in memory. When it does not,
pass `partition = grid(cellsize)` and a streamed `vectra_node` as `y`. Both
inputs are binned to a uniform grid and joined one shard at a time, so neither
side is ever fully resident. A coordinate maps to the cell

```
cell  =  floor( (coord - origin) / cellsize )
```

per axis. Each left feature is assigned to the single cell of its reference
point; each right feature is replicated to every cell its bounding box overlaps.
A left row is therefore emitted exactly once, and the result equals the resident
join for point left geometry.

```{r join-partition}
g_poly <- tempfile(fileext = ".vtr")
write_vtr(data.frame(
  NAME = nc$NAME,
  geometry = st_as_binary(st_geometry(nc), hex = TRUE)
), g_poly)

tagged2 <- tbl(fp) |>
  spatial_join(tbl(g_poly), coords = c("x", "y"), crs = crs_nc,
               partition = grid(80000))
t2 <- collect(tagged2)
sum(!is.na(t2$NAME))
```

All 500 points are tagged, matching the resident join. The `cellsize` is the one
tuning knob: large enough that one cell's features fit in memory, small enough
that the shards stay balanced. For point-in-polygon tagging any cell larger than
the polygons works; for an extended-on-extended join choose it larger than the
left features. The nearest-feature predicate is not local to one cell, so under
`partition` a left point searches its own cell and the eight around it; set
`cellsize` at least as large as the largest expected nearest distance.

### Same answer as resident sf

Streaming changes the memory profile, not the result. The streamed join returns
exactly what `st_join()` returns on the whole layer held in memory, feature for
feature. Run both on the 500-point set and compare.

```{r join-equality}
streamed <- collect(
  tbl(fp) |> spatial_join(nc["NAME"], coords = c("x", "y"), crs = crs_nc))

resident <- st_join(
  st_as_sf(collect(tbl(fp)), coords = c("x", "y"), crs = crs_nc, remove = FALSE),
  nc["NAME"], join = st_intersects)

all.equal(streamed$NAME[order(streamed$id)],
          resident$NAME[order(resident$id)])
```

The county tags match for every point. The equality is the contract: a streamed
verb is the resident sf call run batch by batch, so the choice between them is a
memory decision, not an accuracy trade-off.

### At scale

The point of streaming shows when the layer stops fitting in memory. Scatter
200,000 points across the state's bounding box and filter them to the five-county
region. Only one batch is resident at a time, so the same code that ran on 500
points runs on 200,000 with a flat memory profile.

```{r filter-scale}
set.seed(42)
bb <- st_bbox(nc)
n_big <- 2e5
big <- data.frame(id = seq_len(n_big),
                  x = runif(n_big, bb["xmin"], bb["xmax"]),
                  y = runif(n_big, bb["ymin"], bb["ymax"]))
fbig <- tempfile(fileext = ".vtr")
write_vtr(big, fbig)

kept <- tbl(fbig) |>
  spatial_filter(region, coords = c("x", "y"), crs = crs_nc) |>
  collect()
nrow(kept)
```

A few thousand of the 200,000 land in the five counties. At a billion points the
file would be larger but the resident set would not grow: the engine still pulls
one batch, tests it in C against the region's spatial index, and moves on. That
is the whole proposition. The operation that a desktop GIS can only run on a
layer it can open is the same operation, run past a fixed memory budget.

```{r cleanup-scale, include = FALSE}
unlink(fbig)
```

## Nearest features with spatial_knn

Where `spatial_join()` attaches the single nearest feature, `spatial_knn()`
returns the `k` nearest features per streamed point, one row per pair, each with
the rank (1 is nearest) and the distance. The candidate layer `y` is held
resident while the points stream.

```{r knn}
towns <- suppressWarnings(st_centroid(st_geometry(nc)))[1:5]
towns <- st_sf(town = nc$NAME[1:5], geometry = towns)

set.seed(1)
pts <- suppressWarnings(st_coordinates(st_sample(nc, 100)))
f_pts <- tempfile(fileext = ".vtr")
write_vtr(data.frame(id = seq_len(nrow(pts)), x = pts[, 1], y = pts[, 2]), f_pts)

tbl(f_pts) |>
  spatial_knn(towns, k = 2, coords = c("x", "y"), crs = crs_nc, y_id = "town") |>
  collect() |> head()
unlink(f_pts)
```

## Merge geometries with spatial_dissolve

Dissolve unions the geometries within each group into one feature, the GIS
"Dissolve" tool: counties into states, parcels into ownership blocks. Unlike the
per-batch verbs, it needs every geometry of a group together to union them, so it
rides a different tier. The input is spilled once and routed into one disjoint
shard per group in a single bounded pass; each shard is then read back and
unioned. Peak memory is the routing budget during the pass, then one group's
geometries while that group unions. Partition on a key whose groups fit in
memory.

Split the counties into two bands by their 1974 SIDS count and merge each band
into a single feature, summing births along the way.

```{r dissolve}
nc$band <- ifelse(nc$SID74 > 5, "high", "low")
fb <- tempfile(fileext = ".vtr")
write_vtr(data.frame(
  band = nc$band, BIR74 = nc$BIR74,
  geometry = st_as_binary(st_geometry(nc), hex = TRUE)
), fb)

merged <- tbl(fb) |>
  spatial_dissolve(by = "band", crs = crs_nc,
                   .fun = list(births = function(d) sum(d$BIR74)))
m_sf <- collect(merged)
m_sf
```

Two features come back, one per band, each carrying the summed `births` of its
counties. The `.fun` argument is a named list of summaries; each function takes
the group's data.frame and returns one value, and the list name becomes the
output column. Plot the two bands.

```{r dissolve-plot}
plot(m_sf["band"], main = "Counties dissolved into two SIDS bands")
```

With no `by`, the whole layer dissolves into one feature, which is the fast way
to compute a layer's outline out of core. On planar data each group is unioned
natively off the WKB column; geographic data with s2 on, or an extra `st_union()`
argument such as `is_coverage = TRUE`, unions through sf instead.

## Split overlaps apart with spatial_overlay

Dissolve merges overlapping geometry; overlay splits it. `spatial_overlay()`
takes a polygon layer that overlaps itself and cuts it along every overlap into
disjoint pieces, returning one row per piece per covering polygon. Where `k`
polygons overlap, that piece appears `k` times, each row carrying one source
polygon's attributes. This is the operation a GIS exposes as "Union (single
layer)": the overlap is retained once per contributing feature rather than
dissolved away. It answers questions like "which protected areas, designated in
which years, cover this exact patch of sea".

Three overlapping squares designated in different years make the smallest
meaningful example.

```{r overlay}
sq <- function(a, b) st_polygon(list(rbind(
  c(a, 0), c(b, 0), c(b, 1), c(a, 1), c(a, 0))))
polys <- st_sf(year = c(1990L, 2010L, 2000L),
               geometry = st_sfc(sq(0, 2), sq(1, 3), sq(1.5, 3.5)))

pieces <- collect(spatial_overlay(polys))
nrow(pieces)
length(unique(pieces$piece_id))
```

The three squares decompose into five disjoint pieces, returned as nine rows
because the overlapping pieces are repeated once per covering square. The
`piece_id` column keys the pieces: rows sharing an id are the same patch of
ground seen from different source polygons. Resolve the duplication with a
grouped slice. Earliest designation year wins:

```{r overlay-resolve}
first <- spatial_overlay(polys) |>
  group_by(piece_id) |>
  slice_min(year, n = 1, with_ties = FALSE) |>
  collect()
nrow(first)

plot(first["year"], main = "Overlay pieces, earliest year wins")
```

Five pieces remain, one per disjoint patch, each labelled with the earliest year
that covers it. Swapping `slice_min` for `slice_max` gives the latest; any
grouped reduction works, because the pieces are an ordinary streamed node by this
point.

Three properties make the overlay trustworthy on real data. First, correctness
is checked, not assumed: for a valid decomposition the piece areas covering an
input must sum to that input's area, and the engine verifies this invariant per
batch, warning if coverage drifts past a relative `1e-4` rather than returning
silently wrong. Second, memory is bounded by tiling: a connected cluster of many
large overlapping polygons can node into far more linework than the input, so
clusters too large for the budget are tiled over their own extent and clipped,
with each feature cleaned exactly once. Third, the snapping grid is explicit:
coordinates are snapped to a grid derived from their magnitude so the
near-duplicate boundaries that overlapping polygons share coincide, which is what
keeps the pieces disjoint instead of leaving sliver faces.

Unlike the other verbs, `spatial_overlay()` takes a resident `sf` object, not a
lazy node: building the overlap graph needs the geometries in memory. The
exploded result, typically several times larger than the input, is what streams
to disk. `mem_limit` caps the peak working set and `threads` sets the parallel
overlay width; raise both for speed on a big machine, lower `mem_limit` for
tighter memory.

### A second layer

Passing a second layer `y` nodes two layers into one planar partition and
carries the attributes of the covering `x`-record and `y`-record onto each
piece. A `how` argument selects which pieces to keep: `"intersection"` (covered
by both, the default), `"union"` (every piece of either, the absent side filled
with `NA`), `"identity"` (all of `x` split by `y`), or `"symdiff"` (pieces in
exactly one layer).

```{r overlay-two}
zones <- st_sf(zone = c("A", "B"),
               geometry = st_sfc(sq(0, 1.5), sq(1.5, 3)))

inter <- spatial_overlay(polys, zones, how = "intersection") |> collect()
inter
```

Each piece now carries both its `year` (from the squares) and its `zone` (from
the zone layer). `vars_y` narrows the carried `y` columns, and a name shared
with an `x` column is disambiguated with a `.x` / `.y` suffix. A file-path `y`
is read in batches with `layer_y` / `query_y`, the same way `x` is.

## Materializing and round-tripping

Three exits bring a streamed result out. `collect()` returns the data.frame with
geometry still a WKB string, useful when the next step is a
non-spatial verb or another `.vtr` write. `collect()` decodes the WKB and
reattaches the node's CRS, giving an `sf` object. `write_vtr()` on a node streams
the result straight to a new `.vtr` file without ever holding it whole, so a
multi-step spatial pipeline can land its output on disk under the same fixed
memory budget it ran in.

```{r roundtrip}
out <- tempfile(fileext = ".vtr")
tbl(fp) |>
  spatial_filter(region, coords = c("x", "y"), crs = crs_nc) |>
  write_vtr(out)
nrow(collect(tbl(out)))
```

The CRS lives on the node, not in the file. A pipeline that opens projected data,
maps, filters, and joins keeps the projection because each verb carries it
forward; you state `crs =` once at the first step that needs it, or let it
inherit from the upstream node. The `.vtr` file itself stores only the WKB bytes,
so reopening a written file and calling `collect()` needs the CRS supplied
again if you want it labelled.

## Geometry functions in expressions

vectra carries geometry through the engine as hex-encoded WKB in an ordinary
string column. The verbs above wrap whole sf operations around that column.
This section covers the other half of the spatial surface: a family of `st_*` functions that work inside the expression verbs themselves,
so a measure, a predicate, or a geometry transform is just another term in
`mutate()`, `filter()`, or `summarise()`.

These functions run on the GEOS C library straight off the WKB column, one row
at a time, with no per-batch round-trip through sf. `filter(st_area(geometry) >
1e6)` prunes the stream in C, and `mutate(here = st_centroid(geometry))` adds a
new WKB geometry column that any later verb can read. The per-row decode is
parallelised with OpenMP, so a measure over a large layer uses every core.


### A layer to work on

The examples use the North Carolina counties shipped with sf. Writing the layer
to a `.vtr` is the usual first step: the geometry becomes a hex-WKB string
column (named `geometry` by convention), and the attributes ride alongside it.

```{r gx-data}
nc <- st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)

f <- tempfile(fileext = ".vtr")
write_vtr(data.frame(
  NAME     = nc$NAME,
  BIR74    = nc$BIR74,
  geometry = st_as_binary(st_geometry(nc), hex = TRUE)
), f)

tbl(f)
```

The counties are stored in longitude and latitude, so every measure below is
planar in those units, the same convention the streaming verbs follow. Project
the layer first if you need metric areas or distances.

### Measures

A measure reads a geometry and returns a number, so it drops into `mutate()` as
an ordinary column.

```{r gx-measures}
tbl(f) |>
  mutate(area   = st_area(geometry),
         perim  = st_perimeter(geometry),
         nverts = st_npoints(geometry)) |>
  select(NAME, area, perim, nverts) |>
  collect() |>
  head()
```

`st_length()` returns the boundary length of a polygon (an alias of
`st_perimeter()`) and the line length of a linestring. `st_ngeometries()`
counts the parts of a multi-geometry. `st_x()` and `st_y()` read the coordinate
of a point and return `NA` for anything that is not a point, which makes them
most useful on a centroid:

```{r gx-coords}
tbl(f) |>
  mutate(centroid = st_centroid(geometry),
         cx = st_x(centroid),
         cy = st_y(centroid)) |>
  select(NAME, cx, cy) |>
  collect() |>
  head()
```

A geometry-valued function such as `st_centroid()` produces a new WKB column
(`centroid` above), and the next term reads it like any other column. That is
the whole mechanism: geometry in, geometry or a scalar out, all as columns.

### Predicates

A unary predicate tests one geometry: `st_is_valid()`, `st_is_empty()`,
`st_is_simple()`. A binary predicate tests a topological relation against a
second geometry: `st_intersects()`, `st_within()`, `st_contains()`,
`st_overlaps()`, `st_touches()`, `st_crosses()`, `st_equals()`,
`st_disjoint()`, `st_covers()`, `st_covered_by()`.

In `filter()` a predicate keeps the rows where the relation holds, the
geometry-expression form of select-by-location:

```{r gx-predicate-filter}
aoi <- st_as_sfc(st_bbox(c(xmin = -80, ymin = 35, xmax = -78, ymax = 36.5)),
                 crs = st_crs(nc))

tbl(f) |>
  filter(st_intersects(geometry, aoi)) |>
  collect() |>
  nrow()
```

The second geometry here is a constant `sf` object. It is parsed once and shared
read-only across every row, so testing a whole stream against one area of
interest stays cheap. A multi-feature object is unioned to a single geometry
first.

In `mutate()` the same call returns a logical column, ready for a later verb:

```{r gx-predicate-mutate}
tbl(f) |>
  mutate(near_raleigh = st_intersects(geometry, aoi)) |>
  filter(near_raleigh) |>
  select(NAME) |>
  collect() |>
  head()
```

### Distance

`st_distance()` returns the shortest planar distance to a second geometry,
again a constant or another column:

```{r gx-distance}
raleigh <- st_sfc(st_point(c(-78.64, 35.78)), crs = st_crs(nc))

tbl(f) |>
  mutate(centroid   = st_centroid(geometry),
         d_raleigh  = st_distance(centroid, raleigh)) |>
  select(NAME, d_raleigh) |>
  arrange(d_raleigh) |>
  collect() |>
  head()
```

When the second argument is a geometry column instead of a constant, the
distance is computed row by row between the two columns.

### Aggregating a measure

Because a measure is an ordinary numeric column, it aggregates like one. A
grouped `summarise()` over a measure is a zonal total computed entirely in the
stream:

```{r gx-summarise}
tbl(f) |>
  mutate(area = st_area(geometry)) |>
  summarise(total_area = sum(area), counties = n()) |>
  collect()
```

### Transforms

A transform returns a geometry, so it builds a new WKB column. Materialise it
with `collect()`, which reads the WKB column back into an `sf` object (point
it at the column with `geom =`, and pass the `crs` the layer was stored in).

```{r gx-transforms}
hulls <- tbl(f) |>
  mutate(geometry = st_convex_hull(geometry)) |>
  select(NAME, geometry) |>
  collect(geom = "geometry", crs = st_crs(nc))

hulls
```

The transform set is `st_centroid()`, `st_point_on_surface()` (a point
guaranteed to lie on the geometry), `st_boundary()`, `st_envelope()` (the
bounding rectangle), `st_convex_hull()`, `st_make_valid()` (repair an invalid
geometry), and two parameterised forms: `st_buffer(g, dist)` and
`st_simplify(g, tol)`. Buffering each county and reading the areas back:

```{r gx-buffer}
tbl(f) |>
  mutate(geometry = st_buffer(geometry, 0.1)) |>
  select(NAME, geometry) |>
  collect(geom = "geometry", crs = st_crs(nc)) |>
  st_area() |>
  head()
```

`st_geometry_type()` returns the GEOS type name (`"Point"`, `"Polygon"`,
`"MultiPolygon"`, and so on) as a string column.

### The second geometry of a binary op

For `st_distance()` and the binary predicates, the second argument can be:

- another geometry column, compared row by row;

- a constant `sf` or `sfc` object, parsed once and reused across the stream (a
  multi-feature object is unioned to one geometry);
- a hex-WKB string.

### Missing geometry

A missing (`NA`) or unparseable geometry, or an operation with no answer (the
coordinate of a non-point, the distance to a missing geometry), yields `NA` for
that row rather than stopping the query:

```{r gx-na}
g <- tempfile(fileext = ".vtr")
write_vtr(data.frame(
  id = 1:4,
  geometry = c(st_as_binary(st_geometry(nc)[1:3], hex = TRUE), NA)
), g)

tbl(g) |>
  mutate(area = st_area(geometry)) |>
  collect()
```

### Where this fits

The `st_*` expressions are the scalar, per-row layer of vectra's spatial
surface. They cover measures, predicates, and the common single-geometry
transforms at the price of a column term, with no sf object built per batch. For
an arbitrary per-feature transform that has no `st_*` form, reach for
`spatial_map()`, which runs any sf-in, sf-out function over each feature. For
constructions that read a whole feature set at once (dissolves, overlays, hulls
of a group, planar topology), the set-wise `spatial_*` verbs above and
[Coverage and topology](#coverage-and-topology) are the tools.

```{r gx-cleanup, include = FALSE}
unlink(c(f, g))
```

## Coverage and topology

The per-feature verbs above treat each geometry on its own. A second family of verbs reads a whole feature set at once: the
constructions and topology operations that need every feature in a group
together, not one row at a time. These ride a shared partition tier. A `by`
argument routes the layer into one shard per group, each shard is processed as
an independent coverage, and the shards stream so the resident budget is one
group rather than the whole layer.

This section covers building polygons from lines, cleaning a polygon coverage,
decomposing geometry into its parts, set-wise constructions, and linear
referencing. Every block runs on the optional
[sf](https://r-spatial.github.io/sf/) package.


### From lines to polygons

`spatial_polygonize()` builds the polygonal faces enclosed by a line network,
the inverse of taking polygon boundaries. A group's lines are unioned and noded,
then the faces of that arrangement are returned, one per row.

```{r ct-polygonize}
grid <- st_sfc(
  st_linestring(rbind(c(0, 0), c(2, 0))),
  st_linestring(rbind(c(0, 1), c(2, 1))),
  st_linestring(rbind(c(0, 2), c(2, 2))),
  st_linestring(rbind(c(0, 0), c(0, 2))),
  st_linestring(rbind(c(1, 0), c(1, 2))),
  st_linestring(rbind(c(2, 0), c(2, 2))))

f <- tempfile(fileext = ".vtr")
write_vtr(data.frame(geometry = st_as_binary(grid, hex = TRUE)), f)

tbl(f) |> spatial_polygonize() |> collect()
```

`spatial_line_merge()` is the line counterpart of a dissolve: it sews segments
that meet end to end into maximal linestrings, one row per chain. Segments
meeting at a junction of degree greater than two stay separate.

```{r ct-line-merge}
seg <- st_sfc(
  st_linestring(rbind(c(0, 0), c(1, 0))),
  st_linestring(rbind(c(1, 0), c(2, 0))),
  st_linestring(rbind(c(2, 0), c(3, 0))))

g <- tempfile(fileext = ".vtr")
write_vtr(data.frame(geometry = st_as_binary(seg, hex = TRUE)), g)

tbl(g) |> spatial_line_merge() |> collect()
```

### Cleaning a coverage

Real coverages arrive with jittered shared borders, sub-pixel slivers, and more
detail than a map needs. Three verbs clean them while keeping adjacent polygons
edge-matched.

`spatial_snap_grid()` rounds coordinates to a regular grid and repairs the
result, one batch at a time. It is the fixed-precision snap-rounding that
`spatial_overlay()` applies internally, exposed as a standalone verb, so a layer
can be pre-noded to a common precision without running a full overlay.

```{r ct-snap-grid}
p <- st_polygon(list(rbind(c(0.04, 0.03), c(1.02, 0.01),
                           c(0.98, 1.03), c(0.01, 0.97), c(0.04, 0.03))))
h <- tempfile(fileext = ".vtr")
write_vtr(data.frame(id = 1L, geometry = st_as_binary(st_sfc(p), hex = TRUE)), h)

tbl(h) |> spatial_snap_grid(0.1) |> collect()
```

`spatial_snap()` instead pulls a streamed layer onto a resident reference layer:
vertices and edges within `tolerance` of `y` are moved onto it, so two layers
that should share a boundary line up exactly.

```{r ct-snap}
ref  <- st_sfc(st_linestring(rbind(c(0, 0), c(10, 0))))
line <- st_linestring(rbind(c(0, 0.2), c(5, 0.1), c(10, 0.2)))

sn <- tempfile(fileext = ".vtr")
write_vtr(data.frame(id = 1L, geometry = st_as_binary(st_sfc(line), hex = TRUE)), sn)

tbl(sn) |> spatial_snap(ref, tolerance = 0.5) |> collect()
```

`spatial_eliminate()` cleans a polygon coverage by absorbing every feature
smaller than `max_area` into a neighbour, the QGIS "Eliminate". Each sliver joins
the neighbour it shares the longest border with, or the largest-area neighbour
with `into = "largest_area"`. An area-rooted union-find collapses chains of
slivers so a connected run flows to its single largest member, whose attributes
survive.

```{r ct-eliminate}
big    <- st_polygon(list(rbind(
  c(0, 0), c(10, 0), c(10, 10), c(0, 10), c(0, 0))))
sliver <- st_polygon(list(rbind(
  c(10, 0), c(10.3, 0), c(10.3, 10), c(10, 10), c(10, 0))))

e <- tempfile(fileext = ".vtr")
write_vtr(data.frame(
  id = c("keep", "sliver"),
  geometry = st_as_binary(st_sfc(big, sliver), hex = TRUE)), e)

tbl(e) |> spatial_eliminate(max_area = 5) |> collect()
```

The thin sliver is gone, absorbed into the square it bordered, and the survivor
keeps the `keep` attributes.

`spatial_simplify()` simplifies a polygon coverage without tearing shared edges.
Boundaries are unioned so a shared border is one line, noded into arcs, each arc
simplified once with its junction endpoints pinned, then re-polygonized. Adjacent
polygons stay edge-matched with no slivers. This is the topology-preserving
simplification a per-feature `spatial_map(~ sf::st_simplify())` cannot give,
because that simplifies each polygon's copy of a shared border independently.

```{r ct-simplify}
p1 <- st_polygon(list(rbind(
  c(0, 0), c(1, 0), c(1, 0.5), c(1, 1), c(0, 1), c(0, 0))))
p2 <- st_polygon(list(rbind(
  c(1, 0), c(2, 0), c(2, 1), c(1, 1), c(1, 0.5), c(1, 0))))

s <- tempfile(fileext = ".vtr")
write_vtr(data.frame(
  id = c("a", "b"),
  geometry = st_as_binary(st_sfc(p1, p2), hex = TRUE)), s)

tbl(s) |> spatial_simplify(tolerance = 0.6) |> collect()
```

### Decomposing geometry

`spatial_explode()` splits every multipart geometry into its single-part
components, a `MULTIPOLYGON` into one row per polygon and so on, copying the
source attributes onto each part. An optional `part` argument names a 1-based
part-index column. It is the streaming "multipart to singleparts" tool and the
inverse of `spatial_dissolve()`.

```{r ct-explode}
mp <- st_multipolygon(list(
  list(rbind(c(0, 0), c(1, 0), c(1, 1), c(0, 1), c(0, 0))),
  list(rbind(c(2, 2), c(3, 2), c(3, 3), c(2, 3), c(2, 2)))))

m <- tempfile(fileext = ".vtr")
write_vtr(data.frame(
  id = 1L, geometry = st_as_binary(st_sfc(mp), hex = TRUE)), m)

tbl(m) |> spatial_explode(part = "part_id") |> collect()
```

`spatial_topology()` decomposes a polygon coverage into the arcs of its planar
topology. The unioned boundaries are noded so a shared border becomes one arc,
tagged with the identifiers of the (up to two) polygons on either side: two for
an internal shared edge, one for an outer edge. It is the inverse of
`spatial_polygonize()`.

```{r ct-topology}
q1 <- st_polygon(list(rbind(c(0, 0), c(1, 0), c(1, 1), c(0, 1), c(0, 0))))
q2 <- st_polygon(list(rbind(c(1, 0), c(2, 0), c(2, 1), c(1, 1), c(1, 0))))

tp <- tempfile(fileext = ".vtr")
write_vtr(data.frame(
  id = c("a", "b"),
  geometry = st_as_binary(st_sfc(q1, q2), hex = TRUE)), tp)

tbl(tp) |> spatial_topology(id = "id") |> collect()
```

The shared edge between `a` and `b` appears once, with both face columns filled;
each outer edge appears once with the second face `NA`.

`spatial_centerline()` traces the centerline (medial axis) of each polygon from
the Voronoi diagram of its densified boundary: the Voronoi edges that fall inside
the polygon are its skeleton, merged into maximal lines. `density` sets the
boundary sampling and `prune` drops the short spurs the skeleton grows toward
convex corners. It is the usual approximation for a river or road centerline from
a filled shape; non-polygon geometry passes through unchanged.

```{r ct-centerline}
strip <- st_polygon(list(rbind(
  c(0, 0), c(10, 0), c(10, 2), c(0, 2), c(0, 0))))

cl <- tempfile(fileext = ".vtr")
write_vtr(data.frame(geometry = st_as_binary(st_sfc(strip), hex = TRUE)), cl)

tbl(cl) |> spatial_centerline(density = 0.25, prune = 0.5) |> collect()
```

### Set-wise constructions

`spatial_construct()` builds a geometry construction from a whole feature set,
the constructions a per-feature `spatial_map()` cannot express. A `kind`
argument selects it: `"convex_hull"`, `"concave_hull"`, `"envelope"`,
`"oriented_box"`, `"enclosing_circle"`, `"inscribed_circle"`, `"pole"` (the pole
of inaccessibility, the centre of the maximum inscribed circle), `"voronoi"`,
and `"delaunay"`. A `by` argument builds one construction per group; the
enclosing kinds emit one feature per group, the tessellations one polygon per
cell.

```{r ct-construct-hull}
nc <- st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)
nc$band <- nc$SID74 > 5

n <- tempfile(fileext = ".vtr")
write_vtr(data.frame(
  band = nc$band,
  geometry = st_as_binary(st_geometry(nc), hex = TRUE)), n)

tbl(n) |>
  spatial_construct("convex_hull", by = "band", crs = st_crs(nc)) |>
  collect()
```

A tessellation kind returns one polygon per cell, each carrying the group's `by`
values:

```{r ct-construct-voronoi}
xy <- rbind(c(0, 0), c(1, 0), c(1, 1), c(0, 1), c(0.4, 0.6))
cloud <- st_sfc(lapply(seq_len(nrow(xy)), function(i) st_point(xy[i, ])))

v <- tempfile(fileext = ".vtr")
write_vtr(data.frame(geometry = st_as_binary(cloud, hex = TRUE)), v)

tbl(v) |> spatial_construct("voronoi") |> collect()
```

### Linear referencing

`spatial_locate()` places streamed points along a resident line layer: each
point gets its nearest line's identifier, the measure (distance along that line
to the point's projection), and the perpendicular offset, with an optional `snap`
onto the line. The inverse, a measure back to a point, is
`sf::st_line_interpolate()` through `spatial_map()`.

```{r ct-locate}
roads <- st_sf(road = c("main", "side"), geometry = st_sfc(
  st_linestring(rbind(c(0, 0), c(10, 0))),
  st_linestring(rbind(c(0, 5), c(0, 15)))))

pts <- data.frame(id = 1:2, x = c(3, 1), y = c(1, 9))
l <- tempfile(fileext = ".vtr")
write_vtr(pts, l)

tbl(l) |>
  spatial_locate(roads, coords = c("x", "y"), y_id = "road") |>
  collect()
```

Point 1 is one unit off `main` at measure 3; point 2 is one unit off `side` at
measure 4. With `snap = TRUE` each point's geometry is replaced by its
projection onto the matched line.

```{r ct-cleanup, include = FALSE}
unlink(c(f, g, h, sn, e, s, m, tp, cl, n, v, l))
```

## Network analysis

A line layer becomes a routable graph when its segments are joined at their
endpoints and weighted by a traversal cost. vectra builds that graph once and
holds it resident in a `vectra_network` object, the network counterpart of the
small resident `sf` layer the streaming verbs compare against. Origins then
stream past the resident graph: each batch of start points is routed in native
C, so a national set of origins flows past a fixed memory budget while the graph
stays in memory.

The graph and the shortest-path solver are pure C, a binary-heap Dijkstra over a
compressed adjacency, with no `igraph` dependency. The solver parallelises over a
batch of origins with OpenMP.


### Building a network

`spatial_network()` takes a line layer and returns a `vectra_network`. Endpoints
are snapped to shared nodes (exactly coincident by default, or within
`tolerance` CRS units), and each edge is weighted by its geometry length unless a
`weight` column names another cost.

The examples use a small grid of unit streets.

```{r nw-build}
mk <- function(x1, y1, x2, y2)
  st_linestring(rbind(c(x1, y1), c(x2, y2)))

streets <- st_sfc(
  mk(0, 0, 1, 0), mk(1, 0, 2, 0), mk(2, 0, 3, 0),   # bottom row
  mk(0, 1, 1, 1), mk(1, 1, 2, 1), mk(2, 1, 3, 1),   # middle row
  mk(0, 2, 1, 2), mk(1, 2, 2, 2), mk(2, 2, 3, 2),   # top row
  mk(0, 0, 0, 1), mk(0, 1, 0, 2),                   # left verticals
  mk(1, 0, 1, 1), mk(1, 1, 1, 2),
  mk(2, 0, 2, 1), mk(2, 1, 2, 2),
  mk(3, 0, 3, 1), mk(3, 1, 3, 2))                   # right verticals

net <- spatial_network(streets)
net
```

### Shortest routes

`spatial_route()` streams a layer of origin points past the network and returns
the shortest path from each origin to one or more destinations. With
`geometry = FALSE` it returns only the cost, so a set of destinations per origin
is an origin-destination cost matrix in long form.

```{r nw-route-cost}
f <- tempfile(fileext = ".vtr")
write_vtr(data.frame(id = 1L, x = 0, y = 0), f)

dest <- st_sfc(st_point(c(3, 2)))

tbl(f) |>
  spatial_route(net, to = dest, coords = c("x", "y"), geometry = FALSE) |>
  collect()
```

The cost is the five unit steps from the bottom-left corner to the top-right.
With `geometry = TRUE` (the default) each row also carries the route line, ready
to materialise with `collect()`:

```{r nw-route-geom}
tbl(f) |>
  spatial_route(net, to = dest, coords = c("x", "y")) |>
  collect()
```

Several destinations, named through `to_id`, give one row per origin and
destination, the cost matrix in long form:

```{r nw-od-matrix}
g <- tempfile(fileext = ".vtr")
write_vtr(data.frame(id = 1:2, x = c(0, 3), y = c(0, 0)), g)

dests <- st_sf(
  name = c("top_left", "top_right"),
  geometry = st_sfc(st_point(c(0, 2)), st_point(c(3, 2))))

tbl(g) |>
  spatial_route(net, to = dests, to_id = "name", geometry = FALSE,
                coords = c("x", "y")) |>
  collect()
```

Each origin carries its own attributes (`id` here) onto every route row.
Unreachable pairs return an infinite cost rather than dropping the row.

### Service areas and isochrones

`spatial_service_area()` returns, per origin, the part of the network reachable
within a cost budget. A vector of budgets returns nested travel-cost bands, one
row per origin and band. The `output` argument selects the shape: the reachable
nodes as a multipoint, the reachable edges as lines, or their convex hull as a
service-area polygon.

```{r nw-service-nodes}
tbl(f) |>
  spatial_service_area(net, cost = c(1, 2), output = "nodes",
                       coords = c("x", "y")) |>
  collect()
```

The first band reaches the start node and its two immediate neighbours; the
second band reaches everything within two steps. Asking for `output = "polygon"`
wraps each band in its hull, the usual drive-time catchment:

```{r nw-service-polygon}
tbl(f) |>
  spatial_service_area(net, cost = 2, output = "polygon",
                       coords = c("x", "y")) |>
  collect() |>
  st_area()
```

### Weighted and directed networks

By default every edge is two-way and weighted by length. A `weight` column sets
a different traversal cost (a travel time, a toll), and `directed = TRUE` makes
edges one-way. The `direction` column then carries one-way codes (`"B"`
two-way, `"FT"` along the digitised direction, `"TF"` against it, `"N"` closed),
and `weight_to` gives the reverse cost on a two-way edge.

```{r nw-weighted}
streets_df <- st_sf(
  cost = runif(length(streets), 1, 3),
  geometry = streets)

wnet <- spatial_network(streets_df, weight = "cost")

tbl(f) |>
  spatial_route(wnet, to = dest, coords = c("x", "y"), geometry = FALSE) |>
  collect()
```

The route now minimises the summed `cost` weight instead of the step count.

### The streaming model

The graph is the resident budget: it is built once and held in the
`vectra_network` object for the life of the queries against it. The origin layer
is what streams. Each batch of origins is solved against the resident graph with
one Dijkstra per origin, run in parallel across the batch, so the query side
scales by streaming while memory tracks the graph plus one batch. This is the
network analogue of the resident-`y` pattern the
streaming spatial verbs use for a small locator layer.

```{r nw-cleanup, include = FALSE}
unlink(c(f, g))
```

## The cost model

Each verb states what it holds resident, so the toolbox reads as a cost model.
Three tiers cover everything here.

| Tier | Verbs | Resident set |
|---|---|---|
| Monoid fold | `spatial_map`, `spatial_filter`, `spatial_clip`, `rasterize` | one batch + small `y` |
| Partition | `spatial_dissolve`, partitioned `spatial_join` | routing budget, then one shard |
| All-to-all | `spatial_overlay` | one overlap cluster (tiled) |

The fold tier is the cheapest: one batch at a time, no spill, memory flat across
the whole stream. The partition tier spills the input once and processes it one
group or shard at a time, so its peak is set by the largest group rather than the
whole layer. The all-to-all tier is inherently resident because the operation is
global, and the overlay bounds it by tiling overlap clusters. Operations that are
global by nature and not tiled, such as Voronoi diagrams, Delaunay triangulation,
convex hulls, and neighbour graphs, are not streamed at all: `collect()` the
layer and run sf or terra on it.

A second axis is whether a step runs natively or through sf. The recognised
predicates and operations run in C straight off the WKB column when the data is
planar, which means projected, or geographic with spherical geometry off, or of
unknown CRS. Geographic data with `sf::sf_use_s2()` on keeps the spherical sf
path so the answer matches sf exactly. The native path parses the resident layer
once into a spatial index and tests each batch in C with no decode; the fallback
decodes each batch to `sf`. For the large planar workloads vectra targets, the
native path is the common case, which is why the examples project up front.

Four options tune the machinery:

- `vectra.spatial_flush` (default 500,000): rows buffered before a spill flush.

- `vectra.partition_budget`: rows held while routing a dissolve or partitioned
  join before flushing shards.
- `vectra.overlay_mem_limit` (default 2 GB): the overlay's peak working-set
  budget.
- `vectra.overlay_threads` and `vectra.spatial_threads`: worker counts for the
  overlay and the per-batch GEOS loops.

## Practical guidance

A few rules of thumb, with the numbers that drive them.

**Project before you stream.** A projected CRS puts the recognised verbs on the
native GEOS path and makes distances and areas Euclidean. Geographic data with
s2 on works but decodes every batch to sf and computes on the sphere, which is
correct but slower; reach for it only when spherical accuracy over large extents
matters more than throughput.

**Size the join grid to the data, not the machine.** For point-in-polygon
tagging set `grid(cellsize)` to anything larger than your polygons; the join is
exact for point left geometry at any such size. For an extended-on-extended join
make `cellsize` larger than the left features. For a partitioned nearest-feature
join make it at least the largest expected nearest distance, because the search
only reaches the eight neighbouring cells.

**Match `flush_rows` to feature weight.** The 500,000 default suits points and
small polygons. For heavy geometry, dense coastlines, detailed admin boundaries,
lower it so a spill batch is a workable size; for millions of tiny point rows,
raise it to cut the number of temporary files.

**Partition dissolve on a key whose groups fit.** Dissolve holds one group's
geometry while it unions. Dissolving a country into one outline is fine; grouping
by a key with a single enormous group is not, because that group must be resident
to union. Split the key finer if a group blows the budget.

**Trust the overlay's coverage warning.** If `spatial_overlay()` warns that
coverage drifted past `1e-4`, the named input features are finer than the
snapping grid or invalid after it. Pass a smaller `grid =` for fine geometry, or
inspect the `coverage_offenders` attribute on the result for the worst rows. A
clean run means the pieces reconstruct the inputs exactly.

When the verbs map to tiers and resident sets, the choice of tool is a memory
decision:

| Need | Verb | Holds resident |
|---|---|---|
| Per-feature transform | `spatial_map` | one batch |
| Keep / drop by location | `spatial_filter` | one batch + locator |
| Cut geometry to a mask | `spatial_clip` | one batch + mask |
| Tag with attributes | `spatial_join` | one batch + polygons |
| Tag, both sides huge | `spatial_join(partition=)` | one grid shard |
| Merge by group | `spatial_dissolve` | one group |
| Split self-overlaps | `spatial_overlay` | one overlap cluster |

**When not to stream.** For a layer that already fits in memory, sf is simpler
and faster: there is no gain in spilling 10,000 features through run-files when
`st_join()` runs in a blink. Streaming earns its overhead when the layer is
larger than RAM, or when a step in a longer pipeline would otherwise force a full
materialization. For genuinely global operations, Voronoi, Delaunay, hulls,
neighbour graphs, there is nothing to stream: collect the layer and run the
resident tool. The streaming envelope covers the operations whose memory can be
bounded, and the cost model names the ones it cannot.

```{r cleanup, include = FALSE}
unlink(c(f_poly, cent, fp, g_poly, fb, out))
```

## Where to go next

- [Out-of-core GIS with vectra](spatial.html) for the cost model and the raster
  half of the toolbox: `zonal()`, `focal()`, `terrain()`, `warp()`,
  `polygonize()`, `contours()`, `proximity()`.
- [Larger-than-RAM strategy](large-data.html) for the spill and partition tiers
  the spatial verbs reuse.
- [Species distribution modelling](sdm.html) for an end-to-end ecological
  workflow built on these verbs.
- The function reference for every spatial verb and its arguments.
```