Skip to content

Commit

Permalink
feat: set default CRS as 4326
Browse files Browse the repository at this point in the history
  • Loading branch information
@furqaankhan
furqaankhan committed Dec 1, 2023
1 parent 3e09bbb commit 7d38d3e
Show file tree
Hide file tree
Showing 7 changed files with 77 additions and 42 deletions.
10 changes: 4 additions & 6 deletions common/src/main/java/org/apache/sedona/common/raster/PixelFunctionEditors.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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<GridCoverage2D, Geometry> pair = RasterUtils.setDefaultCRSAndTransform(raster, geom);
raster = pair.getLeft();
geom = pair.getRight();

// checking if the raster contains the geometry
if (!RasterPredicates.rsIntersects(raster, geom)) {
Expand Down
19 changes: 13 additions & 6 deletions common/src/main/java/org/apache/sedona/common/raster/PixelFunctions.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -251,12 +252,18 @@ public static List<Double> values(GridCoverage2D rasterGeom, List<Geometry> 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<GridCoverage2D, Geometry> pair = RasterUtils.setDefaultCRSAndTransform(rasterGeom, geometries.get(i));

geometries.set(i, pair.getRight());

if (i == 0) {
rasterGeom = pair.getLeft();
}
}

Expand Down
10 changes: 4 additions & 6 deletions common/src/main/java/org/apache/sedona/common/raster/RasterBandEditors.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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<GridCoverage2D, Geometry> 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
Expand Down
35 changes: 35 additions & 0 deletions common/src/main/java/org/apache/sedona/common/utils/RasterUtils.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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<GridCoverage2D, Geometry> 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} <br><br>
Expand Down
11 changes: 4 additions & 7 deletions common/src/test/java/org/apache/sedona/common/raster/FunctionEditorsTest.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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);
}
Expand Down
6 changes: 3 additions & 3 deletions common/src/test/java/org/apache/sedona/common/raster/FunctionsTest.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -343,7 +343,7 @@ public void valuesImplicitTransform() throws TransformException, FactoryExceptio
GridCoverage2D raster = rasterFromGeoTiff(resourceFolder + "raster_geotiff_color/FAA_UTM18N_NAD83.tif");

List<Geometry> 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<Double> values = PixelFunctions.values(raster, points, 1);
assertEquals(2, values.size());
Expand Down
28 changes: 14 additions & 14 deletions spark/common/src/test/scala/org/apache/sedona/sql/rasteralgebraTest.scala
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down Expand Up @@ -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)
}

Expand All @@ -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)
Expand All @@ -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)
Expand Down Expand Up @@ -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))
Expand All @@ -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)
Expand All @@ -611,19 +611,19 @@ 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))

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))
Expand Down

0 comments on commit 7d38d3e

Please sign in to comment.