feat(rust/sedona-raster-gdal): add RS_Clip

docs/reference/sql/rs_clip.qmd

@@ -0,0 +1,131 @@
+---
+# Licensed to the Apache Software Foundation (ASF) under one
+# or more contributor license agreements.  See the NOTICE file
+# distributed with this work for additional information
+# regarding copyright ownership.  The ASF licenses this file
+# to you under the Apache License, Version 2.0 (the
+# "License"); you may not use this file except in compliance
+# with the License.  You may obtain a copy of the License at
+#
+#   http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing,
+# software distributed under the License is distributed on an
+# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+# KIND, either express or implied.  See the License for the
+# specific language governing permissions and limitations
+# under the License.
+
+title: RS_Clip
+description: >
+  Clips a raster to a geometry, setting pixels outside the geometry to nodata
+  and (by default) cropping to the geometry's bounding box.
+kernels:
+  - returns: raster
+    args:
+    - raster
+    - {name: band, type: integer}
+    - {name: geom, type: geometry}
+  - returns: raster
+    args:
+    - raster
+    - {name: band, type: integer}
+    - {name: geom, type: geometry}
+    - {name: all_touched, type: boolean}
+  - returns: raster
+    args:
+    - raster
+    - {name: band, type: integer}
+    - {name: geom, type: geometry}
+    - {name: all_touched, type: boolean}
+    - {name: no_data_value, type: double}
+  - returns: raster
+    args:
+    - raster
+    - {name: band, type: integer}
+    - {name: geom, type: geometry}
+    - {name: all_touched, type: boolean}
+    - {name: no_data_value, type: double}
+    - {name: crop, type: boolean}
+  - returns: raster
+    args:
+    - raster
+    - name: band
+      type: integer
+      description: >
+        1-based index of the band to clip, or 0 to clip every band. A NULL
+        band yields a NULL result.
+    - name: geom
+      type: geometry
+      description: >
+        Clip geometry. Reprojected into the raster's CRS when both carry one;
+        it is an error for exactly one side to have a CRS.
+    - name: all_touched
+      type: boolean
+      description: >
+        If true, keep every pixel the geometry touches; otherwise keep only
+        pixels whose center falls inside it. Defaults to false.
+    - name: no_data_value
+      type: double
+      description: >
+        Value written to pixels outside the geometry. Defaults to the band's
+        own nodata value, or the band data type's minimum if it has none. It is
+        an error if the value does not fit the band's data type.
+    - name: crop
+      type: boolean
+      description: Crop the output to the geometry's bounding box. Defaults to true.
+    - name: lenient
+      type: boolean
+      description: >
+        If true, a geometry that does not intersect the raster yields NULL;
+        if false, it raises an error. Defaults to true.
+---
+
+## Description
+
+`RS_Clip` clips a raster to a geometry, in the spirit of PostGIS `ST_Clip`.
+Pixels of the selected band that fall outside the geometry are set to the
+band's nodata value (or `no_data_value` if given); pixels inside are kept. When
+`crop` is true (the default) the output is cropped to the geometry's bounding
+box intersected with the raster extent, snapped to the pixel grid; unselected
+pixels within that window remain as nodata padding. Otherwise the original
+extent is preserved.
+
+The geometry is reprojected into the raster's CRS when both carry one. If
+exactly one side has a CRS, the call errors; if neither does, the geometry is
+used as-is.
+
+A pixel is selected only when its center falls inside the geometry. A geometry
+smaller than a pixel, or a thin sliver, can therefore select no pixels even
+though it overlaps the raster; set `all_touched` to true to keep every pixel the
+geometry touches instead.
+
+The clip is a 2-D `(y, x)` operation. For an N-D raster (bands with extra
+dimensions such as `time`), the same mask and crop window are broadcast across
+every non-spatial plane, so the non-spatial dimensions are preserved and only
+the `(y, x)` extent changes.
+
+## Examples
+
+```sql
+SELECT RS_Width(
+  RS_Clip(
+    RS_Example(),
+    1,
+    ST_GeomFromText('POLYGON ((60 90, 160 90, 160 190, 60 190, 60 90))', 'OGC:CRS84')
+  )
+);
+```
+
+```sql
+SELECT RS_Width(
+  RS_Clip(
+    RS_Example(),
+    1,
+    ST_GeomFromText('POLYGON ((60 90, 160 90, 160 190, 60 190, 60 90))', 'OGC:CRS84'),
+    false,
+    0.0,
+    false
+  )
+);
+```

