diff --git a/CHANGELOG.md b/CHANGELOG.md
index 222120e..9308198 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
@@ -22,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/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/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/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..cc3a0e7 100644
--- a/docs/terrame_fill_correspondence.md
+++ b/docs/terrame_fill_correspondence.md
@@ -44,19 +44,21 @@ 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` 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
@@ -79,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)
@@ -107,40 +120,39 @@ 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) |
-| `population` — `sum`, `area = true` | — | — | — | not supported |
+| `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.**
 
 - **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
 
@@ -161,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.
@@ -197,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.
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}))