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());

    }
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 ) { qcritIter . setSample ( i , j , 0 , 5 ) ; } else { double checkStable = tanPhiValue * ( 1 - 1 / rhoValue ) + cohValue / ( 9810 * rhoValue * hsValue ) * ( 1 + pow ( slopeValue , 2 ) ) ; if ( slopeValue < checkStable ) { 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 ) ; 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 ( ) ; 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 ( ) ) ; }
Calculates the trasmissivity in every pixel of the map.