rust/sedona-raster-gdal/Cargo.toml

@@ -42,6 +42,7 @@ sedona-common = { workspace = true }
 sedona-expr = { workspace = true }
 sedona-functions = { workspace = true }
 sedona-gdal = { workspace = true }
+sedona-proj = { workspace = true }
 sedona-raster = { workspace = true }
 sedona-raster-functions = { workspace = true }
 sedona-schema = { workspace = true }
@@ -54,10 +55,16 @@ bindgen = ["sedona-gdal/bindgen"]
 [dev-dependencies]
 criterion = { workspace = true }
 sedona-gdal = { workspace = true, features = ["gdal-sys"] }
-sedona-testing = { workspace = true }
+sedona-proj = { workspace = true, features = ["proj-sys"] }
+sedona-testing = { workspace = true, features = ["criterion"] }
 tempfile = { workspace = true }
 tokio = { workspace = true, features = ["rt-multi-thread", "macros"] }
 
+[[bench]]
+harness = false
+name = "rs_clip"
+path = "benches/rs_clip.rs"
+
 [[bench]]
 harness = false
 name = "rs_frompath"

rust/sedona-raster-gdal/benches/rs_clip.rs

@@ -0,0 +1,140 @@
+// Licensed to the Apache Software Foundation (ASF) under one
+// or more contributor license agreements.  See the NOTICE file
+// distributed with this work for additional information
+// regarding copyright ownership.  The ASF licenses this file
+// to you under the Apache License, Version 2.0 (the
+// "License"); you may not use this file except in compliance
+// with the License.  You may obtain a copy of the License at
+//
+//   http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing,
+// software distributed under the License is distributed on an
+// "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+// KIND, either express or implied.  See the License for the
+// specific language governing permissions and limitations
+// under the License.
+
+//! Benchmarks for the RS_Clip UDF.
+//!
+//! RS_Clip rasterizes the clip geometry into a mask, sets pixels outside it to
+//! nodata, and (by default) crops to the geometry's bounding box.
+//!
+//! Each case builds a raster whose world extent is exactly the clip-polygon
+//! generator's `[-10, 10]²` bounds at the requested resolution, so every
+//! generated polygon lands on the raster and the full mask/crop path runs.
+//! `all_touched = true` guarantees an overlapping polygon burns at least one
+//! pixel even when it is smaller than a cell (otherwise a sub-pixel polygon
+//! would produce an empty mask and hit the no-intersection early return — which
+//! is what an earlier version of this benchmark accidentally measured, because
+//! it placed the rasters far outside the polygons entirely).
+//!
+//! Axes:
+//! - **Raster resolution** (`64²`, `256²`, `1024²`) with a small polygon: cost
+//!   is dominated by O(width × height) mask handling (rasterize + mask scan).
+//! - **Clip polygon complexity** (vertex count) at a fixed resolution, driving
+//!   the GDAL rasterization cost.
+//! - **Large clip**: a polygon covering most of the raster, so the mask/crop
+//!   copy (`apply_mask_and_crop`) dominates rather than mask rasterization.
+
+use std::sync::Arc;
+
+use arrow_array::{ArrayRef, BinaryArray, BooleanArray, Int32Array};
+use arrow_schema::DataType;
+use criterion::{criterion_group, criterion_main, Criterion};
+use datafusion_expr::ScalarUDF;
+use sedona_schema::datatypes::{SedonaType, RASTER, WKB_GEOMETRY};
+use sedona_testing::{
+    benchmark_util::BenchmarkArgSpec, create::make_wkb, raster_spec::RasterSpec,
+    testers::ScalarUdfTester,
+};
+
+fn criterion_benchmark(c: &mut Criterion) {
+    let f = sedona_raster_gdal::register::default_function_set();
+    let udf: ScalarUDF = f
+        .scalar_udf("rs_clip")
+        .expect("rs_clip is registered")
+        .clone()
+        .into();
+
+    // RS_Clip(raster, band, geom, all_touched).
+    let tester = ScalarUdfTester::new(
+        udf,
+        vec![
+            RASTER,
+            SedonaType::Arrow(DataType::Int32),
+            WKB_GEOMETRY,
+            SedonaType::Arrow(DataType::Boolean),
+        ],
+    );
+
+    let band: ArrayRef = Arc::new(Int32Array::from(vec![1]));
+    let all_touched: ArrayRef = Arc::new(BooleanArray::from(vec![true]));
+
+    // A north-up raster covering exactly the polygon generator's [-10, 10]²
+    // bounds at the requested resolution, so every generated polygon overlaps.
+    let build_raster = |w: i64, h: i64| -> ArrayRef {
+        let transform = [-10.0, 20.0 / w as f64, 0.0, 10.0, 0.0, -20.0 / h as f64];
+        let values = vec![1u8; (w * h) as usize];
+        Arc::new(
+            RasterSpec::d2(w, h)
+                .band_values(&values)
+                .crs(None)
+                .transform(transform)
+                .build(),
+        )
+    };
+
+    // `c` is passed in rather than captured, so several cases can share it.
+    let run = |c: &mut Criterion, label: &str, raster: ArrayRef, geom: ArrayRef| {
+        c.bench_function(label, |b| {
+            b.iter(|| {
+                tester
+                    .invoke_arrays(vec![
+                        raster.clone(),
+                        band.clone(),
+                        geom.clone(),
+                        all_touched.clone(),
+                    ])
+                    .unwrap()
+            })
+        });
+    };
+
+    let gen_polygon = |vertices: usize| -> ArrayRef {
+        BenchmarkArgSpec::Polygon(vertices)
+            .build_arrays(0, 1, 1)
+            .expect("build clip polygon")
+            .remove(0)
+    };
+
+    // Resolution sweep (simple 8-vertex polygon).
+    for (w, h) in [(64i64, 64i64), (256, 256), (1024, 1024)] {
+        let label = format!("raster-gdal rs_clip Clip(Raster({w}x{h}), Polygon(8))");
+        run(c, &label, build_raster(w, h), gen_polygon(8));
+    }
+
+    // Polygon-complexity axis at a fixed 64×64 resolution.
+    run(
+        c,
+        "raster-gdal rs_clip Clip(Raster(64x64), Polygon(50))",
+        build_raster(64, 64),
+        gen_polygon(50),
+    );
+
+    // Large clip: the polygon covers nearly the whole raster, so the crop copy
+    // (apply_mask_and_crop) dominates — this is where its row memcpy matters.
+    let big_geom: ArrayRef = Arc::new(BinaryArray::from_iter_values([make_wkb(
+        "POLYGON ((-9.5 -9.5, 9.5 -9.5, 9.5 9.5, -9.5 9.5, -9.5 -9.5))",
+    )
+    .as_slice()]));
+    run(
+        c,
+        "raster-gdal rs_clip Clip(Raster(1024x1024), Polygon(large))",
+        build_raster(1024, 1024),
+        big_geom,
+    );
+}
+
+criterion_group!(benches, criterion_benchmark);
+criterion_main!(benches);

