From f14c66eef9d5fa664750c3d806682f6c7e4bf050 Mon Sep 17 00:00:00 2001 From: jameswillis Date: Wed, 2 Sep 2026 13:19:11 -0700 Subject: [PATCH 1/2] feat(rust/sedona-raster-functions): add RS_GeoTransform composite accessor Returns the geotransform decomposition Sedona Spark reports: a struct of magnitudeI, magnitudeJ, thetaI, thetaIJ, offsetX, offsetY. The math mirrors Spark's RasterAccessors#getGeoTransform exactly (including the acos sign tests) since Spark is the parity target; thetaI therefore only agrees with the existing RS_Rotation when |skewX| == |skewY|. The parity harness's result_to_tuples now passes struct columns through as dicts (they cannot be cast to string), matching the list-column treatment, and the spark-parity suite gains anchored north-up and 3-4-5-skew cases. --- docs/reference/sql/rs_geotransform.qmd | 37 ++++ integration/spark-parity/test_rs_scalar.py | 55 +++++ python/sedonadb/python/sedonadb/testing.py | 6 +- .../tests/functions/test_raster_functions.py | 30 +++ rust/sedona-raster-functions/src/register.rs | 1 + .../src/rs_geotransform.rs | 194 +++++++++++++++++- 6 files changed, 321 insertions(+), 2 deletions(-) create mode 100644 docs/reference/sql/rs_geotransform.qmd diff --git a/docs/reference/sql/rs_geotransform.qmd b/docs/reference/sql/rs_geotransform.qmd new file mode 100644 index 0000000000..a4d2c904b9 --- /dev/null +++ b/docs/reference/sql/rs_geotransform.qmd @@ -0,0 +1,37 @@ +# 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_GeoTransform +description: > + Returns the raster's geotransform decomposed into a struct of six doubles, + matching Sedona Spark's RS_GeoTransform: magnitudeI and magnitudeJ (pixel + sizes along the transformed i and j axes), thetaI (rotation of the raster in + radians, positive clockwise), thetaIJ (angle from the transformed i axis to + the transformed j axis in radians, positive counter-clockwise), and offsetX + and offsetY (the upper-left corner of the upper-left pixel). For the raw + affine coefficients instead of this decomposition, use the individual + accessors (RS_ScaleX, RS_SkewX, RS_UpperLeftX, ...) or RS_GeoReference. +kernels: + - returns: struct + args: [raster] +--- + +## Examples + +```sql +SELECT RS_GeoTransform(RS_Example()); +``` diff --git a/integration/spark-parity/test_rs_scalar.py b/integration/spark-parity/test_rs_scalar.py index 4018cc8c5a..c60f8cfad8 100644 --- a/integration/spark-parity/test_rs_scalar.py +++ b/integration/spark-parity/test_rs_scalar.py @@ -33,6 +33,8 @@ engine raises, that error is what trips the xfail. """ +import math + import pytest from sedonadb.testing import SedonaDB, compare @@ -156,6 +158,59 @@ def test_rs_band_nodata_null_raster(tmp_path): compare(sql, sedona, spark) +def test_rs_geotransform_north_up(tmp_path): + """Both engines decompose a north-up geotransform into the same struct. + + The default fixture bbox (100, 482, 114, 500) over a 7x6 grid gives + scaleX=2, scaleY=-3 with no skew, so the anchor is hand-derivable: + magnitudes are the pixel sizes, thetaI is acos(1) = 0, and thetaIJ is + acos(0) = pi/2 negated by its sign test (the i-to-j separation of a + y-down raster is -90 degrees).""" + sedona, spark = SedonaDB(), SedonaSpark() + for eng in (sedona, spark): + eng.create_random_raster_view("gt_raster", tmp_path / "gt.tif") + sql = "SELECT RS_GeoTransform(rast) FROM gt_raster" + anchor = { + "magnitudeI": 2.0, + "magnitudeJ": 3.0, + "thetaI": 0.0, + "thetaIJ": -math.pi / 2, + "offsetX": 100.0, + "offsetY": 500.0, + } + compare(sql, sedona, spark, expected=[(anchor,)]) + + +def test_rs_geotransform_skewed(tmp_path): + """A sheared transform exercises the acos sign tests in the decomposition. + + skewX=5 and skewY=3 are deliberately distinct — equal skews are the + degenerate regime where the magnitudes coincide and thetaIJ collapses to + +/-pi/2. The 3-4-5 / 5-12-13 pairs keep the anchor exactly representable: + magnitudeI = sqrt(4^2 + 3^2) = 5, magnitudeJ = sqrt(12^2 + 5^2) = 13, + thetaI = -acos(4/5) (negative because skewY > 0), and thetaIJ = + -acos(-16/65) (products and magnitudes are exact, and the sign test + acos(-63/65) exceeds pi/2). The anchor repeats the implementation's + operation order so every value is bit-exact.""" + sedona, spark = SedonaDB(), SedonaSpark() + for eng in (sedona, spark): + eng.create_random_raster_view( + "gt_skew_raster", + tmp_path / "gt_skew.tif", + gdal_transform=(100.0, 4.0, 5.0, 500.0, 3.0, -12.0), + ) + sql = "SELECT RS_GeoTransform(rast) FROM gt_skew_raster" + anchor = { + "magnitudeI": 5.0, + "magnitudeJ": 13.0, + "thetaI": -math.acos(4.0 / 5.0), + "thetaIJ": -math.acos(-16.0 / 65.0), + "offsetX": 100.0, + "offsetY": 500.0, + } + compare(sql, sedona, spark, expected=[(anchor,)]) + + @pytest.mark.xfail( reason="SedonaDB coalesces a NULL band index to band 1 (unwrap_or(1) in " "rs_band_accessors.rs); Sedona Spark returns NULL" diff --git a/python/sedonadb/python/sedonadb/testing.py b/python/sedonadb/python/sedonadb/testing.py index 15ae3634a8..4c23449689 100644 --- a/python/sedonadb/python/sedonadb/testing.py +++ b/python/sedonadb/python/sedonadb/testing.py @@ -292,7 +292,9 @@ def result_to_tuples( Geometry columns are rendered as WKT strings. List columns (e.g. the `List` returned by `RS_Values`) can't be cast to string, so they pass through as Python lists and are compared by value — assert them with - an expected cell that is itself a list, e.g. ``[([1.0, None],)]``. + an expected cell that is itself a list, e.g. ``[([1.0, None],)]``. Struct + columns (e.g. the one returned by `RS_GeoTransform`) likewise pass + through, as field-name-to-value dicts. """ tab = self.result_to_table(result) columns = [] @@ -302,6 +304,8 @@ def result_to_tuples( columns.append(ga.format_wkt(col, precision=wkt_precision).to_pylist()) elif pa.types.is_list(col.type) or pa.types.is_large_list(col.type): columns.append(col.to_pylist()) + elif pa.types.is_struct(col.type): + columns.append(col.to_pylist()) else: columns.append(col.cast(pa.string()).to_pylist()) diff --git a/python/sedonadb/tests/functions/test_raster_functions.py b/python/sedonadb/tests/functions/test_raster_functions.py index 83b7a8f255..0e8f026c19 100644 --- a/python/sedonadb/tests/functions/test_raster_functions.py +++ b/python/sedonadb/tests/functions/test_raster_functions.py @@ -15,6 +15,8 @@ # specific language governing permissions and limitations # under the License. +import math + import numpy as np import pandas as pd import pytest @@ -241,6 +243,34 @@ def test_rs_value_matches_rasterio(con): assert got == expected +def test_rs_geotransform(): + # RS_Example's geotransform is scaleX=2, skewX=1, skewY=1, scaleY=2 with + # upper-left (43.08, 79.07). The expected struct follows Sedona Spark's + # RS_GeoTransform decomposition: magnitudes sqrt(scaleX^2 + skewY^2) and + # sqrt(scaleY^2 + skewX^2), thetaI = -acos(scaleX / magnitudeI) (negative + # because skewY > 0), and thetaIJ positive because its sign test stays + # under pi/2 for this south-up transform. The thetaIJ expression repeats + # the implementation's operation order: sqrt(5)**2 is 5.000000000000001, + # so this is acos(4 / that), not acos(0.8). + eng = SedonaDB() + magnitude = math.sqrt(5.0) + eng.assert_query_result( + "SELECT RS_GeoTransform(RS_Example())", + [ + ( + { + "magnitudeI": magnitude, + "magnitudeJ": magnitude, + "thetaI": -math.acos(2.0 / magnitude), + "thetaIJ": math.acos(4.0 / (magnitude * magnitude)), + "offsetX": 43.08, + "offsetY": 79.07, + }, + ) + ], + ) + + def test_rs_setgeoreference_roundtrips_with_getter(): # RS_GeoReference emits scaleX, skewY, skewX, scaleY, upperLeftX, upperLeftY; # RS_SetGeoReference accepts the same six values back (GDAL order). diff --git a/rust/sedona-raster-functions/src/register.rs b/rust/sedona-raster-functions/src/register.rs index 09651d6ace..e0305e746b 100644 --- a/rust/sedona-raster-functions/src/register.rs +++ b/rust/sedona-raster-functions/src/register.rs @@ -52,6 +52,7 @@ pub fn default_function_set() -> FunctionSet { crate::rs_example::rs_example_udf, crate::rs_georeference::rs_georeference_udf, crate::rs_isempty::rs_isempty_udf, + crate::rs_geotransform::rs_geotransform_udf, crate::rs_geotransform::rs_rotation_udf, crate::rs_geotransform::rs_scalex_udf, crate::rs_geotransform::rs_scaley_udf, diff --git a/rust/sedona-raster-functions/src/rs_geotransform.rs b/rust/sedona-raster-functions/src/rs_geotransform.rs index c13ad74e4d..eebe635b1a 100644 --- a/rust/sedona-raster-functions/src/rs_geotransform.rs +++ b/rust/sedona-raster-functions/src/rs_geotransform.rs @@ -14,11 +14,14 @@ // KIND, either express or implied. See the License for the // specific language governing permissions and limitations // under the License. +use std::f64::consts::FRAC_PI_2; use std::{sync::Arc, vec}; use crate::executor::RasterExecutor; use arrow_array::builder::Float64Builder; -use arrow_schema::DataType; +use arrow_array::{ArrayRef, StructArray}; +use arrow_buffer::NullBufferBuilder; +use arrow_schema::{DataType, Field, Fields}; use datafusion_common::error::Result; use datafusion_expr::{ColumnarValue, Volatility}; use sedona_expr::scalar_udf::{SedonaScalarKernel, SedonaScalarUDF}; @@ -124,6 +127,66 @@ pub fn rs_rotation_udf() -> SedonaScalarUDF { ) } +/// RS_GeoTransform() scalar UDF implementation +/// +/// Returns the raster's geotransform decomposed into pixel magnitudes along +/// the transformed i/j axes, rotation and axis-separation angles, and the +/// upper-left offsets, as a struct matching Sedona Spark's `RS_GeoTransform`. +pub fn rs_geotransform_udf() -> SedonaScalarUDF { + SedonaScalarUDF::new( + "rs_geotransform", + vec![Arc::new(RsGeoTransformComposite {})], + Volatility::Immutable, + ) +} + +fn geotransform_fields() -> Fields { + Fields::from(vec![ + Field::new("magnitudeI", DataType::Float64, false), + Field::new("magnitudeJ", DataType::Float64, false), + Field::new("thetaI", DataType::Float64, false), + Field::new("thetaIJ", DataType::Float64, false), + Field::new("offsetX", DataType::Float64, false), + Field::new("offsetY", DataType::Float64, false), + ]) +} + +/// Decompose a GDAL-ordered geotransform into the six values Sedona Spark's +/// `RS_GeoTransform` reports. The math (including the sign conventions of the +/// two `acos` sign tests) mirrors Sedona Spark's +/// `RasterAccessors#getGeoTransform`, which is the parity target; note that +/// `thetaI` therefore only agrees with `RS_Rotation` when `|skewX| == |skewY|`. +fn decompose_geotransform(gt: &[f64]) -> [f64; 6] { + let (offset_x, scale_x, skew_x) = (gt[0], gt[1], gt[2]); + let (offset_y, skew_y, scale_y) = (gt[3], gt[4], gt[5]); + + // Pixel sizes along the transformed i (west-east) and j (north-south) axes + let magnitude_i = (scale_x * scale_x + skew_y * skew_y).sqrt(); + let magnitude_j = (scale_y * scale_y + skew_x * skew_x).sqrt(); + + // Rotation of the raster (radians, positive clockwise) + let mut theta_i = (scale_x / magnitude_i).acos(); + if (skew_y / magnitude_i).acos() < FRAC_PI_2 { + theta_i = -theta_i; + } + + // Angle from the transformed i axis to the transformed j axis (radians, + // positive counter-clockwise) + let mut theta_ij = ((scale_x * skew_x + skew_y * scale_y) / (magnitude_i * magnitude_j)).acos(); + if ((-skew_y * skew_x + scale_x * scale_y) / (magnitude_i * magnitude_j)).acos() > FRAC_PI_2 { + theta_ij = -theta_ij; + } + + [ + magnitude_i, + magnitude_j, + theta_i, + theta_ij, + offset_x, + offset_y, + ] +} + #[derive(Debug, Clone)] enum GeoTransformParam { Rotation, @@ -181,13 +244,71 @@ impl SedonaScalarKernel for RsGeoTransform { } } +#[derive(Debug)] +struct RsGeoTransformComposite {} + +impl SedonaScalarKernel for RsGeoTransformComposite { + fn return_type(&self, args: &[SedonaType]) -> Result> { + let matcher = ArgMatcher::new( + vec![ArgMatcher::is_raster()], + SedonaType::Arrow(DataType::Struct(geotransform_fields())), + ); + + matcher.match_args(args) + } + + fn invoke_batch( + &self, + arg_types: &[SedonaType], + args: &[ColumnarValue], + ) -> Result { + let executor = RasterExecutor::new(arg_types, args); + let num_iterations = executor.num_iterations(); + let mut builders: Vec = (0..6) + .map(|_| Float64Builder::with_capacity(num_iterations)) + .collect(); + let mut validity = NullBufferBuilder::new(num_iterations); + + executor.execute_raster_void(|_i, raster_opt| { + match raster_opt { + None => { + validity.append_null(); + // The fields are non-nullable, so null rows carry a + // placeholder in every child under a null struct slot. + for builder in builders.iter_mut() { + builder.append_value(0.0); + } + } + Some(raster) => { + validity.append_non_null(); + let components = decompose_geotransform(raster.transform()); + for (builder, component) in builders.iter_mut().zip(components) { + builder.append_value(component); + } + } + } + Ok(()) + })?; + + let arrays: Vec = builders + .iter_mut() + .map(|builder| Arc::new(builder.finish()) as ArrayRef) + .collect(); + let struct_array = StructArray::try_new(geotransform_fields(), arrays, validity.finish())?; + executor.finish(Arc::new(struct_array)) + } +} + #[cfg(test)] mod tests { use super::*; use arrow_array::Float64Array; + use arrow_buffer::NullBuffer; use datafusion_expr::ScalarUDF; use rstest::rstest; + use sedona_raster::builder::RasterBuilder; use sedona_schema::datatypes::RASTER; + use sedona_schema::raster::BandDataType; use sedona_testing::compare::assert_array_equal; use sedona_testing::rasters::generate_test_rasters; use sedona_testing::testers::ScalarUdfTester; @@ -197,6 +318,9 @@ mod tests { let udf: ScalarUDF = rs_rotation_udf().into(); assert_eq!(udf.name(), "rs_rotation"); + let udf: ScalarUDF = rs_geotransform_udf().into(); + assert_eq!(udf.name(), "rs_geotransform"); + let udf: ScalarUDF = rs_scalex_udf().into(); assert_eq!(udf.name(), "rs_scalex"); @@ -256,4 +380,72 @@ mod tests { let result = tester.invoke_array(Arc::new(rasters)).unwrap(); assert_array_equal(&result, &expected); } + + /// A north-up raster, a null, and a sheared raster with distinct skews + /// (equal skews are the degenerate regime where the magnitudes coincide + /// and thetaIJ collapses to +/-pi/2). The transforms are chosen so every + /// expected component below is exactly representable: the magnitudes are + /// square roots of perfect squares (3-4-5 and 5-12-13) and the `acos` + /// inputs are exact one-rounding quotients. + fn build_composite_test_rasters() -> arrow_array::StructArray { + let mut builder = RasterBuilder::new(3); + for transform in [ + Some([1.0, 3.0, 0.0, 2.0, 0.0, -4.0]), + None, + Some([10.0, 4.0, 5.0, 20.0, 3.0, -12.0]), + ] { + match transform { + None => builder.append_null().unwrap(), + Some([ulx, scale_x, skew_x, uly, skew_y, scale_y]) => { + builder + .start_raster_2d(2, 2, ulx, uly, scale_x, scale_y, skew_x, skew_y, None) + .unwrap(); + builder.start_band_2d(BandDataType::UInt8, None).unwrap(); + builder.band_data_writer().append_value([0u8; 4]); + builder.finish_band().unwrap(); + builder.finish_raster().unwrap(); + } + } + } + builder.finish().unwrap() + } + + #[test] + fn udf_invoke_composite() { + let tester = ScalarUdfTester::new(rs_geotransform_udf().into(), vec![RASTER]); + + // North-up: magnitudes are |scaleX|/|scaleY|, thetaI is acos(1) = 0, + // and thetaIJ is acos(0) = pi/2 negated by its sign test (the + // i-to-j separation of a y-down raster is -90 degrees). + let north_up = [3.0, 4.0, 0.0, -FRAC_PI_2, 1.0, 2.0]; + // Sheared with distinct skews: magnitudeI = sqrt(4^2 + 3^2) = 5, + // magnitudeJ = sqrt(12^2 + 5^2) = 13, thetaI = acos(4/5) negated + // because skewY > 0, and thetaIJ = acos(-16/65) negated because its + // sign test acos(-63/65) exceeds pi/2. + let skewed = [ + 5.0, + 13.0, + -(4.0f64 / 5.0).acos(), + -(-16.0f64 / 65.0).acos(), + 10.0, + 20.0, + ]; + + let columns: Vec = (0..6) + .map(|i| Arc::new(Float64Array::from(vec![north_up[i], 0.0, skewed[i]])) as ArrayRef) + .collect(); + let expected: Arc = Arc::new( + StructArray::try_new( + geotransform_fields(), + columns, + Some(NullBuffer::from(vec![true, false, true])), + ) + .unwrap(), + ); + + let result = tester + .invoke_array(Arc::new(build_composite_test_rasters())) + .unwrap(); + assert_array_equal(&result, &expected); + } } From f94f2a149a56f4835ef480310ae5a0bb39aaa89a Mon Sep 17 00:00:00 2001 From: jameswillis Date: Wed, 2 Sep 2026 13:38:03 -0700 Subject: [PATCH 2/2] test(integration): one spark-parity file per RS_ function Move the RS_GeoTransform cases into their own test_rs_geotransform.py, rename test_rs_scalar.py -> test_rs_bandnodatavalue.py and test_rs_raster_out.py -> test_rs_setbandnodatavalue.py to match their single-function contents, and record the file-per-function convention in the suite README. --- integration/spark-parity/README.md | 3 + ...s_scalar.py => test_rs_bandnodatavalue.py} | 60 +------------- .../spark-parity/test_rs_geotransform.py | 82 +++++++++++++++++++ ...r_out.py => test_rs_setbandnodatavalue.py} | 15 ++-- 4 files changed, 96 insertions(+), 64 deletions(-) rename integration/spark-parity/{test_rs_scalar.py => test_rs_bandnodatavalue.py} (74%) create mode 100644 integration/spark-parity/test_rs_geotransform.py rename integration/spark-parity/{test_rs_raster_out.py => test_rs_setbandnodatavalue.py} (95%) diff --git a/integration/spark-parity/README.md b/integration/spark-parity/README.md index aada3d3e7d..a75de481f6 100644 --- a/integration/spark-parity/README.md +++ b/integration/spark-parity/README.md @@ -65,6 +65,9 @@ Useful knobs, both read by `sedonadb.testing_spark`: ## Conventions +One test file per RS_ function, named `test_rs_.py` — parity coverage +for a function lands in its own file, not appended to a neighbor's. + Where the two engines are known to diverge and we intend to close the gap, mark the case `xfail(reason=...)` rather than deleting or loosening it. The suite then doubles as a catalog of what to fix, and flips to `xpass` the day the fix lands. diff --git a/integration/spark-parity/test_rs_scalar.py b/integration/spark-parity/test_rs_bandnodatavalue.py similarity index 74% rename from integration/spark-parity/test_rs_scalar.py rename to integration/spark-parity/test_rs_bandnodatavalue.py index c60f8cfad8..5b5bfeee08 100644 --- a/integration/spark-parity/test_rs_scalar.py +++ b/integration/spark-parity/test_rs_bandnodatavalue.py @@ -14,7 +14,7 @@ # KIND, either express or implied. See the License for the # specific language governing permissions and limitations # under the License. -"""SedonaDB vs Sedona Spark parity for the scalar raster readers. +"""SedonaDB vs Sedona Spark parity for `RS_BandNoDataValue`. Sedona Spark is the compatibility target, so each test runs one shared SQL string on both engines and asserts they agree with the harness-level @@ -33,8 +33,6 @@ engine raises, that error is what trips the xfail. """ -import math - import pytest from sedonadb.testing import SedonaDB, compare @@ -136,7 +134,8 @@ def test_rs_band_nodata_fractional_on_int_band(tmp_path): ) def test_rs_band_nodata_out_of_range_band(band, tmp_path): """An out-of-range band index gets the same answer from both engines. - Contrast the setter, which both engines refuse (see test_rs_raster_out.py).""" + Contrast the setter, which both engines refuse (see + test_rs_setbandnodatavalue.py).""" sedona, spark = SedonaDB(), SedonaSpark() for eng in (sedona, spark): eng.create_random_raster_view("oob_raster", tmp_path / "oob.tif", nodata=7.0) @@ -158,59 +157,6 @@ def test_rs_band_nodata_null_raster(tmp_path): compare(sql, sedona, spark) -def test_rs_geotransform_north_up(tmp_path): - """Both engines decompose a north-up geotransform into the same struct. - - The default fixture bbox (100, 482, 114, 500) over a 7x6 grid gives - scaleX=2, scaleY=-3 with no skew, so the anchor is hand-derivable: - magnitudes are the pixel sizes, thetaI is acos(1) = 0, and thetaIJ is - acos(0) = pi/2 negated by its sign test (the i-to-j separation of a - y-down raster is -90 degrees).""" - sedona, spark = SedonaDB(), SedonaSpark() - for eng in (sedona, spark): - eng.create_random_raster_view("gt_raster", tmp_path / "gt.tif") - sql = "SELECT RS_GeoTransform(rast) FROM gt_raster" - anchor = { - "magnitudeI": 2.0, - "magnitudeJ": 3.0, - "thetaI": 0.0, - "thetaIJ": -math.pi / 2, - "offsetX": 100.0, - "offsetY": 500.0, - } - compare(sql, sedona, spark, expected=[(anchor,)]) - - -def test_rs_geotransform_skewed(tmp_path): - """A sheared transform exercises the acos sign tests in the decomposition. - - skewX=5 and skewY=3 are deliberately distinct — equal skews are the - degenerate regime where the magnitudes coincide and thetaIJ collapses to - +/-pi/2. The 3-4-5 / 5-12-13 pairs keep the anchor exactly representable: - magnitudeI = sqrt(4^2 + 3^2) = 5, magnitudeJ = sqrt(12^2 + 5^2) = 13, - thetaI = -acos(4/5) (negative because skewY > 0), and thetaIJ = - -acos(-16/65) (products and magnitudes are exact, and the sign test - acos(-63/65) exceeds pi/2). The anchor repeats the implementation's - operation order so every value is bit-exact.""" - sedona, spark = SedonaDB(), SedonaSpark() - for eng in (sedona, spark): - eng.create_random_raster_view( - "gt_skew_raster", - tmp_path / "gt_skew.tif", - gdal_transform=(100.0, 4.0, 5.0, 500.0, 3.0, -12.0), - ) - sql = "SELECT RS_GeoTransform(rast) FROM gt_skew_raster" - anchor = { - "magnitudeI": 5.0, - "magnitudeJ": 13.0, - "thetaI": -math.acos(4.0 / 5.0), - "thetaIJ": -math.acos(-16.0 / 65.0), - "offsetX": 100.0, - "offsetY": 500.0, - } - compare(sql, sedona, spark, expected=[(anchor,)]) - - @pytest.mark.xfail( reason="SedonaDB coalesces a NULL band index to band 1 (unwrap_or(1) in " "rs_band_accessors.rs); Sedona Spark returns NULL" diff --git a/integration/spark-parity/test_rs_geotransform.py b/integration/spark-parity/test_rs_geotransform.py new file mode 100644 index 0000000000..ab81691734 --- /dev/null +++ b/integration/spark-parity/test_rs_geotransform.py @@ -0,0 +1,82 @@ +# 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. +"""SedonaDB vs Sedona Spark parity for `RS_GeoTransform`. + +Each test runs one shared SQL string on both engines and asserts they agree +via the harness-level `compare(sql, sedona, spark, expected=...)`; see the +README for the suite's conventions. The anchors here repeat the +implementation's floating-point operation order (both engines compute e.g. +sqrt(a*a + b*b) then divide before acos), so every expected value is +bit-exact rather than approximate. +""" + +import math + +from sedonadb.testing import SedonaDB, compare +from sedonadb.testing_spark import SedonaSpark + + +def test_rs_geotransform_north_up(tmp_path): + """Both engines decompose a north-up geotransform into the same struct. + + The default fixture bbox (100, 482, 114, 500) over a 7x6 grid gives + scaleX=2, scaleY=-3 with no skew, so the anchor is hand-derivable: + magnitudes are the pixel sizes, thetaI is acos(1) = 0, and thetaIJ is + acos(0) = pi/2 negated by its sign test (the i-to-j separation of a + y-down raster is -90 degrees).""" + sedona, spark = SedonaDB(), SedonaSpark() + for eng in (sedona, spark): + eng.create_random_raster_view("gt_raster", tmp_path / "gt.tif") + sql = "SELECT RS_GeoTransform(rast) FROM gt_raster" + anchor = { + "magnitudeI": 2.0, + "magnitudeJ": 3.0, + "thetaI": 0.0, + "thetaIJ": -math.pi / 2, + "offsetX": 100.0, + "offsetY": 500.0, + } + compare(sql, sedona, spark, expected=[(anchor,)]) + + +def test_rs_geotransform_skewed(tmp_path): + """A sheared transform exercises the acos sign tests in the decomposition. + + skewX=5 and skewY=3 are deliberately distinct — equal skews are the + degenerate regime where the magnitudes coincide and thetaIJ collapses to + +/-pi/2. The 3-4-5 / 5-12-13 pairs keep the anchor exactly representable: + magnitudeI = sqrt(4^2 + 3^2) = 5, magnitudeJ = sqrt(12^2 + 5^2) = 13, + thetaI = -acos(4/5) (negative because skewY > 0), and thetaIJ = + -acos(-16/65) (products and magnitudes are exact, and the sign test + acos(-63/65) exceeds pi/2).""" + sedona, spark = SedonaDB(), SedonaSpark() + for eng in (sedona, spark): + eng.create_random_raster_view( + "gt_skew_raster", + tmp_path / "gt_skew.tif", + gdal_transform=(100.0, 4.0, 5.0, 500.0, 3.0, -12.0), + ) + sql = "SELECT RS_GeoTransform(rast) FROM gt_skew_raster" + anchor = { + "magnitudeI": 5.0, + "magnitudeJ": 13.0, + "thetaI": -math.acos(4.0 / 5.0), + "thetaIJ": -math.acos(-16.0 / 65.0), + "offsetX": 100.0, + "offsetY": 500.0, + } + compare(sql, sedona, spark, expected=[(anchor,)]) diff --git a/integration/spark-parity/test_rs_raster_out.py b/integration/spark-parity/test_rs_setbandnodatavalue.py similarity index 95% rename from integration/spark-parity/test_rs_raster_out.py rename to integration/spark-parity/test_rs_setbandnodatavalue.py index b99ebf6030..ffa915850b 100644 --- a/integration/spark-parity/test_rs_raster_out.py +++ b/integration/spark-parity/test_rs_setbandnodatavalue.py @@ -14,7 +14,7 @@ # KIND, either express or implied. See the License for the # specific language governing permissions and limitations # under the License. -"""SedonaDB vs Sedona Spark parity for raster-in / raster-out functions. +"""SedonaDB vs Sedona Spark parity for `RS_SetBandNoDataValue`. These exercise the raster round-trip *out* of each engine: SedonaDB decodes its native raster column; Sedona Spark transports the result as GeoTIFF bytes @@ -22,17 +22,18 @@ `DecodedRaster` (pixels + geotransform + per-band nodata), which is how the harness-level `sedonadb.testing.compare` compares a raster result — a NULL raster from every engine also counts as agreement. Same -`xfail`-for-known-divergence policy as the scalar suite. +`xfail`-for-known-divergence policy as the rest of the suite. Both engines are constructed directly rather than through fixtures: this suite is only run deliberately, so a missing pyspark, JVM, or Sedona jar should be a failure with a real traceback, not a skip. `SedonaSpark` caches its `SparkSession` on the class, so building one per test reuses the same JVM. -`RS_SetBandNoDataValue` is the first case on purpose: it is raster-in/raster-out -but passes pixels through untouched, so a mismatch is a round-trip bug, not an -operation divergence — it isolates the transport. Pixel-transforming ops -(RS_Resample, RS_MapAlgebra) come next, with their known divergences as xfails. +`RS_SetBandNoDataValue` was the harness's first raster-out case on purpose: it +is raster-in/raster-out but passes pixels through untouched, so a mismatch is a +round-trip bug, not an operation divergence — it isolates the transport. +Pixel-transforming ops (RS_Resample, RS_MapAlgebra) each get their own file, +with their known divergences as xfails. """ import pytest @@ -207,7 +208,7 @@ def test_rs_setbandnodata_negative_on_uint8(tmp_path): @pytest.mark.parametrize("band", [0, 3]) def test_rs_setbandnodata_out_of_range_band_rejected(band, tmp_path): """Both engines refuse an out-of-range band index — unlike the getter, - where SedonaDB returns NULL (see test_rs_scalar.py).""" + where SedonaDB returns NULL (see test_rs_bandnodatavalue.py).""" sedona, spark = SedonaDB(), SedonaSpark() for eng in (sedona, spark): eng.create_random_raster_view(