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

/**
 * Create a new raster using a given raster as template.
 * 
 * @param coverage
 * @param makeNew
 */
public Raster( GridCoverage2D raster, boolean makeNew ) {
    this.makeNew = makeNew;

    crs = raster.getCoordinateReferenceSystem();
    regionMap = getRegionParamsFromGridCoverage(raster);
    cols = regionMap.getCols();
    rows = regionMap.getRows();
    west = regionMap.getWest();
    south = regionMap.getSouth();
    east = regionMap.getEast();
    north = regionMap.getNorth();
    xRes = (east - west) / cols;
    yRes = (north - south) / rows;

    if (makeNew) {
        newWR = createWritableRaster(cols, rows, null, null, HMConstants.doubleNovalue);
        iter = RandomIterFactory.createWritable(newWR, null);
    } else {
        iter = RandomIterFactory.create(raster.getRenderedImage(), null);
    }
}
 

/**
 * Constructor.
 * 
 * @param lasFile the las folder index file.
 * @param inDem a dem to normalize the elevation. If <code>null</code>, the original las elevation is used.
 * @param elevThreshold a threshold to use for the elevation normalization.
 * @param inCrs the data {@link org.opengis.referencing.crs.CoordinateReferenceSystem}. if null, the one of the dem is read, if available.
 */
LasFileDataManager( File lasFile, GridCoverage2D inDem, double elevThreshold, CoordinateReferenceSystem inCrs ) {
    this.lasFile = lasFile;
    this.inDem = inDem;
    this.elevThreshold = elevThreshold;

    fileNamesList.add(lasFile.getName());

    try {
        // prj file rules if available
        inCrs = CrsUtilities.readProjectionFile(lasFile.getAbsolutePath(), "las");
    } catch (Exception e) {
        // ignore and try to read
    }
    if (inCrs != null) {
        crs = inCrs;
    } else if (inDem != null) {
        crs = inDem.getCoordinateReferenceSystem();
    } else {
        throw new IllegalArgumentException("The Crs can't be null.");
    }
}
 

public void testTiff() throws Exception {
    File file = getFile("formats/tiff/test.tif");
    GridCoverage2D raster = OmsRasterReader.readRaster(file.getAbsolutePath());
    RegionMap regionMap = CoverageUtilities.getRegionParamsFromGridCoverage(raster);
    Double novalue = CoverageUtilities.getNovalue(raster);
    assertNull(novalue);

    assertEquals(30, regionMap.getCols());
    assertEquals(26, regionMap.getRows());
    assertEquals(7874.0, regionMap.getXres());
    assertEquals(8120.769230769229580, regionMap.getYres(), DELTA);
    assertEquals(688054.25, regionMap.getWest());
    assertEquals(5472037.36448598, regionMap.getSouth(), DELTA);

    CoordinateReferenceSystem crs = raster.getCoordinateReferenceSystem();
    String epsg = CrsUtilities.getCodeFromCrs(crs);
    assertEquals("EPSG:26921", epsg);

}
 

/**
 * Creates a new {@link GridCoverage2D} using an existing as template.
 * 
 * @param template the template to use.
 * @param value the value to set the new raster to, if not <code>null</code>.
 * @param writableRasterHolder an array of length 1 to place the writable raster in, that 
 *                  was can be used to populate the coverage. If <code>null</code>, it is ignored.
 * @return the new coverage.
 */
public static GridCoverage2D createCoverageFromTemplate( GridCoverage2D template, Double value,
        WritableRaster[] writableRasterHolder ) {
    RegionMap regionMap = getRegionParamsFromGridCoverage(template);

    double west = regionMap.getWest();
    double south = regionMap.getSouth();
    double east = regionMap.getEast();
    double north = regionMap.getNorth();
    int cols = regionMap.getCols();
    int rows = regionMap.getRows();
    ComponentSampleModel sampleModel = new ComponentSampleModel(DataBuffer.TYPE_DOUBLE, cols, rows, 1, cols, new int[]{0});

    WritableRaster raster = RasterFactory.createWritableRaster(sampleModel, null);
    if (value != null) {
        // autobox only once
        double v = value;
        for( int y = 0; y < rows; y++ ) {
            for( int x = 0; x < cols; x++ ) {
                raster.setSample(x, y, 0, v);
            }
        }
    }
    if (writableRasterHolder != null) {
        writableRasterHolder[0] = raster;
    }

    Envelope2D writeEnvelope = new Envelope2D(template.getCoordinateReferenceSystem(), west, south, east - west,
            north - south);
    GridCoverageFactory factory = CoverageFactoryFinder.getGridCoverageFactory(null);
    GridCoverage2D coverage2D = factory.create("newraster", raster, writeEnvelope);
    return coverage2D;
}
 

