feat(rust/sedona-raster-gdal): RS_ZonalStats and RS_ZonalStatsAll - #1066
Conversation
… UDFs Compute summary statistics of the raster pixels covered by a zone geometry. RS_ZonalStats returns one statistic (count, sum, mean, median, mode, stddev, variance, min, max) as a Float64; RS_ZonalStatsAll returns them all as a struct. Optional settings (band, all_touched, exclude_nodata, lenient) travel as a single JSON options argument rather than a positional overload ladder. Nodata is compared in the band's own byte representation to avoid a lossy f64 round-trip, and the band value type is dispatched once outside the pixel loop. Includes a criterion benchmark over raster resolution and zone complexity.
Mark both functions experimental and give runnable examples over RS_Example().
Cross-check every statistic against a rasterio + numpy reference through the dataframe API, with SQL-text smokes and error-path coverage.
…ersection The no-intersection gate now tests a true geometry intersection between the zone and the raster footprint (its convex hull) rather than a zone-envelope intersected with the raster extent (a bounding-box overlap), matching Sedona Spark's rsIntersects gate. A zone whose bounding box overlaps the raster but whose geometry is disjoint is now a no-intersection case (NULL under lenient, an error under strict) instead of count 0. Reuses the RS_Intersects convex-hull predicate through a shared raster_intersects_geom_wkb helper in sedona-raster-functions.
…loads Replace the trailing JSON-options argument with Apache Sedona Spark's positional overload ladder (band, stat_type, all_touched, exclude_nodata, lenient) so Spark SQL tends to run unchanged. The band-less overloads still error on a multiband raster rather than defaulting to band 1. Drop the now-unused serde/serde_json dependencies.
f5071b7 to
1e0fbab
Compare
…een RS_Clip and RS_ZonalStats Extract the pixel-window addressing and geometry rasterization that RS_Clip and RS_ZonalStats each carried into a shared mask module: one PixelWindow, one envelope_window, one rasterize_geometry_mask. RS_Clip now builds its geotransform from raster.transform() instead of an inline coefficient array. Behavior is unchanged; the deliberately-different per-pixel consumers (apply_mask_to_band writes nodata, collect_masked_values reads f64) stay separate.
…, and error paths - Reuse crs_utils::align_wkb_to_crs in both kernels: the equal-CRS fast path now borrows the zone WKB instead of reprojecting and copying it every row, and a CRS on exactly one side still errors (same policy, shared wording). - Guard that the band's nodata sentinel width matches the dtype byte size before the byte-equality compare, so a malformed nodata errors instead of silently disabling nodata exclusion. - Return sedona_internal_err! instead of unwrap() when appending struct fields, so a future field-layout change surfaces an error rather than aborting a Python release build.
…nsistency - Rename the RS_ZonalStats / RS_ZonalStatsAll arguments zone -> roi and exclude_nodata -> exclude_no_data so the whole signature mirrors Sedona Spark's getZonalStats(raster, roi, band, statType, allTouched, excludeNoData, lenient) (band / stat_type / all_touched / lenient already matched). Updates the qmd frontmatter (the generated accessor API), kernels, docs, and tests. - compute_statistics: a NaN pixel now yields NaN for every statistic (numpy semantics) instead of finite min/max (f64::min/f64::max silently skip NaN) alongside NaN sum/mean. Adds a unit test.
paleolimbot
left a comment
There was a problem hiding this comment.
I think some thought about the corner cases of stats handling (and how this is done in non-Sedona Spark places) should be considered here before this merges, but in general this looks great!
| // Rasterize the roi into a window-sized 0/1 mask (moves `geometry`, whose | ||
| // only remaining use is the burn). | ||
| let mask = rasterize_geometry_mask(gdal, geometry, &transform, &window, params.all_touched)?; |
There was a problem hiding this comment.
Is this mask heap allocated / would it benefit from scratch space in a loop since it's temporary?
|
Maybe a helper class in a separate file like struct StatsEvaluator<T> {
values: Vec<T>,
values_sorted: bool,
count: i64,
current_sum: Option<T>,
}
impl <T> for StatsEvaluator<T> {
fn sum(&mut self) -> T {
if let Some(cached_sum) = self.current_sum {
cached_sum
} else { // if floating point, sort values first, then compute sum and cache it }
}
}...or something. That will let you benchmark and test it independently, perhaps against an existing library that has spent some time considering the corner cases of (instead of going on hunches about what will or won't be faster or correct). This approach would also let you be lazy about computing stats (instead of computing all of them for every case and discarding all except one in the case of RS_ZonalStats). You could have one lazily constructed cached StatsEvaluator for each RasterDataType (for most cases there will only be one data type per batch so this won't explode scratch space). |
Co-authored-by: Dewey Dunnington <dewey@dunnington.ca>
# Conflicts: # rust/sedona-raster-gdal/Cargo.toml # rust/sedona-raster-gdal/src/lib.rs # rust/sedona-raster-gdal/src/register.rs
… and output Address three RS_ZonalStats review comments, all behavior-preserving: - Mode: compute it from the already-sorted values in a single longest-run pass (tie -> larger) instead of a HashMap keyed on the value bit pattern. - Single statistic: split the shared masking/collection (collect_zonal_values) from the computation so RS_ZonalStats computes only the requested statistic (compute_single_statistic sorts only for median/mode), while RS_ZonalStatsAll still computes every statistic. - Struct output: assemble the RS_ZonalStatsAll struct from one typed builder per field plus an outer validity buffer (ZonalStatsBuilders), building the StructArray once at the end rather than downcasting a StructBuilder's field builders on every row. The sample (two-pass) variance, median, and NaN/empty semantics are unchanged, factored into shared helpers. A test pins that the single-stat path returns exactly what the full computation does.
…t numpy Add rasterio+numpy parity coverage for the NaN and infinity edge cases Dewey's review flagged. A float64 fixture plants a NaN (then +inf) pixel inside the roi (not the nodata sentinel, so it is not excluded), and the resulting statistics are compared field-by-field against the numpy reference over the same masked selection: - NaN poisons every statistic except count (numpy semantics); - +inf flows through — sum/mean/max/mode become +inf, min/median stay finite, variance/stddev become NaN. A NaN/inf-aware comparator handles the fact that NaN never equals itself. This validates the existing behavior against a trusted reference; no kernel change.
# Conflicts: # python/sedonadb/tests/functions/test_rs_zonalstats.py
paleolimbot
left a comment
There was a problem hiding this comment.
Thank you!
Float64 as a working type works for me (it's what GDAL does), although it would be good to check at least one adversarial example where the input is all INT64_MAX (if there isn't one already that I missed) to make sure we don't get a panic anywhere. I suggested a minor tweak to avoid a few passes over the data that shouldn't be too hard.
…ean for min/max/variance The all-statistics path takes min and max from the single sort it already performs and reuses the computed mean for the variance pass, replacing two folds and a redundant mean computation. The single-statistic path keeps its O(n) folds so the stats that need no ordering do not pay for a sort.
Experimental
RS_ZonalStatsandRS_ZonalStatsAll(ported from #704), matching Sedona Spark's positional overload ladder.