Examples with IDWInterpolator org.hortonmachine.gears.modules.r.interpolation2d.core.IDWInterpolator used on opensource projects

Example 1 with IDWInterpolator

use of org.hortonmachine.gears.modules.r.interpolation2d.core.IDWInterpolator in project hortonmachine by TheHortonMachine.

the class OmsSurfaceInterpolator method process.

@Execute
public void process() throws Exception {
    checkNull(inGrid);
    gridGeometry = inGrid.getGridGeometry();
    RegionMap regionMap = CoverageUtilities.getRegionParamsFromGridCoverage(inGrid);
    final int cols = regionMap.getCols();
    int rows = regionMap.getRows();
    coordinatesSpatialTree = new STRtree();
    if (inVector != null) {
        checkNull(fCat);
        GeometryDescriptor geometryDescriptor = inVector.getSchema().getGeometryDescriptor();
        if (!EGeometryType.isPoint(geometryDescriptor)) {
            throw new ModelsIllegalargumentException("The geometry has to be a point geometry.", this, pm);
        }
        SimpleFeatureIterator featureIterator = inVector.features();
        Coordinate[] coordinates = new Coordinate[inVector.size()];
        int index = 0;
        pm.beginTask("Indexing control points...", coordinates.length);
        while (featureIterator.hasNext()) {
            SimpleFeature feature = featureIterator.next();
            Geometry geometry = (Geometry) feature.getDefaultGeometry();
            coordinates[index] = geometry.getCoordinate();
            double value = ((Number) feature.getAttribute(fCat)).doubleValue();
            coordinates[index].z = value;
            Envelope env = new Envelope(coordinates[index]);
            coordinatesSpatialTree.insert(env, coordinates[index]);
            pm.worked(1);
        }
        pm.done();
        pm.message("Indexed control points: " + coordinates.length);
    } else {
        // create it from grid
        pm.beginTask("Indexing control points...", cols);
        RandomIter inIter = CoverageUtilities.getRandomIterator(inGrid);
        int count = 0;
        for (int r = 0; r < rows; r++) {
            for (int c = 0; c < cols; c++) {
                double value = inIter.getSampleDouble(c, r, 0);
                if (!HMConstants.isNovalue(value)) {
                    Coordinate coordinate = CoverageUtilities.coordinateFromColRow(c, r, gridGeometry);
                    coordinate.z = value;
                    Envelope env = new Envelope(coordinate);
                    coordinatesSpatialTree.insert(env, coordinate);
                    count++;
                }
            }
            pm.worked(1);
        }
        pm.done();
        pm.message("Indexed control points (from input grid): " + count);
    }
    coordinatesSpatialTree.build();
    if (pMode.equals(IDW)) {
        interpolator = new IDWInterpolator(pBuffer);
    } else {
        interpolator = new TPSInterpolator(pBuffer);
    }
    WritableRaster interpolatedWR = CoverageUtilities.createWritableRaster(cols, rows, null, null, HMConstants.doubleNovalue);
    final WritableRandomIter interpolatedIter = RandomIterFactory.createWritable(interpolatedWR, null);
    boolean doMultiThread = pMaxThreads > 1;
    ExecutorService fixedThreadPool = null;
    if (doMultiThread)
        fixedThreadPool = Executors.newFixedThreadPool(pMaxThreads);
    pm.beginTask("Performing interpolation...", rows);
    final double[] eval = new double[1];
    for (int r = 0; r < rows; r++) {
        final int row = r;
        if (doMultiThread) {
            Runnable runner = new Runnable() {

                public void run() {
                    processing(cols, coordinatesSpatialTree, interpolatedIter, eval, row);
                }
            };
            fixedThreadPool.execute(runner);
        } else {
            processing(cols, coordinatesSpatialTree, interpolatedIter, eval, row);
        }
    }
    if (doMultiThread) {
        try {
            fixedThreadPool.shutdown();
            fixedThreadPool.awaitTermination(30, TimeUnit.DAYS);
            fixedThreadPool.shutdownNow();
        } catch (InterruptedException e) {
            e.printStackTrace();
        }
    }
    pm.done();
    outRaster = CoverageUtilities.buildCoverage("interpolatedraster", interpolatedWR, regionMap, inGrid.getCoordinateReferenceSystem());
}

Example 2 with IDWInterpolator

use of org.hortonmachine.gears.modules.r.interpolation2d.core.IDWInterpolator in project hortonmachine by TheHortonMachine.

the class OmsBobTheBuilder method process.

