From d6010cce7581b202340916f1b2b5f89665767977 Mon Sep 17 00:00:00 2001 From: Vincent Date: Wed, 28 Oct 2026 23:37:07 +0800 Subject: [PATCH] [GH-3155] Transform query geometry into raster CRS for RS_ predicates Cross-CRS RS_Intersects/Contains/Within previously reprojected both the raster footprint and the query into WGS84. Corner-only chords under-cover polar and other highly curved footprints (#3323). Prefer transforming the geometry into the raster's native CRS (SedonaDB parity) so the footprint stays exact; fall back to WGS84 when that transform fails. Raster-raster predicates keep the WGS84 path for antimeridian safety. Signed-off-by: Vincent --- .../common/raster/RasterPredicates.java | 64 ++++++++++++------- .../common/raster/RasterPredicatesTest.java | 34 ++++++++++ 2 files changed, 76 insertions(+), 22 deletions(-) diff --git a/common/src/main/java/org/apache/sedona/common/raster/RasterPredicates.java b/common/src/main/java/org/apache/sedona/common/raster/RasterPredicates.java index 433b74f92e2..7cd8a57d5d6 100644 --- a/common/src/main/java/org/apache/sedona/common/raster/RasterPredicates.java +++ b/common/src/main/java/org/apache/sedona/common/raster/RasterPredicates.java @@ -53,10 +53,11 @@ public class RasterPredicates { /** - * Test if a raster intersects a query window. If both the raster and the query window have a CRS, - * the query window and the envelope of the raster will be transformed to a common CRS before - * testing for intersection. Please note that the CRS transformation will be lenient, which means - * that the transformation may not be accurate. + * Test if a raster intersects a query window. When the CRS differs, the query window is + * transformed into the raster's CRS so the raster footprint stays exact (affine pixel→world maps + * lines to lines). Falls back to comparing both sides in WGS84 if that transform fails. Please + * note that the CRS transformation will be lenient, which means that the transformation may not + * be accurate. * * @param raster the raster * @param geometry the query window @@ -307,18 +308,21 @@ private static Pair convertCRSIfNeeded( return Pair.of(rasterGeometry, queryWindow); } - // Raster has a non-authoritative CRS, or the CRS of the raster is different from the - // CRS of the query window. We'll transform both sides to a common CRS (WGS84) before - // testing for relationship. - CoordinateReferenceSystem queryWindowCRS; - queryWindowCRS = FunctionsGeoTools.sridToCRS(queryWindowSRID); - Geometry transformedQueryWindow = transformGeometryToWGS84(queryWindow, queryWindowCRS); - - // Transform the raster envelope. Here we don't use the envelope transformation method - // provided by GeoTools since it performs poorly when the raster envelope crosses the - // anti-meridian. - Geometry transformedRasterGeometry = transformGeometryToWGS84(rasterGeometry, rasterCRS); - return Pair.of(transformedRasterGeometry, transformedQueryWindow); + // Prefer transforming the (usually smaller) geometry into the raster's native CRS so the + // raster footprint stays an exact straight-edged quadrilateral. Matches SedonaDB and avoids + // the corner-chord under-coverage of reprojecting the raster footprint to WGS84 (#3155). + CoordinateReferenceSystem queryWindowCRS = FunctionsGeoTools.sridToCRS(queryWindowSRID); + try { + Geometry transformedQueryWindow = + transformGeometry(queryWindow, queryWindowCRS, rasterCRS, false); + return Pair.of(rasterGeometry, transformedQueryWindow); + } catch (RuntimeException preferredFailed) { + // Fallback: transform both sides to WGS84. Needed when the preferred direction fails + // (e.g. some antimeridian / polar degeneracies) and preserves historical join behavior. + Geometry transformedQueryWindow = transformGeometryToWGS84(queryWindow, queryWindowCRS); + Geometry transformedRasterGeometry = transformGeometryToWGS84(rasterGeometry, rasterCRS); + return Pair.of(transformedRasterGeometry, transformedQueryWindow); + } } private static Pair convertCRSIfNeeded( @@ -345,8 +349,10 @@ private static Pair convertCRSIfNeeded( return Pair.of(leftGeometry, rightGeometry); } - // Transform both sides to WGS84, and then return transformed geometries for evaluating - // predicates. + // Raster-raster: keep comparing in WGS84. Transforming one large footprint into the other's + // projected CRS silently mishandles antimeridian-crossing cases (see + // testRasterRasterPredicates*), whereas the geometry→raster path benefits from an exact + // raster footprint (#3155). Geometry transformedLeftGeometry = transformGeometryToWGS84(leftGeometry, leftCRS); Geometry transformedRightGeometry = transformGeometryToWGS84(rightGeometry, rightCRS); return Pair.of(transformedLeftGeometry, transformedRightGeometry); @@ -401,14 +407,28 @@ public static boolean isCRSMatchesSRID(CoordinateReferenceSystem crs, int srid) private static Geometry transformGeometryToWGS84( Geometry geometry, CoordinateReferenceSystem crs) { - if (crs == DefaultGeographicCRS.WGS84) { + return transformGeometry(geometry, crs, DefaultGeographicCRS.WGS84, true); + } + + /** + * Transform {@code geometry} from {@code sourceCRS} into {@code targetCRS}. When {@code + * antiMeridianSafe} is true and the source is projected, apply antimeridian-safe wrapping after + * transforming into a geographic CRS (historical WGS84 path behavior). + */ + private static Geometry transformGeometry( + Geometry geometry, + CoordinateReferenceSystem sourceCRS, + CoordinateReferenceSystem targetCRS, + boolean antiMeridianSafe) { + if (sourceCRS == targetCRS || CRS.equalsIgnoreMetadata(sourceCRS, targetCRS)) { return geometry; } try { - MathTransform transform = - CachedCRSTransformFinder.findTransform(crs, DefaultGeographicCRS.WGS84); + MathTransform transform = CachedCRSTransformFinder.findTransform(sourceCRS, targetCRS); Geometry transformedGeometry = JTS.transform(geometry, transform); - if (!(crs instanceof GeographicCRS)) { + if (antiMeridianSafe + && !(sourceCRS instanceof GeographicCRS) + && targetCRS instanceof GeographicCRS) { transformedGeometry = GeomUtils.antiMeridianSafeGeom(transformedGeometry); } return transformedGeometry; diff --git a/common/src/test/java/org/apache/sedona/common/raster/RasterPredicatesTest.java b/common/src/test/java/org/apache/sedona/common/raster/RasterPredicatesTest.java index 0f5f22acbb3..6140ce7b46b 100644 --- a/common/src/test/java/org/apache/sedona/common/raster/RasterPredicatesTest.java +++ b/common/src/test/java/org/apache/sedona/common/raster/RasterPredicatesTest.java @@ -241,6 +241,40 @@ public void testIntersectsCrossingAntiMeridian() { Assert.assertFalse(result); } + /** + * Polar rasters reprojected to WGS84 with corner-only chords lose most of the polar cap (#3323). + * Transforming the query point into the raster CRS keeps the footprint exact (#3155). + */ + @Test + public void testIntersectsPolarRasterInNativeCrs() { + // 20x20 EPSG:3413 raster covering (-1e6,-1e6)–(1e6,1e6): a 2000 km square on the north pole. + GridCoverage2D raster = + createRandomRaster( + DataBuffer.TYPE_BYTE, 20, 20, -1_000_000, 1_000_000, 100_000, 1, "EPSG:3413"); + + Geometry inside = GEOMETRY_FACTORY.createPoint(new Coordinate(45, 82)); + inside.setSRID(4326); + Assert.assertTrue(RasterPredicates.rsIntersects(raster, inside)); + Assert.assertTrue(RasterPredicates.rsContains(raster, inside)); + + inside = GEOMETRY_FACTORY.createPoint(new Coordinate(180, 85)); + inside.setSRID(4326); + Assert.assertTrue(RasterPredicates.rsIntersects(raster, inside)); + + inside = GEOMETRY_FACTORY.createPoint(new Coordinate(0, 85)); + inside.setSRID(4326); + Assert.assertTrue(RasterPredicates.rsIntersects(raster, inside)); + + inside = GEOMETRY_FACTORY.createPoint(new Coordinate(90, 85)); + inside.setSRID(4326); + Assert.assertTrue(RasterPredicates.rsIntersects(raster, inside)); + Assert.assertTrue(RasterPredicates.rsContains(raster, inside)); + + Geometry outside = GEOMETRY_FACTORY.createPoint(new Coordinate(0, 60)); + outside.setSRID(4326); + Assert.assertFalse(RasterPredicates.rsIntersects(raster, outside)); + } + @Test public void testContainsNoCrs() throws FactoryException { Geometry geometry = GEOMETRY_FACTORY.toGeometry(new Envelope(5, 10, 5, 10));