Java Code Examples for org.geotools.coverage.grid.GridCoverage2D#getGridGeometry()

/**
 * Get the array of region parameters covered by the {@link GridCoverage2D coverage}. 
 * 
 * @param gridCoverage the coverage.
 * @return the array of region parameters as [n, s, w, e, xres, yres, cols, rows]
 */
public static double[] getRegionArrayFromGridCoverage( GridCoverage2D gridCoverage ) {
    Envelope envelope = gridCoverage.getEnvelope();
    DirectPosition lowerCorner = envelope.getLowerCorner();
    double[] westSouth = lowerCorner.getCoordinate();
    DirectPosition upperCorner = envelope.getUpperCorner();
    double[] eastNorth = upperCorner.getCoordinate();

    GridGeometry2D gridGeometry = gridCoverage.getGridGeometry();
    GridEnvelope2D gridRange = gridGeometry.getGridRange2D();
    int height = gridRange.height;
    int width = gridRange.width;

    AffineTransform gridToCRS = (AffineTransform) gridGeometry.getGridToCRS();
    double xRes = XAffineTransform.getScaleX0(gridToCRS);
    double yRes = XAffineTransform.getScaleY0(gridToCRS);

    double[] params = new double[]{eastNorth[1], westSouth[1], westSouth[0], eastNorth[0], xRes, yRes, width, height};

    return params;
}
 

/**
 * Get the parameters of the region covered by the {@link GridCoverage2D coverage}. 
 * 
 * @param gridCoverage the coverage.
 * @return the {@link HashMap map} of parameters. ( {@link #NORTH} and the 
 *          other static vars can be used to retrieve them.
 */
public static RegionMap getRegionParamsFromGridCoverage( GridCoverage2D gridCoverage ) {
    RegionMap envelopeParams = new RegionMap();

    Envelope envelope = gridCoverage.getEnvelope();

    DirectPosition lowerCorner = envelope.getLowerCorner();
    double[] westSouth = lowerCorner.getCoordinate();
    DirectPosition upperCorner = envelope.getUpperCorner();
    double[] eastNorth = upperCorner.getCoordinate();

    GridGeometry2D gridGeometry = gridCoverage.getGridGeometry();
    GridEnvelope2D gridRange = gridGeometry.getGridRange2D();
    int height = gridRange.height;
    int width = gridRange.width;

    AffineTransform gridToCRS = (AffineTransform) gridGeometry.getGridToCRS();
    double xRes = XAffineTransform.getScaleX0(gridToCRS);
    double yRes = XAffineTransform.getScaleY0(gridToCRS);

    envelopeParams.put(NORTH, eastNorth[1]);
    envelopeParams.put(SOUTH, westSouth[1]);
    envelopeParams.put(WEST, westSouth[0]);
    envelopeParams.put(EAST, eastNorth[0]);
    envelopeParams.put(XRES, xRes);
    envelopeParams.put(YRES, yRes);
    envelopeParams.put(ROWS, (double) height);
    envelopeParams.put(COLS, (double) width);

    return envelopeParams;
}
 

/**
 * Get the array of rows and cols. 
 * 
 * @param gridCoverage the coverage.
 * @return the array as [cols, rows]
 */
public static int[] getRegionColsRows( GridCoverage2D gridCoverage ) {
    GridGeometry2D gridGeometry = gridCoverage.getGridGeometry();
    GridEnvelope2D gridRange = gridGeometry.getGridRange2D();
    int height = gridRange.height;
    int width = gridRange.width;
    int[] params = new int[]{width, height};
    return params;
}
 

/**
 * Get the cols and rows ranges to use to loop the original gridcoverage.
 * 
 * @param gridCoverage the coverage.
 * @param subregion the sub region of the coverage to get the cols and rows to loop on.
 * @return the array of looping values in the form [minCol, maxCol, minRow, maxRow].
 * @throws Exception
 */
public static int[] getLoopColsRowsForSubregion( GridCoverage2D gridCoverage, Envelope2D subregion ) throws Exception {
    GridGeometry2D gridGeometry = gridCoverage.getGridGeometry();
    GridEnvelope2D subRegionGrid = gridGeometry.worldToGrid(subregion);
    int minCol = subRegionGrid.x;
    int maxCol = subRegionGrid.x + subRegionGrid.width;
    int minRow = subRegionGrid.y;
    int maxRow = subRegionGrid.y + subRegionGrid.height;
    return new int[]{minCol, maxCol, minRow, maxRow};
}
 

