Examples with RandomIter javax.media.jai.iterator.RandomIter used on opensource projects

Example 71 with RandomIter

use of javax.media.jai.iterator.RandomIter in project hortonmachine by TheHortonMachine.

the class OmsShalstab method qcrit.

/**
 * Calculates the trasmissivity in every pixel of the map.
 */
private void qcrit(RenderedImage slope, RenderedImage ab, RandomIter trasmissivityRI, RandomIter frictionRI, RandomIter cohesionRI, RandomIter hsIter, RandomIter effectiveRI, RandomIter densityRI) {
    HashMap<String, Double> regionMap = CoverageUtilities.getRegionParamsFromGridCoverage(inSlope);
    int cols = regionMap.get(CoverageUtilities.COLS).intValue();
    int rows = regionMap.get(CoverageUtilities.ROWS).intValue();
    RandomIter slopeRI = RandomIterFactory.create(slope, null);
    RandomIter abRI = RandomIterFactory.create(ab, null);
    WritableRaster qcritWR = CoverageUtilities.createWritableRaster(cols, rows, null, null, null);
    WritableRandomIter qcritIter = RandomIterFactory.createWritable(qcritWR, null);
    WritableRaster classiWR = CoverageUtilities.createWritableRaster(cols, rows, null, null, null);
    WritableRandomIter classiIter = RandomIterFactory.createWritable(classiWR, null);
    pm.beginTask("Creating qcrit map...", rows);
    for (int j = 0; j < rows; j++) {
        pm.worked(1);
        for (int i = 0; i < cols; i++) {
            double slopeValue = slopeRI.getSampleDouble(i, j, 0);
            double tanPhiValue = frictionRI.getSampleDouble(i, j, 0);
            double cohValue = cohesionRI.getSampleDouble(i, j, 0);
            double rhoValue = densityRI.getSampleDouble(i, j, 0);
            double hsValue = hsIter.getSampleDouble(i, j, 0);
            if (!isNovalue(slopeValue) && !isNovalue(tanPhiValue) && !isNovalue(cohValue) && !isNovalue(rhoValue) && !isNovalue(hsValue)) {
                if (hsValue <= EPS || slopeValue > pRock) {
                    qcritIter.setSample(i, j, 0, ROCK);
                } else {
                    double checkUnstable = tanPhiValue + cohValue / (9810.0 * rhoValue * hsValue) * (1 + pow(slopeValue, 2));
                    if (slopeValue >= checkUnstable) {
                        /*
                             * uncond unstable
                             */
                        qcritIter.setSample(i, j, 0, 5);
                    } else {
                        double checkStable = tanPhiValue * (1 - 1 / rhoValue) + cohValue / (9810 * rhoValue * hsValue) * (1 + pow(slopeValue, 2));
                        if (slopeValue < checkStable) {
                            /*
                                 * uncond. stable
                                 */
                            qcritIter.setSample(i, j, 0, 0);
                        } else {
                            double qCrit = trasmissivityRI.getSampleDouble(i, j, 0) * sin(atan(slopeValue)) / abRI.getSampleDouble(i, j, 0) * rhoValue * (1 - slopeValue / tanPhiValue + cohValue / (9810 * rhoValue * hsValue * tanPhiValue) * (1 + pow(slopeValue, 2))) * 1000;
                            qcritIter.setSample(i, j, 0, qCrit);
                            /*
                                 * see the Qcrit (critical effective
                                 * precipitation) that leads the slope to
                                 * instability (see article of Montgomery et Al,
                                 * Hydrological Processes, 12, 943-955, 1998)
                                 */
                            double value = qcritIter.getSampleDouble(i, j, 0);
                            if (value > 0 && value < 50)
                                qcritIter.setSample(i, j, 0, 1);
                            if (value >= 50 && value < 100)
                                qcritIter.setSample(i, j, 0, 2);
                            if (value >= 100 && value < 200)
                                qcritIter.setSample(i, j, 0, 3);
                            if (value >= 200)
                                qcritIter.setSample(i, j, 0, 4);
                        }
                    }
                }
            } else {
                qcritIter.setSample(i, j, 0, doubleNovalue);
            }
        }
    }
    pm.done();
    /*
         * build the class matrix 1=inc inst 2=inc stab 3=stab 4=instab
         * rock=presence of rock
         */
    pm.beginTask("Creating stability map...", rows);
    double Tq = 0;
    for (int j = 0; j < rows; j++) {
        pm.worked(1);
        for (int i = 0; i < cols; i++) {
            Tq = trasmissivityRI.getSampleDouble(i, j, 0) / (effectiveRI.getSampleDouble(i, j, 0) / 1000.0);
            double slopeValue = slopeRI.getSampleDouble(i, j, 0);
            double abValue = abRI.getSampleDouble(i, j, 0);
            double tangPhiValue = frictionRI.getSampleDouble(i, j, 0);
            double cohValue = cohesionRI.getSampleDouble(i, j, 0);
            double rhoValue = densityRI.getSampleDouble(i, j, 0);
            double hsValue = hsIter.getSampleDouble(i, j, 0);
            if (!isNovalue(slopeValue) && !isNovalue(abValue) && !isNovalue(tangPhiValue) && !isNovalue(cohValue) && !isNovalue(rhoValue) && !isNovalue(hsValue)) {
                if (hsValue <= EPS || slopeValue > pRock) {
                    classiIter.setSample(i, j, 0, ROCK);
                } else {
                    double checkUncondUnstable = tangPhiValue + cohValue / (9810 * rhoValue * hsValue) * (1 + pow(slopeValue, 2));
                    double checkUncondStable = tangPhiValue * (1 - 1 / rhoValue) + cohValue / (9810 * rhoValue * hsValue) * (1 + pow(slopeValue, 2));
                    double checkStable = Tq * sin(atan(slopeValue)) * rhoValue * (1 - slopeValue / tangPhiValue + cohValue / (9810 * rhoValue * hsValue * tangPhiValue) * (1 + pow(slopeValue, 2)));
                    if (slopeValue >= checkUncondUnstable) {
                        classiIter.setSample(i, j, 0, 1);
                    } else if (slopeValue < checkUncondStable) {
                        classiIter.setSample(i, j, 0, 2);
                    } else if (abValue < checkStable && classiIter.getSampleDouble(i, j, 0) != 1 && classiIter.getSampleDouble(i, j, 0) != 2) {
                        classiIter.setSample(i, j, 0, 3);
                    } else {
                        classiIter.setSample(i, j, 0, 4);
                    }
                }
            } else {
                classiIter.setSample(i, j, 0, doubleNovalue);
            }
        }
    }
    pm.done();
    outQcrit = CoverageUtilities.buildCoverage("qcrit", qcritWR, regionMap, inSlope.getCoordinateReferenceSystem());
    outShalstab = CoverageUtilities.buildCoverage("classi", classiWR, regionMap, inSlope.getCoordinateReferenceSystem());
}

