Search in sources :

Example 51 with WritableRandomIter

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

the class OmsHillshade method normalVector.

protected WritableRaster normalVector(WritableRaster pitWR, double res) {
    int minX = pitWR.getMinX();
    int minY = pitWR.getMinY();
    int rows = pitWR.getHeight();
    int cols = pitWR.getWidth();
    RandomIter pitIter = RandomIterFactory.create(pitWR, null);
    /*
         * Initialize the Image of the normal vector in the central point of the
         * cells, which have 3 components so the Image have 3 bands..
         */
    SampleModel sm = RasterFactory.createBandedSampleModel(5, cols, rows, 3);
    WritableRaster tmpNormalVectorWR = CoverageUtilities.createWritableRaster(cols, rows, null, sm, 0.0);
    WritableRandomIter tmpNormaIter = RandomIterFactory.createWritable(tmpNormalVectorWR, null);
    /*
         * apply the corripio's formula (is the formula (3) in the article)
         */
    for (int j = minY; j < minX + rows - 1; j++) {
        for (int i = minX; i < minX + cols - 1; i++) {
            double zij = pitIter.getSampleDouble(i, j, 0);
            double zidxj = pitIter.getSampleDouble(i + 1, j, 0);
            double zijdy = pitIter.getSampleDouble(i, j + 1, 0);
            double zidxjdy = pitIter.getSampleDouble(i + 1, j + 1, 0);
            double firstComponent = 0.5 * res * (zij - zidxj + zijdy - zidxjdy);
            double secondComponent = 0.5 * res * (zij + zidxj - zijdy - zidxjdy);
            double thirthComponent = (res * res);
            double den = Math.sqrt(firstComponent * firstComponent + secondComponent * secondComponent + thirthComponent * thirthComponent);
            tmpNormaIter.setPixel(i, j, new double[] { firstComponent / den, secondComponent / den, thirthComponent / den });
        }
    }
    pitIter.done();
    return tmpNormalVectorWR;
}
Also used : SampleModel(java.awt.image.SampleModel) WritableRandomIter(javax.media.jai.iterator.WritableRandomIter) WritableRaster(java.awt.image.WritableRaster) RandomIter(javax.media.jai.iterator.RandomIter) WritableRandomIter(javax.media.jai.iterator.WritableRandomIter)

Example 52 with WritableRandomIter

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

the class OmsDebrisFlow method process.