/**
 * Calculates the profile of a raster map between given {@link Coordinate coordinates}.
 * 
 * @param coverage the coverage from which to extract the profile.
 * @param coordinates the coordinates to use to trace the profile.
 * @return the list of {@link ProfilePoint}s.
 * @throws Exception
 */
public static List<ProfilePoint> doProfile( GridCoverage2D coverage, Coordinate... coordinates ) throws Exception {
    GridGeometry2D gridGeometry = coverage.getGridGeometry();
    RenderedImage renderedImage = coverage.getRenderedImage();
    RandomIter iter = RandomIterFactory.create(renderedImage, null);

    List<ProfilePoint> profilePointsList = doProfile(iter, gridGeometry, coordinates);

    iter.done();
    return profilePointsList;
}
 

public static String getRegionPrint( GridCoverage2D coverage ) {
    org.opengis.geometry.Envelope envelope = coverage.getEnvelope();

    DirectPosition lowerCorner = envelope.getLowerCorner();
    double[] westSouth = lowerCorner.getCoordinate();
    DirectPosition upperCorner = envelope.getUpperCorner();
    double[] eastNorth = upperCorner.getCoordinate();

    GridGeometry2D gridGeometry = coverage.getGridGeometry();
    GridEnvelope2D gridRange = gridGeometry.getGridRange2D();

    String numberSpace = "                ";
    String numberRest = "                  ";

    StringBuilder sb = new StringBuilder();
    sb.append("        +----------------------   ").append(String.format("%17.8f", eastNorth[1]))
            .append("   ----------------------+").append("\n");
    sb.append("        +                         ").append(String.format("%17.8f", gridRange.getMinY()))
            .append("                         +").append("\n");
    sb.append("        +                         ").append(numberSpace).append("                          +").append("\n");
    sb.append("        +                         ").append(numberSpace).append("                          +").append("\n");
    sb.append(String.format("%17.8f", westSouth[0])).append(numberRest).append(numberSpace).append(numberRest)
            .append(String.format("%17.8f", eastNorth[0])).append("\n");
    sb.append(String.format("%17.8f", gridRange.getMinX())).append(numberRest).append(numberSpace).append(numberRest)
            .append(String.format("%17.8f", gridRange.getMaxX())).append("\n");
    sb.append("        +                         ").append(numberSpace).append("                          +").append("\n");
    sb.append("        +                         ").append(numberSpace).append("                          +").append("\n");
    sb.append("        +                         ").append(String.format("%17.8f", gridRange.getMaxY()))
            .append("                         +").append("\n");
    sb.append("        +----------------------   ").append(String.format("%17.8f", westSouth[1]))
            .append("   ----------------------+").append("\n");

    return sb.toString();
}
 

@Execute
public void process() throws Exception {

    OmsTimeSeriesIteratorReader reader = new OmsTimeSeriesIteratorReader();
    reader.file = inData;
    reader.idfield = fStationsid;
    // reader.tStart = "2000-01-01 00:00";
    // reader.tTimestep = 60;
    // reader.tEnd = "2000-01-01 00:00";
    reader.fileNovalue = "-9999";
    reader.initProcess();

    OmsKriging kriging = new OmsKriging();
    kriging.inStations = getVector(inStations);
    kriging.fStationsid = fStationsid;
    kriging.fStationsZ = fStationsZ;
    kriging.inInterpolate = getVector(inInterpolate);
    kriging.fInterpolateid = fInterpolateid;
    kriging.fPointZ = fPointZ;
    kriging.pMode = pMode;
    if (pIntegralscale!=null && pIntegralscale.trim().length()>0) {
        String[] split = pIntegralscale.split(",");
        double[] integralScaleDouble = new double[split.length];
        for( int i = 0; i < split.length; i++ ) {
            try {
                integralScaleDouble[i] = Double.parseDouble(split[i]);
            } catch (Exception e) {
                throw new ModelsIllegalargumentException("Problems with integral scale: " + pIntegralscale, this, pm);
            }
        }
        kriging.pIntegralscale = integralScaleDouble;
    }
    kriging.pVariance = pVariance;
    kriging.doLogarithmic = doLogarithmic;
    GridCoverage2D interpolationGrid = getRaster(inInterpolationGrid);
    if (interpolationGrid != null) {
        kriging.inInterpolationGrid = interpolationGrid.getGridGeometry();
    }
    kriging.defaultVariogramMode = defaultVariogramMode;
    kriging.pSemivariogramType = pSemivariogramType;
    kriging.doIncludezero = doIncludezero;
    kriging.pA = pA;
    kriging.pS = pS;
    kriging.pNug = pNug;
    kriging.pm = pm;
    kriging.doProcess = doProcess;
    kriging.doReset = doReset;

    OmsTimeSeriesIteratorWriter writer = null;

    pm.beginTask("Processing...", IHMProgressMonitor.UNKNOWN);
    while( reader.doProcess ) {
        reader.nextRecord();

        if (writer == null) {
            writer = new OmsTimeSeriesIteratorWriter();
            writer.file = outData;
            writer.tStart = reader.tStart;
            writer.tTimestep = reader.tTimestep;
        }
        pm.message("timestep: " + reader.tCurrent);

        HashMap<Integer, double[]> id2ValueMap = reader.outData;
        kriging.inData = id2ValueMap;
        kriging.process();

        // write csv data
        writer.inData = kriging.outData;
        writer.writeNextLine();

        // write raster data
        dumpRaster(kriging.outGrid, outGrid);
    }
    pm.done();

    reader.close();
    writer.close();
}
 

