james-willis commented on code in PR #1066:
URL: https://github.com/apache/sedona-db/pull/1066#discussion_r3670484904
##########
python/sedonadb/tests/functions/test_rs_zonalstats.py:
##########
@@ -15,195 +15,290 @@
# specific language governing permissions and limitations
# under the License.
-"""RS_ZonalStats parity against geometry_mask + numpy reductions.
-
-The rasterio comparator selects pixels with
-`rasterio.features.geometry_mask`, drops pixels valued at the band nodata,
-and reduces in float64. stddev/variance are the sample (ddof=1) statistics —
-that is what Sedona computes. The diagonal-edged zone under the centroid
-rule is on the Sedona Spark deviation ledger (its scanline rasterizer
-mis-places x-intercepts on non-square pixels and drops some center-inside
-pixels there, apache/sedona#3111). Zones that select no pixels are not
-compared here.
+"""RS_ZonalStats / RS_ZonalStatsAll cross-checked against a numpy reference.
+
+Both functions mirror Apache Sedona Spark's positional overloads, so the tests
+call them positionally: `(raster, roi, stat_type)` /
+`(raster, roi, band, stat_type[, all_touched[, exclude_no_data[, lenient]]])`
+for RS_ZonalStats and the same ladder without `stat_type` for RS_ZonalStatsAll.
+
+The fixture raster is CRS-less (so nothing reprojects and pixel selection is
+bit-comparable). The reference rasterizes the roi with `rasterio.features`
+(the same GDAL rasterizer the kernel uses) and reduces the selected pixels with
+numpy; exact-selection statistics (count, sum, min, max, median, mode) are
+compared exactly and the float-accumulation ones (mean, variance, stddev) with
+a tolerance.
+
+rasterio is required to write the fixture GeoTIFF, so the whole module skips
+when it is unavailable rather than importing it at module scope.
"""
+import numpy as np
import pyarrow as pa
import pytest
-from sedonadb.raster_testing import (
- Deviation,
- SedonaDB,
- SedonaSpark,
- expect_deviations,
- random_raster_data,
- write_geotiff,
-)
-
pytest.importorskip("rasterio")
-pytest.importorskip("shapely")
-pytestmark = pytest.mark.skipif(
- not SedonaDB.implements("zonal_stats"),
- reason="RS_ZonalStats is not implemented in SedonaDB (the parity subject)",
+from sedonadb.raster_testing import ( # noqa: E402
+ random_raster_data,
+ write_geotiff,
)
-# GDAL-order geotransform: origin (100, 500), 2-wide by 3-tall north-up
-# pixels; with a 7x6 raster the extent is x in [100, 114], y in [482, 500].
+# GDAL-order geotransform: origin (100, 500), 2-wide by 3-tall north-up pixels;
+# a 6x7 raster then spans x in [100, 114], y in [482, 500].
GDAL_TRANSFORM = (100.0, 2.0, 0.0, 500.0, 0.0, -3.0)
-HEIGHT, WIDTH = 6, 7
+BANDS, HEIGHT, WIDTH = 1, 6, 7
+NODATA = -9999
+
+# A rectangle well inside the raster that selects a block of pixels.
GEOM_RECT = (
"POLYGON ((102.6 495.8, 109.3 495.8, 109.3 485.9, 102.6 485.9, 102.6
495.8))"
)
-# Diagonal edges make all_touched matter, while staying clear of the corner
-# pixels where the fixture plants the dtype extremes (a float64 extreme in
-# the zone would push the squared-deviation statistics to infinity).
-GEOM_TRIANGLE = "POLYGON ((102.7 497.4, 112.4 496.9, 104.2 483.7, 102.7
497.4))"
-
-STATS = ["count", "sum", "mean", "min", "max", "stddev", "variance", "median"]
-
-DEVIATIONS = [
- Deviation(
- SedonaSpark,
- "zonal_stats",
- matches=lambda p: p.get("wkt") == GEOM_TRIANGLE and not
p.get("all_touched"),
- reason="Sedona's scanline rasterizer mis-places x-intercepts on "
- "non-square pixels and drops some center-inside pixels along "
- "diagonal edges; GDAL selects every center-inside pixel "
- "(https://github.com/apache/sedona/issues/3111)",
- ),
-]
+# Entirely outside the raster extent.
+GEOM_DISJOINT = "POLYGON ((900 900, 910 900, 910 890, 900 890, 900 900))"
+# Bounding box overlaps the raster, but the geometry itself is disjoint: the
+# triangle sits in the far corner of its bounding box, clear of the raster. A
+# bounding-box gate would burn no pixels and report count 0; a true-geometry
gate
+# (matching Sedona Spark's rsIntersects) treats it as a no-intersection case.
+GEOM_DISJOINT_BBOX = "POLYGON ((124 490, 124 510, 108 510, 124 490))"
+# A thin strip crossing the x = 104 pixel boundary but covering no pixel center
+# (centers sit at odd x): selects nothing unless all_touched.
+GEOM_SLIVER = "POLYGON ((103.6 499, 104.4 499, 104.4 483, 103.6 483, 103.6
499))"
+STATS = ["count", "sum", "mean", "median", "mode", "stddev", "variance",
"min", "max"]
+EXACT_STATS = {"count", "sum", "min", "max", "median", "mode"}
[email protected]("stat", STATS)
[email protected](
- ("wkt", "all_touched"),
- [
- (GEOM_RECT, False),
- (GEOM_RECT, True),
- (GEOM_TRIANGLE, False),
- (GEOM_TRIANGLE, True),
- ],
- ids=["rect-centroid", "rect-touched", "triangle-centroid",
"triangle-touched"],
-)
-def test_rs_zonalstats_matches_comparators(
- subject, comparator, request, tmp_path, wkt, all_touched, stat
-):
- """Every statistic over the float64 fixture, on both selection rules.
- The zone stays clear of the corners so the planted dtype extremes don't
- collapse sums to infinity."""
- expect_deviations(request, comparator, "zonal_stats", DEVIATIONS)
- tiff = tmp_path / "zonal.tif"
- write_geotiff(
- tiff,
- random_raster_data("float64", bands=2, height=HEIGHT, width=WIDTH),
- gdal_transform=GDAL_TRANSFORM,
- )
- got = subject.zonal_stats(tiff, wkt, band=2, stat=stat,
all_touched=all_touched)
- expected = comparator.zonal_stats(
- tiff, wkt, band=2, stat=stat, all_touched=all_touched
+def fixture_raster(tmp_path):
+ """A single-band int32 raster with planted nodata and a repeated value.
+
+ Returns `(path, band)` where `band` is the `(HEIGHT, WIDTH)` numpy array.
+ Two interior pixels hold the nodata value and three hold a repeated value
+ (66) so the mode is unambiguous and nodata exclusion is observable.
+ """
+ data = random_raster_data(
+ "int32",
+ bands=BANDS,
+ height=HEIGHT,
+ width=WIDTH,
+ seed=7,
+ plants={(1, 1): NODATA, (2, 2): NODATA, (1, 2): 66, (2, 3): 66, (3,
1): 66},
)
- # Engines reduce in different orders, so exact float equality is not
- # attainable; 1e-9 passes summation noise and still fails any semantic
- # mismatch (selection, nodata handling, ddof).
- assert got == pytest.approx(expected, rel=1e-9), (wkt, all_touched, stat)
-
-
[email protected]("stat", ["count", "sum"])
-def test_rs_zonalstats_excludes_nodata(subject, comparator, tmp_path, stat):
- """A pixel valued at the band nodata inside the zone is excluded from
- the reduction by every engine."""
- tiff = tmp_path / "zonal_nodata.tif"
- write_geotiff(
- tiff,
- random_raster_data(
- "uint8", bands=1, height=HEIGHT, width=WIDTH, plants={(2, 3): 200}
- ),
- gdal_transform=GDAL_TRANSFORM,
- nodata=200.0,
+ path = tmp_path / "zonal.tif"
+ write_geotiff(path, data, gdal_transform=GDAL_TRANSFORM, nodata=NODATA)
+ return path, data[0]
+
+
+def numpy_reference(band, wkt, *, all_touched, exclude_no_data):
+ """Reference statistics over the pixels the roi selects, via
rasterio+numpy.
+
+ Returns a dict of every statistic, or the sentinel string ``"empty"`` when
+ the selection is empty (the caller maps that to count 0 / NULLs).
+ """
+ import rasterio.features
+ import shapely
Review Comment:
done
--
This is an automated message from the Apache Git Service.
To respond to the message, please log on to GitHub and use the
URL above to go to the specific comment.
To unsubscribe, e-mail: [email protected]
For queries about this service, please contact Infrastructure at:
[email protected]
