paleolimbot commented on code in PR #1027:
URL: https://github.com/apache/sedona-db/pull/1027#discussion_r3600359175
##########
python/sedonadb/tests/functions/test_raster_functions.py:
##########
@@ -283,6 +292,183 @@ def
test_rs_setbandnodatavalue_two_arg_requires_single_band():
)
+# RS_Values samples one pixel per sub-point and returns a List<Double> in input
+# order. Same raster facts as `test_rs_value_point`: band `b` is filled with
`b`,
+# (74.58, 110.57) is the centroid of pixel (10, 10), and (44.58, 80.57) is the
+# centroid of the top-left pixel (0, 0), which is set to nodata. A point far
+# outside the footprint yields a NULL element; the whole result is never NULL
+# here because the geometry is non-null.
[email protected](
+ ("expr", "expected"),
+ [
+ # Two in-bounds points + one outside, on band 1.
+ (
+ "RS_Values(RS_Example(), ST_GeomFromText('MULTIPOINT (74.58
110.57, 74.58 110.57, 0 0)', 'OGC:CRS84'), 1)",
+ [1.0, 1.0, None],
+ ),
+ # Explicit band selects the plane; nodata corner and outside are NULL.
+ (
+ "RS_Values(RS_Example(), ST_GeomFromText('MULTIPOINT (74.58
110.57, 44.58 80.57, 0 0)', 'OGC:CRS84'), 2)",
+ [2.0, None, None],
+ ),
+ (
+ "RS_Values(RS_Example(), ST_GeomFromText('MULTIPOINT (74.58
110.57)', 'OGC:CRS84'), 3)",
+ [3.0],
+ ),
+ # A bare Point is accepted and yields a one-element list.
+ (
+ "RS_Values(RS_Example(), ST_Point(74.58, 110.57, 'OGC:CRS84'), 1)",
+ [1.0],
+ ),
+ # An empty MultiPoint is an empty list (not NULL).
+ (
+ "RS_Values(RS_Example(), ST_GeomFromText('MULTIPOINT EMPTY',
'OGC:CRS84'), 1)",
+ [],
+ ),
+ ],
+)
+def test_rs_values_multipoint(expr, expected):
+ SedonaDB().assert_query_result(f"SELECT {expr}", [(expected,)])
+
+
+def test_rs_values_default_band_requires_single_band(con):
+ # RS_Example has 3 bands, so omitting the band is ambiguous and errors.
+ with pytest.raises(Exception, match="specify which band"):
+ con.sql(
+ "SELECT RS_Values(RS_Example(), ST_GeomFromText('MULTIPOINT (74.58
110.57)', 'OGC:CRS84'))"
+ ).to_arrow_table()
+
+
+def test_rs_values_ensureloaded_outdb(con, sedona_testing):
+ """RS_Values over an OutDb raster exercises the needs_pixels ->
+ RS_EnsureLoaded planner path end to end: the raster from RS_FromPath
carries
+ no pixels, so the planner must materialise it before the kernel samples.
+
+ sentinel2.tif's top-left pixel holds 2324 (see test_rs_ensureloaded); its
+ world center is derived from the raster's own georeference so the test does
+ not hard-code the file's geotransform.
+ """
+ path = sedona_testing / "data/raster/sentinel2.tif"
+ t = con.sql("SELECT RS_FromPath($1) AS raster", params=(str(path),))
+ view = "test_rs_values_ensureloaded_outdb_raster"
+ t.to_view(view)
+ try:
+ # `.to_pylist()` converts the whole (one-row) table in one pass;
+ # indexing a column with `[0]` first would force a chunk-combining
+ # copy in pyarrow.
+ meta = (
+ con.sql(
+ f"SELECT RS_GeoReference(raster) AS georef, RS_SRID(raster) AS
srid FROM {view}"
+ )
+ .to_arrow_table()
+ .to_pylist()[0]
+ )
+ georef = [float(v) for v in meta["georef"].split()]
+ scale_x, skew_y, skew_x, scale_y, ul_x, ul_y = georef
+ srid = meta["srid"]
+ # World center of pixel (0, 0): upper-left corner + half a pixel step.
+ cx = ul_x + 0.5 * scale_x + 0.5 * skew_x
+ cy = ul_y + 0.5 * skew_y + 0.5 * scale_y
+
+ values = (
+ con.sql(
+ f"""
+ SELECT RS_Values(
+ raster,
+ ST_GeomFromText('MULTIPOINT ({cx} {cy})', 'EPSG:{srid}'),
+ 1
+ ) AS v FROM {view}
+ """
+ )
+ .to_arrow_table()["v"]
+ .to_pylist()
+ )
+ assert values == [[2324.0]]
+ finally:
+ con.drop_view(view)
+
+
+def test_rs_values_matches_rasterio(con):
+ """Cross-check RS_Values against rasterio on a random raster.
+
+ The plural counterpart of `test_rs_value_matches_rasterio`: the same dense
+ set of sample points is passed as a single MultiPoint, so one `RS_Values`
+ call returns a list that must match rasterio's per-point reads element for
+ element, in order.
+ """
+ import numpy as np
+
+ 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): every pixel center
plus
+ # four off-center positions (kept inside the pixel to avoid floor
ambiguity),
+ # then a batch of random interior points.
+ 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 every point in one MultiPoint via a single RS_Values
call.
+ raster = Raster.from_numpy(data, transform=gdal_transform)
+ wkt = "MULTIPOINT (" + ", ".join(f"{x} {y}" for x, y in zip(xs, ys)) + ")"
+ got = (
+ con.sql(
+ "SELECT RS_Values($1, ST_GeomFromText($2)) AS v",
+ params=(raster, wkt),
+ )
+ .to_arrow_table()["v"]
+ .to_pylist()[0]
+ )
+
+ assert got == pytest.approx(expected)
Review Comment:
Does this need to be approximated or is the difference here an indicator of
different behaviour?
##########
python/sedonadb/tests/functions/test_raster_functions.py:
##########
@@ -283,6 +292,183 @@ def
test_rs_setbandnodatavalue_two_arg_requires_single_band():
)
+# RS_Values samples one pixel per sub-point and returns a List<Double> in input
+# order. Same raster facts as `test_rs_value_point`: band `b` is filled with
`b`,
+# (74.58, 110.57) is the centroid of pixel (10, 10), and (44.58, 80.57) is the
+# centroid of the top-left pixel (0, 0), which is set to nodata. A point far
+# outside the footprint yields a NULL element; the whole result is never NULL
+# here because the geometry is non-null.
[email protected](
+ ("expr", "expected"),
+ [
+ # Two in-bounds points + one outside, on band 1.
+ (
+ "RS_Values(RS_Example(), ST_GeomFromText('MULTIPOINT (74.58
110.57, 74.58 110.57, 0 0)', 'OGC:CRS84'), 1)",
+ [1.0, 1.0, None],
+ ),
+ # Explicit band selects the plane; nodata corner and outside are NULL.
+ (
+ "RS_Values(RS_Example(), ST_GeomFromText('MULTIPOINT (74.58
110.57, 44.58 80.57, 0 0)', 'OGC:CRS84'), 2)",
+ [2.0, None, None],
+ ),
+ (
+ "RS_Values(RS_Example(), ST_GeomFromText('MULTIPOINT (74.58
110.57)', 'OGC:CRS84'), 3)",
+ [3.0],
+ ),
+ # A bare Point is accepted and yields a one-element list.
+ (
+ "RS_Values(RS_Example(), ST_Point(74.58, 110.57, 'OGC:CRS84'), 1)",
+ [1.0],
+ ),
+ # An empty MultiPoint is an empty list (not NULL).
+ (
+ "RS_Values(RS_Example(), ST_GeomFromText('MULTIPOINT EMPTY',
'OGC:CRS84'), 1)",
+ [],
+ ),
+ ],
+)
+def test_rs_values_multipoint(expr, expected):
+ SedonaDB().assert_query_result(f"SELECT {expr}", [(expected,)])
+
+
+def test_rs_values_default_band_requires_single_band(con):
+ # RS_Example has 3 bands, so omitting the band is ambiguous and errors.
+ with pytest.raises(Exception, match="specify which band"):
+ con.sql(
+ "SELECT RS_Values(RS_Example(), ST_GeomFromText('MULTIPOINT (74.58
110.57)', 'OGC:CRS84'))"
+ ).to_arrow_table()
+
+
+def test_rs_values_ensureloaded_outdb(con, sedona_testing):
+ """RS_Values over an OutDb raster exercises the needs_pixels ->
+ RS_EnsureLoaded planner path end to end: the raster from RS_FromPath
carries
+ no pixels, so the planner must materialise it before the kernel samples.
+
+ sentinel2.tif's top-left pixel holds 2324 (see test_rs_ensureloaded); its
+ world center is derived from the raster's own georeference so the test does
+ not hard-code the file's geotransform.
+ """
+ path = sedona_testing / "data/raster/sentinel2.tif"
+ t = con.sql("SELECT RS_FromPath($1) AS raster", params=(str(path),))
+ view = "test_rs_values_ensureloaded_outdb_raster"
+ t.to_view(view)
+ try:
+ # `.to_pylist()` converts the whole (one-row) table in one pass;
+ # indexing a column with `[0]` first would force a chunk-combining
+ # copy in pyarrow.
+ meta = (
+ con.sql(
+ f"SELECT RS_GeoReference(raster) AS georef, RS_SRID(raster) AS
srid FROM {view}"
+ )
+ .to_arrow_table()
+ .to_pylist()[0]
+ )
+ georef = [float(v) for v in meta["georef"].split()]
+ scale_x, skew_y, skew_x, scale_y, ul_x, ul_y = georef
+ srid = meta["srid"]
+ # World center of pixel (0, 0): upper-left corner + half a pixel step.
+ cx = ul_x + 0.5 * scale_x + 0.5 * skew_x
+ cy = ul_y + 0.5 * skew_y + 0.5 * scale_y
+
+ values = (
+ con.sql(
+ f"""
+ SELECT RS_Values(
+ raster,
+ ST_GeomFromText('MULTIPOINT ({cx} {cy})', 'EPSG:{srid}'),
+ 1
+ ) AS v FROM {view}
+ """
+ )
+ .to_arrow_table()["v"]
+ .to_pylist()
+ )
+ assert values == [[2324.0]]
+ finally:
+ con.drop_view(view)
+
+
+def test_rs_values_matches_rasterio(con):
+ """Cross-check RS_Values against rasterio on a random raster.
+
+ The plural counterpart of `test_rs_value_matches_rasterio`: the same dense
+ set of sample points is passed as a single MultiPoint, so one `RS_Values`
+ call returns a list that must match rasterio's per-point reads element for
+ element, in order.
+ """
+ import numpy as np
Review Comment:
This can be module-level
##########
python/sedonadb/tests/functions/test_raster_functions.py:
##########
@@ -283,6 +292,183 @@ def
test_rs_setbandnodatavalue_two_arg_requires_single_band():
)
+# RS_Values samples one pixel per sub-point and returns a List<Double> in input
+# order. Same raster facts as `test_rs_value_point`: band `b` is filled with
`b`,
+# (74.58, 110.57) is the centroid of pixel (10, 10), and (44.58, 80.57) is the
+# centroid of the top-left pixel (0, 0), which is set to nodata. A point far
+# outside the footprint yields a NULL element; the whole result is never NULL
+# here because the geometry is non-null.
[email protected](
+ ("expr", "expected"),
+ [
+ # Two in-bounds points + one outside, on band 1.
+ (
+ "RS_Values(RS_Example(), ST_GeomFromText('MULTIPOINT (74.58
110.57, 74.58 110.57, 0 0)', 'OGC:CRS84'), 1)",
+ [1.0, 1.0, None],
+ ),
+ # Explicit band selects the plane; nodata corner and outside are NULL.
+ (
+ "RS_Values(RS_Example(), ST_GeomFromText('MULTIPOINT (74.58
110.57, 44.58 80.57, 0 0)', 'OGC:CRS84'), 2)",
+ [2.0, None, None],
+ ),
+ (
+ "RS_Values(RS_Example(), ST_GeomFromText('MULTIPOINT (74.58
110.57)', 'OGC:CRS84'), 3)",
+ [3.0],
+ ),
+ # A bare Point is accepted and yields a one-element list.
+ (
+ "RS_Values(RS_Example(), ST_Point(74.58, 110.57, 'OGC:CRS84'), 1)",
+ [1.0],
+ ),
+ # An empty MultiPoint is an empty list (not NULL).
+ (
+ "RS_Values(RS_Example(), ST_GeomFromText('MULTIPOINT EMPTY',
'OGC:CRS84'), 1)",
+ [],
+ ),
+ ],
+)
+def test_rs_values_multipoint(expr, expected):
+ SedonaDB().assert_query_result(f"SELECT {expr}", [(expected,)])
+
+
+def test_rs_values_default_band_requires_single_band(con):
+ # RS_Example has 3 bands, so omitting the band is ambiguous and errors.
+ with pytest.raises(Exception, match="specify which band"):
+ con.sql(
+ "SELECT RS_Values(RS_Example(), ST_GeomFromText('MULTIPOINT (74.58
110.57)', 'OGC:CRS84'))"
+ ).to_arrow_table()
+
+
+def test_rs_values_ensureloaded_outdb(con, sedona_testing):
+ """RS_Values over an OutDb raster exercises the needs_pixels ->
+ RS_EnsureLoaded planner path end to end: the raster from RS_FromPath
carries
+ no pixels, so the planner must materialise it before the kernel samples.
+
+ sentinel2.tif's top-left pixel holds 2324 (see test_rs_ensureloaded); its
+ world center is derived from the raster's own georeference so the test does
+ not hard-code the file's geotransform.
+ """
+ path = sedona_testing / "data/raster/sentinel2.tif"
+ t = con.sql("SELECT RS_FromPath($1) AS raster", params=(str(path),))
+ view = "test_rs_values_ensureloaded_outdb_raster"
+ t.to_view(view)
+ try:
+ # `.to_pylist()` converts the whole (one-row) table in one pass;
+ # indexing a column with `[0]` first would force a chunk-combining
+ # copy in pyarrow.
+ meta = (
+ con.sql(
+ f"SELECT RS_GeoReference(raster) AS georef, RS_SRID(raster) AS
srid FROM {view}"
+ )
+ .to_arrow_table()
+ .to_pylist()[0]
+ )
+ georef = [float(v) for v in meta["georef"].split()]
+ scale_x, skew_y, skew_x, scale_y, ul_x, ul_y = georef
+ srid = meta["srid"]
+ # World center of pixel (0, 0): upper-left corner + half a pixel step.
+ cx = ul_x + 0.5 * scale_x + 0.5 * skew_x
+ cy = ul_y + 0.5 * skew_y + 0.5 * scale_y
+
+ values = (
+ con.sql(
+ f"""
+ SELECT RS_Values(
+ raster,
+ ST_GeomFromText('MULTIPOINT ({cx} {cy})', 'EPSG:{srid}'),
+ 1
+ ) AS v FROM {view}
+ """
+ )
+ .to_arrow_table()["v"]
+ .to_pylist()
+ )
+ assert values == [[2324.0]]
+ finally:
+ con.drop_view(view)
+
+
+def test_rs_values_matches_rasterio(con):
+ """Cross-check RS_Values against rasterio on a random raster.
+
+ The plural counterpart of `test_rs_value_matches_rasterio`: the same dense
+ set of sample points is passed as a single MultiPoint, so one `RS_Values`
+ call returns a list that must match rasterio's per-point reads element for
+ element, in order.
+ """
+ import numpy as np
+
+ 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)
Review Comment:
Something bbox-parameterized here would be helpful as well
##########
rust/sedona-geometry/src/transform.rs:
##########
@@ -415,6 +415,57 @@ where
Ok(())
}
+/// Visit each point of a Point/MultiPoint geometry as `Some((x, y))` — or
+/// `None` for an empty (sub-)point — optionally transforming each coordinate
+/// with `trans` before visiting.
+///
+/// Compared to [`transform`], which writes a transformed WKB copy, this hands
+/// the caller each coordinate directly, so point-sampling consumers can
+/// transform and consume in one pass with no intermediate geometry
+/// materialisation. Any other geometry type is an error.
+pub fn visit_point_coords<F>(
+ geom: impl GeometryTrait<T = f64>,
+ trans: Option<&dyn CrsTransform>,
+ mut visit: F,
+) -> Result<(), SedonaGeometryError>
+where
+ F: FnMut(Option<(f64, f64)>) -> Result<(), SedonaGeometryError>,
+{
+ fn visit_one<P, F>(
+ point: &P,
+ trans: Option<&dyn CrsTransform>,
+ visit: &mut F,
+ ) -> Result<(), SedonaGeometryError>
+ where
+ P: PointTrait<T = f64>,
+ F: FnMut(Option<(f64, f64)>) -> Result<(), SedonaGeometryError>,
+ {
+ match point.coord() {
+ Some(coord) => {
+ let mut xy = (coord.x(), coord.y());
+ if let Some(trans) = trans {
+ trans.transform_coord(&mut xy)?;
+ }
+ visit(Some(xy))
+ }
+ None => visit(None),
+ }
+ }
+
+ match geom.as_type() {
+ GeometryType::Point(point) => visit_one(point, trans, &mut visit),
+ GeometryType::MultiPoint(multi_point) => {
+ for point in multi_point.points() {
+ visit_one(&point, trans, &mut visit)?;
+ }
+ Ok(())
+ }
+ _ => Err(SedonaGeometryError::Invalid(
+ "expected a Point or MultiPoint geometry".to_string(),
+ )),
Review Comment:
I think you could reasonably also accept a GeometryCollection whose members
are points or multipoints?
##########
rust/sedona-raster-functions/src/crs_utils.rs:
##########
@@ -17,11 +17,32 @@
use std::borrow::Cow;
+use datafusion_common::config::ConfigOptions;
use datafusion_common::{exec_datafusion_err, exec_err, DataFusionError,
Result};
+use sedona_common::option::SedonaOptions;
use sedona_geometry::transform::{transform, CrsEngine};
+use sedona_proj::transform::with_global_proj_engine;
use sedona_schema::crs::{deserialize_crs, CoordinateReferenceSystem, Crs,
CrsRef};
use wkb::reader::read_wkb;
+/// 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>(
Review Comment:
Same as the note on this in the RS_Clip PR...if the issue is the tester, we
should fix the tester so it can be invoked and sees a reasonable CRS engine
where it needs to be. I'd like to drop the sedona-proj dependency completely
here and this workaround doesn't quite get us closer to that.
--
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]