private double[] resFromCoverage( GridCoverage2D dem ) {
    GridGeometry2D gridGeometry = dem.getGridGeometry();
    AffineTransform gridToCRS = (AffineTransform) gridGeometry.getGridToCRS();
    double[] res = new double[]{XAffineTransform.getScaleX0(gridToCRS), XAffineTransform.getScaleY0(gridToCRS)};
    return res;
}
 

/**
 * Extracts sections on the dtm from a riverline at regular intervals.
 * 
 * @param riverLine the river line to consider for the cross sections extraction.
 *          <b>The river line has to be oriented from upstream to downstream.</b>
 * @param elevation the elevation {@link GridCoverage2D}.
 * @param sectionsInterval the distance to use between extracted sections. 
 * @param sectionsWidth the width of the extracted sections.
 * @param bridgePoints the list of bridge {@link Point}s. 
 * @param bridgeWidthAttribute the name of the attribute in the bridges feature that defines the width of the bridge.
 * @param bridgeBuffer a buffer to use for the bridge inside which the bridge is considered to be on the river.
 * @param monitor the progress monitor.
 * @throws Exception
 */
public RiverSectionsFromDtmExtractor( //
        LineString riverLine, //
        GridCoverage2D elevation, //
        double sectionsInterval, //
        double sectionsWidth, //
        List<FeatureMate> bridgePoints, //
        String bridgeWidthAttribute, //
        double bridgeBuffer, //
        IHMProgressMonitor monitor //
) throws Exception {
    crs = elevation.getCoordinateReferenceSystem();

    RandomIter elevIter = CoverageUtilities.getRandomIterator(elevation);
    GridGeometry2D gridGeometry = elevation.getGridGeometry();
    RegionMap regionMap = CoverageUtilities.getRegionParamsFromGridCoverage(elevation);
    Envelope envelope = regionMap.toEnvelope();

    riverPointsList = new ArrayList<RiverPoint>();
    LengthIndexedLine indexedLine = new LengthIndexedLine(riverLine);

    double length = riverLine.getLength();
    int totalWork = (int) (length / sectionsInterval);
    monitor.beginTask("Extracting sections...", totalWork);
    double runningLength = 0;
    while( runningLength <= length ) {
        // important to extract from left to right

        // TODO extract with point before and after in order to have more regular sections
        Coordinate leftPoint = indexedLine.extractPoint(runningLength, sectionsWidth);
        Coordinate rightPoint = indexedLine.extractPoint(runningLength, -sectionsWidth);
        if (envelope.intersects(leftPoint) && envelope.intersects(rightPoint)) {
            RiverPoint netPoint = getNetworkPoint(riverLine, elevIter, gridGeometry, runningLength, null, leftPoint,
                    rightPoint);
            if (netPoint != null)
                riverPointsList.add(netPoint);
        }
        runningLength = runningLength + sectionsInterval;
        monitor.worked(1);
    }
    monitor.done();

    process(riverLine, sectionsWidth, bridgePoints, bridgeWidthAttribute, bridgeBuffer, elevIter, gridGeometry, envelope,
            indexedLine);
}
 

