feat(rust/sedona-raster-functions): add RS_Value for spatial 2d rasters

docs/reference/sql/rs_value.qmd

@@ -0,0 +1,56 @@
+---
+# 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_Value
+description: >
+  Returns the value of a single raster pixel as a double, selected by a point
+  geometry. Returns null if the location is outside the raster or the pixel
+  holds the band's nodata value.
+kernels:
+  - returns: float64
+    args:
+    - raster
+    - name: point
+      type: geometry
+      description: >
+        The pixel that contains this point is sampled (no resampling).
+        Reprojected into the raster CRS when both carry one.
+  - returns: float64
+    args:
+    - raster
+    - name: point
+      type: geometry
+    - name: band
+      type: int
+      description: Band index (1-based). Defaults to 1 if not specified.
+---
+
+## Description
+
+`RS_Value` samples one pixel of a raster. The location is given as a point
+geometry — the value of the pixel that contains the point is returned, with no
+interpolation. The band defaults to 1.
+
+The result is `NULL` when the point falls outside the raster, or when the
+sampled pixel equals the band's nodata value. Only 2-D rasters are supported.
+
+## Examples
+
+```sql
+SELECT RS_Value(RS_Example(), ST_Point(74.58, 110.57, 'OGC:CRS84'));
+```

python/sedonadb/pyproject.toml

@@ -42,6 +42,7 @@ test = [
     "polars",
     "pytest",
     "pyyaml",
+    "rasterio",
 ]
 geopandas = [
     "adbc-driver-manager[dbapi]",

python/sedonadb/tests/functions/test_raster_functions.py

@@ -102,3 +102,123 @@ def test_rs_ensureloaded(con, sedona_testing):
     assert arr.shape == (512, 512)
     assert arr.dtype == "uint16"
     assert arr[0, 0] == 2324
+
+
+# Point sampling. RS_Example fills band `b` with the constant value `b`, except
+# the top-left pixel which is set to the nodata value (127). (74.58, 110.57) is
+# the centroid of pixel (10, 10) (0-based) in the raster's OGC:CRS84 space; the
+# point and raster share a CRS so no reprojection happens. A point far outside
+# the footprint yields NULL. (The `needs_pixels` -> RS_EnsureLoaded planner path
+# is covered against a real OutDb raster by `test_rs_ensureloaded`.)
+@pytest.mark.parametrize(
+    ("expr", "expected"),
+    [
+        (
+            "RS_Value(RS_Example(), ST_SetCRS(ST_Point(74.58, 110.57), 'OGC:CRS84'))",
+            1.0,
+        ),
+        (
+            "RS_Value(RS_Example(), ST_SetCRS(ST_Point(74.58, 110.57), 'OGC:CRS84'), 2)",
+            2.0,
+        ),
+        (
+            "RS_Value(RS_Example(), ST_SetCRS(ST_Point(74.58, 110.57), 'OGC:CRS84'), 3)",
+            3.0,
+        ),
+        ("RS_Value(RS_Example(), ST_SetCRS(ST_Point(0.0, 0.0), 'OGC:CRS84'))", None),
+        # POINT EMPTY has no location to sample -> NULL (not an error).
+        (
+            "RS_Value(RS_Example(), ST_SetCRS(ST_GeomFromText('POINT EMPTY'), 'OGC:CRS84'))",
+            None,
+        ),
+    ],
+)
+def test_rs_value_point(expr, expected):
+    SedonaDB().assert_query_result(f"SELECT {expr}", expected)
+
+
+def test_rs_value_matches_rasterio(con):
+    """Cross-check RS_Value against rasterio on a random raster.
+
+    Builds an in-memory raster from a random numpy array with a known
+    geotransform and no CRS (so neither engine reprojects), then samples a dense
+    set of points and asserts RS_Value returns exactly what rasterio reads at the
+    same world coordinates. Points cover every pixel center plus four off-center
+    positions per pixel (toward the corners, kept inside the pixel to avoid floor
+    ambiguity at exact boundaries) and a batch of random interior points.
+    """
+    import numpy as np
+    import pandas as pd
+
+    pytest.importorskip("rasterio")
+    from rasterio.io import MemoryFile
+    from rasterio.transform import Affine
+
+    from sedonadb.raster import Raster
+
+    rng = np.random.default_rng(42)
+    height, width = 7, 5
+    data = rng.random((height, width)) * 1000.0
+
+    # GDAL-order geotransform: origin (100, 500), 2-wide pixels, -3 tall
+    # (north-up), no skew. Shared verbatim by both engines.
+    gdal_transform = (100.0, 2.0, 0.0, 500.0, 0.0, -3.0)
+    affine = Affine.from_gdal(*gdal_transform)
+
+    # Sample points in pixel space (col_frac, row_frac).
+    pixel_points = []
+    for row in range(height):
+        for col in range(width):
+            for du, dv in [
+                (0.5, 0.5),
+                (0.25, 0.25),
+                (0.75, 0.75),
+                (0.25, 0.75),
+                (0.75, 0.25),
+            ]:
+                pixel_points.append((col + du, row + dv))
+    n_random = 150
+    rand_cols = rng.integers(0, width, n_random)
+    rand_rows = rng.integers(0, height, n_random)
+    pixel_points.extend(
+        zip(
+            rand_cols + rng.uniform(0.1, 0.9, n_random),
+            rand_rows + rng.uniform(0.1, 0.9, n_random),
+        )
+    )
+
+    # Map pixel-space positions to world coordinates via the shared affine.
+    xs, ys = zip(*(affine * (u, v) for u, v in pixel_points))
+
+    # rasterio reference: a real GDAL read of the same array (no CRS).
+    with MemoryFile() as mem:
+        with mem.open(
+            driver="GTiff",
+            height=height,
+            width=width,
+            count=1,
+            dtype="float64",
+            transform=affine,
+        ) as dst:
+            dst.write(data, 1)
+        with mem.open() as src:
+            expected = [vals[0] for vals in src.sample(list(zip(xs, ys)))]
+
+    # sedonadb: sample the same points via RS_Value over a scalar raster.
+    raster = Raster.from_numpy(data, transform=gdal_transform)
+    pts = con.create_data_frame(pd.DataFrame({"idx": range(len(xs)), "x": xs, "y": ys}))
+    view = "test_rs_value_matches_rasterio_pts"
+    pts.to_view(view)
+    try:
+        got = (
+            con.sql(
+                f"SELECT RS_Value($1, ST_Point(x, y)) AS v FROM {view} ORDER BY idx",
+                params=(raster,),
+            )
+            .to_arrow_table()["v"]
+            .to_pylist()
+        )
+    finally:
+        con.drop_view(view)
+
+    assert got == pytest.approx(expected)