Add the EDR /area endpoint (Polygon zonal reduction)#60
Conversation
emmanuelmathot
left a comment
emmanuelmathot
left a comment
There was a problem hiding this comment.
Reviewed end to end: 544 tests pass, ruff + mypy --strict clean, and the modeling is spec-correct (the CovJSON schema pins composite.coordinates to ["x","y"], so the divergence from a lat-first CRS's referencing order is required, and it's locked with a test). The pre-allocation cell ceiling closes the reprojection and all-rings bypass paths it claims to. Two documentation-level notes on the reduction's semantics, inline; neither blocks merge.
| "Reduce the dataset over the polygon `coords` (a WKT " | ||
| "`POLYGON((x y, ...))`) to a single value per band by `stat` " | ||
| "(default `mean`) and return a CoverageJSON Polygon coverage. By " | ||
| "default the polygon is interpreted in, and the output labeled " | ||
| "with, CRS84 (longitude/latitude); pass `crs` to override. A " | ||
| "polygon that selects no valid pixels (outside the dataset, or " | ||
| "all nodata) yields a `null` value rather than an error. Vertical " | ||
| "selection (a `z` level, or a 3-D `POLYGON Z`) is rejected: the " | ||
| "2-D raster backing has no vertical dimension to reduce over." |
There was a problem hiding this comment.
The reduction's spatial semantics aren't surfaced to users here. rio-tiler's feature() rasterizes the cutline with all_touched=True, so every pixel the polygon boundary touches is included at full weight, with no area weighting. For small polygons (a few pixels across), boundary over-inclusion can make the result diverge noticeably from an area-weighted zonal statistic (e.g. exactextract) or from rasterstats' default (all_touched=False).
This is a reasonable consequence of the dependency-free (no shapely) design goal, not a bug -- but a user comparing against another tool would be surprised. Suggest one line in this description (and the README) noting the reduction is an unweighted, all-touched pixel statistic, not area-weighted.
| case Stat.MEDIAN: | ||
| reduced = np.ma.median(flat, axis=1) # type: ignore[no-untyped-call] | ||
| case Stat.STD: | ||
| reduced = np.ma.std(flat, axis=1) |
There was a problem hiding this comment.
np.ma.std defaults to ddof=0, i.e. this is the population standard deviation. The endpoint advertises std with no qualifier, and a user might reasonably expect the sample standard deviation (ddof=1). Worth a word confirming the convention -- either in the Stat docstring or the /area stat parameter description -- so it's an explicit choice rather than an inherited numpy default.
|
@emmanuelmathot, I made the explicit wording changes you requested, so please take another look. |
Add `reduce_bands` (reduce a masked (bands, H, W) array to one scalar per band by a `Stat`) and the `Polygon`/`PolygonInput` variant with its `imagedata_to_polygon_input` converter, toward the /area endpoint. The converter takes each band's dtype from the reduced array, not the source raster: the statistic sets the range's value type (mean of an int raster is float, count is integer).
Rename the rio-tiler converters for the domain they produce, now that `ImageData` maps to two domains (Grid and Polygon) and the old per-source naming collides: imagedata_to_coverage_input -> imagedata_to_grid_input pointdata_to_coverage_input -> pointdata_to_point_input
Add the `PolygonInput` arm to `to_coverage` (composite polygon axis, scalar range per band), extend the `CoverageInput` union to match, and share the 0-D range builder as `_create_scalar_ranges`.
Parse a WKT POLYGON, clip the dataset with Reader.feature(), and reduce each band to one scalar by `stat` (default mean) as a Polygon coverage. A polygon selecting no valid pixels yields null; a degenerate (zero-area) or oversized polygon is a 400.
Exercise /area in the container smoke check and bring the README and dev app current for the three shipped endpoints (the README had only /bbox).
A projected-CRS polygon reprojects to a destination grid far larger than its source-pixel footprint (Web Mercator near the poles), so the source-grid ceiling under-counted the read and let a multi-GB allocation through. Measure the ceiling on the destination grid feature() produces, via a shared _output_grid_dimensions extracted from _resolve_grid_dimensions, and reject a degenerate polygon in the route before opening the dataset.
The Content-Crs header and exclude_none serialization were duplicated verbatim across the /bbox, /position, and /area handlers. Route them through one _covjson_response(coverage, label_crs) helper.
Rename _read_feature -> _read_polygon_image (named for its result and domain like its _read_bounded_image / _read_point siblings, not the rio-tiler method) and reduce_bands -> reduce_each_band (it reduces within each band, not the band axis).
Thread each band through the reduction's identity: the description
names the statistic ("<stat> of <band>") and count drops the unit as a
dimensionless pixel count.
The pre-read cell-count ceiling measured only the exterior ring's box, but rio-tiler's feature() bounds every ring. A permissively-accepted hole reaching past the exterior slipped an arbitrarily large read past the guard, allocating it before the post-read backstop could reject (a resource-exhaustion path). Polygon.bounds now spans all rings, so the read is bounded before allocation on the same-CRS and reprojected paths alike.
A present-but-empty ?stat= was rejected with a 400 ("Unsupported stat
''"), unlike the other selector guards, which treat a blank query value
as absent. area_stat now defaults a blank stat to the mean.
875a44f to
5f0907e
Compare
The reduction's conventions were unstated: feature() cuts with all_touched=True, so a grazed boundary pixel counts as much as one wholly inside, and np.ma.std is the population standard deviation. Document both on the endpoint, the stat parameter, and the README. Give Stat.STD the label "population standard deviation", the only description that travels with a saved coverage, as COUNT already does.
5f0907e to
3d73e7b
Compare
2 participants
Adds GET /area: clips a raster to a WKT POLYGON (holes supported) and reduces it to one scalar per band by a statistic (min, max, mean, median, std, sum, count; default mean), returning a CoverageJSON Polygon-domain coverage. An out-of-bounds or all-nodata polygon reduces to null (200). The reduction statistic is recorded in each band's description and unit. The native-resolution read is bounded before allocation by a cell-count ceiling measured on the destination grid across all polygon rings. Dependency-free (no shapely).
Closes #56