@Execute
public void process() throws Exception {
    RegionMap regionMap = CoverageUtilities.getRegionParamsFromGridCoverage(inElev);
    int cols = regionMap.getCols();
    int rows = regionMap.getRows();
    double xRes = regionMap.getXres();
    double yRes = regionMap.getYres();
    double west = regionMap.getWest();
    double east = regionMap.getEast();
    double south = regionMap.getSouth();
    double north = regionMap.getNorth();
    if (!isBetween(pEasting, west, east) || !isBetween(pNorthing, south, north)) {
        throw new ModelsIllegalargumentException("Input coordinates have to be within the map boundaries.", this, pm);
    }
    double thresArea = pMcoeff * pow(pVolume, (2.0 / 3.0));
    GridGeometry2D gridGeometry = inElev.getGridGeometry();
    int[] colRow = CoverageUtilities.colRowFromCoordinate(new Coordinate(pEasting, pNorthing), gridGeometry, null);
    RandomIter elevIter = CoverageUtilities.getRandomIterator(inElev);
    int startCol = colRow[0];
    int startRow = colRow[1];
    double startValue = elevIter.getSampleDouble(startCol, startRow, 0);
    if (isNovalue(startValue)) {
        throw new ModelsIllegalargumentException("Input coordinates are on a novalue elevation point.", this, pm);
    }
    WritableRaster mcsWR = CoverageUtilities.createWritableRaster(cols, rows, null, null, HMConstants.doubleNovalue);
    WritableRandomIter probIter = RandomIterFactory.createWritable(mcsWR, null);
    Random flatRnd = new Random();
    int processedMc = 0;
    for (int mc = 0; mc < pMontecarlo; mc++) {
        pm.message("Montecarlo n." + mc);
        processedMc = mc;
        /*
             * for every Montecarlo loop, get a flow path
             */
        int centerCol = startCol;
        int centerRow = startRow;
        TreeSet<Point> touchedPoints = new TreeSet<Point>();
        touchedPoints.add(new Point(centerCol, centerRow));
        boolean doStop = false;
        // cicle over every cell neighbours along the way
        // -System.currentTimeMillis());
        Random randomGenerator = new Random();
        do {
            // System.out.println(centerCol + "/" + centerRow + " --- " + cols + "/" + rows);
            double centerValue = elevIter.getSampleDouble(centerCol, centerRow, 0);
            List<SlopeProbability> spList = new ArrayList<SlopeProbability>();
            double slopeSum = 0;
            for (int x = -1; x <= 1; x++) {
                for (int y = -1; y <= 1; y++) {
                    if (x == 0 && y == 0) {
                        continue;
                    }
                    int tmpCol = centerCol + x;
                    int tmpRow = centerRow + y;
                    if (touchedPoints.contains(new Point(tmpCol, tmpRow))) {
                        continue;
                    }
                    // if point is outside jump it
                    if (!isBetween(tmpCol, 0, cols - 1) || !isBetween(tmpRow, 0, rows - 1)) {
                        continue;
                    }
                    // if point is novalue, jump it
                    double nextValue = elevIter.getSampleDouble(tmpCol, tmpRow, 0);
                    if (isNovalue(nextValue)) {
                        continue;
                    }
                    double distance = pythagoras(abs((tmpCol - centerCol) * xRes), abs((tmpRow - centerRow) * yRes));
                    double slope = (nextValue - centerValue) / distance;
                    // we take only negative and 0 slope, downhill
                    if (slope > 0) {
                        continue;
                    }
                    slope = abs(slope);
                    SlopeProbability sp = new SlopeProbability();
                    sp.fromCol = centerCol;
                    sp.fromRow = centerRow;
                    sp.fromElev = centerValue;
                    sp.toCol = tmpCol;
                    sp.toRow = tmpRow;
                    sp.toElev = nextValue;
                    sp.slope = slope;
                    slopeSum = slopeSum + slope;
                    spList.add(sp);
                }
            }
            if (spList.size() == 0) {
                /*
                     * touched border or slope is not negative
                     */
                doStop = true;
            } else {
                // get a random number between 0 and 1
                double random = randomGenerator.nextDouble();
                if (spList.size() == 1) {
                    // direction is only one
                    SlopeProbability sp = spList.get(0);
                    centerCol = sp.toCol;
                    centerRow = sp.toRow;
                } else {
                    Collections.sort(spList);
                    /*
                         * case in which the slopes are all 0
                         */
                    if (dEq(slopeSum, 0.0)) {
                        // choose a random and go on
                        int size = spList.size();
                        double rnd = flatRnd.nextDouble();
                        int index = (int) round(rnd * size) - 1;
                        if (index < 0)
                            index = 0;
                        SlopeProbability sp = spList.get(index);
                        centerCol = sp.toCol;
                        centerRow = sp.toRow;
                    } else /*
                         * normal case in which the slopes have a value
                         */
                    {
                        // cumulate the probability
                        for (int i = 0; i < spList.size(); i++) {
                            SlopeProbability sp = spList.get(i);
                            double p = sp.slope / slopeSum;
                            sp.probability = p;
                            if (i != 0) {
                                SlopeProbability tmpSp = spList.get(i - 1);
                                sp.probability = sp.probability + tmpSp.probability;
                            }
                        }
                        for (int i = 1; i < spList.size(); i++) {
                            SlopeProbability sp1 = spList.get(i - 1);
                            SlopeProbability sp2 = spList.get(i);
                            // if (random < sp1.probability) {
                            if (random < sp1.probability) {
                                centerCol = sp1.toCol;
                                centerRow = sp1.toRow;
                                break;
                            // } else if (random >= sp1.probability && random <
                            // sp2.probability)
                            // {
                            } else if (random >= sp1.probability && random < sp2.probability) {
                                centerCol = sp2.toCol;
                                centerRow = sp2.toRow;
                                break;
                            }
                        }
                    }
                }
                touchedPoints.add(new Point(centerCol, centerRow));
                double outValue = probIter.getSampleDouble(centerCol, centerRow, 0);
                if (isNovalue(outValue)) {
                    outValue = 0.0;
                }
                probIter.setSample(centerCol, centerRow, 0, outValue + 1.0);
            }
        } while (!doStop);
        /*
             * check if the max area is flooded
             */
        int floodedCellNum = 0;
        for (int r = 0; r < rows; r++) {
            for (int c = 0; c < cols; c++) {
                double value = probIter.getSampleDouble(c, r, 0);
                if (isNovalue(value)) {
                    continue;
                }
                floodedCellNum++;
            }
        }
        double floodedArea = floodedCellNum * xRes * yRes;
        if (thresArea <= floodedArea) {
            break;
        }
    }
    double probSum = 0.0;
    double validCells = 0.0;
    for (int r = 0; r < rows; r++) {
        for (int c = 0; c < cols; c++) {
            double prob = probIter.getSampleDouble(c, r, 0);
            if (isNovalue(prob)) {
                continue;
            }
            double newProb = prob / (processedMc - 1);
            probIter.setSample(c, r, 0, newProb);
            probSum = probSum + sqrt(newProb);
            validCells++;
        }
    }
    /*
         * calculate deposition
         */
    double avgProb = probSum / validCells;
    double avgHeight = pDcoeff * pow(pVolume, 1.0 / 3.0);
    WritableRaster depoWR = CoverageUtilities.createWritableRaster(cols, rows, null, null, HMConstants.doubleNovalue);
    WritableRandomIter depoIter = RandomIterFactory.createWritable(depoWR, null);
    for (int r = 0; r < rows; r++) {
        for (int c = 0; c < cols; c++) {
            double probValue = probIter.getSampleDouble(c, r, 0);
            if (isNovalue(probValue)) {
                continue;
            }
            double depoValue = avgHeight * sqrt(probValue) / avgProb;
            depoIter.setSample(c, r, 0, depoValue);
        }
    }
    outMcs = CoverageUtilities.buildCoverage("mcs", mcsWR, regionMap, inElev.getCoordinateReferenceSystem());
    outDepo = CoverageUtilities.buildCoverage("depo", depoWR, regionMap, inElev.getCoordinateReferenceSystem());
}
Also used : GridGeometry2D(org.geotools.coverage.grid.GridGeometry2D) RegionMap(org.hortonmachine.gears.utils.RegionMap) ArrayList(java.util.ArrayList) RandomIter(javax.media.jai.iterator.RandomIter) WritableRandomIter(javax.media.jai.iterator.WritableRandomIter) ModelsIllegalargumentException(org.hortonmachine.gears.libs.exceptions.ModelsIllegalargumentException) WritableRandomIter(javax.media.jai.iterator.WritableRandomIter) Random(java.util.Random) Coordinate(org.locationtech.jts.geom.Coordinate) WritableRaster(java.awt.image.WritableRaster) TreeSet(java.util.TreeSet) Execute(oms3.annotations.Execute)