Example 72 with RandomIter

use of javax.media.jai.iterator.RandomIter in project hortonmachine by TheHortonMachine.

the class CoverageUtilities method coverageValuesMapper.

/**
 * Mappes the values of a map (valuesMap) into the valid pixels of the second map (maskMap).
 *
 * @param valuesMap the map holding the values that are needed in the resulting map.
 * @param maskMap the map to use as mask for the values.
 * @return the map containing the values of the valuesMap, but only in the places in which the maskMap is valid.
 */
public static GridCoverage2D coverageValuesMapper(GridCoverage2D valuesMap, GridCoverage2D maskMap) {
    RegionMap valuesRegionMap = getRegionParamsFromGridCoverage(valuesMap);
    int cs = valuesRegionMap.getCols();
    int rs = valuesRegionMap.getRows();
    RegionMap maskRegionMap = getRegionParamsFromGridCoverage(maskMap);
    int tmpcs = maskRegionMap.getCols();
    int tmprs = maskRegionMap.getRows();
    if (cs != tmpcs || rs != tmprs) {
        throw new IllegalArgumentException("The raster maps have to be of equal size to be mapped.");
    }
    RandomIter valuesIter = RandomIterFactory.create(valuesMap.getRenderedImage(), null);
    RandomIter maskIter = RandomIterFactory.create(maskMap.getRenderedImage(), null);
    WritableRaster writableRaster = createWritableRaster(cs, rs, null, null, HMConstants.doubleNovalue);
    WritableRandomIter outIter = RandomIterFactory.createWritable(writableRaster, null);
    for (int r = 0; r < rs; r++) {
        for (int c = 0; c < cs; c++) {
            if (!isNovalue(maskIter.getSampleDouble(c, r, 0))) {
                // if not nv, put the value from the valueMap in the new map
                double value = valuesIter.getSampleDouble(c, r, 0);
                if (!isNovalue(value))
                    outIter.setSample(c, r, 0, value);
            }
        }
    }
    GridCoverage2D outCoverage = buildCoverage(// $NON-NLS-1$
    "mapped", // $NON-NLS-1$
    writableRaster, // $NON-NLS-1$
    maskRegionMap, valuesMap.getCoordinateReferenceSystem());
    return outCoverage;
}

