Expose ground-source and ground-derivation controls on the point-cloud CHM request

[amarcozzi]

Summary

Follow-on to #330. The point-cloud CHM handler decides how to establish ground entirely on its own and offers the user no say in it. Give the create request two controls: which ground source to use, and the parameters of the derivation when one is derived.

Canopy height is height above ground, so the ground surface bounds the quality of everything downstream — the CHM, the tree heights ITD derives from it, and the canopy fuels built from those. Today a user who can see it went wrong (via source.ground.ground_coverage and max_ground_distance_m, which #330 added for exactly this reason) has no lever to do anything about it.

What the handler does today

services/griddle/griddle/handlers/chm_point_cloud.py decides by content, not provenance: if the point cloud document's summary.point_classes contains ASPRS class 2, those returns are the ground surface; otherwise the surface is derived by min-surface → PMF → snap-back. It records which happened as source.ground.ground_source (classification | derived).

The PMF parameters are module constants, stored in metres and converted to cells at run time:

PMF_MAX_WINDOW_M = 33.0
PMF_SLOPE = 1.0
PMF_INITIAL_DISTANCE_M = 0.15
PMF_MAX_DISTANCE_M = 2.5
GROUND_SNAP_TOLERANCE_M = 0.5

These are filters.pmf's published defaults, chosen deliberately — hand-picked values measured 4.8× worse than the published ones during #330 profiling.

Control 1 — which ground source

The automatic rule is right most of the time and wrong in one specific way: it trusts a vendor classification whenever one exists. Some acquisitions classify ground poorly, and a user who has compared the result against known terrain should be able to say "ignore the classification, derive it".

Mirror the response field rather than inventing a parallel vocabulary:

{
  "source_point_cloud_id": "{{POINT_CLOUD_ID}}",
  "ground_source": "derived"
}

• auto (default) — current behavior: classification when class 2 is present, derived otherwise.
• classification — require the cloud's class 2 returns; 422 if the cloud has none, rather than silently falling back.
• derived — always derive, ignoring any classification present.

The response's source.ground.ground_source keeps reporting what actually happened, so auto stays observable.

Note the original plan item called this use_3dep_ground_classification: bool. That name is wrong on both counts — the decision is not 3DEP-specific (it applies equally to uploads, by content), and a bool cannot express "require it, don't fall back". Use the enum.

Control 2 — the derivation parameters

Only meaningful when the surface is derived. The #330 profiling identified the controlling variable precisely:

The controlling variable is void width vs. the largest window, and every local-neighbourhood method fails identically.

So max_window_m is the field that actually matters — a user working a site with wide ground-return-free voids (large buildings, closed evergreen canopy, the failure cases we measured at 7.07 m RMSE on dense urban) wants to raise it. The rest are exposed for completeness and because reproducing a published result sometimes requires them.

{
  "source_point_cloud_id": "{{POINT_CLOUD_ID}}",
  "ground_source": "derived",
  "ground_derivation": {
    "max_window_m": 66.0,
    "slope": 1.0,
    "initial_distance_m": 0.15,
    "max_distance_m": 2.5,
    "snap_tolerance_m": 0.5
  }
}

Every field optional, defaulting to the constant above. Reject ground_derivation when ground_source is classification — it would have no effect, and silently ignoring it is worse than a 422.

Validator worth writing: the window is converted as window_cells = max(3, round(window_m / resolution)), so at a coarse alignment.resolution the requested window silently collapses toward the 3-cell floor. A 33 m window at 11 m cells is 3 cells; at 20 m cells it is still 3 cells but no longer means 33 m. Either reject the combination or warn — do not let the stored parameters imply something the run did not do.

Explicitly out of scope: the pmf | smrf algorithm choice

The original plan item proposed ground_classification_algorithm: Literal["pmf","smrf"] with per-algorithm config objects. Leave it out, because we measured it and there is no second algorithm worth offering.

CHM RMSE against vendor-classified ground, classification discarded, same eight clouds:

Method	RMSE
PDAL filters.pmf	0.190 m
ours, PMF + snap-back	0.238 m
ours, surface only	0.363 m
PDAL filters.smrf	1.147 m
ours at hand-picked slope=0.3	1.149 m

SMRF at its own defaults lost to PMF by ~5×, and #330 concluded replicating it was not worth doing. A discriminated union with one member is ceremony, not design.

When a second algorithm does earn its place, the shape to copy already exists in this codebase: StemIsolationAlgorithm in services/api/api/resources/inventories/tree/chm/schema.py is a Field(discriminator="name") union of StemIsolationLmf | StemIsolationVwf on the adjacent ITD path. Same problem, same resource family — mirror it rather than reinventing.

The more valuable follow-up is not a second local-neighbourhood algorithm at all. Per the #330 profiling, every such method fails the same way, and the real fix is a large-scale term — a coarse-to-fine surface whose coarsest window exceeds the widest expected void, or an object-extent test demoting large flat elevated components before the morphological pass. That deserves its own issue.

Reproducibility

Both controls must be persisted on the grid's source, resolved, the way alignment.resolution already is — a stored grid should record exactly what it was built with, not what was requested. That means writing back the effective PMF parameters even when the user supplied none.

Baseline the derivation against PDAL filters.pmf

The RMSE table above is a one-off measurement from #330 profiling. Nothing checked into
this repo reproduces it. The only PDAL code we have is
services/lakitu/scripts/profile_chm.py — untracked, lazily importing pdal inside a
function, and profiling memory and throughput, not accuracy. The accuracy harness that
produced 0.190 / 0.238 does not exist here.

That is tolerable while the PMF parameters are module constants: one configuration, one
measurement, taken once. This issue is what makes it intolerable. The moment
max_window_m, slope, initial_distance_m, and max_distance_m become user-settable,
"0.238 m RMSE" stops describing our implementation and starts describing one point in a
space we have never measured. We would be shipping knobs whose effect we cannot check.

The validator item under Control 2 is a claim about exactly that space —
window_cells = max(3, round(window_m / resolution)) collapses toward the 3-cell floor at
coarse resolution — and we currently have no way to measure what it costs. A quick sweep on
synthetic terrain (bare 900 m cone, no vegetation at all, unclassified so the derived path
runs; the correct CHM is 0.0 everywhere) shows the cost is real and non-monotonic:

alignment.resolution	window (cells)	window (m)	max false CHM, 9% slope	max false CHM, 60% slope
1 m	33	33	0.00	11.46
2 m	16	32	0.00	4.37
5 m	7	35	0.54	7.00
10 m	3	30	1.07	8.49
30 m	3	90	3.77	25.46

At 30 m cells the floor makes the window ~3× wider in metres than intended, and the
window <= 3 branch pins the tolerance at PMF_INITIAL_DISTANCE_M (0.15 m) instead of the
slope-derived value — so the wide opening is accepted almost unconditionally. _pmf's own
docstring says that tolerance "is what stops a wide window from flattening a real ridge";
the cell-count branch bypasses it at precisely the resolutions where the window is widest.
A bare summit reading as a 25 m crown is something ITD will happily detect.

These numbers are synthetic and indicative, not a validation. Getting real ones is the
point of this section.

What to measure

1. Agreement with filters.pmf at the published defaults. Our _pmf against PDAL's,
   same input surface, 1 m. Tightest possible test: same algorithm, same parameters, so any
   divergence is an implementation bug rather than a design difference. Compare the ground
   surface, not just the end-to-end CHM, so a discrepancy is localized to the filter.
2. Reproduce the RMSE table as an assertion. End-to-end CHM against vendor-classified
   ground on the eight-cloud set, with a stated tolerance. Turns 0.190 / 0.238 from folklore
   into a regression guard.
3. The resolution sweep, 1 → 30 m. PDAL's filters.pmf takes cell_size and
   max_window_size in cells directly, so it can be driven to the intended metre window at
   each resolution — which makes it the reference for what the 3-cell floor actually costs.
   This is the measurement that should decide whether the validator rejects the combination,
   warns, or whether the floor itself should be reworked.
4. A degenerate-terrain case with a known answer. Bare steep terrain, no vegetation, so
   the correct CHM is 0.0 by construction and no vendor reference is needed. Cheap,
   deterministic, and the one item here that can run on every PR.

Where it runs

PDAL is a conda-first native package and is not a dependency of any service. Adding it to
griddle's runtime image to satisfy a test would be backwards — the #330 conclusion was that
PDAL OOMs on large clouds and the streamed numpy path replaces it.

• Items 1–3 are opt-in: pytest.importorskip("pdal") plus a marker, and they need the
  eight-cloud fixture set, which is too large for CI. Run on demand.
• Item 4 needs neither PDAL nor fixtures and belongs in
  services/griddle/tests/handlers/test_chm_point_cloud.py, alongside the existing
  TestGroundGapFill resolution tests — which today assert the conversion arithmetic
  (_fill_cells) but never that the resulting surface is right.
• The accuracy harness itself should be recovered or rewritten and checked in next to
  profile_chm.py, which covers the orthogonal axis.

Ripples

• Breaking change: no. New optional request fields, new optional stored fields. Beta regardless.
• OpenAPI: field descriptions plus a json_schema_extra example showing the derived path. This is where users will actually meet these controls.
• FastFuels-Web: the CHM how-to filed as silvxlabs/FastFuels-Web#256 covers source.ground; it would gain the request-side controls. Worth a comment there once this lands rather than a second issue.
• Python SDK: fastfuels-sdk-python has no canopy grid support at all today, so this rides on that larger decision rather than being a separate SDK task.

Acceptance criteria

• ground_source on CreatePointCloudChmRequest as a typed enum, defaulting to auto
• classification returns 422 when the cloud carries no class 2 rather than falling back
• ground_derivation config with all fields optional and defaulted to the published PMF values
• 422 when ground_derivation is given with ground_source: "classification"
• The resolution/window-collapse interaction handled explicitly, not silently
• Effective values persisted on the grid source, including when defaulted
• Handler reads them instead of its module constants; constants become the defaults in one place
• Router tests for each rejection, against a live server
• OpenAPI descriptions and an example body for the derived path
• No pmf | smrf discriminator in this issue
• Accuracy harness checked in (the existing profile_chm.py covers memory/throughput only)
• _pmf agrees with filters.pmf at the published defaults on a common input surface at 1 m, within a stated tolerance
• The 0.190 / 0.238 m RMSE table reproduces as an assertion rather than prose
• Resolution sweep 1 → 30 m measured against a PDAL reference driven to the same metre window
• A no-PDAL, no-fixture degenerate-terrain test (bare steep terrain → CHM ≈ 0) in the griddle unit suite
• PDAL stays out of every service's runtime dependencies; the baseline tests skip cleanly when it is absent

Related

• Convert a point cloud to a CHM (grids/canopy point_cloud source) #330 — the CHM endpoint these controls attach to
• Investigate a spatially partitioned point-cloud artifact: chunked LAZ vs a custom COPC writer #483 — chunked point-cloud processing; the derived-ground path is where tile seams would show
• silvxlabs/FastFuels-Web#256 — the CHM how-to that documents source.ground