Example 53 with WritableRandomIter

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

the class OmsDebrisTriggerCnr method process.

@Execute
public void process() throws Exception {
    checkNull(inElev, inNet, inTca);
    // calculate gradient map degrees
    OmsGradient gradient = new OmsGradient();
    gradient.inElev = inElev;
    gradient.pMode = Variables.FINITE_DIFFERENCES;
    gradient.doDegrees = true;
    gradient.pm = pm;
    gradient.process();
    GridCoverage2D gradientCoverageDeg = gradient.outSlope;
    // calculate gradient map %
    gradient = new OmsGradient();
    gradient.inElev = inElev;
    gradient.pMode = Variables.FINITE_DIFFERENCES;
    gradient.doDegrees = false;
    gradient.pm = pm;
    gradient.process();
    GridCoverage2D gradientCoverageTan = gradient.outSlope;
    // ritaglio della mappa di gradient lungo il reticolo
    // idrografico ed estrazione delle sole celle con
    // * pendenza minore di 38 gradi
    // * area cumulata minore di 10 km2
    RegionMap regionMap = CoverageUtilities.getRegionParamsFromGridCoverage(inElev);
    int cols = regionMap.getCols();
    int rows = regionMap.getRows();
    double xres = regionMap.getXres();
    double yres = regionMap.getYres();
    RenderedImage netRI = inNet.getRenderedImage();
    RandomIter netIter = RandomIterFactory.create(netRI, null);
    RenderedImage tcaRI = inTca.getRenderedImage();
    RandomIter tcaIter = RandomIterFactory.create(tcaRI, null);
    RenderedImage gradientDegRI = gradientCoverageDeg.getRenderedImage();
    RandomIter gradientDegIter = RandomIterFactory.create(gradientDegRI, null);
    RenderedImage gradientTanRI = gradientCoverageTan.getRenderedImage();
    RandomIter gradientTanIter = RandomIterFactory.create(gradientTanRI, null);
    WritableRaster outputWR = CoverageUtilities.createWritableRaster(cols, rows, null, null, HMConstants.doubleNovalue);
    WritableRandomIter outputIter = RandomIterFactory.createWritable(outputWR, null);
    pm.beginTask("Extracting trigger points...", cols);
    for (int r = 0; r < rows; r++) {
        for (int c = 0; c < cols; c++) {
            double net = netIter.getSampleDouble(c, r, 0);
            // all only along the network
            if (!isNovalue(net)) {
                double tca = tcaIter.getSampleDouble(c, r, 0);
                // tca in km2 along the net
                double tcaKm2 = tca * xres * yres / 1000000;
                // gradient in degrees along the net
                double gradientDeg = gradientDegIter.getSampleDouble(c, r, 0);
                // gradient in tan along the net
                double gradientTan = gradientTanIter.getSampleDouble(c, r, 0);
                /*
                     * calculate the trigger threshold:
                     * 
                     *  S = 0.32 * A^-0.2
                     *  where:
                     *   S = gradient in m/m
                     *   A = tca in km2
                     */
                double triggerThreshold = 0.32 * pow(tcaKm2, -0.2);
                if (// 
                gradientTan > triggerThreshold && // 
                gradientDeg < pGradthres && tcaKm2 < pTcathres) {
                    // we have a trigger point
                    outputIter.setSample(c, r, 0, triggerThreshold);
                }
            }
        }
        pm.worked(1);
    }
    pm.done();
    outTriggers = CoverageUtilities.buildCoverage("triggers", outputWR, regionMap, inElev.getCoordinateReferenceSystem());
}
Also used : GridCoverage2D(org.geotools.coverage.grid.GridCoverage2D) WritableRandomIter(javax.media.jai.iterator.WritableRandomIter) WritableRaster(java.awt.image.WritableRaster) RegionMap(org.hortonmachine.gears.utils.RegionMap) RandomIter(javax.media.jai.iterator.RandomIter) WritableRandomIter(javax.media.jai.iterator.WritableRandomIter) RenderedImage(java.awt.image.RenderedImage) OmsGradient(org.hortonmachine.hmachine.modules.geomorphology.gradient.OmsGradient) Execute(oms3.annotations.Execute)