Example 73 with RandomIter

use of javax.media.jai.iterator.RandomIter in project hortonmachine by TheHortonMachine.

the class CoverageUtilities method calculateHypsographic.

/**
 * Calculates the hypsographic curve for the given raster, using the supplied bins.
 *
 * @param elevationCoverage the elevation raster.
 * @param bins the bins to use.
 * @param pm the monitor.
 * @return the matrix containing the hypsographic curve in [elev, area] pairs per row.
 */
public static double[][] calculateHypsographic(GridCoverage2D elevationCoverage, int bins, IHMProgressMonitor pm) {
    if (pm == null) {
        pm = new DummyProgressMonitor();
    }
    RegionMap regionMap = getRegionParamsFromGridCoverage(elevationCoverage);
    int cols = regionMap.getCols();
    int rows = regionMap.getRows();
    double xres = regionMap.getXres();
    double yres = regionMap.getYres();
    /*
         * calculate the maximum and minimum value of the raster data
         */
    RandomIter elevIter = getRandomIterator(elevationCoverage);
    double maxRasterValue = Double.NEGATIVE_INFINITY;
    double minRasterValue = Double.POSITIVE_INFINITY;
    for (int r = 0; r < rows; r++) {
        for (int c = 0; c < cols; c++) {
            double value = elevIter.getSampleDouble(c, r, 0);
            if (isNovalue(value)) {
                continue;
            }
            maxRasterValue = max(maxRasterValue, value);
            minRasterValue = min(minRasterValue, value);
        }
    }
    /*
         * subdivide the whole value range in bins and count the number of pixels in each bin
         */
    double binWidth = (maxRasterValue - minRasterValue) / (bins);
    double[] pixelPerBinCount = new double[bins];
    double[] areaAtGreaterElevation = new double[bins];
    pm.beginTask("Performing calculation of hypsographic curve...", rows);
    for (int r = 0; r < rows; r++) {
        for (int c = 0; c < cols; c++) {
            double value = elevIter.getSampleDouble(c, r, 0);
            if (isNovalue(value)) {
                continue;
            }
            for (int k = 0; k < pixelPerBinCount.length; k++) {
                double thres = minRasterValue + k * binWidth;
                if (value >= thres) {
                    pixelPerBinCount[k] = pixelPerBinCount[k] + 1;
                    areaAtGreaterElevation[k] = areaAtGreaterElevation[k] + (yres * xres);
                }
            }
        }
        pm.worked(1);
    }
    pm.done();
    double[][] hypso = new double[pixelPerBinCount.length][3];
    for (int j = 0; j < hypso.length; j++) {
        hypso[j][0] = minRasterValue + (j * binWidth) + (binWidth / 2.0);
        hypso[j][1] = areaAtGreaterElevation[j] / 1000000.0;
    }
    return hypso;
}

Example 74 with RandomIter

use of javax.media.jai.iterator.RandomIter in project hortonmachine by TheHortonMachine.

the class CoverageUtilities method gridcoverageToCellPolygons.

/**
 * Extracts a list of polygons from the cell bounds of a given {@link GridCoverage2D coverage}.
 *
 * <p><b>Note that the cells are added in a rows
 * and cols order (for each row evaluate each column).</b></p>
 *
 * <p>The userdata of the geometry contains the value of the raster.
 *
 * @param coverage the coverage to use.
 * @param keepCoordinatePredicate an optional predicate to filter out some of the cells.
 * @return the list of envelope geometries.
 */
