Optional: guardrail + convenience for density (Mg/ha) rasters (mean already works; sum is the trap)

[cboettig]

Optional ergonomics — not a capability gap

(Reframed after discussion in data-workflows.) On reflection this is not a missing reducer: cng-datasets already handles per-area/density rasters correctly via --hex-resampling mean. Filing for the datasets team to assess two optional ergonomics improvements.

Background (data-workflows #202, resolved)

The Conservation International carbon rasters (Noon et al. 2022) store carbon as density, Mg C ha⁻¹, on a ~300 m WGS84 grid. They had been hexed with --hex-resampling sum, which produced a meaningless sum of densities ~6–7× below the true total (irrecoverable-2018: 21 Gt vs published 139 Gt).

The root cause was the reducer choice, not the tool. sum is only for values already integrated per pixel (e.g. GHS-POP persons/cell); a density is intensive and aggregates with mean, with the total being mean_density × cell_area.

What the tool already does right

--hex-resampling mean correctly yields the area-weighted mean density per H3 cell for a density raster. The correct total is then mean × H3 cell area (a query-time derivation, or bake it in downstream). No core capability is missing — we corrected all 10 carbon layers this way, and the decision rule ("read the source units; the reducer follows the units, not the conceptual quantity") is now documented in the data-workflows AGENTS.md.

Optional improvements for the datasets team to weigh

1. sum-on-density guardrail (low effort, high value): when --hex-resampling sum is used and the source band's units/metadata look per-area (ha, /ha, km-2, density, …), emit a warning — "summing a per-area density is usually wrong; did you mean mean (then × cell area for a total)?" This would have caught the carbon mistake automatically.
2. --density-to-total convenience (pure sugar): a flag emitting mean × exact H3 cell ground area, so a density→total product is one step and uses the exact per-cell area instead of a manual latitude approximation. Equivalent to today's mean + a downstream multiply.

Neither blocks anything — the correct path works today. The operator-decision safeguard lives in data-workflows AGENTS.md; this issue is only about optional tooling ergonomics.

Related: #95 (max/min — done), #88/#92 (seam geometry).

[cboettig]

Proposal: operator-declared value semantics → one-step accurate totals

I verified the mechanics against the actual irrecoverable_c_total_2018 COG and the exactextract build we ship. Two findings reshape the two optional items above.

Verified findings

1. The source COG carries no units metadata. gdalinfo on s3://public-carbon/v2/cogs/irrecoverable_c_total_2018.tif:

• Band 1 Type=Float32, Min=0 Max=895.89, NoData 3.4e38, grid 0.00269° (~300 m) global WGS84.
• band unit: None, band description: None. The only band metadata is RepresentationType=ATHEMATIC + statistics. No ha, /ha, density, or any units string anywhere.

➡️ Consequence for item 1 (units guardrail): it would not have caught the carbon mistake — there are no units to sniff, and real-world density COGs routinely omit them. A units heuristic is worth keeping as a cheap bonus, but it cannot be the primary safeguard.

2. exactextract (we ship 0.3.0) computes the exact area-integral in one op. On a constant-density test raster (value = 10 over a 1°×1° equatorial box):

op	result	meaning
sum (= coverage_weight=fraction)	100000	Σ(value·fraction) — the meaningless sum-of-densities
sum(coverage_weight=area_spherical_km2)	123914	Σ(value·pixel_area_km²) — the true integral
count(coverage_weight=area_spherical_km2)	12391.4 km²	confirms it = 10 × ground area

For a density in Mg ha⁻¹ the conserved total is sum(coverage_weight=area_spherical_km2) × 100 (1 km² = 100 ha). Unit factors: ha → ×100, km² → ×1, m² → ×1e6. This is exact per-pixel with true latitude-dependent area — and, crucially, exact at partial-coverage / nodata-edge cells too, where mean × full-cell-area overestimates by filling gaps with the mean density.

Proposed design — declared semantics, accurate by construction

Since the tool cannot infer density-vs-amount (finding 1), let the operator declare it on the invocation and make the accurate total the one-step result:

--value-semantics {amount|density}   # default amount (today's behavior; sum conserves)
--density-area-unit {ha|km2|m2}       # required denominator when density; default ha
--value-units "Mg/ha"                 # optional free-text, recorded to output + STAC for provenance

Behavior when --value-semantics density:

• reducer sum (or a clearer alias total) → internally sum(coverage_weight=area_spherical_km2) × unit_factor. The emitted column is the conserved total; SUM() over the hex layer is the true grand total — no downstream multiply, no latitude approximation.
• reducer mean → area-weighted mean density (unchanged), the right intensive value for choropleths; mean × h3_cell_area documented as an approximate total.

This subsumes item 2 (--density-to-total) but implements it via the exact integral rather than mean × cell_area, so coastal/fragmented layers don't bias high.

Guardrail, restated (item 1)

Because units are usually absent, the effective guardrail is: warn whenever sum is requested and --value-semantics was not declared —

sum conserves only for per-pixel amounts; if this raster is a per-area density, pass --value-semantics density (emits the area-integrated total) or mean (emits density).

Keep the units-string sniff as an additional trigger when metadata happens to exist. This catches the carbon class of error regardless of whether the COG carries units.

Scope / constraints

• Lives entirely in RasterProcessor + the exactextract call site (cog.py) and the CLI; small, localized.
• exact-extract method only — warp-centroid forwards to gdal.Warp and can't area-integrate; that combination must error.
• Tests: synthetic constant-density raster where Σ(value·area) is analytic — assert (a) the integral conserves the total, (b) two cells at different latitudes get correctly different areas (the property raw sum gets wrong), (c) the sum-without-declared-semantics warning fires.

Net: mean stays correct and nothing is forced, but a density raster can now be processed to the physically accurate total in a single declared invocation — the workflow we actually want for the carbon layers.

Aside for the data team: the published irrecoverable-carbon-2018 h9 hex carbon column has max ≈ 1504, above the source COG's pixel max (≈ 896), so its current units aren't self-describing — worth a separate look, independent of this tooling change.

[cboettig]

Refined CLI shape (after discussion)

Putting density on the reducer axis rather than behind a --value-semantics/--value-is-density flag — it reads more naturally and sits where operators already look:

--hex-resampling total --value-units "Mg/ha"

• total = area-integrated sum, implemented as sum(coverage_weight=area_spherical_km2) × unit_factor. The emitted column is the conserved grand total, so SUM() over the hex layer is the true total — exact per-cell area, exact at nodata/coastal edges.
• --value-units is required for total and must carry a per-area denominator. The denominator sets the area unit and the scale factor: /ha → ×100, /km2 → ×1, /m2 → ×1e6. Recognized per-area forms: /ha, ha-1, ha⁻¹, per ha (and km2/m2 equivalents). Units are also recorded to the output + STAC for provenance.
• Reducer vocabulary becomes sum | mean | max | min | mode | total.
• exact-extract method only — warp-centroid can't area-integrate, so total there errors.
• Guardrail retained: when sum is chosen but --value-units is per-area, warn — "summing a per-area density is usually wrong; did you mean total?"

So the carbon layers become a single physically-correct invocation:

cng-datasets workflow --dataset carbon/irrecoverable-2018 \
  --source-url .../irrecoverable_c_total_2018.tif \
  --hex-resampling total --value-units "Mg/ha" ...