You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@sedona.apache.org by ja...@apache.org on 2026/07/29 16:36:42 UTC
(sedona-db) branch main updated: feat(rust/sedona-raster-gdal): add RS_Resample UDF (#1064)
This is an automated email from the ASF dual-hosted git repository.
james-willis 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 1965063ae feat(rust/sedona-raster-gdal): add RS_Resample UDF (#1064)
1965063ae is described below
commit 1965063ae44599eaf5deafa1fc6e57177f7007aa
Author: James Willis <ja...@users.noreply.github.com>
AuthorDate: Wed Jul 29 09:36:37 2026 -0700
feat(rust/sedona-raster-gdal): add RS_Resample UDF (#1064)
---
c/sedona-gdal/src/gdal.rs | 9 +
docs/reference/sql/rs_resample.qmd | 145 +++
python/sedonadb/python/sedonadb/raster_testing.py | 227 +++-
.../sedonadb/tests/functions/test_rs_resample.py | 367 +++++-
rust/sedona-raster-gdal/Cargo.toml | 5 +
rust/sedona-raster-gdal/benches/rs_resample.rs | 130 ++
rust/sedona-raster-gdal/src/lib.rs | 2 +
rust/sedona-raster-gdal/src/register.rs | 1 +
rust/sedona-raster-gdal/src/rs_reproject_match.rs | 113 +-
rust/sedona-raster-gdal/src/rs_resample.rs | 1344 ++++++++++++++++++++
rust/sedona-raster-gdal/src/utils.rs | 247 +++-
11 files changed, 2438 insertions(+), 152 deletions(-)
diff --git a/c/sedona-gdal/src/gdal.rs b/c/sedona-gdal/src/gdal.rs
index 58f94aac0..637b9ad8a 100644
--- a/c/sedona-gdal/src/gdal.rs
+++ b/c/sedona-gdal/src/gdal.rs
@@ -72,6 +72,15 @@ impl Gdal {
self.api.version_info(request)
}
+ /// The linked GDAL library version as an integer, from
+ /// `GDALVersionInfo("VERSION_NUM")`: `major*1_000_000 + minor*10_000 +
+ /// rev*100 + build` (for example `3_080_400` for GDAL 3.8.4, `3_130_000` for
+ /// 3.13.0). Falls back to `0` when the value cannot be parsed, which callers
+ /// treat as "oldest" — the most conservative choice for version gates.
+ pub fn version_num(&self) -> i32 {
+ self.version_info("VERSION_NUM").trim().parse().unwrap_or(0)
+ }
+
// -- Config --------------------------------------------------------------
/// Set a thread-local GDAL configuration option.
diff --git a/docs/reference/sql/rs_resample.qmd b/docs/reference/sql/rs_resample.qmd
new file mode 100644
index 000000000..2efc40510
--- /dev/null
+++ b/docs/reference/sql/rs_resample.qmd
@@ -0,0 +1,145 @@
+---
+# 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_Resample
+description: >
+ Resamples a raster onto a new pixel grid — a new width/height, a new pixel
+ size, and/or a snapped origin — in the same CRS.
+kernels:
+ - returns: raster
+ args:
+ - name: raster
+ type: raster
+ description: The input raster to resample.
+ - name: reference_raster
+ type: raster
+ description: >
+ A raster whose grid — dimensions (or pixel size, when `use_scale` is
+ true) and origin — defines the output grid. Must share the input's CRS.
+ - name: use_scale
+ type: boolean
+ description: >
+ When false, `width_or_scale`/`height_or_scale` are output pixel
+ dimensions and the world extent is preserved; when true, they are the
+ target pixel size.
+ - name: algorithm
+ type: string
+ description: >
+ Resampling algorithm (case-insensitive): `NearestNeighbor` (the
+ default), `Bilinear`, `Cubic` (alias `Bicubic`), `CubicSpline`,
+ `Lanczos`, `Average`, or `Mode`.
+ - returns: raster
+ args:
+ - name: raster
+ type: raster
+ - name: width_or_scale
+ type: double
+ description: >
+ Output width in pixels when `use_scale` is false (a fractional value is
+ an error), or the target x pixel size (`scaleX`) when `use_scale` is true.
+ - name: height_or_scale
+ type: double
+ description: >
+ Output height in pixels when `use_scale` is false (a fractional value is
+ an error), or the target y pixel size (`scaleY`, negative for a north-up
+ raster) when `use_scale` is true.
+ - name: use_scale
+ type: boolean
+ - name: algorithm
+ type: string
+ - returns: raster
+ args:
+ - name: raster
+ type: raster
+ - name: width_or_scale
+ type: double
+ - name: height_or_scale
+ type: double
+ - name: grid_x
+ type: double
+ description: X coordinate of the grid the output origin snaps to.
+ - name: grid_y
+ type: double
+ description: Y coordinate of the grid the output origin snaps to.
+ - name: use_scale
+ type: boolean
+ - name: algorithm
+ type: string
+---
+
+::: callout-warning
+**Experimental.** This function is experimental; its behavior may change without
+notice.
+:::
+
+## Description
+
+`RS_Resample` changes a raster's pixel grid while keeping its band count and
+order and per-band nodata. Pixel values are recomputed by GDAL with the chosen
+`algorithm` (nearest neighbour by default). The output stays in the input's
+CRS — `RS_Resample` never reprojects (use
+[`RS_ReprojectMatch`](rs_reprojectmatch.qmd) for that).
+
+`RS_Resample` has three positional overloads:
+
+- **`RS_Resample(raster, width_or_scale, height_or_scale, use_scale, algorithm)`** —
+ when `use_scale` is false, `width_or_scale`/`height_or_scale` are the output pixel
+ dimensions (a fractional value is an error) and the world extent is preserved.
+ When `use_scale` is true they are the target pixel size (`scaleX`, `scaleY`,
+ with `scaleY` negative for a north-up raster); the pixel size is kept exact and
+ the extent grows by up to one pixel so it tiles into whole pixels. Cells
+ outside the source footprint become nodata.
+- **`RS_Resample(raster, width_or_scale, height_or_scale, grid_x, grid_y, use_scale,
+ algorithm)`** — additionally snaps the output origin to the grid anchored at
+ `(grid_x, grid_y)`, expanding outward so the original extent stays covered.
+- **`RS_Resample(raster, reference_raster, use_scale, algorithm)`** — takes the
+ target dimensions (or pixel size, when `use_scale`) and the origin grid from
+ `reference_raster`, which must share the input's CRS (it is an error otherwise).
+
+`algorithm` (case-insensitive) is one of `NearestNeighbor` (the default),
+`Bilinear`, `Cubic` (alias `Bicubic`), `CubicSpline`, `Lanczos`, `Average`, or
+`Mode`.
+
+Dimensions and pixel size use the rows-then-columns intuition of the raster
+(`height_or_scale`/`scaleY` for rows, `width_or_scale`/`scaleX` for columns). In the
+extent-preserving dimension form, a skewed or rotated raster is re-gridded with
+its footprint preserved (both the scale and skew terms of each axis scale by the
+same factor). A scale change or origin snap on a skewed raster is not supported
+and errors.
+
+For an N-D raster (bands with extra dimensions such as `time`), the resample is
+a 2-D `(y, x)` operation broadcast across every non-spatial plane, so the
+non-spatial dimensions are preserved and only the `(y, x)` extent changes.
+
+::: callout-note
+Extent-growing scale changes and origin snaps read the entire raster into memory.
+:::
+
+## Examples
+
+Downsample the example raster (64 x 32) to half its dimensions:
+
+```sql
+SELECT RS_Width(RS_Resample(RS_Example(), 32, 16, false, 'NearestNeighbor'));
+```
+
+Upsample to twice its dimensions with bilinear resampling:
+
+```sql
+SELECT RS_Height(RS_Resample(RS_Example(), 128, 64, false, 'Bilinear'));
+```
diff --git a/python/sedonadb/python/sedonadb/raster_testing.py b/python/sedonadb/python/sedonadb/raster_testing.py
index a1806d5dd..68825501d 100644
--- a/python/sedonadb/python/sedonadb/raster_testing.py
+++ b/python/sedonadb/python/sedonadb/raster_testing.py
@@ -266,12 +266,22 @@ class RasterEngine:
self._not_implemented("as_raster")
def resample(
- self, path, *, width: int, height: int, algorithm: str = "nearestneighbor"
+ self,
+ path,
+ *,
+ width: Optional[int] = None,
+ height: Optional[int] = None,
+ scale_x: Optional[float] = None,
+ scale_y: Optional[float] = None,
+ algorithm: str = "NearestNeighbor",
) -> DecodedRaster:
- """Resample the raster to `width` x `height` pixels over the same extent.
+ """Resample the GeoTIFF at `path` onto a new grid, in the same CRS.
- `algorithm` is one of nearestneighbor, bilinear, bicubic
- (case-insensitive; each engine maps to its own resampler names).
+ Selects the output grid two ways: target dimensions (`width`/`height`,
+ extent-preserving) or a target pixel size (`scale_x`/`scale_y`, which
+ keeps the size exact and grows the extent to whole pixels). Band count
+ and nodata are preserved; the result is a decoded raster (pixels,
+ geotransform, nodata). `RS_Resample` never reprojects.
"""
self._not_implemented("resample")
@@ -516,6 +526,81 @@ class SedonaDB(RasterEngine):
).to_arrow_table()["r"]
return decode_raster(result[0])
+ def resample(
+ self,
+ path,
+ *,
+ width=None,
+ height=None,
+ scale_x=None,
+ scale_y=None,
+ algorithm="NearestNeighbor",
+ ):
+ # useScale selects between the dimension and pixel-size forms of the
+ # 5-argument overload
+ # RS_Resample(raster, widthOrScale, heightOrScale, useScale, algorithm).
+ if scale_x is not None:
+ use_scale, width_or_scale, height_or_scale = (
+ True,
+ float(scale_x),
+ float(scale_y),
+ )
+ else:
+ use_scale, width_or_scale, height_or_scale = (
+ False,
+ float(width),
+ float(height),
+ )
+
+ # Arguments travel as table columns (not literals) so the kernel runs
+ # its real array path.
+ table = pa.table(
+ {
+ "path": pa.array([str(path)], type=pa.utf8()),
+ "width_or_scale": pa.array([width_or_scale], type=pa.float64()),
+ "height_or_scale": pa.array([height_or_scale], type=pa.float64()),
+ "use_scale": pa.array([use_scale], type=pa.bool_()),
+ "algorithm": pa.array([algorithm], type=pa.utf8()),
+ }
+ )
+ df = self._con.create_data_frame(table)
+ result = df.select(
+ r=df.path.funcs.rs_frompath().funcs.rs_resample(
+ df.width_or_scale,
+ df.height_or_scale,
+ df.use_scale,
+ df.algorithm,
+ ),
+ ).to_arrow_table()["r"]
+ return decode_raster(result[0])
+
+ def resample_to_reference(
+ self, path, reference_path, *, use_scale=False, algorithm="NearestNeighbor"
+ ) -> Optional[DecodedRaster]:
+ """Run the 4-argument reference overload
+ `RS_Resample(raster, referenceRaster, useScale, algorithm)`.
+
+ The target dimensions (or pixel size, when `use_scale`) and origin come
+ from `referenceRaster`, which must share the input's CRS.
+ """
+ table = pa.table(
+ {
+ "path": pa.array([str(path)], type=pa.utf8()),
+ "ref": pa.array([str(reference_path)], type=pa.utf8()),
+ "use_scale": pa.array([use_scale], type=pa.bool_()),
+ "algorithm": pa.array([algorithm], type=pa.utf8()),
+ }
+ )
+ df = self._con.create_data_frame(table)
+ result = df.select(
+ r=df.path.funcs.rs_frompath().funcs.rs_resample(
+ df.ref.funcs.rs_frompath(),
+ df.use_scale,
+ df.algorithm,
+ ),
+ ).to_arrow_table()["r"]
+ return decode_raster(result[0])
+
def reproject_match(self, path, reference_path, *, algorithm="NearestNeighbor"):
# Both rasters travel as table columns (not literals) so the kernel runs
# its real array path.
@@ -1066,39 +1151,87 @@ class Rasterio(RasterEngine):
)
return DecodedRaster(pixels[np.newaxis], tuple(transform.to_gdal()), [nodata])
- def resample(self, path, *, width, height, algorithm="nearestneighbor"):
+ def resample(
+ self,
+ path,
+ *,
+ width=None,
+ height=None,
+ scale_x=None,
+ scale_y=None,
+ algorithm="NearestNeighbor",
+ ):
import rasterio
- import rasterio.crs
- from rasterio.enums import Resampling
- from rasterio.warp import reproject
- resampling = {
- "nearestneighbor": Resampling.nearest,
- "bilinear": Resampling.bilinear,
- "bicubic": Resampling.cubic,
- }[algorithm.lower()]
+ resampling = _RASTERIO_RESAMPLING.get(algorithm.lower())
+ if resampling is None:
+ raise ValueError(f"unsupported resampling algorithm {algorithm!r}")
+
with rasterio.open(str(path)) as src:
- source = src.read()
- dst_transform = src.transform * src.transform.scale(
- src.width / width, src.height / height
- )
- # reproject() insists on a CRS; with the same one on both sides
- # nothing reprojects, so an arbitrary CRS keeps a CRS-less
- # fixture a pure grid resample.
- crs = src.crs or rasterio.crs.CRS.from_epsg(3857)
- destination = np.zeros((src.count, height, width), dtype=source.dtype)
- reproject(
- source,
- destination,
- src_transform=src.transform,
- src_crs=crs,
- dst_transform=dst_transform,
- dst_crs=crs,
- resampling=resampling,
+ if scale_x is not None:
+ return self._scale_grow(src, scale_x, scale_y, resampling)
+
+ # Dimension mode, same CRS: RS_Resample reads the source into the
+ # target grid (GDAL RasterIO), preserving the extent.
+ pixels = src.read(
+ out_shape=(src.count, height, width), resampling=resampling
)
- return DecodedRaster(
- destination, tuple(dst_transform.to_gdal()), list(src.nodatavals)
+ # Scale both axes of the source transform by the pixel-count ratio —
+ # the independent affine equivalent of RS_Resample's footprint-
+ # preserving grid (scale and skew terms scale together).
+ transform = src.transform * src.transform.scale(
+ src.width / width, src.height / height
)
+ nodata = [src.nodata] * src.count
+ return DecodedRaster(pixels, tuple(transform.to_gdal()), nodata)
+
+ def _scale_grow(self, src, scale_x, scale_y, resampling):
+ """Reference for scale mode: keep the pixel size exact and grow the extent.
+
+ Dimensions are `ceil(extent / pixel_size)` (so the last row/column may
+ extend past the source), the origin and skew are the source's, and the
+ grown border reads back as nodata. Same CRS on both sides — a pure
+ regrid, not a reprojection.
+ """
+ import math
+
+ from rasterio.transform import Affine
+
+ left, bottom, right, top = src.bounds
+ out_w = math.ceil(abs(right - left) / abs(scale_x))
+ out_h = math.ceil(abs(top - bottom) / abs(scale_y))
+ a = src.transform
+ # Affine(scale_x, skew_x, origin_x, skew_y, scale_y, origin_y): exact
+ # pixel size, source origin and skew.
+ dst_transform = Affine(scale_x, a.b, a.c, a.d, scale_y, a.f)
+ return self._warp_into(
+ src, dst_transform, out_w, out_h, src.crs, src.crs, resampling
+ )
+
+ @staticmethod
+ def _warp_into(src, dst_transform, dst_w, dst_h, src_crs, dst_crs, resampling):
+ import numpy as np
+ from rasterio.warp import reproject
+
+ src_data = src.read()
+ dst = np.zeros((src.count, dst_h, dst_w), dtype=src_data.dtype)
+ # `dst_nodata` fills cells outside the source footprint with the band
+ # nodata (0 when the band has none) — the destination pre-fill
+ # RS_Resample uses. Source nodata is deliberately not passed, matching
+ # RS_Resample's warp (source nodata values pass through, not masked).
+ reproject(
+ source=src_data,
+ destination=dst,
+ src_transform=src.transform,
+ src_crs=src_crs,
+ dst_transform=dst_transform,
+ dst_crs=dst_crs,
+ dst_nodata=src.nodata,
+ resampling=resampling,
+ )
+ return DecodedRaster(
+ dst, tuple(dst_transform.to_gdal()), [src.nodata] * src.count
+ )
def reproject_match(self, path, reference_path, *, algorithm="NearestNeighbor"):
"""Reference for `RS_ReprojectMatch`: warp the input onto the reference grid.
@@ -1219,6 +1352,36 @@ class Rasterio(RasterEngine):
return decode_geotiff_bytes(data)
+# RS_Resample algorithm names (lower-cased) to rasterio's Resampling enum, for
+# the reference engine. Only the algorithms the parity tests exercise are
+# mapped; unknown names raise in `Rasterio.resample`.
+def _rasterio_resampling_map():
+ from rasterio.enums import Resampling
+
+ return {
+ "nearestneighbor": Resampling.nearest,
+ "nearestneighbour": Resampling.nearest,
+ "bilinear": Resampling.bilinear,
+ "cubic": Resampling.cubic,
+ "average": Resampling.average,
+ }
+
+
+class _LazyResamplingMap:
+ """Builds the resampling map on first use so importing this module does not
+ require rasterio (it is an optional reference dependency)."""
+
+ _map = None
+
+ def get(self, key):
+ if _LazyResamplingMap._map is None:
+ _LazyResamplingMap._map = _rasterio_resampling_map()
+ return _LazyResamplingMap._map.get(key)
+
+
+_RASTERIO_RESAMPLING = _LazyResamplingMap()
+
+
def _is_nodata(sampled, nodata) -> bool:
"""Whether a sampled value equals the band nodata, NaN-aware (a NaN
sentinel matches NaN pixels, which bare `==` never would)."""
diff --git a/python/sedonadb/tests/functions/test_rs_resample.py b/python/sedonadb/tests/functions/test_rs_resample.py
index 13b2d88f9..82caa21bf 100644
--- a/python/sedonadb/tests/functions/test_rs_resample.py
+++ b/python/sedonadb/tests/functions/test_rs_resample.py
@@ -15,57 +15,352 @@
# specific language governing permissions and limitations
# under the License.
-"""RS_Resample parity against `rasterio.warp.reproject` on the same grid.
-
-Only nearest-neighbor resampling is compared: it is a pure pixel-pick, so
-engines that agree on grid alignment agree bit for bit, extremes included.
-The interpolating algorithms are not asserted against rasterio because the
-implementations differ at the raster border (geotools yields NaN outside the
-interpolation support where GDAL clamps to the edge), so there is no
-engine-independent expected value there.
+"""RS_Resample cross-checked against a rasterio reference implementation.
+
+SedonaDB is the subject; rasterio is the reference the results are checked
+against (no Sedona Spark here — Spark parity for RS_Resample is a later
+follow-up). Each test writes a raster once (numpy array + GDAL geotransform)
+and resamples it through both engines from `sedonadb.raster_testing`.
+`RS_Resample` is same-CRS only (it never reprojects), so:
+
+- **Dimension mode** (useScale=false) shares GDAL's RasterIO path with rasterio's
+ `Dataset.read(out_shape=..., resampling=...)`. Nearest-neighbour resampling is
+ a pure pixel-pick, so pixels are compared exactly (upsampling by an integer
+ factor is block replication; downsampling decimates on the same grid) —
+ including the planted dtype extremes; bilinear accumulates over neighbours and
+ rasterio may bundle a different GDAL build, so it uses a small tolerance.
+- **Scale mode** (useScale=true) keeps the pixel size exact and grows the extent
+ to whole pixels (Sedona Spark semantics), filling the grown border with nodata.
+ SedonaDB reads the covered sub-window with a fractional-window RasterIO; the
+ reference reproduces the same nearest sampling and nodata fill with
+ `rasterio.warp.reproject` (same CRS on both sides).
+- **The reference-raster overload** takes the target grid and origin from a
+ reference raster in the same CRS.
+
+Arguments travel as table columns (not literals) so the kernel runs its real
+array path. Error pathways are asserted directly against RS_Resample's messages.
"""
+import numpy as np
+import pyarrow as pa
import pytest
from sedonadb.raster_testing import (
+ Rasterio,
SedonaDB,
- assert_decoded_equal,
- random_raster_data,
- write_geotiff,
+ assert_transform_and_nodata,
+ decode_raster,
+ write_grid_geotiff,
+ write_random_geotiff,
)
-pytest.importorskip("rasterio")
-pytest.importorskip("shapely")
+# North-up GDAL-order geotransform: origin (100, 500), 2x2 pixels.
+GDAL_TRANSFORM = (100.0, 2.0, 0.0, 500.0, 0.0, -2.0)
+
+DTYPES = ["uint8", "uint16", "int16", "int32", "float32", "float64"]
-pytestmark = pytest.mark.skipif(
- not SedonaDB.implements("resample"),
- reason="RS_Resample is not implemented in SedonaDB (the parity subject)",
-)
-# GDAL-order geotransform: origin (100, 500), 2-wide by 3-tall north-up
-# pixels; with a 7x6 raster the extent is x in [100, 114], y in [482, 500].
-GDAL_TRANSFORM = (100.0, 2.0, 0.0, 500.0, 0.0, -3.0)
-HEIGHT, WIDTH = 6, 7
+@pytest.fixture()
+def sedona(con):
+ return SedonaDB(con)
+
+
+@pytest.fixture()
+def reference():
+ return Rasterio.create_or_skip()
+
+
+def _write(
+ tmp_path,
+ dtype,
+ *,
+ bands,
+ height,
+ width,
+ nodata=None,
+ transform=GDAL_TRANSFORM,
+ crs=None,
+ name="in",
+):
+ path = tmp_path / f"resample_{name}_{dtype}_{width}x{height}.tif"
+ write_random_geotiff(
+ path,
+ dtype,
+ bands=bands,
+ height=height,
+ width=width,
+ gdal_transform=transform,
+ crs=crs,
+ nodata=nodata,
+ )
+ return path
+
+
+@pytest.mark.parametrize("dtype", DTYPES)
+def test_rs_resample_nearest_upsample_matches_rasterio(
+ sedona, reference, tmp_path, dtype
+):
+ """A 2x integer nearest-neighbour upsample replicates each source pixel into
+ a 2x2 block. Because block replication is unambiguous, pixels are bit-exact
+ against rasterio; the transform halves both pixel dimensions while keeping
+ the extent, and nodata is preserved."""
+ tiff = _write(tmp_path, dtype, bands=2, height=3, width=4, nodata=None)
+
+ got = sedona.resample(tiff, width=8, height=6)
+ expected = reference.resample(tiff, width=8, height=6)
+
+ np.testing.assert_array_equal(got.pixels, expected.pixels)
+ assert got.pixels.dtype == expected.pixels.dtype
+ assert got.pixels.shape == (2, 6, 8)
+ assert_transform_and_nodata(got, expected)
@pytest.mark.parametrize("dtype", ["uint8", "float64"])
+def test_rs_resample_nearest_downsample_matches_rasterio(
+ sedona, reference, tmp_path, dtype
+):
+ """Nearest-neighbour downsampling to a non-integer pixel ratio still picks
+ source pixels verbatim (GDAL RasterIO decimation on both sides), so SedonaDB
+ and rasterio agree bit-for-bit, planted dtype extremes included. A 7x6 source
+ to 4x3 gives a 7/4 column ratio; the output transform is the unchanged extent
+ divided by the new shape."""
+ tiff = _write(tmp_path, dtype, bands=2, height=6, width=7, nodata=None)
+
+ got = sedona.resample(tiff, width=4, height=3)
+ expected = reference.resample(tiff, width=4, height=3)
+
+ np.testing.assert_array_equal(got.pixels, expected.pixels)
+ assert got.pixels.dtype == expected.pixels.dtype
+ assert got.pixels.shape == (2, 3, 4)
+ assert_transform_and_nodata(got, expected)
+
+
+@pytest.mark.parametrize("dtype", ["int64", "uint64"])
+def test_rs_resample_nearest_preserves_64bit_ints_exactly(sedona, tmp_path, dtype):
+ """A plain nearest width/height change preserves 64-bit integers bit-exactly:
+ the RasterIO decimation copies source samples in their native type, so values
+ above 2**53 (which GDAL's floating working type would round) survive. A 2x
+ upsample replicates each source pixel into a 2x2 block, compared against the
+ source read back at native resolution — all through SedonaDB, since rasterio
+ need not support 64-bit dtypes."""
+ pytest.importorskip("rasterio") # _write needs rasterio to author the input
+ tiff = _write(tmp_path, dtype, bands=1, height=3, width=4, nodata=None)
+
+ # An identity resample (same grid, nearest) reads the source pixels back.
+ src = sedona.resample(tiff, width=4, height=3).pixels
+ # random_raster_data plants the dtype extremes in opposite corners, so the
+ # source really does carry a value a double cannot represent exactly.
+ assert int(src.max()) > 2**53 or int(src.min()) < -(2**53)
+
+ got = sedona.resample(tiff, width=8, height=6)
+ expected = np.repeat(np.repeat(src, 2, axis=1), 2, axis=2)
+ np.testing.assert_array_equal(got.pixels, expected)
+ assert got.pixels.dtype == np.dtype(dtype)
+
+
+@pytest.mark.parametrize(
+ "kwargs",
+ [
+ {"width": 4, "height": 4, "algorithm": "Bilinear"}, # interpolating
+ {"scale_x": 4.0, "scale_y": -4.0}, # scale mode regrids through warp
+ ],
+)
+def test_rs_resample_rejects_64bit_ints_off_the_nearest_fast_path(
+ sedona, tmp_path, kwargs
+):
+ """Only a plain nearest width/height change preserves 64-bit integers. Any
+ interpolation or regrid routes pixels through GDAL's floating working type, so
+ RS_Resample rejects Int64/UInt64 there rather than silently corrupting them."""
+ pytest.importorskip("rasterio") # _write needs rasterio to author the input
+ tiff = _write(tmp_path, "int64", bands=1, height=8, width=8, nodata=None)
+ with pytest.raises(Exception, match="Int64/UInt64"):
+ sedona.resample(tiff, **kwargs)
+
+
+def test_rs_resample_nodata_preserved(sedona, reference, tmp_path):
+ """The per-band nodata survives a resample unchanged (metadata is copied,
+ not recomputed)."""
+ tiff = _write(tmp_path, "uint8", bands=1, height=3, width=4, nodata=200.0)
+
+ got = sedona.resample(tiff, width=8, height=6)
+ expected = reference.resample(tiff, width=8, height=6)
+
+ np.testing.assert_array_equal(got.pixels, expected.pixels)
+ assert got.pixels.dtype == expected.pixels.dtype
+ assert_transform_and_nodata(got, expected)
+ assert got.nodata == [200]
+
+
+def test_rs_resample_bilinear_matches_rasterio(sedona, reference, tmp_path):
+ """Bilinear resampling blends neighbouring pixels, so it is compared to
+ rasterio with a small tolerance (rasterio may use a different GDAL build).
+ Uses a float band so the interpolated values are not integer-truncated."""
+ tiff = _write(tmp_path, "float64", bands=1, height=8, width=8, nodata=None)
+
+ got = sedona.resample(tiff, width=4, height=4, algorithm="Bilinear")
+ expected = reference.resample(tiff, width=4, height=4, algorithm="Bilinear")
+
+ np.testing.assert_allclose(got.pixels, expected.pixels, rtol=1e-6, atol=1e-6)
+ assert got.pixels.shape == (1, 4, 4)
+ assert_transform_and_nodata(got, expected)
+
+
+def test_rs_resample_scale_mode_matches_dimension_mode(sedona, reference, tmp_path):
+ """Scale mode picks the output dimensions as ceil(extent / scale) and then
+ settles the pixel size to preserve the extent. For a 4x3 raster of 2x2
+ pixels (extent 8x6), a target pixel size of 1 gives an 8x6 output — the same
+ grid as the explicit dimensions, so the two modes agree pixel-for-pixel."""
+ tiff = _write(tmp_path, "uint8", bands=1, height=3, width=4, nodata=None)
+
+ by_scale = sedona.resample(tiff, scale_x=1.0, scale_y=-1.0)
+ expected = reference.resample(tiff, width=8, height=6)
+
+ assert by_scale.pixels.shape == (1, 6, 8)
+ np.testing.assert_array_equal(by_scale.pixels, expected.pixels)
+ assert by_scale.pixels.dtype == expected.pixels.dtype
+ assert_transform_and_nodata(by_scale, expected)
+
+
+def test_rs_resample_scale_mode_grows_extent_matches_rasterio(
+ sedona, reference, tmp_path
+):
+ """Scale mode keeps the requested pixel size exact and grows the extent to
+ whole pixels (Sedona Spark semantics). A 4x3 raster of 2x2 pixels (extent
+ 8x6) at pixel size 5 gives ceil(8/5)=2 x ceil(6/5)=2 output pixels spanning
+ 10x10 — one grown pixel past the source on each axis — so the right/bottom
+ border reads back as the band nodata. SedonaDB regrids with a fractional
+ RasterIO window; the reference reproduces the same nearest sampling and
+ nodata fill with a same-CRS `rasterio.warp.reproject`, so the covered pixels
+ (and the nodata fill) match exactly.
+
+ The source carries a CRS because the rasterio reference reprojects (same CRS
+ on both sides) to reproduce the grid; RS_Resample handles the CRS-less case
+ too (covered by the Rust tests)."""
+ tiff = _write(
+ tmp_path, "uint8", bands=2, height=3, width=4, nodata=200.0, crs="EPSG:4326"
+ )
+
+ got = sedona.resample(tiff, scale_x=5.0, scale_y=-5.0)
+ expected = reference.resample(tiff, scale_x=5.0, scale_y=-5.0)
+
+ assert got.pixels.shape == (2, 2, 2)
+ np.testing.assert_array_equal(got.pixels, expected.pixels)
+ assert got.pixels.dtype == expected.pixels.dtype
+ assert_transform_and_nodata(got, expected)
+
+
+def test_rs_resample_reference_raster_matches_dimension_mode(
+ sedona, reference, tmp_path
+):
+ """The 4-argument reference overload (useScale=false) takes the reference's
+ dimensions and origin. With the reference sharing the source's origin and
+ extent, this is the same grid as an explicit dimension resample, so it agrees
+ with the reference engine's dimension mode pixel-for-pixel."""
+ tiff = _write(
+ tmp_path, "uint8", bands=2, height=3, width=4, nodata=None, crs="EPSG:4326"
+ )
+ # Reference: 8x6 grid over the same extent/origin as the 4x3 source.
+ ref = tmp_path / "ref.tif"
+ write_grid_geotiff(
+ ref,
+ gdal_transform=(100.0, 1.0, 0.0, 500.0, 0.0, -1.0),
+ width=8,
+ height=6,
+ crs="EPSG:4326",
+ )
+
+ got = sedona.resample_to_reference(tiff, ref)
+ expected = reference.resample(tiff, width=8, height=6)
+
+ assert got.pixels.shape == (2, 6, 8)
+ np.testing.assert_array_equal(got.pixels, expected.pixels)
+ assert got.pixels.dtype == expected.pixels.dtype
+ assert_transform_and_nodata(got, expected)
+
+
+def test_rs_resample_reference_raster_crs_mismatch_errors(con, tmp_path):
+ """A reference raster in a different CRS errors — RS_Resample never
+ reprojects."""
+ pytest.importorskip("rasterio") # _write / write_grid_geotiff need rasterio
+ tiff = _write(
+ tmp_path, "uint8", bands=1, height=3, width=4, crs="EPSG:4326", name="in"
+ )
+ ref = tmp_path / "ref.tif"
+ write_grid_geotiff(
+ ref,
+ gdal_transform=(0.0, 1.0, 0.0, 6.0, 0.0, -1.0),
+ width=8,
+ height=6,
+ crs="EPSG:3857",
+ )
+ with pytest.raises(Exception, match="does not reproject"):
+ SedonaDB(con).resample_to_reference(tiff, ref)
+
+
+def test_rs_resample_null_raster_is_null(con):
+ """A NULL raster row yields a NULL result rather than erroring."""
+ table = pa.table(
+ {
+ "path": pa.array([None], type=pa.utf8()),
+ "width_or_scale": pa.array([8.0], type=pa.float64()),
+ "height_or_scale": pa.array([6.0], type=pa.float64()),
+ "use_scale": pa.array([False], type=pa.bool_()),
+ "algorithm": pa.array(["NearestNeighbor"], type=pa.utf8()),
+ }
+ )
+ df = con.create_data_frame(table)
+ result = df.select(
+ r=df.path.funcs.rs_frompath().funcs.rs_resample(
+ df.width_or_scale, df.height_or_scale, df.use_scale, df.algorithm
+ )
+ ).to_arrow_table()["r"]
+ assert decode_raster(result[0]) is None
+
+
@pytest.mark.parametrize(
- ("width", "height"), [(14, 12), (4, 3)], ids=["upsample", "downsample"]
+ ("width_or_scale", "height_or_scale", "use_scale", "algorithm", "match"),
+ [
+ (3.5, 6.0, False, "NearestNeighbor", "whole number"),
+ (0.0, 6.0, False, "NearestNeighbor", "positive"),
+ (0.0, -1.0, True, "NearestNeighbor", "non-zero"),
+ (8.0, 6.0, False, "sinc", "unknown algorithm"),
+ ],
)
-def test_rs_resample_nearest_matches_comparators(
- subject, comparator, tmp_path, dtype, width, height
+def test_rs_resample_argument_errors(
+ con, tmp_path, width_or_scale, height_or_scale, use_scale, algorithm, match
):
- """Nearest-neighbor picks source pixels verbatim over the unchanged
- extent — the planted dtype extremes must survive when selected, and the
- output transform is the extent divided by the new shape (including the
- non-integer 7/4 column scale)."""
- tiff = tmp_path / f"resample_{dtype}.tif"
- write_geotiff(
- tiff,
- random_raster_data(dtype, bands=2, height=HEIGHT, width=WIDTH),
- gdal_transform=GDAL_TRANSFORM,
+ """Invalid argument combinations raise with a clear message (executed through
+ the real array path, not constant-folded)."""
+ tiff = _write_uint8(tmp_path)
+ table = pa.table(
+ {
+ "path": pa.array([str(tiff)], type=pa.utf8()),
+ "width_or_scale": pa.array([width_or_scale], type=pa.float64()),
+ "height_or_scale": pa.array([height_or_scale], type=pa.float64()),
+ "use_scale": pa.array([use_scale], type=pa.bool_()),
+ "algorithm": pa.array([algorithm], type=pa.utf8()),
+ }
)
+ df = con.create_data_frame(table)
+ with pytest.raises(Exception, match=match):
+ df.select(
+ r=df.path.funcs.rs_frompath().funcs.rs_resample(
+ df.width_or_scale, df.height_or_scale, df.use_scale, df.algorithm
+ )
+ ).to_arrow_table()
+
+
+def test_rs_resample_sql_smoke(con, tmp_path):
+ """One SQL-text invocation keeps the parser path covered (everything else
+ routes through the expression API)."""
+ tiff = _write_uint8(tmp_path)
+ tab = con.sql(
+ "SELECT RS_Width(RS_Resample(RS_FromPath($1), 8, 6, false, 'NearestNeighbor')) AS w",
+ params=(str(tiff),),
+ ).to_arrow_table()
+ assert tab["w"][0].as_py() == 8
+
- got = subject.resample(tiff, width=width, height=height)
- expected = comparator.resample(tiff, width=width, height=height)
- assert_decoded_equal(got, expected, context=(dtype, width, height))
+def _write_uint8(tmp_path):
+ pytest.importorskip("rasterio") # write_random_geotiff needs rasterio
+ return _write(tmp_path, "uint8", bands=1, height=3, width=4, nodata=None)
diff --git a/rust/sedona-raster-gdal/Cargo.toml b/rust/sedona-raster-gdal/Cargo.toml
index 7c40103bf..591eaef94 100644
--- a/rust/sedona-raster-gdal/Cargo.toml
+++ b/rust/sedona-raster-gdal/Cargo.toml
@@ -95,6 +95,11 @@ harness = false
name = "rs_as_raster"
path = "benches/rs_as_raster.rs"
+[[bench]]
+harness = false
+name = "rs_resample"
+path = "benches/rs_resample.rs"
+
[[bench]]
harness = false
name = "rs_tile"
diff --git a/rust/sedona-raster-gdal/benches/rs_resample.rs b/rust/sedona-raster-gdal/benches/rs_resample.rs
new file mode 100644
index 000000000..035f8d212
--- /dev/null
+++ b/rust/sedona-raster-gdal/benches/rs_resample.rs
@@ -0,0 +1,130 @@
+// 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.
+
+//! Benchmarks for the RS_Resample UDF.
+//!
+//! RS_Resample reads the full source raster through GDAL's RasterIO resampling
+//! into a new grid. The cost is dominated by the resampling kernel over the
+//! larger of the source/output pixel counts, so the axes are:
+//!
+//! - **Resolution sweep** downsampling 2:1 with nearest neighbour, over
+//! `256²`, `512²`, `1024²` sources — the kernel scans O(source pixels).
+//! - **Upsample 1:2** with nearest neighbour, where the output pixel count
+//! dominates.
+//! - **Bilinear vs nearest** at a fixed resolution, isolating the kernel cost.
+//!
+//! The ground-truth comparison for correctness lives in the Python parity test
+//! (`test_rs_resample.py`), which cross-checks these same operations against
+//! rasterio's `Dataset.read(out_shape=..., resampling=...)` — the identical
+//! GDAL RasterIO path — so nearest is bit-exact and bilinear matches within a
+//! float tolerance.
+
+use std::sync::Arc;
+
+use arrow_array::{ArrayRef, BooleanArray, Float64Array, StringArray};
+use arrow_schema::DataType;
+use criterion::{criterion_group, criterion_main, Criterion};
+use datafusion_expr::ScalarUDF;
+use sedona_schema::datatypes::{SedonaType, RASTER};
+use sedona_testing::{raster_spec::RasterSpec, testers::ScalarUdfTester};
+
+fn criterion_benchmark(c: &mut Criterion) {
+ let f = sedona_raster_gdal::register::default_function_set();
+ let udf: ScalarUDF = f
+ .scalar_udf("rs_resample")
+ .expect("rs_resample is registered")
+ .clone()
+ .into();
+
+ // The 5-argument dimension overload:
+ // RS_Resample(raster, widthOrScale, heightOrScale, useScale, algorithm).
+ let tester = ScalarUdfTester::new(
+ udf,
+ vec![
+ RASTER,
+ SedonaType::Arrow(DataType::Float64),
+ SedonaType::Arrow(DataType::Float64),
+ SedonaType::Arrow(DataType::Boolean),
+ SedonaType::Arrow(DataType::Utf8),
+ ],
+ );
+
+ // A north-up single-band raster covering [0, size]² at unit pixels.
+ let build_raster = |size: i64| -> ArrayRef {
+ let values = vec![1u8; (size * size) as usize];
+ Arc::new(
+ RasterSpec::d2(size, size)
+ .band_values(&values)
+ .crs(None)
+ .bbox(0.0, 0.0, size as f64, size as f64)
+ .build(),
+ )
+ };
+
+ // Dimension mode (useScale=false), so widthOrScale/heightOrScale are the
+ // output pixel dimensions.
+ let run =
+ |c: &mut Criterion, label: &str, raster: ArrayRef, width: f64, height: f64, alg: &str| {
+ let width: ArrayRef = Arc::new(Float64Array::from(vec![width]));
+ let height: ArrayRef = Arc::new(Float64Array::from(vec![height]));
+ let use_scale: ArrayRef = Arc::new(BooleanArray::from(vec![false]));
+ let algorithm: ArrayRef = Arc::new(StringArray::from(vec![alg.to_string()]));
+ c.bench_function(label, |b| {
+ b.iter(|| {
+ tester
+ .invoke_arrays(vec![
+ raster.clone(),
+ width.clone(),
+ height.clone(),
+ use_scale.clone(),
+ algorithm.clone(),
+ ])
+ .unwrap()
+ })
+ });
+ };
+
+ // Resolution sweep: 2:1 nearest downsample.
+ for size in [256i64, 512, 1024] {
+ let half = (size / 2) as f64;
+ let label = format!("raster-gdal rs_resample Downsample2x(Nearest, Raster({size}x{size}))");
+ run(c, &label, build_raster(size), half, half, "NearestNeighbor");
+ }
+
+ // 1:2 nearest upsample: output pixel count dominates.
+ run(
+ c,
+ "raster-gdal rs_resample Upsample2x(Nearest, Raster(512x512))",
+ build_raster(512),
+ 1024.0,
+ 1024.0,
+ "NearestNeighbor",
+ );
+
+ // Bilinear vs nearest at a fixed resolution isolates the kernel cost.
+ run(
+ c,
+ "raster-gdal rs_resample Downsample2x(Bilinear, Raster(1024x1024))",
+ build_raster(1024),
+ 512.0,
+ 512.0,
+ "Bilinear",
+ );
+}
+
+criterion_group!(benches, criterion_benchmark);
+criterion_main!(benches);
diff --git a/rust/sedona-raster-gdal/src/lib.rs b/rust/sedona-raster-gdal/src/lib.rs
index 376e5566e..28bb22a92 100644
--- a/rust/sedona-raster-gdal/src/lib.rs
+++ b/rust/sedona-raster-gdal/src/lib.rs
@@ -41,6 +41,7 @@ mod rs_frompath;
mod rs_metadata;
mod rs_polygonize;
mod rs_reproject_match;
+mod rs_resample;
mod rs_tile;
mod source_uri;
mod utils;
@@ -59,6 +60,7 @@ pub use rs_frompath::rs_frompath_udf;
pub use rs_metadata::rs_metadata_udf;
pub use rs_polygonize::rs_polygonize_udf;
pub use rs_reproject_match::rs_reproject_match_udf;
+pub use rs_resample::rs_resample_udf;
pub use rs_tile::rs_tile_udf;
pub use utils::{
append_as_indb_raster, append_as_outdb_raster, append_nd_from_dataset, dataset_to_indb_raster,
diff --git a/rust/sedona-raster-gdal/src/register.rs b/rust/sedona-raster-gdal/src/register.rs
index f0d0c20b5..c2fb8d5ad 100644
--- a/rust/sedona-raster-gdal/src/register.rs
+++ b/rust/sedona-raster-gdal/src/register.rs
@@ -28,6 +28,7 @@ pub fn default_function_set() -> FunctionSet {
function_set.insert_scalar_udf(crate::rs_metadata::rs_metadata_udf());
function_set.insert_scalar_udf(crate::rs_polygonize::rs_polygonize_udf());
function_set.insert_scalar_udf(crate::rs_reproject_match::rs_reproject_match_udf());
+ function_set.insert_scalar_udf(crate::rs_resample::rs_resample_udf());
function_set.insert_scalar_udf(crate::rs_tile::rs_tile_udf());
function_set
}
diff --git a/rust/sedona-raster-gdal/src/rs_reproject_match.rs b/rust/sedona-raster-gdal/src/rs_reproject_match.rs
index 2ce57b466..e3310c6a1 100644
--- a/rust/sedona-raster-gdal/src/rs_reproject_match.rs
+++ b/rust/sedona-raster-gdal/src/rs_reproject_match.rs
@@ -55,11 +55,13 @@ use sedona_raster_functions::rs_ensure_loaded::{
use sedona_raster_functions::RasterExecutor;
use sedona_schema::datatypes::{SedonaType, RASTER};
use sedona_schema::matchers::ArgMatcher;
-use sedona_schema::raster::BandDataType;
use crate::gdal_common::{raster_ref_to_gdal_mem, with_gdal, GdalBandLayout};
use crate::gdal_dataset_provider::configure_thread_local_options;
-use crate::utils::{append_warped_nd_from_dataset, WarpGrid};
+use crate::utils::{
+ append_warped_nd_from_dataset, parse_resample_algorithm, reject_lossy_resample_dtypes, Grid,
+ OutputGrid,
+};
/// RS_ReprojectMatch() scalar UDF implementation.
///
@@ -168,7 +170,7 @@ impl SedonaScalarKernel for RsReprojectMatch {
return Ok(());
};
- let alg = parse_algorithm(algorithm)?;
+ let alg = parse_resample_algorithm(algorithm, "RS_ReprojectMatch")?;
reproject_match(
gdal,
raster,
@@ -255,20 +257,11 @@ fn reproject_match(
let band_count = raster.bands().len();
- // GDAL's warp routes 64-bit integer pixels through a floating working type
- // (a double for nearest/interpolation, a 32-bit float for mode on GDAL <
- // 3.13), so no resampling method preserves Int64/UInt64 exactly. Reject
- // these dtypes up front rather than mis-warping silently.
- for i in 0..band_count {
- if let Some(BandDataType::Int64 | BandDataType::UInt64) = raster.band_data_type(i) {
- return exec_err!(
- "RS_ReprojectMatch does not support Int64/UInt64 rasters: GDAL's warp routes \
- 64-bit integer pixels through a floating working type (a double for \
- nearest/interpolation, a 32-bit float for mode on GDAL < 3.13), which cannot \
- represent them exactly. Cast to a supported type (e.g. Int32/Float64) first."
- );
- }
- }
+ // A warp always routes pixels through a floating working type, so no
+ // resampling method preserves Int64/UInt64 exactly (nor Int32/UInt32 + Mode on
+ // GDAL < 3.13). Reject those up front rather than mis-warping silently.
+ // RS_Resample shares this check.
+ reject_lossy_resample_dtypes(gdal, raster, band_count, alg, true, "RS_ReprojectMatch")?;
let band_indices: Vec<usize> = (1..=band_count).collect();
let layout = GdalBandLayout::from_raster(raster, &band_indices)?;
@@ -277,45 +270,23 @@ fn reproject_match(
// returns, so `raster` outlives it.
let src_dataset = unsafe { raster_ref_to_gdal_mem(gdal, raster, &band_indices)? };
- let warp_grid = WarpGrid {
- transform,
- width: ref_width as usize,
- height: ref_height as usize,
+ let output = OutputGrid {
+ grid: Grid {
+ transform,
+ width: ref_width,
+ height: ref_height,
+ },
crs: reference.crs(),
alg,
- warp_memory_limit_bytes,
};
- append_warped_nd_from_dataset(gdal, &src_dataset, &layout, builder, &warp_grid)
-}
-
-/// Map an algorithm name (case-insensitive) to a GDAL [`ResampleAlg`]. An empty
-/// string defaults to nearest neighbour (matching Sedona Spark).
-///
-/// Accepts the GDAL names plus the Spark spellings (American `NearestNeighbor`
-/// and `Bicubic`, which GDAL calls `Cubic`). The algorithms without a warp
-/// equivalent (Gauss) are not accepted.
-fn parse_algorithm(name: &str) -> Result<ResampleAlg> {
- if name.is_empty() {
- return Ok(ResampleAlg::NearestNeighbour);
- }
- let alg = match name.to_ascii_lowercase().as_str() {
- "nearestneighbor" | "nearestneighbour" | "nearest" | "near" => {
- ResampleAlg::NearestNeighbour
- }
- "bilinear" => ResampleAlg::Bilinear,
- "cubic" | "bicubic" => ResampleAlg::Cubic,
- "cubicspline" => ResampleAlg::CubicSpline,
- "lanczos" => ResampleAlg::Lanczos,
- "average" => ResampleAlg::Average,
- "mode" => ResampleAlg::Mode,
- _ => {
- return exec_err!(
- "RS_ReprojectMatch: unknown algorithm {name:?}; expected one of \
- NearestNeighbor, Bilinear, Cubic, CubicSpline, Lanczos, Average, Mode"
- );
- }
- };
- Ok(alg)
+ append_warped_nd_from_dataset(
+ gdal,
+ &src_dataset,
+ &layout,
+ builder,
+ &output,
+ warp_memory_limit_bytes,
+ )
}
#[cfg(test)]
@@ -610,4 +581,40 @@ mod tests {
}
}
}
+
+ #[test]
+ fn int32_mode_gated_on_gdal_version() {
+ // Mode uses a 32-bit float working type on GDAL < 3.13 (a double from 3.13
+ // on), which cannot represent Int32/UInt32 above 2^24. Mode is a selection,
+ // so silently rounding a category code is wrong: reject Int32 + Mode on
+ // older GDAL, accept it (double working type) on GDAL >= 3.13. Other
+ // algorithms are unaffected — a double is exact for 32-bit ints.
+ let reference = RasterSpec::d2(2, 2)
+ .band_values(&[0u8; 4])
+ .bbox(0.0, 0.0, 2.0, 2.0)
+ .crs(Some("EPSG:4326"));
+ let big: i32 = (1i32 << 24) + 1;
+ let int32_source = RasterSpec::d2(2, 2)
+ .band_values(&[big, big, big, big])
+ .bbox(0.0, 0.0, 2.0, 2.0)
+ .crs(Some("EPSG:4326"));
+
+ // Bilinear routes through a double (exact for 32-bit ints), so Int32 is
+ // accepted on every GDAL version — the gate is Mode-specific.
+ tester()
+ .invoke_scalar_scalar_scalar(&int32_source, &reference, "Bilinear")
+ .unwrap();
+
+ let gdal_version = with_gdal(|gdal| Ok(gdal.version_num())).unwrap();
+ let mode = tester().invoke_scalar_scalar_scalar(&int32_source, &reference, "Mode");
+ if gdal_version >= 3_130_000 {
+ mode.unwrap();
+ } else {
+ let err = mode.unwrap_err().to_string();
+ assert!(
+ err.contains("Mode resampling of Int32/UInt32"),
+ "got: {err}"
+ );
+ }
+ }
}
diff --git a/rust/sedona-raster-gdal/src/rs_resample.rs b/rust/sedona-raster-gdal/src/rs_resample.rs
new file mode 100644
index 000000000..0020102ac
--- /dev/null
+++ b/rust/sedona-raster-gdal/src/rs_resample.rs
@@ -0,0 +1,1344 @@
+// 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_Resample UDF - Resample a raster to a new resolution/grid.
+//!
+//! Changes a raster's pixel grid — its width/height or its pixel size — and,
+//! optionally, snaps the output origin to a grid (`gridX`/`gridY`). Band
+//! count/order and per-band nodata are preserved. Pixel values are recomputed by
+//! GDAL: a plain resolution change reads the source into the new grid with
+//! RasterIO resampling, while a scale change that grows the extent or an origin
+//! snap go through GDAL's warp (which maps every output pixel centre exactly),
+//! leaving the grown/shifted border at nodata.
+//!
+//! `RS_Resample` is a **same-CRS** operation — it never reprojects. The
+//! reference-raster overload requires the reference to share the input's CRS and
+//! errors otherwise (matching Apache Sedona (Spark)).
+//!
+//! Dimension mode preserves the world extent (the pixel size is derived from
+//! it). Scale mode keeps the requested pixel size exact and grows the extent up
+//! to one pixel so it tiles into whole pixels — matching Sedona Spark's
+//! `RS_Resample`.
+//!
+//! The argument surface mirrors Sedona Spark's positional overloads verbatim, so
+//! Spark SQL tends to run unchanged: a reference-raster form, a `(widthOrScale,
+//! heightOrScale)` form, and a `(widthOrScale, heightOrScale, gridX, gridY)`
+//! form, each ending in `(useScale, algorithm)`.
+
+use std::sync::Arc;
+
+use arrow_array::ArrayRef;
+use arrow_schema::DataType;
+use datafusion_common::cast::{as_boolean_array, as_float64_array, as_string_array};
+use datafusion_common::config::ConfigOptions;
+use datafusion_common::error::Result;
+use datafusion_common::exec_err;
+use datafusion_expr::{ColumnarValue, Volatility};
+
+use sedona_common::sedona_internal_err;
+use sedona_expr::scalar_udf::{SedonaScalarKernel, SedonaScalarUDF};
+use sedona_gdal::gdal::Gdal;
+use sedona_gdal::geo_transform::{GeoTransform, GeoTransformEx};
+use sedona_gdal::raster::types::ResampleAlg;
+use sedona_raster::array::RasterRefImpl;
+use sedona_raster::builder::RasterBuilder;
+use sedona_raster::traits::RasterRef;
+use sedona_raster_functions::crs_utils::{crs_transform_required, resolve_crs};
+use sedona_raster_functions::rs_ensure_loaded::{
+ NEEDS_PIXELS_METADATA_KEY, RETURNS_BYTES_METADATA_KEY,
+};
+use sedona_raster_functions::RasterExecutor;
+use sedona_schema::datatypes::{SedonaType, RASTER};
+use sedona_schema::matchers::ArgMatcher;
+
+use crate::gdal_common::{raster_ref_to_gdal_mem, with_gdal, GdalBandLayout};
+use crate::gdal_dataset_provider::configure_thread_local_options;
+use crate::utils::{
+ append_resampled_nd_from_dataset, append_warped_nd_from_dataset, parse_resample_algorithm,
+ reject_lossy_resample_dtypes, Grid, OutputGrid,
+};
+
+/// RS_Resample() scalar UDF implementation.
+///
+/// Resamples a raster onto a new pixel grid in the same CRS. The positional
+/// overloads match Apache Sedona (Spark):
+/// - `RS_Resample(raster, referenceRaster, useScale, algorithm)` — 4 args
+/// - `RS_Resample(raster, widthOrScale, heightOrScale, useScale, algorithm)` — 5 args
+/// - `RS_Resample(raster, widthOrScale, heightOrScale, gridX, gridY, useScale, algorithm)` — 7 args
+pub fn rs_resample_udf() -> SedonaScalarUDF {
+ SedonaScalarUDF::new(
+ "rs_resample",
+ vec![
+ Arc::new(RsResample { arg_count: 4 }),
+ Arc::new(RsResample { arg_count: 5 }),
+ Arc::new(RsResample { arg_count: 7 }),
+ ],
+ Volatility::Immutable,
+ )
+ // Reads band pixels (so the planner materializes OutDb rasters via
+ // RS_EnsureLoaded first) and emits a fresh InDb raster (so its output is
+ // already loaded and isn't wrapped again).
+ .with_metadata(NEEDS_PIXELS_METADATA_KEY, "true")
+ .with_metadata(RETURNS_BYTES_METADATA_KEY, "true")
+}
+
+/// Kernel implementation for RS_Resample.
+///
+/// The `useScale` boolean and the 4/5/7-argument overload ladder below are
+/// Sedona Spark's `RS_Resample` signature reproduced verbatim, so Spark SQL
+/// runs unchanged. That parity is why the modes are one function switching on
+/// `arg_count` / `useScale` rather than separately-named functions: the shared
+/// SQL surface is deliberate, not an accident of implementation.
+#[derive(Debug)]
+struct RsResample {
+ /// Number of arguments in the matched signature (4, 5, or 7).
+ arg_count: usize,
+}
+
+impl SedonaScalarKernel for RsResample {
+ fn return_type(&self, args: &[SedonaType]) -> Result<Option<SedonaType>> {
+ let matchers = match self.arg_count {
+ // (raster, referenceRaster, useScale, algorithm)
+ 4 => vec![
+ ArgMatcher::is_raster(),
+ ArgMatcher::is_raster(),
+ ArgMatcher::is_boolean(),
+ ArgMatcher::is_string(),
+ ],
+ // (raster, widthOrScale, heightOrScale, useScale, algorithm)
+ 5 => vec![
+ ArgMatcher::is_raster(),
+ ArgMatcher::is_numeric(),
+ ArgMatcher::is_numeric(),
+ ArgMatcher::is_boolean(),
+ ArgMatcher::is_string(),
+ ],
+ // (raster, widthOrScale, heightOrScale, gridX, gridY, useScale, algorithm)
+ 7 => vec![
+ ArgMatcher::is_raster(),
+ ArgMatcher::is_numeric(),
+ ArgMatcher::is_numeric(),
+ ArgMatcher::is_numeric(),
+ ArgMatcher::is_numeric(),
+ ArgMatcher::is_boolean(),
+ ArgMatcher::is_string(),
+ ],
+ _ => {
+ return sedona_internal_err!(
+ "RS_Resample: unexpected arg_count {}",
+ self.arg_count
+ );
+ }
+ };
+ ArgMatcher::new(matchers, RASTER).match_args(args)
+ }
+
+ fn invoke_batch(
+ &self,
+ arg_types: &[SedonaType],
+ args: &[ColumnarValue],
+ ) -> Result<ColumnarValue> {
+ self.invoke_batch_from_args(arg_types, args, &SedonaType::Arrow(DataType::Null), 0, None)
+ }
+
+ fn invoke_batch_from_args(
+ &self,
+ arg_types: &[SedonaType],
+ args: &[ColumnarValue],
+ _return_type: &SedonaType,
+ _num_rows: usize,
+ config_options: Option<&ConfigOptions>,
+ ) -> Result<ColumnarValue> {
+ match self.arg_count {
+ 4 => self.invoke_reference(arg_types, args, config_options),
+ 5 | 7 => self.invoke_dims_or_scale(arg_types, args, config_options),
+ other => sedona_internal_err!("RS_Resample: unexpected arg_count {other}"),
+ }
+ }
+}
+
+impl RsResample {
+ /// `RS_Resample(raster, referenceRaster, useScale, algorithm)`: the target
+ /// grid (dimensions or pixel size) and origin come from the reference raster,
+ /// which must share the input's CRS.
+ fn invoke_reference(
+ &self,
+ arg_types: &[SedonaType],
+ args: &[ColumnarValue],
+ config_options: Option<&ConfigOptions>,
+ ) -> Result<ColumnarValue> {
+ let num_iterations = RasterExecutor::num_iterations_over(args);
+
+ let use_scale_array = boolean_column(&args[2], num_iterations)?;
+ let mut use_scale_iter = use_scale_array.iter();
+ let algorithm_array = string_column(&args[3], num_iterations)?;
+ let mut algorithm_iter = algorithm_array.iter();
+
+ let mut builder = RasterBuilder::new(num_iterations);
+
+ // The executor iterates the (raster, referenceRaster) pair; useScale and
+ // algorithm are advanced in lockstep below.
+ let exec_arg_types = vec![arg_types[0].clone(), arg_types[1].clone()];
+ let exec_args = vec![args[0].clone(), args[1].clone()];
+ let executor =
+ RasterExecutor::new_with_num_iterations(&exec_arg_types, &exec_args, num_iterations);
+
+ with_gdal(|gdal| {
+ configure_thread_local_options(gdal, config_options)?;
+ executor.execute_raster_raster_void(|_i, raster_opt, reference_opt| {
+ let use_scale = use_scale_iter.next().flatten();
+ let algorithm = algorithm_iter.next().flatten();
+ let (Some(raster), Some(reference), Some(use_scale), Some(algorithm)) =
+ (raster_opt, reference_opt, use_scale, algorithm)
+ else {
+ builder.append_null()?;
+ return Ok(());
+ };
+
+ // Same-CRS only: RS_Resample never reprojects (matching Spark,
+ // which throws when the reference SRID differs). Compare CRSes
+ // semantically rather than by raw string, so two equivalent
+ // encodings of the same CRS are not spuriously rejected; a CRS on
+ // exactly one side is likewise an error (via crs_transform_required).
+ let raster_crs = resolve_crs(raster.crs())?;
+ let reference_crs = resolve_crs(reference.crs())?;
+ if crs_transform_required(
+ raster_crs.as_deref(),
+ reference_crs.as_deref(),
+ "input raster",
+ "referenceRaster",
+ )? {
+ return exec_err!(
+ "RS_Resample: referenceRaster CRS differs from the input raster CRS; \
+ RS_Resample does not reproject"
+ );
+ }
+
+ let plan = plan_from_reference(reference, use_scale, algorithm)?;
+ resample_raster(gdal, raster, &plan, &mut builder)?;
+ Ok(())
+ })?;
+
+ let out: ArrayRef = Arc::new(builder.finish()?);
+ RasterExecutor::finish_over(args, out)
+ })
+ }
+
+ /// The 5- and 7-argument overloads: the target grid comes from
+ /// `(widthOrScale, heightOrScale)`, with an optional origin snap to
+ /// `(gridX, gridY)` in the 7-argument form.
+ fn invoke_dims_or_scale(
+ &self,
+ arg_types: &[SedonaType],
+ args: &[ColumnarValue],
+ config_options: Option<&ConfigOptions>,
+ ) -> Result<ColumnarValue> {
+ let num_iterations = RasterExecutor::num_iterations_over(args);
+
+ // Column indices differ between the 5- and 7-argument overloads.
+ let has_grid = self.arg_count == 7;
+ let (use_scale_idx, algorithm_idx) = if has_grid { (5, 6) } else { (3, 4) };
+
+ let width_array = float_column(&args[1], num_iterations)?;
+ let height_array = float_column(&args[2], num_iterations)?;
+ let mut width_iter = width_array.iter();
+ let mut height_iter = height_array.iter();
+
+ let grid_x_array = if has_grid {
+ Some(float_column(&args[3], num_iterations)?)
+ } else {
+ None
+ };
+ let grid_y_array = if has_grid {
+ Some(float_column(&args[4], num_iterations)?)
+ } else {
+ None
+ };
+ let mut grid_x_iter = grid_x_array.as_ref().map(|a| a.iter());
+ let mut grid_y_iter = grid_y_array.as_ref().map(|a| a.iter());
+
+ let use_scale_array = boolean_column(&args[use_scale_idx], num_iterations)?;
+ let mut use_scale_iter = use_scale_array.iter();
+ let algorithm_array = string_column(&args[algorithm_idx], num_iterations)?;
+ let mut algorithm_iter = algorithm_array.iter();
+
+ let mut builder = RasterBuilder::new(num_iterations);
+
+ let exec_arg_types = vec![arg_types[0].clone()];
+ let exec_args = vec![args[0].clone()];
+ let executor =
+ RasterExecutor::new_with_num_iterations(&exec_arg_types, &exec_args, num_iterations);
+
+ with_gdal(|gdal| {
+ configure_thread_local_options(gdal, config_options)?;
+ executor.execute_raster_void(|_i, raster_opt| {
+ let width_or_scale = width_iter.next().flatten();
+ let height_or_scale = height_iter.next().flatten();
+ let grid_x = grid_x_iter.as_mut().and_then(|it| it.next().flatten());
+ let grid_y = grid_y_iter.as_mut().and_then(|it| it.next().flatten());
+ let use_scale = use_scale_iter.next().flatten();
+ let algorithm = algorithm_iter.next().flatten();
+
+ let (
+ Some(raster),
+ Some(width_or_scale),
+ Some(height_or_scale),
+ Some(use_scale),
+ Some(algorithm),
+ ) = (
+ raster_opt,
+ width_or_scale,
+ height_or_scale,
+ use_scale,
+ algorithm,
+ )
+ else {
+ builder.append_null()?;
+ return Ok(());
+ };
+
+ let grid = match (has_grid, grid_x, grid_y) {
+ (false, _, _) => None,
+ (true, Some(gx), Some(gy)) => Some((gx, gy)),
+ // A NULL gridX/gridY in the 7-arg form yields a NULL row.
+ (true, _, _) => {
+ builder.append_null()?;
+ return Ok(());
+ }
+ };
+
+ let plan = plan_from_dims_or_scale(
+ width_or_scale,
+ height_or_scale,
+ grid,
+ use_scale,
+ algorithm,
+ )?;
+ resample_raster(gdal, raster, &plan, &mut builder)?;
+ Ok(())
+ })?;
+
+ let out: ArrayRef = Arc::new(builder.finish()?);
+ RasterExecutor::finish_over(args, out)
+ })
+ }
+}
+
+/// Cast one argument to a `Float64` array (expanding a scalar to `num_iterations`
+/// rows), so a per-row column and a scalar are iterated the same way.
+fn float_column(arg: &ColumnarValue, num_iterations: usize) -> Result<arrow_array::Float64Array> {
+ let array = arg
+ .clone()
+ .cast_to(&DataType::Float64, None)?
+ .into_array(num_iterations)?;
+ Ok(as_float64_array(&array)?.clone())
+}
+
+/// Cast one argument to a `Boolean` array (expanding a scalar to `num_iterations`
+/// rows).
+fn boolean_column(arg: &ColumnarValue, num_iterations: usize) -> Result<arrow_array::BooleanArray> {
+ let array = arg
+ .clone()
+ .cast_to(&DataType::Boolean, None)?
+ .into_array(num_iterations)?;
+ Ok(as_boolean_array(&array)?.clone())
+}
+
+/// Cast one argument to a `Utf8` array (expanding a scalar to `num_iterations`
+/// rows).
+fn string_column(arg: &ColumnarValue, num_iterations: usize) -> Result<arrow_array::StringArray> {
+ let array = arg
+ .clone()
+ .cast_to(&DataType::Utf8, None)?
+ .into_array(num_iterations)?;
+ Ok(as_string_array(&array)?.clone())
+}
+
+/// Which axis the output grid is pinned on.
+#[derive(Debug, Clone, Copy, PartialEq)]
+enum TargetGrid {
+ /// Fixed output dimensions; the pixel size is derived so the world extent
+ /// is preserved.
+ Dims { width: i64, height: i64 },
+ /// Fixed target pixel size; the output dimensions are `ceil(extent / scale)`
+ /// and the extent grows up to one pixel so it tiles into whole pixels
+ /// (matching Sedona Spark).
+ Scale { scale_x: f64, scale_y: f64 },
+}
+
+/// A validated resample request: the target grid, an optional origin snap, and
+/// the resampling algorithm.
+#[derive(Debug, Clone, PartialEq)]
+struct ResamplePlan {
+ target: TargetGrid,
+ /// `(grid_x, grid_y)` if the output origin is snapped to a grid.
+ grid_snap: Option<(f64, f64)>,
+ algorithm: ResampleAlg,
+}
+
+/// Build a plan from a `(widthOrScale, heightOrScale)` pair and optional origin
+/// snap. When `use_scale` is false the pair is target dimensions (rejecting a
+/// fractional width/height); when true it is the target pixel size.
+fn plan_from_dims_or_scale(
+ width_or_scale: f64,
+ height_or_scale: f64,
+ grid: Option<(f64, f64)>,
+ use_scale: bool,
+ algorithm: &str,
+) -> Result<ResamplePlan> {
+ let algorithm = parse_resample_algorithm(algorithm, "RS_Resample")?;
+
+ let target = if use_scale {
+ if !width_or_scale.is_finite()
+ || !height_or_scale.is_finite()
+ || width_or_scale == 0.0
+ || height_or_scale == 0.0
+ {
+ return exec_err!(
+ "RS_Resample: scaleX and scaleY must be finite and non-zero \
+ (got {width_or_scale}, {height_or_scale})"
+ );
+ }
+ TargetGrid::Scale {
+ scale_x: width_or_scale,
+ scale_y: height_or_scale,
+ }
+ } else {
+ let width = whole_dimension(width_or_scale, "width")?;
+ let height = whole_dimension(height_or_scale, "height")?;
+ TargetGrid::Dims { width, height }
+ };
+
+ let grid_snap = match grid {
+ Some((grid_x, grid_y)) => {
+ if !grid_x.is_finite() || !grid_y.is_finite() {
+ return exec_err!("RS_Resample: gridX and gridY must be finite");
+ }
+ Some((grid_x, grid_y))
+ }
+ None => None,
+ };
+
+ Ok(ResamplePlan {
+ target,
+ grid_snap,
+ algorithm,
+ })
+}
+
+/// Build a plan from a reference raster: its dimensions (or pixel size, when
+/// `use_scale`) and its upper-left origin as the snap grid. The core grid math is
+/// shared with [`plan_from_dims_or_scale`].
+fn plan_from_reference(
+ reference: &RasterRefImpl<'_>,
+ use_scale: bool,
+ algorithm: &str,
+) -> Result<ResamplePlan> {
+ let m = reference.metadata();
+ let (width_or_scale, height_or_scale) = if use_scale {
+ (m.scale_x(), m.scale_y())
+ } else {
+ (m.width() as f64, m.height() as f64)
+ };
+ // The reference's upper-left corner is the grid the output origin snaps to.
+ let grid = Some((m.upper_left_x(), m.upper_left_y()));
+ plan_from_dims_or_scale(width_or_scale, height_or_scale, grid, use_scale, algorithm)
+}
+
+/// A whole-number dimension from a `Double`: rejects a non-finite, fractional, or
+/// non-positive value rather than silently truncating it.
+fn whole_dimension(value: f64, name: &str) -> Result<i64> {
+ if !value.is_finite() {
+ return exec_err!("RS_Resample: {name} must be finite (got {value})");
+ }
+ if value.fract() != 0.0 {
+ return exec_err!(
+ "RS_Resample: {name} must be a whole number when useScale is false (got {value})"
+ );
+ }
+ let dim = value as i64;
+ if dim <= 0 {
+ return exec_err!("RS_Resample: {name} must be positive (got {dim})");
+ }
+ if dim > i32::MAX as i64 {
+ return exec_err!("RS_Resample: {name} {dim} exceeds the maximum raster dimension");
+ }
+ Ok(dim)
+}
+
+/// Spark's `approximateEquals` tolerance for deciding whether the grid snap
+/// actually moved the origin. This is a world-coordinate distance, so on a
+/// lon/lat raster `1e-6` degrees is ~10 cm at the equator — below the accuracy
+/// anyone resampling in lon/lat is likely to care about; for projected rasters
+/// it is well below one pixel.
+const SNAP_EPS: f64 = 1e-6;
+
+/// Resample one raster into `builder` according to `plan`.
+///
+/// Band count/order and per-band nodata are preserved. A plain dimension change
+/// takes the extent-preserving RasterIO path; a scale change (which grows the
+/// extent) or an origin snap go through GDAL's warp, which maps every output
+/// pixel centre exactly and fills any output cells outside the source footprint
+/// with nodata. The resample is a 2-D `(y, x)` operation broadcast across every
+/// non-spatial
+/// plane of an N-D band via [`GdalBandLayout`].
+fn resample_raster(
+ gdal: &Gdal,
+ raster: &RasterRefImpl<'_>,
+ plan: &ResamplePlan,
+ builder: &mut RasterBuilder,
+) -> Result<()> {
+ let metadata = raster.metadata();
+ let src = Grid::from_metadata(&metadata);
+ if src.width <= 0 || src.height <= 0 {
+ return exec_err!(
+ "RS_Resample: source raster has non-positive dimensions {}x{}",
+ src.width,
+ src.height
+ );
+ }
+
+ let band_count = raster.bands().len();
+
+ // The regrid path — a scale change, origin snap, or reference grid — moves
+ // the output off the source grid and/or grows the extent past the source
+ // footprint, which a single axis-aligned RasterIO read cannot express, so it
+ // goes through GDAL's warp (which maps every output pixel centre exactly).
+ // A plain dimension change stays on the extent-preserving RasterIO
+ // decimation.
+ let is_regrid = matches!(plan.target, TargetGrid::Scale { .. }) || plan.grid_snap.is_some();
+
+ // Int64/UInt64 are preserved only on the nearest-neighbour RasterIO fast path
+ // (a plain dimension change), where each output pixel is a copied source
+ // sample — a value selection, not a computation — so it is bit-exact for any
+ // dtype whether up- or down-sampling (pinned by the
+ // nearest_fast_path_preserves_* tests). Every other case routes pixels through
+ // a floating working type — the warp path always, and RasterIO for any
+ // interpolating algorithm — so those dtypes are rejected (Int32/UInt32 + Mode
+ // only on GDAL < 3.13). RS_ReprojectMatch shares this check.
+ let routes_through_float = is_regrid || plan.algorithm != ResampleAlg::NearestNeighbour;
+ reject_lossy_resample_dtypes(
+ gdal,
+ raster,
+ band_count,
+ plan.algorithm,
+ routes_through_float,
+ "RS_Resample",
+ )?;
+
+ // A scale change or origin snap on a skewed (non-north-up) raster is an
+ // explicit unsupported case: the extent-preserving dimension path carries the
+ // shear through unchanged, but growing the extent or snapping the origin
+ // against a sheared grid is not supported.
+ if is_regrid && (src.transform[2] != 0.0 || src.transform[4] != 0.0) {
+ return exec_err!(
+ "RS_Resample: a scale change or grid snap on a skewed (non-north-up) raster is not \
+ supported (only an extent-preserving width/height change is)"
+ );
+ }
+
+ let band_indices: Vec<usize> = (1..=band_count).collect();
+ let layout = GdalBandLayout::from_raster(raster, &band_indices)?;
+ // SAFETY: the returned dataset references `raster`'s band bytes zero-copy;
+ // it is only read below and dropped before `resample_raster` returns, so
+ // `raster` outlives it.
+ let src_dataset = unsafe { raster_ref_to_gdal_mem(gdal, raster, &band_indices)? };
+
+ // Output geotransform matching Sedona Spark: the skew is carried through
+ // unchanged and the dimensions / pixel size follow the requested mode. The
+ // output stays in the source CRS — RS_Resample never reprojects.
+ let output = OutputGrid {
+ grid: compute_output_grid(&src, plan),
+ crs: raster.crs(),
+ alg: plan.algorithm,
+ };
+
+ if !is_regrid {
+ // Extent-preserving RasterIO decimation supplies the pixels (bit-exact
+ // for nearest; interpolating algorithms blend as GDAL RasterIO does).
+ return append_resampled_nd_from_dataset(&src_dataset, &layout, builder, &output);
+ }
+
+ // Warp path: scale-mode grow-extent, origin snap, and reference-grid
+ // resampling. Same CRS on both sides (RS_Resample never reprojects), so this
+ // is a pure regrid onto `output` — GDAL maps each output pixel centre exactly
+ // (no covered-window clamp) and fills the grown/shifted border, uncovered by
+ // the source, with each band's nodata. `0.0` uses GDAL's own default
+ // working-memory cache size.
+ append_warped_nd_from_dataset(gdal, &src_dataset, &layout, builder, &output, 0.0)
+}
+
+/// Resolve the output grid, mirroring Apache Sedona (Spark)'s
+/// `RasterEditors.resample`: dimension mode preserves the extent, scale mode
+/// keeps the pixel size exact and grows the extent, and `gridX`/`gridY` snap the
+/// origin outward to the grid. The CRS is unchanged from the source.
+fn compute_output_grid(src: &Grid, plan: &ResamplePlan) -> Grid {
+ let env = src.envelope();
+ let [src_ul_x, _, src_skew_x, src_ul_y, src_skew_y, src_scale_y] = src.transform;
+
+ // Dimensions and pixel size per the requested mode. Skew is carried through
+ // unchanged (Spark keeps the original shear in both modes).
+ let (mut width, mut height, mut scale_x, mut scale_y) = match plan.target {
+ TargetGrid::Dims { width, height } => {
+ let sx = (env.max_x - env.min_x) / width as f64;
+ let sy = src_scale_y.signum() * (env.max_y - env.min_y) / height as f64;
+ (width, height, sx, sy)
+ }
+ TargetGrid::Scale { scale_x, scale_y } => {
+ let width = ((env.max_x - env.min_x) / scale_x.abs()).ceil().max(1.0) as i64;
+ let height = ((env.max_y - env.min_y) / scale_y.abs()).ceil().max(1.0) as i64;
+ (width, height, scale_x, scale_y)
+ }
+ };
+ let skew_x = src_skew_x;
+ let skew_y = src_skew_y;
+ let mut upper_left_x = src_ul_x;
+ let mut upper_left_y = src_ul_y;
+
+ // Origin snap: move the origin to the grid-cell corner containing it
+ // (expanding outward so the original extent stays covered), then regrow the
+ // dimension (scale mode) or re-fit the pixel size (dimension mode).
+ if let Some((grid_x, grid_y)) = plan.grid_snap {
+ let current = [upper_left_x, scale_x, skew_x, upper_left_y, skew_y, scale_y];
+ let (snap_x, snap_y) = snap_origin(grid_x, grid_y, ¤t);
+ if (snap_x - upper_left_x).abs() > SNAP_EPS {
+ match plan.target {
+ TargetGrid::Dims { .. } => scale_x = (env.max_x - snap_x).abs() / width as f64,
+ _ => {
+ width = ((env.max_x - snap_x).abs() / scale_x.abs()).ceil().max(1.0) as i64;
+ }
+ }
+ upper_left_x = snap_x;
+ }
+ if (snap_y - upper_left_y).abs() > SNAP_EPS {
+ match plan.target {
+ TargetGrid::Dims { .. } => {
+ scale_y = scale_y.signum() * (env.min_y - snap_y).abs() / height as f64;
+ }
+ _ => {
+ height = ((env.min_y - snap_y).abs() / scale_y.abs()).ceil().max(1.0) as i64;
+ }
+ }
+ upper_left_y = snap_y;
+ }
+ }
+
+ Grid {
+ transform: [upper_left_x, scale_x, skew_x, upper_left_y, skew_y, scale_y],
+ width,
+ height,
+ }
+}
+
+/// Snap the origin of `transform` (its `[0]`/`[3]` upper-left) to the grid
+/// anchored at `(grid_x, grid_y)` with `transform`'s pixel size and skew: the
+/// returned point is the grid-cell corner that contains the origin
+/// (inverse-affine then floor, matching Spark's `getGridCoordinatesFromWorld` +
+/// `getWorldCornerCoordinates`).
+///
+/// The math is the shared affine machinery: build the grid's own geotransform
+/// (the source pixel size/skew anchored at `(grid_x, grid_y)`), invert it to
+/// find the origin's fractional cell, floor to the containing corner, and map
+/// that corner back. A degenerate (uninvertible) transform leaves the origin
+/// unmoved.
+fn snap_origin(grid_x: f64, grid_y: f64, transform: &GeoTransform) -> (f64, f64) {
+ let [ul_x, scale_x, skew_x, ul_y, skew_y, scale_y] = *transform;
+ let anchored: GeoTransform = [grid_x, scale_x, skew_x, grid_y, skew_y, scale_y];
+ let Ok(inverse) = anchored.invert() else {
+ return (ul_x, ul_y);
+ };
+ let (col, row) = inverse.apply(ul_x, ul_y);
+ anchored.apply(col.floor(), row.floor())
+}
+
+#[cfg(test)]
+mod tests {
+ use super::*;
+ use datafusion_common::ScalarValue;
+ use sedona_raster::array::RasterStructArray;
+ use sedona_testing::raster_spec::{
+ assert_raster_scalar_equals, assert_rasters_equal, raster_array, RasterSpec,
+ };
+ use sedona_testing::testers::ScalarUdfTester;
+
+ /// Tester for the 5-argument `(raster, widthOrScale, heightOrScale, useScale,
+ /// algorithm)` overload.
+ fn tester5() -> ScalarUdfTester {
+ ScalarUdfTester::new(
+ rs_resample_udf().into(),
+ vec![
+ RASTER,
+ SedonaType::Arrow(DataType::Float64),
+ SedonaType::Arrow(DataType::Float64),
+ SedonaType::Arrow(DataType::Boolean),
+ SedonaType::Arrow(DataType::Utf8),
+ ],
+ )
+ }
+
+ /// Tester for the 7-argument `(..., gridX, gridY, useScale, algorithm)`
+ /// overload.
+ fn tester7() -> ScalarUdfTester {
+ ScalarUdfTester::new(
+ rs_resample_udf().into(),
+ vec![
+ RASTER,
+ SedonaType::Arrow(DataType::Float64),
+ SedonaType::Arrow(DataType::Float64),
+ SedonaType::Arrow(DataType::Float64),
+ SedonaType::Arrow(DataType::Float64),
+ SedonaType::Arrow(DataType::Boolean),
+ SedonaType::Arrow(DataType::Utf8),
+ ],
+ )
+ }
+
+ /// Tester for the 4-argument `(raster, referenceRaster, useScale, algorithm)`
+ /// overload.
+ fn tester_ref() -> ScalarUdfTester {
+ ScalarUdfTester::new(
+ rs_resample_udf().into(),
+ vec![
+ RASTER,
+ RASTER,
+ SedonaType::Arrow(DataType::Boolean),
+ SedonaType::Arrow(DataType::Utf8),
+ ],
+ )
+ }
+
+ fn raster_scalar(spec: &RasterSpec) -> ColumnarValue {
+ ColumnarValue::Scalar(ScalarValue::Struct(Arc::new(spec.build())))
+ }
+
+ fn f64_scalar(v: f64) -> ColumnarValue {
+ ColumnarValue::Scalar(ScalarValue::Float64(Some(v)))
+ }
+
+ fn bool_scalar(v: bool) -> ColumnarValue {
+ ColumnarValue::Scalar(ScalarValue::Boolean(Some(v)))
+ }
+
+ fn str_scalar(v: &str) -> ColumnarValue {
+ ColumnarValue::Scalar(ScalarValue::Utf8(Some(v.to_string())))
+ }
+
+ fn unwrap_scalar(result: Result<ColumnarValue>) -> ScalarValue {
+ match result.unwrap() {
+ ColumnarValue::Scalar(sv) => sv,
+ other => panic!("expected a scalar result, got {other:?}"),
+ }
+ }
+
+ /// Invoke the 5-argument dimension form `(width, height)` with nearest.
+ fn resample_dims(source: &RasterSpec, width: f64, height: f64) -> ScalarValue {
+ unwrap_scalar(tester5().invoke(vec![
+ raster_scalar(source),
+ f64_scalar(width),
+ f64_scalar(height),
+ bool_scalar(false),
+ str_scalar("NearestNeighbor"),
+ ]))
+ }
+
+ /// Invoke the 5-argument scale form `(scaleX, scaleY)` with nearest.
+ fn resample_scale(source: &RasterSpec, scale_x: f64, scale_y: f64) -> ScalarValue {
+ unwrap_scalar(tester5().invoke(vec![
+ raster_scalar(source),
+ f64_scalar(scale_x),
+ f64_scalar(scale_y),
+ bool_scalar(true),
+ str_scalar("NearestNeighbor"),
+ ]))
+ }
+
+ #[test]
+ fn resample_to_same_grid_is_identity() {
+ // Resampling to the source dimensions is a straight copy: pixels,
+ // transform, CRS, and nodata are all preserved exactly.
+ let source = RasterSpec::d2(2, 1)
+ .band_values(&[10u8, 20])
+ .nodata(255u8)
+ .bbox(0.0, 0.0, 4.0, 2.0)
+ .crs(Some("EPSG:4326"));
+
+ let result = resample_dims(&source, 2.0, 1.0);
+ assert_raster_scalar_equals(&result, &source);
+ }
+
+ #[test]
+ fn nearest_upsample_replicates_pixels() {
+ // A 2x1 -> 4x2 nearest upsample by an integer factor replicates each
+ // source pixel into a 2x2 block (block replication is unambiguous, so
+ // this is bit-exact and stable across GDAL versions). The world extent
+ // x[0,4], y[0,2] is preserved, halving the pixel size.
+ let source = RasterSpec::d2(2, 1)
+ .band_values(&[10u8, 20])
+ .bbox(0.0, 0.0, 4.0, 2.0)
+ .crs(Some("EPSG:4326"));
+
+ let result = resample_dims(&source, 4.0, 2.0);
+
+ let expected = RasterSpec::d2(4, 2)
+ .band_values(&[10u8, 10, 20, 20, 10, 10, 20, 20])
+ .bbox(0.0, 0.0, 4.0, 2.0)
+ .crs(Some("EPSG:4326"));
+ assert_raster_scalar_equals(&result, &expected);
+ }
+
+ #[test]
+ fn nearest_downsample_block_constant_is_stable() {
+ // Nearest downsampling picks source pixels; which column of a 2:1
+ // decimation it selects is GDAL-version-dependent, so the source here
+ // is block-constant per output cell (cols 0,1 -> 10; cols 2,3 -> 20;
+ // both rows identical) making the selected value unambiguous: [10, 20].
+ let source = RasterSpec::d2(4, 2)
+ .band_values(&[10u8, 10, 20, 20, 10, 10, 20, 20])
+ .bbox(0.0, 0.0, 4.0, 2.0)
+ .crs(Some("EPSG:4326"));
+
+ let result = resample_dims(&source, 2.0, 1.0);
+
+ let expected = RasterSpec::d2(2, 1)
+ .band_values(&[10u8, 20])
+ .bbox(0.0, 0.0, 4.0, 2.0)
+ .crs(Some("EPSG:4326"));
+ assert_raster_scalar_equals(&result, &expected);
+ }
+
+ #[test]
+ fn scale_mode_derives_dimensions() {
+ // Source pixel size is 2 (extent 4 over 2 columns, extent 2 over 1 row);
+ // requesting pixel size 1 gives dims ceil(4/1)=4 x ceil(2/1)=2. The extent
+ // divides evenly, so the extent does not grow: the output pixel size is
+ // exactly the requested 1 / -1 and the bbox is unchanged.
+ let source = RasterSpec::d2(2, 1)
+ .band_values(&[10u8, 20])
+ .bbox(0.0, 0.0, 4.0, 2.0)
+ .crs(Some("EPSG:4326"));
+
+ let result = resample_scale(&source, 1.0, -1.0);
+
+ let expected = RasterSpec::d2(4, 2)
+ .band_values(&[10u8, 10, 20, 20, 10, 10, 20, 20])
+ .bbox(0.0, 0.0, 4.0, 2.0)
+ .crs(Some("EPSG:4326"));
+ assert_raster_scalar_equals(&result, &expected);
+ }
+
+ #[test]
+ fn scale_mode_grows_extent_to_whole_pixels() {
+ // Extent 4 wide with pixel size 3 does not tile evenly: Spark keeps the
+ // scale exact (3) and grows to ceil(4/3)=2 columns, so the output spans
+ // x[0,6] (2*3) — one grown column past the source right edge (x=4) — and
+ // that grown column, uncovered by the source, is nodata. Height: extent 2
+ // with scale -3 gives ceil(2/3)=1 row spanning y[-1,2].
+ //
+ // Nearest maps each output column center to a source column: col 0 center
+ // x=1.5 -> source col 0 (10); col 1 center x=4.5 is past the source, so it
+ // falls to the nodata fill. Row 0 center y=0.5 -> source row 0.
+ let source = RasterSpec::d2(2, 1)
+ .band_values(&[10u8, 20])
+ .nodata(255u8)
+ .bbox(0.0, 0.0, 4.0, 2.0)
+ .crs(Some("EPSG:4326"));
+
+ let result = resample_scale(&source, 3.0, -3.0);
+
+ // Origin unchanged (0, 2); output 2x1 at pixel size 3/-3 -> bbox x[0,6],
+ // y[-1,2]. The source only covers x[0,4], so the second column is nodata.
+ let expected = RasterSpec::d2(2, 1)
+ .band_values(&[10u8, 255])
+ .nodata(255u8)
+ .bbox(0.0, -1.0, 6.0, 2.0)
+ .crs(Some("EPSG:4326"));
+ assert_raster_scalar_equals(&result, &expected);
+ }
+
+ #[test]
+ fn scale_mode_grows_extent_without_crs() {
+ // The regrid path must also work on a CRS-less raster (a same-pixel-space
+ // regrid) — the Python parity fixtures write CRS-less GeoTIFFs. Same grow
+ // arithmetic as the CRS'd case above.
+ let source = RasterSpec::d2(2, 1)
+ .band_values(&[10u8, 20])
+ .nodata(255u8)
+ .bbox(0.0, 0.0, 4.0, 2.0)
+ .crs(None);
+
+ let result = resample_scale(&source, 3.0, -3.0);
+
+ let expected = RasterSpec::d2(2, 1)
+ .band_values(&[10u8, 255])
+ .nodata(255u8)
+ .bbox(0.0, -1.0, 6.0, 2.0)
+ .crs(None);
+ assert_raster_scalar_equals(&result, &expected);
+ }
+
+ #[test]
+ fn grid_snap_shifts_origin_and_regrows() {
+ // Scale-mode pixel size 3/-3 snapped to a grid anchored at (1, 3): the
+ // source origin (0, 2) snaps outward to the grid corner containing it —
+ // x: floor((0-1)/3) = -1 -> 1 + (-1)*3 = -2; y: floor((2-3)/-3) = 0 ->
+ // 3 + 0*(-3) = 3. New origin (-2, 3) (a corner of the (1,3)/3 grid), and
+ // the source extent stays covered. Width regrows to reach the source
+ // right edge x=4: ceil(|4-(-2)|/3) = 2; height ceil(|0-3|/3) = 1.
+ let source = RasterSpec::d2(2, 1)
+ .band_values(&[10u8, 20])
+ .nodata(255u8)
+ .bbox(0.0, 0.0, 4.0, 2.0)
+ .crs(Some("EPSG:4326"));
+
+ // 7-arg overload: (raster, scaleX=3, scaleY=-3, gridX=1, gridY=3,
+ // useScale=true, algorithm).
+ let result = unwrap_scalar(tester7().invoke(vec![
+ raster_scalar(&source),
+ f64_scalar(3.0),
+ f64_scalar(-3.0),
+ f64_scalar(1.0),
+ f64_scalar(3.0),
+ bool_scalar(true),
+ str_scalar("NearestNeighbor"),
+ ]));
+
+ // Origin (-2, 3), 2x1 at 3/-3 -> bbox x[-2,4], y[0,3]. Cell centers are
+ // (-0.5, 1.5) — left of the source (x<0), nodata — and (2.5, 1.5) —
+ // inside source column 1, value 20.
+ let expected = RasterSpec::d2(2, 1)
+ .band_values(&[255u8, 20])
+ .nodata(255u8)
+ .bbox(-2.0, 0.0, 4.0, 3.0)
+ .crs(Some("EPSG:4326"));
+ assert_raster_scalar_equals(&result, &expected);
+ }
+
+ #[test]
+ fn reference_raster_matches_its_grid() {
+ // The 4-arg reference overload with useScale=false takes the reference's
+ // dimensions and origin. Here the reference is a 4x2 grid over the same
+ // extent as the 2x1 source, so the result is the 2x upsample.
+ let source = RasterSpec::d2(2, 1)
+ .band_values(&[10u8, 20])
+ .bbox(0.0, 0.0, 4.0, 2.0)
+ .crs(Some("EPSG:4326"));
+ let reference = RasterSpec::d2(4, 2)
+ .band_values(&[0u8; 8])
+ .bbox(0.0, 0.0, 4.0, 2.0)
+ .crs(Some("EPSG:4326"));
+
+ let result = unwrap_scalar(tester_ref().invoke(vec![
+ raster_scalar(&source),
+ raster_scalar(&reference),
+ bool_scalar(false),
+ str_scalar("NearestNeighbor"),
+ ]));
+
+ let expected = RasterSpec::d2(4, 2)
+ .band_values(&[10u8, 10, 20, 20, 10, 10, 20, 20])
+ .bbox(0.0, 0.0, 4.0, 2.0)
+ .crs(Some("EPSG:4326"));
+ assert_raster_scalar_equals(&result, &expected);
+ }
+
+ #[test]
+ fn reference_raster_crs_mismatch_errors() {
+ // A reference raster in a different CRS must error — RS_Resample never
+ // reprojects.
+ let source = RasterSpec::d2(2, 1)
+ .band_values(&[10u8, 20])
+ .bbox(0.0, 0.0, 4.0, 2.0)
+ .crs(Some("EPSG:4326"));
+ let reference = RasterSpec::d2(4, 2)
+ .band_values(&[0u8; 8])
+ .bbox(0.0, 0.0, 4.0, 2.0)
+ .crs(Some("EPSG:3857"));
+
+ let err = tester_ref()
+ .invoke(vec![
+ raster_scalar(&source),
+ raster_scalar(&reference),
+ bool_scalar(false),
+ str_scalar("NearestNeighbor"),
+ ])
+ .unwrap_err()
+ .to_string();
+ assert!(err.contains("does not reproject"), "got: {err}");
+ }
+
+ #[test]
+ fn multiband_nearest_upsample() {
+ // Every band is resampled and the band order/count preserved.
+ let source = RasterSpec::d2(2, 1)
+ .band_values(&[10u8, 20])
+ .nodata(0u8)
+ .band_values(&[1u8, 2])
+ .bbox(0.0, 0.0, 4.0, 2.0)
+ .crs(Some("EPSG:4326"));
+
+ let result = resample_dims(&source, 4.0, 2.0);
+
+ let expected = RasterSpec::d2(4, 2)
+ .band_values(&[10u8, 10, 20, 20, 10, 10, 20, 20])
+ .nodata(0u8)
+ .band_values(&[1u8, 1, 2, 2, 1, 1, 2, 2])
+ .bbox(0.0, 0.0, 4.0, 2.0)
+ .crs(Some("EPSG:4326"));
+ assert_raster_scalar_equals(&result, &expected);
+ }
+
+ #[test]
+ fn nd_band_broadcasts_across_planes() {
+ // A [time=2, y=1, x=2] band upsampled to (y=2, x=4): the time dimension
+ // is preserved and each time plane is resampled independently. Time 0 is
+ // [10, 20], time 1 is [30, 40].
+ let source = RasterSpec::nd(&["time", "y", "x"], &[2, 1, 2])
+ .band_values(&[10u8, 20, 30, 40])
+ .transform([0.0, 2.0, 0.0, 2.0, 0.0, -2.0])
+ .crs(Some("EPSG:4326"));
+
+ let result = resample_dims(&source, 4.0, 2.0);
+
+ let expected = RasterSpec::nd(&["time", "y", "x"], &[2, 2, 4])
+ .band_values(&[
+ 10u8, 10, 20, 20, 10, 10, 20, 20, // time 0
+ 30, 30, 40, 40, 30, 30, 40, 40, // time 1
+ ])
+ .transform([0.0, 1.0, 0.0, 2.0, 0.0, -1.0])
+ .crs(Some("EPSG:4326"));
+ assert_raster_scalar_equals(&result, &expected);
+ }
+
+ #[test]
+ fn null_raster_yields_null() {
+ let result = tester5()
+ .invoke(vec![
+ ColumnarValue::Array(Arc::new(raster_array(vec![
+ Some(
+ RasterSpec::d2(2, 1)
+ .band_values(&[10u8, 20])
+ .bbox(0.0, 0.0, 4.0, 2.0)
+ .crs(Some("EPSG:4326")),
+ ),
+ None,
+ ]))),
+ f64_scalar(4.0),
+ f64_scalar(2.0),
+ bool_scalar(false),
+ str_scalar("NearestNeighbor"),
+ ])
+ .unwrap();
+ let ColumnarValue::Array(arr) = result else {
+ panic!("expected an array result");
+ };
+
+ let expected = vec![
+ Some(
+ RasterSpec::d2(4, 2)
+ .band_values(&[10u8, 10, 20, 20, 10, 10, 20, 20])
+ .bbox(0.0, 0.0, 4.0, 2.0)
+ .crs(Some("EPSG:4326")),
+ ),
+ None,
+ ];
+ assert_rasters_equal(&arr, &expected);
+ }
+
+ #[test]
+ fn skewed_raster_keeps_skew_and_derives_scale_from_envelope() {
+ // A skewed (non-north-up) raster resamples with the shear carried through
+ // UNCHANGED and the new pixel size derived from the axis-aligned envelope
+ // divided across the new dimensions, matching Sedona Spark (which keeps
+ // the original skew and derives the scale from the envelope, ignoring the
+ // skew for the scale). Source transform [ulx=0, scale_x=2, skew_x=0.5,
+ // uly=2, skew_y=0.5, scale_y=-2] over 2x1 has envelope x in [0, 4.5],
+ // y in [0, 3], so upsampling to 4x2 gives scale_x = 4.5/4 = 1.125,
+ // scale_y = -3/2 = -1.5, with skew unchanged at 0.5. The pixels are the
+ // extent-preserving RasterIO decimation.
+ let source = RasterSpec::d2(2, 1)
+ .band_values(&[10u8, 20])
+ .transform([0.0, 2.0, 0.5, 2.0, 0.5, -2.0])
+ .crs(Some("EPSG:4326"));
+
+ let result = resample_dims(&source, 4.0, 2.0);
+
+ let expected = RasterSpec::d2(4, 2)
+ .band_values(&[10u8, 10, 20, 20, 10, 10, 20, 20])
+ .transform([0.0, 1.125, 0.5, 2.0, 0.5, -1.5])
+ .crs(Some("EPSG:4326"));
+ assert_raster_scalar_equals(&result, &expected);
+ }
+
+ #[test]
+ fn skewed_scale_change_errors() {
+ // A skewed raster in the grow/snap regrid path (scale mode) is an
+ // explicit unsupported case and errors, rather than silently producing a
+ // wrong grid.
+ let source = RasterSpec::d2(2, 1)
+ .band_values(&[10u8, 20])
+ .transform([0.0, 2.0, 0.5, 2.0, 0.5, -2.0])
+ .crs(Some("EPSG:4326"));
+
+ let err = tester5()
+ .invoke(vec![
+ raster_scalar(&source),
+ f64_scalar(3.0),
+ f64_scalar(-3.0),
+ bool_scalar(true),
+ str_scalar("NearestNeighbor"),
+ ])
+ .unwrap_err()
+ .to_string();
+ assert!(err.contains("skewed"), "got: {err}");
+ }
+
+ #[test]
+ fn int64_allowed_only_on_the_nearest_fast_path() {
+ // 64-bit ints are exact only on the nearest RasterIO fast path (a plain
+ // dimension change). Any warp-bound case — an interpolating algorithm, or
+ // a scale change / origin snap / reference grid even with nearest —
+ // routes pixels through a floating working type that cannot represent
+ // Int64/UInt64 exactly, so RS_Resample rejects it (matching
+ // RS_ReprojectMatch).
+ let source = RasterSpec::d2(2, 1)
+ .band_values(&[10i64, 20])
+ .bbox(0.0, 0.0, 4.0, 2.0)
+ .crs(Some("EPSG:4326"));
+
+ // Non-nearest dimension change: warp-bound, rejected.
+ let err = tester5()
+ .invoke(vec![
+ raster_scalar(&source),
+ f64_scalar(4.0),
+ f64_scalar(2.0),
+ bool_scalar(false),
+ str_scalar("Bilinear"),
+ ])
+ .unwrap_err()
+ .to_string();
+ assert!(err.contains("Int64/UInt64"), "got: {err}");
+
+ // Nearest *scale* change also warps (grow-extent / origin snap), so it is
+ // rejected too despite the nearest algorithm.
+ let err = tester5()
+ .invoke(vec![
+ raster_scalar(&source),
+ f64_scalar(3.0),
+ f64_scalar(-3.0),
+ bool_scalar(true),
+ str_scalar("NearestNeighbor"),
+ ])
+ .unwrap_err()
+ .to_string();
+ assert!(err.contains("Int64/UInt64"), "got: {err}");
+
+ // Nearest dimension change is the bit-exact fast path and succeeds.
+ tester5()
+ .invoke(vec![
+ raster_scalar(&source),
+ f64_scalar(4.0),
+ f64_scalar(2.0),
+ bool_scalar(false),
+ str_scalar("NearestNeighbor"),
+ ])
+ .unwrap();
+ }
+
+ #[test]
+ fn nearest_fast_path_preserves_int64_exactly() {
+ // The premise for allowing Int64/UInt64 only on the nearest fast path is
+ // that RasterIO nearest copies a source sample rather than computing it, so
+ // it survives values that f64 cannot represent (2^53+1 rounds to 2^53, and
+ // 2^53+3 to 2^53+4, if routed through a floating working type). Verified for
+ // both up- and down-sampling.
+ let a: i64 = (1i64 << 53) + 1;
+ let b: i64 = (1i64 << 53) + 3;
+
+ // 2x1 -> 4x2 nearest upsample replicates each pixel.
+ let up_src = RasterSpec::d2(2, 1)
+ .band_values(&[a, b])
+ .bbox(0.0, 0.0, 4.0, 2.0)
+ .crs(Some("EPSG:4326"));
+ let up_expected = RasterSpec::d2(4, 2)
+ .band_values(&[a, a, b, b, a, a, b, b])
+ .bbox(0.0, 0.0, 4.0, 2.0)
+ .crs(Some("EPSG:4326"));
+ assert_raster_scalar_equals(&resample_dims(&up_src, 4.0, 2.0), &up_expected);
+
+ // 4x1 -> 2x1 nearest decimation (block-constant so the selected column is
+ // unambiguous, as in nearest_downsample_block_constant_is_stable).
+ let down_src = RasterSpec::d2(4, 1)
+ .band_values(&[a, a, b, b])
+ .bbox(0.0, 0.0, 4.0, 1.0)
+ .crs(Some("EPSG:4326"));
+ let down_expected = RasterSpec::d2(2, 1)
+ .band_values(&[a, b])
+ .bbox(0.0, 0.0, 4.0, 1.0)
+ .crs(Some("EPSG:4326"));
+ assert_raster_scalar_equals(&resample_dims(&down_src, 2.0, 1.0), &down_expected);
+ }
+
+ #[test]
+ fn nearest_fast_path_preserves_uint64_exactly() {
+ // UInt64 above i64::MAX and above 2^63 (neither representable in f64) also
+ // survives the nearest fast path unchanged.
+ let a: u64 = u64::MAX - 1;
+ let b: u64 = (1u64 << 63) + 5;
+ let source = RasterSpec::d2(2, 1)
+ .band_values(&[a, b])
+ .bbox(0.0, 0.0, 4.0, 2.0)
+ .crs(Some("EPSG:4326"));
+ let expected = RasterSpec::d2(4, 2)
+ .band_values(&[a, a, b, b, a, a, b, b])
+ .bbox(0.0, 0.0, 4.0, 2.0)
+ .crs(Some("EPSG:4326"));
+ assert_raster_scalar_equals(&resample_dims(&source, 4.0, 2.0), &expected);
+ }
+
+ #[test]
+ fn int32_mode_gated_on_gdal_version() {
+ // Mode selects the most-common source value, so it should return an exact
+ // source value — but GDAL routes Mode through a 32-bit float working type
+ // on GDAL < 3.13 (a double from 3.13 on), and a 32-bit float cannot
+ // represent Int32 above 2^24 (2^24+1 -> 2^24, 2^24+3 -> 2^24+4). So on older
+ // GDAL RS_Resample rejects Int32/UInt32 + Mode rather than silently
+ // corrupting a category code; on GDAL >= 3.13 the value survives exactly.
+ let a: i32 = (1i32 << 24) + 1;
+ let b: i32 = (1i32 << 24) + 3;
+ // 4x1 -> 2x1 block-constant windows so Mode's selection is unambiguous.
+ let source = RasterSpec::d2(4, 1)
+ .band_values(&[a, a, b, b])
+ .bbox(0.0, 0.0, 4.0, 1.0)
+ .crs(Some("EPSG:4326"));
+ let invoke = || {
+ tester5().invoke(vec![
+ raster_scalar(&source),
+ f64_scalar(2.0),
+ f64_scalar(1.0),
+ bool_scalar(false),
+ str_scalar("Mode"),
+ ])
+ };
+
+ let gdal_version = with_gdal(|gdal| Ok(gdal.version_num())).unwrap();
+ if gdal_version >= 3_130_000 {
+ let expected = RasterSpec::d2(2, 1)
+ .band_values(&[a, b])
+ .bbox(0.0, 0.0, 4.0, 1.0)
+ .crs(Some("EPSG:4326"));
+ assert_raster_scalar_equals(&unwrap_scalar(invoke()), &expected);
+ } else {
+ let err = invoke().unwrap_err().to_string();
+ assert!(err.contains("Int32/UInt32"), "got: {err}");
+ assert!(err.contains("Mode"), "got: {err}");
+ }
+ }
+
+ #[test]
+ fn bilinear_scale_change_of_constant_is_constant_via_warp() {
+ // A non-nearest algorithm on the regrid (warp) path — here a scale change.
+ // Interpolating a constant band yields that constant everywhere — exact
+ // and GDAL-version-independent — which pins the warp path end-to-end (the
+ // blended values near edges are pinned against rasterio in the Python
+ // parity tests). Pixel size 1 over the 4x4 extent divides evenly (no grown
+ // border), so the whole 4x4 output is the constant.
+ let source = RasterSpec::d2(2, 2)
+ .band_values(&[7.0f64, 7.0, 7.0, 7.0])
+ .bbox(0.0, 0.0, 4.0, 4.0)
+ .crs(Some("EPSG:4326"));
+
+ // 5-arg overload: (raster, scaleX=1, scaleY=-1, useScale=true, algorithm).
+ let result = unwrap_scalar(tester5().invoke(vec![
+ raster_scalar(&source),
+ f64_scalar(1.0),
+ f64_scalar(-1.0),
+ bool_scalar(true),
+ str_scalar("Bilinear"),
+ ]));
+
+ let expected = RasterSpec::d2(4, 4)
+ .band_values(&[7.0f64; 16])
+ .bbox(0.0, 0.0, 4.0, 4.0)
+ .crs(Some("EPSG:4326"));
+ assert_raster_scalar_equals(&result, &expected);
+ }
+
+ #[test]
+ fn fractional_dimension_errors() {
+ // useScale=false with a fractional width is rejected rather than
+ // truncated.
+ let source = RasterSpec::d2(2, 1)
+ .band_values(&[10u8, 20])
+ .bbox(0.0, 0.0, 4.0, 2.0)
+ .crs(Some("EPSG:4326"));
+
+ let err = tester5()
+ .invoke(vec![
+ raster_scalar(&source),
+ f64_scalar(3.5),
+ f64_scalar(2.0),
+ bool_scalar(false),
+ str_scalar("NearestNeighbor"),
+ ])
+ .unwrap_err()
+ .to_string();
+ assert!(err.contains("whole number"), "got: {err}");
+ }
+
+ #[test]
+ fn unknown_algorithm_errors_through_udf() {
+ let source = RasterSpec::d2(2, 1)
+ .band_values(&[10u8, 20])
+ .bbox(0.0, 0.0, 4.0, 2.0)
+ .crs(Some("EPSG:4326"));
+
+ let err = tester5()
+ .invoke(vec![
+ raster_scalar(&source),
+ f64_scalar(4.0),
+ f64_scalar(2.0),
+ bool_scalar(false),
+ str_scalar("sinc"),
+ ])
+ .unwrap_err()
+ .to_string();
+ assert!(err.contains("unknown algorithm"), "got: {err}");
+ }
+
+ #[test]
+ fn output_keeps_input_crs_and_origin() {
+ // Structural pin: the output CRS always equals the input CRS (RS_Resample
+ // never reprojects). A 4x4 EPSG:4326 raster resampled to 8x8 keeps its
+ // CRS and origin.
+ let source = RasterSpec::d2(4, 4)
+ .band_values(&(0u8..16).collect::<Vec<_>>())
+ .bbox(10.0, 40.0, 14.0, 44.0)
+ .crs(Some("EPSG:4326"));
+
+ let result = resample_dims(&source, 8.0, 8.0);
+
+ let ScalarValue::Struct(struct_array) = &result else {
+ panic!("expected a raster struct scalar, got {result:?}");
+ };
+ let rasters = RasterStructArray::try_new(struct_array).unwrap();
+ let raster = rasters.get(0).unwrap();
+ assert_eq!(raster.crs(), Some("EPSG:4326"));
+ let m = raster.metadata();
+ assert_eq!(m.upper_left_x(), 10.0);
+ assert_eq!(m.upper_left_y(), 44.0);
+ }
+}
diff --git a/rust/sedona-raster-gdal/src/utils.rs b/rust/sedona-raster-gdal/src/utils.rs
index 3bd5e41a9..f0f7cc451 100644
--- a/rust/sedona-raster-gdal/src/utils.rs
+++ b/rust/sedona-raster-gdal/src/utils.rs
@@ -20,19 +20,20 @@
use arrow_array::StructArray;
use arrow_buffer::Buffer;
use datafusion_common::error::Result;
-use datafusion_common::exec_datafusion_err;
+use datafusion_common::{exec_datafusion_err, exec_err};
use sedona_gdal::dataset::Dataset;
use sedona_gdal::gdal::Gdal;
use sedona_gdal::gdal_dyn_bindgen::{GDAL_OF_RASTER, GDAL_OF_READONLY};
-use sedona_gdal::geo_transform::GeoTransform;
+use sedona_gdal::geo_transform::{GeoTransform, GeoTransformEx};
use sedona_gdal::mem::MemDatasetBuilder;
use sedona_gdal::raster::types::DatasetOptions;
use sedona_gdal::raster::types::ResampleAlg;
use sedona_gdal::spatial_ref::SpatialRef;
+use sedona_raster::array::RasterRefImpl;
use sedona_raster::builder::RasterBuilder;
-use sedona_raster::traits::{BandMetadata, RasterMetadata};
-use sedona_schema::raster::StorageType;
+use sedona_raster::traits::{BandMetadata, RasterMetadata, RasterRef};
+use sedona_schema::raster::{BandDataType, StorageType};
use crate::gdal_common::{
band_data_type_to_gdal, band_nodata_to_bytes, convert_gdal_err, gdal_to_band_data_type,
@@ -208,23 +209,195 @@ pub fn append_nd_from_dataset(
append_nd_from_dataset_inner(dataset, layout, builder, &metadata, crs.as_deref(), None)
}
-/// The output grid a warp writes into: geotransform, pixel dimensions, target
-/// CRS, and the resampling algorithm.
+/// A raster's affine grid: its geotransform and pixel dimensions.
///
-/// The warp path can move the output origin off the source grid, grow the
-/// extent past the source footprint, and change the CRS — so `crs` here is the
-/// *destination* CRS (equal to the source CRS for a same-CRS regrid, or the
-/// reprojection target). `crs` is carried through rather than read back from
-/// GDAL so a caller can preserve the exact CRS string it was given.
-pub struct WarpGrid<'a> {
+/// This is the geometric core shared by a source raster's grid (via
+/// [`Grid::from_metadata`]) and every resample/warp output target, replacing the
+/// several near-identical transform-plus-dimensions structs that had accumulated
+/// here.
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub struct Grid {
pub transform: GeoTransform,
- pub width: usize,
- pub height: usize,
+ pub width: i64,
+ pub height: i64,
+}
+
+/// Axis-aligned world bounding box of a [`Grid`]'s four corners.
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub struct Envelope {
+ pub min_x: f64,
+ pub max_x: f64,
+ pub min_y: f64,
+ pub max_y: f64,
+}
+
+impl Grid {
+ /// The grid of an existing raster, read from its [`RasterMetadata`].
+ pub fn from_metadata(m: &RasterMetadata) -> Self {
+ Self {
+ transform: [
+ m.upper_left_x(),
+ m.scale_x(),
+ m.skew_x(),
+ m.upper_left_y(),
+ m.skew_y(),
+ m.scale_y(),
+ ],
+ width: m.width(),
+ height: m.height(),
+ }
+ }
+
+ /// Bounding box of the four grid corners (handles skew; reduces to
+ /// `width * |scale_x|` etc. for a north-up grid). Corners are mapped with the
+ /// shared [`GeoTransformEx::apply`] — the same affine used everywhere else.
+ pub fn envelope(&self) -> Envelope {
+ let w = self.width as f64;
+ let h = self.height as f64;
+ let mut min_x = f64::INFINITY;
+ let mut max_x = f64::NEG_INFINITY;
+ let mut min_y = f64::INFINITY;
+ let mut max_y = f64::NEG_INFINITY;
+ for (col, row) in [(0.0, 0.0), (w, 0.0), (0.0, h), (w, h)] {
+ let (x, y) = self.transform.apply(col, row);
+ min_x = min_x.min(x);
+ max_x = max_x.max(x);
+ min_y = min_y.min(y);
+ max_y = max_y.max(y);
+ }
+ Envelope {
+ min_x,
+ max_x,
+ min_y,
+ max_y,
+ }
+ }
+}
+
+/// The output grid a resample or warp writes into: a [`Grid`] plus the CRS to
+/// record on the output and the resampling algorithm.
+///
+/// `crs` is the CRS stamped on the output — the source CRS for a same-CRS
+/// resample or regrid, the reprojection target for a warp that changes CRS. It
+/// is carried through verbatim rather than read back from GDAL so a caller can
+/// preserve the exact CRS string it was given.
+pub struct OutputGrid<'a> {
+ pub grid: Grid,
pub crs: Option<&'a str>,
pub alg: ResampleAlg,
- /// GDAL's working-memory cache size (bytes) for the warp; `0.0` uses
- /// GDAL's own default.
- pub warp_memory_limit_bytes: f64,
+}
+
+/// GDAL 3.13 is the first release that resamples Mode through a double working
+/// type instead of a 32-bit float.
+const GDAL_3_13_0: i32 = 3_130_000;
+
+/// Map an algorithm name (case-insensitive) to a GDAL [`ResampleAlg`]. An empty
+/// string defaults to nearest neighbour (matching Sedona Spark). `func_name`
+/// names the calling SQL function for the error message.
+///
+/// Accepts the GDAL names plus the Spark spellings (American `NearestNeighbor`
+/// and `Bicubic`, which GDAL calls `Cubic`). Shared by the RasterIO-resampling
+/// and warping UDFs so they accept an identical algorithm surface.
+pub fn parse_resample_algorithm(name: &str, func_name: &str) -> Result<ResampleAlg> {
+ if name.is_empty() {
+ return Ok(ResampleAlg::NearestNeighbour);
+ }
+ let alg = match name.to_ascii_lowercase().as_str() {
+ "nearestneighbor" | "nearestneighbour" | "nearest" | "near" => {
+ ResampleAlg::NearestNeighbour
+ }
+ "bilinear" => ResampleAlg::Bilinear,
+ "cubic" | "bicubic" => ResampleAlg::Cubic,
+ "cubicspline" => ResampleAlg::CubicSpline,
+ "lanczos" => ResampleAlg::Lanczos,
+ "average" => ResampleAlg::Average,
+ "mode" => ResampleAlg::Mode,
+ _ => {
+ return exec_err!(
+ "{func_name}: unknown algorithm {name:?}; expected one of \
+ NearestNeighbor, Bilinear, Cubic, CubicSpline, Lanczos, Average, Mode"
+ );
+ }
+ };
+ Ok(alg)
+}
+
+/// Reject the integer dtypes GDAL's floating resample/warp working type cannot
+/// represent exactly.
+///
+/// `routes_through_float` is true when the pixels will pass through a floating
+/// working type: always for a warp/reproject, and for RasterIO resampling only
+/// on a regrid or an interpolating algorithm — a plain nearest-neighbour
+/// decimation is an exact native-type copy, so it passes `false` and every dtype
+/// is allowed. `func_name` names the calling SQL function for the error message.
+///
+/// `Int64`/`UInt64` are never representable in a floating working type (a double
+/// is exact only to 2^53). `Mode` additionally uses a 32-bit float working type
+/// on GDAL < 3.13, which cannot represent `Int32`/`UInt32` above 2^24; since Mode
+/// is a value selection (the most-common source value), silently rounding a
+/// category code is especially wrong, so those are rejected on older GDAL too.
+pub fn reject_lossy_resample_dtypes(
+ gdal: &Gdal,
+ raster: &RasterRefImpl<'_>,
+ band_count: usize,
+ alg: ResampleAlg,
+ routes_through_float: bool,
+ func_name: &str,
+) -> Result<()> {
+ if !routes_through_float {
+ return Ok(());
+ }
+ let mode_uses_float32 = alg == ResampleAlg::Mode && gdal.version_num() < GDAL_3_13_0;
+ for i in 0..band_count {
+ match raster.band_data_type(i) {
+ Some(BandDataType::Int64 | BandDataType::UInt64) => {
+ return exec_err!(
+ "{func_name} does not support Int64/UInt64 rasters: GDAL routes 64-bit \
+ integer pixels through a floating working type that cannot represent them \
+ exactly; cast to Int32/Float64 first."
+ );
+ }
+ Some(BandDataType::Int32 | BandDataType::UInt32) if mode_uses_float32 => {
+ return exec_err!(
+ "{func_name} does not support Mode resampling of Int32/UInt32 rasters on GDAL \
+ {} (before 3.13): Mode routes pixels through a 32-bit float working type that \
+ cannot represent 32-bit integers above 2^24 exactly. Upgrade to GDAL >= 3.13, \
+ which resamples Mode through a double, or cast to Float64 first.",
+ gdal.version_info("RELEASE_NAME")
+ );
+ }
+ _ => {}
+ }
+ }
+ Ok(())
+}
+
+/// Append a GDAL dataset as a single **N-D** in-db raster, resampling every band
+/// to `output`'s dimensions with `output.alg`.
+///
+/// The spatial analog of [`append_nd_from_dataset`]: the full source window of
+/// each GDAL band is read into an `output.grid.width` x `output.grid.height`
+/// buffer using GDAL's RasterIO resampling, so band count/order and the
+/// non-spatial structure in `layout` are preserved and only the trailing
+/// `(y, x)` extent changes.
+pub fn append_resampled_nd_from_dataset(
+ dataset: &Dataset,
+ layout: &GdalBandLayout,
+ builder: &mut RasterBuilder,
+ output: &OutputGrid<'_>,
+) -> Result<()> {
+ let metadata = output
+ .grid
+ .transform
+ .to_raster_metadata(output.grid.width as usize, output.grid.height as usize);
+ append_nd_from_dataset_inner(
+ dataset,
+ layout,
+ builder,
+ &metadata,
+ output.crs,
+ Some(output.alg),
+ )
}
/// Warp every band of `src_dataset` into `grid`, appending the result as a
@@ -236,17 +409,21 @@ pub struct WarpGrid<'a> {
/// source, or shift under a reprojection — read back as nodata rather than as an
/// accidental zero (`GDALReprojectImage` writes only covered pixels). GDAL warps
/// into those buffers; the warped bytes are then regrouped into N-D bands via
-/// [`append_nd_from_dataset_inner`], carrying `grid.crs` through verbatim rather
-/// than round-tripping it back out of GDAL.
+/// [`append_nd_from_dataset_inner`], carrying `output.crs` through verbatim
+/// rather than round-tripping it back out of GDAL.
+///
+/// `warp_memory_limit_bytes` is GDAL's working-memory cache size for the warp;
+/// pass `0.0` for GDAL's own default.
pub fn append_warped_nd_from_dataset(
gdal: &Gdal,
src_dataset: &Dataset,
layout: &GdalBandLayout,
builder: &mut RasterBuilder,
- grid: &WarpGrid<'_>,
+ output: &OutputGrid<'_>,
+ warp_memory_limit_bytes: f64,
) -> Result<()> {
- let out_width = grid.width;
- let out_height = grid.height;
+ let out_width = output.grid.width as usize;
+ let out_height = output.grid.height as usize;
// One owned, nodata-pre-filled buffer per source band, its planes
// concatenated plane-major; the destination's DATAPOINTER bands point into
@@ -278,9 +455,9 @@ pub fn append_warped_nd_from_dataset(
// SAFETY: the DATAPOINTER buffers in `band_buffers` outlive `dst_dataset`.
let dst_dataset = unsafe { dst_builder.build(gdal).map_err(convert_gdal_err)? };
dst_dataset
- .set_geo_transform(&grid.transform)
+ .set_geo_transform(&output.grid.transform)
.map_err(convert_gdal_err)?;
- if let Some(crs) = grid.crs {
+ if let Some(crs) = output.crs {
dst_dataset.set_projection(crs).map_err(convert_gdal_err)?;
}
@@ -306,15 +483,18 @@ pub fn append_warped_nd_from_dataset(
gdal.reproject_image(
src_dataset,
&dst_dataset,
- grid.alg,
- grid.warp_memory_limit_bytes,
+ output.alg,
+ warp_memory_limit_bytes,
)
.map_err(convert_gdal_err)?;
- let metadata = grid.transform.to_raster_metadata(out_width, out_height);
+ let metadata = output
+ .grid
+ .transform
+ .to_raster_metadata(out_width, out_height);
// Read the warped buffers back out (native read, no further resampling),
// regrouping the flat GDAL bands into N-D raster bands.
- append_nd_from_dataset_inner(&dst_dataset, layout, builder, &metadata, grid.crs, None)
+ append_nd_from_dataset_inner(&dst_dataset, layout, builder, &metadata, output.crs, None)
}
/// Allocate a `total`-byte buffer pre-filled with the little-endian `nodata`
@@ -383,9 +563,14 @@ fn append_nd_from_dataset_inner(
)
.map_err(|e| exec_datafusion_err!("Failed to start band: {}", e))?;
- let mut band_data: Vec<u8> = Vec::with_capacity(
- plan.plane_count * out_width * out_height * plan.data_type.byte_size(),
- );
+ let capacity = out_width
+ .checked_mul(out_height)
+ .and_then(|a| a.checked_mul(plan.data_type.byte_size()))
+ .and_then(|a| a.checked_mul(plan.plane_count))
+ .ok_or_else(|| {
+ exec_datafusion_err!("resampled band size overflow ({out_width}x{out_height})")
+ })?;
+ let mut band_data: Vec<u8> = Vec::with_capacity(capacity);
for _ in 0..plan.plane_count {
let band = dataset
.rasterband(gdal_band)