@Execute
public void process() throws Exception {
    checkNull(inRaster, inArea, inElevations, fElevation);
    RegionMap regionMap = CoverageUtilities.getRegionParamsFromGridCoverage(inRaster);
    int cols = regionMap.getCols();
    int rows = regionMap.getRows();
    double west = regionMap.getWest();
    double east = regionMap.getEast();
    double south = regionMap.getSouth();
    double north = regionMap.getNorth();
    ReferencedEnvelope vectorBounds = inArea.getBounds();
    if (!isBetween(vectorBounds.getMaxX(), west, east) || !isBetween(vectorBounds.getMinX(), west, east) || !isBetween(vectorBounds.getMaxY(), south, north) || !isBetween(vectorBounds.getMinY(), south, north)) {
        throw new ModelsIllegalargumentException("The vector map has to be within the raster map boundaries.", this, pm);
    }
    List<FeatureMate> polygonMates = FeatureUtilities.featureCollectionToMatesList(inArea);
    String polygonMessage = "This operation can be applied only to a single polygon.";
    if (polygonMates.size() != 1) {
        throw new ModelsIllegalargumentException(polygonMessage, this, pm);
    }
    FeatureMate polygonMate = polygonMates.get(0);
    Geometry polygon = polygonMate.getGeometry();
    if (polygon instanceof MultiPolygon) {
        polygon = polygon.getGeometryN(0);
    }
    if (!(polygon instanceof Polygon)) {
        throw new ModelsIllegalargumentException(polygonMessage, this, pm);
    }
    List<FeatureMate> pointsMates = FeatureUtilities.featureCollectionToMatesList(inElevations);
    if (doUseOnlyInternal) {
        PreparedGeometry preparedPolygon = PreparedGeometryFactory.prepare(polygon);
        List<FeatureMate> tmpPointsMates = new ArrayList<FeatureMate>();
        for (FeatureMate pointMate : pointsMates) {
            Geometry geometry = pointMate.getGeometry();
            if (preparedPolygon.covers(geometry)) {
                tmpPointsMates.add(pointMate);
            }
        }
        pointsMates = tmpPointsMates;
    }
    if (pointsMates.size() < 4) {
        throw new ModelsIllegalargumentException("You need at least 4 elevation points (the more, the better) to gain a decent interpolation.", this, pm);
    }
    List<Coordinate> controlPointsList = new ArrayList<Coordinate>();
    if (doPolygonborder) {
        pm.beginTask("Extract polygon border...", IHMProgressMonitor.UNKNOWN);
        Coordinate[] polygonCoordinates = polygon.getCoordinates();
        List<ProfilePoint> profile = CoverageUtilities.doProfile(inRaster, polygonCoordinates);
        for (ProfilePoint profilePoint : profile) {
            Coordinate position = profilePoint.getPosition();
            double elevation = profilePoint.getElevation();
            Coordinate coord = new Coordinate(position.x, position.y, elevation);
            controlPointsList.add(coord);
        }
        pm.done();
    }
    for (FeatureMate pointsMate : pointsMates) {
        Coordinate coordinate = pointsMate.getGeometry().getCoordinate();
        double elev = pointsMate.getAttribute(fElevation, Double.class);
        Coordinate coord = new Coordinate(coordinate.x, coordinate.y, elev);
        controlPointsList.add(coord);
    }
    Coordinate[] controlPoints = controlPointsList.toArray(new Coordinate[0]);
    GridGeometry2D gridGeometry = inRaster.getGridGeometry();
    RandomIter elevIter = CoverageUtilities.getRandomIterator(inRaster);
    WritableRaster outputWR = CoverageUtilities.createWritableRaster(cols, rows, null, null, HMConstants.doubleNovalue);
    WritableRandomIter outputIter = RandomIterFactory.createWritable(outputWR, null);
    DefaultFeatureCollection newCollection = new DefaultFeatureCollection();
    newCollection.add(polygonMate.getFeature());
    OmsScanLineRasterizer slRasterizer = new OmsScanLineRasterizer();
    slRasterizer.pm = pm;
    slRasterizer.inVector = newCollection;
    slRasterizer.pCols = cols;
    slRasterizer.pRows = rows;
    slRasterizer.pNorth = north;
    slRasterizer.pSouth = south;
    slRasterizer.pEast = east;
    slRasterizer.pWest = west;
    slRasterizer.pValue = 1.0;
    slRasterizer.process();
    GridCoverage2D outRasterized = slRasterizer.outRaster;
    if (pMaxbuffer < 0)
        pMaxbuffer = Math.max(vectorBounds.getWidth(), vectorBounds.getHeight());
    IDWInterpolator interpolator = new IDWInterpolator(pMaxbuffer);
    final GridCoordinates2D gridCoord = new GridCoordinates2D();
    RandomIter rasterizedIter = CoverageUtilities.getRandomIterator(outRasterized);
    pm.beginTask("Interpolating...", cols);
    for (int r = 0; r < rows; r++) {
        for (int c = 0; c < cols; c++) {
            double probValue = rasterizedIter.getSampleDouble(c, r, 0);
            if (isNovalue(probValue)) {
                continue;
            }
            gridCoord.setLocation(c, r);
            DirectPosition world = gridGeometry.gridToWorld(gridCoord);
            double[] coordinate = world.getCoordinate();
            double interpolated = interpolator.getValue(controlPoints, new Coordinate(coordinate[0], coordinate[1]));
            outputIter.setSample(c, r, 0, interpolated);
        }
        pm.worked(1);
    }
    pm.done();
    pm.beginTask("Merging with original raster...", cols);
    for (int r = 0; r < rows; r++) {
        for (int c = 0; c < cols; c++) {
            double interpolatedValue = outputIter.getSampleDouble(c, r, 0);
            double rasterValue = elevIter.getSampleDouble(c, r, 0);
            if (isNovalue(interpolatedValue)) {
                if (!isNovalue(rasterValue))
                    outputIter.setSample(c, r, 0, rasterValue);
            } else {
                if (doErode) {
                    // any value generated is ok
                    outputIter.setSample(c, r, 0, interpolatedValue);
                } else {
                    // only values higher than the raster are ok
                    if (!isNovalue(rasterValue)) {
                        if (rasterValue < interpolatedValue) {
                            outputIter.setSample(c, r, 0, interpolatedValue);
                        } else {
                            outputIter.setSample(c, r, 0, rasterValue);
                        }
                    } else {
                        outputIter.setSample(c, r, 0, interpolatedValue);
                    }
                }
            }
        }
        pm.worked(1);
    }
    pm.done();
    outRaster = CoverageUtilities.buildCoverage("raster", outputWR, regionMap, inRaster.getCoordinateReferenceSystem());
}

