james-willis commented on code in PR #1027:
URL: https://github.com/apache/sedona-db/pull/1027#discussion_r3606256847
##########
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:
removed. wasn't needed
##########
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:
done
--
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]
