This is an automated email from the ASF dual-hosted git repository.
jiayu pushed a commit to branch master
in repository https://gitbox.apache.org/repos/asf/sedona.git
The following commit(s) were added to refs/heads/master by this push:
new f5d0202dae [SEDONA-707] Refactor rasterization process and add
allTouched parameter for rasterization functions (#1788)
f5d0202dae is described below
commit f5d0202daeae5523b56924d5bd2c958c4ad90806
Author: Pranav Toggi <[email protected]>
AuthorDate: Mon Feb 24 18:18:57 2025 -0800
[SEDONA-707] Refactor rasterization process and add allTouched parameter
for rasterization functions (#1788)
* draft 1
* draft 2 - support search grid for sub-geoms in geom collections.
* fix useGeometryExtent in rasterizeLinestring
* profiling1
* profiling1
* profiling2
* profiling3
* profiling4
* profiling4
* profiling5
* fix rasterizeGeomExtent
* profiling6
* Fix rasterExtent bug; Update scala tests; Add rasterization specific tests
* Fix scala test; undo POM changes
* surface allTouched for RS_Clip;clean up
* spotless fix
* profiling...
* temp pom changes
* Update docs
* temp pom changes
* profiling 8
* test2
* test2
* test3
* refactor polygon rasterization
* fix lint
* Fix floating point issues; Handle multipolygons in rasterizePolygon
* Fix floating point issues; Handle multipolygons in rasterizePolygon
* Avoid slope calculation for vertical segments
* address comments
---
.../sedona/common/raster/PixelFunctionEditors.java | 103 ++--
.../sedona/common/raster/RasterBandAccessors.java | 86 ++-
.../sedona/common/raster/RasterBandEditors.java | 40 +-
.../sedona/common/raster/RasterConstructors.java | 155 ++---
.../apache/sedona/common/raster/Rasterization.java | 662 +++++++++++++++++++++
.../sedona/common/raster/FunctionEditorsTest.java | 74 ++-
.../common/raster/RasterBandAccessorsTest.java | 88 +--
.../common/raster/RasterBandEditorsTest.java | 11 +-
.../common/raster/RasterConstructorsTest.java | 174 +++---
docs/api/sql/Raster-operators.md | 80 ++-
docs/api/sql/Raster-writer.md | 27 +-
.../expressions/raster/PixelFunctionEditors.scala | 1 +
.../expressions/raster/RasterBandAccessors.scala | 1 +
.../expressions/raster/RasterBandEditors.scala | 1 +
.../expressions/raster/RasterConstructors.scala | 5 +-
.../expressions/raster/RasterFunctions.scala | 19 +-
.../scala/org/apache/sedona/sql/rasterIOTest.scala | 34 +-
.../org/apache/sedona/sql/rasteralgebraTest.scala | 81 +--
18 files changed, 1269 insertions(+), 373 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 e4c748e9d3..5a5db76cf7 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
@@ -26,8 +26,6 @@ import javax.media.jai.RasterFactory;
import org.apache.commons.lang3.tuple.Pair;
import org.apache.sedona.common.utils.RasterUtils;
import org.geotools.coverage.grid.GridCoverage2D;
-import org.geotools.geometry.DirectPosition2D;
-import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.Geometry;
import org.opengis.referencing.FactoryException;
import org.opengis.referencing.operation.TransformException;
@@ -121,13 +119,19 @@ public class PixelFunctionEditors {
* @param band Band of the raster to be edited
* @param geom Geometry region to update
* @param value Value to updated in the said region
+ * @param allTouched When set to true, sets value for all pixels touched by
geom
* @param keepNoData To keep no data value or not
* @return An updated raster
* @throws FactoryException
* @throws TransformException
*/
public static GridCoverage2D setValues(
- GridCoverage2D raster, int band, Geometry geom, double value, boolean
keepNoData)
+ GridCoverage2D raster,
+ int band,
+ Geometry geom,
+ double value,
+ boolean allTouched,
+ boolean keepNoData)
throws FactoryException, TransformException {
RasterUtils.ensureBand(raster, band);
@@ -148,9 +152,10 @@ public class PixelFunctionEditors {
if (keepNoData) {
noDataValue = RasterBandAccessors.getBandNoDataValue(raster, band);
- rasterizedGeom = RasterConstructors.asRaster(geom, raster, bandDataType,
value, noDataValue);
+ rasterizedGeom =
+ RasterConstructors.asRaster(geom, raster, bandDataType, allTouched,
value, noDataValue);
} else {
- rasterizedGeom = RasterConstructors.asRaster(geom, raster, bandDataType,
value);
+ rasterizedGeom = RasterConstructors.asRaster(geom, raster, bandDataType,
allTouched, value);
}
Raster rasterizedGeomData =
RasterUtils.getRaster(rasterizedGeom.getRenderedImage());
@@ -162,55 +167,36 @@ public class PixelFunctionEditors {
widthOriginalRaster = RasterAccessors.getWidth(raster);
WritableRaster rasterCopied = makeCopiedRaster(raster);
- String geometryType = geom.getGeometryType();
-
- // Taking a shortcut if the provided geometry is a Point or MultiPoint,
then taking the lat lon
- if (geometryType.equalsIgnoreCase(Geometry.TYPENAME_POINT)
- || geometryType.equalsIgnoreCase(Geometry.TYPENAME_MULTIPOINT)) {
- Coordinate[] coordinates = geom.getCoordinates();
- for (Coordinate pointCoordinate : coordinates) {
- int[] pointLocation =
- raster
- .getGridGeometry()
- .worldToGrid(new DirectPosition2D(pointCoordinate.x,
pointCoordinate.y))
- .getCoordinateValues();
- double[] pixel = rasterCopied.getPixel(pointLocation[0],
pointLocation[1], (double[]) null);
- pixel[band - 1] = rasterizedGeomData.getPixel(0, 0, (double[])
null)[0];
- rasterCopied.setPixel(pointLocation[0], pointLocation[1], pixel);
- }
- }
// Converting geometry to raster and then iterating through them
- else {
- int[] pixelLocation = RasterUtils.getGridCoordinatesFromWorld(raster,
colX, rowY);
- int x = pixelLocation[0], y = pixelLocation[1];
+ int[] pixelLocation = RasterUtils.getGridCoordinatesFromWorld(raster,
colX, rowY);
+ int x = pixelLocation[0], y = pixelLocation[1];
- // rasterX & rasterY are the starting pixels for the target raster
- int rasterX = Math.max(x, 0);
- int rasterY = Math.max(y, 0);
- // geometryX & geometryY are the starting pixels for the geometry raster
- int geometryX = rasterX - x;
- int geometryY = rasterY - y;
- // widthRegion & heightRegion are the size of the region to update
- int widthRegion = Math.min(widthGeometryRaster - geometryX,
widthOriginalRaster - rasterX);
- int heightRegion = Math.min(heightGeometryRaster - geometryY,
heightOriginalRaster - rasterY);
+ // rasterX & rasterY are the starting pixels for the target raster
+ int rasterX = Math.max(x, 0);
+ int rasterY = Math.max(y, 0);
+ // geometryX & geometryY are the starting pixels for the geometry raster
+ int geometryX = rasterX - x;
+ int geometryY = rasterY - y;
+ // widthRegion & heightRegion are the size of the region to update
+ int widthRegion = Math.min(widthGeometryRaster - geometryX,
widthOriginalRaster - rasterX);
+ int heightRegion = Math.min(heightGeometryRaster - geometryY,
heightOriginalRaster - rasterY);
- for (int j = 0; j < heightRegion; j++) {
- for (int i = 0; i < widthRegion; i++) {
- double[] pixel = rasterCopied.getPixel(rasterX + i, rasterY + j,
(double[]) null);
- // [0] as only one band in the rasterized Geometry
- double pixelNew =
- rasterizedGeomData.getPixel(geometryX + i, geometryY + j,
(double[]) null)[0];
- // skipping 0 from the rasterized geometry as
- if (pixelNew == 0
- || keepNoData && noDataValue != null && noDataValue ==
pixel[band - 1]) {
- continue;
- } else {
- pixel[band - 1] = pixelNew;
- }
- rasterCopied.setPixel(rasterX + i, rasterY + j, pixel);
+ for (int j = 0; j < heightRegion; j++) {
+ for (int i = 0; i < widthRegion; i++) {
+ double[] pixel = rasterCopied.getPixel(rasterX + i, rasterY + j,
(double[]) null);
+ // [0] as only one band in the rasterized Geometry
+ double pixelNew =
+ rasterizedGeomData.getPixel(geometryX + i, geometryY + j,
(double[]) null)[0];
+ // skipping 0 from the rasterized geometry as
+ if (pixelNew == 0 || keepNoData && noDataValue != null && noDataValue
== pixel[band - 1]) {
+ continue;
+ } else {
+ pixel[band - 1] = pixelNew;
}
+ rasterCopied.setPixel(rasterX + i, rasterY + j, pixel);
}
}
+
return RasterUtils.clone(
rasterCopied,
raster.getSampleDimensions(),
@@ -219,6 +205,25 @@ public class PixelFunctionEditors {
true); // keep metadata since this is essentially the same raster
}
+ /**
+ * Returns a raster by replacing the values of pixels in a specified
geometry region. It converts
+ * the Geometry to a raster using RS_AsRaster. A convenience function with
keepNoData as false.
+ *
+ * @param raster Input raster to be updated
+ * @param band Band of the raster to be edited
+ * @param geom Geometry region to update
+ * @param value Value to updated in the said region
+ * @param allTouched When set to true, sets value for all pixels touched by
geom
+ * @return An updated raster
+ * @throws FactoryException
+ * @throws TransformException
+ */
+ public static GridCoverage2D setValues(
+ GridCoverage2D raster, int band, Geometry geom, double value, boolean
allTouched)
+ throws FactoryException, TransformException {
+ return setValues(raster, band, geom, value, allTouched, false);
+ }
+
/**
* Returns a raster by replacing the values of pixels in a specified
geometry region. It converts
* the Geometry to a raster using RS_AsRaster. A convenience function with
keepNoData as false.
@@ -234,7 +239,7 @@ public class PixelFunctionEditors {
public static GridCoverage2D setValues(
GridCoverage2D raster, int band, Geometry geom, double value)
throws FactoryException, TransformException {
- return setValues(raster, band, geom, value, false);
+ return setValues(raster, band, geom, value, false, false);
}
/**
diff --git
a/common/src/main/java/org/apache/sedona/common/raster/RasterBandAccessors.java
b/common/src/main/java/org/apache/sedona/common/raster/RasterBandAccessors.java
index 4ed8e592eb..7f816d8b52 100644
---
a/common/src/main/java/org/apache/sedona/common/raster/RasterBandAccessors.java
+++
b/common/src/main/java/org/apache/sedona/common/raster/RasterBandAccessors.java
@@ -92,6 +92,7 @@ public class RasterBandAccessors {
* @param raster Raster to use for computing stats
* @param roi Geometry to define the region of interest
* @param band Band to be used for computation
+ * @param allTouched Include pixels touched by roi geometry
* @param excludeNoData Specifies whether to exclude no-data value or not
* @param lenient Return null if the raster and roi do not intersect when
set to true, otherwise
* will throw an exception
@@ -99,9 +100,14 @@ public class RasterBandAccessors {
* @throws FactoryException
*/
public static double[] getZonalStatsAll(
- GridCoverage2D raster, Geometry roi, int band, boolean excludeNoData,
boolean lenient)
+ GridCoverage2D raster,
+ Geometry roi,
+ int band,
+ boolean allTouched,
+ boolean excludeNoData,
+ boolean lenient)
throws FactoryException {
- List<Object> objects = getStatObjects(raster, roi, band, excludeNoData,
lenient);
+ List<Object> objects = getStatObjects(raster, roi, band, allTouched,
excludeNoData, lenient);
if (objects == null) {
return null;
}
@@ -128,14 +134,28 @@ public class RasterBandAccessors {
* @param raster Raster to use for computing stats
* @param roi Geometry to define the region of interest
* @param band Band to be used for computation
+ * @param allTouched Include pixels touched by roi geometry
* @param excludeNoData Specifies whether to exclude no-data value or not
* @return An array with all the stats for the region
* @throws FactoryException
*/
public static double[] getZonalStatsAll(
- GridCoverage2D raster, Geometry roi, int band, boolean excludeNoData)
+ GridCoverage2D raster, Geometry roi, int band, boolean allTouched,
boolean excludeNoData)
throws FactoryException {
- return getZonalStatsAll(raster, roi, band, excludeNoData, true);
+ return getZonalStatsAll(raster, roi, band, allTouched, excludeNoData,
true);
+ }
+
+ /**
+ * @param raster Raster to use for computing stats
+ * @param roi Geometry to define the region of interest
+ * @param band Band to be used for computation
+ * @param allTouched Include pixels touched by roi geometry
+ * @return An array with all the stats for the region, excludeNoData is set
to true
+ * @throws FactoryException
+ */
+ public static double[] getZonalStatsAll(
+ GridCoverage2D raster, Geometry roi, int band, boolean allTouched)
throws FactoryException {
+ return getZonalStatsAll(raster, roi, band, allTouched, true);
}
/**
@@ -147,7 +167,7 @@ public class RasterBandAccessors {
*/
public static double[] getZonalStatsAll(GridCoverage2D raster, Geometry roi,
int band)
throws FactoryException {
- return getZonalStatsAll(raster, roi, band, true);
+ return getZonalStatsAll(raster, roi, band, false);
}
/**
@@ -159,7 +179,7 @@ public class RasterBandAccessors {
*/
public static double[] getZonalStatsAll(GridCoverage2D raster, Geometry roi)
throws FactoryException {
- return getZonalStatsAll(raster, roi, 1, true);
+ return getZonalStatsAll(raster, roi, 1);
}
/**
@@ -167,6 +187,7 @@ public class RasterBandAccessors {
* @param roi Geometry to define the region of interest
* @param band Band to be used for computation
* @param statType Define the statistic to be computed
+ * @param allTouched Include pixels touched by roi geometry
* @param excludeNoData Specifies whether to exclude no-data value or not
* @param lenient Return null if the raster and roi do not intersect when
set to true, otherwise
* will throw an exception
@@ -179,10 +200,11 @@ public class RasterBandAccessors {
Geometry roi,
int band,
String statType,
+ boolean allTouched,
boolean excludeNoData,
boolean lenient)
throws FactoryException {
- List<Object> objects = getStatObjects(raster, roi, band, excludeNoData,
lenient);
+ List<Object> objects = getStatObjects(raster, roi, band, allTouched,
excludeNoData, lenient);
if (objects == null) {
return null;
}
@@ -217,10 +239,42 @@ public class RasterBandAccessors {
}
}
+ /**
+ * @param raster Raster to use for computing stats
+ * @param roi Geometry to define the region of interest
+ * @param band Band to be used for computation
+ * @param statType Define the statistic to be computed
+ * @param allTouched Include pixels touched by roi geometry
+ * @param excludeNoData Specifies whether to exclude no-data value or not
+ * @return A double precision floating point number representing the
requested statistic
+ * calculated over the specified region.
+ * @throws FactoryException
+ */
public static Double getZonalStats(
- GridCoverage2D raster, Geometry roi, int band, String statType, boolean
excludeNoData)
+ GridCoverage2D raster,
+ Geometry roi,
+ int band,
+ String statType,
+ boolean allTouched,
+ boolean excludeNoData)
throws FactoryException {
- return getZonalStats(raster, roi, band, statType, excludeNoData, true);
+ return getZonalStats(raster, roi, band, statType, allTouched,
excludeNoData, true);
+ }
+
+ /**
+ * @param raster Raster to use for computing stats
+ * @param roi Geometry to define the region of interest
+ * @param band Band to be used for computation
+ * @param statType Define the statistic to be computed
+ * @param allTouched Include pixels touched by roi geometry
+ * @return A double precision floating point number representing the
requested statistic
+ * calculated over the specified region. The excludeNoData is set to
true.
+ * @throws FactoryException
+ */
+ public static Double getZonalStats(
+ GridCoverage2D raster, Geometry roi, int band, String statType, boolean
allTouched)
+ throws FactoryException {
+ return getZonalStats(raster, roi, band, statType, allTouched, true);
}
/**
@@ -234,7 +288,7 @@ public class RasterBandAccessors {
*/
public static Double getZonalStats(GridCoverage2D raster, Geometry roi, int
band, String statType)
throws FactoryException {
- return getZonalStats(raster, roi, band, statType, true);
+ return getZonalStats(raster, roi, band, statType, false);
}
/**
@@ -248,7 +302,7 @@ public class RasterBandAccessors {
*/
public static Double getZonalStats(GridCoverage2D raster, Geometry roi,
String statType)
throws FactoryException {
- return getZonalStats(raster, roi, 1, statType, true);
+ return getZonalStats(raster, roi, 1, statType, false);
}
/**
@@ -267,6 +321,7 @@ public class RasterBandAccessors {
* @param raster Raster to use for computing stats
* @param roi Geometry to define the region of interest
* @param band Band to be used for computation
+ * @param allTouched Include pixels touched by roi geometry
* @param excludeNoData Specifies whether to exclude no-data value or not
* @param lenient Return null if the raster and roi do not intersect when
set to true, otherwise
* will throw an exception
@@ -274,7 +329,12 @@ public class RasterBandAccessors {
* @throws FactoryException
*/
private static List<Object> getStatObjects(
- GridCoverage2D raster, Geometry roi, int band, boolean excludeNoData,
boolean lenient)
+ GridCoverage2D raster,
+ Geometry roi,
+ int band,
+ boolean allTouched,
+ boolean excludeNoData,
+ boolean lenient)
throws FactoryException {
RasterUtils.ensureBand(raster, band);
@@ -301,7 +361,7 @@ public class RasterBandAccessors {
Double noDataValue = RasterBandAccessors.getBandNoDataValue(raster, band);
// Adding an arbitrary value '150' for the pixels that are under the
geometry.
GridCoverage2D rasterizedGeom =
- RasterConstructors.asRasterWithRasterExtent(roi, raster, datatype,
150, null);
+ RasterConstructors.asRasterWithRasterExtent(roi, raster, datatype,
allTouched, 150, null);
Raster rasterziedData =
RasterUtils.getRaster(rasterizedGeom.getRenderedImage());
int width = RasterAccessors.getWidth(rasterizedGeom),
height = RasterAccessors.getHeight(rasterizedGeom);
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 781ef7bed3..a05872722b 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
@@ -272,6 +272,7 @@ public class RasterBandEditors {
* @param raster Raster to clip
* @param band Band number to perform clipping
* @param geometry Specify ROI
+ * @param allTouched Include all pixels touched by roi geometry
* @param noDataValue no-Data value for empty cells
* @param crop Specifies to keep the original extent or not
* @param lenient Return null if the raster and geometry do not intersect
when set to true,
@@ -282,6 +283,7 @@ public class RasterBandEditors {
GridCoverage2D raster,
int band,
Geometry geometry,
+ boolean allTouched,
double noDataValue,
boolean crop,
boolean lenient)
@@ -343,7 +345,8 @@ public class RasterBandEditors {
}
// rasterize the geometry to iterate over the clipped raster
- GridCoverage2D rasterized = RasterConstructors.asRaster(geometry,
raster, bandType, 150);
+ GridCoverage2D rasterized =
+ RasterConstructors.asRaster(geometry, raster, bandType, allTouched,
150);
Raster rasterizedData =
RasterUtils.getRaster(rasterized.getRenderedImage());
double[] metadataRasterized = RasterAccessors.metadata(rasterized);
int widthRasterized = (int) metadataRasterized[2],
@@ -406,14 +409,20 @@ public class RasterBandEditors {
* @param raster Raster to clip
* @param band Band number to perform clipping
* @param geometry Specify ROI
+ * @param allTouched Include all pixels touched by roi geometry
* @param noDataValue no-Data value for empty cells
* @param crop Specifies to keep the original extent or not
* @return A clip Raster with defined ROI by the geometry
*/
public static GridCoverage2D clip(
- GridCoverage2D raster, int band, Geometry geometry, double noDataValue,
boolean crop)
+ GridCoverage2D raster,
+ int band,
+ Geometry geometry,
+ boolean allTouched,
+ double noDataValue,
+ boolean crop)
throws FactoryException, TransformException {
- return clip(raster, band, geometry, noDataValue, crop, true);
+ return clip(raster, band, geometry, allTouched, noDataValue, crop, true);
}
/**
@@ -422,13 +431,14 @@ public class RasterBandEditors {
* @param raster Raster to clip
* @param band Band number to perform clipping
* @param geometry Specify ROI
+ * @param allTouched Include all pixels touched by roi geometry
* @param noDataValue no-Data value for empty cells
* @return A clip Raster with defined ROI by the geometry
*/
public static GridCoverage2D clip(
- GridCoverage2D raster, int band, Geometry geometry, double noDataValue)
+ GridCoverage2D raster, int band, Geometry geometry, boolean allTouched,
double noDataValue)
throws FactoryException, TransformException {
- return clip(raster, band, geometry, noDataValue, true);
+ return clip(raster, band, geometry, allTouched, noDataValue, true);
}
/**
@@ -438,9 +448,11 @@ public class RasterBandEditors {
* @param raster Raster to clip
* @param band Band number to perform clipping
* @param geometry Specify ROI
+ * @param allTouched Include all pixels touched by roi geometry
* @return A clip Raster with defined ROI by the geometry
*/
- public static GridCoverage2D clip(GridCoverage2D raster, int band, Geometry
geometry)
+ public static GridCoverage2D clip(
+ GridCoverage2D raster, int band, Geometry geometry, boolean allTouched)
throws FactoryException, TransformException {
boolean isDataTypeIntegral =
RasterUtils.isDataTypeIntegral(
@@ -452,6 +464,20 @@ public class RasterBandEditors {
} else {
noDataValue = Double.MIN_VALUE;
}
- return clip(raster, band, geometry, noDataValue, true);
+ return clip(raster, band, geometry, allTouched, noDataValue, true);
+ }
+
+ /**
+ * Return a clipped raster with the specified ROI by the geometry. No-data
value will be taken as
+ * the lowest possible value for the data type and crop will be `true`.
+ *
+ * @param raster Raster to clip
+ * @param band Band number to perform clipping
+ * @param geometry Specify ROI
+ * @return A clip Raster with defined ROI by the geometry
+ */
+ public static GridCoverage2D clip(GridCoverage2D raster, int band, Geometry
geometry)
+ throws FactoryException, TransformException {
+ return clip(raster, band, geometry, false);
}
}
diff --git
a/common/src/main/java/org/apache/sedona/common/raster/RasterConstructors.java
b/common/src/main/java/org/apache/sedona/common/raster/RasterConstructors.java
index 9e16be4d43..c43596c6dd 100644
---
a/common/src/main/java/org/apache/sedona/common/raster/RasterConstructors.java
+++
b/common/src/main/java/org/apache/sedona/common/raster/RasterConstructors.java
@@ -18,12 +18,11 @@
*/
package org.apache.sedona.common.raster;
-import java.awt.Rectangle;
+import java.awt.*;
import java.awt.image.Raster;
import java.awt.image.RenderedImage;
import java.awt.image.WritableRaster;
import java.io.IOException;
-import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
import java.util.Map;
@@ -42,10 +41,6 @@ import org.geotools.feature.simple.SimpleFeatureBuilder;
import org.geotools.feature.simple.SimpleFeatureTypeBuilder;
import org.geotools.gce.arcgrid.ArcGridReader;
import org.geotools.gce.geotiff.GeoTiffReader;
-import org.geotools.geometry.Envelope2D;
-import org.geotools.geometry.jts.JTS;
-import org.geotools.geometry.jts.ReferencedEnvelope;
-import org.geotools.process.vector.VectorToRasterProcess;
import org.geotools.referencing.CRS;
import org.geotools.referencing.crs.DefaultEngineeringCRS;
import org.geotools.referencing.operation.transform.AffineTransform2D;
@@ -106,6 +101,7 @@ public class RasterConstructors {
* @param geom The geometry to convert
* @param raster The reference raster
* @param pixelType The data type of pixel/cell of resultant raster
+ * @param allTouched When set to true, rasterizes all pixels touched by geom
* @param value The value of the pixel of the resultant raster
* @param noDataValue The noDataValue of the resultant raster
* @param useGeometryExtent The way to generate extent of the resultant
raster. Use the extent of
@@ -117,21 +113,31 @@ public class RasterConstructors {
Geometry geom,
GridCoverage2D raster,
String pixelType,
+ boolean allTouched,
double value,
Double noDataValue,
boolean useGeometryExtent)
throws FactoryException {
+
List<Object> objects =
- rasterization(geom, raster, pixelType, value, noDataValue,
useGeometryExtent);
+ Rasterization.rasterize(geom, raster, pixelType, value,
useGeometryExtent, allTouched);
+
WritableRaster writableRaster = (WritableRaster) objects.get(0);
GridCoverage2D rasterized = (GridCoverage2D) objects.get(1);
- return RasterUtils.clone(
- writableRaster,
- rasterized.getSampleDimensions(),
- rasterized,
- noDataValue,
- false); // no need to original raster metadata since this is a new
raster.
+ GridCoverage2D resultRaster =
+ RasterUtils.clone(
+ writableRaster,
+ rasterized.getSampleDimensions(),
+ rasterized,
+ noDataValue,
+ false); // no need to original raster metadata since this is a new
raster.
+
+ if (noDataValue != null) {
+ resultRaster = RasterBandEditors.setBandNoDataValue(resultRaster, 1,
noDataValue);
+ }
+
+ return resultRaster;
}
/**
@@ -141,15 +147,22 @@ public class RasterConstructors {
* @param geom The geometry to convert
* @param raster The reference raster
* @param pixelType The data type of pixel/cell of resultant raster
+ * @param allTouched When set to true, rasterizes all pixels touched by geom
* @param value The value of the pixel of the resultant raster
* @param noDataValue The noDataValue of the resultant raster
* @return Rasterized Geometry
* @throws FactoryException
*/
public static GridCoverage2D asRaster(
- Geometry geom, GridCoverage2D raster, String pixelType, double value,
Double noDataValue)
+ Geometry geom,
+ GridCoverage2D raster,
+ String pixelType,
+ boolean allTouched,
+ double value,
+ Double noDataValue)
throws FactoryException {
- return asRaster(geom, raster, pixelType, value, noDataValue, true);
+ GridCoverage2D result = asRaster(geom, raster, pixelType, allTouched,
value, noDataValue, true);
+ return result;
}
/**
@@ -164,7 +177,7 @@ public class RasterConstructors {
*/
public static GridCoverage2D asRaster(Geometry geom, GridCoverage2D raster,
String pixelType)
throws FactoryException {
- return asRaster(geom, raster, pixelType, 1, null);
+ return asRaster(geom, raster, pixelType, false, 1, null);
}
/**
@@ -174,97 +187,32 @@ public class RasterConstructors {
* @param geom The geometry to convert
* @param raster The reference raster
* @param pixelType The data type of pixel/cell of resultant raster.
- * @param value The value of the pixel of the resultant raster
+ * @param allTouched When set to true, rasterizes all pixels touched by geom
* @return Rasterized Geometry
* @throws FactoryException
*/
public static GridCoverage2D asRaster(
- Geometry geom, GridCoverage2D raster, String pixelType, double value)
+ Geometry geom, GridCoverage2D raster, String pixelType, boolean
allTouched)
throws FactoryException {
- return asRaster(geom, raster, pixelType, value, null);
+ return asRaster(geom, raster, pixelType, allTouched, 1, null);
}
- private static List<Object> rasterization(
- Geometry geom,
- GridCoverage2D raster,
- String pixelType,
- double value,
- Double noDataValue,
- boolean useGeometryExtent)
+ /**
+ * Returns a raster that is converted from the geometry provided. A
convenience function for
+ * asRaster.
+ *
+ * @param geom The geometry to convert
+ * @param raster The reference raster
+ * @param pixelType The data type of pixel/cell of resultant raster.
+ * @param allTouched When set to true, rasterizes all pixels touched by geom
+ * @param value The value of the pixel of the resultant raster
+ * @return Rasterized Geometry
+ * @throws FactoryException
+ */
+ public static GridCoverage2D asRaster(
+ Geometry geom, GridCoverage2D raster, String pixelType, boolean
allTouched, double value)
throws FactoryException {
- DefaultFeatureCollection featureCollection =
- getFeatureCollection(geom, raster.getCoordinateReferenceSystem());
-
- double[] metadata = RasterAccessors.metadata(raster);
- // The current implementation doesn't support rasters with properties below
- // It is not a problem as most rasters don't have these properties
- // ScaleX < 0
- if (metadata[4] < 0) {
- throw new IllegalArgumentException(
- String.format("ScaleX %f of the raster is negative, it should be
positive", metadata[4]));
- }
- // ScaleY > 0
- if (metadata[5] > 0) {
- throw new IllegalArgumentException(
- String.format(
- "ScaleY %f of the raster is positive. It should be negative.",
metadata[5]));
- }
- // SkewX should be zero
- if (metadata[6] != 0) {
- throw new IllegalArgumentException(
- String.format("SkewX %d of the raster is not zero.", metadata[6]));
- }
- // SkewY should be zero
- if (metadata[7] != 0) {
- throw new IllegalArgumentException(
- String.format("SkewY %d of the raster is not zero.", metadata[7]));
- }
-
- Envelope2D bound = null;
-
- int width, height;
- if (useGeometryExtent) {
- bound =
- JTS.getEnvelope2D(geom.getEnvelopeInternal(),
raster.getCoordinateReferenceSystem2D());
- double scaleX = Math.abs(metadata[4]), scaleY = Math.abs(metadata[5]);
- width = Math.max((int) Math.ceil(bound.getWidth() / scaleX), 1);
- height = Math.max((int) Math.ceil(bound.getHeight() / scaleY), 1);
- bound =
- new Envelope2D(
- bound.getCoordinateReferenceSystem(),
- bound.getMinX(),
- bound.getMinY(),
- width * scaleX,
- height * scaleY);
- } else {
- ReferencedEnvelope envelope =
- ReferencedEnvelope.create(raster.getEnvelope(),
raster.getCoordinateReferenceSystem());
- bound = JTS.getEnvelope2D(envelope,
raster.getCoordinateReferenceSystem2D());
- GridEnvelope2D gridRange = raster.getGridGeometry().getGridRange2D();
- width = gridRange.width;
- height = gridRange.height;
- }
-
- VectorToRasterProcess rasterProcess = new VectorToRasterProcess();
- GridCoverage2D rasterized =
- rasterProcess.execute(
- featureCollection, width, height, "value", Double.toString(value),
bound, null);
- if (noDataValue != null) {
- rasterized = RasterBandEditors.setBandNoDataValue(rasterized, 1,
noDataValue);
- }
- WritableRaster writableRaster =
- RasterFactory.createBandedRaster(
- RasterUtils.getDataTypeCode(pixelType), width, height, 1, null);
- double[] samples =
- RasterUtils.getRaster(rasterized.getRenderedImage())
- .getSamples(0, 0, width, height, 0, (double[]) null);
- writableRaster.setSamples(0, 0, width, height, 0, samples);
-
- List<Object> objects = new ArrayList<>();
- objects.add(writableRaster);
- objects.add(rasterized);
-
- return objects;
+ return asRaster(geom, raster, pixelType, allTouched, value, null);
}
/**
@@ -280,9 +228,14 @@ public class RasterConstructors {
* @throws FactoryException
*/
public static GridCoverage2D asRasterWithRasterExtent(
- Geometry geom, GridCoverage2D raster, String pixelType, double value,
Double noDataValue)
+ Geometry geom,
+ GridCoverage2D raster,
+ String pixelType,
+ boolean allTouched,
+ double value,
+ Double noDataValue)
throws FactoryException {
- return asRaster(geom, raster, pixelType, value, noDataValue, false);
+ return asRaster(geom, raster, pixelType, allTouched, value, noDataValue,
false);
}
public static DefaultFeatureCollection getFeatureCollection(
diff --git
a/common/src/main/java/org/apache/sedona/common/raster/Rasterization.java
b/common/src/main/java/org/apache/sedona/common/raster/Rasterization.java
new file mode 100644
index 0000000000..652e822e92
--- /dev/null
+++ b/common/src/main/java/org/apache/sedona/common/raster/Rasterization.java
@@ -0,0 +1,662 @@
+/*
+ * 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 java.awt.image.WritableRaster;
+import java.math.BigDecimal;
+import java.math.RoundingMode;
+import java.util.*;
+import javax.media.jai.RasterFactory;
+import org.apache.sedona.common.Functions;
+import org.apache.sedona.common.utils.RasterUtils;
+import org.geotools.coverage.grid.GridCoverage2D;
+import org.geotools.coverage.grid.GridCoverageFactory;
+import org.geotools.geometry.Envelope2D;
+import org.geotools.geometry.jts.JTS;
+import org.locationtech.jts.geom.*;
+import org.opengis.geometry.BoundingBox;
+import org.opengis.referencing.FactoryException;
+
+public class Rasterization {
+ protected static List<Object> rasterize(
+ Geometry geom,
+ GridCoverage2D raster,
+ String pixelType,
+ double value,
+ boolean useGeometryExtent,
+ boolean allTouched)
+ throws FactoryException {
+
+ // Validate the input geometry and raster metadata
+ double[] metadata = RasterAccessors.metadata(raster);
+ validateRasterMetadata(metadata);
+ if (!RasterPredicates.rsIntersects(raster, geom)) {
+ throw new IllegalArgumentException("Geometry does not intersect
Raster.");
+ }
+
+ Envelope2D rasterExtent = raster.getEnvelope2D();
+
+ Envelope2D geomExtent = rasterizeGeomExtent(geom, raster, metadata,
allTouched);
+
+ RasterizationParams params =
+ calculateRasterizationParams(raster, useGeometryExtent, metadata,
geomExtent, pixelType);
+
+ rasterizeGeometry(raster, metadata, geom, params, geomExtent, value,
allTouched);
+
+ // Create a GridCoverage2D for the rasterized result
+ GridCoverageFactory coverageFactory = new GridCoverageFactory();
+ GridCoverage2D rasterized;
+
+ if (useGeometryExtent) {
+ rasterized = coverageFactory.create("rasterized", params.writableRaster,
geomExtent);
+ } else {
+ rasterized = coverageFactory.create("rasterized", params.writableRaster,
rasterExtent);
+ }
+
+ // Return results compatible with the original function
+ List<Object> objects = new ArrayList<>();
+ objects.add(params.writableRaster);
+ objects.add(rasterized);
+
+ return objects;
+ }
+
+ private static void rasterizeGeometry(
+ GridCoverage2D raster,
+ double[] metadata,
+ Geometry geom,
+ RasterizationParams params,
+ Envelope2D geomExtent,
+ double value,
+ boolean allTouched)
+ throws FactoryException {
+
+ switch (geom.getGeometryType()) {
+ case "GeometryCollection":
+ case "MultiPolygon":
+ case "MultiPoint":
+ rasterizeGeometryCollection(raster, metadata, geom, params, value,
allTouched);
+ break;
+ case "Point":
+ rasterizePoint(geom, params, geomExtent, value);
+ break;
+ case "LineString":
+ case "MultiLineString":
+ case "LinearRing":
+ rasterizeLineString(geom, params, value, geomExtent);
+ break;
+ default:
+ rasterizePolygon(geom, params, geomExtent, value, allTouched);
+ break;
+ }
+ }
+
+ private static void rasterizeGeometryCollection(
+ GridCoverage2D raster,
+ double[] metadata,
+ Geometry geom,
+ RasterizationParams params,
+ double value,
+ boolean allTouched)
+ throws FactoryException {
+
+ for (int i = 0; i < geom.getNumGeometries(); i++) {
+ Geometry subGeom = geom.getGeometryN(i);
+ Envelope2D subGeomExtent = rasterizeGeomExtent(subGeom, raster,
metadata, allTouched);
+ rasterizeGeometry(raster, metadata, subGeom, params, subGeomExtent,
value, allTouched);
+ }
+ }
+
+ private static void rasterizePoint(
+ Geometry geom, RasterizationParams params, Envelope2D geomExtent, double
value) {
+
+ int startX = (int) Math.round(((geomExtent.getMinX() - params.upperLeftX)
/ params.scaleX));
+ int startY = (int) Math.round(((geomExtent.getMinY() - params.upperLeftY)
/ params.scaleY));
+ int x = startX;
+ int y = startY;
+
+ for (double worldY = geomExtent.getMinY();
+ worldY < geomExtent.getMaxY();
+ worldY += params.scaleY, y++) {
+ x = startX;
+ for (double worldX = geomExtent.getMinX();
+ worldX < geomExtent.getMaxX();
+ worldX += params.scaleX, x++) {
+
+ // Flip y-axis (since raster Y starts from top-left)
+ int yIndex = -y - 1;
+
+ // Create envelope for this pixel
+ double cellMaxX = worldX + params.scaleX;
+ double cellMaxY = worldY + params.scaleY;
+
+ boolean intersects =
+ geom.intersects(JTS.toGeometry(new Envelope(worldX, cellMaxX,
worldY, cellMaxY)));
+
+ if (intersects) {
+ params.writableRaster.setSample(x, yIndex, 0, value);
+ }
+ }
+ }
+ }
+
+ private static void rasterizeLineString(
+ Geometry geom, RasterizationParams params, double value, Envelope2D
geomExtent) {
+ for (int i = 0; i < geom.getNumGeometries(); i++) {
+ LineString line = (LineString) geom.getGeometryN(i);
+ Coordinate[] coords = line.getCoordinates();
+
+ for (int j = 0; j < coords.length - 1; j++) {
+ // Extract start and end points for the segment
+ LineSegment clippedSegment =
+ clipSegmentToRasterBounds(coords[j], coords[j + 1], geomExtent);
+
+ Coordinate start;
+ Coordinate end;
+ if (clippedSegment != null) {
+ start = new Coordinate(clippedSegment.p0.x, clippedSegment.p0.y);
+ end = new Coordinate(clippedSegment.p1.x, clippedSegment.p1.y);
+ } else {
+ continue; // Skip case where segment is completely outside geomExtent
+ }
+
+ double x0 = (start.x - params.upperLeftX) / params.scaleX;
+ double y0 = (params.upperLeftY - start.y) / params.scaleY;
+ double x1 = (end.x - params.upperLeftX) / params.scaleX;
+ double y1 = (params.upperLeftY - end.y) / params.scaleY;
+
+ // Apply Bresenham for this segment
+ drawLineBresenham(params, x0, y0, x1, y1, value, 0.2);
+ }
+ }
+ }
+
+ // Modified Bresenham with Fractional Steps
+ private static void drawLineBresenham(
+ RasterizationParams params,
+ double x0,
+ double y0,
+ double x1,
+ double y1,
+ double value,
+ double stepSize) {
+
+ double dx = x1 - x0;
+ double dy = y1 - y0;
+
+ // Compute the number of steps based on the larger of dx or dy
+ double distance = Math.sqrt(dx * dx + dy * dy);
+ int steps = (int) Math.ceil(distance / stepSize);
+
+ // Calculate the step increment for each axis
+ double stepX = dx / steps;
+ double stepY = dy / steps;
+
+ // Start stepping through the line
+ double x = x0;
+ double y = y0;
+
+ for (int i = 0; i <= steps; i++) {
+ int rasterX = (int) (Math.floor(x));
+ int rasterY = (int) (Math.floor(y));
+
+ // Only write if within raster bounds
+ if (rasterX >= 0
+ && rasterX < params.writableRaster.getWidth()
+ && rasterY >= 0
+ && rasterY < params.writableRaster.getHeight()) {
+ params.writableRaster.setSample(rasterX, rasterY, 0, value);
+ }
+
+ // Increment by fractional steps
+ x += stepX;
+ y += stepY;
+ }
+ }
+
+ private static LineSegment clipSegmentToRasterBounds(
+ Coordinate p1, Coordinate p2, Envelope2D geomExtent) {
+ double minX = geomExtent.getMinX();
+ double maxX = geomExtent.getMaxX();
+ double minY = geomExtent.getMinY();
+ double maxY = geomExtent.getMaxY();
+
+ double x1 = p1.x, y1 = p1.y;
+ double x2 = p2.x, y2 = p2.y;
+
+ boolean p1Inside = isInsideBounds(x1, y1, minX, maxX, minY, maxY);
+ boolean p2Inside = isInsideBounds(x2, y2, minX, maxX, minY, maxY);
+
+ if (p1Inside && p2Inside) {
+ // Both points inside: no clipping needed
+ return new LineSegment(p1, p2);
+ }
+
+ if (!p1Inside && !p2Inside) {
+ // Both points outside: no clipping needed
+ return null;
+ }
+
+ double dx = x2 - x1;
+ double dy = y2 - y1;
+
+ // Clip using parametric line equation
+ double[] tValues = {0, 1}; // Stores valid segment proportions
+
+ // Clip against minX and maxX
+ if (dx != 0) {
+ double tMin = (minX - x1) / dx;
+ double tMax = (maxX - x1) / dx;
+ if (dx < 0) {
+ double temp = tMin;
+ tMin = tMax;
+ tMax = temp;
+ }
+ tValues[0] = Math.max(tValues[0], tMin);
+ tValues[1] = Math.min(tValues[1], tMax);
+ }
+
+ // Clip against minY and maxY
+ if (dy != 0) {
+ double tMin = (minY - y1) / dy;
+ double tMax = (maxY - y1) / dy;
+ if (dy < 0) {
+ double temp = tMin;
+ tMin = tMax;
+ tMax = temp;
+ }
+ tValues[0] = Math.max(tValues[0], tMin);
+ tValues[1] = Math.min(tValues[1], tMax);
+ }
+
+ // If tValues are invalid (segment is outside), return null
+ if (tValues[0] > tValues[1]) {
+ return null; // No valid clipped segment
+ }
+
+ // Compute new clipped endpoints
+ Coordinate newP1 = new Coordinate(x1 + tValues[0] * dx, y1 + tValues[0] *
dy);
+ Coordinate newP2 = new Coordinate(x1 + tValues[1] * dx, y1 + tValues[1] *
dy);
+
+ return new LineSegment(newP1, newP2);
+ }
+
+ // Helper function to check if a point is inside the bounding box
+ private static boolean isInsideBounds(
+ double x, double y, double minX, double maxX, double minY, double maxY) {
+ return x >= minX && x <= maxX && y >= minY && y <= maxY;
+ }
+
+ private static Envelope2D rasterizeGeomExtent(
+ Geometry geom, GridCoverage2D raster, double[] metadata, boolean
allTouched) {
+
+ if (Objects.equals(geom.getGeometryType(), "MultiLineString")) {
+ allTouched = true;
+ }
+ if (Objects.equals(geom.getGeometryType(), "MultiPoint")) {
+ allTouched = true;
+ }
+
+ Envelope2D rasterExtent =
+ JTS.getEnvelope2D(geom.getEnvelopeInternal(),
raster.getCoordinateReferenceSystem2D());
+
+ // Using BigDecimal to avoid floating point errors
+ BigDecimal upperLeftX = BigDecimal.valueOf(metadata[0]);
+ BigDecimal upperLeftY = BigDecimal.valueOf(metadata[1]);
+ BigDecimal scaleX = BigDecimal.valueOf(metadata[4]);
+ BigDecimal scaleY = BigDecimal.valueOf(metadata[5]);
+
+ // Compute the aligned min/max values
+ double alignedMinX =
+ (scaleX
+ .multiply(
+ BigDecimal.valueOf(
+ Math.floor((rasterExtent.getMinX() + metadata[0]) /
metadata[4])))
+ .subtract(upperLeftX))
+ .doubleValue();
+
+ double alignedMinY =
+ (scaleY
+ .multiply(
+ BigDecimal.valueOf(
+ Math.ceil((rasterExtent.getMinY() + metadata[1]) /
metadata[5])))
+ .subtract(upperLeftY))
+ .doubleValue();
+
+ double alignedMaxX =
+ (scaleX
+ .multiply(
+ BigDecimal.valueOf(
+ Math.ceil((rasterExtent.getMaxX() + metadata[0]) /
metadata[4])))
+ .subtract(upperLeftX))
+ .doubleValue();
+
+ double alignedMaxY =
+ (scaleY
+ .multiply(
+ BigDecimal.valueOf(
+ Math.floor((rasterExtent.getMaxY() + metadata[1]) /
metadata[5])))
+ .subtract(upperLeftY))
+ .doubleValue();
+
+ // For points and LineStrings at intersection of 2 or more pixels,
+ // extend search grid by 1 pixel in each direction
+ if (alignedMaxX == alignedMinX) {
+ alignedMinX -= metadata[4];
+ alignedMaxX += metadata[4];
+ }
+ if (alignedMaxY == alignedMinY) {
+ alignedMinY += metadata[5];
+ alignedMaxY -= metadata[5];
+ }
+
+ // Get the extent of the original raster
+ double originalMinX = raster.getEnvelope().getMinimum(0);
+ double originalMinY = raster.getEnvelope().getMinimum(1);
+ double originalMaxX = raster.getEnvelope().getMaximum(0);
+ double originalMaxY = raster.getEnvelope().getMaximum(1);
+
+ // Extend the aligned extent by 1 pixel if allTouched is true and if any
geometry extent meets
+ // the aligned extent
+ if (allTouched) {
+ alignedMinX -= (rasterExtent.getMinX() == alignedMinX) ? metadata[4] : 0;
+ alignedMinY += (rasterExtent.getMinY() == alignedMinY) ? metadata[5] : 0;
+ alignedMaxX += (rasterExtent.getMaxX() == alignedMaxX) ? metadata[4] : 0;
+ alignedMaxY -= (rasterExtent.getMaxY() == alignedMaxY) ? metadata[5] : 0;
+ }
+
+ alignedMinX = Math.max(alignedMinX, originalMinX);
+ alignedMinY = Math.max(alignedMinY, originalMinY);
+ alignedMaxX = Math.min(alignedMaxX, originalMaxX);
+ alignedMaxY = Math.min(alignedMaxY, originalMaxY);
+
+ // Create the aligned raster extent
+ Envelope2D alignedRasterExtent =
+ new Envelope2D(
+ rasterExtent.getCoordinateReferenceSystem(),
+ alignedMinX,
+ alignedMinY,
+ alignedMaxX - alignedMinX,
+ alignedMaxY - alignedMinY);
+
+ return alignedRasterExtent;
+ }
+
+ private static RasterizationParams calculateRasterizationParams(
+ GridCoverage2D raster,
+ boolean useGeometryExtent,
+ double[] metadata,
+ Envelope2D geomExtent,
+ String pixelType) {
+
+ double upperLeftX = 0;
+ double upperLeftY = 0;
+ if (useGeometryExtent) {
+ upperLeftX = geomExtent.getMinX();
+ upperLeftY = geomExtent.getMaxY();
+ } else {
+ upperLeftX = metadata[0];
+ upperLeftY = metadata[1];
+ }
+
+ WritableRaster writableRaster;
+ if (useGeometryExtent) {
+ int geomExtentWidth = (int) (Math.round(geomExtent.getWidth() /
metadata[4]));
+ int geomExtentHeight = (int) (Math.round(geomExtent.getHeight() /
-metadata[5]));
+
+ writableRaster =
+ RasterFactory.createBandedRaster(
+ RasterUtils.getDataTypeCode(pixelType), geomExtentWidth,
geomExtentHeight, 1, null);
+ } else {
+ int rasterWidth = (int)
raster.getGridGeometry().getGridRange2D().getWidth();
+ int rasterHeight = (int)
raster.getGridGeometry().getGridRange2D().getHeight();
+ writableRaster =
+ RasterFactory.createBandedRaster(
+ RasterUtils.getDataTypeCode(pixelType), rasterWidth,
rasterHeight, 1, null);
+ }
+
+ return new RasterizationParams(
+ writableRaster, pixelType, metadata[4], -metadata[5], upperLeftX,
upperLeftY);
+ }
+
+ private static void validateRasterMetadata(double[] rasterMetadata) {
+ if (rasterMetadata[4] < 0) {
+ throw new IllegalArgumentException("ScaleX cannot be negative");
+ }
+ if (rasterMetadata[5] > 0) {
+ throw new IllegalArgumentException("ScaleY must be negative");
+ }
+ if (rasterMetadata[6] != 0 || rasterMetadata[7] != 0) {
+ throw new IllegalArgumentException("SkewX and SkewY must be zero");
+ }
+ }
+
+ private static class RasterizationParams {
+ WritableRaster writableRaster;
+ String pixelType;
+ double scaleX;
+ double scaleY;
+ double upperLeftX;
+ double upperLeftY;
+
+ RasterizationParams(
+ WritableRaster writableRaster,
+ String pixelType,
+ double scaleX,
+ double scaleY,
+ double upperLeftX,
+ double upperLeftY) {
+ this.writableRaster = writableRaster;
+ this.pixelType = pixelType;
+ this.scaleX = scaleX;
+ this.scaleY = scaleY;
+ this.upperLeftX = upperLeftX;
+ this.upperLeftY = upperLeftY;
+ }
+ }
+
+ public static void rasterizePolygon(
+ Geometry geom,
+ RasterizationParams params,
+ Envelope2D geomExtent,
+ double value,
+ boolean allTouched) {
+ if (!(geom instanceof Polygon)) {
+ throw new IllegalArgumentException("Only Polygon geometry is supported");
+ }
+
+ // Clip polygon to rasterExtent
+ // Using buffer(0) to avoid TopologyException : found non-noded
intersection.
+ Geometry clippedGeom =
+ Functions.intersection(JTS.toGeometry((BoundingBox) geomExtent),
Functions.buffer(geom, 0));
+
+ if (Objects.equals(clippedGeom.getGeometryType(), "MultiPolygon")) {
+ for (int i = 0; i < clippedGeom.getNumGeometries(); i++) {
+ Geometry subGeom = clippedGeom.getGeometryN(i);
+ rasterizePolygon(subGeom, params, geomExtent, value, allTouched);
+ }
+ return;
+ }
+
+ Polygon polygon = (Polygon) clippedGeom;
+
+ // Compute scanline X-intercepts
+ Map<Double, TreeSet<Double>> scanlineIntersections =
+ computeScanlineIntersections(polygon, params, value, geomExtent,
allTouched);
+
+ // Process intersections to get startXs and endXs for each scanline
+ Map<Integer, List<int[]>> scanlineFillRanges =
computeFillRanges(scanlineIntersections);
+
+ // Burn values between startX and endX pairs
+ fillPolygon(scanlineFillRanges, params, value);
+ }
+
+ /** Computes scanline intersections by iterating over polygon edges. */
+ private static Map<Double, TreeSet<Double>> computeScanlineIntersections(
+ Polygon polygon,
+ RasterizationParams params,
+ double value,
+ Envelope2D geomExtent,
+ boolean allTouched) {
+
+ Map<Double, TreeSet<Double>> scanlineIntersections = new HashMap<>();
+ List<LineString> allRings = new ArrayList<>();
+ allRings.add(polygon.getExteriorRing());
+
+ for (int i = 0; i < polygon.getNumInteriorRing(); i++) {
+ allRings.add(polygon.getInteriorRingN(i));
+ }
+
+ for (LineString ring : allRings) {
+ Coordinate[] coords = ring.getCoordinates();
+ int numPoints = coords.length;
+
+ if (allTouched) {
+ rasterizeLineString(ring, params, value, geomExtent);
+ }
+
+ for (int i = 0; i < numPoints - 1; i++) {
+ Coordinate worldP1 = coords[i];
+ Coordinate worldP2 = coords[i + 1];
+
+ // Ensure worldP1.y ≤ worldP2.y
+ if (worldP1.y > worldP2.y) {
+ Coordinate temp = worldP1;
+ worldP1 = worldP2;
+ worldP2 = temp;
+ }
+
+ if (worldP1.y == worldP2.y) {
+ continue; // Skip horizontal edges
+ }
+
+ // Calculate scan line limits to iterate between for each segment
+ // Using BigDecimal to avoid floating point errors
+ double yStart =
+ Math.ceil(
+ (BigDecimal.valueOf(params.upperLeftY)
+ .subtract(BigDecimal.valueOf(worldP1.y))
+ .divide(BigDecimal.valueOf(params.scaleY),
RoundingMode.CEILING))
+ .doubleValue());
+
+ double yEnd =
+ Math.floor(
+ (BigDecimal.valueOf(params.upperLeftY)
+ .subtract(BigDecimal.valueOf(worldP2.y))
+ .divide(BigDecimal.valueOf(params.scaleY),
RoundingMode.FLOOR))
+ .doubleValue());
+
+ // Contain y range within geomExtent; Use centroid y line as scan line
+ yEnd = Math.max(0.5, Math.abs(yEnd) + 0.5);
+ yStart = Math.min((params.writableRaster.getHeight() - 0.5),
Math.abs(yStart) - 0.5);
+
+ double p1X = (worldP1.x - params.upperLeftX) / params.scaleX;
+ double p1Y = (params.upperLeftY - worldP1.y) / params.scaleY;
+
+ if (worldP1.x == worldP2.x) {
+ // Vertical line case: directly set xIntercept. Avoid divide by zero
error when
+ // calculating slope
+ for (double y = yStart; y >= yEnd; y--) {
+ double xIntercept = p1X; // Vertical line, xIntercept is constant
+ if (xIntercept < 0 || xIntercept >
params.writableRaster.getWidth()) {
+ continue; // Skip xIntercepts outside geomExtent
+ }
+ scanlineIntersections.computeIfAbsent(y, k -> new
TreeSet<>()).add(xIntercept);
+ }
+ } else {
+ double slope = (worldP2.y - worldP1.y) / (worldP2.x - worldP1.x);
+ // System.out.println("slope: " + slope);
+
+ for (double y = yStart; y >= yEnd; y--) {
+ double xIntercept = p1X + ((p1Y - y) / slope);
+ if ((xIntercept < 0) || (xIntercept >
params.writableRaster.getWidth())) {
+ continue; // skip xIntercepts outside geomExtent
+ }
+ if (!scanlineIntersections.containsKey(y)) {
+ scanlineIntersections.put(y, new TreeSet<>());
+ }
+ scanlineIntersections.get(y).add(xIntercept);
+ }
+ }
+ }
+ }
+ return scanlineIntersections;
+ }
+
+ /** Computes startX and endX pairs for each scanline by leveraging sorted
X-intercepts. */
+ private static Map<Integer, List<int[]>> computeFillRanges(
+ Map<Double, TreeSet<Double>> scanlineIntersections) {
+
+ Map<Integer, List<int[]>> scanlineFillRanges = new HashMap<>();
+
+ for (Map.Entry<Double, TreeSet<Double>> entry :
scanlineIntersections.entrySet()) {
+ double y = entry.getKey();
+ TreeSet<Double> xIntercepts = entry.getValue();
+ List<int[]> ranges = new ArrayList<>();
+
+ Iterator<Double> it = xIntercepts.iterator();
+ while (it.hasNext()) {
+ double x1 = it.next();
+ if (!it.hasNext()) {
+ ranges.add(new int[] {(int) Math.floor(x1)});
+ continue;
+ }
+ double x2 = it.next();
+
+ // Compute centroids
+ double firstCentroid = Math.ceil(x1 - 0.5) + 0.5;
+ double lastCentroid = Math.floor(x2 + 0.5) - 0.5;
+
+ if (lastCentroid < firstCentroid) {
+ continue; // No valid pixel to include
+ }
+ if (lastCentroid == firstCentroid) {
+ ranges.add(new int[] {(int) Math.floor(firstCentroid)});
+ } else {
+ ranges.add(new int[] {(int) Math.floor(firstCentroid), (int)
Math.floor(lastCentroid)});
+ }
+ }
+
+ scanlineFillRanges.put((int) Math.floor(y), ranges);
+ }
+
+ return scanlineFillRanges;
+ }
+
+ /** Burns pixels between startX and endX for each scanline. */
+ private static void fillPolygon(
+ Map<Integer, List<int[]>> scanlineFillRanges, RasterizationParams
params, double value) {
+
+ for (Map.Entry<Integer, List<int[]>> entry :
scanlineFillRanges.entrySet()) {
+ int y = entry.getKey();
+ for (int[] range : entry.getValue()) {
+ if (range.length == 1) {
+ params.writableRaster.setSample(range[0], y, 0, value);
+ continue;
+ }
+ int xStart = range[0];
+ int xEnd = range[1];
+
+ for (int x = xStart; x <= xEnd; x++) {
+ params.writableRaster.setSample(x, y, 0, value);
+ }
+ }
+ }
+ }
+}
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 d3c1ed5c7e..ef34a5060d 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
@@ -32,7 +32,8 @@ import org.opengis.referencing.operation.TransformException;
public class FunctionEditorsTest extends RasterTestBase {
@Test
- public void testSetValuesWithEmptyRaster() throws FactoryException {
+ public void testSetValuesWithEmptyRaster()
+ throws FactoryException, ParseException, TransformException {
GridCoverage2D emptyRaster = RasterConstructors.makeEmptyRaster(1, 5, 5,
0, 0, 1, -1, 0, 0, 0);
double[] values =
new double[] {1, 1, 1, 0, 0, 0, 1, 2, 3, 3, 5, 6, 7, 0, 0, 3, 0, 0, 3,
0, 0, 0, 0, 0, 0};
@@ -56,6 +57,41 @@ public class FunctionEditorsTest extends RasterTestBase {
17.0, 18.0, 19.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
};
assertArrayEquals(actual, expected, 0.0);
+
+ // Test Rasterization with GeometryCollection
+ emptyRaster = RasterConstructors.makeEmptyRaster(1, 10, 10, 0, 0, 1, -1,
0, 0, 0);
+ String polygon =
+ "GEOMETRYCOLLECTION(POINT(1.5 -4.5), POINT(2 -2), MULTIPOINT((4 -4),
(6 -6)), LINESTRING(1 -1, 8 -1), MULTILINESTRING((2 -8, 7 -8), (3 -6, 5 -6)),
POLYGON((3 -3, 5 -3, 5 -5, 3 -5, 3 -3)), MULTIPOLYGON(((6 -2, 8 -2, 8 -4, 6 -4,
6 -2)),((1 -7, 2 -7, 2 -9, 1 -9, 1 -7))))";
+ Geometry geom = Constructors.geomFromWKT(polygon, 0);
+ GridCoverage2D result = PixelFunctionEditors.setValues(emptyRaster, 1,
geom, 10, false, false);
+ actual = MapAlgebra.bandAsArray(result, 1);
+ expected =
+ new double[] {
+ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 10.0, 10.0,
10.0, 10.0, 10.0, 10.0,
+ 10.0, 0.0, 0.0, 0.0, 10.0, 10.0, 0.0, 0.0, 0.0, 10.0, 10.0, 0.0,
0.0, 0.0, 0.0, 0.0, 10.0,
+ 10.0, 0.0, 10.0, 10.0, 0.0, 0.0, 0.0, 10.0, 0.0, 10.0, 10.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0,
+ 0.0, 0.0, 0.0, 0.0, 10.0, 10.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 10.0,
10.0, 10.0, 10.0, 0.0,
+ 0.0, 0.0, 0.0, 10.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
10.0, 10.0, 10.0, 10.0,
+ 10.0, 10.0, 10.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0
+ };
+ assertArrayEquals(actual, expected, 0.0);
+
+ // Test Rasterization on Polygons with holes
+ polygon =
+ "POLYGON ((2.5 -2.5, 7.5 -2.5, 7.5 -7.5, 2.5 -7.5, 2.5 -2.5), (4 -4, 6
-4, 6 -6, 4 -6, 4 -4))";
+ geom = Constructors.geomFromWKT(polygon, 0);
+ result = PixelFunctionEditors.setValues(emptyRaster, 1, geom, 10, false,
false);
+ actual = MapAlgebra.bandAsArray(result, 1);
+ expected =
+ new double[] {
+ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0,
+ 0.0, 0.0, 0.0, 0.0, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0, 0.0, 0.0,
0.0, 0.0, 10.0, 10.0,
+ 10.0, 10.0, 10.0, 10.0, 0.0, 0.0, 0.0, 0.0, 10.0, 10.0, 0.0, 0.0,
10.0, 10.0, 0.0, 0.0,
+ 0.0, 0.0, 10.0, 10.0, 0.0, 0.0, 10.0, 10.0, 0.0, 0.0, 0.0, 0.0,
10.0, 10.0, 10.0, 10.0,
+ 10.0, 10.0, 0.0, 0.0, 0.0, 0.0, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0,
0.0, 0.0, 0.0, 0.0,
+ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0
+ };
+ assertArrayEquals(actual, expected, 0.0);
}
@Test
@@ -67,13 +103,20 @@ public class FunctionEditorsTest extends RasterTestBase {
"POLYGON ((236722 4204770, 243900 4204770, 243900 4197590, 236722
4197590, 236722 4204770))";
Geometry geom = Constructors.geomFromWKT(polygon, 26918);
- GridCoverage2D result = PixelFunctionEditors.setValues(raster, 1, geom,
10, false);
+ GridCoverage2D result = PixelFunctionEditors.setValues(raster, 1, geom,
10, true, false);
- Geometry point = Constructors.geomFromWKT("POINT (243700 4197797)", 26918);
+ Geometry point = Constructors.geomFromWKT("POINT (236698.272 4204794.5)",
26918);
double actual = PixelFunctions.value(result, point, 1);
double expected = 10.0;
assertEquals(expected, actual, 0d);
+ result = PixelFunctionEditors.setValues(raster, 1, geom, 10, false);
+
+ point = Constructors.geomFromWKT("POINT (243700 4197797)", 26918);
+ actual = PixelFunctions.value(result, point, 1);
+ expected = 10.0;
+ assertEquals(expected, actual, 0d);
+
point = Constructors.geomFromWKT("POINT (240311 4202806)", 26918);
actual = PixelFunctions.value(result, point, 1);
assertEquals(expected, actual, 0d);
@@ -92,13 +135,20 @@ public class FunctionEditorsTest extends RasterTestBase {
"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);
- GridCoverage2D result = PixelFunctionEditors.setValues(raster, 1, geom,
10, false);
+ GridCoverage2D result = PixelFunctionEditors.setValues(raster, 1, geom,
10, true, false);
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);
+ result = PixelFunctionEditors.setValues(raster, 1, geom, 10, false);
+
+ point = Constructors.geomFromWKT("POINT (-77.9146 37.8916)", 0);
+ actual = PixelFunctions.value(result, point, 1);
+ expected = 10.0;
+ assertEquals(expected, actual, 0d);
+
point = Constructors.geomFromWKT("POINT (-77.9549 37.9357)", 0);
actual = PixelFunctions.value(result, point, 1);
assertEquals(expected, actual, 0d);
@@ -113,12 +163,18 @@ public class FunctionEditorsTest extends RasterTestBase {
"POLYGON ((-8682522.873537656 4572703.890837922, -8673439.664183248
4572993.532747675, -8673155.57366801 4563873.2099182755, -8701890.325907696
4562931.7093397, -8682522.873537656 4572703.890837922))";
Geometry geom = Constructors.geomFromWKT(polygon, 3857);
- GridCoverage2D result = PixelFunctionEditors.setValues(raster, 1, geom,
10, false);
+ GridCoverage2D result = PixelFunctionEditors.setValues(raster, 1, geom,
10, true, false);
Geometry point = Constructors.geomFromWKT("POINT (243700 4197797)", 26918);
double actual = PixelFunctions.value(result, point, 1);
double expected = 10.0;
assertEquals(expected, actual, 0d);
+ result = PixelFunctionEditors.setValues(raster, 1, geom, 10, false);
+
+ point = Constructors.geomFromWKT("POINT (243700 4197797)", 26918);
+ actual = PixelFunctions.value(result, point, 1);
+ expected = 10.0;
+ assertEquals(expected, actual, 0d);
point = Constructors.geomFromWKT("POINT (235749.0869 4200557.7397)",
26918);
actual = PixelFunctions.value(result, point, 1);
@@ -150,7 +206,7 @@ public class FunctionEditorsTest extends RasterTestBase {
double[] actual = MapAlgebra.bandAsArray(raster, 1);
double[] expected =
new double[] {
- 4235.0, 0.0, 0.0, 0.0, 4235.0, 4235.0, 0.0, 0.0, 4235.0, 4235.0,
4235.0, 0.0, 4235.0, 0.0,
+ 4235.0, 0.0, 0.0, 0.0, 4235.0, 4235.0, 0.0, 0.0, 4235.0, 0.0,
4235.0, 0.0, 4235.0, 0.0,
0.0, 4235.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
};
assertArrayEquals(expected, actual, 0.1d);
@@ -161,7 +217,7 @@ public class FunctionEditorsTest extends RasterTestBase {
actual = MapAlgebra.bandAsArray(raster, 1);
expected =
new double[] {
- 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0,
+ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 35.0, 0.0, 0.0,
0.0, 35.0, 0.0, 0.0, 0.0, 0.0
};
assertArrayEquals(expected, actual, 0.1d);
@@ -172,8 +228,8 @@ public class FunctionEditorsTest extends RasterTestBase {
actual = MapAlgebra.bandAsArray(raster, 1);
expected =
new double[] {
- 0.0, 400.0, 0.0, 0.0, 0.0, 400.0, 0.0, 0.0, 0.0, 0.0, 400.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0,
- 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 400.0
+ 400.0, 400.0, 0.0, 0.0, 400.0, 400.0, 400.0, 0.0, 0.0, 400.0, 400.0,
0.0, 0.0, 0.0, 0.0,
+ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 400.0, 400.0
};
assertArrayEquals(expected, actual, 0.1d);
diff --git
a/common/src/test/java/org/apache/sedona/common/raster/RasterBandAccessorsTest.java
b/common/src/test/java/org/apache/sedona/common/raster/RasterBandAccessorsTest.java
index e4fe01d575..94986d039a 100644
---
a/common/src/test/java/org/apache/sedona/common/raster/RasterBandAccessorsTest.java
+++
b/common/src/test/java/org/apache/sedona/common/raster/RasterBandAccessorsTest.java
@@ -92,39 +92,40 @@ public class RasterBandAccessorsTest extends RasterTestBase
{
"POLYGON ((236722 4204770, 243900 4204770, 243900 4197590, 221170
4197590, 236722 4204770))";
Geometry geom = Constructors.geomFromWKT(polygon,
RasterAccessors.srid(raster));
- double actual = RasterBandAccessors.getZonalStats(raster, geom, 1, "sum",
false);
- double expected = 1.0719726E7;
+ double actual = RasterBandAccessors.getZonalStats(raster, geom, 1, "sum",
false, false);
+ double expected = 1.0896994E7;
assertEquals(expected, actual, 0d);
- actual = RasterBandAccessors.getZonalStats(raster, geom, 2, "mean", false);
- expected = 220.7527;
+ actual = RasterBandAccessors.getZonalStats(raster, geom, 2, "mean", false,
false);
+ expected = 220.76055239764008;
assertEquals(expected, actual, FP_TOLERANCE);
actual = RasterBandAccessors.getZonalStats(raster, geom, 1, "count");
- expected = 184792.0;
+ expected = 185953.0;
assertEquals(expected, actual, 0.1d);
- actual = RasterBandAccessors.getZonalStats(raster, geom, 3, "variance",
false);
- expected = 13549.6263;
+ actual = RasterBandAccessors.getZonalStats(raster, geom, 3, "variance",
false, false);
+ expected = 13560.056183580346;
assertEquals(expected, actual, FP_TOLERANCE);
actual = RasterBandAccessors.getZonalStats(raster, geom, "max");
expected = 255.0;
assertEquals(expected, actual, 1E-1);
- actual = RasterBandAccessors.getZonalStats(raster, geom, 1, "min", false);
+ actual = RasterBandAccessors.getZonalStats(raster, geom, 1, "min", false,
false);
expected = 0.0;
assertEquals(expected, actual, 1E-1);
- actual = RasterBandAccessors.getZonalStats(raster, geom, 1, "sd", false);
- expected = 92.1500;
+ actual = RasterBandAccessors.getZonalStats(raster, geom, 1, "sd", false,
false);
+ expected = 92.53811912983977;
assertEquals(expected, actual, FP_TOLERANCE);
geom =
Constructors.geomFromWKT(
"POLYGON ((-77.96672569800863073 37.91971182746296876,
-77.9688630154902711 37.89620133516485367, -77.93936803424354309
37.90517806858776595, -77.96672569800863073 37.91971182746296876))",
0);
- Double statValue = RasterBandAccessors.getZonalStats(raster, geom, 1,
"sum", false, true);
+ Double statValue =
+ RasterBandAccessors.getZonalStats(raster, geom, 1, "sum", false,
false, true);
assertNotNull(statValue);
Geometry nonIntersectingGeom =
@@ -132,12 +133,14 @@ public class RasterBandAccessorsTest extends
RasterTestBase {
"POLYGON ((-78.22106647832458748 37.76411511479908967,
-78.20183062098976734 37.72863564460374874, -78.18088490966962922
37.76753482276972562, -78.22106647832458748 37.76411511479908967))",
0);
statValue =
- RasterBandAccessors.getZonalStats(raster, nonIntersectingGeom, 1,
"sum", false, true);
+ RasterBandAccessors.getZonalStats(
+ raster, nonIntersectingGeom, 1, "sum", false, false, true);
assertNull(statValue);
assertThrows(
IllegalArgumentException.class,
() ->
- RasterBandAccessors.getZonalStats(raster, nonIntersectingGeom, 1,
"sum", false, false));
+ RasterBandAccessors.getZonalStats(
+ raster, nonIntersectingGeom, 1, "sum", false, false, false));
}
@Test
@@ -149,32 +152,32 @@ public class RasterBandAccessorsTest extends
RasterTestBase {
// Testing implicit CRS transformation
Geometry geom = Constructors.geomFromWKT(polygon, 0);
- double actual = RasterBandAccessors.getZonalStats(raster, geom, 1, "sum",
true);
- double expected = 3213526.0;
+ double actual = RasterBandAccessors.getZonalStats(raster, geom, 1, "sum",
false, true);
+ double expected = 3229013.0;
assertEquals(expected, actual, 0d);
- actual = RasterBandAccessors.getZonalStats(raster, geom, 1, "mean", true);
- expected = 226.5599;
+ actual = RasterBandAccessors.getZonalStats(raster, geom, 1, "mean", false,
true);
+ expected = 226.61330619692416;
assertEquals(expected, actual, FP_TOLERANCE);
actual = RasterBandAccessors.getZonalStats(raster, geom, 1, "count");
- expected = 14184.0;
+ expected = 14249.0;
assertEquals(expected, actual, 0.1d);
actual = RasterBandAccessors.getZonalStats(raster, geom, "variance");
- expected = 5606.4233;
+ expected = 5596.966403503485;
assertEquals(expected, actual, FP_TOLERANCE);
actual = RasterBandAccessors.getZonalStats(raster, geom, "max");
expected = 255.0;
assertEquals(expected, actual, 1E-1);
- actual = RasterBandAccessors.getZonalStats(raster, geom, 1, "min", true);
+ actual = RasterBandAccessors.getZonalStats(raster, geom, 1, "min", false,
true);
expected = 1.0;
assertEquals(expected, actual, 1E-1);
- actual = RasterBandAccessors.getZonalStats(raster, geom, 1, "sd", true);
- expected = 74.8760;
+ actual = RasterBandAccessors.getZonalStats(raster, geom, 1, "sd", false,
true);
+ expected = 74.81287592054916;
assertEquals(expected, actual, FP_TOLERANCE);
}
@@ -187,16 +190,16 @@ public class RasterBandAccessorsTest extends
RasterTestBase {
"POLYGON ((-8673439.6642 4572993.5327, -8673155.5737 4563873.2099,
-8701890.3259 4562931.7093, -8682522.8735 4572703.8908, -8673439.6642
4572993.5327))";
Geometry geom = Constructors.geomFromWKT(polygon, 3857);
- double[] actual = RasterBandAccessors.getZonalStatsAll(raster, geom, 1,
false);
+ double[] actual = RasterBandAccessors.getZonalStatsAll(raster, geom, 1,
false, false, false);
double[] expected =
new double[] {
- 184792.0,
- 1.0719726E7,
- 58.00968656653401,
+ 185953.0,
+ 1.0896994E7,
+ 58.600796975566816,
0.0,
0.0,
- 92.15004748703687,
- 8491.631251863151,
+ 92.53811912983977,
+ 8563.303492088418,
0.0,
255.0
};
@@ -206,31 +209,46 @@ public class RasterBandAccessorsTest extends
RasterTestBase {
Constructors.geomFromWKT(
"POLYGON ((-77.96672569800863073 37.91971182746296876,
-77.9688630154902711 37.89620133516485367, -77.93936803424354309
37.90517806858776595, -77.96672569800863073 37.91971182746296876))",
0);
- actual = RasterBandAccessors.getZonalStatsAll(raster, geom, 1, false);
+ actual = RasterBandAccessors.getZonalStatsAll(raster, geom, 1, false,
false, false);
assertNotNull(actual);
Geometry nonIntersectingGeom =
Constructors.geomFromWKT(
"POLYGON ((-78.22106647832458748 37.76411511479908967,
-78.20183062098976734 37.72863564460374874, -78.18088490966962922
37.76753482276972562, -78.22106647832458748 37.76411511479908967))",
0);
- actual = RasterBandAccessors.getZonalStatsAll(raster, nonIntersectingGeom,
1, false, true);
+ actual =
+ RasterBandAccessors.getZonalStatsAll(raster, nonIntersectingGeom, 1,
false, false, true);
assertNull(actual);
assertThrows(
IllegalArgumentException.class,
- () -> RasterBandAccessors.getZonalStatsAll(raster,
nonIntersectingGeom, 1, false, false));
+ () ->
+ RasterBandAccessors.getZonalStatsAll(
+ raster, nonIntersectingGeom, 1, false, false, false));
}
@Test
- public void testZonalStatsAllWithNoData() throws IOException,
FactoryException, ParseException {
+ public void testZonalStatsAllWithNoData()
+ throws IOException, FactoryException, ParseException, TransformException
{
GridCoverage2D raster =
rasterFromGeoTiff(resourceFolder +
"raster/raster_with_no_data/test5.tiff");
String polygon =
"POLYGON((-167.750000 87.750000, -155.250000 87.750000, -155.250000
40.250000, -180.250000 40.250000, -167.750000 87.750000))";
Geometry geom = Constructors.geomFromWKT(polygon,
RasterAccessors.srid(raster));
- double[] actual = RasterBandAccessors.getZonalStatsAll(raster, geom, 1,
true);
+ double[] actual = RasterBandAccessors.getZonalStatsAll(raster, geom, 1,
false, true);
double[] expected =
- new double[] {14184.0, 3213526.0, 226.5599, 255.0, 255.0, 74.8760,
5606.4233, 1.0, 255.0};
+ new double[] {
+ 14249.0,
+ 3229013.0,
+ 226.61330619692416,
+ 255.0,
+ 255.0,
+ 74.81287592054916,
+ 5596.966403503485,
+ 1.0,
+ 255.0
+ };
+
assertArrayEquals(expected, actual, FP_TOLERANCE);
}
@@ -247,7 +265,7 @@ public class RasterBandAccessorsTest extends RasterTestBase
{
// Testing implicit CRS transformation
Geometry geom = Constructors.geomFromWKT("POLYGON((2 -2, 2 -6, 6 -6, 6 -2,
2 -2))", 0);
- double[] actual = RasterBandAccessors.getZonalStatsAll(raster, geom, 1,
true);
+ double[] actual = RasterBandAccessors.getZonalStatsAll(raster, geom, 1,
false, true);
double[] expected = new double[] {13.0, 114.0, 8.7692, 9.0, 11.0, 4.7285,
22.3589, 1.0, 16.0};
assertArrayEquals(expected, actual, FP_TOLERANCE);
}
diff --git
a/common/src/test/java/org/apache/sedona/common/raster/RasterBandEditorsTest.java
b/common/src/test/java/org/apache/sedona/common/raster/RasterBandEditorsTest.java
index 6a560b7280..73f81e8ed1 100644
---
a/common/src/test/java/org/apache/sedona/common/raster/RasterBandEditorsTest.java
+++
b/common/src/test/java/org/apache/sedona/common/raster/RasterBandEditorsTest.java
@@ -184,7 +184,7 @@ public class RasterBandEditorsTest extends RasterTestBase {
"POLYGON ((-8682522.873537656 4572703.890837922, -8673439.664183248
4572993.532747675, -8673155.57366801 4563873.2099182755, -8701890.325907696
4562931.7093397, -8682522.873537656 4572703.890837922))";
Geometry geom = Constructors.geomFromWKT(polygon, 3857);
- GridCoverage2D clippedRaster = RasterBandEditors.clip(raster, 1, geom,
200, false);
+ GridCoverage2D clippedRaster = RasterBandEditors.clip(raster, 1, geom,
false, 200, false);
double[] clippedMetadata =
Arrays.stream(RasterAccessors.metadata(clippedRaster), 0, 9).toArray();
double[] originalMetadata =
Arrays.stream(RasterAccessors.metadata(raster), 0, 9).toArray();
@@ -211,7 +211,7 @@ public class RasterBandEditorsTest extends RasterTestBase {
"POLYGON ((236722 4204770, 243900 4204770, 243900 4197590, 221170
4197590, 236722 4204770))";
Geometry geom = Constructors.geomFromWKT(polygon,
RasterAccessors.srid(raster));
- GridCoverage2D clippedRaster = RasterBandEditors.clip(raster, 1, geom,
200, false);
+ GridCoverage2D clippedRaster = RasterBandEditors.clip(raster, 1, geom,
false, 200, false);
double[] clippedMetadata =
Arrays.stream(RasterAccessors.metadata(clippedRaster), 0, 9).toArray();
double[] originalMetadata =
Arrays.stream(RasterAccessors.metadata(raster), 0, 9).toArray();
@@ -232,7 +232,7 @@ public class RasterBandEditorsTest extends RasterTestBase {
Double[] expectedValues = new Double[] {null, null, 0.0, 0.0, null};
assertTrue(Arrays.equals(expectedValues, actualValues));
- GridCoverage2D croppedRaster = RasterBandEditors.clip(raster, 1, geom,
200, true);
+ GridCoverage2D croppedRaster = RasterBandEditors.clip(raster, 1, geom,
false, 200, true);
assertEquals(0, croppedRaster.getRenderedImage().getMinX());
assertEquals(0, croppedRaster.getRenderedImage().getMinY());
GridCoverage2D croppedRaster2 =
Serde.deserialize(Serde.serialize(croppedRaster));
@@ -263,10 +263,11 @@ public class RasterBandEditorsTest extends RasterTestBase
{
// Throws an exception in non-lenient mode
assertThrows(
CannotCropException.class,
- () -> RasterBandEditors.clip(raster, 1, nonIntersectingGeom, 200,
false, false));
+ () -> RasterBandEditors.clip(raster, 1, nonIntersectingGeom, false,
200, false, false));
// Returns null in lenient mode
- GridCoverage2D result = RasterBandEditors.clip(raster, 1,
nonIntersectingGeom, 200, false);
+ GridCoverage2D result =
+ RasterBandEditors.clip(raster, 1, nonIntersectingGeom, false, 200,
false);
assertNull(result);
raster.dispose(true);
}
diff --git
a/common/src/test/java/org/apache/sedona/common/raster/RasterConstructorsTest.java
b/common/src/test/java/org/apache/sedona/common/raster/RasterConstructorsTest.java
index a79165b4ee..e5b7c2d2bf 100644
---
a/common/src/test/java/org/apache/sedona/common/raster/RasterConstructorsTest.java
+++
b/common/src/test/java/org/apache/sedona/common/raster/RasterConstructorsTest.java
@@ -81,131 +81,157 @@ public class RasterConstructorsTest extends
RasterTestBase {
// Polygon
GridCoverage2D raster =
RasterConstructors.makeEmptyRaster(2, 255, 255, 1, -1, 2, -2, 0, 0,
4326);
- Geometry geom = Constructors.geomFromWKT("POLYGON((15 15, 18 20, 15 24, 24
25, 15 15))", 0);
- GridCoverage2D rasterized = RasterConstructors.asRaster(geom, raster, "d",
3093151, 3d);
+
+ Geometry geom =
+ Constructors.geomFromWKT("POLYGON((15 -15, 18 -20, 15 -24, 24 -25, 15
-15))", 0);
+ GridCoverage2D rasterized = RasterConstructors.asRaster(geom, raster, "d",
true, 3093151, 3d);
+
double[] actual = MapAlgebra.bandAsArray(rasterized, 1);
+
double[] expected =
new double[] {
- 3093151.0, 3093151.0, 3093151.0, 3093151.0, 0.0, 0.0, 3093151.0,
3093151.0, 0.0, 0.0, 0.0,
- 3093151.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0
+ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3093151.0, 3093151.0, 0.0, 0.0, 0.0,
0.0, 0.0, 3093151.0,
+ 3093151.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3093151.0, 3093151.0, 0.0, 0.0,
0.0, 3093151.0,
+ 3093151.0, 3093151.0, 3093151.0, 0.0, 3093151.0, 3093151.0,
3093151.0, 3093151.0,
+ 3093151.0, 3093151.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3093151.0
};
+
assertArrayEquals(expected, actual, 0.1d);
rasterized = RasterConstructors.asRaster(geom, raster, "d");
actual = MapAlgebra.bandAsArray(rasterized, 1);
expected =
new double[] {
- 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0,
- 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
+ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0,
0.0, 0.0, 0.0, 1.0, 1.0,
+ 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 0.0
};
+
assertArrayEquals(expected, actual, 0.1d);
// MultiPolygon
geom =
Constructors.geomFromWKT(
- "MULTIPOLYGON (((15 15, 1.5 5.5, 3.5 5.5, 3.5 1.5, 15 15)), ((4.4
2.4, 4.4 6.4, 6.4 6.4, 6.4 2.4, 4.4 2.4)))",
+ "MULTIPOLYGON ( ((2 -2, 4 -2, 4 -4, 2 -4, 2 -2)), ((4 -4, 6 -4, 6
-6, 5 -7, 4 -6, 4 -4)), ((6 -6, 8 -6, 8 -8, 6 -8, 6 -6)), ((8 -6, 10 -6, 10 -4,
9 -3, 8 -4, 8 -6)) )",
0);
- rasterized = RasterConstructors.asRaster(geom, raster, "d", 3093151, 3d);
+
+ rasterized = RasterConstructors.asRaster(geom, raster, "d", true, 3093151,
3d, true);
actual = MapAlgebra.bandAsArray(rasterized, 1);
expected =
new double[] {
- 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0,
- 0.0, 0.0, 0.0, 0.0, 0.0, 3093151.0, 3093151.0, 0.0, 0.0, 0.0, 0.0,
3093151.0, 3093151.0,
- 0.0, 0.0, 0.0, 0.0, 0.0, 3093151.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0,
- 0.0
+ 3093151.0, 3093151.0, 0.0, 0.0, 0.0, 3093151.0, 3093151.0,
3093151.0, 3093151.0,
+ 3093151.0, 0.0, 3093151.0, 3093151.0, 3093151.0, 3093151.0, 0.0,
0.0, 3093151.0,
+ 3093151.0, 0.0
};
assertArrayEquals(expected, actual, 0.1d);
rasterized = RasterConstructors.asRaster(geom, raster, "d");
actual = MapAlgebra.bandAsArray(rasterized, 1);
+
expected =
new double[] {
- 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0,
- 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0,
0.0, 0.0, 0.0, 0.0, 0.0,
- 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
+ 1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 1.0, 1.0,
1.0, 1.0, 0.0, 0.0, 1.0,
+ 1.0, 0.0
};
assertArrayEquals(expected, actual, 0.1d);
// MultiLineString
- geom = Constructors.geomFromWKT("MULTILINESTRING ((5 5, 10 10), (10 10, 15
15, 20 20))", 0);
- rasterized = RasterConstructors.asRaster(geom, raster, "d", 3093151, 3d);
+ geom =
+ Constructors.geomFromWKT("MULTILINESTRING ((5 -5, 10 -10), (10 -10, 15
-15, 20 -20))", 0);
+ rasterized = RasterConstructors.asRaster(geom, raster, "d", false,
3093151, 3d);
+
actual = MapAlgebra.bandAsArray(rasterized, 1);
expected =
new double[] {
- 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3093151.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 3093151.0,
- 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3093151.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 3093151.0,
- 0.0, 0.0, 0.0, 0.0, 0.0, 3093151.0, 3093151.0, 0.0, 0.0, 0.0, 0.0,
0.0, 3093151.0, 0.0,
- 0.0, 0.0, 0.0, 0.0, 0.0, 3093151.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
3093151.0, 0.0, 0.0,
- 0.0, 0.0, 0.0, 0.0, 0.0
+ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3093151.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0,
+ 0.0, 0.0, 0.0, 3093151.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 3093151.0, 0.0,
+ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3093151.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0,
+ 0.0, 3093151.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
3093151.0, 0.0, 0.0, 0.0,
+ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3093151.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0,
+ 3093151.0
};
assertArrayEquals(expected, actual, 0.1d);
rasterized = RasterConstructors.asRaster(geom, raster, "d");
+
actual = MapAlgebra.bandAsArray(rasterized, 1);
+
expected =
new double[] {
- 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 1.0, 0.0, 0.0, 0.0,
- 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0,
0.0, 0.0, 0.0, 0.0, 1.0,
- 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
1.0, 0.0, 0.0, 0.0, 0.0,
- 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
+ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0,
+ 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0,
0.0, 0.0, 0.0, 0.0, 0.0,
+ 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 1.0, 0.0, 0.0, 0.0,
+ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 1.0, 0.0,
+ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0
};
assertArrayEquals(expected, actual, 0.1d);
// LinearRing
- geom = Constructors.geomFromWKT("LINEARRING (10 10, 18 20, 15 24, 24 25,
10 10)", 0);
- rasterized = RasterConstructors.asRaster(geom, raster, "d", 3093151, 3d);
+ geom = Constructors.geomFromWKT("LINEARRING (10 -10, 18 -20, 15 -24, 24
-25, 10 -10)", 0);
+ rasterized = RasterConstructors.asRaster(geom, raster, "d", false,
3093151, 3d);
actual = MapAlgebra.bandAsArray(rasterized, 1);
+
expected =
new double[] {
- 0.0, 0.0, 0.0, 0.0, 3093151.0, 3093151.0, 3093151.0, 0.0, 0.0,
3093151.0, 3093151.0, 0.0,
- 3093151.0, 0.0, 0.0, 0.0, 0.0, 3093151.0, 3093151.0, 0.0, 0.0, 0.0,
0.0, 0.0, 3093151.0,
- 3093151.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3093151.0, 0.0, 0.0, 0.0, 0.0,
0.0, 3093151.0, 0.0,
- 0.0, 0.0, 0.0, 0.0, 3093151.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3093151.0,
0.0, 0.0, 0.0, 0.0,
- 0.0, 0.0
+ 3093151.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3093151.0, 3093151.0,
0.0, 0.0, 0.0, 0.0,
+ 0.0, 0.0, 0.0, 3093151.0, 3093151.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 3093151.0,
+ 3093151.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3093151.0, 3093151.0,
0.0, 0.0, 0.0, 0.0,
+ 0.0, 0.0, 0.0, 3093151.0, 3093151.0, 0.0, 0.0, 0.0, 0.0, 0.0,
3093151.0, 3093151.0,
+ 3093151.0, 3093151.0, 0.0, 0.0, 0.0, 0.0, 3093151.0, 3093151.0,
3093151.0, 3093151.0,
+ 3093151.0
};
assertArrayEquals(expected, actual, 0.1d);
rasterized = RasterConstructors.asRaster(geom, raster, "d");
actual = MapAlgebra.bandAsArray(rasterized, 1);
+
expected =
new double[] {
- 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0,
0.0, 0.0, 0.0, 0.0, 1.0,
- 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0,
1.0, 0.0, 0.0, 0.0, 0.0,
- 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0,
1.0, 0.0, 0.0, 0.0, 0.0,
- 0.0, 0.0
+ 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 1.0,
+ 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 1.0,
+ 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0,
0.0, 0.0, 1.0, 1.0, 1.0,
+ 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0
};
+
assertArrayEquals(expected, actual, 0.1d);
// MultiPoints
- geom = Constructors.geomFromWKT("MULTIPOINT ((5 5), (10 10), (15 15))", 0);
- rasterized = RasterConstructors.asRaster(geom, raster, "d", 3093151, 3d);
+ geom = Constructors.geomFromWKT("MULTIPOINT ((5 -5), (10 -10), (15 -15))",
0);
+ rasterized = RasterConstructors.asRaster(geom, raster, "d", false,
3093151, 3d);
actual = MapAlgebra.bandAsArray(rasterized, 1);
+
expected =
new double[] {
- 0.0, 0.0, 0.0, 0.0, 3093151.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
3093151.0, 0.0, 0.0,
- 0.0, 0.0, 0.0, 0.0, 0.0, 3093151.0, 0.0, 0.0, 0.0, 0.0
+ 3093151.0, 3093151.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3093151.0, 3093151.0,
0.0, 0.0, 0.0, 0.0,
+ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3093151.0,
0.0, 0.0, 0.0, 0.0, 0.0,
+ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3093151.0,
3093151.0, 0.0, 0.0, 0.0,
+ 0.0, 0.0, 3093151.0, 3093151.0
};
assertArrayEquals(expected, actual, 0.1d);
rasterized = RasterConstructors.asRaster(geom, raster, "d");
actual = MapAlgebra.bandAsArray(rasterized, 1);
+
expected =
new double[] {
- 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0,
0.0, 0.0, 0.0, 0.0, 0.0,
- 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0
+ 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0,
+ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0,
+ 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0
};
+
assertArrayEquals(expected, actual, 0.1d);
// Point
- geom = Constructors.geomFromWKT("POINT (5 5)", 0);
- rasterized = RasterConstructors.asRaster(geom, raster, "d", 3093151, 3d);
+ geom = Constructors.geomFromWKT("POINT (5 -5)", 0);
+ rasterized = RasterConstructors.asRaster(geom, raster, "d", false,
3093151, 3d);
+
actual = MapAlgebra.bandAsArray(rasterized, 1);
- expected = new double[] {3093151.0};
+
+ expected = new double[] {3093151.0, 3093151.0, 3093151.0, 3093151.0};
assertArrayEquals(expected, actual, 0.1d);
rasterized = RasterConstructors.asRaster(geom, raster, "d");
actual = MapAlgebra.bandAsArray(rasterized, 1);
- expected = new double[] {1.0};
+ expected = new double[] {1.0, 1.0, 1.0, 1.0};
assertArrayEquals(expected, actual, 0.1d);
}
@@ -214,15 +240,15 @@ public class RasterConstructorsTest extends
RasterTestBase {
// Horizontal LineString
GridCoverage2D raster =
RasterConstructors.makeEmptyRaster(2, 255, 255, 1, -1, 2, -2, 0, 0,
4326);
- Geometry geom = Constructors.geomFromEWKT("LINESTRING(1 1, 2 1, 10 1)");
- GridCoverage2D rasterized = RasterConstructors.asRaster(geom, raster, "d",
3093151, 0d);
+ Geometry geom = Constructors.geomFromEWKT("LINESTRING(1 -1, 2 -1, 10 -1)");
+ GridCoverage2D rasterized = RasterConstructors.asRaster(geom, raster, "d",
false, 3093151, 0d);
double[] actual = MapAlgebra.bandAsArray(rasterized, 1);
double[] expected = new double[] {3093151.0, 3093151.0, 3093151.0,
3093151.0, 3093151.0};
assertArrayEquals(expected, actual, 0.1d);
// Vertical LineString
- geom = Constructors.geomFromEWKT("LINESTRING(1 1, 1 2, 1 10)");
- rasterized = RasterConstructors.asRaster(geom, raster, "d", 3093151, 0d);
+ geom = Constructors.geomFromEWKT("LINESTRING(1 -1, 1 -2, 1 -10)");
+ rasterized = RasterConstructors.asRaster(geom, raster, "d", false,
3093151, 0d);
actual = MapAlgebra.bandAsArray(rasterized, 1);
expected = new double[] {3093151.0, 3093151.0, 3093151.0, 3093151.0,
3093151.0};
assertArrayEquals(expected, actual, 0.1d);
@@ -235,34 +261,38 @@ public class RasterConstructorsTest extends
RasterTestBase {
rasterFromGeoTiff(resourceFolder +
"raster/raster_with_no_data/test5.tiff");
Geometry geom =
Constructors.geomFromWKT("POLYGON((1.5 1.5, 3.8 3.0, 4.5 4.4, 3.4 3.5,
1.5 1.5))", 0);
- GridCoverage2D rasterized = RasterConstructors.asRaster(geom, raster, "d",
612028, 5d);
+ GridCoverage2D rasterized = RasterConstructors.asRaster(geom, raster, "d",
true, 612028, 5d);
double[] actual = Arrays.stream(MapAlgebra.bandAsArray(rasterized,
1)).toArray();
double[] expected = {
- 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0,
- 0.0, 0.0, 0.0, 612028.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 612028.0, 612028.0,
- 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 612028.0, 612028.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0,
- 0.0, 0.0, 0.0, 612028.0, 612028.0, 612028.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 612028.0,
- 612028.0, 612028.0, 612028.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
612028.0, 612028.0, 612028.0,
- 612028.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 612028.0, 612028.0,
612028.0, 0.0, 0.0, 0.0, 0.0,
- 0.0, 0.0, 0.0, 0.0, 612028.0, 612028.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 612028.0,
- 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 612028.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0,
- 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
+ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 612028.0,
612028.0, 0.0, 0.0, 0.0,
+ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 612028.0, 612028.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0,
+ 0.0, 0.0, 0.0, 612028.0, 612028.0, 612028.0, 612028.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0,
+ 0.0, 612028.0, 612028.0, 612028.0, 612028.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0,
+ 612028.0, 612028.0, 612028.0, 612028.0, 612028.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0,
+ 612028.0, 612028.0, 612028.0, 612028.0, 612028.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0,
+ 612028.0, 612028.0, 612028.0, 612028.0, 612028.0, 612028.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0,
+ 612028.0, 612028.0, 612028.0, 612028.0, 612028.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0,
+ 612028.0, 612028.0, 612028.0, 612028.0, 612028.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0,
+ 612028.0, 612028.0, 612028.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 612028.0,
+ 612028.0, 612028.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
612028.0, 612028.0, 0.0,
+ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 612028.0, 612028.0,
0.0, 0.0, 0.0, 0.0, 0.0,
+ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
};
assertArrayEquals(expected, actual, 0.1d);
- rasterized = RasterConstructors.asRaster(geom, raster, "d", 5484);
+ rasterized = RasterConstructors.asRaster(geom, raster, "d", false, 5484);
actual = Arrays.stream(MapAlgebra.bandAsArray(rasterized, 1)).toArray();
expected =
new double[] {
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0,
- 0.0, 0.0, 0.0, 0.0, 5484.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 5484.0,
- 5484.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 5484.0, 5484.0,
0.0, 0.0, 0.0, 0.0,
- 0.0, 0.0, 0.0, 0.0, 0.0, 5484.0, 5484.0, 5484.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0,
- 5484.0, 5484.0, 5484.0, 5484.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
5484.0, 5484.0, 5484.0,
- 5484.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 5484.0, 5484.0, 5484.0,
0.0, 0.0, 0.0, 0.0,
- 0.0, 0.0, 0.0, 0.0, 5484.0, 5484.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 5484.0,
- 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 5484.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0,
- 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0
+ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 5484.0, 5484.0,
+ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 5484.0, 5484.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0,
+ 0.0, 0.0, 0.0, 5484.0, 5484.0, 5484.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 5484.0,
+ 5484.0, 5484.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 5484.0,
5484.0, 5484.0, 0.0, 0.0,
+ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 5484.0, 5484.0, 5484.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0,
+ 0.0, 5484.0, 5484.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
5484.0, 5484.0, 0.0,
+ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 5484.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0,
+ 0.0, 0.0, 5484.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0
};
assertArrayEquals(expected, actual, 0.1d);
}
@@ -275,7 +305,7 @@ public class RasterConstructorsTest extends RasterTestBase {
Constructors.geomFromWKT(
"POLYGON((1.5 1.5, 3.8 3.0, 4.5 4.4, 3.4 3.5, 1.5 1.5))",
RasterAccessors.srid(raster));
GridCoverage2D rasterized =
- RasterConstructors.asRasterWithRasterExtent(geom, raster, "d", 612028,
5d);
+ RasterConstructors.asRasterWithRasterExtent(geom, raster, "d", false,
612028, 5d);
int widthActual = RasterAccessors.getWidth(rasterized);
int widthExpected = RasterAccessors.getWidth(raster);
@@ -293,7 +323,7 @@ public class RasterConstructorsTest extends RasterTestBase {
Geometry geom =
Constructors.geomFromWKT("POLYGON((0.1 0.1, 0.1 0.4, 0.4 0.4, 0.4 0.1,
0.1 0.1))", 0);
GridCoverage2D rasterized =
- RasterConstructors.asRasterWithRasterExtent(geom, raster, "d", 100d,
0d);
+ RasterConstructors.asRasterWithRasterExtent(geom, raster, "d", false,
100d, 0d);
assertEquals(0, rasterized.getEnvelope2D().x, 1e-6);
assertEquals(0, rasterized.getEnvelope2D().y, 1e-6);
assertEquals(0.5, rasterized.getEnvelope2D().width, 1e-6);
diff --git a/docs/api/sql/Raster-operators.md b/docs/api/sql/Raster-operators.md
index e874888b91..e456cc0693 100644
--- a/docs/api/sql/Raster-operators.md
+++ b/docs/api/sql/Raster-operators.md
@@ -1160,6 +1160,11 @@ Output:
Introduction: This returns a statistic value specified by `statType` over the
region of interest defined by `zone`. It computes the statistic from the pixel
values within the ROI geometry and returns the result. If the `excludeNoData`
parameter is not specified, it will default to `true`. This excludes NoData
values from the statistic calculation. Additionally, if the `band` parameter is
not provided, band 1 will be used by default for the statistic computation. The
valid options for `st [...]
+The `allTouched` parameter (Since `v1.7.1`) determines how pixels are selected:
+
+- When true, any pixel touched by the geometry will be included.
+- When false (default), only pixels whose centroid intersects with the
geometry will be included.
+
- `count`: Number of pixels in the region.
- `sum`: Sum of pixel values.
- `mean|average|avg`: Arithmetic mean.
@@ -1183,15 +1188,19 @@ Introduction: This returns a statistic value specified
by `statType` over the re
Format:
```
-RS_ZonalStats(raster: Raster, zone: Geometry, band: Integer, statType: String,
excludeNoData: Boolean, lenient: Boolean)
+RS_ZonalStats(raster: Raster, zone: Geometry, band: Integer, statType: String,
allTouched: Boolean, excludeNoData: Boolean, lenient: Boolean)
```
```
-RS_ZonalStats(raster: Raster, zone: Geometry, band: Integer, statType: String,
excludeNoData: Boolean)
+RS_ZonalStats(raster: Raster, zone: Geometry, band: Integer, statType: String,
allTouched: Boolean, excludeNoData: Boolean)
```
```
-RS_ZonalStats(raster: Raster, zone: Geometry, band: Integer, statType: String)
+RS_ZonalStats(raster: Raster, zone: Geometry, band: Integer, statType: String,
allTouched: Boolean)
+```
+
+```
+RS_ZonalStats(raster: Raster, zone: Geometry, statType: String, allTouched:
Boolean)
```
```
@@ -1203,7 +1212,7 @@ Since: `v1.5.1`
SQL Example
```sql
-RS_ZonalStats(rast1, geom1, 1, 'sum', false)
+RS_ZonalStats(rast1, geom1, 1, 'sum', true, false)
```
Output:
@@ -1215,7 +1224,7 @@ Output:
SQL Example
```sql
-RS_ZonalStats(rast2, geom2, 1, 'mean', true)
+RS_ZonalStats(rast2, geom2, 1, 'mean', false, true)
```
Output:
@@ -1228,6 +1237,11 @@ Output:
Introduction: Returns a struct of statistic values, where each statistic is
computed over a region defined by the `zone` geometry. The struct has the
following schema:
+The `allTouched` parameter (Since `v1.7.1`) determines how pixels are selected:
+
+- When true, any pixel touched by the geometry will be included.
+- When false (default), only pixels whose centroid intersects with the
geometry will be included.
+
- count: Count of the pixels.
- sum: Sum of the pixel values.
- mean: Arithmetic mean.
@@ -1251,11 +1265,15 @@ Introduction: Returns a struct of statistic values,
where each statistic is comp
Format:
```
-RS_ZonalStatsAll(raster: Raster, zone: Geometry, band: Integer, excludeNodata:
Boolean, lenient: Boolean)
+RS_ZonalStatsAll(raster: Raster, zone: Geometry, band: Integer, allTouched:
Boolean, excludeNodata: Boolean, lenient: Boolean)
```
```
-RS_ZonalStatsAll(raster: Raster, zone: Geometry, band: Integer, excludeNodata:
Boolean)
+RS_ZonalStatsAll(raster: Raster, zone: Geometry, band: Integer, allTouched:
Boolean, excludeNodata: Boolean)
+```
+
+```
+RS_ZonalStatsAll(raster: Raster, zone: Geometry, band: Integer, allTouched:
Boolean)
```
```
@@ -1271,7 +1289,7 @@ Since: `v1.5.1`
SQL Example
```sql
-RS_ZonalStatsAll(rast1, geom1, 1, false)
+RS_ZonalStatsAll(rast1, geom1, 1, true, false)
```
Output:
@@ -1283,7 +1301,7 @@ Output:
SQL Example
```sql
-RS_ZonalStatsAll(rast2, geom2, 1, true)
+RS_ZonalStatsAll(rast2, geom2, 1, false, true)
```
Output:
@@ -1446,6 +1464,11 @@ Introduction: Returns a raster that is clipped by the
given geometry.
If `crop` is not specified then it will default to `true`, meaning it will
make the resulting raster shrink to the geometry's extent and if `noDataValue`
is not specified then the resulting raster will have the minimum possible value
for the band pixel data type.
+The `allTouched` parameter (Since `v1.7.1`) determines how pixels are selected:
+
+- When true, any pixel touched by the geometry will be included.
+- When false (default), only pixels whose centroid intersects with the
geometry will be included.
+
!!!Note
- Since `v1.5.1`, if the coordinate reference system (CRS) of the input
`geom` geometry differs from that of the `raster`, then `geom` will be
transformed to match the CRS of the `raster`. If the `raster` or `geom` doesn't
have a CRS then it will default to `4326/WGS84`.
- Since `v1.7.0`, `RS_Clip` function will return `null` if the `raster`
and `geometry` geometry do not intersect. If you want to throw an exception in
this case, you can set the `lenient` parameter to `false`.
@@ -1453,15 +1476,19 @@ If `crop` is not specified then it will default to
`true`, meaning it will make
Format:
```
-RS_Clip(raster: Raster, band: Integer, geom: Geometry, noDataValue: Double,
crop: Boolean, lenient: Boolean)
+RS_Clip(raster: Raster, band: Integer, geom: Geometry, allTouched: Boolean,
noDataValue: Double, crop: Boolean, lenient: Boolean)
+```
+
+```
+RS_Clip(raster: Raster, band: Integer, geom: Geometry, allTouched: Boolean,
noDataValue: Double, crop: Boolean)
```
```
-RS_Clip(raster: Raster, band: Integer, geom: Geometry, noDataValue: Double,
crop: Boolean)
+RS_Clip(raster: Raster, band: Integer, geom: Geometry, allTouched: Boolean,
noDataValue: Double)
```
```
-RS_Clip(raster: Raster, band: Integer, geom: Geometry, noDataValue: Double)
+RS_Clip(raster: Raster, band: Integer, geom: Geometry, allTouched: Boolean)
```
```
@@ -1480,7 +1507,7 @@ SQL Example
SELECT RS_Clip(
RS_FromGeoTiff(content), 1,
ST_GeomFromWKT('POLYGON ((236722 4204770, 243900 4204770, 243900
4197590, 221170 4197590, 236722 4204770))'),
- 200, true
+ false, 200, true
)
```
@@ -1494,7 +1521,7 @@ SQL Example
SELECT RS_Clip(
RS_FromGeoTiff(content), 1,
ST_GeomFromWKT('POLYGON ((236722 4204770, 243900 4204770, 243900
4197590, 221170 4197590, 236722 4204770))'),
- 200, false
+ false, 200, false
)
```
@@ -2039,14 +2066,14 @@ to this function.
!!!Note
Since `v1.5.1`, if the coordinate reference system (CRS) of the input
`geom` geometry differs from that of the `raster`, then `geom` will be
transformed to match the CRS of the `raster`. If the `raster` or `geom` doesn't
have a CRS then it will default to `4326/WGS84`.
-Format:
+Format without ROI `geom`:
```
RS_SetValues(raster: Raster, bandIndex: Integer, colX: Integer, rowY: Integer,
width: Integer, height: Integer, newValues: ARRAY[Double], keepNoData: Boolean
= false)
```
```
-RS_SetValues(raster: Raster, bandIndex: Integer, geom: Geometry, newValue:
Double, keepNoData: Boolean = false)
+RS_SetValues(raster: Raster, bandIndex: Integer, colX: Integer, rowY: Integer,
width: Integer, height: Integer, newValues: ARRAY[Double])
```
Since: `v1.5.0`
@@ -2056,12 +2083,31 @@ set to the corresponding value in `newValues`. The
`newValues` should be provide
The geometry variant of this function accepts all types of Geometries, and it
sets the `newValue` in the specified region under the `geom`.
+The `allTouched` parameter (Since `v1.7.1`) determines how pixels are selected:
+
+- When true, any pixel touched by the geometry will be included.
+- When false (default), only pixels whose centroid intersects with the
geometry will be included.
+
!!!note
If the shape of `newValues` doesn't match with provided `width` and
`height`, `IllegalArgumentException` is thrown.
!!!Note
If the mentioned `bandIndex` doesn't exist, this will throw an
`IllegalArgumentException`.
+Format with ROI `geom`:
+
+```
+RS_SetValues(raster: Raster, bandIndex: Integer, geom: Geometry, newValue:
Double, allTouched: Boolean = false, keepNoData: Boolean = false)
+```
+
+```
+RS_SetValues(raster: Raster, bandIndex: Integer, geom: Geometry, newValue:
Double, allTouched: Boolean = false)
+```
+
+```
+RS_SetValues(raster: Raster, bandIndex: Integer, geom: Geometry, newValue:
Double)
+```
+
SQL Example
```sql
@@ -2091,7 +2137,7 @@ SELECT RS_BandAsArray(
RS_MakeEmptyRaster(1, 5, 5, 1, -1, 1, -1, 0, 0, 0),
Array(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), 1
),
- 1, ST_GeomFromWKT('POLYGON((1 -1, 3 -3, 6 -6, 4 -1, 1 -1))'), 255,
false
+ 1, ST_GeomFromWKT('POLYGON((1 -1, 3 -3, 6 -6, 4 -1, 1 -1))'), 255,
false, false
)
)
```
diff --git a/docs/api/sql/Raster-writer.md b/docs/api/sql/Raster-writer.md
index c4c1d57195..df709ebc96 100644
--- a/docs/api/sql/Raster-writer.md
+++ b/docs/api/sql/Raster-writer.md
@@ -245,21 +245,32 @@ The newly created DataFrame can be written to disk again
but must be under a dif
#### RS_AsRaster
-Introduction: Converts a Geometry to a Raster dataset. Defaults to using `1.0`
for cell `value` and `null` for `noDataValue` if not provided. Supports all
geometry types.
-The `pixelType` argument defines data type of the output raster. This can be
one of the following, D (double), F (float), I (integer), S (short), US
(unsigned short) or B (byte).
-The `useGeometryExtent` argument defines the extent of the resultant raster.
When set to `true`, it corresponds to the extent of `geom`, and when set to
false, it corresponds to the extent of `raster`. Default value is `true` if not
set.
+Introduction: `RS_AsRaster` converts a vector geometry into a raster dataset
by assigning a specified value to all pixels covered by the geometry. Unlike
`RS_Clip`, which extracts a subset of an existing raster while preserving its
original values, `RS_AsRaster` generates a new raster where the geometry is
rasterized onto a raster grid. The function supports all geometry types and
takes the following parameters:
+
+* `geom`: The geometry to be rasterized.
+* `raster`: The reference raster to be used for overlaying the `geom` on.
+* `pixelType`: Defines data type of the output raster. This can be one of the
following, D (double), F (float), I (integer), S (short), US (unsigned short)
or B (byte).
+* `allTouched` (Since: `v1.7.1`): Decides the pixel selection criteria. If set
to `true`, the function selects all pixels touched by the geometry, else,
selects only pixels who's centroids intersect the geometry. Defaults to `false`.
+* `Value`: The value to be used for assigning pixels covered by the geometry.
Defaults to using `1.0` for cell `value` if not provided.
+* `noDataValue`: Used for assigning the no data value of the resultant raster.
Defaults to `null` if not provided.
+* `useGeometryExtent`: Defines the extent of the resultant raster. When set to
`true`, it corresponds to the extent of `geom`, and when set to false, it
corresponds to the extent of `raster`. Default value is `true` if not set.
+
Format:
```
-RS_AsRaster(geom: Geometry, raster: Raster, pixelType: String, value: Double,
noDataValue: Double, useGeometryExtent: Boolean)
+RS_AsRaster(geom: Geometry, raster: Raster, pixelType: String, allTouched:
Boolean, value: Double, noDataValue: Double, useGeometryExtent: Boolean)
+```
+
+```
+RS_AsRaster(geom: Geometry, raster: Raster, pixelType: String, allTouched:
Boolean, value: Double, noDataValue: Double)
```
```
-RS_AsRaster(geom: Geometry, raster: Raster, pixelType: String, value: Double,
noDataValue: Double)
+RS_AsRaster(geom: Geometry, raster: Raster, pixelType: String, allTouched:
Boolean, value: Double)
```
```
-RS_AsRaster(geom: Geometry, raster: Raster, pixelType: String, value: Double)
+RS_AsRaster(geom: Geometry, raster: Raster, pixelType: String, allTouched:
Boolean)
```
```
@@ -286,7 +297,7 @@ SQL Example
SELECT RS_AsRaster(
ST_GeomFromWKT('POLYGON((15 15, 18 20, 15 24, 24 25, 15 15))'),
RS_MakeEmptyRaster(2, 255, 255, 3, -215, 2, -2, 0, 0, 4326),
- 'D', 255.0, 0d
+ 'D', false, 255.0, 0d
)
```
@@ -316,7 +327,7 @@ GridCoverage2D["g...
SELECT RS_AsRaster(
ST_GeomFromWKT('POLYGON((15 15, 18 20, 15 24, 24 25, 15 15))'),
RS_MakeEmptyRaster(2, 255, 255, 3, 215, 2, -2, 0, 0, 0),
- 'D',255, 0d, false
+ 'D', true, 255, 0d, false
)
```
diff --git
a/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/PixelFunctionEditors.scala
b/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/PixelFunctionEditors.scala
index 086213dd7a..330b62e7f9 100644
---
a/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/PixelFunctionEditors.scala
+++
b/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/PixelFunctionEditors.scala
@@ -28,6 +28,7 @@ case class RS_SetValues(inputExpressions: Seq[Expression])
extends InferredExpression(
inferrableFunction8(PixelFunctionEditors.setValues),
inferrableFunction7(PixelFunctionEditors.setValues),
+ inferrableFunction6(PixelFunctionEditors.setValues),
inferrableFunction5(PixelFunctionEditors.setValues),
inferrableFunction4(PixelFunctionEditors.setValues)) {
protected def withNewChildrenInternal(newChildren: IndexedSeq[Expression]) =
{
diff --git
a/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterBandAccessors.scala
b/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterBandAccessors.scala
index 3108be64f9..6ac6c1770e 100644
---
a/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterBandAccessors.scala
+++
b/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterBandAccessors.scala
@@ -50,6 +50,7 @@ case class RS_Count(inputExpressions: Seq[Expression])
case class RS_ZonalStats(inputExpressions: Seq[Expression])
extends InferredExpression(
+ inferrableFunction7(RasterBandAccessors.getZonalStats),
inferrableFunction6(RasterBandAccessors.getZonalStats),
inferrableFunction5(RasterBandAccessors.getZonalStats),
inferrableFunction4(RasterBandAccessors.getZonalStats),
diff --git
a/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterBandEditors.scala
b/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterBandEditors.scala
index b9690115f0..ae39b404b1 100644
---
a/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterBandEditors.scala
+++
b/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterBandEditors.scala
@@ -59,6 +59,7 @@ case class RS_Union(inputExpressions: Seq[Expression])
case class RS_Clip(inputExpressions: Seq[Expression])
extends InferredExpression(
+ inferrableFunction7(RasterBandEditors.clip),
inferrableFunction6(RasterBandEditors.clip),
inferrableFunction5(RasterBandEditors.clip),
inferrableFunction4(RasterBandEditors.clip),
diff --git
a/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterConstructors.scala
b/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterConstructors.scala
index 6119719272..b033c29a67 100644
---
a/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterConstructors.scala
+++
b/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterConstructors.scala
@@ -42,10 +42,11 @@ case class RS_FromArcInfoAsciiGrid(inputExpressions:
Seq[Expression])
case class RS_AsRaster(inputExpressions: Seq[Expression])
extends InferredExpression(
- inferrableFunction5(RasterConstructors.asRaster),
inferrableFunction3(RasterConstructors.asRaster),
inferrableFunction4(RasterConstructors.asRaster),
- inferrableFunction6(RasterConstructors.asRaster)) {
+ inferrableFunction5(RasterConstructors.asRaster),
+ inferrableFunction6(RasterConstructors.asRaster),
+ inferrableFunction7(RasterConstructors.asRaster)) {
protected def withNewChildrenInternal(newChildren: IndexedSeq[Expression]) =
{
copy(inputExpressions = newChildren)
}
diff --git
a/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterFunctions.scala
b/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterFunctions.scala
index 04da5704ad..6b59ff9ff0 100644
---
a/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterFunctions.scala
+++
b/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterFunctions.scala
@@ -182,23 +182,28 @@ case class RS_ZonalStatsAll(inputExpressions:
Seq[Expression])
} else {
1
}
- val noData = if (inputExpressions.length >= 4) {
+ val allTouched = if (inputExpressions.length >= 4) {
inputExpressions(3).eval(input).asInstanceOf[Boolean]
} else {
- true
+ false
}
- val lenient = if (inputExpressions.length >= 5) {
+ val noData = if (inputExpressions.length >= 5) {
inputExpressions(4).eval(input).asInstanceOf[Boolean]
} else {
true
}
+ val lenient = if (inputExpressions.length >= 6) {
+ inputExpressions(5).eval(input).asInstanceOf[Boolean]
+ } else {
+ true
+ }
// Check if the raster geometry is null
if (rasterGeom == null) {
null
} else {
val zonalStatsAll =
- RasterBandAccessors.getZonalStatsAll(rasterGeom, roi, band, noData,
lenient)
+ RasterBandAccessors.getZonalStatsAll(rasterGeom, roi, band,
allTouched, noData, lenient)
// Create an InternalRow with the zonalStatsAll
if (zonalStatsAll == null) {
return null
@@ -220,8 +225,10 @@ case class RS_ZonalStatsAll(inputExpressions:
Seq[Expression])
Seq(RasterUDT, GeometryUDT, IntegerType)
} else if (inputExpressions.length == 4) {
Seq(RasterUDT, GeometryUDT, IntegerType, BooleanType)
- } else if (inputExpressions.length >= 5) {
- Seq(RasterUDT, GeometryUDT, IntegerType, BooleanType)
+ } else if (inputExpressions.length == 5) {
+ Seq(RasterUDT, GeometryUDT, IntegerType, BooleanType, BooleanType)
+ } else if (inputExpressions.length >= 6) {
+ Seq(RasterUDT, GeometryUDT, IntegerType, BooleanType, BooleanType,
BooleanType)
} else {
Seq(RasterUDT, GeometryUDT)
}
diff --git
a/spark/common/src/test/scala/org/apache/sedona/sql/rasterIOTest.scala
b/spark/common/src/test/scala/org/apache/sedona/sql/rasterIOTest.scala
index d256d6a31a..a62b975454 100644
--- a/spark/common/src/test/scala/org/apache/sedona/sql/rasterIOTest.scala
+++ b/spark/common/src/test/scala/org/apache/sedona/sql/rasterIOTest.scala
@@ -91,9 +91,10 @@ class rasterIOTest extends TestBaseScala with BeforeAndAfter
with GivenWhenThen
}
it("Passed RS_AsRaster with empty raster") {
- var df = sparkSession.sql(
- "SELECT RS_MakeEmptyRaster(2, 255, 255, 3, -215, 2, -2, 0, 0, 4326) as
raster, ST_GeomFromWKT('POLYGON((15 15, 18 20, 15 24, 24 25, 15 15))') as geom")
- var rasterized = df.selectExpr("RS_AsRaster(geom, raster, 'd', 255, 0d)
as rasterized")
+ val df = sparkSession.sql(
+ "SELECT RS_MakeEmptyRaster(2, 255, 255, 3, 215, 2, -2, 0, 0, 4326) as
raster, ST_GeomFromWKT('POLYGON((15 15, 18 20, 15 24, 24 25, 15 15))') as geom")
+ var rasterized =
+ df.selectExpr("RS_AsRaster(geom, raster, 'd', false, 255, 0d) as
rasterized")
var actual = rasterized
.selectExpr("RS_BandAsArray(rasterized, 1)")
.first()
@@ -103,7 +104,7 @@ class rasterIOTest extends TestBaseScala with
BeforeAndAfter with GivenWhenThen
"Array(255.0, 255.0, 255.0, 255.0, 0.0, 0.0, 255.0, 255.0, 0.0, 0.0,
0.0, 255.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)"
assertEquals(expected, actual)
- rasterized = df.selectExpr("RS_AsRaster(geom, raster, 'd', 3093151) as
rasterized")
+ rasterized = df.selectExpr("RS_AsRaster(geom, raster, 'd', false,
3093151) as rasterized")
actual = rasterized
.selectExpr("RS_BandAsArray(rasterized, 1)")
.first()
@@ -125,19 +126,20 @@ class rasterIOTest extends TestBaseScala with
BeforeAndAfter with GivenWhenThen
}
it("Passed RS_AsRaster LineString") {
- var df = sparkSession.sql(
- "SELECT RS_MakeEmptyRaster(2, 255, 255, 3, -215, 2, -2, 0, 0, 4326) as
raster, ST_GeomFromWKT('LINESTRING(1 1, 2 1, 10 1)') as geom")
- var rasterized = df.selectExpr("RS_AsRaster(geom, raster, 'd', 255, 0d)
as rasterized")
+ val df = sparkSession.sql(
+ "SELECT RS_MakeEmptyRaster(2, 255, 255, 3, 215, 2, -2, 0, 0, 4326) as
raster, ST_GeomFromWKT('LINESTRING(1 1, 2 1, 10 1)') as geom")
+ var rasterized =
+ df.selectExpr("RS_AsRaster(geom, raster, 'd', false, 255, 0d) as
rasterized")
var actual = rasterized
.selectExpr("RS_BandAsArray(rasterized, 1)")
.first()
.getSeq(0)
.mkString("Array(", ", ", ")")
- var expected = "Array(255.0, 255.0, 255.0, 255.0, 255.0)"
+ var expected = "Array(0.0, 0.0, 0.0, 0.0, 255.0, 255.0, 255.0, 255.0)"
assertEquals(expected, actual)
rasterized = df.selectExpr(
- "RS_AsRaster(ST_GeomFromWKT('LINESTRING(1 1, 1 2, 1 10)'), raster,
'd', 255, 0d) as rasterized")
+ "RS_AsRaster(ST_GeomFromWKT('LINESTRING(4 1, 4 2, 4 10)'), raster,
'd', false, 255, 0d) as rasterized")
actual = rasterized
.selectExpr("RS_BandAsArray(rasterized, 1)")
.first()
@@ -149,9 +151,11 @@ class rasterIOTest extends TestBaseScala with
BeforeAndAfter with GivenWhenThen
it("Passed RS_AsRaster with raster") {
var df = sparkSession.read.format("binaryFile").load(resourceFolder +
"raster/test1.tiff")
- df =
- df.selectExpr("ST_GeomFromWKT('POINT(5 5)') as geom",
"RS_FromGeoTiff(content) as raster")
- var rasterized = df.selectExpr("RS_AsRaster(geom, raster, 'd', 61784,
0d) as rasterized")
+ df = df.selectExpr(
+ "ST_GeomFromText('POINT (-13085817.809482181 3993868.8560156375)',
3857) as geom",
+ "RS_FromGeoTiff(content) as raster")
+ var rasterized =
+ df.selectExpr("RS_AsRaster(geom, raster, 'd', false, 61784, 0d) as
rasterized")
var actual = rasterized
.selectExpr("RS_BandAsArray(rasterized, 1)")
.first()
@@ -160,7 +164,7 @@ class rasterIOTest extends TestBaseScala with
BeforeAndAfter with GivenWhenThen
var expected = "Array(61784.0)"
assertEquals(expected, actual)
- rasterized = df.selectExpr("RS_AsRaster(geom, raster, 'd', 255) as
rasterized")
+ rasterized = df.selectExpr("RS_AsRaster(geom, raster, 'd', false, 255)
as rasterized")
actual = rasterized
.selectExpr("RS_BandAsArray(rasterized, 1)")
.first()
@@ -183,13 +187,13 @@ class rasterIOTest extends TestBaseScala with
BeforeAndAfter with GivenWhenThen
var df = sparkSession.sql(
"SELECT RS_MakeEmptyRaster(2, 255, 255, 3, 215, 2, -2, 0, 0, 0) as
raster, ST_GeomFromWKT('POLYGON((15 15, 18 20, 15 24, 24 25, 15 15))') as geom")
var rasterized =
- df.selectExpr("RS_AsRaster(geom, raster, 'd', 255, 0d, false) as
rasterized")
+ df.selectExpr("RS_AsRaster(geom, raster, 'd', false, 255, 0d, false)
as rasterized")
var actualSeq =
rasterized.selectExpr("RS_BandAsArray(rasterized,
1)").first().getSeq[Double](0)
var actualMax = actualSeq.max
var actualSum = actualSeq.sum
var expectedMax = 255.0d
- var expectedSum = 255.0 * 9
+ var expectedSum = 255.0 * 7
assertEquals(expectedMax, actualMax, 1e-5)
assertEquals(expectedSum, actualSum, 1e-5)
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 ecdf257de9..1f7142a73d 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
@@ -681,7 +681,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")
+ val resultDf = df.selectExpr("RS_SetValues(raster, 1, geom, 10, false,
false) as result")
var actual = resultDf
.selectExpr("RS_Value(result, ST_GeomFromWKT('POINT (-77.9146
37.8916)'))")
@@ -733,7 +733,7 @@ class rasteralgebraTest extends TestBaseScala with
BeforeAndAfter with GivenWhen
"RS_BandAsArray(RS_SetValues(raster, 1, ST_GeomFromWKT('LINESTRING(1
-1, 1 -4, 2 -2, 3 -3, 4 -4, 5 -4, 6 -6)'), 255), 1)")
.first()
.getSeq(0)
- var expected = Seq(255.0, 0.0, 0.0, 0.0, 0.0, 255.0, 255.0, 0.0, 0.0,
0.0, 255.0, 255.0,
+ var expected = Seq(255.0, 0.0, 0.0, 0.0, 0.0, 255.0, 255.0, 0.0, 0.0,
0.0, 255.0, 0.0,
255.0, 0.0, 0.0, 255.0, 0.0, 0.0, 255.0, 255.0, 0.0, 0.0, 0.0, 0.0,
255.0)
assert(expected.equals(actual))
@@ -742,8 +742,8 @@ class rasteralgebraTest extends TestBaseScala with
BeforeAndAfter with GivenWhen
"RS_BandAsArray(RS_SetValues(raster, 1, ST_GeomFromWKT('POINT(2
-2)'), 25), 1)")
.first()
.getSeq(0)
- expected = Seq(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 25.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0,
- 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
+ expected = Seq(25.0, 25.0, 0.0, 0.0, 0.0, 25.0, 25.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0,
+ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
assert(expected.equals(actual))
actual = inputDf
@@ -751,8 +751,8 @@ class rasteralgebraTest extends TestBaseScala with
BeforeAndAfter with GivenWhen
"RS_BandAsArray(RS_SetValues(raster, 1,
ST_GeomFromWKT('MULTIPOINT((2 -2), (2 -1), (3 -3))'), 400), 1)")
.first()
.getSeq(0)
- expected = Seq(0.0, 400.0, 0.0, 0.0, 0.0, 0.0, 400.0, 0.0, 0.0, 0.0,
0.0, 0.0, 400.0, 0.0,
- 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
+ expected = Seq(400.0, 400.0, 0.0, 0.0, 0.0, 400.0, 400.0, 400.0, 0.0,
0.0, 0.0, 400.0,
+ 400.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
assert(expected.equals(actual))
actual = inputDf
@@ -761,7 +761,16 @@ class rasteralgebraTest extends TestBaseScala with
BeforeAndAfter with GivenWhen
.first()
.getSeq(0)
expected = Seq(150.0, 150.0, 150.0, 0.0, 0.0, 0.0, 150.0, 150.0, 150.0,
0.0, 0.0, 0.0,
- 150.0, 150.0, 0.0, 0.0, 0.0, 0.0, 150.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
+ 150.0, 150.0, 0.0, 0.0, 0.0, 0.0, 150.0, 0.0, 0.0, 0.0, 0.0, 0.0,
150.0)
+ assert(expected.equals(actual))
+
+ actual = inputDf
+ .selectExpr(
+ "RS_BandAsArray(RS_SetValues(raster, 1, ST_GeomFromWKT('POLYGON((1
-1, 3 -3, 6 -6, 4 -1, 1 -1))'), 150, true), 1)")
+ .first()
+ .getSeq(0)
+ expected = Seq(150.0, 150.0, 150.0, 150.0, 0.0, 0.0, 150.0, 150.0,
150.0, 0.0, 0.0, 0.0,
+ 150.0, 150.0, 150.0, 0.0, 0.0, 0.0, 150.0, 150.0, 0.0, 0.0, 0.0, 0.0,
150.0)
assert(expected.equals(actual))
}
@@ -955,7 +964,7 @@ class rasteralgebraTest extends TestBaseScala with
BeforeAndAfter with GivenWhen
.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 clippedDf = df.selectExpr("RS_Clip(raster, 1, geom, 200, false) as
clipped")
+ val clippedDf = df.selectExpr("RS_Clip(raster, 1, geom, false, 200,
false) as clipped")
// Extract the metadata for the original raster
val clippedMetadata = df
@@ -996,7 +1005,7 @@ class rasteralgebraTest extends TestBaseScala with
BeforeAndAfter with GivenWhen
.selectExpr(
"RS_FromGeoTiff(content) as raster",
"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 clippedDf = df.selectExpr("RS_Clip(raster, 1, geom, false, 200,
false) as clipped")
val clippedMetadata =
df.selectExpr("RS_Metadata(raster)").first().getStruct(0).toSeq.slice(0, 9)
@@ -1015,7 +1024,7 @@ class rasteralgebraTest extends TestBaseScala with
BeforeAndAfter with GivenWhen
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")
+ val croppedDf = df.selectExpr("RS_Clip(raster, 1, geom, false, 200,
false) as cropped")
actualValues = croppedDf
.selectExpr(
"RS_Values(cropped, " +
@@ -1621,20 +1630,22 @@ class rasteralgebraTest extends TestBaseScala with
BeforeAndAfter with GivenWhen
df = df.selectExpr(
"RS_FromGeoTiff(content) as raster",
"ST_GeomFromWKT('POLYGON ((236722 4204770, 243900 4204770, 243900
4197590, 221170 4197590, 236722 4204770))', 26918) as geom")
- var actual = df.selectExpr("RS_ZonalStats(raster, geom, 1, 'sum',
true)").first().get(0)
- assertEquals(1.0719726e7, actual)
+ var actual =
+ df.selectExpr("RS_ZonalStats(raster, geom, 1, 'sum', false,
true)").first().get(0)
+ assertEquals(1.0896994e7, actual)
- actual = df.selectExpr("RS_ZonalStats(raster, geom, 1, 'count',
false)").first().get(0)
- assertEquals(184792.0, actual)
+ actual =
+ df.selectExpr("RS_ZonalStats(raster, geom, 1, 'count', false,
false)").first().get(0)
+ assertEquals(185953.0, actual)
- actual = df.selectExpr("RS_ZonalStats(raster, geom, 1, 'mean',
false)").first().get(0)
- assertEquals(58.00968656653401, actual)
+ actual = df.selectExpr("RS_ZonalStats(raster, geom, 1, 'mean', true,
false)").first().get(0)
+ assertEquals(58.650240700685295, actual)
actual = df.selectExpr("RS_ZonalStats(raster, geom, 1,
'variance')").first().get(0)
- assertEquals(8491.631251863151, actual)
+ assertEquals(8563.303492088418, actual)
actual = df.selectExpr("RS_ZonalStats(raster, geom,
'sd')").first().get(0)
- assertEquals(92.15004748703687, actual)
+ assertEquals(92.53811912983977, actual)
// Test with a polygon in EPSG:4326
actual = df
@@ -1660,20 +1671,22 @@ class rasteralgebraTest extends TestBaseScala with
BeforeAndAfter with GivenWhen
df = df.selectExpr(
"RS_FromGeoTiff(content) as raster",
"ST_GeomFromWKT('POLYGON((-167.750000 87.750000, -155.250000
87.750000, -155.250000 40.250000, -180.250000 40.250000, -167.750000
87.750000))', 4326) as geom")
- var actual = df.selectExpr("RS_ZonalStats(raster, geom, 1, 'sum',
true)").first().get(0)
- assertEquals(3213526.0, actual)
+ var actual =
+ df.selectExpr("RS_ZonalStats(raster, geom, 1, 'sum', false,
true)").first().get(0)
+ assertEquals(3229013.0, actual)
- actual = df.selectExpr("RS_ZonalStats(raster, geom, 1, 'count',
true)").first().get(0)
- assertEquals(14184.0, actual)
+ actual = df.selectExpr("RS_ZonalStats(raster, geom, 1, 'count', true,
true)").first().get(0)
+ assertEquals(14648.0, actual)
- actual = df.selectExpr("RS_ZonalStats(raster, geom, 1, 'mean',
false)").first().get(0)
- assertEquals(226.49605300253515, actual)
+ actual =
+ df.selectExpr("RS_ZonalStats(raster, geom, 1, 'mean', false,
false)").first().get(0)
+ assertEquals(226.54970883322787, actual)
- actual = df.selectExpr("RS_ZonalStats(raster, geom, 1,
'variance')").first().get(0)
- assertEquals(5606.423398599913, actual)
+ actual = df.selectExpr("RS_ZonalStats(raster, geom, 1, 'variance',
false)").first().get(0)
+ assertEquals(5596.966403503485, actual)
actual = df.selectExpr("RS_ZonalStats(raster, geom,
'sd')").first().get(0)
- assertEquals(74.87605357255357, actual)
+ assertEquals(74.81287592054916, actual)
}
it("Passed RS_ZonalStatsAll") {
@@ -1684,19 +1697,19 @@ class rasteralgebraTest extends TestBaseScala with
BeforeAndAfter with GivenWhen
"RS_FromGeoTiff(content) as raster",
"ST_GeomFromWKT('POLYGON ((-8673439.6642 4572993.5327, -8673155.5737
4563873.2099, -8701890.3259 4562931.7093, -8682522.8735 4572703.8908,
-8673439.6642 4572993.5327))', 3857) as geom")
val actual = df
- .selectExpr("RS_ZonalStatsAll(raster, geom, 1, true)")
+ .selectExpr("RS_ZonalStatsAll(raster, geom, 1, false, true)")
.first()
.getStruct(0)
.toSeq
.slice(0, 9)
- val expected = Seq(184792.0, 1.0719726e7, 58.00968656653401, 0.0, 0.0,
92.15004748703687,
- 8491.631251863151, 0.0, 255.0)
+ val expected = Seq(185953.0, 1.0896994e7, 58.600796975566816, 0.0, 0.0,
92.53811912983977,
+ 8563.303492088418, 0.0, 255.0)
assertTrue(expected.equals(actual))
// Test with a polygon that does not intersect the raster in lenient mode
val actual2 = df
.selectExpr(
- "RS_ZonalStatsAll(raster, ST_GeomFromWKT('POLYGON
((-78.22106647832458748 37.76411511479908967, -78.20183062098976734
37.72863564460374874, -78.18088490966962922 37.76753482276972562,
-78.22106647832458748 37.76411511479908967))'), 1, false)")
+ "RS_ZonalStatsAll(raster, ST_GeomFromWKT('POLYGON
((-78.22106647832458748 37.76411511479908967, -78.20183062098976734
37.72863564460374874, -78.18088490966962922 37.76753482276972562,
-78.22106647832458748 37.76411511479908967))'), 1, true, false)")
.first()
.getStruct(0)
assertNull(actual2)
@@ -1710,13 +1723,13 @@ class rasteralgebraTest extends TestBaseScala with
BeforeAndAfter with GivenWhen
"RS_FromGeoTiff(content) as raster",
"ST_GeomFromWKT('POLYGON((-167.750000 87.750000, -155.250000
87.750000, -155.250000 40.250000, -180.250000 40.250000, -167.750000
87.750000))', 0) as geom")
val actual = df
- .selectExpr("RS_ZonalStatsAll(raster, geom, 1, true)")
+ .selectExpr("RS_ZonalStatsAll(raster, geom, 1, false, true)")
.first()
.getStruct(0)
.toSeq
.slice(0, 9)
- val expected = Seq(14184.0, 3213526.0, 226.55992667794473, 255.0, 255.0,
74.87605357255357,
- 5606.423398599913, 1.0, 255.0)
+ val expected = Seq(14249.0, 3229013.0, 226.61330619692416, 255.0, 255.0,
74.81287592054916,
+ 5596.966403503485, 1.0, 255.0)
assertTrue(expected.equals(actual))
}