/**
 * Create a subcoverage given a template coverage and an envelope.
 * 
 * @param template the template coverage used for the resolution.
 * @param subregion the envelope to extract to the new coverage. This should
 *                  be snapped on the resolution of the coverage, in order to avoid 
 *                  shifts.
 * @param value the value to set the new raster to, if not <code>null</code>. 
 * @param writableRasterHolder an array of length 1 to place the writable raster in, that 
 *                  was can be used to populate the coverage. If <code>null</code>, it is ignored.
 * @return the new coverage.
 */
public static GridCoverage2D createSubCoverageFromTemplate( GridCoverage2D template, Envelope2D subregion, Double value,
        WritableRaster[] writableRasterHolder ) {
    RegionMap regionMap = getRegionParamsFromGridCoverage(template);
    double xRes = regionMap.getXres();
    double yRes = regionMap.getYres();

    double west = subregion.getMinX();
    double south = subregion.getMinY();
    double east = subregion.getMaxX();
    double north = subregion.getMaxY();

    int cols = (int) ((east - west) / xRes);
    int rows = (int) ((north - south) / yRes);
    ComponentSampleModel sampleModel = new ComponentSampleModel(DataBuffer.TYPE_DOUBLE, cols, rows, 1, cols, new int[]{0});

    WritableRaster writableRaster = RasterFactory.createWritableRaster(sampleModel, null);
    if (value != null) {
        // autobox only once
        double v = value;
        for( int y = 0; y < rows; y++ ) {
            for( int x = 0; x < cols; x++ ) {
                writableRaster.setSample(x, y, 0, v);
            }
        }
    }
    if (writableRasterHolder != null)
        writableRasterHolder[0] = writableRaster;
    Envelope2D writeEnvelope = new Envelope2D(template.getCoordinateReferenceSystem(), west, south, east - west,
            north - south);
    GridCoverageFactory factory = CoverageFactoryFinder.getGridCoverageFactory(null);
    GridCoverage2D coverage2D = factory.create("newraster", writableRaster, writeEnvelope);
    return coverage2D;
}
 

/**
 * Constructor.
 * 
 * @param lasFolderIndexFile the las folder index file.
 * @param inDem a dem to normalize the elevation. If <code>null</code>, the original las elevation is used.
 * @param elevThreshold a threshold to use for the elevation normalization.
 * @param inCrs the data {@link org.opengis.referencing.crs.CoordinateReferenceSystem}. if null, the one of the dem is read, if available.
 */
LasFolderIndexDataManager( File lasFolderIndexFile, GridCoverage2D inDem, double elevThreshold,
        CoordinateReferenceSystem inCrs ) {
    this.lasFolderIndexFile = lasFolderIndexFile;
    this.inDem = inDem;
    this.elevThreshold = elevThreshold;
    lasFolder = lasFolderIndexFile.getParentFile();

    try {
        // prj file rules if available
        inCrs = CrsUtilities.readProjectionFile(lasFolderIndexFile.getAbsolutePath(), "lasfolder");
    } catch (Exception e) {
        // ignore and try to read
    }
    if (inCrs != null) {
        crs = inCrs;
    } else if (inDem != null) {
        crs = inDem.getCoordinateReferenceSystem();
    } else {
        throw new IllegalArgumentException("The Crs can't be null.");
    }
    // indexMap = new LinkedHashMap<String, Pair>() {
    // @Override
    // protected boolean removeEldestEntry(Entry<String, Pair> eldest) {
    // return size() > READERCACHE;
    // }
    // };
    fileName2LasReaderMap = new WeakValueHashMap<String, Pair>();
    fileName2IndexMap = new WeakValueHashMap<String, STRtreeJGT>();
    fileName4LasReaderMapSupport = new ArrayList<String>();
}
 

