You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@sedona.apache.org by pa...@apache.org on 2026/06/24 21:43:32 UTC
(sedona-db) branch main updated: feat(rust/sedona-raster-functions): add RS_Value for spatial 2d rasters (#974)
This is an automated email from the ASF dual-hosted git repository.
paleolimbot pushed a commit to branch main
in repository https://gitbox.apache.org/repos/asf/sedona-db.git
The following commit(s) were added to refs/heads/main by this push:
new daf5fd48f feat(rust/sedona-raster-functions): add RS_Value for spatial 2d rasters (#974)
daf5fd48f is described below
commit daf5fd48f5b7bf7b1c480f9b834b32056f946de0
Author: James Willis <ja...@users.noreply.github.com>
AuthorDate: Wed Jun 24 14:43:26 2026 -0700
feat(rust/sedona-raster-functions): add RS_Value for spatial 2d rasters (#974)
---
docs/reference/sql/rs_value.qmd | 56 ++
python/sedonadb/pyproject.toml | 1 +
.../tests/functions/test_raster_functions.py | 120 +++
rust/sedona-geometry/src/wkb_header.rs | 177 +++-
.../benches/native-raster-functions.rs | 37 +
rust/sedona-raster-functions/src/executor.rs | 11 +-
rust/sedona-raster-functions/src/lib.rs | 1 +
rust/sedona-raster-functions/src/register.rs | 1 +
rust/sedona-raster-functions/src/rs_value.rs | 992 +++++++++++++++++++++
9 files changed, 1381 insertions(+), 15 deletions(-)
diff --git a/docs/reference/sql/rs_value.qmd b/docs/reference/sql/rs_value.qmd
new file mode 100644
index 000000000..f99d7efc0
--- /dev/null
+++ b/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'));
+```
diff --git a/python/sedonadb/pyproject.toml b/python/sedonadb/pyproject.toml
index a610d1641..b0f83c6a8 100644
--- a/python/sedonadb/pyproject.toml
+++ b/python/sedonadb/pyproject.toml
@@ -42,6 +42,7 @@ test = [
"polars",
"pytest",
"pyyaml",
+ "rasterio",
]
geopandas = [
"adbc-driver-manager[dbapi]",
diff --git a/python/sedonadb/tests/functions/test_raster_functions.py b/python/sedonadb/tests/functions/test_raster_functions.py
index fd927f0fd..37640bcf8 100644
--- a/python/sedonadb/tests/functions/test_raster_functions.py
+++ b/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)
diff --git a/rust/sedona-geometry/src/wkb_header.rs b/rust/sedona-geometry/src/wkb_header.rs
index b1cdfcac0..2e1b1c8e0 100644
--- a/rust/sedona-geometry/src/wkb_header.rs
+++ b/rust/sedona-geometry/src/wkb_header.rs
@@ -152,6 +152,77 @@ impl WkbHeader {
}
}
+/// Read the `(x, y)` of a WKB Point without parsing the geometry generally.
+///
+/// This is a fast path for the common case of sampling/indexing by a single
+/// point: a Point WKB has a fixed layout (`byte order`, `type code`, optional
+/// SRID, then the coordinate), so the coordinate is at a known offset and can
+/// be read with a couple of fixed-offset loads. It deliberately handles only
+/// Points — any other geometry is an error.
+///
+/// Returns `Ok(None)` for `POINT EMPTY` (encoded as `NaN`/`NaN`), `Ok(Some((x, y)))`
+/// for a finite point, and an error if the bytes are not a Point or are truncated.
+/// Accepts both ISO (with Z/M dimension flags) and EWKB (with Z/M/SRID flags)
+/// encodings; only the first two ordinates are read.
+pub fn read_point_xy(buf: &[u8]) -> Result<Option<(f64, f64)>, SedonaGeometryError> {
+ let little_endian = match buf.first() {
+ Some(1) => true,
+ Some(0) => false,
+ _ => {
+ return Err(SedonaGeometryError::Invalid(
+ "WKB: missing or invalid byte order".to_string(),
+ ))
+ }
+ };
+
+ let read_u32 = |offset: usize| -> Result<u32, SedonaGeometryError> {
+ let bytes: [u8; 4] = buf
+ .get(offset..offset + 4)
+ .ok_or_else(|| SedonaGeometryError::Invalid("WKB: truncated header".to_string()))?
+ .try_into()
+ .unwrap();
+ Ok(if little_endian {
+ u32::from_le_bytes(bytes)
+ } else {
+ u32::from_be_bytes(bytes)
+ })
+ };
+
+ let read_f64 = |offset: usize| -> Result<f64, SedonaGeometryError> {
+ let bytes: [u8; 8] = buf
+ .get(offset..offset + 8)
+ .ok_or_else(|| SedonaGeometryError::Invalid("WKB: truncated coordinate".to_string()))?
+ .try_into()
+ .unwrap();
+ Ok(if little_endian {
+ f64::from_le_bytes(bytes)
+ } else {
+ f64::from_be_bytes(bytes)
+ })
+ };
+
+ let type_code = read_u32(1)?;
+ // The low three bits of the base type identify a Point (1) regardless of the
+ // ISO dimension digits or EWKB high-bit flags.
+ if type_code & 0x7 != 1 {
+ return Err(SedonaGeometryError::Invalid(format!(
+ "WKB: expected a Point (type code {type_code})"
+ )));
+ }
+
+ // EWKB carries a 4-byte SRID between the type code and the coordinate.
+ let has_srid = type_code & SRID_FLAG_BIT != 0;
+ let coord_offset = if has_srid { 9 } else { 5 };
+
+ let x = read_f64(coord_offset)?;
+ let y = read_f64(coord_offset + 8)?;
+ if x.is_nan() && y.is_nan() {
+ Ok(None)
+ } else {
+ Ok(Some((x, y)))
+ }
+}
+
// A helper struct for calculating the WKBHeader
struct WkbBuffer<'a> {
buf: &'a [u8],
@@ -377,6 +448,7 @@ fn calc_dimensions(code: u32) -> Result<Dimensions, SedonaGeometryError> {
#[cfg(test)]
mod tests {
use super::*;
+ use sedona_testing::fixtures::*;
use std::str::FromStr;
use wkb::writer::{write_geometry, WriteOptions};
use wkt::Wkt;
@@ -514,8 +586,6 @@ mod tests {
#[test]
fn ewkb() {
- use sedona_testing::fixtures::*;
-
// Test POINT with SRID 4326
let header = WkbHeader::try_new(&POINT_WITH_SRID_4326_EWKB).unwrap();
assert_eq!(header.srid(), 4326);
@@ -598,8 +668,6 @@ mod tests {
#[test]
fn srid_linestring() {
- use sedona_testing::fixtures::*;
-
let header = WkbHeader::try_new(&LINESTRING_WITH_SRID_4326_EWKB).unwrap();
assert_eq!(header.srid(), 4326);
assert_eq!(
@@ -613,8 +681,6 @@ mod tests {
#[test]
fn srid_polygon() {
- use sedona_testing::fixtures::*;
-
let header = WkbHeader::try_new(&POLYGON_WITH_SRID_4326_EWKB).unwrap();
assert_eq!(header.srid(), 4326);
assert_eq!(header.geometry_type_id().unwrap(), GeometryTypeId::Polygon);
@@ -625,8 +691,6 @@ mod tests {
#[test]
fn multipoint_with_srid() {
- use sedona_testing::fixtures::*;
-
let header = WkbHeader::try_new(&MULTIPOINT_WITH_SRID_4326_EWKB).unwrap();
assert_eq!(header.srid(), 4326);
assert_eq!(
@@ -640,8 +704,6 @@ mod tests {
#[test]
fn srid_empty_geometries_with_srid() {
- use sedona_testing::fixtures::*;
-
// Test POINT EMPTY with SRID
let header = WkbHeader::try_new(&POINT_EMPTY_WITH_SRID_4326_EWKB).unwrap();
assert_eq!(header.srid(), 4326);
@@ -944,7 +1006,6 @@ mod tests {
#[test]
fn incomplete_ewkb_buffers() {
- use sedona_testing::fixtures::*;
// Test incomplete EWKB buffers
// 1 byte_order + 4 geometry type + 4 srid + 8 x + 8 y
@@ -1059,4 +1120,98 @@ mod tests {
let header = WkbHeader::try_new(&wkb).unwrap();
assert!(header.is_empty().unwrap());
}
+
+ #[test]
+ fn read_point_xy_iso_little_endian() {
+ // 2-D / 3-D / 4-D ISO points all expose the same leading (x, y), and the
+ // result must agree with the general parser's `first_xy`.
+ for wkt_value in [
+ "POINT (1 2)",
+ "POINT Z (1 2 3)",
+ "POINT M (1 2 3)",
+ "POINT ZM (1 2 3 4)",
+ ] {
+ let wkb = make_wkb(wkt_value);
+ assert_eq!(
+ read_point_xy(&wkb).unwrap(),
+ Some((1.0, 2.0)),
+ "{wkt_value}"
+ );
+ assert_eq!(
+ read_point_xy(&wkb).unwrap(),
+ Some(WkbHeader::try_new(&wkb).unwrap().first_xy()),
+ "disagrees with first_xy for {wkt_value}"
+ );
+ }
+ }
+
+ #[test]
+ fn read_point_xy_big_endian() {
+ // Byte order 0x00, type 1 (BE), x = 1.0, y = 2.0 — all big-endian.
+ let wkb: [u8; 21] = [
+ 0x00, // big-endian
+ 0x00, 0x00, 0x00, 0x01, // type code 1 (Point)
+ 0x3F, 0xF0, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, // x = 1.0
+ 0x40, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, // y = 2.0
+ ];
+ assert_eq!(read_point_xy(&wkb).unwrap(), Some((1.0, 2.0)));
+ }
+
+ #[test]
+ fn read_point_xy_ewkb_with_srid() {
+ // The SRID sits between the type code and the coordinate; the reader must
+ // skip it for every dimensionality.
+ for wkb in [
+ &POINT_WITH_SRID_4326_EWKB[..],
+ &POINT_Z_WITH_SRID_3857_EWKB[..],
+ &POINT_M_WITH_SRID_4326_EWKB[..],
+ &POINT_ZM_WITH_SRID_4326_EWKB[..],
+ ] {
+ assert_eq!(read_point_xy(wkb).unwrap(), Some((1.0, 2.0)));
+ assert_eq!(
+ read_point_xy(wkb).unwrap(),
+ Some(WkbHeader::try_new(wkb).unwrap().first_xy())
+ );
+ }
+ }
+
+ #[test]
+ fn read_point_xy_empty_is_none() {
+ // POINT EMPTY is encoded as NaN/NaN and reads back as "no coordinate".
+ assert_eq!(read_point_xy(&make_wkb("POINT EMPTY")).unwrap(), None);
+ assert_eq!(read_point_xy(&make_wkb("POINT Z EMPTY")).unwrap(), None);
+ assert_eq!(
+ read_point_xy(&POINT_EMPTY_WITH_SRID_4326_EWKB).unwrap(),
+ None
+ );
+ }
+
+ #[test]
+ fn read_point_xy_rejects_non_point() {
+ for wkt_value in [
+ "LINESTRING (1 2, 3 4)",
+ "POLYGON ((0 0, 0 1, 1 0, 0 0))",
+ "MULTIPOINT ((1 2))",
+ "GEOMETRYCOLLECTION (POINT (1 2))",
+ ] {
+ let err = read_point_xy(&make_wkb(wkt_value)).unwrap_err().to_string();
+ assert!(err.contains("expected a Point"), "{wkt_value}: {err}");
+ }
+ }
+
+ #[test]
+ fn read_point_xy_rejects_truncated() {
+ let wkb = make_wkb("POINT (1 2)");
+ // Empty input has no byte order; everything shorter than a full 2-D point
+ // header + coordinate must error rather than read past the end.
+ assert!(read_point_xy(&[]).is_err());
+ for i in 1..wkb.len() {
+ assert!(
+ read_point_xy(&wkb[0..i]).is_err(),
+ "0..{i} unexpectedly succeeded"
+ );
+ }
+ // Invalid byte-order flag.
+ assert!(read_point_xy(&[0x02, 0, 0, 0, 1]).is_err());
+ }
}
diff --git a/rust/sedona-raster-functions/benches/native-raster-functions.rs b/rust/sedona-raster-functions/benches/native-raster-functions.rs
index 2b6ffeef1..156ae4c0f 100644
--- a/rust/sedona-raster-functions/benches/native-raster-functions.rs
+++ b/rust/sedona-raster-functions/benches/native-raster-functions.rs
@@ -214,6 +214,43 @@ fn criterion_benchmark(c: &mut Criterion) {
),
);
+ // RS_Value(raster, point) — array raster (a distinct raster per row).
+ benchmark::scalar(
+ c,
+ &f,
+ "native-raster",
+ "rs_value",
+ BenchmarkArgs::ArrayArray(
+ Raster(64, 64),
+ Transformed(Box::new(Point), sd_apply_default_crs_udf().into()),
+ ),
+ );
+ // RS_Value(raster, point) — scalar raster (one raster, many points). This
+ // is the case where per-point band/nodata/buffer resolution can be hoisted
+ // out of the sample loop, so it is the primary target for optimization.
+ benchmark::scalar(
+ c,
+ &f,
+ "native-raster",
+ "rs_value",
+ BenchmarkArgs::ScalarArray(
+ Raster(64, 64),
+ Transformed(Box::new(Point), sd_apply_default_crs_udf().into()),
+ ),
+ );
+ // RS_Value(raster, point, band)
+ benchmark::scalar(
+ c,
+ &f,
+ "native-raster",
+ "rs_value",
+ BenchmarkArgs::ArrayArrayScalar(
+ Raster(64, 64),
+ Transformed(Box::new(Point), sd_apply_default_crs_udf().into()),
+ Int32(1, 2),
+ ),
+ );
+
// RS_Intersects(raster, geometry) - point
benchmark::scalar(
c,
diff --git a/rust/sedona-raster-functions/src/executor.rs b/rust/sedona-raster-functions/src/executor.rs
index 34959bf7d..39cbb2573 100644
--- a/rust/sedona-raster-functions/src/executor.rs
+++ b/rust/sedona-raster-functions/src/executor.rs
@@ -47,7 +47,7 @@ pub struct RasterExecutor<'a, 'b> {
// operations are expensive relative to per-element dispatch overhead, the cost
// of matching on each access is negligible in practice.
#[derive(Clone)]
-enum ItemWkbAccessor {
+pub(crate) enum ItemWkbAccessor {
Binary(BinaryArray),
BinaryView(BinaryViewArray),
}
@@ -76,7 +76,7 @@ impl ItemWkbAccessor {
// Same enum-dispatch rationale as `ItemWkbAccessor` above: the per-element
// match cost is dwarfed by the raster and CRS operations performed on each row.
-enum GeomWkbCrsAccessor {
+pub(crate) enum GeomWkbCrsAccessor {
WkbArray {
wkb: ItemWkbAccessor,
static_crs: Crs,
@@ -104,7 +104,7 @@ enum GeomWkbCrsAccessor {
impl GeomWkbCrsAccessor {
#[inline]
- fn get(&mut self, i: usize) -> Result<(Option<&[u8]>, CrsRef<'_>)> {
+ pub(crate) fn get(&mut self, i: usize) -> Result<(Option<&[u8]>, CrsRef<'_>)> {
match self {
Self::Null => Ok((None, None)),
Self::WkbArray { wkb, static_crs } => {
@@ -527,7 +527,10 @@ impl<'a, 'b> RasterExecutor<'a, 'b> {
}
}
- fn make_geom_wkb_crs_accessor(&self, arg_index: usize) -> Result<GeomWkbCrsAccessor> {
+ pub(crate) fn make_geom_wkb_crs_accessor(
+ &self,
+ arg_index: usize,
+ ) -> Result<GeomWkbCrsAccessor> {
let sedona_type = self
.arg_types
.get(arg_index)
diff --git a/rust/sedona-raster-functions/src/lib.rs b/rust/sedona-raster-functions/src/lib.rs
index 3a5baa8df..4de2d5975 100644
--- a/rust/sedona-raster-functions/src/lib.rs
+++ b/rust/sedona-raster-functions/src/lib.rs
@@ -38,4 +38,5 @@ pub mod rs_size;
pub mod rs_slice;
pub mod rs_spatial_predicates;
pub mod rs_srid;
+pub mod rs_value;
pub mod rs_worldcoordinate;
diff --git a/rust/sedona-raster-functions/src/register.rs b/rust/sedona-raster-functions/src/register.rs
index cfb2f008f..e3d118da4 100644
--- a/rust/sedona-raster-functions/src/register.rs
+++ b/rust/sedona-raster-functions/src/register.rs
@@ -77,6 +77,7 @@ pub fn default_function_set() -> FunctionSet {
crate::rs_spatial_predicates::rs_within_udf,
crate::rs_srid::rs_crs_udf,
crate::rs_srid::rs_srid_udf,
+ crate::rs_value::rs_value_udf,
crate::rs_worldcoordinate::rs_rastertoworldcoord_udf,
crate::rs_worldcoordinate::rs_rastertoworldcoordx_udf,
crate::rs_worldcoordinate::rs_rastertoworldcoordy_udf,
diff --git a/rust/sedona-raster-functions/src/rs_value.rs b/rust/sedona-raster-functions/src/rs_value.rs
new file mode 100644
index 000000000..b28b2a59a
--- /dev/null
+++ b/rust/sedona-raster-functions/src/rs_value.rs
@@ -0,0 +1,992 @@
+// 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.
+
+//! `RS_Value` — sample a raster's pixel value at a point.
+//!
+//! ```text
+//! RS_Value(raster, point) -> Double -- band defaults to 1
+//! RS_Value(raster, point, band) -> Double
+//! ```
+//!
+//! Returns the value of the pixel that contains the point (no resampling). The
+//! result is `NULL` when the raster/arguments are null, the point is empty, the
+//! point is out of bounds, or the value equals the band's nodata.
+//!
+//! The function is tagged [`NEEDS_PIXELS_METADATA_KEY`], so the planner wraps
+//! its raster argument in `RS_EnsureLoaded`; by the time a kernel runs the band
+//! bytes are materialised InDb and a value is read directly from the band's
+//! [`NdBuffer`](sedona_raster::traits::NdBuffer) — no GDAL involved. Only 2-D
+//! rasters are supported; a band with extra (non-spatial) dimensions errors.
+
+use std::sync::Arc;
+
+use arrow_array::{builder::Float64Builder, Array, ArrayRef, Float64Array, StructArray};
+use arrow_schema::DataType;
+use datafusion_common::cast::as_int32_array;
+use datafusion_common::{exec_datafusion_err, exec_err, Result, ScalarValue};
+use datafusion_expr::{ColumnarValue, Volatility};
+use sedona_expr::scalar_udf::{SedonaScalarKernel, SedonaScalarUDF};
+use sedona_geometry::wkb_header::read_point_xy;
+use sedona_proj::transform::with_global_proj_engine;
+use sedona_raster::affine_transformation::AffineMatrix;
+use sedona_raster::array::RasterStructArray;
+use sedona_raster::traits::{nodata_bytes_to_f64_lossless, NdBuffer, RasterRef};
+use sedona_schema::crs::CrsRef;
+use sedona_schema::{datatypes::SedonaType, matchers::ArgMatcher};
+
+use crate::crs_utils::{crs_transform_wkb, resolve_crs};
+use crate::executor::RasterExecutor;
+use crate::rs_ensure_loaded::NEEDS_PIXELS_METADATA_KEY;
+
+/// `RS_Value()` scalar UDF — sample a pixel value at a point.
+pub fn rs_value_udf() -> SedonaScalarUDF {
+ SedonaScalarUDF::new(
+ "rs_value",
+ vec![
+ Arc::new(RsValuePoint { with_band: false }), // RS_Value(raster, point)
+ Arc::new(RsValuePoint { with_band: true }), // RS_Value(raster, point, band)
+ ],
+ Volatility::Immutable,
+ )
+ // The kernels read pixel bytes, so the raster argument must be materialised
+ // InDb first; the planner injects RS_EnsureLoaded based on this flag.
+ .with_metadata(NEEDS_PIXELS_METADATA_KEY, "true")
+}
+
+/// Kernel for `RS_Value(raster, point[, band])`.
+#[derive(Debug)]
+struct RsValuePoint {
+ with_band: bool,
+}
+
+impl SedonaScalarKernel for RsValuePoint {
+ fn return_type(&self, args: &[SedonaType]) -> Result<Option<SedonaType>> {
+ let mut matchers = vec![
+ ArgMatcher::is_raster(),
+ ArgMatcher::is_geometry_or_geography(),
+ ];
+ if self.with_band {
+ matchers.push(ArgMatcher::is_integer());
+ }
+ let matcher = ArgMatcher::new(matchers, SedonaType::Arrow(DataType::Float64));
+ matcher.match_args(args)
+ }
+
+ fn invoke_batch(
+ &self,
+ arg_types: &[SedonaType],
+ args: &[ColumnarValue],
+ ) -> Result<ColumnarValue> {
+ // Fast path: a constant (scalar) raster lets us resolve the affine
+ // transform and CRS once for the whole batch instead of per point — the
+ // common RS_Value(raster_expr, point_column[, band]) shape. The band
+ // argument does not change this: a constant band also hoists its buffer,
+ // and only a band *column* falls back to per-row band resolution.
+ if let ColumnarValue::Scalar(ScalarValue::Struct(raster_struct)) = &args[0] {
+ return self.invoke_scalar_raster(arg_types, args, raster_struct.as_ref());
+ }
+
+ let executor = RasterExecutor::new(arg_types, args);
+ let num_iterations = executor.num_iterations();
+ let mut builder = Float64Builder::with_capacity(num_iterations);
+
+ // The optional band argument, materialised once as an Int32 array. Held
+ // as an `ArrayRef` so the typed view below borrows it instead of cloning
+ // the typed `Int32Array`.
+ let band_arr = if self.with_band {
+ Some(int32_array_arg(&args[2], num_iterations)?)
+ } else {
+ None
+ };
+ let band_array = band_arr.as_ref().map(|a| as_int32_array(a)).transpose()?;
+ let mut band_iter = band_array.map(|a| a.iter());
+
+ // Reprojecting the point into the raster CRS needs a PROJ engine.
+ with_global_proj_engine(|engine| {
+ executor.execute_raster_wkb_crs_void(|raster_opt, wkb_opt, point_crs| {
+ let (raster, point_wkb, band_num) =
+ match (raster_opt, wkb_opt, next_band(&mut band_iter)) {
+ (Some(raster), Some(point_wkb), Some(band_num)) => {
+ (raster, point_wkb, band_num)
+ }
+ _ => {
+ builder.append_null();
+ return Ok(());
+ }
+ };
+
+ // Parse the point and bring it into the raster's CRS. Null/empty
+ // points (and non-finite coordinates) have no location to sample.
+ let raster_crs = resolve_crs(raster.crs())?;
+ let Some((x, y)) =
+ resolve_point_xy(point_wkb, point_crs, raster_crs.as_deref(), false, engine)?
+ else {
+ builder.append_null();
+ return Ok(());
+ };
+
+ let affine = AffineMatrix::from_metadata(&raster.metadata());
+ match xy_to_pixel(&affine, x, y)? {
+ Some((col, row)) => match sample_pixel(raster, col, row, band_num)? {
+ Some(value) => builder.append_value(value),
+ None => builder.append_null(),
+ },
+ None => builder.append_null(),
+ }
+ Ok(())
+ })
+ })?;
+
+ executor.finish(Arc::new(builder.finish()))
+ }
+}
+
+impl RsValuePoint {
+ /// Optimized path for a constant (scalar) raster: the affine transform and
+ /// raster CRS are resolved once for the whole batch, then a selection vector
+ /// drives a tight sample loop. This serves every band shape:
+ /// - no band argument or a constant band → the band buffer is hoisted too;
+ /// - a band *column* → the band buffer is resolved per row (its `NdBuffer`
+ /// borrows from the band, which can't be cached across distinct bands),
+ /// but the affine/CRS/reproject work is still hoisted.
+ ///
+ /// Sampling behaviour matches the general path: it uses the same
+ /// [`resolve_point_xy`]/[`xy_to_pixel`] helpers, so null/empty/non-finite
+ /// points yield NULL identically. Band resolution (and its 2-D check) is
+ /// deferred until at least one point needs sampling, so an all-null/all-empty
+ /// batch returns NULL without touching the band — as the general path does.
+ fn invoke_scalar_raster(
+ &self,
+ arg_types: &[SedonaType],
+ args: &[ColumnarValue],
+ raster_struct: &StructArray,
+ ) -> Result<ColumnarValue> {
+ let executor = RasterExecutor::new(arg_types, args);
+ let n = executor.num_iterations();
+
+ let rasters = RasterStructArray::try_new(raster_struct)?;
+ if rasters.is_null(0) {
+ // A NULL raster makes every output NULL.
+ return executor.finish(Arc::new(Float64Array::from(vec![None; n])));
+ }
+ let raster = rasters.get(0)?;
+
+ // Band selection: a missing band argument (default 1) or a scalar band is
+ // constant for the batch and lets us hoist the band buffer; a band column
+ // is resolved per row. A NULL scalar band makes every output NULL.
+ let mut const_band: Option<usize> = None;
+ let mut band_values: Option<ArrayRef> = None;
+ if self.with_band {
+ match &args[2] {
+ ColumnarValue::Scalar(scalar) => {
+ let arr = ColumnarValue::Scalar(scalar.clone())
+ .cast_to(&DataType::Int32, None)?
+ .into_array(1)?;
+ let arr = as_int32_array(&arr)?;
+ if arr.is_null(0) {
+ return executor.finish(Arc::new(Float64Array::from(vec![None; n])));
+ }
+ // Match `next_band`: clamp to 0 so band 0/negative surface as a
+ // not-1-based error from `Bands::band` rather than being coerced.
+ const_band = Some(arr.value(0).max(0) as usize);
+ }
+ other => band_values = Some(int32_array_arg(other, n)?),
+ }
+ } else {
+ const_band = Some(1);
+ }
+ let band_array = band_values
+ .as_ref()
+ .map(|a| as_int32_array(a))
+ .transpose()?;
+
+ // Affine transform and raster CRS, resolved once for all points.
+ let affine = AffineMatrix::from_metadata(&raster.metadata());
+ let raster_crs = resolve_crs(raster.crs())?;
+ // Decide reprojection once when the point CRS is column-level (the common
+ // case). This skips a per-point `crs_equals`, which allocates a String —
+ // ~15 ns/point, uniform across CRS types.
+ let needs_reproject = column_reproject_decision(&arg_types[1], raster_crs.as_deref())?;
+ let skip_reproject = needs_reproject == Some(false);
+
+ let mut geom = executor.make_geom_wkb_crs_accessor(1)?;
+
+ // Phase 1 — selection vector: collect (row, x, y, band) for the points
+ // worth sampling (non-null, non-empty Point with a non-null band),
+ // reprojected into the raster CRS. Skipped rows stay NULL in the output.
+ let mut selection: Vec<(usize, f64, f64, usize)> = Vec::with_capacity(n);
+ let mut band_iter = band_array.map(|a| a.iter());
+ with_global_proj_engine(|engine| {
+ for i in 0..n {
+ // Advance the band column in lockstep with the row index; a NULL
+ // band element leaves the row NULL (matching the general path).
+ let band_num = match const_band {
+ Some(b) => b,
+ None => match next_band(&mut band_iter) {
+ Some(b) => b,
+ None => continue,
+ },
+ };
+ let (wkb_opt, point_crs) = geom.get(i)?;
+ let Some(point_wkb) = wkb_opt else {
+ continue;
+ };
+ if let Some((x, y)) = resolve_point_xy(
+ point_wkb,
+ point_crs,
+ raster_crs.as_deref(),
+ skip_reproject,
+ engine,
+ )? {
+ selection.push((i, x, y, band_num));
+ }
+ }
+ Ok(())
+ })?;
+
+ let mut out: Vec<Option<f64>> = vec![None; n];
+ if selection.is_empty() {
+ return executor.finish(Arc::new(Float64Array::from(out)));
+ }
+
+ // Phase 2 — sample. A constant band resolves its buffer/nodata once (now
+ // that we know a point needs sampling); a band column resolves per row.
+ match const_band {
+ Some(band_num) => {
+ let band = raster
+ .bands()
+ .band(band_num)
+ .map_err(|e| exec_datafusion_err!("RS_Value: {e}"))?;
+ if !band.is_spatial_2d() {
+ return exec_err!(
+ "RS_Value supports 2-D rasters only; band is not a 2-D (y, x) grid"
+ );
+ }
+ let buffer = band
+ .nd_buffer()
+ .map_err(|e| exec_datafusion_err!("RS_Value: {e}"))?;
+ let nodata = band
+ .nodata_as_f64()
+ .map_err(|e| exec_datafusion_err!("RS_Value: {e}"))?;
+ for (i, x, y, _band) in selection {
+ if let Some((col, row)) = xy_to_pixel(&affine, x, y)? {
+ out[i] = read_pixel(&buffer, nodata, col, row)?;
+ }
+ }
+ }
+ None => {
+ for (i, x, y, band_num) in selection {
+ if let Some((col, row)) = xy_to_pixel(&affine, x, y)? {
+ out[i] = sample_pixel(&raster, col, row, band_num)?;
+ }
+ }
+ }
+ }
+
+ executor.finish(Arc::new(Float64Array::from(out)))
+ }
+}
+
+/// Materialise an integer argument as an owned `Int32` [`ArrayRef`] for the
+/// batch. Callers keep the returned `ArrayRef` alive and borrow a typed
+/// `&Int32Array` view from it (via `as_int32_array`) rather than cloning the
+/// typed array.
+fn int32_array_arg(arg: &ColumnarValue, num_iterations: usize) -> Result<ArrayRef> {
+ arg.clone()
+ .cast_to(&DataType::Int32, None)?
+ .into_array(num_iterations)
+}
+
+/// Advance the optional band-number iterator one row, yielding the 1-based band
+/// to sample. A missing band argument defaults to band 1; a NULL band element
+/// returns `None`, which the caller propagates to a NULL result. Band 0 and
+/// negative values map to 0 so [`Bands::band`](sedona_raster::traits::Bands::band)
+/// rejects them as not 1-based rather than being silently coerced.
+fn next_band(
+ band_iter: &mut Option<arrow_array::iterator::ArrayIter<&arrow_array::Int32Array>>,
+) -> Option<usize> {
+ match band_iter.as_mut() {
+ None => Some(1),
+ Some(iter) => iter.next().flatten().map(|b| b.max(0) as usize),
+ }
+}
+
+/// Reproject `point_wkb` from its CRS into the raster CRS, returning the
+/// transformed WKB only when a reprojection actually happened (so the caller
+/// can sample the original bytes otherwise — no allocation in the common case).
+///
+/// Errors if exactly one of the point / raster carries a CRS: sampling across a
+/// known and an unknown CRS would silently mislocate the point.
+///
+/// Unlike the spatial predicates (`RS_Intersects` et al.), which fall back to a
+/// WGS84 pivot when a direct transform between two CRSes fails, a failed
+/// transform here is propagated as an error. Sampling has to land the point in
+/// the raster's own CRS — that is the only space its affine/pixel grid is
+/// defined in — so there is no neutral CRS to fall back to: a WGS84 pivot would
+/// silently sample the wrong pixel rather than compare geometries in a shared
+/// space.
+fn reproject_point(
+ point_wkb: &[u8],
+ point_crs: CrsRef<'_>,
+ raster_crs: CrsRef<'_>,
+ engine: &dyn sedona_geometry::transform::CrsEngine,
+) -> Result<Option<Vec<u8>>> {
+ match (point_crs, raster_crs) {
+ (Some(point_crs), Some(raster_crs)) => {
+ if point_crs.crs_equals(raster_crs) {
+ Ok(None)
+ } else {
+ Ok(Some(crs_transform_wkb(
+ point_wkb, point_crs, raster_crs, engine,
+ )?))
+ }
+ }
+ (None, None) => Ok(None),
+ (Some(_), None) => {
+ exec_err!("RS_Value: point has a CRS but the raster does not")
+ }
+ (None, Some(_)) => {
+ exec_err!("RS_Value: raster has a CRS but the point does not")
+ }
+ }
+}
+
+/// For a **column-level** point CRS, decide once whether sampling needs a
+/// reprojection into the raster CRS:
+/// - `Some(true)` — CRSes differ, reproject each point;
+/// - `Some(false)` — CRSes match (or both absent), sample original coordinates;
+/// - `None` — the point CRS is carried per row, so decide per point.
+///
+/// Errors when exactly one side carries a CRS. Hoisting this out of the per-point
+/// loop avoids a per-point `crs_equals`, whose `to_authority_code()` allocates a
+/// `String` every call (~15 ns/point, uniform across CRS types).
+fn column_reproject_decision(
+ point_type: &SedonaType,
+ raster_crs: CrsRef<'_>,
+) -> Result<Option<bool>> {
+ let point_crs = match point_type {
+ SedonaType::Wkb(_, c) | SedonaType::WkbView(_, c) => c.as_deref(),
+ // A per-item CRS varies by row; the caller decides per point.
+ _ => return Ok(None),
+ };
+ match (point_crs, raster_crs) {
+ (Some(p), Some(r)) => Ok(Some(!p.crs_equals(r))),
+ (None, None) => Ok(Some(false)),
+ (Some(_), None) => exec_err!("RS_Value: point has a CRS but the raster does not"),
+ (None, Some(_)) => exec_err!("RS_Value: raster has a CRS but the point does not"),
+ }
+}
+
+/// Parse a Point WKB and return its `(x, y)` in the raster's CRS, or `None` when
+/// there is nothing to sample (the point is empty — both ordinates NaN).
+///
+/// `skip_reproject` is the hoisted column-level decision: when `true` the
+/// original coordinates are returned without consulting [`reproject_point`] (and
+/// its per-point `crs_equals`); when `false` reprojection is decided per point.
+fn resolve_point_xy(
+ point_wkb: &[u8],
+ point_crs: CrsRef<'_>,
+ raster_crs: CrsRef<'_>,
+ skip_reproject: bool,
+ engine: &dyn sedona_geometry::transform::CrsEngine,
+) -> Result<Option<(f64, f64)>> {
+ let Some((px, py)) =
+ read_point_xy(point_wkb).map_err(|e| exec_datafusion_err!("RS_Value: {e}"))?
+ else {
+ return Ok(None);
+ };
+ if skip_reproject {
+ return Ok(Some((px, py)));
+ }
+ // A reprojection only happens when both sides carry a (differing) CRS;
+ // otherwise the original coordinates are sampled directly.
+ match reproject_point(point_wkb, point_crs, raster_crs, engine)? {
+ Some(reprojected) => {
+ let xy = read_point_xy(&reprojected)
+ .map_err(|e| exec_datafusion_err!("RS_Value: {e}"))?
+ .ok_or_else(|| exec_datafusion_err!("RS_Value: reprojected point is empty"))?;
+ Ok(Some(xy))
+ }
+ None => Ok(Some((px, py))),
+ }
+}
+
+/// Map a coordinate in the raster's CRS to a 0-based `(col, row)` pixel index,
+/// or `None` when the point has no location to sample.
+///
+/// A non-finite coordinate (a Point with a NaN/Inf ordinate, e.g. `POINT(NaN 5)`)
+/// returns `None`: without this guard a NaN would survive `inv_transform` and the
+/// saturating `f64 -> i64` cast would turn it into 0 (in bounds), silently
+/// sampling pixel column 0 rather than yielding NULL.
+fn xy_to_pixel(affine: &AffineMatrix, x: f64, y: f64) -> Result<Option<(i64, i64)>> {
+ if !x.is_finite() || !y.is_finite() {
+ return Ok(None);
+ }
+ let (raster_x, raster_y) = affine
+ .inv_transform(x, y)
+ .map_err(|e| exec_datafusion_err!("RS_Value: {e}"))?;
+ // Floor (not truncate toward zero) so a point just outside the top/left edge
+ // maps to a negative index and is rejected as out of bounds, rather than
+ // truncating to 0 and sampling an edge pixel. The `f64 -> i64` cast saturates
+ // (never panics); the bounds check in `read_pixel`/`sample_pixel` rejects an
+ // out-of-range index as NULL.
+ Ok(Some((raster_x.floor() as i64, raster_y.floor() as i64)))
+}
+
+/// Sample band `band_num` (1-based) at 0-based pixel `(col, row)` as `f64`.
+///
+/// Returns `None` when the pixel is out of bounds or equals the band's nodata.
+/// Reads exactly one pixel by computing its byte offset from the band's
+/// [`NdBuffer`](sedona_raster::traits::NdBuffer) strides — zero-copy and O(1),
+/// no whole-band materialisation. Errors if the band index is out of range or
+/// the band is not 2-D.
+fn sample_pixel(
+ raster: &dyn RasterRef,
+ col: i64,
+ row: i64,
+ band_num: usize,
+) -> Result<Option<f64>> {
+ let band = raster
+ .bands()
+ .band(band_num)
+ .map_err(|e| exec_datafusion_err!("RS_Value: {e}"))?;
+
+ // 2-D only: the band must be a recognized spatial (y, x) grid, not just any
+ // two-axis band (e.g. (time, band) would have len 2 but no spatial meaning).
+ if !band.is_spatial_2d() {
+ return exec_err!("RS_Value supports 2-D rasters only; band is not a 2-D (y, x) grid");
+ }
+ let buffer = band
+ .nd_buffer()
+ .map_err(|e| exec_datafusion_err!("RS_Value: {e}"))?;
+ let nodata = band
+ .nodata_as_f64()
+ .map_err(|e| exec_datafusion_err!("RS_Value: {e}"))?;
+ read_pixel(&buffer, nodata, col, row)
+}
+
+/// Read pixel `(col, row)` from an already-resolved band buffer and nodata
+/// value. Returns `None` for an out-of-bounds pixel or one that equals nodata.
+/// Hoisting the band resolution out of this function lets callers sample many
+/// points from one buffer without re-resolving per point.
+fn read_pixel(buffer: &NdBuffer, nodata: Option<f64>, col: i64, row: i64) -> Result<Option<f64>> {
+ let (height, width) = (buffer.shape[0], buffer.shape[1]);
+ if row < 0 || row >= height || col < 0 || col >= width {
+ return Ok(None);
+ }
+
+ // Byte offset of the (row, col) pixel via the band's own strides, so the
+ // read stays correct for any layout the producer hands us. Checked
+ // arithmetic throughout: `row`/`col` are already in bounds, but a corrupt
+ // stride or offset must surface as an error, never an i64 overflow panic.
+ let size = buffer.data_type.byte_size() as i64;
+ let byte_offset = row
+ .checked_mul(buffer.strides[0])
+ .zip(col.checked_mul(buffer.strides[1]))
+ .and_then(|(r, c)| r.checked_add(c))
+ .and_then(|rc| rc.checked_add(buffer.offset as i64))
+ .ok_or_else(|| exec_datafusion_err!("RS_Value: pixel byte offset overflow"))?;
+ let end_offset = byte_offset
+ .checked_add(size)
+ .ok_or_else(|| exec_datafusion_err!("RS_Value: pixel byte offset overflow"))?;
+ let start = usize::try_from(byte_offset)
+ .map_err(|_| exec_datafusion_err!("RS_Value: negative pixel byte offset"))?;
+ let end = usize::try_from(end_offset)
+ .map_err(|_| exec_datafusion_err!("RS_Value: pixel byte offset overflow"))?;
+ let bytes = buffer.buffer.get(start..end).ok_or_else(|| {
+ exec_datafusion_err!("RS_Value: pixel is out of the band's buffer bounds")
+ })?;
+
+ // Decode the pixel to f64. The lossless converter errors (rather than
+ // silently rounding) on Int64/UInt64 values beyond f64's exact-integer
+ // range (2^53) — RS_Value returns a Double, so such a pixel can't be
+ // represented faithfully; failing loudly is preferred over a wrong value.
+ let value = nodata_bytes_to_f64_lossless(bytes, &buffer.data_type)
+ .map_err(|e| exec_datafusion_err!("RS_Value: {e}"))?;
+
+ if let Some(nodata) = nodata {
+ if value == nodata || (value.is_nan() && nodata.is_nan()) {
+ return Ok(None);
+ }
+ }
+
+ Ok(Some(value))
+}
+
+#[cfg(test)]
+mod tests {
+ use super::*;
+ use arrow_array::{Array, Float64Array, Int32Array};
+ use datafusion_expr::ScalarUDF;
+ use sedona_raster::array::RasterStructArray;
+ use sedona_schema::crs::lnglat;
+ use sedona_schema::datatypes::{Edges, RASTER};
+ use sedona_schema::raster::BandDataType;
+ use sedona_testing::create::create_array as create_geom_array;
+ use sedona_testing::raster_spec::RasterSpec;
+ use sedona_testing::rasters::generate_test_rasters;
+ use sedona_testing::testers::ScalarUdfTester;
+
+ /// Resolve a single `RasterRefImpl` from a one-row spec for direct
+ /// `sample_pixel` exercise.
+ fn sample(spec: RasterSpec, col: i64, row: i64, band: usize) -> Result<Option<f64>> {
+ let array = spec.build();
+ let rasters = RasterStructArray::try_new(&array).unwrap();
+ let raster = rasters.get(0).unwrap();
+ sample_pixel(&raster, col, row, band)
+ }
+
+ #[test]
+ fn udf_metadata() {
+ let udf: ScalarUDF = rs_value_udf().into();
+ assert_eq!(udf.name(), "rs_value");
+ }
+
+ #[test]
+ fn udf_marks_needs_pixels() {
+ assert_eq!(
+ rs_value_udf()
+ .metadata()
+ .get(NEEDS_PIXELS_METADATA_KEY)
+ .map(String::as_str),
+ Some("true")
+ );
+ }
+
+ #[test]
+ fn return_type_is_float64() {
+ // (raster, point) resolves to a Float64 output.
+ let return_type = RsValuePoint { with_band: false }
+ .return_type(&[RASTER, SedonaType::Wkb(Edges::Planar, lnglat())])
+ .unwrap();
+ assert_eq!(return_type, Some(SedonaType::Arrow(DataType::Float64)));
+ }
+
+ #[test]
+ fn samples_2d_pixels_row_major() {
+ // 3x2 raster, row-major pixels:
+ // row0 = [10, 20, 30], row1 = [40, 50, 60]
+ let spec = || RasterSpec::d2(3, 2).band_values(&[10u8, 20, 30, 40, 50, 60]);
+ assert_eq!(sample(spec(), 0, 0, 1).unwrap(), Some(10.0)); // top-left
+ assert_eq!(sample(spec(), 2, 0, 1).unwrap(), Some(30.0)); // top-right
+ assert_eq!(sample(spec(), 0, 1, 1).unwrap(), Some(40.0)); // bottom-left
+ assert_eq!(sample(spec(), 2, 1, 1).unwrap(), Some(60.0)); // bottom-right
+ }
+
+ #[test]
+ fn out_of_bounds_pixel_is_none() {
+ let spec = || RasterSpec::d2(3, 2).band_values(&[10u8, 20, 30, 40, 50, 60]);
+ assert_eq!(sample(spec(), 3, 0, 1).unwrap(), None); // col == width
+ assert_eq!(sample(spec(), 0, 2, 1).unwrap(), None); // row == height
+ assert_eq!(sample(spec(), -1, 0, 1).unwrap(), None); // negative
+ }
+
+ #[test]
+ fn nodata_pixel_is_none() {
+ let spec = RasterSpec::d2(2, 1).band_values(&[7u8, 9]).nodata(9u8);
+ assert_eq!(sample(spec.clone(), 0, 0, 1).unwrap(), Some(7.0));
+ assert_eq!(sample(spec, 1, 0, 1).unwrap(), None);
+ }
+
+ #[test]
+ fn second_band_is_addressable() {
+ let spec = RasterSpec::d2(2, 1)
+ .band_values(&[1u8, 2])
+ .band_values(&[30u8, 40]);
+ assert_eq!(sample(spec.clone(), 1, 0, 1).unwrap(), Some(2.0));
+ assert_eq!(sample(spec, 1, 0, 2).unwrap(), Some(40.0));
+ }
+
+ #[test]
+ fn float_band_values_round_trip() {
+ let spec = RasterSpec::d2(2, 1).band_values(&[1.5f32, -2.5]);
+ assert_eq!(sample(spec.clone(), 0, 0, 1).unwrap(), Some(1.5));
+ assert_eq!(sample(spec, 1, 0, 1).unwrap(), Some(-2.5));
+ }
+
+ #[test]
+ fn band_out_of_range_errors() {
+ let spec = RasterSpec::d2(2, 1).band_values(&[1u8, 2]);
+ let err = sample(spec, 0, 0, 2).unwrap_err().to_string();
+ assert!(err.contains("RS_Value"), "unexpected error: {err}");
+ }
+
+ #[test]
+ fn band_zero_errors() {
+ // Band 0 is not coerced to band 1 — it surfaces as a 1-based error.
+ let spec = RasterSpec::d2(2, 1).band_values(&[1u8, 2]);
+ let err = sample(spec, 0, 0, 0).unwrap_err().to_string();
+ assert!(err.contains("1-based"), "unexpected error: {err}");
+ }
+
+ #[test]
+ fn nan_nodata_pixel_is_none() {
+ // A float band whose nodata is NaN: a NaN pixel reads as NULL (NaN != NaN
+ // makes the `==` check insufficient), a normal pixel reads as its value.
+ let spec = RasterSpec::d2(2, 1)
+ .band_values(&[f64::NAN, 1.0])
+ .nodata(f64::NAN);
+ assert_eq!(sample(spec.clone(), 0, 0, 1).unwrap(), None);
+ assert_eq!(sample(spec, 1, 0, 1).unwrap(), Some(1.0));
+ }
+
+ #[test]
+ fn non_2d_band_errors() {
+ // A band with a leading non-spatial dimension is rejected.
+ let spec = RasterSpec::nd(&["time", "y", "x"], &[2, 2, 1]).band(BandDataType::UInt8);
+ let err = sample(spec, 0, 0, 1).unwrap_err().to_string();
+ assert!(err.contains("2-D"), "unexpected error: {err}");
+ }
+
+ #[test]
+ fn point_crs_mismatch_errors() {
+ let udf: ScalarUDF = rs_value_udf().into();
+
+ // Raster has a CRS (generate_test_rasters sets OGC:CRS84), point does not.
+ let geom_type = SedonaType::Wkb(Edges::Planar, None);
+ let tester = ScalarUdfTester::new(udf.clone(), vec![RASTER, geom_type.clone()]);
+ let rasters = generate_test_rasters(1, None).unwrap();
+ let geoms = create_geom_array(&[Some("POINT (2.1 2.6)")], &geom_type);
+ let err = tester
+ .invoke_arrays(vec![Arc::new(rasters), geoms])
+ .unwrap_err()
+ .to_string();
+ assert!(
+ err.contains("raster has a CRS but the point does not"),
+ "unexpected error: {err}"
+ );
+
+ // Point has a CRS, raster does not.
+ let geom_type = SedonaType::Wkb(Edges::Planar, lnglat());
+ let tester = ScalarUdfTester::new(udf, vec![RASTER, geom_type.clone()]);
+ let rasters = RasterSpec::d2(2, 2)
+ .band(BandDataType::UInt8)
+ .crs(None)
+ .build();
+ let geoms = create_geom_array(&[Some("POINT (0 0)")], &geom_type);
+ let err = tester
+ .invoke_arrays(vec![Arc::new(rasters), geoms])
+ .unwrap_err()
+ .to_string();
+ assert!(
+ err.contains("point has a CRS but the raster does not"),
+ "unexpected error: {err}"
+ );
+ }
+
+ #[test]
+ fn non_point_geometry_errors() {
+ let udf: ScalarUDF = rs_value_udf().into();
+ let geom_type = SedonaType::Wkb(Edges::Planar, lnglat());
+ let tester = ScalarUdfTester::new(udf, vec![RASTER, geom_type.clone()]);
+ let rasters = generate_test_rasters(1, None).unwrap();
+ let geoms = create_geom_array(&[Some("LINESTRING (0 0, 1 1)")], &geom_type);
+ let err = tester
+ .invoke_arrays(vec![Arc::new(rasters), geoms])
+ .unwrap_err()
+ .to_string();
+ assert!(err.contains("expected a Point"), "unexpected error: {err}");
+ }
+
+ #[test]
+ fn empty_point_is_none() {
+ // POINT EMPTY has no location to sample, so the result is NULL rather
+ // than an error. The empty check short-circuits before CRS resolution,
+ // so a missing/again-mismatched point CRS does not matter here.
+ let udf: ScalarUDF = rs_value_udf().into();
+ let geom_type = SedonaType::Wkb(Edges::Planar, lnglat());
+ let tester = ScalarUdfTester::new(udf, vec![RASTER, geom_type.clone()]);
+ let rasters = generate_test_rasters(1, None).unwrap();
+ let geoms = create_geom_array(&[Some("POINT EMPTY")], &geom_type);
+ let result = tester
+ .invoke_arrays(vec![Arc::new(rasters), geoms])
+ .unwrap();
+ let arr = result.as_any().downcast_ref::<Float64Array>().unwrap();
+ assert!(arr.is_null(0), "POINT EMPTY should sample to NULL");
+ }
+
+ #[test]
+ fn point_just_outside_top_edge_is_none() {
+ // North-up raster: origin (0, 10), 1x1 pixels (geotransform
+ // [c, a, b, f, d, e] = [0, 1, 0, 10, 0, -1]), so world y decreases down
+ // the rows. A point at y = 10.5 is just *above* the top edge: its pixel
+ // row is -0.5, which must floor to -1 (out of bounds -> NULL), not
+ // truncate toward zero to 0 (the top row).
+ let udf: ScalarUDF = rs_value_udf().into();
+ let geom_type = SedonaType::Wkb(Edges::Planar, lnglat());
+ let tester = ScalarUdfTester::new(udf, vec![RASTER, geom_type.clone()]);
+ let raster = || {
+ RasterSpec::d2(2, 2)
+ .band_values(&[1u8, 2, 3, 4])
+ .transform([0.0, 1.0, 0.0, 10.0, 0.0, -1.0])
+ .build()
+ };
+
+ // Just above the top edge -> NULL.
+ let geoms = create_geom_array(&[Some("POINT (0.5 10.5)")], &geom_type);
+ let result = tester
+ .invoke_arrays(vec![Arc::new(raster()), geoms])
+ .unwrap();
+ let arr = result.as_any().downcast_ref::<Float64Array>().unwrap();
+ assert!(arr.is_null(0), "point above the top edge should be NULL");
+
+ // Just inside the top row -> the top-left value (1).
+ let geoms = create_geom_array(&[Some("POINT (0.5 9.5)")], &geom_type);
+ let result = tester
+ .invoke_arrays(vec![Arc::new(raster()), geoms])
+ .unwrap();
+ let arr = result.as_any().downcast_ref::<Float64Array>().unwrap();
+ assert_eq!(arr.value(0), 1.0);
+ }
+
+ #[test]
+ fn scalar_raster_samples_via_fast_path() {
+ // A constant (scalar) raster with the default band takes the optimized
+ // hoisted path; verify it samples the right pixel and rejects
+ // out-of-bounds, matching the general (array-raster) path.
+ let udf: ScalarUDF = rs_value_udf().into();
+ let geom_type = SedonaType::Wkb(Edges::Planar, lnglat());
+ let tester = ScalarUdfTester::new(udf, vec![RASTER, geom_type.clone()]);
+
+ // North-up 2x2 raster, origin (0, 10), 1x1 pixels; row-major [1,2 / 3,4].
+ let raster = RasterSpec::d2(2, 2)
+ .band_values(&[1u8, 2, 3, 4])
+ .transform([0.0, 1.0, 0.0, 10.0, 0.0, -1.0])
+ .build();
+
+ let sample = |wkt: &str| -> Option<f64> {
+ let geoms = create_geom_array(&[Some(wkt)], &geom_type);
+ let result = tester
+ .invoke(vec![
+ ColumnarValue::Scalar(ScalarValue::Struct(Arc::new(raster.clone()))),
+ ColumnarValue::Array(geoms),
+ ])
+ .unwrap();
+ let arr = result.into_array(1).unwrap();
+ let arr = arr.as_any().downcast_ref::<Float64Array>().unwrap();
+ (!arr.is_null(0)).then(|| arr.value(0))
+ };
+
+ assert_eq!(sample("POINT (0.5 9.5)"), Some(1.0)); // pixel (0, 0)
+ assert_eq!(sample("POINT (1.5 8.5)"), Some(4.0)); // pixel (1, 1)
+ assert_eq!(sample("POINT (100 100)"), None); // outside the footprint
+ }
+
+ #[test]
+ fn scalar_raster_constant_band_uses_fast_path() {
+ // A scalar raster with a constant (scalar) band argument takes the same
+ // hoisted fast path as the 2-arg default-band case, just on that band.
+ let udf: ScalarUDF = rs_value_udf().into();
+ let geom_type = SedonaType::Wkb(Edges::Planar, lnglat());
+ let tester = ScalarUdfTester::new(
+ udf,
+ vec![
+ RASTER,
+ geom_type.clone(),
+ SedonaType::Arrow(DataType::Int32),
+ ],
+ );
+ // North-up 2x2, two bands: band 1 [1,2,3,4], band 2 [10,20,30,40].
+ let raster = RasterSpec::d2(2, 2)
+ .band_values(&[1u8, 2, 3, 4])
+ .band_values(&[10u8, 20, 30, 40])
+ .transform([0.0, 1.0, 0.0, 10.0, 0.0, -1.0])
+ .build();
+ let geoms = create_geom_array(&[Some("POINT (0.5 9.5)")], &geom_type); // pixel (0, 0)
+ let result = tester
+ .invoke(vec![
+ ColumnarValue::Scalar(ScalarValue::Struct(Arc::new(raster))),
+ ColumnarValue::Array(geoms),
+ ColumnarValue::Scalar(ScalarValue::Int32(Some(2))),
+ ])
+ .unwrap();
+ let arr = result.into_array(1).unwrap();
+ let arr = arr.as_any().downcast_ref::<Float64Array>().unwrap();
+ assert_eq!(arr.value(0), 10.0); // band 2, pixel (0, 0)
+ }
+
+ #[test]
+ fn scalar_raster_band_column_resolves_per_row() {
+ // A scalar raster with a band *column* still hoists affine/CRS/reproject,
+ // but resolves the band buffer per row. A NULL band element -> NULL.
+ let udf: ScalarUDF = rs_value_udf().into();
+ let geom_type = SedonaType::Wkb(Edges::Planar, lnglat());
+ let tester = ScalarUdfTester::new(
+ udf,
+ vec![
+ RASTER,
+ geom_type.clone(),
+ SedonaType::Arrow(DataType::Int32),
+ ],
+ );
+ let raster = RasterSpec::d2(2, 2)
+ .band_values(&[1u8, 2, 3, 4])
+ .band_values(&[10u8, 20, 30, 40])
+ .transform([0.0, 1.0, 0.0, 10.0, 0.0, -1.0])
+ .build();
+ // Three points at pixel (0, 0), sampling band 1, band 2, then a NULL band.
+ let geoms = create_geom_array(
+ &[
+ Some("POINT (0.5 9.5)"),
+ Some("POINT (0.5 9.5)"),
+ Some("POINT (0.5 9.5)"),
+ ],
+ &geom_type,
+ );
+ let bands = Arc::new(Int32Array::from(vec![Some(1), Some(2), None]));
+ let result = tester
+ .invoke(vec![
+ ColumnarValue::Scalar(ScalarValue::Struct(Arc::new(raster))),
+ ColumnarValue::Array(geoms),
+ ColumnarValue::Array(bands),
+ ])
+ .unwrap();
+ let arr = result.into_array(3).unwrap();
+ let arr = arr.as_any().downcast_ref::<Float64Array>().unwrap();
+ assert_eq!(arr.value(0), 1.0); // band 1
+ assert_eq!(arr.value(1), 10.0); // band 2
+ assert!(arr.is_null(2), "NULL band element should be NULL");
+ }
+
+ #[test]
+ fn scalar_raster_null_scalar_band_is_all_null() {
+ // A NULL scalar band makes every output NULL without touching the band.
+ let udf: ScalarUDF = rs_value_udf().into();
+ let geom_type = SedonaType::Wkb(Edges::Planar, lnglat());
+ let tester = ScalarUdfTester::new(
+ udf,
+ vec![
+ RASTER,
+ geom_type.clone(),
+ SedonaType::Arrow(DataType::Int32),
+ ],
+ );
+ let raster = RasterSpec::d2(2, 2).band_values(&[1u8, 2, 3, 4]).build();
+ let geoms = create_geom_array(&[Some("POINT (0.5 0.5)")], &geom_type);
+ let result = tester
+ .invoke(vec![
+ ColumnarValue::Scalar(ScalarValue::Struct(Arc::new(raster))),
+ ColumnarValue::Array(geoms),
+ ColumnarValue::Scalar(ScalarValue::Int32(None)),
+ ])
+ .unwrap();
+ let arr = result.into_array(1).unwrap();
+ let arr = arr.as_any().downcast_ref::<Float64Array>().unwrap();
+ assert!(arr.is_null(0), "NULL scalar band should sample to NULL");
+ }
+
+ #[test]
+ fn non_finite_point_ordinate_is_none() {
+ // A point with a NaN/Inf ordinate (e.g. POINT(NaN 5)) has no location to
+ // sample. Without the finite guard a NaN survives `inv_transform` and the
+ // saturating cast turns it into pixel column 0, silently sampling a value.
+ let raster = RasterSpec::d2(2, 2)
+ .band_values(&[1u8, 2, 3, 4])
+ .transform([0.0, 1.0, 0.0, 10.0, 0.0, -1.0])
+ .build();
+ let rasters = RasterStructArray::try_new(&raster).unwrap();
+ let affine = AffineMatrix::from_metadata(&rasters.get(0).unwrap().metadata());
+
+ assert_eq!(xy_to_pixel(&affine, f64::NAN, 5.0).unwrap(), None);
+ assert_eq!(xy_to_pixel(&affine, 5.0, f64::NAN).unwrap(), None);
+ assert_eq!(xy_to_pixel(&affine, f64::INFINITY, 5.0).unwrap(), None);
+ // A finite in-bounds point still maps to a real pixel.
+ assert_eq!(xy_to_pixel(&affine, 0.5, 9.5).unwrap(), Some((0, 0)));
+ }
+
+ #[test]
+ fn scalar_all_null_points_defer_band_resolution() {
+ // The scalar fast path defers band/2-D validation until a point needs
+ // sampling: an all-null point column over an unsupported (non-2-D) raster
+ // returns all-NULL rather than erroring, matching the general path. A
+ // valid point over the same raster still surfaces the 2-D error.
+ let udf: ScalarUDF = rs_value_udf().into();
+ let geom_type = SedonaType::Wkb(Edges::Planar, None);
+ let tester = ScalarUdfTester::new(udf, vec![RASTER, geom_type.clone()]);
+ let raster = RasterSpec::nd(&["time", "y", "x"], &[2, 2, 1])
+ .band(BandDataType::UInt8)
+ .crs(None)
+ .build();
+
+ let invoke = |wkt: Option<&str>| {
+ let geoms = create_geom_array(&[wkt], &geom_type);
+ tester.invoke(vec![
+ ColumnarValue::Scalar(ScalarValue::Struct(Arc::new(raster.clone()))),
+ ColumnarValue::Array(geoms),
+ ])
+ };
+
+ // All-null point -> NULL, no band resolution, no error.
+ let result = invoke(None).unwrap();
+ let arr = result.into_array(1).unwrap();
+ let arr = arr.as_any().downcast_ref::<Float64Array>().unwrap();
+ assert!(
+ arr.is_null(0),
+ "all-null point over a non-2-D raster should be NULL"
+ );
+
+ // A real point forces band resolution, which rejects the non-2-D raster.
+ let err = invoke(Some("POINT (0 0)")).unwrap_err().to_string();
+ assert!(err.contains("2-D"), "unexpected error: {err}");
+ }
+
+ #[test]
+ fn crs_decision_equal_crs_skips_reproject() {
+ // The common case: a lng/lat point CRS and a lng/lat raster are detected
+ // as equal, so the per-point reproject (and its crs_equals) is skipped.
+ // This is what the B1 hoist relies on — if it returned Some(true) the
+ // optimization would silently no-op.
+ let raster = RasterSpec::d2(2, 2).band(BandDataType::UInt8).build(); // default lng/lat
+ let rasters = RasterStructArray::try_new(&raster).unwrap();
+ let raster_crs = resolve_crs(rasters.get(0).unwrap().crs()).unwrap();
+ let point_type = SedonaType::Wkb(Edges::Planar, lnglat());
+ assert_eq!(
+ column_reproject_decision(&point_type, raster_crs.as_deref()).unwrap(),
+ Some(false),
+ "lng/lat point + lng/lat raster must be detected as equal (skip reproject)"
+ );
+ }
+
+ #[test]
+ fn crs_decision_differing_crs_reprojects() {
+ use sedona_schema::crs::deserialize_crs;
+ let raster = RasterSpec::d2(2, 2)
+ .crs(Some("EPSG:4326"))
+ .band(BandDataType::UInt8)
+ .build();
+ let rasters = RasterStructArray::try_new(&raster).unwrap();
+ let raster_crs = resolve_crs(rasters.get(0).unwrap().crs()).unwrap();
+ let point_type = SedonaType::Wkb(Edges::Planar, deserialize_crs("EPSG:3857").unwrap());
+ assert_eq!(
+ column_reproject_decision(&point_type, raster_crs.as_deref()).unwrap(),
+ Some(true),
+ "EPSG:3857 point + EPSG:4326 raster must require reprojection"
+ );
+ }
+
+ #[test]
+ fn crs_decision_one_sided_crs_errors() {
+ let raster = RasterSpec::d2(2, 2)
+ .crs(None)
+ .band(BandDataType::UInt8)
+ .build();
+ let rasters = RasterStructArray::try_new(&raster).unwrap();
+ let raster_crs = resolve_crs(rasters.get(0).unwrap().crs()).unwrap();
+ let point_type = SedonaType::Wkb(Edges::Planar, lnglat());
+ assert!(column_reproject_decision(&point_type, raster_crs.as_deref()).is_err());
+ }
+}