Example 54 with WritableRandomIter

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

the class OmsInsolation method normalVector.

protected WritableRaster normalVector(WritableRaster pitWR, double res) {
    int minX = pitWR.getMinX();
    int minY = pitWR.getMinY();
    int rows = pitWR.getHeight();
    int cols = pitWR.getWidth();
    RandomIter pitIter = RandomIterFactory.create(pitWR, null);
    /*
         * Initializa the Image of the normal vector in the central point of the
         * cells, which have 3 components so the Image have 3 bands..
         */
    SampleModel sm = RasterFactory.createBandedSampleModel(5, cols, rows, 3);
    WritableRaster tmpNormalVectorWR = CoverageUtilities.createWritableRaster(cols, rows, null, sm, 0.0);
    WritableRandomIter tmpNormalIter = RandomIterFactory.createWritable(tmpNormalVectorWR, null);
    /*
         * apply the corripio's formula (is the formula (3) in the article)
         */
    for (int j = minY; j < minX + rows - 1; j++) {
        for (int i = minX; i < minX + cols - 1; i++) {
            double zij = pitIter.getSampleDouble(i, j, 0);
            double zidxj = pitIter.getSampleDouble(i + 1, j, 0);
            double zijdy = pitIter.getSampleDouble(i, j + 1, 0);
            double zidxjdy = pitIter.getSampleDouble(i + 1, j + 1, 0);
            double firstComponent = res * (zij - zidxj + zijdy - zidxjdy);
            double secondComponent = res * (zij + zidxj - zijdy - zidxjdy);
            double thirthComponent = 2 * (res * res);
            double den = Math.sqrt(firstComponent * firstComponent + secondComponent * secondComponent + thirthComponent * thirthComponent);
            tmpNormalIter.setPixel(i, j, new double[] { firstComponent / den, secondComponent / den, thirthComponent / den });
        }
    }
    pitIter.done();
    return tmpNormalVectorWR;
}
Also used : SampleModel(java.awt.image.SampleModel) WritableRandomIter(javax.media.jai.iterator.WritableRandomIter) WritableRaster(java.awt.image.WritableRaster) RandomIter(javax.media.jai.iterator.RandomIter) WritableRandomIter(javax.media.jai.iterator.WritableRandomIter) Point(org.locationtech.jts.geom.Point)

