[GEOS-9068] Allow SLDService to limit classification to a given numbe…

…r of standard deviations from the average

doc/en/user/source/extensions/sldservice/index.rst

@@ -205,6 +205,10 @@ The parameters usable to customize the ColorMap are:
      - allows to run the classification on a specific bounding box. Recommended when the overall dataset is too big, and the classification can be performed on a smaller dataset, or to enhance the visualization of a particular subset of data
      - same syntax as WMS/WFS, expected axis order is east/north unless the spatial reference system is explicitly provided, ``minx,miny,max,maxy[,srsName]``
      - 
+   * - stddevs
+     - limits the data the classifier is working on to a range of "stddevs" standard deviations around the mean value. 
+     - a positive floating point number (e.g., '1', '2.5', '3').
+     -
 
 Examples
 ~~~~~~~~~~

src/extension/sldService/src/main/java/org/geoserver/sldservice/rest/ClassifierController.java

@@ -54,6 +54,8 @@
 import org.geotools.data.util.NullProgressListener;
 import org.geotools.factory.CommonFactoryFinder;
 import org.geotools.feature.FeatureCollection;
+import org.geotools.feature.visitor.CalcResult;
+import org.geotools.feature.visitor.StandardDeviationVisitor;
 import org.geotools.filter.function.RangedClassifier;
 import org.geotools.geometry.jts.ReferencedEnvelope;
 import org.geotools.image.util.ImageUtilities;
@@ -68,6 +70,7 @@
 import org.geotools.styling.Style;
 import org.geotools.styling.StyledLayerDescriptor;
 import org.geotools.util.Converters;
+import org.geotools.util.NumberRange;
 import org.geotools.util.factory.Hints;
 import org.geotools.util.logging.Logging;
 import org.geotools.xml.styling.SLDTransformer;
@@ -81,6 +84,7 @@
 import org.opengis.feature.type.PropertyDescriptor;
 import org.opengis.filter.Filter;
 import org.opengis.filter.FilterFactory2;
+import org.opengis.filter.PropertyIsBetween;
 import org.opengis.filter.spatial.BBOX;
 import org.opengis.referencing.FactoryException;
 import org.opengis.referencing.operation.TransformException;
