org.jgrasstools.gears.modules.r.summary.OmsZonalStatsIM Maven / Gradle / Ivy

/*
 * This file is part of JGrasstools (http://www.jgrasstools.org)
 * (C) HydroloGIS - www.hydrologis.com 
 * 
 * JGrasstools is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see .
 */
package org.jgrasstools.gears.modules.r.summary;

import static org.jgrasstools.gears.i18n.GearsMessages.OMSHYDRO_AUTHORCONTACTS;
import static org.jgrasstools.gears.i18n.GearsMessages.OMSHYDRO_AUTHORNAMES;
import static org.jgrasstools.gears.i18n.GearsMessages.OMSHYDRO_LICENSE;

import java.awt.image.Raster;
import java.io.File;
import java.util.List;
import java.util.concurrent.ConcurrentLinkedQueue;

import javax.media.jai.iterator.RandomIter;
import javax.media.jai.iterator.RandomIterFactory;

import oms3.annotations.Author;
import oms3.annotations.Description;
import oms3.annotations.Execute;
import oms3.annotations.In;
import oms3.annotations.Keywords;
import oms3.annotations.Label;
import oms3.annotations.License;
import oms3.annotations.Name;
import oms3.annotations.Out;
import oms3.annotations.Status;

import org.geotools.coverage.grid.GridCoverage2D;
import org.geotools.coverage.grid.GridGeometry2D;
import org.geotools.data.simple.SimpleFeatureCollection;
import org.geotools.feature.DefaultFeatureCollection;
import org.geotools.feature.simple.SimpleFeatureBuilder;
import org.geotools.geometry.jts.ReferencedEnvelope;
import org.jgrasstools.gears.libs.modules.JGTConstants;
import org.jgrasstools.gears.libs.modules.JGTModelIM;
import org.jgrasstools.gears.utils.features.FeatureUtilities;
import org.jgrasstools.gears.utils.math.NumericsUtilities;
import org.opengis.feature.simple.SimpleFeature;

import com.vividsolutions.jts.geom.Envelope;
import com.vividsolutions.jts.geom.Geometry;

@Description("Calculate zonal stats on image mosaic datasets.")
@Author(name = OMSHYDRO_AUTHORNAMES, contact = OMSHYDRO_AUTHORCONTACTS)
@Keywords("zonalstats, image mosaic")
@Label(JGTConstants.RASTERPROCESSING)
@Name("zonalstats_im")
@Status(Status.EXPERIMENTAL)
@License(OMSHYDRO_LICENSE)
public class OmsZonalStatsIM extends JGTModelIM {

    @Description("The image mosaic map to process..")
    @In
    public String inRaster = null;

    @Description("The polygons map on which to do the stats.")
    @In
    public SimpleFeatureCollection inVector = null;

    @Description("The size of the bins to use for pretiling the region.")
    @In
    public double pBinSize = 10000;

    @Description("Percentage of minimum active cells to have a valid stat.")
    @In
    public double pPercentageThres = 20.0;

    @Description("Total meanvalue (also produced by this module) for the calculation of the mean absolute deviation.")
    @In
    public Double pTotalMean = null;

    @Description("The input polygons with the added stats values.")
    @Out
    public SimpleFeatureCollection outVector;

    /**
     * The array holding:
     * 
     *  totalMean
     *  userTotalMean
     *  totalActiveCells
     * 
     * if {@link #pTotalMean} is != null.
     */
    double[] tm_usertm_tactivecells = new double[3];