Example 3 with IDWInterpolator

use of org.hortonmachine.gears.modules.r.interpolation2d.core.IDWInterpolator in project hortonmachine by TheHortonMachine.

the class OmsHoleFiller method process.

@Execute
public void process() throws Exception {
    checkNull(inRaster);
    ISurfaceInterpolator interpolator;
    if (pMode.equals(IDW)) {
        interpolator = new IDWInterpolator(pBuffer);
    } else {
        interpolator = new TPSInterpolator(pBuffer);
    }
    RegionMap regionMap = CoverageUtilities.getRegionParamsFromGridCoverage(inRaster);
    int rows = regionMap.getRows();
    int cols = regionMap.getCols();
    WritableRaster outWR = CoverageUtilities.renderedImage2WritableRaster(inRaster.getRenderedImage(), false);
    WritableRandomIter outIter = CoverageUtilities.getWritableRandomIterator(outWR);
    GridGeometry2D gridGeometry = inRaster.getGridGeometry();
    PreparedGeometry preparedRoi = null;
    if (inROI != null) {
        List<Geometry> roiList = FeatureUtilities.featureCollectionToGeometriesList(inROI, false, null);
        GeometryCollection gc = new GeometryCollection(roiList.toArray(GeometryUtilities.TYPE_GEOMETRY), gf);
        preparedRoi = PreparedGeometryFactory.prepare(gc);
    }
    pm.beginTask("Filling holes...", cols - 2);
    for (int r = 1; r < rows - 1; r++) {
        for (int c = 1; c < cols - 1; c++) {
            if (pm.isCanceled()) {
                return;
            }
            double value = outIter.getSampleDouble(c, r, 0);
            if (isNovalue(value)) {
                DirectPosition worldPosition = gridGeometry.gridToWorld(new GridCoordinates2D(c, r));
                double[] coordinate = worldPosition.getCoordinate();
                Coordinate pointCoordinate = new Coordinate(coordinate[0], coordinate[1]);
                Point point = gf.createPoint(pointCoordinate);
                if (preparedRoi == null || preparedRoi.intersects(point)) {
                    // TODO this could be done considering more points and more far away points.
                    // For now, this works.
                    List<Coordinate> surroundingValids = getValidSurroundingPoints(outIter, gridGeometry, c, r);
                    if (surroundingValids.size() > 3) {
                        double newValue = interpolator.getValue(surroundingValids.toArray(new Coordinate[0]), pointCoordinate);
                        outIter.setSample(c, r, 0, newValue);
                    }
                }
            }
        }
        pm.worked(1);
    }
    pm.done();
    outIter.done();
    outRaster = CoverageUtilities.buildCoverage("nulled", outWR, regionMap, inRaster.getCoordinateReferenceSystem());
}