rust/sedona-raster-gdal/src/lib.rs

@@ -33,6 +33,7 @@ mod gdal_common;
 mod gdal_dataset_provider;
 
 mod raster_loader;
+mod rs_clip;
 mod rs_frompath;
 mod rs_metadata;
 mod source_uri;
@@ -44,6 +45,7 @@ pub use gdal_common::{
     nodata_bytes_to_f64, nodata_f64_to_bytes, GdalBandLayout, GdalBandPlan,
 };
 pub use raster_loader::{GdalLoader, GDAL_FORMAT};
+pub use rs_clip::rs_clip_udf;
 pub use rs_frompath::rs_frompath_udf;
 pub use rs_metadata::rs_metadata_udf;
 pub use utils::{

rust/sedona-raster-gdal/src/register.rs

@@ -20,6 +20,7 @@ use sedona_expr::function_set::FunctionSet;
 /// Export the set of GDAL-backed functions defined in this crate.
 pub fn default_function_set() -> FunctionSet {
     let mut function_set = FunctionSet::new();
+    function_set.insert_scalar_udf(crate::rs_clip::rs_clip_udf());
     function_set.insert_scalar_udf(crate::rs_frompath::rs_frompath_udf());
     function_set.insert_scalar_udf(crate::rs_metadata::rs_metadata_udf());
     function_set