    @Execute
    public void process() throws Exception {
        checkNull(inVector);

        addSource(new File(inRaster));

        boolean hasUserTotalMean = false;
        if (pTotalMean != null) {
            hasUserTotalMean = true;
            tm_usertm_tactivecells[1] = pTotalMean;
        }

        ReferencedEnvelope bounds = inVector.getBounds();
        double[] xBins = NumericsUtilities.range2Bins(bounds.getMinX(), bounds.getMaxX(), pBinSize, false);
        double[] yBins = NumericsUtilities.range2Bins(bounds.getMinY(), bounds.getMaxY(), pBinSize, false);

        SimpleFeatureBuilder featureBuilder = OmsZonalStats.createFeatureBuilder(bounds.getCoordinateReferenceSystem(),
                hasUserTotalMean);

        outVector = new DefaultFeatureCollection();
        List geometriesList = FeatureUtilities.featureCollectionToGeometriesList(inVector, true, null);

        ConcurrentLinkedQueue allGeometriesQueue = new ConcurrentLinkedQueue();
        allGeometriesQueue.addAll(geometriesList);
        ConcurrentLinkedQueue keepGeometriesQueue;
        ConcurrentLinkedQueue removeGeometriesQueue;
        pm.beginTask("Processing polygons...", xBins.length - 1);
        for( int x = 0; x < xBins.length - 1; x++ ) {
            for( int y = 0; y < yBins.length - 1; y++ ) {
                Envelope envelope = new Envelope(xBins[x], xBins[x + 1], yBins[y], yBins[y + 1]);
                Envelope readEnvelope = null;
                keepGeometriesQueue = new ConcurrentLinkedQueue();
                removeGeometriesQueue = new ConcurrentLinkedQueue();
                for( Geometry geometry : allGeometriesQueue ) {
                    Envelope geometryenvelope = geometry.getEnvelopeInternal();
                    if (geometryenvelope.intersects(envelope)) {
                        removeGeometriesQueue.add(geometry);
                        if (readEnvelope == null) {
                            readEnvelope = new Envelope(geometryenvelope);
                        } else {
                            readEnvelope.expandToInclude(geometryenvelope);
                        }
                    } else {
                        keepGeometriesQueue.add(geometry);
                    }
                }
                allGeometriesQueue = keepGeometriesQueue;
                if (removeGeometriesQueue.size() == 0) {
                    continue;
                }

                // pm.message("" + readEnvelope);
                GridCoverage2D readGC = getGridCoverage(0, readEnvelope);
                GridGeometry2D gridGeometry = readGC.getGridGeometry();
                Raster readRaster = readGC.getRenderedImage().getData();
                RandomIter readIter = RandomIterFactory.create(readRaster, null);
                for( Geometry geometry : removeGeometriesQueue ) {
                    double[] polygonStats = OmsZonalStats.polygonStats(geometry, gridGeometry, readIter, hasUserTotalMean,
                            tm_usertm_tactivecells, pPercentageThres, pm);
                    if (polygonStats == null) {
                        continue;
                    }

                    Object[] values;
                    if (pTotalMean == null) {
                        values = new Object[]{geometry, //
                                polygonStats[0], //
                                polygonStats[1], //
                                polygonStats[2], //
                                polygonStats[3], //
                                polygonStats[4], //
                                (int) polygonStats[5], //
                                (int) polygonStats[6] //
                        };
                    } else {
                        values = new Object[]{geometry, //
                                polygonStats[0], //
                                polygonStats[1], //
                                polygonStats[2], //
                                polygonStats[3], //
                                polygonStats[4], //
                                polygonStats[5], //
                                (int) polygonStats[6], //
                                (int) polygonStats[7] //
                        };
                    }

                    featureBuilder.addAll(values);
                    SimpleFeature feature = featureBuilder.buildFeature(null);
                    ((DefaultFeatureCollection) outVector).add(feature);
                }
            }
            pm.worked(1);
        }
        pm.done();

        if (!hasUserTotalMean) {
            tm_usertm_tactivecells[0] = tm_usertm_tactivecells[0] / tm_usertm_tactivecells[2];
            pm.message("Total mean: " + tm_usertm_tactivecells[0]);
        }

        dispose();
    }

    protected void processCell( int readCol, int readRow, int writeCol, int writeRow, int readCols, int readRows, int writeCols,
            int writeRows ) {
        // not used in this case
    }
}