/**
     * Extracts sections on the dtm on predefined points with a given id and KS.
     * 
     * @param riverLine the river line to consider for the cross sections extraction.
     *          <b>The river line has to be oriented from upstream to downstream.</b>
     * @param riverPointCoordinates the points coordinates.
     * @param riverPointIds the points ids.
     * @param riverPointKs the points KS.
     * @param elevation the elevation {@link GridCoverage2D}.
     * @param sectionsInterval the distance to use between extracted sections. 
     * @param sectionsWidth the width of the extracted sections.
     * @param bridgePoints the list of bridge {@link Point}s. 
     * @param bridgeWidthAttribute the name of the attribute in the bridges feature that defines the width of the bridge.
     * @param bridgeBuffer a buffer to use for the bridge inside which the bridge is considered to be on the river.
     * @param monitor the progress monitor.
     * @throws Exception
     */
    public RiverSectionsFromDtmExtractor( //
            LineString riverLine, //
            Coordinate[] riverPointCoordinates, //
            int[] riverPointIds, //
            double[] riverPointKs, GridCoverage2D elevation, //
            double sectionsInterval, //
            double sectionsWidth, //
            List<FeatureMate> bridgePoints, //
            String bridgeWidthAttribute, //
            double bridgeBuffer, //
            IHMProgressMonitor monitor //
    ) throws Exception {
        crs = elevation.getCoordinateReferenceSystem();

        RandomIter elevIter = CoverageUtilities.getRandomIterator(elevation);
        GridGeometry2D gridGeometry = elevation.getGridGeometry();
        RegionMap regionMap = CoverageUtilities.getRegionParamsFromGridCoverage(elevation);
        Envelope envelope = regionMap.toEnvelope();

        riverPointsList = new ArrayList<RiverPoint>();
        LengthIndexedLine indexedLine = new LengthIndexedLine(riverLine);

        monitor.beginTask("Extracting sections in supplied net points...", riverPointCoordinates.length);
        for( int i = 0; i < riverPointCoordinates.length; i++ ) {
            double pointIndex = indexedLine.indexOf(riverPointCoordinates[i]);
            // important to extract from left to right
            Coordinate leftPoint = indexedLine.extractPoint(pointIndex, sectionsWidth);
            Coordinate rightPoint = indexedLine.extractPoint(pointIndex, -sectionsWidth);
            if (envelope.intersects(leftPoint) && envelope.intersects(rightPoint)) {
                // TODO add ks of net point for section

                RiverPoint netPoint = getNetworkPoint(riverLine, elevIter, gridGeometry, pointIndex, riverPointKs[i], leftPoint,
                        rightPoint);
                netPoint.setSectionId(riverPointIds[i]);
//                netPoint.setSectionGaukler(riverPointKs[i]);
                
                if (netPoint != null)
                    riverPointsList.add(netPoint);
            }
            monitor.worked(1);
        }
        monitor.done();

        process(riverLine, sectionsWidth, bridgePoints, bridgeWidthAttribute, bridgeBuffer, elevIter, gridGeometry, envelope,
                indexedLine);
    }
 

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

    SimpleFeatureCollection buildingsFC = null;
    GridCoverage2D inDtmGC = null;
    try (ALasDataManager lasHandler = ALasDataManager.getDataManager(new File(inLas), inDtmGC, 0, null)) {
        lasHandler.open();

        ReferencedEnvelope3D e = lasHandler.getEnvelope3D();

        if (pRasterResolution != null) {
            OmsRasterReader reader = new OmsRasterReader();
            reader.file = inDtm;
            reader.pNorth = e.getMaxY();
            reader.pSouth = e.getMinY();
            reader.pWest = e.getMinX();
            reader.pEast = e.getMaxX();
            reader.pXres = pRasterResolution;
            reader.pYres = pRasterResolution;
            reader.process();
            inDtmGC = reader.outRaster;
        } else {
            inDtmGC = getRaster(inDtm);
            RegionMap regionMap = CoverageUtilities.getRegionParamsFromGridCoverage(inDtmGC);
            pRasterResolution = regionMap.getXres();
        }

        Envelope2D envelope2d = inDtmGC.getEnvelope2D();
        Polygon regionPolygon = FeatureUtilities.envelopeToPolygon(envelope2d);

        WritableRaster[] buildingsHolder = new WritableRaster[1];
        GridCoverage2D newBuildingsRaster = CoverageUtilities.createCoverageFromTemplate(inDtmGC, 1.0, buildingsHolder);
        GridGeometry2D gridGeometry = newBuildingsRaster.getGridGeometry();

        java.awt.Point p = new java.awt.Point();
        IHMProgressMonitor pm = new LogProgressMonitor();
        pm.beginTask("Reading points...", IHMProgressMonitor.UNKNOWN);
        List<LasRecord> pointsInGeometry = lasHandler.getPointsInGeometry(regionPolygon, true);
        pm.done();
        pm.beginTask("Buildings filtering...", (int) pointsInGeometry.size());
        for( LasRecord dot : pointsInGeometry ) {
            Coordinate c = new Coordinate(dot.x, dot.y);
            if (!envelope2d.contains(dot.x, dot.y)) {
                continue;
            }

            double dtmValue = CoverageUtilities.getValue(inDtmGC, c);
            double height = abs(dot.z - dtmValue);
            CoverageUtilities.colRowFromCoordinate(c, gridGeometry, p);
            if (height < pDtmThres) {
                buildingsHolder[0].setSample(p.x, p.y, 0, doubleNovalue);
            }
            pm.worked(1);
        }
        pm.done();

        OmsMorpher morpher = new OmsMorpher();
        morpher.pm = pm;
        morpher.inMap = newBuildingsRaster;
        morpher.pMode = Variables.OPEN;
        morpher.process();
        newBuildingsRaster = morpher.outMap;

        OmsVectorizer buildingsVectorizer = new OmsVectorizer();
        buildingsVectorizer.pm = pm;
        buildingsVectorizer.inRaster = newBuildingsRaster;
        buildingsVectorizer.doRemoveHoles = true;
        double minAreaInCells = pMinArea / pRasterResolution / pRasterResolution;
        buildingsVectorizer.pThres = minAreaInCells;
        buildingsVectorizer.fDefault = "rast";
        buildingsVectorizer.process();

        buildingsFC = buildingsVectorizer.outVector;
        dumpVector(buildingsFC, outBuildings);
        if (doSmoothing && outCleanBuildings != null) {
            // smooth
            SimpleFeatureCollection smoothedFC = smoothBuildings(pDensifyResolution, buildingsFC);

            // remove trees
            DefaultFeatureCollection removedAndSmoothedFC = removeNonBuildings(lasHandler, smoothedFC, inDtmGC,
                    pBuildingsBuffer);

            dumpVector(removedAndSmoothedFC, outCleanBuildings);
        }
    }

}
 

