From de305db0a00135607c66dcdd8b190ec0514cbb82 Mon Sep 17 00:00:00 2001 From: jameswillis Date: Fri, 17 Jul 2026 01:05:43 -0700 Subject: [PATCH 1/5] [GH-3111] Fix scanline x-intercepts on rasters with non-square pixels --- .../sedona/common/raster/Rasterization.java | 6 +- .../common/raster/RasterConstructorsTest.java | 81 +++++++++++++++++++ 2 files changed, 86 insertions(+), 1 deletion(-) diff --git a/common/src/main/java/org/apache/sedona/common/raster/Rasterization.java b/common/src/main/java/org/apache/sedona/common/raster/Rasterization.java index 1449700b609..24b72b49e98 100644 --- a/common/src/main/java/org/apache/sedona/common/raster/Rasterization.java +++ b/common/src/main/java/org/apache/sedona/common/raster/Rasterization.java @@ -669,7 +669,11 @@ private static Map> computeScanlineIntersections( double xMax = (geomExtent.getMaxX() - params.upperLeftX) / params.scaleX; for (double y = yStart; y >= yEnd; y--) { - double xIntercept = p1X + ((p1Y - y) / slope); + // p1X, p1Y and y are in pixel units while slope is world-units dy + // over dx, so converting the pixel dy to world units (scaleY) and + // the resulting world dx back to pixels (scaleX) keeps the + // intercept in pixel space for any pixel aspect ratio. + double xIntercept = p1X + ((y - p1Y) * params.scaleY / slope / params.scaleX); if ((xIntercept < xMin) || (xIntercept >= xMax)) { continue; // Skip xIntercepts outside geomExtent } diff --git a/common/src/test/java/org/apache/sedona/common/raster/RasterConstructorsTest.java b/common/src/test/java/org/apache/sedona/common/raster/RasterConstructorsTest.java index 7bfb37c30dc..cb556f6a59c 100644 --- a/common/src/test/java/org/apache/sedona/common/raster/RasterConstructorsTest.java +++ b/common/src/test/java/org/apache/sedona/common/raster/RasterConstructorsTest.java @@ -388,6 +388,87 @@ public void testAsRasterWithEmptyRaster() throws FactoryException, ParseExceptio assertArrayEquals(expected, actual, 0.1d); } + @Test + public void testAsRasterWithNonSquarePixels() throws FactoryException, ParseException { + // Pixels are 2 world units wide and 3 tall, so the x-intercept math cannot + // silently conflate pixel-space and world-space slopes the way square + // pixels allow. Expected matrices are produced by GDAL + // (rasterio.features.rasterize) on the same grid; pixels whose centers are + // inside the geometry must be burned. + GridCoverage2D raster = + RasterConstructors.makeEmptyRaster(1, "d", 7, 6, 100, 500, 2, -3, 0, 0, 0); + Geometry geom = + Constructors.geomFromWKT( + "POLYGON ((102.7 497.4, 112.4 496.9, 104.2 483.7, 102.7 497.4))", 0); + + GridCoverage2D rasterized = + RasterConstructors.asRaster(geom, raster, "d", false, 1d, 0d, false); + double[] actual = MapAlgebra.bandAsArray(rasterized, 1); + double[] expected = + new double[] { + 0, 0, 0, 0, 0, 0, 0, + 0, 1, 1, 1, 1, 1, 0, + 0, 0, 1, 1, 1, 0, 0, + 0, 0, 1, 1, 0, 0, 0, + 0, 0, 1, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0 + }; + assertArrayEquals(expected, actual, 0.1d); + + // allTouched adds every boundary-touched pixel on the same grid + rasterized = RasterConstructors.asRaster(geom, raster, "d", true, 1d, 0d, false); + actual = MapAlgebra.bandAsArray(rasterized, 1); + expected = + new double[] { + 0, 1, 1, 1, 1, 1, 0, + 0, 1, 1, 1, 1, 1, 1, + 0, 1, 1, 1, 1, 1, 0, + 0, 1, 1, 1, 1, 0, 0, + 0, 1, 1, 1, 0, 0, 0, + 0, 0, 1, 0, 0, 0, 0 + }; + assertArrayEquals(expected, actual, 0.1d); + + // An interior ring with diagonal edges exercises hole intercepts too + geom = + Constructors.geomFromWKT( + "POLYGON ((101.1 498.9, 113.5 497.8, 112.9 483.2, 102.3 484.6, 101.1 498.9), " + + "(104.1 496.9, 110.2 496.9, 106.3 487.8, 104.1 496.9))", + 0); + rasterized = RasterConstructors.asRaster(geom, raster, "d", false, 1d, 0d, false); + actual = MapAlgebra.bandAsArray(rasterized, 1); + expected = + new double[] { + 0, 1, 1, 0, 0, 0, 0, + 0, 1, 0, 0, 0, 1, 1, + 0, 1, 1, 0, 1, 1, 1, + 0, 1, 1, 0, 1, 1, 1, + 0, 1, 1, 1, 1, 1, 1, + 0, 0, 0, 0, 0, 1, 0 + }; + assertArrayEquals(expected, actual, 0.1d); + + // Bottom-up variant of the same grid: the same triangle, vertically + // flipped in pixel space + GridCoverage2D rasterBottomUp = + RasterConstructors.makeEmptyRaster(1, "d", 7, 6, 100, 482, 2, 3, 0, 0, 0); + geom = + Constructors.geomFromWKT( + "POLYGON ((102.7 497.4, 112.4 496.9, 104.2 483.7, 102.7 497.4))", 0); + rasterized = RasterConstructors.asRaster(geom, rasterBottomUp, "d", false, 1d, 0d, false); + actual = MapAlgebra.bandAsArray(rasterized, 1); + expected = + new double[] { + 0, 0, 0, 0, 0, 0, 0, + 0, 0, 1, 0, 0, 0, 0, + 0, 0, 1, 1, 0, 0, 0, + 0, 0, 1, 1, 1, 0, 0, + 0, 1, 1, 1, 1, 1, 0, + 0, 0, 0, 0, 0, 0, 0 + }; + assertArrayEquals(expected, actual, 0.1d); + } + @Test public void testAsRasterLingString() throws FactoryException, ParseException { // Horizontal LineString From 877d6ff3ff907edcdbd138219f6e28227c354fb1 Mon Sep 17 00:00:00 2001 From: jameswillis Date: Fri, 17 Jul 2026 01:37:05 -0700 Subject: [PATCH 2/5] [GH-3111] Add seeded RS_AsRaster fuzz test against rasterio --- python/tests/sql/test_rasterize_parity.py | 114 ++++++++++++++++++++++ 1 file changed, 114 insertions(+) create mode 100644 python/tests/sql/test_rasterize_parity.py diff --git a/python/tests/sql/test_rasterize_parity.py b/python/tests/sql/test_rasterize_parity.py new file mode 100644 index 00000000000..2d547f8a724 --- /dev/null +++ b/python/tests/sql/test_rasterize_parity.py @@ -0,0 +1,114 @@ +# 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. + +import random + +import numpy as np +import rasterio.features +from affine import Affine +from shapely import from_wkt +from shapely.geometry import Polygon + +from tests.test_base import TestBase + + +def _rasterize_cases(count=100, seed=31113): + """Seeded random polygons over anisotropic north-up and south-up grids. + + The fixed seed makes the corpus deterministic, so a failure is + reproducible from the case id alone. Pixel width and height are drawn + independently, so the grids exercise the non-square pixel aspect ratios + where scanline arithmetic errors hide (square unit grids make pixel-space + and world-space slopes coincide). + """ + rng = random.Random(seed) + cases = [] + while len(cases) < count: + width, height = rng.randint(4, 12), rng.randint(4, 12) + scale_x = round(rng.uniform(0.3, 5.0), 3) + scale_y = round(rng.uniform(0.3, 5.0), 3) * rng.choice([-1, 1]) + upper_left_x = round(rng.uniform(-1000, 1000), 2) + upper_left_y = round(rng.uniform(-1000, 1000), 2) + xs = sorted([upper_left_x, upper_left_x + width * scale_x]) + ys = sorted([upper_left_y, upper_left_y + height * scale_y]) + margin_x = (xs[1] - xs[0]) * 0.05 + margin_y = (ys[1] - ys[0]) * 0.05 + + def random_point(): + return ( + round(rng.uniform(xs[0] + margin_x, xs[1] - margin_x), 3), + round(rng.uniform(ys[0] + margin_y, ys[1] - margin_y), 3), + ) + + num_points = rng.choice([3, 3, 4, 5]) + for _ in range(50): + candidate = Polygon([random_point() for _ in range(num_points)]).buffer(0) + grid_area = (xs[1] - xs[0]) * (ys[1] - ys[0]) + if candidate.geom_type == "Polygon" and candidate.area > grid_area * 0.02: + cases.append( + ( + len(cases), + width, + height, + upper_left_x, + upper_left_y, + scale_x, + scale_y, + candidate.wkt, + ) + ) + break + return cases + + +class TestRasterizeParity(TestBase): + def test_as_raster_centroid_rule_matches_gdal(self): + """RS_AsRaster with allTouched=false must burn exactly the pixels whose + centers fall inside the polygon, matching GDAL + (rasterio.features.rasterize) on the same grid. Only the centroid rule + is compared: allTouched boundary selection intentionally differs from + GDAL's. + """ + cases = _rasterize_cases() + df = self.spark.createDataFrame( + cases, + "id INT, width INT, height INT, ulx DOUBLE, uly DOUBLE, " + "sx DOUBLE, sy DOUBLE, wkt STRING", + ) + rows = df.selectExpr( + "id", + "RS_BandAsArray(RS_AsRaster(ST_GeomFromWKT(wkt), " + "RS_MakeEmptyRaster(1, 'd', width, height, ulx, uly, sx, sy, 0, 0, 0), " + "'d', false, 1.0, 0.0, false), 1) as band", + ).collect() + bands = {row["id"]: row["band"] for row in rows} + assert len(bands) == len(cases) + + for case_id, width, height, ulx, uly, sx, sy, wkt in cases: + actual = np.array(bands[case_id], dtype=float).reshape(height, width) + expected = rasterio.features.rasterize( + [(from_wkt(wkt), 1)], + out_shape=(height, width), + fill=0, + transform=Affine(sx, 0, ulx, 0, sy, uly), + all_touched=False, + dtype="uint8", + ) + assert np.array_equal(actual, expected), ( + f"case {case_id}: grid {width}x{height}, " + f"scale ({sx}, {sy}), origin ({ulx}, {uly}), {wkt}" + ) From 54c8cd59d59ada7627d01a0c639ae5433c51b2c2 Mon Sep 17 00:00:00 2001 From: jameswillis Date: Fri, 17 Jul 2026 10:33:09 -0700 Subject: [PATCH 3/5] [GH-3111] Use shapely.wkt.loads for shapely 1.x compatibility --- python/tests/sql/test_rasterize_parity.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/tests/sql/test_rasterize_parity.py b/python/tests/sql/test_rasterize_parity.py index 2d547f8a724..de56931ae4d 100644 --- a/python/tests/sql/test_rasterize_parity.py +++ b/python/tests/sql/test_rasterize_parity.py @@ -20,8 +20,8 @@ import numpy as np import rasterio.features from affine import Affine -from shapely import from_wkt from shapely.geometry import Polygon +from shapely.wkt import loads as wkt_loads from tests.test_base import TestBase @@ -101,7 +101,7 @@ def test_as_raster_centroid_rule_matches_gdal(self): for case_id, width, height, ulx, uly, sx, sy, wkt in cases: actual = np.array(bands[case_id], dtype=float).reshape(height, width) expected = rasterio.features.rasterize( - [(from_wkt(wkt), 1)], + [(wkt_loads(wkt), 1)], out_shape=(height, width), fill=0, transform=Affine(sx, 0, ulx, 0, sy, uly), From bf49afe92cfdfe1879ad37facdab086305f425fa Mon Sep 17 00:00:00 2001 From: jameswillis Date: Fri, 17 Jul 2026 10:52:35 -0700 Subject: [PATCH 4/5] [GH-3111] Cover allTouched in the fuzz corpus as a strict expected failure --- python/tests/sql/test_rasterize_parity.py | 35 +++++++++++++++++------ 1 file changed, 26 insertions(+), 9 deletions(-) diff --git a/python/tests/sql/test_rasterize_parity.py b/python/tests/sql/test_rasterize_parity.py index de56931ae4d..c5eb65b3d4a 100644 --- a/python/tests/sql/test_rasterize_parity.py +++ b/python/tests/sql/test_rasterize_parity.py @@ -17,6 +17,8 @@ import random +import pytest + import numpy as np import rasterio.features from affine import Affine @@ -76,13 +78,7 @@ def random_point(): class TestRasterizeParity(TestBase): - def test_as_raster_centroid_rule_matches_gdal(self): - """RS_AsRaster with allTouched=false must burn exactly the pixels whose - centers fall inside the polygon, matching GDAL - (rasterio.features.rasterize) on the same grid. Only the centroid rule - is compared: allTouched boundary selection intentionally differs from - GDAL's. - """ + def _assert_matches_gdal(self, all_touched): cases = _rasterize_cases() df = self.spark.createDataFrame( cases, @@ -93,7 +89,7 @@ def test_as_raster_centroid_rule_matches_gdal(self): "id", "RS_BandAsArray(RS_AsRaster(ST_GeomFromWKT(wkt), " "RS_MakeEmptyRaster(1, 'd', width, height, ulx, uly, sx, sy, 0, 0, 0), " - "'d', false, 1.0, 0.0, false), 1) as band", + f"'d', {str(all_touched).lower()}, 1.0, 0.0, false), 1) as band", ).collect() bands = {row["id"]: row["band"] for row in rows} assert len(bands) == len(cases) @@ -105,10 +101,31 @@ def test_as_raster_centroid_rule_matches_gdal(self): out_shape=(height, width), fill=0, transform=Affine(sx, 0, ulx, 0, sy, uly), - all_touched=False, + all_touched=all_touched, dtype="uint8", ) assert np.array_equal(actual, expected), ( f"case {case_id}: grid {width}x{height}, " f"scale ({sx}, {sy}), origin ({ulx}, {uly}), {wkt}" ) + + def test_as_raster_centroid_rule_matches_gdal(self): + """RS_AsRaster with allTouched=false must burn exactly the pixels whose + centers fall inside the polygon, matching GDAL + (rasterio.features.rasterize) on the same grid. + """ + self._assert_matches_gdal(all_touched=False) + + @pytest.mark.xfail( + strict=True, + raises=AssertionError, + reason="allTouched samples ring segments at fixed 0.2-pixel steps " + "(Rasterization.drawLineBresenham), so pixels the boundary clips with " + "a chord shorter than the step fall between samples and are never " + "burned where GDAL burns them. This strict xfail starts failing the " + "day the sampling is replaced with exact cell traversal.", + ) + def test_as_raster_all_touched_matches_gdal(self): + """The same corpus under allTouched=true, kept as a strict expected + failure so the known sampling divergence stays visible.""" + self._assert_matches_gdal(all_touched=True) From b7208b34fbdf21e0f38caaf9a531ff39f77ecbe1 Mon Sep 17 00:00:00 2001 From: jameswillis Date: Fri, 17 Jul 2026 11:38:13 -0700 Subject: [PATCH 5/5] [GH-3111] Add LineString sampler divergence as a strict expected failure --- python/tests/sql/test_rasterize_parity.py | 30 +++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/python/tests/sql/test_rasterize_parity.py b/python/tests/sql/test_rasterize_parity.py index c5eb65b3d4a..be9b0c2791a 100644 --- a/python/tests/sql/test_rasterize_parity.py +++ b/python/tests/sql/test_rasterize_parity.py @@ -129,3 +129,33 @@ def test_as_raster_all_touched_matches_gdal(self): """The same corpus under allTouched=true, kept as a strict expected failure so the known sampling divergence stays visible.""" self._assert_matches_gdal(all_touched=True) + + @pytest.mark.xfail( + strict=True, + raises=AssertionError, + reason="LineString rasterization uses the same fixed-step sampler " + "(Rasterization.drawLineBresenham) regardless of allTouched, so a line " + "that grazes a pixel over a chord shorter than the 0.2-pixel step is " + "not burned where GDAL burns it. This strict xfail starts failing the " + "day the sampling is replaced with exact cell traversal.", + ) + def test_as_raster_linestring_matches_gdal(self): + """A line crossing a 6x6 unit grid: every pixel the line touches must + be burned, matching GDAL. Independent of allTouched (a line has no + interior) and of #3111 (square unit pixels).""" + wkt = "LINESTRING (3.97 1.57, 0.31 3.24)" + band = self.spark.sql( + "SELECT RS_BandAsArray(RS_AsRaster(ST_GeomFromWKT('" + wkt + "'), " + "RS_MakeEmptyRaster(1, 'd', 6, 6, 0.0, 6.0, 1.0, -1.0, 0.0, 0.0, 0), " + "'d', false, 1.0, 0.0, false), 1)" + ).first()[0] + actual = np.array(band, dtype=float).reshape(6, 6) + expected = rasterio.features.rasterize( + [(wkt_loads(wkt), 1)], + out_shape=(6, 6), + fill=0, + transform=Affine(1, 0, 0, 0, -1, 6), + all_touched=True, + dtype="uint8", + ) + assert np.array_equal(actual, expected)