/**
 * 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, outRaster);

    GridCoverage2D inDtmGC = getRaster(inDtm);
    Polygon polygon = CoverageUtilities.getRegionPolygon(inDtmGC);
    CoordinateReferenceSystem crs = inDtmGC.getCoordinateReferenceSystem();

    List<Coordinate> lasCoordinates = new ArrayList<Coordinate>();
    pm.beginTask("Preparing triangulation...", -1);
    try (ALasDataManager lasData = ALasDataManager.getDataManager(new File(inLas), null, 0.0, crs)) {
        lasData.open();
        List<LasRecord> lasPoints = lasData.getPointsInGeometry(polygon, false);
        for( LasRecord lasRecord : lasPoints ) {
            lasCoordinates.add(new Coordinate(lasRecord.x, lasRecord.y, lasRecord.z));
        }
    }

    DelaunayTriangulationBuilder triangulationBuilder = new DelaunayTriangulationBuilder();
    triangulationBuilder.setSites(lasCoordinates);
    Geometry triangles = triangulationBuilder.getTriangles(gf);
    pm.done();

    int numTriangles = triangles.getNumGeometries();
    pm.beginTask("Extracting triangles based on threshold...", numTriangles);
    ArrayList<Geometry> trianglesList = new ArrayList<Geometry>();
    for( int i = 0; i < numTriangles; i++ ) {
        pm.worked(1);
        Geometry geometryN = triangles.getGeometryN(i);
        Coordinate[] coordinates = geometryN.getCoordinates();
        double diff1 = abs(coordinates[0].z - coordinates[1].z);
        if (diff1 > pElevThres) {
            continue;
        }
        double diff2 = abs(coordinates[0].z - coordinates[2].z);
        if (diff2 > pElevThres) {
            continue;
        }
        double diff3 = abs(coordinates[1].z - coordinates[2].z);
        if (diff3 > pElevThres) {
            continue;
        }
        trianglesList.add(geometryN);
    }
    pm.done();

    int newNumTriangles = trianglesList.size();
    int removedNum = numTriangles - newNumTriangles;
    pm.message("Original triangles: " + numTriangles);
    pm.message("New triangles: " + newNumTriangles);
    pm.message("Removed triangles: " + removedNum);

    pm.beginTask("Create triangles index...", newNumTriangles);
    final STRtree tree = new STRtree(trianglesList.size());
    for( Geometry triangle : trianglesList ) {
        Envelope env = triangle.getEnvelopeInternal();
        tree.insert(env, triangle);
        pm.worked(1);
    }
    pm.done();

    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();
    }
    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.createWritableRaster(newCols, newRows, null, null,
            HMConstants.doubleNovalue);

    ThreadedRunnable< ? > runner = new ThreadedRunnable(getDefaultThreadsNum(), null);
    pm.beginTask("Setting raster points...", newCols);
    for( int c = 0; c < newCols; c++ ) {
        final int fCol = c;
        runner.executeRunnable(new Runnable(){
            public void run() {
                try {
                    makeRow(tree, newRows, newGridGeometry2D, newWR, fCol);
                } catch (TransformException e) {
                    e.printStackTrace();
                }
                pm.worked(1);
            }
        });
    }
    runner.waitAndClose();
    pm.done();

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

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

    GridCoverage2D inDtmGC = getRaster(inDtm);
    Polygon polygon = CoverageUtilities.getRegionPolygon(inDtmGC);
    CoordinateReferenceSystem crs = inDtmGC.getCoordinateReferenceSystem();

    try (ALasDataManager lasData = ALasDataManager.getDataManager(new File(inIndexFile), inDtmGC, pThreshold, crs)) {
        lasData.open();
        if (pImpulse != null) {
            lasData.setImpulsesConstraint(new double[]{pImpulse});
        }

        List<LasRecord> lasPoints = lasData.getPointsInGeometry(polygon, false);
        if (lasPoints.size() == 0) {
            pm.message("No points foudn in the given area. Check your input.");
            return;
        }

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

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

        DefaultFeatureCollection newCollection = new DefaultFeatureCollection();
        final SimpleFeatureTypeBuilder b = new SimpleFeatureTypeBuilder();
        b.setName("lasdata");
        b.setCRS(crs);
        b.add("the_geom", Point.class);
        b.add("elev", Double.class);
        final SimpleFeatureType featureType = b.buildFeatureType();
        SimpleFeatureBuilder builder = new SimpleFeatureBuilder(featureType);

        pm.beginTask("Prepare points collection for interpolation...", lasPoints.size());
        for( LasRecord r : lasPoints ) {
            final Point point = gf.createPoint(new Coordinate(r.x, r.y));
            final Object[] values = new Object[]{point, r.z,};
            builder.addAll(values);
            final SimpleFeature feature = builder.buildFeature(null);
            newCollection.add(feature);
            pm.worked(1);
        }
        pm.done();

        OmsRasterGenerator omsRasterGenerator = new OmsRasterGenerator();
        omsRasterGenerator.pNorth = north;
        omsRasterGenerator.pSouth = south;
        omsRasterGenerator.pWest = west;
        omsRasterGenerator.pEast = east;
        omsRasterGenerator.pXres = (east - west) / newCols;
        omsRasterGenerator.pYres = (north - south) / newRows;
        omsRasterGenerator.inCrs = crs;
        omsRasterGenerator.doRandom = false;
        omsRasterGenerator.process();

        OmsSurfaceInterpolator idwInterpolator = new OmsSurfaceInterpolator();
        idwInterpolator.inVector = newCollection;
        idwInterpolator.inGrid = omsRasterGenerator.outRaster;
        idwInterpolator.inMask = null;
        idwInterpolator.fCat = "elev";
        idwInterpolator.pMode = Variables.IDW;
        idwInterpolator.pMaxThreads = getDefaultThreadsNum();
        idwInterpolator.pBuffer = 10.0;
        idwInterpolator.pm = pm;
        idwInterpolator.process();
        GridCoverage2D outRasterGC = idwInterpolator.outRaster;
        dumpRaster(outRasterGC, outRaster);

    }
}