@@ -163,6 +167,7 @@ public Object classify(
             @RequestParam(value = "continuous", required = false, defaultValue = "false")
                     boolean continuous,
             @RequestParam(value = "bbox", required = false) ReferencedEnvelope bbox,
+            @RequestParam(value = "stddevs", required = false) Double stddevs,
             final HttpServletResponse response)
             throws Exception {
         LayerInfo layerInfo = catalog.getLayerByName(layerName);
@@ -173,6 +178,10 @@ public Object classify(
         if (bbox != null && bbox.getCoordinateReferenceSystem() == null) {
             bbox = new ReferencedEnvelope(bbox, layerInfo.getResource().getCRS());
         }
+        if (stddevs != null && stddevs <= 0) {
+            throw new RestException(
+                    "stddevs must be a positive floating point number", HttpStatus.BAD_REQUEST);
+        }
         if (cachingTime > 0) {
             response.setHeader(
                     "cache-control",
@@ -208,7 +217,8 @@ public Object classify(
                                 pointSize,
                                 (FeatureTypeInfo) obj,
                                 ramp,
-                                bbox);
+                                bbox,
+                                stddevs);
             } else if (obj instanceof CoverageInfo) {
                 rules =
                         getRasterRules(
@@ -223,7 +233,8 @@ public Object classify(
                                 (CoverageInfo) obj,
                                 ramp,
                                 continuous,
-                                bbox);
+                                bbox,
+                                stddevs);
             } else {
                 throw new RestException(
                         "The classifier can only work against vector or raster data, "
@@ -369,7 +380,8 @@ private List<Rule> getRasterRules(
             CoverageInfo coverageInfo,
             ColorRamp ramp,
             boolean continuous,
-            ReferencedEnvelope bbox)
+            ReferencedEnvelope bbox,
+            Double stddevs)
             throws Exception {
         int selectedBand = getRequestedBand(property); // one based band name
         // read the image to be classified
@@ -384,6 +396,7 @@ private List<Rule> getRasterRules(
         RenderedImage image = imageReader.getImage();
 
         RasterSymbolizerBuilder builder = new RasterSymbolizerBuilder();
+        builder.setStandardDeviations(stddevs);
         ColorMap colorMap;
         try {
             if (customClasses.isEmpty()) {
@@ -477,7 +490,8 @@ private List<Rule> getVectorRules(
             int pointSize,
             FeatureTypeInfo obj,
             ColorRamp ramp,
-            ReferencedEnvelope bbox)
+            ReferencedEnvelope bbox,
+            Double stddevs)
             throws IOException, TransformException, FactoryException {
         if (property == null || property.isEmpty()) {
             throw new IllegalArgumentException(
@@ -502,6 +516,25 @@ private List<Rule> getVectorRules(
             query.setHints(getQueryHints(viewParams));
             ftCollection =
                     obj.getFeatureSource(new NullProgressListener(), null).getFeatures(query);
+
+            if (stddevs != null) {
+                NumberRange stdDevRange =
+                        getStandardDeviationsRange(property, ftCollection, stddevs);
+                PropertyIsBetween between =
+                        FF.between(
+                                FF.property(property),
+                                FF.literal(stdDevRange.getMinimum()),
+                                FF.literal(stdDevRange.getMaximum()));
+                if (query.getFilter() == Filter.INCLUDE) {
+                    query.setFilter(between);
+                } else {
+                    query.setFilter(FF.and(query.getFilter(), between));
+                }
+
+                // re-query
+                ftCollection =
+                        obj.getFeatureSource(new NullProgressListener(), null).getFeatures(query);
+            }
         }
 
         List<Rule> rules = null;
@@ -586,6 +619,31 @@ else if (geomT == MultiPolygon.class || geomT == Polygon.class) {
         return rules;
     }
 
+    /**
+     * Returns a range of N standard deviations around the mean for the given attribute and
+     * collection
+     */
+    private NumberRange getStandardDeviationsRange(
+            String property, FeatureCollection features, double numStandardDeviations)
+            throws IOException {
+        final StandardDeviationVisitor standardDeviationVisitor =
+                new StandardDeviationVisitor(FF.property(property));
+        features.accepts(standardDeviationVisitor, null);
+        final double mean = standardDeviationVisitor.getMean();
+        final CalcResult result = standardDeviationVisitor.getResult();
+        if (result.getValue() == null) {
+            throw new RestException(
+                    "The standard deviation visit did not find any value, the dataset is empty or previous filters removed all values",
+                    HttpStatus.BAD_REQUEST);
+        }
+        final double standardDeviation = standardDeviationVisitor.getResult().toDouble();
+
+        return new NumberRange(
+                Double.class,
+                mean - standardDeviation * numStandardDeviations,
+                mean + standardDeviation * numStandardDeviations);
+    }
+
     private ColorRamp getColorRamp(
             String customClasses,
             String colorRamp,

src/extension/sldService/src/main/java/org/geoserver/sldservice/utils/classifier/RasterSymbolizerBuilder.java

@@ -6,10 +6,13 @@
 
 import static java.util.Locale.ENGLISH;
 
+import it.geosolutions.jaiext.JAIExt;
 import it.geosolutions.jaiext.classbreaks.ClassBreaksDescriptor;
 import it.geosolutions.jaiext.classbreaks.ClassBreaksRIF;
 import it.geosolutions.jaiext.classbreaks.Classification;
 import it.geosolutions.jaiext.classbreaks.ClassificationMethod;
+import it.geosolutions.jaiext.stats.Statistics;
+import it.geosolutions.jaiext.stats.Statistics.StatsType;
 import java.awt.*;
 import java.awt.image.DataBuffer;
 import java.awt.image.RenderedImage;
@@ -19,14 +22,17 @@
 import java.util.List;
 import java.util.Optional;
 import javax.media.jai.Histogram;
+import javax.media.jai.JAI;
 import javax.media.jai.ParameterBlockJAI;
+import javax.media.jai.RenderedOp;
 import org.geotools.factory.CommonFactoryFinder;
 import org.geotools.filter.function.RangedClassifier;
 import org.geotools.image.ImageWorker;
 import org.geotools.styling.ColorMap;
 import org.geotools.styling.ColorMapEntry;
 import org.geotools.styling.StyleFactory;
 import org.geotools.util.Converters;
+import org.geotools.util.NumberRange;
 import org.geotools.util.factory.GeoTools;
 import org.opengis.filter.FilterFactory2;
 
@@ -57,6 +63,7 @@ public class RasterSymbolizerBuilder {
             Integer.getInteger("org.geoserver.sldService.maxPixels", 4194304);
 
     private long maxPixels;
+    private Double standardDeviations;
 
     /**
      * Builds the {@link RasterSymbolizerBuilder} with a given pixel reading threshold before
@@ -83,29 +90,26 @@ public RasterSymbolizerBuilder() {
      *     exception should be thrown
      */
     public ColorMap uniqueIntervalClassification(RenderedImage image, Integer maxIntervals) {
-        // compute min and max, for the common case avoid doing an extrema a use a pre-defined range
-        // instead
         int low, high;
         int dataType = image.getSampleModel().getDataType();
         ImageWorker iw = getImageWorker(image);
-        switch (dataType) {
-            case DataBuffer.TYPE_BYTE:
-                low = 0;
-                high = 255;
-                break;
-                // The histogram can be very expensive memory wise as it's backed by a
-                // AtomicDouble[],
-                // check how many
-            case DataBuffer.TYPE_USHORT:
-            case DataBuffer.TYPE_SHORT:
-            case DataBuffer.TYPE_INT:
-                low = (int) iw.getMinimums()[0];
-                high = (int) iw.getMaximums()[0];
-                break;
-            default:
-                throw new IllegalArgumentException(
-                        "Cannot perform unique value classification over rasters of float type, only integer numbers are supported. Try a classification by intervals or quantiles instead");
+
+        // compute min and max, for the common case avoid doing an extrema or use a pre-defined
+        // range instead
+        if (dataType == DataBuffer.TYPE_BYTE && standardDeviations == null) {
+            low = 0;
+            high = 255;
+        } else if (dataType == DataBuffer.TYPE_DOUBLE || dataType == DataBuffer.TYPE_FLOAT) {
+            throw new IllegalArgumentException(
+                    "Cannot perform unique value classification over rasters of float type, only integer numbers are supported. Try a classification by intervals or quantiles instead");
+        } else {
+            final NumberRange range = getOperationRange(iw);
+            low = (int) range.getMinimum();
+            high = (int) range.getMaximum();
         }
+
+        // The histogram can be very expensive memory wise as it's backed by a
+        // AtomicDouble[], check how many are they going to be
         if (high - low > MAX_UNIQUE_VALUES) {
             throw new IllegalArgumentException(
                     "Cannot perform unique value classification over rasters with a potential range of values greater than "
@@ -181,8 +185,9 @@ ImageWorker getImageWorker(RenderedImage image) {
     public ColorMap equalIntervalClassification(
             RenderedImage image, int intervals, boolean open, boolean continuous) {
         ImageWorker iw = getImageWorker(image);
-        double low = iw.getMinimums()[0];
-        double high = iw.getMaximums()[0];
+        final NumberRange range = getOperationRange(iw);
+        double low = (int) range.getMinimum();
+        double high = (int) range.getMaximum();
 
         Number[] breaks = new Number[continuous ? intervals : intervals + 1];
         double step = (high - low) / (continuous ? (intervals - 1) : intervals);
@@ -324,9 +329,10 @@ private Number[] getClassificationBreaks(
         pb.set(iw.getYPeriod(), 6);
         pb.set(noData, 7);
         if (numHistogramBins > 0) {
+            final NumberRange range = getOperationRange(iw);
             Double[][] extrema = new Double[2][1];
-            extrema[0][0] = iw.getMinimums()[0];
-            extrema[1][0] = iw.getMaximums()[0];
+            extrema[0][0] = range.getMinimum();
+            extrema[1][0] = range.getMaximum();
             pb.set(extrema, 2);
             pb.set(true, 8);
             pb.set(numHistogramBins, 9);
@@ -379,4 +385,50 @@ public ColorMap createCustomColorMap(
 
         return getColorMapFromBreaks(breaks, open, continuous);
     }
+
+    public void setStandardDeviations(Double standardDeviations) {
+        this.standardDeviations = standardDeviations;
+    }
+
+    private NumberRange getOperationRange(ImageWorker iw) {
+        if (standardDeviations == null) {
+            double min = iw.getMinimums()[0];
+            double max = iw.getMaximums()[0];
+            return new NumberRange(Double.class, min, max);
+        } else {
+            // Create the parameterBlock
+            ParameterBlock pb = new ParameterBlock();
+            pb.setSource(iw.getRenderedImage(), 0);
+            if (JAIExt.isJAIExtOperation("Stats")) {
+                StatsType[] stats =
+                        new StatsType[] {StatsType.MEAN, StatsType.DEV_STD, StatsType.EXTREMA};
+
+                // Image parameters
+                pb.set(iw.getXPeriod(), 0); // xPeriod
+                pb.set(iw.getYPeriod(), 1); // yPeriod
+                pb.set(iw.getROI(), 2); // ROI
+                pb.set(iw.getNoData(), 3); // NoData
+                pb.set(stats, 6); // statistic operation
+                final RenderedOp statsImage = JAI.create("Stats", pb, iw.getRenderingHints());
+                // Retrieving the statistics
+                Statistics[][] results =
+                        (Statistics[][]) statsImage.getProperty(Statistics.STATS_PROPERTY);
+                double mean = (double) results[0][0].getResult();
+                double stddev = (double) results[0][1].getResult();
+                double[] extrema = (double[]) results[0][2].getResult();
+                double min = extrema[0];
+                double max = extrema[1];
+                // return a range centered in the mean with the desired number of standard
+                // deviations, but make sure it does not exceed the data minimim and maximums
+                return new NumberRange(
+                        Double.class,
+                        Math.max(mean - stddev * standardDeviations, min),
+                        Math.min(mean + stddev * standardDeviations, max));
+            } else {
+                // the op should be unique to jai-ext but best be careful
+                throw new IllegalArgumentException(
+                        "Stats image operation is not backed by JAIExt, please enable JAIExt");
+            }
+        }
+    }
 }