feat(rust/sedona-raster-gdal): add RS_Clip - #1000
Conversation
|
@paleolimbot should we hoist the scalar geometry case? Is that a common case we want to optimize? |
|
I think it would be rather difficult to have enough raster tiles that further optimizing the geometry iteration would make a difference given that this is reaching all the way into GDAL and pulling pixels. For what it's worth, I believe the trivial version of this (north up) can be zero copy using views. (Whether that is the best choice for any particular chain of operations is a different story, so maybe not now). |
I was thinking this but never finished the views feature |
paleolimbot
left a comment
paleolimbot
left a comment
There was a problem hiding this comment.
Thank you for working on this!
I think this needs a Python integration test to check that this version of clip (and all its permutations of optional arguments) match what Rasterio would give on real and/or non-trivial input. Or if we are prioritizing PostGIS compatibility, building out the testing framework to be able to run tests against it.
I'm hesitant to go further without solving the nodata value rounding issues, which have also come up in other PRs and apply to all functions.
| // Decide array-vs-scalar over *all* args, not just the raster/geom | ||
| // the executor was given: a per-row band/option column over a scalar | ||
| // raster+geom must still yield an N-row array. (`executor.finish` | ||
| // only inspects its two exec args and would collapse to row 0.) | ||
| let out: ArrayRef = Arc::new(builder.finish()?); |
There was a problem hiding this comment.
Should you just put all the args into the exec and make the exec flexible enough to handle the arg order here?
Faithful copy of RS_Clip (src + bench) from the GDAL raster draft PR apache#704, wired into the module list, register, and bench manifest. This commit does NOT compile against current main — it predates several API changes; the reconciliation is isolated into the next commit so the diff there shows exactly what changed from Kristin's original.
Make Kristin's RS_Clip build and pass on today's main: - Add the sedona-proj dependency (CRS reprojection of the clip geometry), plus proj-sys in dev-dependencies for the transform tests. - Drop the removed RasterBandReader; read band bytes via BandRef::nd_buffer().as_contiguous(). - RasterMetadata width/height are now i64 (were u64). - Rewrite the tests onto the RasterSpec harness with self-contained synthetic rasters (no test4.tiff fixture) and RasterStructArray::try_new. Still 2-D only; N-D plane broadcasting follows in the next commit.
A clip is a 2-D (y, x) operation. For an N-D band (extra leading dims such as time), compute the mask and crop window once and apply them to every non-spatial plane, preserving the non-spatial dims and shrinking only the (y, x) extent. Emit the result through the N-D builder (start_raster_nd / start_band_nd); a plain 2-D raster is just the ["y", "x"] case. Adds an end-to-end test over a [time, y, x] raster.
- docs/reference/sql/rs_clip.qmd documenting all five overloads (optional all_touched / no_data_value / crop / lenient) with runnable examples. - Replace the placeholder bench with one that clips a generated raster by a reprojected polygon, over two polygon complexities, via the standard benchmark::scalar harness (needs sedona-testing's criterion feature).
- Output packaging now considers all args, not just [raster, geom]: a per-row band/option column over a scalar raster+geom yields the full N-row array instead of silently collapsing to row 0 (apache#1). - Band 0 (= all bands) reaches clip_raster instead of being clamped to 1 by band.max(1); a negative band errors rather than silently clipping band 1 (apache#2). - Default nodata is the band's own nodata, else the band data type's minimum, never a silent 0.0 that collides with real zero pixels (apache#3). - 'lenient' softens only the no-intersection case; genuine errors (malformed WKB, GDAL failures) always propagate. Dropped the eprintln! (apache#4). - Collapsed the three duplicate nodata/byte-size encoders onto gdal_common::nodata_f64_to_bytes and BandDataType::byte_size() (apache#6). - Dropped the redundant mask zero-init write (GDAL MEM bands are zero-filled) and the double rasterband(1) fetch (apache#8). Adds regression tests for the band-column row-drop, band-0 all-bands, negative-band error, and error-vs-lenient cases.
RS_Clip reads band pixels but wasn't tagged needs_pixels, so the planner never injected RS_EnsureLoaded ahead of it — an OutDb (e.g. Zarr/GeoTIFF-backed) raster would reach the clip with empty band data. Tag needs_pixels (materialize inputs first) and returns_bytes (output is already InDb, don't re-wrap).
…g + real bench - NULL band (or all_touched/crop/lenient) now yields a NULL row instead of silently clipping all bands; no_data_value NULL stays the 'unset' sentinel - reject a no_data_value that can't be represented in the band type instead of saturating it (e.g. -9999 -> 0 on UInt8) - guard the trailing (y, x) dims with is_spatial_dim_pair before clipping - preserve source band names on the clipped output - borrow the band buffer instead of copying it; reserve the plane accumulator - rebuild the benchmark so the raster covers the clip polygon's extent (the clip path is now actually exercised) and drop the copied CRS-tag shim - export rs_clip_udf from lib.rs; document band 0 = all bands and NULL handling
…ll_touched When a strict (lenient=false) clip selects no pixels, the message now depends on all_touched: a plain 'do not intersect' when every touched pixel was already considered (all_touched=true), and a sub-pixel hint pointing at all_touched otherwise. Documents that pixel-center selection can miss sub-pixel geometries.
apply_mask_and_crop copied one pixel at a time via a dynamic-width copy_from_slice, which the compiler can't vectorize. Copy each contiguous crop-window row in one bulk memcpy, then overwrite only the masked-out pixels with nodata (matching the pattern the crop=false path already uses). Add a large-clip benchmark config where this copy dominates.
…omputed up front The crop window is now the geometry envelope intersected with the raster extent, snapped outward to the pixel grid (PostGIS ST_Clip / gdalwarp -crop_to_cutline semantics), instead of the tight bbox of selected mask pixels. Computing the window before rasterizing rejects disjoint geometries with a cheap envelope check, shrinks the mask to the window, and removes the full-raster crop-window scan.
…s; shared executor helpers Review follow-ups: - nodata resolution works in the band's own byte representation end to end: an explicit no_data_value goes through the validating sedona_raster::traits::nodata_f64_to_bytes (fractional/out-of-range/ beyond-2^53 values error instead of silently saturating), the band's own nodata bytes are used verbatim (no lossy f64 round-trip for Int64/UInt64), and the default sentinel comes from the new BandDataType::min_value_le_bytes - RasterExecutor grows num_iterations_over / finish_over so kernels that project a subset of their arguments into the executor no longer hand-roll the row count and array-vs-scalar decision - the CRS engine comes from SedonaOptions when config options are available (crs_utils::with_crs_engine, falling back to the global PROJ engine) - docs page marked experimental; tests use sedona-testing's make_wkb instead of a GDAL WKT round-trip helper
raster_parity.py: numpy array + GDAL geotransform -> CRS-less GeoTIFF via rasterio -> RS_FromPath on the sedonadb side (exercising the needs_pixels -> RS_EnsureLoaded planner path); the reference is composed from rasterio primitives (geometry_window + geometry_mask) rather than rasterio.mask.mask, which additionally remaps source pixels valued at the band nodata. Option permutations run as rows of one query so the kernel executes its real array path instead of constant-folding literals. test_rs_clip.py: permutation matrix over band/all_touched/nodata/crop x two geometries chosen so every axis discriminates, geometry-shape sweep, default-sentinel sweep including int8/int64/uint64 (exact, no f64 round-trip), signature-defaults equivalence, empty-mask NULL semantics, and error pathways.
# Conflicts: # rust/sedona-raster-functions/src/crs_utils.rs # rust/sedona-raster-gdal/Cargo.toml # rust/sedona-raster-gdal/src/lib.rs # rust/sedona-raster-gdal/src/register.rs
Band.nodata unpacks the sentinel in the band's dtype (exact across the full Int64/UInt64 range). Raster.to_numpy() stacks bands into a (band, ...) array — zero-copy for single-band rasters, a documented copy for multiband since each band owns its own buffer; mixed-dtype rasters error rather than silently promoting.
sedonadb.raster_testing is the raster sibling of testing.py's DBEngine. Raster parity is operation-level rather than shared-SQL (rasterio is not a SQL engine, and raster SQL dialects disagree on names and argument order): RasterEngine exposes one method per operation returning decoded ClipResults, implemented by SedonaDB (the dialect under test; batches option combos as table rows so the kernel runs its array path) and Rasterio (a reference composed from geometry_window + geometry_mask, deliberately not mask.mask, which remaps source-nodata-valued pixels). Dialect engines compare strictly; reference engines take the resolved nodata fill from the caller. Replaces the test-local raster_parity helpers; RS_Clip tests now invoke through the dataframe API instead of f-string SQL.
…cessor The rst.clip() accessor with named parameters is code-generated from rs_clip.qmd's frontmatter at sedonadb-expr build time, replacing the f-string SQL + CAST(NULL AS DOUBLE) construction.
|
Heads up: I've pushed the review follow-ups (Band.nodata / Raster.to_numpy accessors, the sedonadb.raster_testing engine harness, expr-accessor tests) but GitHub is being laggy syncing the branch into the PR — jw/rs-clip is at 171fd09 if the Commits/Files tabs look stale. |
|
merging to continue working on the parity harness, happy to take more feedback into a new PR. |
paleolimbot
left a comment
paleolimbot
left a comment
There was a problem hiding this comment.
Thank you! The parity work is important too...feel free to open issues for these or incorporate them in one of your other PRs.
| /// Run `f` with the session's CRS engine: the [`SedonaOptions`] runtime engine | ||
| /// when config options are available (the query path), falling back to the | ||
| /// process-global PROJ engine otherwise (e.g. direct `invoke_batch` calls). | ||
| pub fn with_crs_engine<T>( |
There was a problem hiding this comment.
I assume the direct calls are coming from the SedonaUDF tester or tests generally, which would use a default config options. We should ensure our testing config options have SedonaOptions in them (I think this is easy) and perhaps a tester.with_crs_engine() so that UDFs that use that particular aspect can add sedona-proj as a dev-dependencies and populate it for their tests.
| // noDataValue at index 4 (when arg_count >= 5) | ||
| let nodata_array = if self.arg_count >= 5 { | ||
| args[4] | ||
| .clone() | ||
| .cast_to(&arrow_schema::DataType::Float64, None)? | ||
| .into_array(num_iterations)? | ||
| } else { | ||
| ScalarValue::Float64(None).to_array_of_size(num_iterations)? | ||
| }; | ||
| let nodata_array = as_float64_array(&nodata_array)?.clone(); |
There was a problem hiding this comment.
We can file a follow-up to fix this cast, which will loose precision for i64 and u64 nodata values and is handled differently in a few PRs. The solution is to resolve the arraydata and add add an accessor for the appropriate slice of the data buffer (which will be lossless little endian bytes regardless of type), and use a NullsBuffer to check for nulls.
| for (x, y) in corners { | ||
| let (col, row) = inverse.apply(x, y); | ||
| min_col = min_col.min(col); | ||
| max_col = max_col.max(col); | ||
| min_row = min_row.min(row); | ||
| max_row = max_row.max(row); | ||
| } |
There was a problem hiding this comment.
Are you sure this logic survives a rotated and/or skewed raster? Is there a test that confirms this works?
| [new_ulx, src[1], src[2], new_uly, src[4], src[5]] | ||
| } else { | ||
| [src[0], src[1], src[2], src[3], src[4], src[5]] | ||
| }; |
There was a problem hiding this comment.
I'm still not sure how this works with rotated input but if there's a test that confirms it works I trust you
| /// Masked/cropped bytes, plane-major in the band's dim order. | ||
| data: Vec<u8>, |
There was a problem hiding this comment.
Here, and in a few other places, there's an allocation of a Vec<u8> of the size of the input that isn't reused. If possible we should allocate one Vec<u8> as scratch space and reuse it for all clip operations in a single batch. I'm not sure if there's an opportunity to save the final output and wrap it as a view that we append to the builder...I didn't see that here but I may have missed it.
3 participants
Migrates RS_Clip out of the GDAL raster draft (#704) into its own module, with N-D plane-broadcast support, SQL reference docs, and a benchmark.
Also adds
sedonadb.raster_testing— a cross-engine raster parity harness (RasterEngine, implemented for SedonaDB and rasterio) used to validate RS_Clip against rasterio across the option permutations — plusBand.nodataandRaster.to_numpy()accessors.Deviations from Sedona Spark's RS_Clip
Checked against
RasterBandEditors.clipin apache/sedona master; recorded here so future contributors can tell which differences are intentional.Integer.MIN_VALUEfor every integral band type (silently truncated on narrower bands, e.g. byte) and Java'sDouble.MIN_VALUEfor floating bands — the smallest positive subnormal (4.9e-324), not the most-negative double. This PR resolves the fill as: explicit argument → the band's own nodata → the exact minimum of the band's data type, produced in band-native bytes so Int64/UInt64 sentinels are exact. That matches PostGIS's precedence (band nodata first) and avoids the JavaMIN_VALUEquirks.(int) noDataValueplus the AWT sample write), so e.g. -9999 on a byte band wraps. Here a value that is not exactly representable in the band's data type is an error.lenient => false. In Spark,lenientonly guards the raster/geometry intersection check; a geometry that intersects the raster but selects no pixels (a sliver passing between pixel centers withallTouched => false) returns an all-nodata raster. Here it yields NULL under the defaultlenient => true, and when strict the error message points atallTouched => true.One behavior confirmed to match Spark, noted because rasterio differs: pixels inside the geometry are copied verbatim, including pixels whose value equals the source band's nodata — they are not remapped to the output nodata the way
rasterio.mask.maskdoes (Spark's copy loop writes rawRaster.getPixelvalues with no source-nodata check; the Python parity tests pin this deliberately).