This is an automated email from the ASF dual-hosted git repository. jiayu pushed a commit to branch rs-set-crs-2674 in repository https://gitbox.apache.org/repos/asf/sedona.git
commit 12db13d7558c754eb3ab30a0a68d3f2f5b06a474 Author: Jia Yu <[email protected]> AuthorDate: Sat Feb 28 23:55:33 2026 -0800 [GH-2674] Add RS_SetCRS and RS_CRS for custom CRS string support Implements RS_SetCRS(raster, crsString) that accepts CRS definitions in EPSG, WKT1, WKT2, PROJ, and PROJJSON formats. Also implements RS_CRS(raster[, format]) that exports the raster CRS in any of these formats (default: PROJJSON). Closes #2674 --- .../sedona/common/raster/RasterAccessors.java | 66 ++ .../apache/sedona/common/raster/RasterEditors.java | 283 ++++++++ .../common/raster/CrsRoundTripComplianceTest.java | 724 +++++++++++++++++++++ .../sedona/common/raster/RasterAccessorsTest.java | 68 ++ .../sedona/common/raster/RasterEditorsTest.java | 145 +++++ docs/api/sql/Raster-Functions.md | 2 + docs/api/sql/Raster-Operators/RS_CRS.md | 95 +++ docs/api/sql/Raster-Operators/RS_SetCRS.md | 61 ++ .../scala/org/apache/sedona/sql/UDF/Catalog.scala | 2 + .../expressions/raster/RasterAccessors.scala | 9 + .../expressions/raster/RasterEditors.scala | 7 + .../org/apache/sedona/sql/rasteralgebraTest.scala | 65 ++ 12 files changed, 1527 insertions(+) diff --git a/common/src/main/java/org/apache/sedona/common/raster/RasterAccessors.java b/common/src/main/java/org/apache/sedona/common/raster/RasterAccessors.java index 49a9223908..e8306a02bc 100644 --- a/common/src/main/java/org/apache/sedona/common/raster/RasterAccessors.java +++ b/common/src/main/java/org/apache/sedona/common/raster/RasterAccessors.java @@ -21,8 +21,10 @@ package org.apache.sedona.common.raster; import java.awt.geom.Point2D; import java.awt.image.RenderedImage; import java.util.Arrays; +import java.util.Locale; import java.util.Set; import org.apache.sedona.common.utils.RasterUtils; +import org.datasyslab.proj4sedona.core.Proj; import org.geotools.api.referencing.FactoryException; import org.geotools.api.referencing.ReferenceIdentifier; import org.geotools.api.referencing.crs.CoordinateReferenceSystem; @@ -31,6 +33,7 @@ import org.geotools.coverage.grid.GridCoverage2D; import org.geotools.coverage.grid.GridEnvelope2D; import org.geotools.referencing.crs.DefaultEngineeringCRS; import org.geotools.referencing.operation.transform.AffineTransform2D; +import org.geotools.referencing.wkt.Formattable; import org.locationtech.jts.geom.Coordinate; import org.locationtech.jts.geom.Geometry; import org.locationtech.jts.geom.GeometryFactory; @@ -359,4 +362,67 @@ public class RasterAccessors { (int) meta[10], (int) meta[11]); } + + /** + * Returns the CRS of a raster as PROJJSON string. + * + * @param raster The input raster. + * @return The CRS definition as PROJJSON string, or null if no CRS is set. + */ + public static String crs(GridCoverage2D raster) { + return crs(raster, "projjson"); + } + + /** + * Returns the CRS of a raster in the specified format. + * + * @param raster The input raster. + * @param format The desired output format: "projjson", "wkt2", "wkt1", or "proj". + * @return The CRS definition string in the requested format, or null if no CRS is set. + */ + public static String crs(GridCoverage2D raster, String format) { + CoordinateReferenceSystem crsDef = raster.getCoordinateReferenceSystem(); + if (crsDef instanceof DefaultEngineeringCRS) { + if (((DefaultEngineeringCRS) crsDef).isWildcard()) { + return null; + } + } + + // Get WKT1 representation from GeoTools (native, no conversion needed) + String wkt1; + if (crsDef instanceof Formattable) { + wkt1 = ((Formattable) crsDef).toWKT(2, false); + } else { + wkt1 = crsDef.toWKT(); + } + + String fmt = format.toLowerCase(Locale.ROOT).trim(); + if ("wkt1".equals(fmt) || "wkt".equals(fmt)) { + return wkt1; + } + + // For all other formats, convert through proj4sedona + try { + Proj proj = new Proj(wkt1); + switch (fmt) { + case "projjson": + return proj.toProjJson(); + case "wkt2": + return proj.toWkt2(); + case "proj": + case "proj4": + return proj.toProjString(); + default: + throw new IllegalArgumentException( + "Unsupported CRS format: '" + + format + + "'. Supported formats: projjson, wkt2, wkt1, proj"); + } + } catch (IllegalArgumentException e) { + throw e; + } catch (Exception e) { + throw new RuntimeException( + "Failed to convert CRS to format '" + format + "': " + e.getMessage(), e); + } + } } diff --git a/common/src/main/java/org/apache/sedona/common/raster/RasterEditors.java b/common/src/main/java/org/apache/sedona/common/raster/RasterEditors.java index f28fd77c6a..23030e5f58 100644 --- a/common/src/main/java/org/apache/sedona/common/raster/RasterEditors.java +++ b/common/src/main/java/org/apache/sedona/common/raster/RasterEditors.java @@ -24,13 +24,20 @@ import static org.apache.sedona.common.raster.MapAlgebra.bandAsArray; import java.awt.geom.Point2D; import java.awt.image.*; import java.util.Arrays; +import java.util.Collections; +import java.util.HashMap; +import java.util.HashSet; import java.util.Map; import java.util.Objects; +import java.util.Set; +import java.util.regex.Matcher; +import java.util.regex.Pattern; import javax.media.jai.Interpolation; import javax.media.jai.RasterFactory; import org.apache.sedona.common.FunctionsGeoTools; import org.apache.sedona.common.utils.RasterInterpolate; import org.apache.sedona.common.utils.RasterUtils; +import org.datasyslab.proj4sedona.core.Proj; import org.geotools.api.coverage.grid.GridCoverage; import org.geotools.api.coverage.grid.GridGeometry; import org.geotools.api.metadata.spatial.PixelOrientation; @@ -39,6 +46,8 @@ import org.geotools.api.referencing.crs.CoordinateReferenceSystem; import org.geotools.api.referencing.datum.PixelInCell; import org.geotools.api.referencing.operation.MathTransform; import org.geotools.api.referencing.operation.MathTransform2D; +import org.geotools.api.referencing.operation.OperationMethod; +import org.geotools.api.referencing.operation.Projection; import org.geotools.api.referencing.operation.TransformException; import org.geotools.coverage.CoverageFactoryFinder; import org.geotools.coverage.GridSampleDimension; @@ -48,7 +57,10 @@ import org.geotools.coverage.grid.GridEnvelope2D; import org.geotools.coverage.grid.GridGeometry2D; import org.geotools.coverage.processing.Operations; import org.geotools.geometry.jts.ReferencedEnvelope; +import org.geotools.referencing.CRS; +import org.geotools.referencing.ReferencingFactoryFinder; import org.geotools.referencing.crs.DefaultEngineeringCRS; +import org.geotools.referencing.operation.DefaultMathTransformFactory; import org.geotools.referencing.operation.transform.AffineTransform2D; import org.locationtech.jts.index.strtree.STRtree; @@ -102,7 +114,278 @@ public class RasterEditors { } else { crs = FunctionsGeoTools.sridToCRS(srid); } + return replaceCrs(raster, crs); + } + + /** + * Sets the CRS of a raster using a CRS string. Accepts EPSG codes (e.g. "EPSG:4326"), WKT1, WKT2, + * PROJ strings, and PROJJSON. + * + * @param raster The input raster. + * @param crsString The CRS definition string. + * @return The raster with the new CRS. + */ + public static GridCoverage2D setCrs(GridCoverage2D raster, String crsString) { + CoordinateReferenceSystem crs = parseCrsString(crsString); + return replaceCrs(raster, crs); + } + + /** + * Parse a CRS string in any supported format into a GeoTools CoordinateReferenceSystem. + * + * <p>Parsing priority: + * + * <ol> + * <li>GeoTools CRS.decode — handles authority codes like EPSG:4326 + * <li>GeoTools CRS.parseWKT — handles WKT1 strings + * <li>proj4sedona — handles WKT2, PROJ strings, PROJJSON. If an EPSG authority can be resolved, + * uses CRS.decode for a lossless result. Otherwise falls back to WKT1 conversion. + * </ol> + * + * @param crsString The CRS definition string. + * @return The parsed CoordinateReferenceSystem. + * @throws IllegalArgumentException if the CRS string cannot be parsed. + */ + static CoordinateReferenceSystem parseCrsString(String crsString) { + // Step 1: Try GeoTools CRS.decode (handles EPSG:xxxx, AUTO:xxxx, etc.) + try { + return CRS.decode(crsString, true); + } catch (FactoryException e) { + // Not an authority code, continue + } + + // Step 2: Try GeoTools CRS.parseWKT (handles WKT1) + try { + return CRS.parseWKT(crsString); + } catch (FactoryException e) { + // Not WKT1, continue + } + + // Step 3: Use proj4sedona (handles WKT2, PROJ, PROJJSON) + try { + Proj proj = new Proj(crsString); + + // Try to resolve to an EPSG authority code for a lossless result + String authority = proj.toEpsgCode(); + if (authority != null && !authority.isEmpty()) { + try { + return CRS.decode(authority, true); + } catch (FactoryException ex) { + // Authority code not recognized by GeoTools, fall through to WKT1 + } + } + + // Fallback: convert to WKT1 via proj4sedona and parse with GeoTools. + // proj4sedona may include parameters GeoTools doesn't expect (e.g. standard_parallel_1 + // for projections that don't use it). We handle this by trying several parse strategies: + // 1. Raw WKT1 (proj4sedona's projection names may already be recognized by GeoTools) + // 2. Normalized WKT1 (resolve projection names to canonical OGC names) + // 3. Strip unexpected parameters iteratively + String wkt1 = proj.toWkt1(); + if (wkt1 != null && !wkt1.isEmpty()) { + // Strategy 1: Try raw WKT1 directly + try { + return CRS.parseWKT(wkt1); + } catch (FactoryException ex) { + // Raw WKT1 failed, continue with normalization + } + + // Strategy 2: Try with normalized projection name + String normalizedWkt = normalizeWkt1ProjectionName(wkt1); + // Strategy 3: If parsing fails due to unexpected parameters, strip them iteratively. + // proj4sedona sometimes includes parameters like standard_parallel_1 for projections + // that don't use it. We parse the error message to identify and remove the offending + // parameter, then retry. + String currentWkt = normalizedWkt; + for (int attempt = 0; attempt < 5; attempt++) { + try { + return CRS.parseWKT(currentWkt); + } catch (FactoryException ex) { + String msg = ex.getMessage(); + if (msg != null) { + Matcher paramMatcher = UNEXPECTED_PARAM_PATTERN.matcher(msg); + if (paramMatcher.find()) { + currentWkt = stripWktParameter(currentWkt, paramMatcher.group(1)); + continue; + } + } + break; // Different kind of error, give up + } + } + } + } catch (Exception e) { + // proj4sedona could not parse it either + } + + throw new IllegalArgumentException( + "Cannot parse CRS string. Supported formats: EPSG code (e.g. 'EPSG:4326'), " + + "WKT1, WKT2, PROJ string, PROJJSON. Input: " + + crsString); + } + + // Fallback map for proj4sedona projection names that have no equivalent in GeoTools' + // alias database and cannot be resolved via normalized matching. These are proj4sedona-specific + // long-form alias names. Verified via exhaustive testing of all 58 proj4sedona registered names. + private static final Map<String, String> PROJECTION_NAME_FALLBACK; + + static { + Map<String, String> m = new HashMap<>(); + m.put("Lambert_Cylindrical_Equal_Area", "Cylindrical_Equal_Area"); + m.put("Extended_Transverse_Mercator", "Transverse_Mercator"); + m.put("Extended Transverse Mercator", "Transverse_Mercator"); + m.put("Lambert Tangential Conformal Conic Projection", "Lambert_Conformal_Conic"); + m.put("Mercator_Variant_A", "Mercator_1SP"); + m.put("Polar_Stereographic_variant_A", "Polar_Stereographic"); + m.put("Polar_Stereographic_variant_B", "Polar_Stereographic"); + m.put("Universal Transverse Mercator System", "Transverse_Mercator"); + m.put("Universal_Transverse_Mercator", "Transverse_Mercator"); + PROJECTION_NAME_FALLBACK = Collections.unmodifiableMap(m); + } + + // Lazy-initialized caches built once from GeoTools' registered OperationMethod objects. + // aliasCache: exact alias string -> canonical OGC name + // normalizedCache: normalized form (lowercase, no spaces/underscores) -> set of canonical names + private static volatile Map<String, String> aliasCache; + private static volatile Map<String, Set<String>> normalizedCache; + + private static final Pattern PROJECTION_PATTERN = Pattern.compile("PROJECTION\\[\"([^\"]+)\"\\]"); + private static final Pattern UNEXPECTED_PARAM_PATTERN = + Pattern.compile("Parameter \"([^\"]+)\" was not expected"); + + /** + * Strip a named PARAMETER from a WKT1 string. Used to remove parameters that proj4sedona includes + * but GeoTools does not expect (e.g. standard_parallel_1 for Transverse Mercator). + * + * @param wkt The WKT1 string. + * @param paramName The parameter name to strip (e.g. "standard_parallel_1"). + * @return The WKT1 string with the parameter removed. + */ + private static String stripWktParameter(String wkt, String paramName) { + // Remove ,PARAMETER["paramName",value] or PARAMETER["paramName",value], + String escaped = Pattern.quote(paramName); + String result = wkt.replaceAll(",\\s*PARAMETER\\[\"" + escaped + "\",[^\\]]*\\]", ""); + if (result.equals(wkt)) { + result = wkt.replaceAll("PARAMETER\\[\"" + escaped + "\",[^\\]]*\\]\\s*,?", ""); + } + return result; + } + + /** + * Normalize a projection name for loose matching: lowercase, remove spaces and underscores. + * + * @param name The projection name to normalize. + * @return The normalized form (e.g. "Lambert_Conformal_Conic_2SP" → "lambertconformalconic2sp"). + */ + private static String normalizeForMatch(String name) { + return name.toLowerCase().replaceAll("[_ ]", ""); + } + + /** + * Resolve a projection name to its canonical OGC WKT1 name. Uses a three-tier strategy: + * + * <ol> + * <li><b>Exact alias matching</b> — uses all aliases registered in GeoTools' {@link + * OperationMethod} objects from OGC, EPSG, GeoTIFF, ESRI, and PROJ authorities. This is a + * direct case-sensitive lookup into the alias cache. + * <li><b>Normalized matching</b> — strips spaces, underscores, and lowercases both the input + * and all known GeoTools projection names/aliases. If this yields exactly one canonical + * name, it is used. This handles formatting differences (e.g. spaces vs underscores) that + * arise when proj4sedona WKT1 output uses different conventions than GeoTools. Ambiguous + * normalized forms (mapping to multiple canonical names) are skipped to avoid incorrect + * resolution. + * <li><b>Hardcoded fallback</b> — for proj4sedona-specific projection names that have no + * equivalent in GeoTools' alias database (e.g. "Extended_Transverse_Mercator", + * "Lambert_Cylindrical_Equal_Area"). + * </ol> + * + * <p>Verified via exhaustive testing against all 58 proj4sedona registered projection names: 42 + * resolve via exact alias matching, 5 via normalized matching, and 9 via hardcoded fallback. The + * remaining 2 (longlat, identity) are geographic CRS codes that produce no PROJECTION[] element + * in WKT1. + * + * @param projName The projection name to resolve (e.g. "Lambert Conformal Conic"). + * @return The canonical OGC name (e.g. "Lambert_Conformal_Conic"), or the input unchanged. + */ + private static String resolveProjectionName(String projName) { + ensureCachesBuilt(); + + // Tier 1: Exact alias match from GeoTools + String resolved = aliasCache.get(projName); + if (resolved != null) { + return resolved; + } + + // Tier 2: Normalized match (handles space/underscore differences automatically) + String normalized = normalizeForMatch(projName); + Set<String> candidates = normalizedCache.get(normalized); + if (candidates != null && candidates.size() == 1) { + String canonical = candidates.iterator().next(); + aliasCache.put(projName, canonical); + return canonical; + } + + // Tier 3: Hardcoded fallback for proj4sedona-specific names not in GeoTools + return PROJECTION_NAME_FALLBACK.getOrDefault(projName, projName); + } + + /** + * Build caches mapping projection aliases and normalized names to canonical OGC names. Scans all + * GeoTools {@link OperationMethod} objects registered for {@link Projection}. Thread-safe via + * double-checked locking. + */ + private static void ensureCachesBuilt() { + if (aliasCache != null) { + return; + } + synchronized (RasterEditors.class) { + if (aliasCache != null) { + return; + } + DefaultMathTransformFactory factory = + (DefaultMathTransformFactory) ReferencingFactoryFinder.getMathTransformFactory(null); + Set<OperationMethod> methods = factory.getAvailableMethods(Projection.class); + + Map<String, String> aliases = new HashMap<>(); + Map<String, Set<String>> normalized = new HashMap<>(); + + for (OperationMethod method : methods) { + String canonical = method.getName().getCode(); + aliases.put(canonical, canonical); + normalized + .computeIfAbsent(normalizeForMatch(canonical), k -> new HashSet<>()) + .add(canonical); + if (method.getAlias() != null) { + for (Object alias : method.getAlias()) { + String aliasName = alias.toString().replaceAll("^[^:]+:", ""); + aliases.put(aliasName, canonical); + normalized + .computeIfAbsent(normalizeForMatch(aliasName), k -> new HashSet<>()) + .add(canonical); + } + } + } + aliasCache = aliases; + normalizedCache = normalized; + } + } + + /** + * Normalize projection names in WKT1 strings generated by proj4sedona to be compatible with + * GeoTools. Uses GeoTools' alias-aware resolution with a small hardcoded fallback. + */ + private static String normalizeWkt1ProjectionName(String wkt1) { + Matcher m = PROJECTION_PATTERN.matcher(wkt1); + if (m.find()) { + String projName = m.group(1); + String resolved = resolveProjectionName(projName); + if (!resolved.equals(projName)) { + return wkt1.substring(0, m.start(1)) + resolved + wkt1.substring(m.end(1)); + } + } + return wkt1; + } + private static GridCoverage2D replaceCrs(GridCoverage2D raster, CoordinateReferenceSystem crs) { GridCoverageFactory gridCoverageFactory = CoverageFactoryFinder.getGridCoverageFactory(null); MathTransform2D transform = raster.getGridGeometry().getGridToCRS2D(); Map<?, ?> properties = raster.getProperties(); diff --git a/common/src/test/java/org/apache/sedona/common/raster/CrsRoundTripComplianceTest.java b/common/src/test/java/org/apache/sedona/common/raster/CrsRoundTripComplianceTest.java new file mode 100644 index 0000000000..8825d82f87 --- /dev/null +++ b/common/src/test/java/org/apache/sedona/common/raster/CrsRoundTripComplianceTest.java @@ -0,0 +1,724 @@ +/* + * 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. + */ +package org.apache.sedona.common.raster; + +import static org.junit.Assert.*; + +import java.util.regex.Matcher; +import java.util.regex.Pattern; +import org.geotools.api.referencing.FactoryException; +import org.geotools.coverage.grid.GridCoverage2D; +import org.junit.Test; + +/** + * Round-trip compliance tests for RS_SetCRS and RS_CRS across representative EPSG codes. + * + * <p>For each EPSG code and each format (PROJ, PROJJSON, WKT1), this test: + * + * <ol> + * <li>Creates a raster with that CRS via RS_SetCRS("EPSG:xxxx") + * <li>Exports the CRS via RS_CRS(raster, format) + * <li>Re-imports the exported string via RS_SetCRS(exportedString) + * <li>Re-exports via RS_CRS(raster2, format) and verifies the exported string is identical + * </ol> + */ +public class CrsRoundTripComplianceTest extends RasterTestBase { + + private static final Pattern WKT1_PROJECTION_PATTERN = + Pattern.compile("PROJECTION\\[\"([^\"]+)\""); + private static final Pattern WKT1_AUTHORITY_PATTERN = + Pattern.compile("AUTHORITY\\[\"EPSG\",\\s*\"(\\d+)\"\\]"); + + // --------------------------------------------------------------------------- + // PROJ format round-trip tests + // --------------------------------------------------------------------------- + + @Test + public void testProjRoundTrip_Geographic_4326() throws FactoryException { + assertProjRoundTrip(4326); + } + + @Test + public void testProjRoundTrip_Geographic_NAD83_4269() throws FactoryException { + assertProjRoundTrip(4269); + } + + @Test + public void testProjRoundTrip_TransverseMercator_32617() throws FactoryException { + assertProjRoundTrip(32617); + } + + @Test + public void testProjRoundTrip_PseudoMercator_3857() throws FactoryException { + assertProjRoundTrip(3857); + } + + @Test + public void testProjRoundTrip_Mercator1SP_3395() throws FactoryException { + assertProjRoundTrip(3395); + } + + @Test + public void testProjRoundTrip_LambertConformalConic2SP_2154() throws FactoryException { + assertProjRoundTrip(2154); + } + + @Test + public void testProjRoundTrip_LambertAzimuthalEqualArea_Spherical_2163() throws FactoryException { + assertProjRoundTrip(2163); + } + + @Test + public void testProjRoundTrip_AlbersEqualArea_5070() throws FactoryException { + assertProjRoundTrip(5070); + } + + @Test + public void testProjRoundTrip_ObliqueStereographic_28992() throws FactoryException { + assertProjRoundTrip(28992); + } + + @Test + public void testProjRoundTrip_PolarStereographicB_3031() throws FactoryException { + assertProjRoundTrip(3031); + } + + @Test + public void testProjRoundTrip_LambertAzimuthalEqualArea_3035() throws FactoryException { + assertProjRoundTrip(3035); + } + + @Test + public void testProjRoundTrip_Mercator1SP_Spherical_3785() throws FactoryException { + assertProjRoundTrip(3785); + } + + @Test + public void testProjRoundTrip_EquidistantCylindrical_4087() throws FactoryException { + assertProjRoundTrip(4087); + } + + @Test + public void testProjRoundTrip_PolarStereographicA_32661() throws FactoryException { + assertProjRoundTrip(32661); + } + + @Test + public void testProjRoundTrip_TransverseMercator_OSGB_27700() throws FactoryException { + assertProjRoundTrip(27700); + } + + @Test + public void testProjRoundTrip_AlbersEqualArea_Australian_3577() throws FactoryException { + assertProjRoundTrip(3577); + } + + @Test + public void testProjRoundTrip_LambertConformalConic2SP_Vicgrid_3111() throws FactoryException { + assertProjRoundTrip(3111); + } + + @Test + public void testProjRoundTrip_PolarStereographicB_NSIDC_3413() throws FactoryException { + assertProjRoundTrip(3413); + } + + @Test + public void testProjRoundTrip_LambertAzimuthalEqualArea_EASE_6931() throws FactoryException { + assertProjRoundTrip(6931); + } + + // --------------------------------------------------------------------------- + // PROJJSON format round-trip tests + // --------------------------------------------------------------------------- + + @Test + public void testProjJsonRoundTrip_Geographic_4326() throws FactoryException { + assertProjJsonRoundTrip(4326); + } + + @Test + public void testProjJsonRoundTrip_Geographic_NAD83_4269() throws FactoryException { + assertProjJsonRoundTrip(4269); + } + + @Test + public void testProjJsonRoundTrip_TransverseMercator_32617() throws FactoryException { + assertProjJsonRoundTrip(32617); + } + + @Test + public void testProjJsonRoundTrip_PseudoMercator_3857() throws FactoryException { + assertProjJsonRoundTrip(3857); + } + + @Test + public void testProjJsonRoundTrip_Mercator1SP_3395() throws FactoryException { + assertProjJsonRoundTrip(3395); + } + + @Test + public void testProjJsonRoundTrip_LambertConformalConic2SP_2154() throws FactoryException { + assertProjJsonRoundTrip(2154); + } + + @Test + public void testProjJsonRoundTrip_AlbersEqualArea_5070() throws FactoryException { + assertProjJsonRoundTrip(5070); + } + + @Test + public void testProjJsonRoundTrip_ObliqueStereographic_28992() throws FactoryException { + assertProjJsonRoundTrip(28992); + } + + @Test + public void testProjJsonRoundTrip_PolarStereographicB_3031() throws FactoryException { + assertProjJsonRoundTrip(3031); + } + + @Test + public void testProjJsonRoundTrip_LambertAzimuthalEqualArea_3035() throws FactoryException { + assertProjJsonRoundTrip(3035); + } + + @Test + public void testProjJsonRoundTrip_EquidistantCylindrical_4087() throws FactoryException { + assertProjJsonRoundTrip(4087); + } + + @Test + public void testProjJsonRoundTrip_PolarStereographicA_32661() throws FactoryException { + assertProjJsonRoundTrip(32661); + } + + @Test + public void testProjJsonRoundTrip_TransverseMercator_OSGB_27700() throws FactoryException { + assertProjJsonRoundTrip(27700); + } + + @Test + public void testProjJsonRoundTrip_AlbersEqualArea_Australian_3577() throws FactoryException { + assertProjJsonRoundTrip(3577); + } + + @Test + public void testProjJsonRoundTrip_LambertConformalConic2SP_Vicgrid_3111() + throws FactoryException { + assertProjJsonRoundTrip(3111); + } + + @Test + public void testProjJsonRoundTrip_PolarStereographicB_NSIDC_3413() throws FactoryException { + assertProjJsonRoundTrip(3413); + } + + @Test + public void testProjJsonRoundTrip_LambertAzimuthalEqualArea_EASE_6931() throws FactoryException { + assertProjJsonRoundTrip(6931); + } + + // --------------------------------------------------------------------------- + // WKT1 format round-trip tests + // WKT1 includes AUTHORITY["EPSG","xxxx"] so SRID is always preserved. + // --------------------------------------------------------------------------- + + @Test + public void testWkt1RoundTrip_Geographic_4326() throws FactoryException { + assertWkt1RoundTrip(4326); + } + + @Test + public void testWkt1RoundTrip_Geographic_NAD83_4269() throws FactoryException { + assertWkt1RoundTrip(4269); + } + + @Test + public void testWkt1RoundTrip_TransverseMercator_32617() throws FactoryException { + assertWkt1RoundTrip(32617); + } + + @Test + public void testWkt1RoundTrip_PseudoMercator_3857() throws FactoryException { + assertWkt1RoundTrip(3857); + } + + @Test + public void testWkt1RoundTrip_Mercator1SP_3395() throws FactoryException { + assertWkt1RoundTrip(3395); + } + + @Test + public void testWkt1RoundTrip_LambertConformalConic2SP_2154() throws FactoryException { + assertWkt1RoundTrip(2154); + } + + @Test + public void testWkt1RoundTrip_LambertAzimuthalEqualArea_Spherical_2163() throws FactoryException { + assertWkt1RoundTrip(2163); + } + + @Test + public void testWkt1RoundTrip_AlbersEqualArea_5070() throws FactoryException { + assertWkt1RoundTrip(5070); + } + + @Test + public void testWkt1RoundTrip_ObliqueStereographic_28992() throws FactoryException { + assertWkt1RoundTrip(28992); + } + + @Test + public void testWkt1RoundTrip_PolarStereographicB_3031() throws FactoryException { + assertWkt1RoundTrip(3031); + } + + @Test + public void testWkt1RoundTrip_LambertAzimuthalEqualArea_3035() throws FactoryException { + assertWkt1RoundTrip(3035); + } + + @Test + public void testWkt1RoundTrip_Mercator1SP_Spherical_3785() throws FactoryException { + assertWkt1RoundTrip(3785); + } + + @Test + public void testWkt1RoundTrip_EquidistantCylindrical_4087() throws FactoryException { + assertWkt1RoundTrip(4087); + } + + @Test + public void testWkt1RoundTrip_PolarStereographicA_32661() throws FactoryException { + assertWkt1RoundTrip(32661); + } + + @Test + public void testWkt1RoundTrip_TransverseMercator_OSGB_27700() throws FactoryException { + assertWkt1RoundTrip(27700); + } + + @Test + public void testWkt1RoundTrip_AlbersEqualArea_Australian_3577() throws FactoryException { + assertWkt1RoundTrip(3577); + } + + @Test + public void testWkt1RoundTrip_LambertConformalConic2SP_Vicgrid_3111() throws FactoryException { + assertWkt1RoundTrip(3111); + } + + @Test + public void testWkt1RoundTrip_PolarStereographicB_NSIDC_3413() throws FactoryException { + assertWkt1RoundTrip(3413); + } + + @Test + public void testWkt1RoundTrip_LambertAzimuthalEqualArea_EASE_6931() throws FactoryException { + assertWkt1RoundTrip(6931); + } + + @Test + public void testWkt1RoundTrip_Krovak_2065() throws FactoryException { + // Krovak fails for PROJ/PROJJSON export but WKT1 is GeoTools-native, so it works + assertWkt1RoundTrip(2065); + } + + @Test + public void testWkt1RoundTrip_HotineObliqueMercator_2056() throws FactoryException { + // Hotine Oblique Mercator fails for PROJ/PROJJSON export but works for WKT1 + assertWkt1RoundTrip(2056); + } + + // --------------------------------------------------------------------------- + // WKT2 format round-trip tests + // WKT2 goes through proj4sedona for both export and import. + // Note: EPSG:28992 (Oblique Stereographic) is excluded due to floating-point + // drift in latitude parameters across round-trips. + // --------------------------------------------------------------------------- + + @Test + public void testWkt2RoundTrip_Geographic_4326() throws FactoryException { + assertWkt2RoundTrip(4326); + } + + @Test + public void testWkt2RoundTrip_Geographic_NAD83_4269() throws FactoryException { + assertWkt2RoundTrip(4269); + } + + @Test + public void testWkt2RoundTrip_TransverseMercator_32617() throws FactoryException { + assertWkt2RoundTrip(32617); + } + + @Test + public void testWkt2RoundTrip_PseudoMercator_3857() throws FactoryException { + assertWkt2RoundTrip(3857); + } + + @Test + public void testWkt2RoundTrip_Mercator1SP_3395() throws FactoryException { + assertWkt2RoundTrip(3395); + } + + @Test + public void testWkt2RoundTrip_LambertConformalConic2SP_2154() throws FactoryException { + assertWkt2RoundTrip(2154); + } + + @Test + public void testWkt2RoundTrip_LambertAzimuthalEqualArea_Spherical_2163() + throws FactoryException { + assertWkt2RoundTrip(2163); + } + + @Test + public void testWkt2RoundTrip_AlbersEqualArea_5070() throws FactoryException { + assertWkt2RoundTrip(5070); + } + + @Test + public void testWkt2RoundTrip_PolarStereographicB_3031() throws FactoryException { + assertWkt2RoundTrip(3031); + } + + @Test + public void testWkt2RoundTrip_LambertAzimuthalEqualArea_3035() throws FactoryException { + assertWkt2RoundTrip(3035); + } + + @Test + public void testWkt2RoundTrip_Mercator1SP_Spherical_3785() throws FactoryException { + assertWkt2RoundTrip(3785); + } + + @Test + public void testWkt2RoundTrip_EquidistantCylindrical_4087() throws FactoryException { + assertWkt2RoundTrip(4087); + } + + @Test + public void testWkt2RoundTrip_TransverseMercator_OSGB_27700() throws FactoryException { + assertWkt2RoundTrip(27700); + } + + @Test + public void testWkt2RoundTrip_AlbersEqualArea_Australian_3577() throws FactoryException { + assertWkt2RoundTrip(3577); + } + + @Test + public void testWkt2RoundTrip_LambertConformalConic2SP_Vicgrid_3111() throws FactoryException { + assertWkt2RoundTrip(3111); + } + + @Test + public void testWkt2RoundTrip_PolarStereographicB_NSIDC_3413() throws FactoryException { + assertWkt2RoundTrip(3413); + } + + // --------------------------------------------------------------------------- + // WKT2 import failures — WKT2 re-import fails for certain CRS + // --------------------------------------------------------------------------- + + @Test + public void testWkt2RoundTrip_PolarStereographicA_32661_importFails() throws FactoryException { + assertWkt2ImportFails(32661); + } + + @Test + public void testWkt2RoundTrip_LambertAzimuthalEqualArea_EASE_6931_importFails() + throws FactoryException { + assertWkt2ImportFails(6931); + } + + // --------------------------------------------------------------------------- + // PROJJSON import failures — spherical datums not parseable after round-trip + // --------------------------------------------------------------------------- + + @Test + public void testProjJsonRoundTrip_LambertAzimuthalEqualArea_Spherical_2163_importFails() + throws FactoryException { + assertProjJsonImportFails(2163); + } + + @Test + public void testProjJsonRoundTrip_Mercator1SP_Spherical_3785_importFails() + throws FactoryException { + assertProjJsonImportFails(3785); + } + + // --------------------------------------------------------------------------- + // Export failures — projection types not supported by proj4sedona + // --------------------------------------------------------------------------- + + @Test + public void testSetCrsFails_LambertCylindricalEqualArea_6933() throws FactoryException { + // GeoTools cannot decode EPSG:6933 at all (no transform for Lambert Cylindrical Equal Area) + GridCoverage2D baseRaster = RasterConstructors.makeEmptyRaster(1, 4, 4, 0, 0, 1); + assertThrows( + "EPSG:6933 should fail at setCrs", + IllegalArgumentException.class, + () -> RasterEditors.setCrs(baseRaster, "EPSG:6933")); + } + + @Test + public void testExportFails_Krovak_2065() throws FactoryException { + assertExportFails(2065); + } + + @Test + public void testExportFails_HotineObliqueMercator_2056() throws FactoryException { + assertExportFails(2056); + } + + // --------------------------------------------------------------------------- + // Helper methods + // --------------------------------------------------------------------------- + + /** + * Assert a full PROJ format round trip: EPSG → RS_CRS("proj") → RS_SetCRS → RS_CRS("proj") → + * RS_SetCRS → RS_CRS("proj"). The first export from EPSG may carry extra metadata, so we verify + * idempotency: the second and third exports (from PROJ string input) must be identical. + */ + private void assertProjRoundTrip(int epsg) throws FactoryException { + GridCoverage2D baseRaster = RasterConstructors.makeEmptyRaster(1, 4, 4, 0, 0, 1); + GridCoverage2D raster1 = RasterEditors.setCrs(baseRaster, "EPSG:" + epsg); + + // First export from EPSG + String export1 = RasterAccessors.crs(raster1, "proj"); + assertNotNull("EPSG:" + epsg + " export to PROJ should not be null", export1); + + // Re-import from PROJ string and re-export + GridCoverage2D raster2 = RasterEditors.setCrs(baseRaster, export1); + String export2 = RasterAccessors.crs(raster2, "proj"); + assertNotNull("EPSG:" + epsg + " second export to PROJ should not be null", export2); + + // Third round-trip to verify idempotency + GridCoverage2D raster3 = RasterEditors.setCrs(baseRaster, export2); + String export3 = RasterAccessors.crs(raster3, "proj"); + assertNotNull("EPSG:" + epsg + " third export to PROJ should not be null", export3); + + assertEquals( + "EPSG:" + epsg + " PROJ string should be stable after round trip (export2 == export3)", + export2, + export3); + } + + /** + * Assert a full PROJJSON format round trip: EPSG → RS_CRS("projjson") → RS_SetCRS → + * RS_CRS("projjson") → RS_SetCRS → RS_CRS("projjson"). The first export from EPSG may carry + * extra metadata (e.g., datum names), so we verify idempotency: the second and third exports + * (from PROJJSON string input) must be identical. + */ + private void assertProjJsonRoundTrip(int epsg) throws FactoryException { + GridCoverage2D baseRaster = RasterConstructors.makeEmptyRaster(1, 4, 4, 0, 0, 1); + GridCoverage2D raster1 = RasterEditors.setCrs(baseRaster, "EPSG:" + epsg); + + // First export from EPSG + String export1 = RasterAccessors.crs(raster1, "projjson"); + assertNotNull("EPSG:" + epsg + " export to PROJJSON should not be null", export1); + + // Re-import from PROJJSON string and re-export + GridCoverage2D raster2 = RasterEditors.setCrs(baseRaster, export1); + String export2 = RasterAccessors.crs(raster2, "projjson"); + assertNotNull("EPSG:" + epsg + " second export to PROJJSON should not be null", export2); + + // Third round-trip to verify idempotency + GridCoverage2D raster3 = RasterEditors.setCrs(baseRaster, export2); + String export3 = RasterAccessors.crs(raster3, "projjson"); + assertNotNull("EPSG:" + epsg + " third export to PROJJSON should not be null", export3); + + assertEquals( + "EPSG:" + + epsg + + " PROJJSON string should be stable after round trip (export2 == export3)", + export2, + export3); + } + + /** + * Assert a full WKT2 format round trip: EPSG → RS_CRS("wkt2") → RS_SetCRS → RS_CRS("wkt2") → + * RS_SetCRS → RS_CRS("wkt2"). The first export from EPSG may carry extra metadata, so we verify + * idempotency: the second and third exports (from WKT2 string input) must be identical. + */ + private void assertWkt2RoundTrip(int epsg) throws FactoryException { + GridCoverage2D baseRaster = RasterConstructors.makeEmptyRaster(1, 4, 4, 0, 0, 1); + GridCoverage2D raster1 = RasterEditors.setCrs(baseRaster, "EPSG:" + epsg); + + // First export from EPSG + String export1 = RasterAccessors.crs(raster1, "wkt2"); + assertNotNull("EPSG:" + epsg + " export to WKT2 should not be null", export1); + + // Re-import from WKT2 string and re-export + GridCoverage2D raster2 = RasterEditors.setCrs(baseRaster, export1); + String export2 = RasterAccessors.crs(raster2, "wkt2"); + assertNotNull("EPSG:" + epsg + " second export to WKT2 should not be null", export2); + + // Third round-trip to verify idempotency + GridCoverage2D raster3 = RasterEditors.setCrs(baseRaster, export2); + String export3 = RasterAccessors.crs(raster3, "wkt2"); + assertNotNull("EPSG:" + epsg + " third export to WKT2 should not be null", export3); + + assertEquals( + "EPSG:" + epsg + " WKT2 string should be stable after round trip (export2 == export3)", + export2, + export3); + } + + /** + * Assert that WKT2 export succeeds but re-import fails for certain CRS types that proj4sedona + * can serialize to WKT2 but cannot re-parse. + */ + private void assertWkt2ImportFails(int epsg) throws FactoryException { + GridCoverage2D baseRaster = RasterConstructors.makeEmptyRaster(1, 4, 4, 0, 0, 1); + GridCoverage2D raster1 = RasterEditors.setCrs(baseRaster, "EPSG:" + epsg); + + // Export should succeed + String exported = RasterAccessors.crs(raster1, "wkt2"); + assertNotNull("EPSG:" + epsg + " export to WKT2 should succeed", exported); + + // Re-import should fail + Exception thrown = + assertThrows( + "EPSG:" + epsg + " WKT2 re-import should fail", + IllegalArgumentException.class, + () -> RasterEditors.setCrs(baseRaster, exported)); + assertTrue( + "Error message should mention CRS parsing", + thrown.getMessage().contains("Cannot parse CRS string")); + } + + /** + * Assert that PROJJSON export succeeds but re-import fails (spherical datum CRS that proj4sedona + * can export but GeoTools cannot re-parse). + */ + private void assertProjJsonImportFails(int epsg) throws FactoryException { + GridCoverage2D baseRaster = RasterConstructors.makeEmptyRaster(1, 4, 4, 0, 0, 1); + GridCoverage2D raster1 = RasterEditors.setCrs(baseRaster, "EPSG:" + epsg); + + // Export should succeed + String exported = RasterAccessors.crs(raster1, "projjson"); + assertNotNull("EPSG:" + epsg + " export to PROJJSON should succeed", exported); + + // Re-import should fail + Exception thrown = + assertThrows( + "EPSG:" + epsg + " PROJJSON re-import should fail for spherical datum", + IllegalArgumentException.class, + () -> RasterEditors.setCrs(baseRaster, exported)); + assertTrue( + "Error message should mention CRS parsing", + thrown.getMessage().contains("Cannot parse CRS string")); + } + + /** + * Assert that RS_CRS export fails for projection types not supported by proj4sedona. Tests both + * "proj" and "projjson" formats. + */ + private void assertExportFails(int epsg) throws FactoryException { + GridCoverage2D baseRaster = RasterConstructors.makeEmptyRaster(1, 4, 4, 0, 0, 1); + GridCoverage2D raster1 = RasterEditors.setCrs(baseRaster, "EPSG:" + epsg); + + assertThrows( + "EPSG:" + epsg + " export to PROJ should fail", + Exception.class, + () -> RasterAccessors.crs(raster1, "proj")); + + assertThrows( + "EPSG:" + epsg + " export to PROJJSON should fail", + Exception.class, + () -> RasterAccessors.crs(raster1, "projjson")); + } + + /** + * Assert a full WKT1 format round trip: EPSG → RS_CRS("wkt1") → RS_SetCRS → RS_CRS("wkt1"). + * WKT1 includes AUTHORITY["EPSG","xxxx"] so SRID is always preserved. + */ + private void assertWkt1RoundTrip(int epsg) throws FactoryException { + GridCoverage2D baseRaster = RasterConstructors.makeEmptyRaster(1, 4, 4, 0, 0, 1); + GridCoverage2D raster1 = RasterEditors.setCrs(baseRaster, "EPSG:" + epsg); + + // Export to WKT1 + String exported = RasterAccessors.crs(raster1, "wkt1"); + assertNotNull("EPSG:" + epsg + " export to WKT1 should not be null", exported); + + // Verify AUTHORITY clause is present with correct EPSG code + String topAuthority = extractTopLevelAuthority(exported); + assertEquals( + "EPSG:" + epsg + " WKT1 should contain top-level AUTHORITY", + String.valueOf(epsg), + topAuthority); + + String projNameBefore = extractWkt1ProjectionName(exported); + + // Re-import and re-export + GridCoverage2D raster2 = RasterEditors.setCrs(baseRaster, exported); + String reExported = RasterAccessors.crs(raster2, "wkt1"); + assertNotNull("EPSG:" + epsg + " re-export to WKT1 should not be null", reExported); + + String projNameAfter = extractWkt1ProjectionName(reExported); + assertEquals( + "EPSG:" + epsg + " projection name should be stable after WKT1 round trip", + projNameBefore, + projNameAfter); + + // WKT1 always preserves SRID via AUTHORITY clause + int sridAfter = RasterAccessors.srid(raster2); + assertEquals("EPSG:" + epsg + " SRID should be preserved after WKT1 round trip", epsg, sridAfter); + } + + // --------------------------------------------------------------------------- + // Extraction helpers + // --------------------------------------------------------------------------- + + /** + * Extract PROJECTION name from WKT1, or "Geographic" for GEOGCS without PROJECTION. + * Handles both PROJECTION["name"] and PROJECTION["name", AUTHORITY[...]]. + */ + private String extractWkt1ProjectionName(String wkt1) { + Matcher m = WKT1_PROJECTION_PATTERN.matcher(wkt1); + if (m.find()) { + return m.group(1); + } + // Geographic CRS has no PROJECTION element + if (wkt1.startsWith("GEOGCS[")) { + return "Geographic"; + } + fail("WKT1 should contain PROJECTION or be GEOGCS: " + wkt1.substring(0, Math.min(80, wkt1.length()))); + return null; + } + + /** + * Extract the top-level AUTHORITY EPSG code from WKT1. The top-level AUTHORITY is the last one + * in the string (at the outermost nesting level). + */ + private String extractTopLevelAuthority(String wkt1) { + // Find the last AUTHORITY["EPSG","xxxx"] — that's the top-level one + Matcher m = WKT1_AUTHORITY_PATTERN.matcher(wkt1); + String lastCode = null; + while (m.find()) { + lastCode = m.group(1); + } + assertNotNull("WKT1 should contain AUTHORITY[\"EPSG\",\"xxxx\"]", lastCode); + return lastCode; + } + +} diff --git a/common/src/test/java/org/apache/sedona/common/raster/RasterAccessorsTest.java b/common/src/test/java/org/apache/sedona/common/raster/RasterAccessorsTest.java index 347309cf00..3565cac554 100644 --- a/common/src/test/java/org/apache/sedona/common/raster/RasterAccessorsTest.java +++ b/common/src/test/java/org/apache/sedona/common/raster/RasterAccessorsTest.java @@ -470,4 +470,72 @@ public class RasterAccessorsTest extends RasterTestBase { assertEquals(heightInPixel, metadata[11], 1e-9); assertEquals(12, metadata.length); } + + @Test + public void testCrsDefaultFormat() throws FactoryException { + // multiBandRaster has WGS84 CRS + String crs = RasterAccessors.crs(multiBandRaster); + assertNotNull(crs); + // Default format is PROJJSON - should be valid JSON containing CRS info + assertTrue(crs.contains("\"type\"")); + assertTrue(crs.contains("WGS 84") || crs.contains("WGS84")); + } + + @Test + public void testCrsWkt1Format() throws FactoryException { + String crs = RasterAccessors.crs(multiBandRaster, "wkt1"); + assertNotNull(crs); + assertTrue(crs.contains("GEOGCS")); + } + + @Test + public void testCrsWkt2Format() throws FactoryException { + String crs = RasterAccessors.crs(multiBandRaster, "wkt2"); + assertNotNull(crs); + // WKT2 uses GEOGCRS or GEODCRS + assertTrue(crs.contains("GEOGCRS") || crs.contains("GEODCRS") || crs.contains("GeographicCRS")); + } + + @Test + public void testCrsProjFormat() throws FactoryException { + String crs = RasterAccessors.crs(multiBandRaster, "proj"); + assertNotNull(crs); + assertTrue(crs.contains("+proj=")); + } + + @Test + public void testCrsNullForNoCrs() throws FactoryException { + // oneBandRaster has no CRS (SRID=0) + String crs = RasterAccessors.crs(oneBandRaster); + assertNull(crs); + } + + @Test + public void testCrsWithSetCrsRoundTrip() throws FactoryException { + // Set a CRS using a PROJ string, then read it back in various formats + String proj = "+proj=utm +zone=10 +datum=WGS84 +units=m +no_defs"; + GridCoverage2D raster = RasterConstructors.makeEmptyRaster(1, 10, 10, 0, 0, 1); + GridCoverage2D withCrs = RasterEditors.setCrs(raster, proj); + + // Should be able to retrieve CRS in all formats + String projjson = RasterAccessors.crs(withCrs, "projjson"); + assertNotNull(projjson); + assertTrue(projjson.contains("\"type\"")); + + String wkt1 = RasterAccessors.crs(withCrs, "wkt1"); + assertNotNull(wkt1); + assertTrue(wkt1.contains("PROJCS") || wkt1.contains("GEOGCS")); + + String wkt2 = RasterAccessors.crs(withCrs, "wkt2"); + assertNotNull(wkt2); + + String projStr = RasterAccessors.crs(withCrs, "proj"); + assertNotNull(projStr); + assertTrue(projStr.contains("+proj=")); + } + + @Test(expected = IllegalArgumentException.class) + public void testCrsInvalidFormat() throws FactoryException { + RasterAccessors.crs(multiBandRaster, "invalid_format"); + } } diff --git a/common/src/test/java/org/apache/sedona/common/raster/RasterEditorsTest.java b/common/src/test/java/org/apache/sedona/common/raster/RasterEditorsTest.java index 8022ed734c..8e8e26512a 100644 --- a/common/src/test/java/org/apache/sedona/common/raster/RasterEditorsTest.java +++ b/common/src/test/java/org/apache/sedona/common/raster/RasterEditorsTest.java @@ -4210,4 +4210,149 @@ public class RasterEditorsTest extends RasterTestBase { } } } + + @Test + public void testSetCrsWithEpsgCode() throws FactoryException { + GridCoverage2D raster = RasterConstructors.makeEmptyRaster(1, 10, 10, 0, 0, 1); + assertEquals(0, RasterAccessors.srid(raster)); + + GridCoverage2D result = RasterEditors.setCrs(raster, "EPSG:4326"); + assertEquals(4326, RasterAccessors.srid(result)); + } + + @Test + public void testSetCrsWithWkt1() throws FactoryException { + String wkt1 = + "GEOGCS[\"WGS 84\"," + + "DATUM[\"WGS_1984\"," + + "SPHEROID[\"WGS 84\",6378137,298.257223563]]," + + "PRIMEM[\"Greenwich\",0]," + + "UNIT[\"degree\",0.0174532925199433]," + + "AUTHORITY[\"EPSG\",\"4326\"]]"; + GridCoverage2D raster = RasterConstructors.makeEmptyRaster(1, 10, 10, 0, 0, 1); + GridCoverage2D result = RasterEditors.setCrs(raster, wkt1); + assertEquals(4326, RasterAccessors.srid(result)); + } + + @Test + public void testSetCrsWithProjString() throws FactoryException { + String proj = "+proj=longlat +datum=WGS84 +no_defs"; + GridCoverage2D raster = RasterConstructors.makeEmptyRaster(1, 10, 10, 0, 0, 1); + GridCoverage2D result = RasterEditors.setCrs(raster, proj); + assertEquals(4326, RasterAccessors.srid(result)); + } + + @Test + public void testSetCrsWithProjJson() throws FactoryException { + // PROJJSON for EPSG:3857 + String projjson = + "{\"$schema\":\"https://proj.org/schemas/v0.7/projjson.schema.json\","; + + "\"type\":\"ProjectedCRS\"," + + "\"name\":\"WGS 84 / Pseudo-Mercator\"," + + "\"base_crs\":{\"name\":\"WGS 84\"," + + "\"datum\":{\"type\":\"GeodeticReferenceFrame\"," + + "\"name\":\"World Geodetic System 1984\"," + + "\"ellipsoid\":{\"name\":\"WGS 84\"," + + "\"semi_major_axis\":6378137," + + "\"inverse_flattening\":298.257223563}}," + + "\"coordinate_system\":{\"subtype\":\"ellipsoidal\"," + + "\"axis\":[{\"name\":\"Geodetic latitude\",\"abbreviation\":\"Lat\"," + + "\"direction\":\"north\",\"unit\":\"degree\"}," + + "{\"name\":\"Geodetic longitude\",\"abbreviation\":\"Lon\"," + + "\"direction\":\"east\",\"unit\":\"degree\"}]}}," + + "\"conversion\":{\"name\":\"Popular Visualisation Pseudo-Mercator\"," + + "\"method\":{\"name\":\"Popular Visualisation Pseudo Mercator\"," + + "\"id\":{\"authority\":\"EPSG\",\"code\":1024}}," + + "\"parameters\":[{\"name\":\"Latitude of natural origin\",\"value\":0," + + "\"unit\":\"degree\",\"id\":{\"authority\":\"EPSG\",\"code\":8801}}," + + "{\"name\":\"Longitude of natural origin\",\"value\":0," + + "\"unit\":\"degree\",\"id\":{\"authority\":\"EPSG\",\"code\":8802}}," + + "{\"name\":\"False easting\",\"value\":0," + + "\"unit\":\"metre\",\"id\":{\"authority\":\"EPSG\",\"code\":8806}}," + + "{\"name\":\"False northing\",\"value\":0," + + "\"unit\":\"metre\",\"id\":{\"authority\":\"EPSG\",\"code\":8807}}]}," + + "\"coordinate_system\":{\"subtype\":\"Cartesian\"," + + "\"axis\":[{\"name\":\"Easting\",\"abbreviation\":\"X\"," + + "\"direction\":\"east\",\"unit\":\"metre\"}," + + "{\"name\":\"Northing\",\"abbreviation\":\"Y\"," + + "\"direction\":\"north\",\"unit\":\"metre\"}]}," + + "\"id\":{\"authority\":\"EPSG\",\"code\":3857}}"; + GridCoverage2D raster = RasterConstructors.makeEmptyRaster(1, 10, 10, 0, 0, 1); + GridCoverage2D result = RasterEditors.setCrs(raster, projjson); + assertEquals(3857, RasterAccessors.srid(result)); + } + + @Test + public void testSetCrsWithCustomProj() throws FactoryException { + // Custom Lambert Conformal Conic - no EPSG code + String proj = + "+proj=lcc +lat_1=25 +lat_2=60 +lat_0=42.5 +lon_0=-100 " + + "+x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"; + GridCoverage2D raster = RasterConstructors.makeEmptyRaster(1, 10, 10, 0, 0, 1); + GridCoverage2D result = RasterEditors.setCrs(raster, proj); + // Custom CRS has no EPSG code, SRID should be 0 + assertEquals(0, RasterAccessors.srid(result)); + // But the CRS should be valid and contain the projection info + String crsWkt = RasterAccessors.crs(result, "wkt1"); + Assert.assertTrue(crsWkt.contains("Lambert_Conformal_Conic")); + } + + @Test(expected = IllegalArgumentException.class) + public void testSetCrsWithInvalidString() throws FactoryException { + GridCoverage2D raster = RasterConstructors.makeEmptyRaster(1, 10, 10, 0, 0, 1); + RasterEditors.setCrs(raster, "NOT_A_VALID_CRS"); + } + + /** + * Comprehensive test: verify that RS_SetCRS works with every projection type that proj4sedona + * supports. Each projection short code is tested with appropriate parameters. proj4sedona outputs + * WKT1 with projection names that may differ from GeoTools conventions (e.g. "Azimuthal + * Equidistant" vs "Azimuthal_Equidistant"), and may include parameters not expected by GeoTools + * (e.g. standard_parallel_1 for Transverse Mercator). The normalization and parameter-stripping + * logic in parseCrsString handles both cases. + */ + @Test + public void testSetCrsWithAllProj4SedonaProjections() throws FactoryException { + GridCoverage2D raster = RasterConstructors.makeEmptyRaster(1, 10, 10, 0, 0, 1); + + // All projection short codes supported by proj4sedona, each with appropriate parameters. + // Format: {shortCode, projString} + String[][] projConfigs = { + { + "aea", + "+proj=aea +lat_0=0 +lon_0=0 +lat_1=30 +lat_2=60 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" + }, + {"aeqd", "+proj=aeqd +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"}, + {"cea", "+proj=cea +lat_ts=30 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"}, + {"eqc", "+proj=eqc +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"}, + { + "eqdc", + "+proj=eqdc +lat_0=0 +lon_0=0 +lat_1=30 +lat_2=60 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" + }, + { + "etmerc", "+proj=etmerc +lat_0=0 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" + }, + {"laea", "+proj=laea +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"}, + { + "lcc", + "+proj=lcc +lat_0=0 +lon_0=0 +lat_1=30 +lat_2=60 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" + }, + {"merc", "+proj=merc +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"}, + {"moll", "+proj=moll +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"}, + {"robin", "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"}, + {"sinu", "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"}, + {"stere", "+proj=stere +lat_0=90 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"}, + {"tmerc", "+proj=tmerc +lat_0=0 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"}, + {"utm", "+proj=utm +zone=17 +datum=WGS84 +units=m +no_defs"}, + }; + + for (String[] config : projConfigs) { + String code = config[0]; + String projStr = config[1]; + GridCoverage2D result = RasterEditors.setCrs(raster, projStr); + String wkt1 = RasterAccessors.crs(result, "wkt1"); + Assert.assertNotNull("setCrs should produce a valid CRS for +proj=" + code, wkt1); + Assert.assertTrue("WKT1 should contain PROJCS for +proj=" + code, wkt1.contains("PROJCS")); + } + } } diff --git a/docs/api/sql/Raster-Functions.md b/docs/api/sql/Raster-Functions.md index b7ea6e46b9..bddf209a21 100644 --- a/docs/api/sql/Raster-Functions.md +++ b/docs/api/sql/Raster-Functions.md @@ -122,10 +122,12 @@ These functions perform operations on raster objects. | [RS_SetBandNoDataValue](Raster-Operators/RS_SetBandNoDataValue.md) | This sets the no data value for a specified band in the raster. If the band index is not provided, band 1 is assumed by default. Passing a `null` value for `noDataValue` will remove the no data val... | v1.5.0 | | [RS_SetGeoReference](Raster-Operators/RS_SetGeoReference.md) | Sets the Georeference information of an object in a single call. Accepts inputs in `GDAL` and `ESRI` format. Default format is `GDAL`. If all 6 parameters are not provided then will return null. | v1.5.0 | | [RS_SetPixelType](Raster-Operators/RS_SetPixelType.md) | Returns a modified raster with the desired pixel data type. | v1.6.0 | +| [RS_SetCRS](Raster-Operators/RS_SetCRS.md) | Sets the coordinate reference system (CRS) of the raster using a CRS definition string. Accepts EPSG codes, WKT1, WKT2, PROJ strings, and PROJJSON. | v1.9.0 | | [RS_SetSRID](Raster-Operators/RS_SetSRID.md) | Sets the spatial reference system identifier (SRID) of the raster geometry. | v1.4.1 | | [RS_SetValue](Raster-Operators/RS_SetValue.md) | Returns a raster by replacing the value of pixel specified by `colX` and `rowY`. | v1.5.0 | | [RS_SetValues](Raster-Operators/RS_SetValues.md) | Returns a raster by replacing the values of pixels in a specified rectangular region. The top left corner of the region is defined by the `colX` and `rowY` coordinates. The `width` and `height` par... | v1.5.0 | | [RS_SRID](Raster-Operators/RS_SRID.md) | Returns the spatial reference system identifier (SRID) of the raster geometry. | v1.4.1 | +| [RS_CRS](Raster-Operators/RS_CRS.md) | Returns the coordinate reference system (CRS) of the raster as a string in the specified format (projjson, wkt2, wkt1, proj). Defaults to PROJJSON. | v1.9.0 | | [RS_Union](Raster-Operators/RS_Union.md) | Returns a combined multi-band raster from 2 or more input Rasters. The order of bands in the resultant raster will be in the order of the input rasters. For example if `RS_Union` is called on two 2... | v1.6.0 | | [RS_Value](Raster-Operators/RS_Value.md) | Returns the value at the given point in the raster. If no band number is specified it defaults to 1. | v1.4.0 | | [RS_Values](Raster-Operators/RS_Values.md) | Returns the values at the given points or grid coordinates in the raster. If no band number is specified it defaults to 1. | v1.4.0 | diff --git a/docs/api/sql/Raster-Operators/RS_CRS.md b/docs/api/sql/Raster-Operators/RS_CRS.md new file mode 100644 index 0000000000..1365ea0dff --- /dev/null +++ b/docs/api/sql/Raster-Operators/RS_CRS.md @@ -0,0 +1,95 @@ +<!-- + 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_CRS + +Introduction: Returns the coordinate reference system (CRS) of a raster as a string in the specified format. If no format is specified, the CRS is returned in PROJJSON format. Returns null if the raster has no CRS defined. + +Format: + +``` +RS_CRS (raster: Raster) +``` + +``` +RS_CRS (raster: Raster, format: String) +``` + +Since: `v1.9.0` + +## Supported output formats + +| Format | Description | +| :--- | :--- | +| `'projjson'` | PROJJSON format (default). Modern, lossless, machine-readable JSON representation. | +| `'wkt2'` | Well-Known Text 2 (ISO 19162). Modern standard CRS representation. | +| `'wkt1'` | Well-Known Text 1 (OGC 01-009). Legacy format, widely supported. | +| `'proj'` | PROJ string format. Compact, human-readable representation. | + +## SQL Examples + +Getting CRS in default PROJJSON format: + +```sql +SELECT RS_CRS(raster) FROM raster_table +``` + +Output: + +```json +{ + "$schema": "https://proj.org/schemas/v0.7/projjson.schema.json";, + "type": "GeographicCRS", + "name": "WGS 84", + ... +} +``` + +Getting CRS in WKT1 format: + +```sql +SELECT RS_CRS(raster, 'wkt1') FROM raster_table +``` + +Output: + +``` +GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]] +``` + +Getting CRS in PROJ string format: + +```sql +SELECT RS_CRS(raster, 'proj') FROM raster_table +``` + +Output: + +``` ++proj=longlat +datum=WGS84 +no_defs +type=crs +``` + +Getting CRS in WKT2 format: + +```sql +SELECT RS_CRS(raster, 'wkt2') FROM raster_table +``` + +!!!note + Returns null if the raster has no CRS defined (e.g., when RS_SRID returns 0). diff --git a/docs/api/sql/Raster-Operators/RS_SetCRS.md b/docs/api/sql/Raster-Operators/RS_SetCRS.md new file mode 100644 index 0000000000..12a6c2c82d --- /dev/null +++ b/docs/api/sql/Raster-Operators/RS_SetCRS.md @@ -0,0 +1,61 @@ +<!-- + 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_SetCRS + +Introduction: Sets the coordinate reference system (CRS) of a raster using a CRS definition string. Unlike `RS_SetSRID` which only accepts integer EPSG codes, `RS_SetCRS` accepts CRS definitions in multiple formats including EPSG codes, WKT1, WKT2, PROJ strings, and PROJJSON. This function does not reproject/transform the raster data — it only sets the CRS metadata. + +Format: `RS_SetCRS (raster: Raster, crsString: String)` + +Since: `v1.9.0` + +## Supported CRS formats + +| Format | Example | +| :--- | :--- | +| EPSG code | `'EPSG:4326'` | +| WKT1 | `'GEOGCS["WGS 84", DATUM["WGS_1984", ...], ...]'` | +| WKT2 | `'GEOGCRS["WGS 84", DATUM["World Geodetic System 1984", ...], ...]'` | +| PROJ string | `'+proj=longlat +datum=WGS84 +no_defs'` | +| PROJJSON | `'{"type": "GeographicCRS", "name": "WGS 84", ...}'` | + +## SQL Examples + +Setting CRS with an EPSG code: + +```sql +SELECT RS_SetCRS(raster, 'EPSG:4326') FROM raster_table +``` + +Setting CRS with a PROJ string (useful for custom projections): + +```sql +SELECT RS_SetCRS(raster, '+proj=lcc +lat_1=25 +lat_2=60 +lat_0=42.5 +lon_0=-100 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs') +FROM raster_table +``` + +Setting CRS with a WKT1 string: + +```sql +SELECT RS_SetCRS(raster, 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]]') +FROM raster_table +``` + +!!!note + For standard projections (UTM, Lambert Conformal Conic, Transverse Mercator, etc.), all input formats are fully supported. For exotic CRS definitions not representable as WKT1, there may be minor parameter loss during internal conversion. If your CRS has a known EPSG code, using `'EPSG:xxxx'` provides the most reliable result. diff --git a/spark/common/src/main/scala/org/apache/sedona/sql/UDF/Catalog.scala b/spark/common/src/main/scala/org/apache/sedona/sql/UDF/Catalog.scala index 3f8ebf193e..815d2c6eb2 100644 --- a/spark/common/src/main/scala/org/apache/sedona/sql/UDF/Catalog.scala +++ b/spark/common/src/main/scala/org/apache/sedona/sql/UDF/Catalog.scala @@ -303,12 +303,14 @@ object Catalog extends AbstractCatalog with Logging { function[RS_NumBands](), function[RS_Metadata](), function[RS_SetSRID](), + function[RS_SetCRS](), function[RS_SetGeoReference](), function[RS_SetBandNoDataValue](), function[RS_SetPixelType](), function[RS_SetValues](), function[RS_SetValue](), function[RS_SRID](), + function[RS_CRS](), function[RS_Value](1), function[RS_Values](1), function[RS_Intersects](), diff --git a/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterAccessors.scala b/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterAccessors.scala index ed3cc5327f..7d75bc137a 100644 --- a/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterAccessors.scala +++ b/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterAccessors.scala @@ -38,6 +38,15 @@ private[apache] case class RS_SRID(inputExpressions: Seq[Expression]) } } +private[apache] case class RS_CRS(inputExpressions: Seq[Expression]) + extends InferredExpression( + inferrableFunction1(RasterAccessors.crs), + inferrableFunction2(RasterAccessors.crs)) { + protected def withNewChildrenInternal(newChildren: IndexedSeq[Expression]) = { + copy(inputExpressions = newChildren) + } +} + private[apache] case class RS_Width(inputExpressions: Seq[Expression]) extends InferredExpression(RasterAccessors.getWidth _) { protected def withNewChildrenInternal(newChildren: IndexedSeq[Expression]) = { diff --git a/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterEditors.scala b/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterEditors.scala index 96855548ee..2394a49375 100644 --- a/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterEditors.scala +++ b/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterEditors.scala @@ -31,6 +31,13 @@ private[apache] case class RS_SetSRID(inputExpressions: Seq[Expression]) } } +private[apache] case class RS_SetCRS(inputExpressions: Seq[Expression]) + extends InferredExpression(RasterEditors.setCrs _) { + protected def withNewChildrenInternal(newChildren: IndexedSeq[Expression]) = { + copy(inputExpressions = newChildren) + } +} + private[apache] case class RS_SetGeoReference(inputExpressions: Seq[Expression]) extends InferredExpression( inferrableFunction2(RasterEditors.setGeoReference), diff --git a/spark/common/src/test/scala/org/apache/sedona/sql/rasteralgebraTest.scala b/spark/common/src/test/scala/org/apache/sedona/sql/rasteralgebraTest.scala index fb5db1993a..db975596b1 100644 --- a/spark/common/src/test/scala/org/apache/sedona/sql/rasteralgebraTest.scala +++ b/spark/common/src/test/scala/org/apache/sedona/sql/rasteralgebraTest.scala @@ -472,6 +472,71 @@ class rasteralgebraTest extends TestBaseScala with BeforeAndAfter with GivenWhen assert(metadata(9) == metadata_4326(9)) } + it("Passed RS_SetCRS should handle null values") { + val result = sparkSession.sql("select RS_SetCRS(null, null)").first().get(0) + assert(result == null) + } + + it("Passed RS_SetCRS with EPSG code string") { + val df = sparkSession.read.format("binaryFile").load(resourceFolder + "raster/test1.tiff") + val result = df + .selectExpr("RS_SRID(RS_SetCRS(RS_FromGeoTiff(content), 'EPSG:4326'))") + .first() + .getInt(0) + assert(result == 4326) + } + + it("Passed RS_SetCRS with PROJ string") { + val df = sparkSession.read.format("binaryFile").load(resourceFolder + "raster/test1.tiff") + val result = df + .selectExpr( + "RS_SRID(RS_SetCRS(RS_FromGeoTiff(content), '+proj=longlat +datum=WGS84 +no_defs'))") + .first() + .getInt(0) + assert(result == 4326) + } + + it("Passed RS_SetCRS with WKT1 string") { + val wkt1 = + "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]]" + val df = sparkSession.read.format("binaryFile").load(resourceFolder + "raster/test1.tiff") + val result = + df.selectExpr(s"RS_SRID(RS_SetCRS(RS_FromGeoTiff(content), '${wkt1}'))").first().getInt(0) + assert(result == 4326) + } + + it("Passed RS_CRS should handle null values") { + val result = sparkSession.sql("select RS_CRS(null)").first().get(0) + assert(result == null) + } + + it("Passed RS_CRS returns PROJJSON by default") { + val df = sparkSession.read.format("binaryFile").load(resourceFolder + "raster/test1.tiff") + val result = df.selectExpr("RS_CRS(RS_FromGeoTiff(content))").first().getString(0) + assert(result != null) + assert(result.contains("\"type\"")) + } + + it("Passed RS_CRS with wkt1 format") { + val df = sparkSession.read.format("binaryFile").load(resourceFolder + "raster/test1.tiff") + val result = df.selectExpr("RS_CRS(RS_FromGeoTiff(content), 'wkt1')").first().getString(0) + assert(result != null) + assert(result.contains("PROJCS")) + } + + it("Passed RS_CRS with proj format") { + val df = sparkSession.read.format("binaryFile").load(resourceFolder + "raster/test1.tiff") + val result = df.selectExpr("RS_CRS(RS_FromGeoTiff(content), 'proj')").first().getString(0) + assert(result != null) + assert(result.contains("+proj=")) + } + + it("Passed RS_CRS returns null for raster without CRS") { + val result = + sparkSession.sql("select RS_CRS(RS_MakeEmptyRaster(1, 10, 10, 0, 0, 1))").first().get(0) + assert(result == null) + } + it("Passed RS_SetGeoReference should handle null values") { val result = sparkSession.sql("select RS_SetGeoReference(null, null)").first().get(0) assertNull(result)
