From 4c8bd95d16d126c4431e196955213483bc398879 Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Fri, 2 Oct 2026 22:07:33 +0000
Subject: [PATCH 1/2] Add areal-weighted vector sum, median, and std ddof

- sum accepts vector sources. With params area = true each polygon's
  attribute is shared among the cells in proportion to the intersected
  area (TerraME's sum with area = true), conserving the total; without it
  every feature adds its whole value to each cell it touches. The column
  defaults to the target name (params column overrides). Vector sources
  are not clipped to the grid for sum, so a polygon crossing the border
  keeps its full denominator.
- median operator (fine-aligned, like std).
- std takes params ddof (0 default, 1 sample); defaults and spec_hash
  are unchanged.

Verified against TerraME on Itaituba: population within 4e-7 in all 620
cells, total conserved. Tests and docs updated.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 CHANGELOG.md                        |  10 ++
 README.md                           |   1 +
 disscube/operators/zonal.py         | 138 ++++++++++++++++++++++-
 docs/architecture/operators.md      |   7 +-
 docs/guides/pipeline_files.md       |   6 +-
 docs/terrame_fill_correspondence.md |   8 +-
 tests/test_sum_area_median.py       | 164 ++++++++++++++++++++++++++++
 7 files changed, 323 insertions(+), 11 deletions(-)
 create mode 100644 tests/test_sum_area_median.py

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 222120e..23e869e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,6 +9,16 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## [Unreleased]
 
 ### Added
+- **`sum` over vector sources**, with `params = {area = true}` for areal weighting: each polygon's
+  attribute (e.g. census population) is shared among the cells in proportion to the intersected
+  area, conserving the total. Reproduces TerraME's `sum` with `area = true` on Itaituba in all 620
+  cells (max error 4×10⁻⁷). Without `area`, every feature adds its whole value to each cell it
+  touches (a point to its own cell). `params = {column = …}` picks the column (default: the target
+  name). Vector sources are no longer clipped to the grid for `sum`, so a polygon crossing the border
+  keeps its full denominator.
+- **`median` operator**: per-cell median of the valid pixels, from the fine-aligned source.
+- **`std` `ddof` param** (`0` default, `1` for the sample estimate). The default is unchanged, so
+  existing `spec_hash`es hold.
 - **`type = "dem"` source**: elevation, or slope in degrees or percent, from SRTM, Copernicus GLO-30
   or TOPODATA (or your own `tiles`), read over the grid plus a `margin` and registered as a raster
   source with checksum and provenance. The slope is computed on a metric (UTM) grid by central
diff --git a/README.md b/README.md
index fc92ac6..f4f1307 100644
--- a/README.md
+++ b/README.md
@@ -180,6 +180,7 @@ for how DisSCube's operators relate to TerraME's *Fill*.
 | `mean` | zonal | average | no |
 | `sum` | zonal | sum | no |
 | `std` | zonal | nearest | no |
+| `median` | zonal | nearest¹ | no |
 | `min` | zonal | min | no |
 | `max` | zonal | max | no |
 | `majority` | zonal | nearest¹ | no |
diff --git a/disscube/operators/zonal.py b/disscube/operators/zonal.py
index 5e1f820..ed06fef 100644
--- a/disscube/operators/zonal.py
+++ b/disscube/operators/zonal.py
@@ -192,6 +192,7 @@ def _continuous_reduce(
     nodata: float | None,
     grid: GridSpec,
     stat: str,
+    ddof: int = 0,
 ) -> tuple[np.ndarray, np.ndarray]:
     """
     Reduce a fine continuous array into the target grid by real windows.
@@ -208,8 +209,11 @@ def _continuous_reduce(
         Sentinel marking invalid fine pixels.
     grid : GridSpec
         Target grid.
-    stat : {"std", "mean", "sum", "min", "max"}
+    stat : {"std", "mean", "sum", "min", "max", "median"}
         Statistic to compute per target cell over valid pixels.
+    ddof : int
+        Delta degrees of freedom for ``std`` (0 = population, 1 = sample). A
+        cell with ``n_valid <= ddof`` valid pixels gets NaN.
 
     Returns
     -------
@@ -244,7 +248,10 @@ def _continuous_reduce(
                 continue
             vals = block[mask]
             if stat == "std":
-                value[ti, tj] = float(np.std(vals))
+                if n_valid > ddof:
+                    value[ti, tj] = float(np.std(vals, ddof=ddof))
+            elif stat == "median":
+                value[ti, tj] = float(np.median(vals))
             elif stat == "mean":
                 value[ti, tj] = float(np.mean(vals))
             elif stat == "sum":
@@ -272,13 +279,107 @@ def compute(self, data, var: Variable, grid: GridSpec) -> xr.DataArray:
 
 
 class SumOperator(Operator):
+    """
+    Sum per cell.
+
+    *Raster source*: the sum of the pixels in the cell.
+
+    *Vector source* (TerraME's ``sum``): the numeric column named like the
+    target (or ``params = {column = ...}``) is summed over the features that
+    reach each cell.
+
+    - ``area = false`` (default): every feature adds its whole value to each
+      cell it touches (points: to the cell containing them).
+    - ``area = true``: areal weighting for polygons — each cell receives
+      ``value × area(cell ∩ polygon) / area(polygon)``, so the total is
+      conserved (census counts, herds, GDP). Polygon area is measured in the
+      grid CRS; use a projected grid for exact ratios. The source is not
+      clipped to the grid, so a polygon crossing the border keeps its full
+      denominator and only the part inside the grid is distributed.
+    """
+
     name = "sum"
     _resampling = Resampling.sum
+    clip_to_grid = False   # area weighting needs the whole polygon; rasters are unaffected
+    params: ClassVar[dict[str, str]] = {
+        "area": "vector polygons: distribute the value in proportion to the intersected area (default false)",
+        "column": "vector: numeric column to sum (default: the target name)",
+    }
 
     def compute(self, data, var: Variable, grid: GridSpec) -> xr.DataArray:
         if isinstance(data, xr.DataArray):
             return _passthrough(data)
-        raise TypeError(f"'sum' requires a raster source, got {type(data).__name__}")
+        if isinstance(data, gpd.GeoDataFrame):
+            return _wrap(_vector_sum(data, var, grid), grid)
+        raise TypeError(f"'sum' got unexpected type {type(data).__name__}")
+
+
+def _vector_sum(data: gpd.GeoDataFrame, var: Variable, grid: GridSpec) -> np.ndarray:
+    import shapely
+
+    column = var.params.get("column", var.name)
+    if column not in data.columns:
+        raise ValueError(f"'sum' needs a numeric column {column!r} in the vector source; it has {sorted(data.columns)}")
+    area = bool(var.params.get("area", False))
+    values = np.asarray(data[column], dtype=np.float64)
+    ok = data.geometry.notna().to_numpy() & ~data.geometry.is_empty.to_numpy() & np.isfinite(values)
+    geoms = shapely.make_valid(np.asarray(data.geometry.values[ok], dtype=object))
+    values = values[ok]
+    res = grid.resolution
+    total = np.zeros(grid.rows * grid.cols)
+    kinds = shapely.get_type_id(geoms)
+
+    if area:
+        # Only polygons have an area to share; ignoring others silently would lose mass.
+        if np.any(~np.isin(kinds, (3, 6))):
+            raise ValueError("'sum' with area = true needs polygon features only")
+        areas = shapely.area(geoms)
+        keep = areas > 0
+        geoms, values, areas = geoms[keep], values[keep], areas[keep]
+        pieces, owner = [], []
+        for k, g in enumerate(geoms):
+            for part in shapely.get_parts(g):
+                for q in _quarters(part, res):
+                    pieces.append(q)
+                    owner.append(k)
+        if not pieces:
+            return total.reshape(grid.rows, grid.cols)
+        pieces_arr = np.empty(len(pieces), dtype=object)
+        pieces_arr[:] = pieces
+        owner_arr = np.asarray(owner)
+        cells = _cell_boxes(grid)
+        tree = shapely.STRtree(cells)
+        piece_idx, cell_idx = tree.query(pieces_arr, predicate="intersects")
+        inter = shapely.area(shapely.intersection(pieces_arr[piece_idx], cells[cell_idx]))
+        k = owner_arr[piece_idx]
+        total += np.bincount(cell_idx, weights=values[k] * inter / areas[k], minlength=total.size)
+        return total.reshape(grid.rows, grid.cols)
+
+    # area = false: whole value to every cell reached
+    is_point = kinds == 0
+    if is_point.any():
+        px, py = shapely.get_x(geoms[is_point]), shapely.get_y(geoms[is_point])
+        col = np.floor((px - grid.bbox[0]) / res).astype(int)
+        row = np.floor((grid.bbox[3] - py) / res).astype(int)
+        inb = (col >= 0) & (col < grid.cols) & (row >= 0) & (row < grid.rows)
+        total += np.bincount(row[inb] * grid.cols + col[inb], weights=values[is_point][inb], minlength=total.size)
+    rest = ~is_point
+    if rest.any():
+        cells = _cell_boxes(grid)
+        tree = shapely.STRtree(cells)
+        g_idx, cell_idx = tree.query(geoms[rest], predicate="intersects")
+        total += np.bincount(cell_idx, weights=values[rest][g_idx], minlength=total.size)
+    return total.reshape(grid.rows, grid.cols)
+
+
+def _cell_boxes(grid: GridSpec) -> np.ndarray:
+    """The grid cells as boxes, row-major (index = row * cols + col)."""
+    import shapely
+
+    res = grid.resolution
+    xmin, ymax = np.meshgrid(grid.bbox[0] + np.arange(grid.cols) * res,
+                             grid.bbox[3] - np.arange(grid.rows) * res)
+    return shapely.box(xmin.ravel(), (ymax - res).ravel(), (xmin + res).ravel(), ymax.ravel())
 
 
 class StdOperator(Operator):
@@ -287,17 +388,44 @@ class StdOperator(Operator):
     # real per-cell window pass, so this operator uses fine alignment.
     _resampling = Resampling.nearest
     needs_fine_alignment = True
-    params: ClassVar[dict[str, str]] = {"subcells": "at most this many fine pixels per cell along each axis (memory bound)"}
+    params: ClassVar[dict[str, str]] = {
+        "subcells": "at most this many fine pixels per cell along each axis (memory bound)",
+        "ddof": "delta degrees of freedom: 0 (default, population) or 1 (sample standard deviation)",
+    }
 
     def compute(self, data, var: Variable, grid: GridSpec) -> xr.DataArray:
         if isinstance(data, xr.DataArray):
+            ddof = var.params.get("ddof", 0)
+            if ddof not in (0, 1):
+                raise ValueError(f"'std' param ddof must be 0 or 1, got {ddof!r}")
             fine, nodata = _fine_array(data)
-            value, cov = _continuous_reduce(fine, nodata, grid, "std")
+            value, cov = _continuous_reduce(fine, nodata, grid, "std", ddof=ddof)
             da = xr.DataArray(value, dims=("y", "x"), coords={"y": grid.ys, "x": grid.xs})
             return da.assign_coords(coverage_purity=(("y", "x"), cov))
         raise TypeError(f"'std' requires a raster source, got {type(data).__name__}")
 
 
+class MedianOperator(Operator):
+    """
+    Median of the valid pixels in each cell (robust to outliers and unmasked
+    clouds). Like ``std`` it needs the sub-cell pixels, so the source is
+    aligned to a fine grid first: the median of per-cell means would not be the
+    median of the pixels.
+    """
+    name = "median"
+    _resampling = Resampling.nearest
+    needs_fine_alignment = True
+    params: ClassVar[dict[str, str]] = {"subcells": "at most this many fine pixels per cell along each axis (memory bound)"}
+
+    def compute(self, data, var: Variable, grid: GridSpec) -> xr.DataArray:
+        if isinstance(data, xr.DataArray):
+            fine, nodata = _fine_array(data)
+            value, cov = _continuous_reduce(fine, nodata, grid, "median")
+            da = xr.DataArray(value, dims=("y", "x"), coords={"y": grid.ys, "x": grid.xs})
+            return da.assign_coords(coverage_purity=(("y", "x"), cov))
+        raise TypeError(f"'median' requires a raster source, got {type(data).__name__}")
+
+
 class MinOperator(Operator):
     name = "min"
     _resampling = Resampling.min
diff --git a/docs/architecture/operators.md b/docs/architecture/operators.md
index eb21c2c..78b7a80 100644
--- a/docs/architecture/operators.md
+++ b/docs/architecture/operators.md
@@ -49,12 +49,15 @@ These operators receive a high-resolution array snapped to the target grid origi
 
 | Operator | Raster | Vector | `requires_class_code` |
 |---|---|---|---|
-| `std` | true per-cell standard deviation | — | no |
+| `std` | true per-cell standard deviation (`ddof` 0 or 1) | — | no |
+| `median` | median of the valid pixels | — | no |
 | `majority` | dominant class by count | rasterizes with `class_code` (or 1) | no |
 | `minority` | least-frequent class by count | rasterizes with `class_code` (or 1) | no |
 | `percentage` | fraction of pixels of the target class | rasterizes with `class_code` | **yes** |
 
-Each accepts `params = {"subcells": n}`, a cap on the fine pixels per cell
+`sum` also accepts a vector source: `params = {area = true}` shares each polygon's value among the cells in proportion to the intersected area (total conserved); without it every feature adds its whole value to each cell it touches. The column is the one named like the target, or `params = {column = "…"}`.
+
+Each of the operators above accepts `params = {"subcells": n}`, a cap on the fine pixels per cell
 along each axis (see `GridAligner._align_fine`).
 
 The last three also produce `coverage_purity` and `dominance_purity` as coordinates of the output `DataArray` (persisted in the Zarr alongside the variable).
diff --git a/docs/guides/pipeline_files.md b/docs/guides/pipeline_files.md
index 890522f..27077b4 100644
--- a/docs/guides/pipeline_files.md
+++ b/docs/guides/pipeline_files.md
@@ -51,6 +51,7 @@ operator = "min_distance"
 | `distance` | `distance` (exact, from the cell centre) or `min_distance` (raster approximation) |
 | `area` | `area` (share of the cell covered by polygons) |
 | `presence`, `count`, `sum`, `minimum`, `maximum`, `stdev` | `presence`, `count`, `sum`, `min`, `max`, `std` |
+| `sum` with `area = true` | `sum` with `params = {area = true}` |
 
 How faithful each operator is to TerraME — and what is not supported yet — is
 measured in [TerraME Fill Cells Correspondence](../terrame_fill_correspondence.md).
@@ -145,7 +146,10 @@ of = ["state_paved", "federal_paved"]
 | Operator | `params` | Meaning |
 |---|---|---|
 | `distance` | `crs` | measure in this CRS — e.g. a projected one, for metres on a geographic grid |
-| `percentage`, `majority`, `minority`, `std` | `subcells` | at most this many fine pixels per cell along each axis: bounds the memory of a fine source over a large grid (a 100 m raster on a 1/12° grid would give ~90 × 90) |
+| `percentage`, `majority`, `minority`, `std`, `median` | `subcells` | at most this many fine pixels per cell along each axis: bounds the memory of a fine source over a large grid (a 100 m raster on a 1/12° grid would give ~90 × 90) |
+| `sum` | `area` | vector polygons: share the value in proportion to the intersected area (default `false`) |
+| `sum` | `column` | vector: numeric column to sum (default: the target name) |
+| `std` | `ddof` | `0` (default, population) or `1` (sample standard deviation) |
 
 An unknown key fails the plan. `fill = "nearest"` gives the cells an
 operator leaves without a value (NaN — e.g. a coastal cell a raster does not
diff --git a/docs/terrame_fill_correspondence.md b/docs/terrame_fill_correspondence.md
index 046f8db..56c0b76 100644
--- a/docs/terrame_fill_correspondence.md
+++ b/docs/terrame_fill_correspondence.md
@@ -53,10 +53,12 @@ and the parity cases in the recipes repository (`cases/terrame_fill`).
 | `distance` | `min_distance` | **approximation — semantics differ** | Rasterizes the features on the target grid and takes the Euclidean distance transform between cell centres (EDT × resolution). TerraME measures the distance from each cell polygon to the nearest feature, so `min_distance` overestimates it by up to about one cell (see the benchmarks). |
 | `average` / `mean` | `mean` | implemented; **parity verified** | Mean value per cell (continuous, area-weighted resampling). |
 | `sum` (raster) | `sum` | implemented | Sum per cell (continuous). |
-| `sum` with `area = true` (polygons) | — | **not implemented** | Distributes a polygon attribute (e.g. census population) over cells in proportion to the intersected area. `sum` accepts raster sources only. |
+| `sum` (vector, `area = false`) | `sum` | implemented | Adds the numeric column named like the target (or `params = {column = …}`) over the features that reach each cell: a point to the cell containing it, any other geometry to every cell it touches. |
+| `sum` with `area = true` (polygons) | `sum`, `params = {area = true}` | implemented; **parity verified** (Itaituba) | Distributes a polygon attribute (e.g. census population) over cells in proportion to the intersected area, so the total is conserved. Reproduces TerraME's `population` within 4×10⁻⁷ in all 620 cells. Polygons are not clipped to the grid: one that crosses the border distributes only its inside share. |
 | `minimum` | `min` | implemented; parity measured (Emas) | Minimum per cell. Matches TerraME in 99.0 % of cells: pixels that straddle a cell border count for both cells, while TerraME assigns each pixel to the cell containing its centre. |
 | `maximum` | `max` | implemented; parity measured (Emas) | Maximum per cell. Matches TerraME in 98.7 % of cells, for the same reason as `min`. |
-| `stdev` / `standardDeviation` | `std` | implemented (window-based) | True per-cell standard deviation over valid pixels. |
+| `stdev` / `standardDeviation` | `std` | implemented (window-based) | True per-cell standard deviation over valid pixels. Population (`ddof = 0`) by default; `params = {ddof = 1}` gives the sample estimate (NaN for a cell with a single valid pixel). |
+| `median` | `median` | implemented (window-based) | Median of the valid pixels per cell; not a TerraME fill strategy, offered for outlier-robust aggregation. |
 | `attribute` (value copy) | `attribute` | implemented (vector) | Rasterize a numeric vector column whose name matches the variable. |
 
 "Parity verified" means the operator reproduces TerraME's own output cell by
@@ -114,7 +116,7 @@ compared in percent (DisSCube fraction × 100).
 | `defor_255` | `percentage × coverage_purity` | 0.001 pp | 0.02 pp | **100 %** within 1 pp |
 | `distroad` — `distance` (lines) | `min_distance` | 1 782 m | 5 891 m | biased +1 782 m (r = 0.983) |
 | `distlocal` — `distance` (points) | `min_distance` | 2 499 m | 6 871 m | biased +2 497 m (r = 0.986) |
-| `population` — `sum`, `area = true` | — | — | — | not supported |
+| `population` — `sum`, `area = true` | `sum`, `area = true` | 0.000 | 0.000 (4×10⁻⁷) | **100 %** within 0.01; total 60 693 conserved |
 
 **Reading the results.**
 
diff --git a/tests/test_sum_area_median.py b/tests/test_sum_area_median.py
new file mode 100644
index 0000000..5950dba
--- /dev/null
+++ b/tests/test_sum_area_median.py
@@ -0,0 +1,164 @@
+"""
+Tests for ``sum`` over vectors (plain and areal-weighted), ``median`` and the
+``ddof`` parameter of ``std``.
+
+Areal weighting must conserve the polygon's value: the cells of a polygon that
+lies wholly inside the grid add up to its attribute, and a polygon crossing the
+border distributes only the share that falls inside.
+"""
+
+import geopandas as gpd
+import numpy as np
+import pytest
+import rasterio
+import rioxarray  # noqa: F401
+from rasterio.transform import from_bounds
+from shapely.geometry import LineString, Point, box
+
+from disscube import CubeClient, Derivation, GridSpec, SpatialSource
+from disscube.models import SpatialDerivation, Variable
+from disscube.models import SpatialSource as RasterSource
+from disscube.pipeline import PipelineContext
+from disscube.pipeline.aggregator import Aggregator
+from disscube.pipeline.aligner import GridAligner
+
+CRS = "EPSG:31982"
+GRID = GridSpec(id="g", type="local", crs=CRS, resolution=100, bbox=[0, 0, 400, 300])  # 3 rows x 4 cols
+
+
+@pytest.fixture
+def cube(tmp_path):
+    c = CubeClient(catalog=str(tmp_path / "c.db"), store=str(tmp_path / "s"))
+    c.register_grid(GRID)
+    return c
+
+
+def _vector(cube, tmp_path, sid, geoms, **columns):
+    path = tmp_path / f"{sid}.geojson"
+    gpd.GeoDataFrame(columns, geometry=geoms, crs=CRS).to_file(path, driver="GeoJSON")
+    cube.register_spatial_source(SpatialSource(id=sid, name=sid, format="vector", crs=CRS, asset_url=str(path)))
+
+
+def _sum(cube, target="pop", source="s", **params):
+    cube.derive_declarative(Derivation(target=target, source_id=source, operator="sum", params=params), grid_id="g")
+    return cube.load(target, grid_id="g").values
+
+
+# --------------------------------------------------------------------------- #
+# sum, area = true
+# --------------------------------------------------------------------------- #
+
+def test_area_sum_splits_a_polygon_in_proportion_to_the_area(cube, tmp_path):
+    # 200 x 100 polygon: covers cell (row 1, col 0) fully and (row 1, col 1) half
+    _vector(cube, tmp_path, "s", [box(0, 100, 150, 200)], pop=[300.0])
+    got = _sum(cube, area=True)
+    assert got[1, 0] == pytest.approx(200.0)       # 100 x 100 of 150 x 100
+    assert got[1, 1] == pytest.approx(100.0)       # 50 x 100 of 150 x 100
+    assert got.sum() == pytest.approx(300.0)
+
+
+def test_area_sum_conserves_the_total_over_many_polygons(cube, tmp_path):
+    polys = [box(10, 10, 390, 150), box(40, 160, 260, 290), box(100, 100, 300, 250)]  # the last overlaps both
+    _vector(cube, tmp_path, "s", polys, pop=[1000.0, 250.0, 40.0])
+    assert _sum(cube, area=True).sum() == pytest.approx(1290.0)
+
+
+def test_area_sum_polygon_crossing_the_border_keeps_its_full_denominator(cube, tmp_path):
+    # half of the polygon (x in 300..500) lies outside the grid (x <= 400)
+    _vector(cube, tmp_path, "s", [box(300, 0, 500, 100)], pop=[80.0])
+    got = _sum(cube, area=True)
+    assert got.sum() == pytest.approx(40.0)         # only the inside half is distributed
+    assert got[2, 3] == pytest.approx(40.0)
+
+
+def test_area_sum_with_a_custom_column(cube, tmp_path):
+    _vector(cube, tmp_path, "s", [box(0, 0, 100, 100)], pop=[5.0], herd=[7.0])
+    got = _sum(cube, target="cattle", area=True, column="herd")
+    assert got[2, 0] == pytest.approx(7.0)
+
+
+def test_area_sum_rejects_non_polygons(cube, tmp_path):
+    _vector(cube, tmp_path, "s", [LineString([(0, 0), (300, 0)])], pop=[1.0])
+    with pytest.raises(ValueError, match="polygon"):
+        _sum(cube, area=True)
+
+
+def test_vector_sum_without_the_column_is_an_error(cube, tmp_path):
+    _vector(cube, tmp_path, "s", [box(0, 0, 100, 100)], other=[1.0])
+    with pytest.raises(ValueError, match="pop"):
+        _sum(cube)
+
+
+# --------------------------------------------------------------------------- #
+# sum, area = false
+# --------------------------------------------------------------------------- #
+
+def test_plain_sum_adds_points_in_their_cell(cube, tmp_path):
+    _vector(cube, tmp_path, "s", [Point(50, 250), Point(60, 260), Point(350, 50)], pop=[1.0, 2.0, 10.0])
+    got = _sum(cube)
+    assert got[0, 0] == 3.0 and got[2, 3] == 10.0
+    assert got.sum() == 13.0
+
+
+def test_plain_sum_gives_a_polygon_to_every_cell_it_touches(cube, tmp_path):
+    _vector(cube, tmp_path, "s", [box(50, 50, 150, 150)], pop=[4.0])   # straddles 4 cells
+    got = _sum(cube)
+    assert got.sum() == 16.0 and (got == 4.0).sum() == 4
+
+
+# --------------------------------------------------------------------------- #
+# median and std ddof (raster, via aligner + aggregator)
+# --------------------------------------------------------------------------- #
+
+def _raster_run(tmp_path, array, var, nodata=None):
+    rows, cols = array.shape
+    path = tmp_path / "r.tif"
+    with rasterio.open(path, "w", driver="GTiff", height=rows, width=cols, count=1, dtype="float32",
+                       crs=CRS, transform=from_bounds(0, 0, 100, 100, cols, rows), nodata=nodata) as dst:
+        dst.write(array.astype(np.float32), 1)
+    grid = GridSpec(id="G1", type="local", crs=CRS, resolution=100, bbox=[0, 0, 100, 100])
+    src = RasterSource(id="S1", name="S1", format="raster", asset_url=str(path), crs=CRS)
+    ctx = PipelineContext(source=src, grid=grid, derivation=SpatialDerivation(
+        source_id="S1", grid_id="G1", role="t", variables=[var]))
+    GridAligner().execute(ctx)
+    Aggregator().execute(ctx)
+    return ctx.data[var.name]
+
+
+def test_median_odd_and_even(tmp_path):
+    odd = _raster_run(tmp_path, np.array([[1, 2, 100]], dtype=np.float32), Variable(name="m", operator="median"))
+    assert odd.values[0, 0] == pytest.approx(2.0)
+    even = _raster_run(tmp_path, np.array([[1, 2], [4, 100]], dtype=np.float32), Variable(name="m", operator="median"))
+    assert even.values[0, 0] == pytest.approx(3.0)
+
+
+def test_median_is_robust_to_an_outlier_where_the_mean_is_not(tmp_path):
+    a = np.array([[10, 11], [12, 5000]], dtype=np.float32)
+    med = _raster_run(tmp_path, a, Variable(name="m", operator="median")).values[0, 0]
+    mean = _raster_run(tmp_path, a, Variable(name="m", operator="mean")).values[0, 0]
+    assert med == pytest.approx(11.5) and mean > 1000
+
+
+def test_median_skips_nodata_and_reports_coverage(tmp_path):
+    out = _raster_run(tmp_path, np.array([[1, -1], [3, 9]], dtype=np.float32), Variable(name="m", operator="median"), nodata=-1)
+    assert out.values[0, 0] == pytest.approx(3.0)
+    assert out["coverage_purity"].values[0, 0] == pytest.approx(0.75)
+
+
+def test_std_ddof_default_is_population_and_1_is_sample(tmp_path):
+    a = np.array([[0, 2], [4, 6]], dtype=np.float32)
+    pop = _raster_run(tmp_path, a, Variable(name="s", operator="std")).values[0, 0]
+    smp = _raster_run(tmp_path, a, Variable(name="s", operator="std", params={"ddof": 1})).values[0, 0]
+    assert pop == pytest.approx(np.std([0, 2, 4, 6]))
+    assert smp == pytest.approx(np.std([0, 2, 4, 6], ddof=1))
+
+
+def test_std_sample_of_a_single_pixel_is_nan(tmp_path):
+    out = _raster_run(tmp_path, np.array([[5, -1], [-1, -1]], dtype=np.float32),
+                      Variable(name="s", operator="std", params={"ddof": 1}), nodata=-1)
+    assert np.isnan(out.values[0, 0])
+
+
+def test_std_rejects_a_bad_ddof(tmp_path):
+    with pytest.raises(ValueError, match="ddof"):
+        _raster_run(tmp_path, np.zeros((2, 2), dtype=np.float32), Variable(name="s", operator="std", params={"ddof": 2}))

From e6347677ac0702296b9f46d5cf6f5d9f87f2eb26 Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Fri, 2 Oct 2026 22:40:14 +0000
Subject: [PATCH 2/2] docs: document TerraME parity against the 2.0.1 goldens
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

The Fill benchmarks were described against a reference file that TerraME
2.0.1 does not reproduce. Against the goldens of LambdaGeo/luccme-goldens:

- coverage divides by the valid pixels, so `percentage` matches with no
  `coverage_purity` correction (Itaituba max 0.0064, Amazônia identical);
- distance is measured from the cell centre to the nearest vertex of the
  feature: `distance` is identical for points and smaller for lines where a
  line passes between vertices (mean 24 m on Itaituba, 300 m on Amazônia);
  `min_distance` errs by up to about a cell (1.2 km and 12 km mean).

Rewrite the correspondence tables, the Itaituba and Amazônia sections, the
known gaps and the positioning statement accordingly; drop the stale "sum with
area = true is not implemented" gap; fix the DistanceOperator docstring, which
claimed TerraME is smaller by half a cell diagonal.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 CHANGELOG.md                        |   9 +++
 disscube/operators/proximity.py     |   6 +-
 docs/terrame_fill_correspondence.md | 101 ++++++++++++++--------------
 3 files changed, 65 insertions(+), 51 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 23e869e..9308198 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -32,6 +32,15 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   (`date = "{year}-07-01"` with `years`); several public servers are tried with retries; HTTP 406
   points to `OSM_CONTACT`. Replaces the OpenStreetMap download script of the Lab15 reconstruction.
 
+### Changed
+- **TerraME parity documented against the TerraME 2.0.1 goldens** (`LambdaGeo/luccme-goldens`)
+  instead of an earlier reference of unrecorded provenance. TerraME 2.0.1 divides coverage by the
+  valid pixels, so `percentage` matches it with no `coverage_purity` correction (Itaituba max 0.0064,
+  Amazônia identical), and it measures distance from the cell centre to the nearest *vertex*:
+  `distance` is identical for points and smaller for lines where a line passes between vertices.
+  The notes saying TerraME measures from the cell polygon and divides by the whole cell were removed,
+  and so was the "half a cell diagonal" remark in the `DistanceOperator` docstring.
+
 ## [0.4.0] - 2026-09-30
 
 ### Changed
diff --git a/disscube/operators/proximity.py b/disscube/operators/proximity.py
index 7bd98c8..1f9cc2a 100644
--- a/disscube/operators/proximity.py
+++ b/disscube/operators/proximity.py
@@ -70,8 +70,10 @@ class DistanceOperator(Operator):
     geometry, and the source is not clipped to the grid, so features outside
     it — a town 50 km away, a road beyond the edge — count. Distances are in
     the grid CRS units (degrees on a geographic grid), measured from cell
-    centres; TerraME's ``distance`` fill measures from the cell polygon, so
-    it is smaller by up to half a cell diagonal.
+    centres. TerraME 2.0.1's ``distance`` fill measures from the cell centre to
+    the nearest *vertex* of the feature: identical to this for points, and
+    larger for lines wherever a line passes between vertices (this measures to
+    the segment).
 
     With ``params = {"crs": ...}`` the cell centres and the features are
     projected to that CRS and the distance is measured there — metres on a
diff --git a/docs/terrame_fill_correspondence.md b/docs/terrame_fill_correspondence.md
index 56c0b76..cc3a0e7 100644
--- a/docs/terrame_fill_correspondence.md
+++ b/docs/terrame_fill_correspondence.md
@@ -44,13 +44,13 @@ and the parity cases in the recipes repository (`cases/terrame_fill`).
 | TerraME fill strategy | DisSCube operator (`name`) | Status | Notes |
 |---|---|---|---|
 | `presence` | `presence` | implemented; parity measured (Emas) | Binary mask: 1 where any feature is present. Matches TerraME in 98.4–99.7 % of cells: lines are rasterized through cell centres, while TerraME marks every cell a line touches. |
-| `coverage` / `percentage` (raster) | `percentage` | implemented (window-based); **parity verified** (Itaituba, Amazônia) | Fraction (0..1) of the target class per cell, **over valid pixels**. TerraME divides by the whole cell instead; `percentage × coverage_purity` reproduces TerraME's value (see the benchmarks). Requires `class_code`. |
+| `coverage` / `percentage` (raster) | `percentage` | implemented (window-based); **parity verified** (Itaituba, Amazônia) | Fraction (0..1) of the target class per cell, **over valid pixels** — the denominator TerraME 2.0.1 uses too, so the value reproduces its goldens with no correction (Itaituba: max 0.0064; Amazônia: identical). `coverage_purity` stays available as metadata. Requires `class_code`. |
 | `area` (polygons) | `area` | implemented (exact); **parity verified** (Amazônia) | Fraction (0..1) of each cell covered by polygons (e.g. protected areas): intersection area / cell area, overlapping polygons counted once. Reproduces TerraME's `protected` within 0.01 in at least 99 % of the cells. |
 | `majority` / `mode` | `majority` | implemented (window-based) | Dominant class per cell; ties resolve to the smallest class value. |
 | `minority` | `minority` | implemented (window-based) | Least-frequent class per cell. |
 | `count` | `count` | implemented | Count of features per cell (proximity operator). |
-| `distance` | `distance` | implemented (exact, from the cell centre) | Euclidean distance from each cell centre to the nearest feature, in CRS units (or in the CRS given as `params = {crs = …}`, e.g. metres on a geographic grid), without clipping the source to the grid (features outside it count). TerraME measures from the cell polygon, so `distance` is larger by at most half a cell diagonal. The LuccME Lab15 cellular space was built with centre distances, and `distance` reproduces its fields. |
-| `distance` | `min_distance` | **approximation — semantics differ** | Rasterizes the features on the target grid and takes the Euclidean distance transform between cell centres (EDT × resolution). TerraME measures the distance from each cell polygon to the nearest feature, so `min_distance` overestimates it by up to about one cell (see the benchmarks). |
+| `distance` | `distance` | implemented (exact, from the cell centre); **parity verified** for points, lines differ by the vertex rule | Euclidean distance from each cell centre to the nearest feature, in CRS units (or in the CRS given as `params = {crs = …}`, e.g. metres on a geographic grid), without clipping the source to the grid (features outside it count). TerraME 2.0.1 measures from the cell centre to the nearest *vertex* of the feature, so for points the two are identical (error < 1 mm on Itaituba and Amazônia) and for lines `distance`, which measures to the segment, is smaller where a line passes between vertices (Itaituba: mean 24 m, 73.7 % of cells within 1 m; Amazônia: mean 300 m, 63.8 %). The LuccME Lab15 cellular space was built with centre distances, and `distance` reproduces its fields. |
+| `distance` | `min_distance` | **approximation — semantics differ** | Rasterizes the features on the target grid and takes the Euclidean distance transform between cell centres (EDT × resolution). TerraME measures from the cell centre to the nearest vertex, so `min_distance` differs from it by up to about one cell: against the 2.0.1 goldens the mean error is 1.2 km on Itaituba (5 km cells) and 11–12 km on Amazônia (50 km cells), biased low because `all_touched` marks the cells a line only touches. Prefer `distance`. |
 | `average` / `mean` | `mean` | implemented; **parity verified** | Mean value per cell (continuous, area-weighted resampling). |
 | `sum` (raster) | `sum` | implemented | Sum per cell (continuous). |
 | `sum` (vector, `area = false`) | `sum` | implemented | Adds the numeric column named like the target (or `params = {column = …}`) over the features that reach each cell: a point to the cell containing it, any other geometry to every cell it touches. |
@@ -81,9 +81,20 @@ design that have not yet been measured against TerraME.
 
 TerraME's `gis` package ships three *Fill* examples — Itaituba (the
 [Fill tutorial](https://github.com/TerraME/terrame/wiki/Fill)), Emas and
-Amazônia — each with its input layers, the Lua script and the cellular space
-TerraME produced. Those outputs are the reference: the same inputs are derived
-with DisSCube on the same grid and compared cell by cell.
+Amazônia — each with its input layers and the Lua script. The reference is the
+cellular space that TerraME 2.0.1 produces from them, archived as the goldens of
+[`LambdaGeo/luccme-goldens`](https://github.com/LambdaGeo/luccme-goldens)
+(v1.0.0, DOI [10.5281/zenodo.23107748](https://doi.org/10.5281/zenodo.23107748)):
+the same inputs are derived with DisSCube on the same grid and compared cell by
+cell.
+
+Earlier versions of this page compared against a different reference file, whose
+provenance was not recorded and which TerraME 2.0.1 does not reproduce (coverage as
+a percentage divided by the whole cell, distance 0 in every cell that contains a
+feature). The explanations drawn from it — that TerraME measures distance from the
+cell polygon and divides coverage by the whole cell — do not hold for the goldens
+and were removed. `elevation`, `population` and all of Emas are identical in both
+files.
 
 The full parity benchmark suite is maintained and executed on CI in the
 [DisSCube Case Studies and Recipes](https://github.com/LambdaGeo/disscube-recipes)
@@ -109,13 +120,13 @@ compared in percent (DisSCube fraction × 100).
 | TerraME fill | DisSCube | Mean abs. error | Max abs. error | Cells within tolerance |
 |---|---|---|---|---|
 | `elevation` — `average` (923 m raster) | `mean` | 0.89 m | 10.13 m | 72 % within 1 m (r = 0.9995) |
-| `defor_7` — `coverage` (60 m raster) | `percentage` | 2.17 pp | 50.13 pp | 92 % within 1 pp |
-| `defor_7` | `percentage × coverage_purity` | 0.09 pp | 0.64 pp | **100 %** within 1 pp |
-| `defor_87` | `percentage × coverage_purity` | 0.06 pp | 0.64 pp | **100 %** within 1 pp |
-| `defor_167` | `percentage × coverage_purity` | 0.01 pp | 0.39 pp | **100 %** within 1 pp |
-| `defor_255` | `percentage × coverage_purity` | 0.001 pp | 0.02 pp | **100 %** within 1 pp |
-| `distroad` — `distance` (lines) | `min_distance` | 1 782 m | 5 891 m | biased +1 782 m (r = 0.983) |
-| `distlocal` — `distance` (points) | `min_distance` | 2 499 m | 6 871 m | biased +2 497 m (r = 0.986) |
+| `defor_7` — `coverage` (60 m raster) | `percentage` | 0.054 pp | 0.64 pp | **100 %** within 1 pp |
+| `defor_87` | `percentage` | 0.054 pp | 0.64 pp | **100 %** within 1 pp |
+| `defor_167` | `percentage` | 0.006 pp | 0.39 pp | **100 %** within 1 pp |
+| `defor_255` | `percentage` | 0.000 pp | 0.002 pp | **100 %** within 1 pp |
+| `distlocal` — `distance` (points) | `distance` | 1×10⁻⁶ m | 5×10⁻⁶ m | **100 %** within 1 mm |
+| `distroad` — `distance` (lines) | `distance` | 24 m | 1 885 m | 73.7 % within 1 m; TerraME measures to the nearest vertex |
+| `distroad`, `distlocal` | `min_distance` | 1 180 m, 1 244 m | 3 449 m, 3 396 m | biased −552 m, −273 m (not recommended) |
 | `population` — `sum`, `area = true` | `sum`, `area = true` | 0.000 | 0.000 (4×10⁻⁷) | **100 %** within 0.01; total 60 693 conserved |
 
 **Reading the results.**
@@ -123,26 +134,25 @@ compared in percent (DisSCube fraction × 100).
 - **Continuous averages agree.** The residual in `elevation` comes from the
   averaging rule (area-weighted resampling vs. TerraME's per-pixel average) on
   a coarse 923 m source.
-- **Coverage agrees once the denominator is made explicit.** Raw
-  `percentage` differs only in the 50 border cells (last column and top row)
-  that the raster covers partially: TerraME divides the class area by the
-  *whole cell*, so its classes sum to less than 100 % there, while DisSCube
-  divides by the *valid pixels* and reports the covered share separately as
-  `coverage_purity`. Their product reproduces TerraME within 0.64 pp in every
-  cell. The difference is a design choice, not an error: DisSCube keeps "how
-  much of the cell is class *k*" apart from "how much of the cell has data",
-  and TerraME's value is recoverable exactly.
-- **Distances differ by construction.** Recomputing the exact distance from
-  each cell polygon to the nearest road reproduces `distroad` exactly (±0.5 m)
-  in 80 % of the cells (mean error 52 m), which identifies TerraME's semantics.
-  `min_distance` instead measures between rasterized cell centres, so at 5 km
-  it overestimates by about a third to a half of a cell on average.
+- **Coverage agrees with no correction.** TerraME 2.0.1 divides the class area
+  by the *valid pixels* of the cell, as DisSCube does; the residual is at most
+  0.64 pp. Multiplying by
+  `coverage_purity` — right for the earlier reference — would now be wrong.
+- **TerraME measures distance from the cell centre to the nearest vertex of the
+  feature.** For the localities (points) that is the feature itself, so
+  `distance` is identical. For roads TerraME ignores the middle of the
+  segments, so `distance` (to the segment) is smaller wherever a road passes
+  between vertices; recomputing the distance to the nearest vertex reproduces
+  `distroad` in 100 % of the cells, to the millimetre. The explanation is
+  inferred from the output; TerraME's source was not read.
+- **`min_distance` is an approximation.** It rasterizes the features with
+  `all_touched`, takes the Euclidean distance transform between cell centres
+  and multiplies by the resolution, so its error reaches about one cell.
 - **Rasters without a declared nodata.** The deforestation raster declares
   none, and its classes include 255. Earlier versions reprojected it with the
   `uint8` default fill value (255) and then treated that value as nodata,
-  dropping the legitimate class 255 (94 % of cells within 1 pp). Since the
-  fix, the fill value can no longer collide with the data, and `defor_255`
-  agrees with TerraME in every cell.
+  dropping the legitimate class 255. Since the fix, the fill value can no longer
+  collide with the data, and `defor_255` agrees with TerraME in every cell.
 
 ### Emas — 5 514 cells, 500 m
 
@@ -163,29 +173,22 @@ the 5 km PRODES raster, roads, ports and indigenous lands.
 
 | TerraME fill | DisSCube | Result |
 |---|---|---|
-| `prodes_10`, `prodes_208` — `coverage` | `percentage × coverage_purity` | **identical** in every cell with PRODES data; in the 55 cells without any, DisSCube reports NaN (purity 0) where TerraME reports 0 |
-| `distroads` — `distance` (lines) | `min_distance` | mean error 17 km; the exact polygon distance matches TerraME in 73 % of cells (83 % within 100 m) |
-| `distports` — `distance` (points) | `min_distance` | mean error 28 km; the exact polygon distance matches in 89 % of cells (91 % within 100 m) |
+| `prodes_10`, `prodes_208` — `coverage` | `percentage` | **identical** (max 5×10⁻¹¹) in every cell with PRODES data; in the 55 cells without any, DisSCube reports NaN where TerraME reports 0 |
+| `distports` — `distance` (points) | `distance` | **identical** (max 5×10⁻⁴ m) |
+| `distroads` — `distance` (lines) | `distance` | mean 300 m, max 18 km; 63.8 % of cells within 1 m. TerraME measures to the nearest vertex (reproduced in 100 % of cells) |
+| `distroads`, `distports` | `min_distance` | mean 12.0 km, 11.2 km; biased low (−6.8 km, −2.7 km); not recommended |
 | `protected` — `area` (polygons) | `area` | intersection area / cell area reproduces TerraME within 0.01 in ≥ 99 % of cells (`tests/test_terrame_parity.py`) |
 
-The exact polygon distance explains most but not all of TerraME's `distance`
-values on Amazônia (the largest residuals reach 13–16 km), so TerraME's rule
-needs to be pinned down before an exact operator is implemented.
-
 ## Known gaps relative to TerraME
 
-- **Exact vector distance.** TerraME's `distance` is measured from the cell
-  polygon to the nearest feature; `min_distance` is a raster approximation
-  between cell centres (see the benchmarks). An exact vector operator is
-  needed for parity; the Amazônia residuals show TerraME's rule is not only
-  the plain geometric distance.
+- **Distance to lines.** TerraME 2.0.1 measures from the cell centre to the
+  nearest *vertex* of a line; `distance` measures to the segment, which is the
+  geometrically correct value and is smaller where a line passes between
+  vertices. There is no vertex mode yet; for points the two coincide.
 - **Cell assignment of pixels and lines.** `min`/`max` count pixels that
   straddle a cell border, and `presence` rasterizes lines through cell
   centres; TerraME assigns each pixel to the cell containing its centre and
   marks every cell a line touches (Emas).
-- **Area-weighted vector aggregation.** TerraME's `sum` with `area = true`
-  has no DisSCube equivalent yet: `sum` accepts raster sources only. (`area`,
-  the fraction of the cell covered by polygons, is exact.)
 - **In-memory, single-tile by design.** The fine-alignment path materializes a
   fine array in memory; very large tiles at a high fine/target ratio are bounded
   by available memory — `params = {subcells = n}` caps the ratio.
@@ -199,7 +202,7 @@ needs to be pinned down before an exact operator is implemented.
 > over a catalogued data cube, with aggregation on windows aligned to the target
 > grid and explicit control of cell purity. On the three Fill examples shipped
 > with TerraME (Itaituba, Emas, Amazônia), raster averages and class coverage
-> reproduce TerraME's output cell by cell once cell purity is applied, and
-> `presence`, `min` and `max` agree in 98–99.7 % of cells, and polygon
-> `area` in ≥ 99 %; exact vector distance and area-weighted sums are the
-> remaining gaps.
+> reproduce TerraME's output cell by cell, distance to points is identical and
+> the areal-weighted `sum` conserves TerraME's totals, while `presence`, `min`
+> and `max` agree in 98–99.7 % of cells and polygon `area` in ≥ 99 %; distance
+> to lines follows TerraME's vertex rule only approximately.