Example 55 with WritableRandomIter

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

the class OmsHazardClassifier method process.

@Execute
public void process() throws Exception {
    checkNull(inIntensityTr100, inIntensityTr200, inIntensityTr30);
    RegionMap regionMap = CoverageUtilities.getRegionParamsFromGridCoverage(inIntensityTr30);
    int nCols = regionMap.getCols();
    int nRows = regionMap.getRows();
    RandomIter tr200Iter = CoverageUtilities.getRandomIterator(inIntensityTr200);
    RandomIter tr100Iter = CoverageUtilities.getRandomIterator(inIntensityTr100);
    RandomIter tr30Iter = CoverageUtilities.getRandomIterator(inIntensityTr30);
    WritableRaster outIP1WR = CoverageUtilities.createWritableRaster(nCols, nRows, null, null, doubleNovalue);
    WritableRandomIter outIP1Iter = RandomIterFactory.createWritable(outIP1WR, null);
    WritableRaster outIP2WR = CoverageUtilities.createWritableRaster(nCols, nRows, null, null, doubleNovalue);
    WritableRandomIter outIP2Iter = RandomIterFactory.createWritable(outIP2WR, null);
    pm.beginTask("Processing map...", nRows);
    for (int r = 0; r < nRows; r++) {
        if (isCanceled(pm)) {
            return;
        }
        for (int c = 0; c < nCols; c++) {
            double tr30 = tr30Iter.getSampleDouble(c, r, 0);
            double tr100 = tr100Iter.getSampleDouble(c, r, 0);
            double tr200 = tr200Iter.getSampleDouble(c, r, 0);
            if (isNovalue(tr30) && isNovalue(tr100) && isNovalue(tr200)) {
                continue;
            }
            double tmpTr30;
            if (isNovalue(tr30)) {
                tmpTr30 = Double.NEGATIVE_INFINITY;
            } else if (dEq(tr30, INTENSITY_LOW)) {
                tmpTr30 = 3.0;
            } else if (dEq(tr30, INTENSITY_MEDIUM)) {
                tmpTr30 = 6.0;
            } else if (dEq(tr30, INTENSITY_HIGH)) {
                tmpTr30 = 9.0;
            } else {
                throw new ModelsIllegalargumentException("Unknown tr30 value: " + tr30, this, pm);
            }
            double tmpTr100;
            if (isNovalue(tr100)) {
                tmpTr100 = Double.NEGATIVE_INFINITY;
            } else if (dEq(tr100, INTENSITY_LOW)) {
                tmpTr100 = 2.0;
            } else if (dEq(tr100, INTENSITY_MEDIUM)) {
                tmpTr100 = 5.0;
            } else if (dEq(tr100, INTENSITY_HIGH)) {
                tmpTr100 = 8.0;
            } else {
                throw new ModelsIllegalargumentException("Unknown tr100 value: " + tr100, this, pm);
            }
            double tmpTr200;
            if (isNovalue(tr200)) {
                tmpTr200 = Double.NEGATIVE_INFINITY;
            } else if (dEq(tr200, INTENSITY_LOW)) {
                tmpTr200 = 1.0;
            } else if (dEq(tr200, INTENSITY_MEDIUM)) {
                tmpTr200 = 4.0;
            } else if (dEq(tr200, INTENSITY_HIGH)) {
                tmpTr200 = 7.0;
            } else {
                throw new ModelsIllegalargumentException("Unknown tr200 value: " + tr200, this, pm);
            }
            int maxValue = (int) max(tmpTr30, max(tmpTr100, tmpTr200));
            double[] reclassIP1 = { // 
            Double.NaN, // 1
            2, // 2
            3, // 3
            3, // 4
            3, // 5
            3, // 6
            4, // 7
            4, // 8
            4, // 9
            4 };
            double[] reclassIP2 = { // 
            Double.NaN, // 1
            2, // 2
            2, // 3
            3, // 4
            3, // 5
            3, // 6
            3, // 7
            4, // 8
            4, // 9
            4 };
            if (maxValue < 1 || maxValue > (reclassIP1.length - 1)) {
                throw new ModelsIllegalargumentException("Unknown max value from tr30/100/200: " + maxValue, this, pm);
            }
            double ip1 = reclassIP1[(int) maxValue];
            double ip2 = reclassIP2[(int) maxValue];
            outIP1Iter.setSample(c, r, 0, ip1);
            outIP2Iter.setSample(c, r, 0, ip2);
        }
        pm.worked(1);
    }
    pm.done();
    outHazardIP1 = CoverageUtilities.buildCoverage("ip1", outIP1WR, regionMap, inIntensityTr100.getCoordinateReferenceSystem());
    outHazardIP2 = CoverageUtilities.buildCoverage("ip2", outIP2WR, regionMap, inIntensityTr100.getCoordinateReferenceSystem());
}
Also used : ModelsIllegalargumentException(org.hortonmachine.gears.libs.exceptions.ModelsIllegalargumentException) WritableRandomIter(javax.media.jai.iterator.WritableRandomIter) WritableRaster(java.awt.image.WritableRaster) RegionMap(org.hortonmachine.gears.utils.RegionMap) RandomIter(javax.media.jai.iterator.RandomIter) WritableRandomIter(javax.media.jai.iterator.WritableRandomIter) Execute(oms3.annotations.Execute)

