From 7d38d3e8482a3a6e1f79ad10c8523180b743dbde Mon Sep 17 00:00:00 2001 From: Furqaanahmed Khan Date: Fri, 1 Dec 2023 15:47:53 -0500 Subject: [PATCH] feat: set default CRS as 4326 --- .../common/raster/PixelFunctionEditors.java | 10 +++--- .../sedona/common/raster/PixelFunctions.java | 19 ++++++---- .../common/raster/RasterBandEditors.java | 10 +++--- .../sedona/common/utils/RasterUtils.java | 35 +++++++++++++++++++ .../common/raster/FunctionEditorsTest.java | 11 +++--- .../sedona/common/raster/FunctionsTest.java | 6 ++-- .../apache/sedona/sql/rasteralgebraTest.scala | 28 +++++++-------- 7 files changed, 77 insertions(+), 42 deletions(-) diff --git a/common/src/main/java/org/apache/sedona/common/raster/PixelFunctionEditors.java b/common/src/main/java/org/apache/sedona/common/raster/PixelFunctionEditors.java index 66d72fdddc..004ae04a8a 100644 --- a/common/src/main/java/org/apache/sedona/common/raster/PixelFunctionEditors.java +++ b/common/src/main/java/org/apache/sedona/common/raster/PixelFunctionEditors.java @@ -18,6 +18,7 @@ */ package org.apache.sedona.common.raster; +import org.apache.commons.lang3.tuple.Pair; import org.apache.sedona.common.Functions; import org.apache.sedona.common.utils.RasterUtils; import org.geotools.coverage.grid.GridCoverage2D; @@ -109,12 +110,9 @@ public static GridCoverage2D setValues(GridCoverage2D raster, int band, int colX public static GridCoverage2D setValues(GridCoverage2D raster, int band, Geometry geom, double value, boolean keepNoData) throws FactoryException, TransformException { RasterUtils.ensureBand(raster, band); - if(RasterAccessors.srid(raster) != geom.getSRID()) { - // implicitly converting geometry CRS to raster CRS - geom = RasterUtils.convertCRSIfNeeded(geom, raster.getCoordinateReferenceSystem()); - // have to set the SRID as RasterUtils.convertCRSIfNeeded doesn't set it even though the geometry is in raster's CRS - geom = Functions.setSRID(geom, RasterAccessors.srid(raster)); - } + Pair pair = RasterUtils.setDefaultCRSAndTransform(raster, geom); + raster = pair.getLeft(); + geom = pair.getRight(); // checking if the raster contains the geometry if (!RasterPredicates.rsIntersects(raster, geom)) { diff --git a/common/src/main/java/org/apache/sedona/common/raster/PixelFunctions.java b/common/src/main/java/org/apache/sedona/common/raster/PixelFunctions.java index 9cc493abe0..a5486f47d2 100644 --- a/common/src/main/java/org/apache/sedona/common/raster/PixelFunctions.java +++ b/common/src/main/java/org/apache/sedona/common/raster/PixelFunctions.java @@ -18,6 +18,7 @@ */ package org.apache.sedona.common.raster; +import org.apache.commons.lang3.tuple.Pair; import org.apache.sedona.common.Functions; import org.apache.sedona.common.utils.RasterUtils; import org.geotools.coverage.grid.GridCoordinates2D; @@ -251,12 +252,18 @@ public static List values(GridCoverage2D rasterGeom, List geom RasterUtils.ensureBand(rasterGeom, band); // Check for invalid band index for (int i = 0; i < geometries.size(); i++) { - if (geometries.get(i) != null && RasterAccessors.srid(rasterGeom) != geometries.get(i).getSRID()) { - // implicitly converting roi geometry CRS to raster CRS - Geometry point = RasterUtils.convertCRSIfNeeded(geometries.get(i), rasterGeom.getCoordinateReferenceSystem()); - // have to set the SRID as RasterUtils.convertCRSIfNeeded doesn't set it even though the geometry is in raster's CRS - point = Functions.setSRID(point, RasterAccessors.srid(rasterGeom)); - geometries.set(i, point); + + if (geometries.get(i) == null) { + continue; + } + + + Pair pair = RasterUtils.setDefaultCRSAndTransform(rasterGeom, geometries.get(i)); + + geometries.set(i, pair.getRight()); + + if (i == 0) { + rasterGeom = pair.getLeft(); } } diff --git a/common/src/main/java/org/apache/sedona/common/raster/RasterBandEditors.java b/common/src/main/java/org/apache/sedona/common/raster/RasterBandEditors.java index a4ef2084ee..27fd410e05 100644 --- a/common/src/main/java/org/apache/sedona/common/raster/RasterBandEditors.java +++ b/common/src/main/java/org/apache/sedona/common/raster/RasterBandEditors.java @@ -19,6 +19,7 @@ package org.apache.sedona.common.raster; import org.apache.commons.lang3.ArrayUtils; +import org.apache.commons.lang3.tuple.Pair; import org.apache.sedona.common.Functions; import org.apache.sedona.common.utils.RasterUtils; import org.geotools.coverage.GridSampleDimension; @@ -190,12 +191,9 @@ public static GridCoverage2D clip(GridCoverage2D raster, int band, Geometry geom RasterUtils.ensureBand(raster, band); GridCoverage2D singleBandRaster = RasterBandAccessors.getBand(raster, new int[]{band}); - if(RasterAccessors.srid(raster) != geometry.getSRID()) { - // implicitly converting geometry CRS to raster CRS - geometry = RasterUtils.convertCRSIfNeeded(geometry, raster.getCoordinateReferenceSystem()); - // have to set the SRID as RasterUtils.convertCRSIfNeeded doesn't set it even though the geometry is in raster's CRS - geometry = Functions.setSRID(geometry, RasterAccessors.srid(raster)); - } + Pair pair = RasterUtils.setDefaultCRSAndTransform(singleBandRaster, geometry); + singleBandRaster = pair.getLeft(); + geometry = pair.getRight(); // Crop the raster // this will shrink the extent of the raster to the geometry diff --git a/common/src/main/java/org/apache/sedona/common/utils/RasterUtils.java b/common/src/main/java/org/apache/sedona/common/utils/RasterUtils.java index ec0103e727..d5d4f1dbf5 100644 --- a/common/src/main/java/org/apache/sedona/common/utils/RasterUtils.java +++ b/common/src/main/java/org/apache/sedona/common/utils/RasterUtils.java @@ -19,8 +19,11 @@ package org.apache.sedona.common.utils; import com.sun.media.imageioimpl.common.BogusColorSpace; +import org.apache.commons.lang3.tuple.Pair; +import org.apache.sedona.common.Functions; import org.apache.sedona.common.FunctionsGeoTools; import org.apache.sedona.common.raster.RasterAccessors; +import org.apache.sedona.common.raster.RasterEditors; import org.geotools.coverage.Category; import org.geotools.coverage.CoverageFactoryFinder; import org.geotools.coverage.GridSampleDimension; @@ -436,6 +439,38 @@ public static Geometry convertCRSIfNeeded(Geometry geometry, CoordinateReference return geometry; } + /** + * If the raster has a CRS, then it transforms the geom to the raster's CRS. + * If any of the inputs, raster or geom doesn't have a CRS, it defaults to 4326. + * @param raster + * @param geom + * @return + * @throws FactoryException + */ + public static Pair setDefaultCRSAndTransform(GridCoverage2D raster, Geometry geom) throws FactoryException { + int rasterSRID = RasterAccessors.srid(raster); + int geomSRID = Functions.getSRID(geom); + + if (rasterSRID == 0) { + raster = RasterEditors.setSrid(raster, 4326); + rasterSRID = RasterAccessors.srid(raster); + } + + if (geomSRID == 0) { + geom = Functions.setSRID(geom, 4326); + geomSRID = Functions.getSRID(geom); + } + + + if (rasterSRID != geomSRID) { + // implicitly converting roi geometry CRS to raster CRS + geom = convertCRSIfNeeded(geom, raster.getCoordinateReferenceSystem()); + // have to set the SRID as RasterUtils.convertCRSIfNeeded doesn't set it even though the geometry is in raster's CRS + geom = Functions.setSRID(geom, RasterAccessors.srid(raster)); + } + return Pair.of(raster, geom); + } + /** * Converts data types to data type codes * @param s pixel data type possible values {D, I, B, F, S, US}

diff --git a/common/src/test/java/org/apache/sedona/common/raster/FunctionEditorsTest.java b/common/src/test/java/org/apache/sedona/common/raster/FunctionEditorsTest.java index 46011868ce..b324e54ebd 100644 --- a/common/src/test/java/org/apache/sedona/common/raster/FunctionEditorsTest.java +++ b/common/src/test/java/org/apache/sedona/common/raster/FunctionEditorsTest.java @@ -20,6 +20,7 @@ import org.apache.sedona.common.Constructors; import org.apache.sedona.common.Functions; +import org.apache.sedona.common.FunctionsGeoTools; import org.apache.sedona.common.utils.RasterUtils; import org.geotools.coverage.grid.GridCoverage2D; import org.junit.Test; @@ -77,21 +78,17 @@ public void testSetValuesWithGeomInRaster() throws IOException, ParseException, @Test public void testSetValuesWithRasterNoSRID() throws IOException, ParseException, FactoryException, TransformException { GridCoverage2D raster = rasterFromGeoTiff(resourceFolder + "raster_geotiff_color/FAA_UTM18N_NAD83.tif"); - String polygon = "POLYGON ((236722 4204770, 243900 4204770, 243900 4197590, 236722 4197590, 236722 4204770))"; + String polygon = "POLYGON ((-77.9148 37.9545, -77.9123 37.8898, -77.9938 37.8878, -77.9964 37.9524, -77.9148 37.9545))"; Geometry geom = Constructors.geomFromWKT(polygon, 0); - raster = RasterEditors.setSrid(raster, 0); - - int rasterSRID = RasterAccessors.srid(raster); - assertEquals(0, rasterSRID); GridCoverage2D result = PixelFunctionEditors.setValues(raster, 1, geom, 10, false); - Geometry point = Constructors.geomFromWKT("POINT (243700 4197797)", 0); + Geometry point = Constructors.geomFromWKT("POINT (-77.9146 37.8916)", 0); double actual = PixelFunctions.value(result, point, 1); double expected = 10.0; assertEquals(expected, actual, 0d); - point = Constructors.geomFromWKT("POINT (240311 4202806)", 0); + point = Constructors.geomFromWKT("POINT (-77.9549 37.9357)", 0); actual = PixelFunctions.value(result, point, 1); assertEquals(expected, actual, 0d); } diff --git a/common/src/test/java/org/apache/sedona/common/raster/FunctionsTest.java b/common/src/test/java/org/apache/sedona/common/raster/FunctionsTest.java index fad42de0ec..2bf3d491bb 100644 --- a/common/src/test/java/org/apache/sedona/common/raster/FunctionsTest.java +++ b/common/src/test/java/org/apache/sedona/common/raster/FunctionsTest.java @@ -76,10 +76,10 @@ public void value() throws TransformException, FactoryException { public void valueImplicitTransform() throws TransformException, FactoryException, IOException, ParseException { GridCoverage2D raster = rasterFromGeoTiff(resourceFolder + "raster_geotiff_color/FAA_UTM18N_NAD83.tif"); - Geometry point = Constructors.geomFromWKT("POINT (-100 100)", 4326); + Geometry point = Constructors.geomFromWKT("POINT (-100 100)", 0); assertNull("Points outside of the envelope should return null.", PixelFunctions.value(raster, point, 1)); - point = Constructors.geomFromWKT("POINT (-77.9146 37.8916)", 4326); + point = Constructors.geomFromWKT("POINT (-77.9146 37.8916)", 0); Double value = PixelFunctions.value(raster, point, 1); assertNotNull(value); assertEquals(99.0d, value, 0.1d); @@ -343,7 +343,7 @@ public void valuesImplicitTransform() throws TransformException, FactoryExceptio GridCoverage2D raster = rasterFromGeoTiff(resourceFolder + "raster_geotiff_color/FAA_UTM18N_NAD83.tif"); List points = Arrays.asList(Constructors.geomFromWKT("POINT (-77.9146 37.8916)", 4326), - Constructors.geomFromWKT("POINT (-78.0059 37.9142)", 4326)); + Constructors.geomFromWKT("POINT (-78.0059 37.9142)", 0)); List values = PixelFunctions.values(raster, points, 1); assertEquals(2, values.size()); 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 8308f65139..61ebc13163 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 @@ -405,7 +405,7 @@ class rasteralgebraTest extends TestBaseScala with BeforeAndAfter with GivenWhen val df = dfFile.selectExpr("RS_FromGeoTiff(content) as raster", "ST_GeomFromWKT('POLYGON ((-8682522.873537656 4572703.890837922, -8673439.664183248 4572993.532747675, -8673155.57366801 4563873.2099182755, -8701890.325907696 4562931.7093397, -8682522.873537656 4572703.890837922))', 3857) as geom") val resultDf = df.selectExpr("RS_SetValues(raster, 1, geom, 10, false) as result") - var actual = resultDf.selectExpr("RS_Value(result, ST_GeomFromWKT('POINT (243700 4197797)', 26918))").first().get(0) + var actual = resultDf.selectExpr("RS_Value(result, ST_GeomFromWKT('POINT (-77.9146 37.8916)'))").first().get(0) val expected = 10.0 assertEquals(expected, actual) @@ -506,7 +506,7 @@ class rasteralgebraTest extends TestBaseScala with BeforeAndAfter with GivenWhen it("Passed RS_Value with raster") { val df = sparkSession.read.format("binaryFile").load(resourceFolder + "raster/test1.tiff") - val result = df.selectExpr("RS_Value(RS_FromGeoTiff(content), ST_Point(-13077301.685, 4002565.802))").first().getDouble(0) + val result = df.selectExpr("RS_Value(RS_FromGeoTiff(content), ST_SetSRID(ST_Point(-13077301.685, 4002565.802), 3857))").first().getDouble(0) assert(result == 255d) } @@ -529,7 +529,7 @@ class rasteralgebraTest extends TestBaseScala with BeforeAndAfter with GivenWhen it("Passed RS_Values with raster") { val df = sparkSession.read.format("binaryFile").load(resourceFolder + "raster/test1.tiff") - val result = df.selectExpr("RS_Values(RS_FromGeoTiff(content), array(ST_Point(-13077301.685, 4002565.802), null))").first().getList[Any](0) + val result = df.selectExpr("RS_Values(RS_FromGeoTiff(content), array(ST_SetSRID(ST_Point(-13077301.685, 4002565.802), 3857), null))").first().getList[Any](0) assert(result.size() == 2) assert(result.get(0) == 255d) assert(result.get(1) == null) @@ -542,7 +542,7 @@ class rasteralgebraTest extends TestBaseScala with BeforeAndAfter with GivenWhen val df = sparkSession.read.format("binaryFile").load(resourceFolder + "raster/test1.tiff") val points = sparkSession.createDataFrame(Seq(("POINT (-13077301.685 4002565.802)",1), ("POINT (0 0)",2))) .toDF("point", "id") - .withColumn("point", expr("ST_GeomFromText(point)")) + .withColumn("point", expr("ST_GeomFromText(point, 3857)")) .groupBy().agg(collect_list("point").alias("point")) val result = df.crossJoin(points).selectExpr("RS_Values(RS_FromGeoTiff(content), point, 1)").first().getList[Any](0) @@ -591,9 +591,9 @@ class rasteralgebraTest extends TestBaseScala with BeforeAndAfter with GivenWhen var actualValues = clippedDf.selectExpr( "RS_Values(clipped, " + - "Array(ST_GeomFromWKT('POINT(223802 4.21769e+06)'),ST_GeomFromWKT('POINT(224759 4.20453e+06)')," + - "ST_GeomFromWKT('POINT(237201 4.20429e+06)'),ST_GeomFromWKT('POINT(237919 4.20357e+06)')," + - "ST_GeomFromWKT('POINT(254668 4.21769e+06)')), 1)" + "Array(ST_GeomFromWKT('POINT(223802 4.21769e+06)', 26918),ST_GeomFromWKT('POINT(224759 4.20453e+06)', 26918)," + + "ST_GeomFromWKT('POINT(237201 4.20429e+06)', 26918),ST_GeomFromWKT('POINT(237919 4.20357e+06)', 26918)," + + "ST_GeomFromWKT('POINT(254668 4.21769e+06)', 26918)), 1)" ).first().get(0) var expectedValues = Seq(null, null, 0.0, 0.0, null) assertTrue(expectedValues.equals(actualValues)) @@ -602,7 +602,7 @@ class rasteralgebraTest extends TestBaseScala with BeforeAndAfter with GivenWhen it("Passed RS_Clip with raster") { val df = sparkSession.read.format("binaryFile").load(resourceFolder + "raster_geotiff_color/FAA_UTM18N_NAD83.tif") .selectExpr("RS_FromGeoTiff(content) as raster", - "ST_GeomFromWKT('POLYGON ((236722 4204770, 243900 4204770, 243900 4197590, 221170 4197590, 236722 4204770))') as geom") + "ST_GeomFromWKT('POLYGON ((236722 4204770, 243900 4204770, 243900 4197590, 221170 4197590, 236722 4204770))', 26918) as geom") val clippedDf = df.selectExpr("RS_Clip(raster, 1, geom, 200, false) as clipped") val clippedMetadata = df.selectExpr("RS_Metadata(raster)").first().getSeq(0).slice(0, 9) @@ -611,9 +611,9 @@ class rasteralgebraTest extends TestBaseScala with BeforeAndAfter with GivenWhen var actualValues = clippedDf.selectExpr( "RS_Values(clipped, " + - "Array(ST_GeomFromWKT('POINT(223802 4.21769e+06)'),ST_GeomFromWKT('POINT(224759 4.20453e+06)')," + - "ST_GeomFromWKT('POINT(237201 4.20429e+06)'),ST_GeomFromWKT('POINT(237919 4.20357e+06)')," + - "ST_GeomFromWKT('POINT(254668 4.21769e+06)')), 1)" + "Array(ST_GeomFromWKT('POINT(223802 4.21769e+06)', 26918),ST_GeomFromWKT('POINT(224759 4.20453e+06)', 26918)," + + "ST_GeomFromWKT('POINT(237201 4.20429e+06)', 26918),ST_GeomFromWKT('POINT(237919 4.20357e+06)', 26918)," + + "ST_GeomFromWKT('POINT(254668 4.21769e+06)', 26918)), 1)" ).first().get(0) var expectedValues = Seq(null, null, 0.0, 0.0, null) assertTrue(expectedValues.equals(actualValues)) @@ -621,9 +621,9 @@ class rasteralgebraTest extends TestBaseScala with BeforeAndAfter with GivenWhen val croppedDf = df.selectExpr("RS_Clip(raster, 1, geom, 200, false) as cropped") actualValues = croppedDf.selectExpr( "RS_Values(cropped, " + - "Array(ST_GeomFromWKT('POINT(236842 4.20465e+06)'),ST_GeomFromWKT('POINT(236961 4.20453e+06)')," + - "ST_GeomFromWKT('POINT(237201 4.20429e+06)'),ST_GeomFromWKT('POINT(237919 4.20357e+06)')," + - "ST_GeomFromWKT('POINT(223802 4.20465e+06)')), 1)" + "Array(ST_GeomFromWKT('POINT(236842 4.20465e+06)', 26918),ST_GeomFromWKT('POINT(236961 4.20453e+06)', 26918)," + + "ST_GeomFromWKT('POINT(237201 4.20429e+06)', 26918),ST_GeomFromWKT('POINT(237919 4.20357e+06)', 26918)," + + "ST_GeomFromWKT('POINT(223802 4.20465e+06)', 26918)), 1)" ).first().get(0) expectedValues = Seq(0.0, 0.0, 0.0, 0.0, null) assertTrue(expectedValues.equals(actualValues))