public static List<Polygon> gridcoverageToCellPolygons(GridCoverage2D coverage, Predicate<Coordinate> keepCoordinatePredicate, boolean doIncrementalMerge, IHMProgressMonitor pm) {
    if (pm == null) {
        pm = new DummyProgressMonitor();
    }
    RegionMap regionMap = CoverageUtilities.getRegionParamsFromGridCoverage(coverage);
    double west = regionMap.getWest();
    double north = regionMap.getNorth();
    double xres = regionMap.getXres();
    double yres = regionMap.getYres();
    int cols = regionMap.getCols();
    int rows = regionMap.getRows();
    int everyRows = 1000;
    GeometryFactory gf = GeometryUtilities.gf();
    RandomIter iter = CoverageUtilities.getRandomIterator(coverage);
    List<Geometry> bigPolygons = new ArrayList<Geometry>();
    List<Polygon> polygons = new ArrayList<Polygon>();
    pm.beginTask("Vectorizing raster cells...", rows);
    for (int r = 0; r < rows; r++) {
        for (int c = 0; c < cols; c++) {
            double w = west + xres * c;
            double e = w + xres;
            double n = north - yres * r;
            double s = n - yres;
            if (keepCoordinatePredicate != null && !keepCoordinatePredicate.test(new Coordinate(w + xres / 2, s + yres / 2))) {
                continue;
            }
            Coordinate[] coords = new Coordinate[5];
            coords[0] = new Coordinate(w, n);
            coords[1] = new Coordinate(e, n);
            coords[2] = new Coordinate(e, s);
            coords[3] = new Coordinate(w, s);
            coords[4] = new Coordinate(w, n);
            LinearRing linearRing = gf.createLinearRing(coords);
            Polygon polygon = gf.createPolygon(linearRing, null);
            polygons.add(polygon);
            double value = iter.getSampleDouble(c, r, 0);
            polygon.setUserData(value);
        }
        if (doIncrementalMerge && ((r > 0 && r % everyRows == 0) || r == rows - 1)) {
            // merge same value to single geometries to save memory
            Map<Integer, List<Geometry>> collected = polygons.parallelStream().filter(poly -> ((Number) poly.getUserData()).doubleValue() != HMConstants.doubleNovalue).collect(Collectors.groupingBy(poly -> ((Number) poly.getUserData()).intValue()));
            pm.message(r + " of " + rows);
            polygons = new ArrayList<Polygon>();
            int size = collected.size();
            int count = 0;
            for (Entry<Integer, List<Geometry>> entry : collected.entrySet()) {
                count++;
                if (count % 10 == 0)
                    pm.message("   -> " + count + " of " + size);
                Integer basinId = entry.getKey();
                List<Geometry> value = entry.getValue();
                Geometry tmpGeom = CascadedPolygonUnion.union(value);
                for (int i = 0; i < tmpGeom.getNumGeometries(); i++) {
                    Polygon geometryN = (Polygon) tmpGeom.getGeometryN(i);
                    geometryN.setUserData(basinId);
                    bigPolygons.add(geometryN);
                }
            }
        }
        pm.worked(1);
    }
    pm.done();
    if (doIncrementalMerge) {
        polygons = new ArrayList<Polygon>();
        for (Geometry tmpGeom : bigPolygons) {
            Object userData = tmpGeom.getUserData();
            for (int i = 0; i < tmpGeom.getNumGeometries(); i++) {
                Polygon geometryN = (Polygon) tmpGeom.getGeometryN(i);
                geometryN.setUserData(userData);
                polygons.add(geometryN);
            }
        }
    }
    return polygons;
}

Example 75 with RandomIter

use of javax.media.jai.iterator.RandomIter in project hortonmachine by TheHortonMachine.

the class CoverageUtilities method invert.

public static GridCoverage2D invert(GridCoverage2D raster, double max) {
    RegionMap regionMap = CoverageUtilities.getRegionParamsFromGridCoverage(raster);
    int nCols = regionMap.getCols();
    int nRows = regionMap.getRows();
    RandomIter rasterIter = CoverageUtilities.getRandomIterator(raster);
    WritableRaster outWR = CoverageUtilities.createWritableRaster(nCols, nRows, null, null, doubleNovalue);
    WritableRandomIter outIter = RandomIterFactory.createWritable(outWR, null);
    for (int r = 0; r < nRows; r++) {
        for (int c = 0; c < nCols; c++) {
            double value = rasterIter.getSampleDouble(c, r, 0);
            if (!isNovalue(value)) {
                outIter.setSample(c, r, 0, max - value);
            }
        }
    }
    GridCoverage2D invertedRaster = CoverageUtilities.buildCoverage("inverted", outWR, regionMap, raster.getCoordinateReferenceSystem());
    return invertedRaster;
}