james-willis opened a new issue, #3111: URL: https://github.com/apache/sedona/issues/3111
## Expected behavior With `allTouched = false`, `RS_AsRaster` is documented to select "only pixels whose centroids intersect the geometry" ([docs](https://sedona.apache.org/latest/api/sql/Raster-Operators/RS_AsRaster/)). This is also the rule GDAL (`gdal_rasterize`, `rasterio.features.rasterize`) and PostGIS `ST_AsRaster` implement, so the three should agree on which pixels are burned. ## Actual behavior On rasters whose pixels are not square in world units (`|scaleX| != |scaleY|`), pixels whose centers are clearly inside the polygon are not burned. Repro on a 7x6 grid with 2x3 world units per pixel (Sedona 1.9.0, Spark 4.0): ```sql SELECT RS_BandAsArray(RS_AsRaster( ST_GeomFromWKT('POLYGON ((102.7 497.4, 112.4 496.9, 104.2 483.7, 102.7 497.4))'), RS_MakeEmptyRaster(1, 'B', 7, 6, 100.0, 500.0, 2.0, -3.0, 0.0, 0.0, 0), 'B', false, 1.0, 0.0, false), 1) ``` Sedona returns (reshaped to 6x7): ``` 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 1 1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 ``` GDAL/rasterio and PostGIS return, for the same geometry and grid: ``` 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 ``` The dropped pixels are not boundary ties. For example pixel (row 2, col 4) has its center at world (109.0, 492.5), which is 0.566 world units inside the polygon boundary; pixel (1, 5) center (111.0, 495.5) is 0.450 units inside. Both are unambiguously "centroid intersects the geometry". ## Root cause `Rasterization.computeScanlineIntersections` computes edge x-intercepts by mixing pixel-space and world-space quantities ([Rasterization.java#L655-L677 @ 0299fa8](https://github.com/apache/sedona/blob/0299fa8a380a/common/src/main/java/org/apache/sedona/common/raster/Rasterization.java#L655-L677)): ```java double p1X = (worldP1.x - params.upperLeftX) / params.scaleX; // pixel space double p1Y = (worldP1.y - params.upperLeftY) / params.scaleY; // pixel space ... double slope = (worldP2.y - worldP1.y) / (worldP2.x - worldP1.x); // WORLD-space slope ... double xIntercept = p1X + ((p1Y - y) / slope); // pixel dy / world slope ``` `(p1Y - y)` is a pixel-space delta but `slope` is in world units, so the intercept is missing the scale factors. The correct pixel-space intercept is: ```java double xIntercept = p1X + ((y - p1Y) * params.scaleY / slope / params.scaleX); ``` I verified this both ways: a line-for-line transcription of the current Java reproduces Sedona's output matrix exactly, and the same transcription with the scale factors restored reproduces the GDAL/PostGIS matrix exactly. Why the test suite doesn't catch it: - For square north-up pixels (`scaleX = s`, `scaleY = -s`) the missing factor `-scaleY/scaleX` equals 1, so the formula is coincidentally correct — unit-scale fixtures can't expose it. - Vertical edges take the constant-x branch (pure pixel space, correct) and horizontal edges are skipped, so axis-aligned rectangles also rasterize correctly. Only non-axis-aligned edges on non-square pixels are affected. ## Impact Correctness. `Rasterization` is reached via `RasterConstructors.asRaster` from `RS_AsRaster`, `RS_Clip` (mask selection), `RS_ZonalStats` / `RS_ZonalStatsAll` (zone pixel selection), and `RS_SetValues` (geometry variant). On any raster with rectangular non-square pixels, all of these operate on a wrong pixel set — e.g. zonal statistics silently aggregate a different set of pixels than the documented centroid rule (and than GDAL/PostGIS) for the same zone polygon. Affected versions: reproduced on 1.9.0 (released jars) and confirmed present in current master (`0299fa8`); the formula predates and was unchanged by SEDONA-740/SEDONA-754/SEDONA-755. -- This is an automated message from the Apache Git Service. To respond to the message, please log on to GitHub and use the URL above to go to the specific comment. To unsubscribe, e-mail: [email protected] For queries about this service, please contact Infrastructure at: [email protected]
