[SEDONA-428] Add RS_ZonalStats & RS_ZonalStatsAll

common/pom.xml

@@ -37,6 +37,10 @@
     </properties>
 
     <dependencies>
+        <dependency>
+            <groupId>org.apache.commons</groupId>
+            <artifactId>commons-math3</artifactId>
+        </dependency>
         <dependency>
             <groupId>org.slf4j</groupId>
             <artifactId>slf4j-api</artifactId>

common/src/main/java/org/apache/sedona/common/raster/RasterBandAccessors.java

@@ -18,16 +18,23 @@
  */
 package org.apache.sedona.common.raster;
 
+import org.apache.commons.lang3.tuple.Pair;
+import org.apache.commons.math3.stat.StatUtils;
+import org.apache.commons.math3.stat.descriptive.DescriptiveStatistics;
+import org.apache.sedona.common.Functions;
 import org.apache.sedona.common.utils.RasterUtils;
 import org.geotools.coverage.GridSampleDimension;
 import org.geotools.coverage.grid.GridCoverage2D;
+import org.locationtech.jts.geom.Geometry;
 import org.opengis.referencing.FactoryException;
 
 import javax.media.jai.RasterFactory;
 import java.awt.image.Raster;
 import java.awt.image.WritableRaster;
+import java.util.ArrayList;
 import java.util.Arrays;
 import java.util.HashMap;