public void process() throws Exception {
    checkNull(inMaxima, inDsm, outRaster);

    GridCoverage2D inDsmGC = getRaster(inDsm);
    RegionMap regionMap = CoverageUtilities.getRegionParamsFromGridCoverage(inDsmGC);
    int cols = regionMap.getCols();
    int rows = regionMap.getRows();
    double xRes = regionMap.getXres();
    double yRes = regionMap.getYres();

    GridGeometry2D gridGeometry = inDsmGC.getGridGeometry();

    RandomIter dsmIter = CoverageUtilities.getRandomIterator(inDsmGC);
    RandomIter dtmIter;
    if (inDtm != null) {
        GridCoverage2D inDtmGC = getRaster(inDtm);
        dtmIter = CoverageUtilities.getRandomIterator(inDtmGC);
    } else {
        WritableRaster dtmWR = CoverageUtilities.createWritableRaster(cols, rows, null, null, 0.0);
        dtmIter = RandomIterFactory.createWritable(dtmWR, null);
    }

    WritableRaster outWR = CoverageUtilities.createWritableRaster(cols, rows, null, null, HMConstants.doubleNovalue);
    outIter = RandomIterFactory.createWritable(outWR, null);

    SimpleFeatureCollection inMaximaFC = getVector(inMaxima);
    List<SimpleFeature> maximaList = FeatureUtilities.featureCollectionToList(inMaximaFC);

    pm.beginTask("Processing region growing...", maximaList.size());
    int index = 0;
    for( SimpleFeature maximaFeature : maximaList ) {
        Coordinate coordinate = ((Geometry) maximaFeature.getDefaultGeometry()).getCoordinate();
        int[] colRow = CoverageUtilities.colRowFromCoordinate(coordinate, gridGeometry, null);

        GridNode startDsmNode = new GridNode(dsmIter, cols, rows, xRes, yRes, colRow[0], colRow[1]);
        GridNode startDtmNode = new GridNode(dtmIter, cols, rows, xRes, yRes, colRow[0], colRow[1]);
        growRegion(startDsmNode, startDtmNode, index, startDsmNode, startDtmNode);

        index++;
        pm.worked(1);
    }
    pm.done();

    dtmIter.done();
    outIter.done();

    GridCoverage2D outRasterGC = CoverageUtilities.buildCoverage("filtered", outWR, regionMap,
            inDsmGC.getCoordinateReferenceSystem());
    dumpRaster(outRasterGC, outRaster);
}