[GH-3111] Fix scanline x-intercepts on rasters with non-square pixels

common/src/main/java/org/apache/sedona/common/raster/Rasterization.java

@@ -669,7 +669,11 @@ private static Map<Double, TreeSet<Double>> 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
             }

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

python/tests/sql/test_rasterize_parity.py

@@ -0,0 +1,161 @@
+# 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 pytest
+
+import numpy as np
+import rasterio.features
+from affine import Affine
+from shapely.geometry import Polygon
+from shapely.wkt import loads as wkt_loads
+
+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 _assert_matches_gdal(self, all_touched):
+        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), "
+            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)
+
+        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(
+                [(wkt_loads(wkt), 1)],
+                out_shape=(height, width),
+                fill=0,
+                transform=Affine(sx, 0, ulx, 0, sy, uly),
+                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)
+
+    @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)