+import java.util.List;
 
 public class RasterBandAccessors {
 
@@ -81,6 +88,190 @@ public static long getCount(GridCoverage2D raster, int band) {
 //        return getCount(raster, 1, excludeNoDataValue);
 //    }
 
+    /**
+     * @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 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) throws FactoryException {
+        List<Object> objects = getStatObjects(raster, roi, band, excludeNoData);
+        DescriptiveStatistics stats = (DescriptiveStatistics) objects.get(0);
+        double[] pixelData = (double[]) objects.get(1);
+
+        // order of stats
+        // count, sum, mean, median, mode, stddev, variance, min, max
+        double[] result = new double[9];
+        result[0] = stats.getN();
+        result[1] = stats.getSum();
+        result[2] = stats.getMean();
+        result[3] = stats.getPercentile(50);
+        result[4] = zonalMode(pixelData);
+        result[5] = stats.getStandardDeviation();
+        result[6] = stats.getVariance();
+        result[7] = stats.getMin();
+        result[8] = stats.getMax();
+
+        return result;
+    }
+
+    /**
+     * @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
+     * @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) throws FactoryException {
+        return getZonalStatsAll(raster, roi, band, true);
+    }
+
+    /**
+     * @param raster Raster to use for computing stats
+     * @param roi Geometry to define the region of interest
+     * @return An array with all the stats for the region, excludeNoData is set to true and band is set to 1
+     * @throws FactoryException
+     */
+    public static double[] getZonalStatsAll(GridCoverage2D raster, Geometry roi) throws FactoryException {
+        return getZonalStatsAll(raster, roi, 1, 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 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) throws FactoryException {
+
+        List<Object> objects = getStatObjects(raster, roi, band, excludeNoData);
+        DescriptiveStatistics stats = (DescriptiveStatistics) objects.get(0);
+        double[] pixelData = (double[]) objects.get(1);
+
+        switch (statType.toLowerCase()) {
+            case "sum":
+                return stats.getSum();
+            case "average":
+            case "avg":
+            case "mean":
+                return stats.getMean();
+            case "count":
+                return stats.getN();
+            case "max":
+                return stats.getMax();
+            case "min":
+                return stats.getMin();
+            case "stddev":
+            case "sd":
+                return stats.getStandardDeviation();
+            case "median":
+                return stats.getPercentile(50);
+            case "mode":
+                return zonalMode(pixelData);
+            case "variance":
+                return stats.getVariance();
+            default:
+                throw new IllegalArgumentException("Please select from the accepted options. Some of the valid options are sum, mean, stddev, etc.");
+        }
+    }
+
+    /**
+     * @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
+     * @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) throws FactoryException {
+        return getZonalStats(raster, roi, band, statType, true);
+    }
+
+    /**
+     * @param raster Raster to use for computing stats
+     * @param roi Geometry to define the region of interest
+     * @param statType Define the statistic to be computed
+     * @return A double precision floating point number representing the requested statistic calculated over the specified region. The excludeNoData is set to true and band is set to 1.
+     * @throws FactoryException
+     */
+    public static double getZonalStats(GridCoverage2D raster, Geometry roi, String statType) throws FactoryException {
+        return getZonalStats(raster, roi, 1, statType, true);
+    }
+
+    /**
+     * @param pixelData An array of double with pixel values
+     * @return Mode of the pixel values. If there is multiple with same occurrence, then the largest value will be returned.
+     */
+    private static double zonalMode(double[] pixelData) {
+        double[] modes = StatUtils.mode(pixelData);
+        return modes[modes.length - 1];
+    }
+
+    /**
+     * An intermediate function to compute zonal statistics
+     * @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 excludeNoData Specifies whether to exclude no-data value or not
+     * @return an object of DescriptiveStatistics and an array of double with pixel data.
+     * @throws FactoryException
+     */
+    private static List<Object> getStatObjects(GridCoverage2D raster, Geometry roi, int band, boolean excludeNoData) throws FactoryException {
+        RasterUtils.ensureBand(raster, band);
+
+        if(RasterAccessors.srid(raster) != roi.getSRID()) {
+            // implicitly converting roi geometry CRS to raster CRS
+            roi = RasterUtils.convertCRSIfNeeded(roi, raster.getCoordinateReferenceSystem());
+            // have to set the SRID as RasterUtils.convertCRSIfNeeded doesn't set it even though the geometry is in raster's CRS
+            roi = Functions.setSRID(roi, RasterAccessors.srid(raster));
+        }
+
+        // checking if the raster contains the geometry
+        if (!RasterPredicates.rsIntersects(raster, roi)) {
+            throw new IllegalArgumentException("The provided geometry is not intersecting the raster. Please provide a geometry that is in the raster's extent.");
+        }
+
+        Raster rasterData = RasterUtils.getRaster(raster.getRenderedImage());
+        String datatype = RasterBandAccessors.getBandType(raster, band);
+        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);
+        Raster rasterziedData = RasterUtils.getRaster(rasterizedGeom.getRenderedImage());
+        int width = RasterAccessors.getWidth(rasterizedGeom), height = RasterAccessors.getHeight(rasterizedGeom);
+        double[] rasterizedPixelData = rasterziedData.getSamples(0, 0, width, height, 0, (double[]) null);
+        double[] rasterPixelData = rasterData.getSamples(0, 0, width, height, band - 1, (double[]) null);
+
+        List<Double> pixelData = new ArrayList<>();
+
+        for (int k = 0; k < rasterPixelData.length; k++) {
+
+            // Pixels with a value of 0 in the rasterizedPixelData are skipped during computation,
+            // as they fall outside the geometry specified by 'roi'.
+            // The region of interest defined by 'roi' contains pixel values of 150,
+            // as initialized when constructing the raster via RasterConstructors.asRasterWithRasterExtent.
+            if (rasterizedPixelData[k] == 0 || excludeNoData && noDataValue != null && rasterPixelData[k] == noDataValue) {
+                continue;
+            } else {
+                pixelData.add(rasterPixelData[k]);
+            }
+        }
+
+        double[] pixelsArray = pixelData.stream().mapToDouble(d -> d).toArray();
+
+        DescriptiveStatistics stats = new DescriptiveStatistics(pixelsArray);
+
+        List<Object> statObjects = new ArrayList<>();
+
+        statObjects.add(stats);
+        statObjects.add(pixelsArray);
+        return statObjects;
+    }
+
     public static double[] getSummaryStats(GridCoverage2D rasterGeom, int band, boolean excludeNoDataValue) {
         RasterUtils.ensureBand(rasterGeom, band);
         Raster raster = RasterUtils.getRaster(rasterGeom.getRenderedImage());

common/src/main/java/org/apache/sedona/common/raster/RasterConstructors.java

@@ -26,6 +26,7 @@
 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.DefaultEngineeringCRS;
 import org.geotools.referencing.operation.transform.AffineTransform2D;
@@ -41,6 +42,9 @@
 import javax.media.jai.RasterFactory;
 import java.awt.image.WritableRaster;
 import java.io.IOException;
+import java.util.ArrayList;
+import java.util.Arrays;
+import java.util.List;
 
 public class RasterConstructors
 {
@@ -67,6 +71,43 @@ public static GridCoverage2D fromGeoTiff(byte[] bytes) throws IOException {
     */
     public static GridCoverage2D asRaster(Geometry geom, GridCoverage2D raster, String pixelType, double value, Double noDataValue) throws FactoryException {
 
+        List<Object> objects = rasterization(geom, raster, pixelType, value, noDataValue, true);
+        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.
+    }
+
+    /**
+     * 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.
+     *
+     * @return Rasterized Geometry
+     * @throws FactoryException
+     */
+    public static GridCoverage2D asRaster(Geometry geom, GridCoverage2D raster, String pixelType) throws FactoryException {
+        return asRaster(geom, raster, pixelType, 1, null);
+    }
+
+    /**
+     * 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 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, double value) throws FactoryException {
+        return asRaster(geom, raster, pixelType, value, null);
+    }
+
+    private static List<Object> rasterization(Geometry geom, GridCoverage2D raster, String pixelType, double value, Double noDataValue, boolean useGeometryExtent) throws FactoryException {
         DefaultFeatureCollection featureCollection = getFeatureCollection(geom, raster.getCoordinateReferenceSystem());
 
         double[] metadata = RasterAccessors.metadata(raster);
@@ -89,7 +130,14 @@ public static GridCoverage2D asRaster(Geometry geom, GridCoverage2D raster, Stri
             throw new IllegalArgumentException(String.format("SkewY %d of the raster is not zero.", metadata[7]));
         }
 
-        Envelope2D bound = JTS.getEnvelope2D(geom.getEnvelopeInternal(), raster.getCoordinateReferenceSystem2D());
+        Envelope2D bound = null;
+
+        if (useGeometryExtent) {
+            bound = JTS.getEnvelope2D(geom.getEnvelopeInternal(), raster.getCoordinateReferenceSystem2D());
+        } else {
+            ReferencedEnvelope envelope = ReferencedEnvelope.create(raster.getEnvelope(), raster.getCoordinateReferenceSystem());
+            bound = JTS.getEnvelope2D(envelope, raster.getCoordinateReferenceSystem2D());
+        }
 
         double scaleX = Math.abs(metadata[4]), scaleY = Math.abs(metadata[5]);
         int width = (int) bound.getWidth(), height = (int) bound.getHeight();
@@ -117,37 +165,32 @@ public static GridCoverage2D asRaster(Geometry geom, GridCoverage2D raster, Stri
         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 RasterUtils.clone(writableRaster, rasterized.getSampleDimensions(), rasterized, noDataValue, false); //no need to original raster metadata since this is a new raster.
-    }
-
-    /**
-     * 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.
-     *
-     * @return Rasterized Geometry
-     * @throws FactoryException
-     */
-    public static GridCoverage2D asRaster(Geometry geom, GridCoverage2D raster, String pixelType) throws FactoryException {
-        return asRaster(geom, raster, pixelType, 1, null);
+        return objects;
     }
 
     /**
-     * Returns a raster that is converted from the geometry provided. A convenience function for asRaster.
-     *
+     * For internal use only!
+     * Returns a raster that is converted from the geometry provided with the extent of the reference raster.
      * @param geom The geometry to convert
      * @param raster The reference raster
-     * @param pixelType The data type of pixel/cell of resultant raster.
+     * @param pixelType The data type of pixel/cell of resultant raster
      * @param value The value of the pixel of the resultant raster
+     * @param noDataValue The noDataValue of the resultant raster
      *
-     * @return Rasterized Geometry
+     * @return Rasterized Geometry with reference raster's extent
      * @throws FactoryException
      */
-    public static GridCoverage2D asRaster(Geometry geom, GridCoverage2D raster, String pixelType, double value) throws FactoryException {
-        return asRaster(geom, raster, pixelType, value, null);
+    public static GridCoverage2D asRasterWithRasterExtent(Geometry geom, GridCoverage2D raster, String pixelType, double value, Double noDataValue) throws FactoryException {
+        List<Object> objects = rasterization(geom, raster, pixelType, value, noDataValue, false);
+        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.
     }
 
     public static DefaultFeatureCollection getFeatureCollection(Geometry geom, CoordinateReferenceSystem crs) {