Aggregations

WritableRandomIter (javax.media.jai.iterator.WritableRandomIter)99 WritableRaster (java.awt.image.WritableRaster)85 RandomIter (javax.media.jai.iterator.RandomIter)60 RegionMap (org.hortonmachine.gears.utils.RegionMap)59 Execute (oms3.annotations.Execute)49 Coordinate (org.locationtech.jts.geom.Coordinate)20 GridCoverage2D (org.geotools.coverage.grid.GridCoverage2D)18 ModelsIllegalargumentException (org.hortonmachine.gears.libs.exceptions.ModelsIllegalargumentException)18 GridGeometry2D (org.geotools.coverage.grid.GridGeometry2D)17 RenderedImage (java.awt.image.RenderedImage)16 Point (java.awt.Point)11 ArrayList (java.util.ArrayList)11 Geometry (org.locationtech.jts.geom.Geometry)11 GridCoordinates2D (org.geotools.coverage.grid.GridCoordinates2D)10 FlowNode (org.hortonmachine.gears.libs.modules.FlowNode)10 DirectPosition2D (org.geotools.geometry.DirectPosition2D)9 CoordinateReferenceSystem (org.opengis.referencing.crs.CoordinateReferenceSystem)8 Envelope (org.locationtech.jts.geom.Envelope)7 DirectPosition (org.opengis.geometry.DirectPosition)7 SampleModel (java.awt.image.SampleModel)6