org.jgrasstools.lesto.modules.raster.LasOnRasterMapper Maven / Gradle / Ivy

Go to download
Show more of this group Show more artifacts with this name
Show all versions of jgt-lesto Show documentation
The newest version!
/*
 * 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.lesto.modules.raster;

import static java.lang.Math.round;
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_DRAFT;
import static org.jgrasstools.gears.i18n.GearsMessages.OMSHYDRO_LICENSE;

import java.awt.Point;
import java.awt.image.WritableRaster;
import java.io.File;
import java.util.List;

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

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.Status;
import oms3.annotations.UI;

import org.geotools.coverage.grid.GridCoverage2D;
import org.geotools.coverage.grid.GridGeometry2D;
import org.jgrasstools.gears.io.las.ALasDataManager;
import org.jgrasstools.gears.io.las.core.LasRecord;
import org.jgrasstools.gears.libs.modules.JGTConstants;
import org.jgrasstools.gears.libs.modules.JGTModel;
import org.jgrasstools.gears.libs.monitor.IJGTProgressMonitor;
import org.jgrasstools.gears.utils.RegionMap;
import org.jgrasstools.gears.utils.coverage.CoverageUtilities;
import org.opengis.referencing.crs.CoordinateReferenceSystem;

import com.vividsolutions.jts.geom.Coordinate;
import com.vividsolutions.jts.geom.Polygon;

@Description("Module that creates a raster by mapping maximum elevation las points.")
@Author(name = OMSHYDRO_AUTHORNAMES, contact = OMSHYDRO_AUTHORCONTACTS)
@Keywords("raster, lidar")
@Label(JGTConstants.LESTO + "/raster")
@Name("lasonrastermapper")
@Status(OMSHYDRO_DRAFT)
@License(OMSHYDRO_LICENSE)
@SuppressWarnings("nls")
public class LasOnRasterMapper extends JGTModel {

    @Description("Las file path.")
    @UI(JGTConstants.FILEIN_UI_HINT)
    @In
    public String inLas = null;

    @Description("A dtm raster to use for the area of interest and lower threshold.")
    @UI(JGTConstants.FILEIN_UI_HINT)
    @In
    public String inDtm;

    @Description("New x resolution (if null, the dtm is used).")
    @In
    public Double pXres;

    @Description("New y resolution (if null, the dtm is used).")
    @In
    public Double pYres;

    @Description("If true, the maxima are mapped else the minima.")
    @In
    public boolean doMax = true;

    @Description("The output raster.")
    @UI(JGTConstants.FILEOUT_UI_HINT)
    @In
    public String outRaster = null;

    @Description("The input dtm, resampled on the new resolution.")
    @UI(JGTConstants.FILEOUT_UI_HINT)
    @In
    public String outDtm = null;

    @Execute
    public void process() throws Exception {
        checkNull(inLas, inDtm, outRaster);

        GridCoverage2D inDtmGC = getRaster(inDtm);
        RegionMap regionMap = CoverageUtilities.getRegionParamsFromGridCoverage(inDtmGC);
        double north = regionMap.getNorth();
        double south = regionMap.getSouth();
        double east = regionMap.getEast();
        double west = regionMap.getWest();
        if (pXres == null || pYres == null) {
            pXres = regionMap.getXres();
            pYres = regionMap.getYres();
        }

        CoordinateReferenceSystem crs = null;
        Polygon polygon = CoverageUtilities.getRegionPolygon(inDtmGC);
        crs = inDtmGC.getCoordinateReferenceSystem();
        GridGeometry2D dtmGridGeometry = inDtmGC.getGridGeometry();

        final int newRows = (int) round((north - south) / pYres);
        int newCols = (int) round((east - west) / pXres);

        final GridGeometry2D newGridGeometry2D = CoverageUtilities.gridGeometryFromRegionValues(north, south, east, west,
                newCols, newRows, crs);
        RegionMap newRegionMap = CoverageUtilities.gridGeometry2RegionParamsMap(newGridGeometry2D);
        final WritableRaster newWR = CoverageUtilities.createDoubleWritableRaster(newCols, newRows, null, null,
                JGTConstants.doubleNovalue);

        RandomIter dtmIter = CoverageUtilities.getRandomIterator(inDtmGC);

        try (ALasDataManager lasData = ALasDataManager.getDataManager(new File(inLas), null, 0.0, crs)) {
            lasData.open();
            pm.beginTask("Reading points on region...", IJGTProgressMonitor.UNKNOWN);
            List lasPoints = lasData.getPointsInGeometry(polygon, false);
            pm.done();

            pm.beginTask("Setting raster points...", lasPoints.size());
            final Point gridPoint = new Point();
            final Point dtmPoint = new Point();
            for( LasRecord lasRecord : lasPoints ) {
                double dotZ = lasRecord.z;
                Coordinate coordinate = new Coordinate(lasRecord.x, lasRecord.y, dotZ);
                CoverageUtilities.colRowFromCoordinate(coordinate, newGridGeometry2D, gridPoint);

                double newRasterValue = newWR.getSampleDouble(gridPoint.x, gridPoint.y, 0);

                CoverageUtilities.colRowFromCoordinate(coordinate, dtmGridGeometry, dtmPoint);
                double dtmValue = dtmIter.getSampleDouble(dtmPoint.x, dtmPoint.y, 0);

                if (doMax) {
                    if (JGTConstants.isNovalue(newRasterValue) || newRasterValue < dotZ) {
                        if (!JGTConstants.isNovalue(dtmValue) && dtmValue > dotZ) {
                            dotZ = dtmValue;
                        }
                        newWR.setSample(gridPoint.x, gridPoint.y, 0, dotZ);
                    }
                } else {
                    if (JGTConstants.isNovalue(newRasterValue) || newRasterValue > dotZ) {
                        if (!JGTConstants.isNovalue(dtmValue) && dtmValue > dotZ) {
                            dotZ = dtmValue;
                        }
                        newWR.setSample(gridPoint.x, gridPoint.y, 0, dotZ);
                    }
                }
                pm.worked(1);
            }
            pm.done();
        }

        GridCoverage2D outRasterGC = CoverageUtilities.buildCoverage("outraster", newWR, newRegionMap, crs);
        dumpRaster(outRasterGC, outRaster);

        if (pXres != null && pYres != null && outDtm != null) {
            WritableRaster[] holder = new WritableRaster[1];
            GridCoverage2D outDtmGC = CoverageUtilities.createCoverageFromTemplate(outRasterGC, JGTConstants.doubleNovalue,
                    holder);

            GridGeometry2D outGridGeometry = outDtmGC.getGridGeometry();

            WritableRandomIter outIter = CoverageUtilities.getWritableRandomIterator(holder[0]);

            RegionMap outRegionMap = CoverageUtilities.getRegionParamsFromGridCoverage(outDtmGC);
            int cols = outRegionMap.getCols();
            int rows = outRegionMap.getRows();

            final Point p = new Point();
            for( int c = 0; c < cols; c++ ) {
                for( int r = 0; r < rows; r++ ) {
                    Coordinate coordinate = CoverageUtilities.coordinateFromColRow(c, r, outGridGeometry);
                    CoverageUtilities.colRowFromCoordinate(coordinate, dtmGridGeometry, p);
                    double dtmValue = dtmIter.getSampleDouble(p.x, p.y, 0);
                    outIter.setSample(c, r, 0, dtmValue);
                }
            }
            outIter.done();

            dumpRaster(outDtmGC, outDtm);
        }
        dtmIter.done();

    }

}