Search in sources :

Example 1 with StationsSelection

use of org.hortonmachine.hmachine.modules.statistics.kriging.utils.StationsSelection in project hortonmachine by TheHortonMachine.

the class OmsKrigingRasterMode method executeKriging.

/**
 * Executing ordinary kriging.
 * <p>
 * <li>Verify if the parameters are correct.
 * <li>Calculating the matrix of the covariance (a).
 * <li>For each point to interpolated, evalutate the know term vector (b)
 * and solve the system (a x)=b where x is the weight.
 * </p>
 *
 * @throws Exception the exception
 */
@Execute
public void executeKriging() throws Exception {
    inInterpolationGrid = inGridCoverage2D.getGridGeometry();
    verifyInput();
    demWR = mapsTransform(inGridCoverage2D);
    LinkedHashMap<Integer, Coordinate> pointsToInterpolateId2Coordinates = null;
    pointsToInterpolateId2Coordinates = getCoordinate(inInterpolationGrid);
    Set<Integer> pointsToInterpolateIdSet = pointsToInterpolateId2Coordinates.keySet();
    Iterator<Integer> idIterator = pointsToInterpolateIdSet.iterator();
    int j = 0;
    double[] result = new double[pointsToInterpolateId2Coordinates.size()];
    int[] idArray = new int[pointsToInterpolateId2Coordinates.size()];
    final DirectPosition gridPoint = new DirectPosition2D();
    MathTransform transf = inInterpolationGrid.getCRSToGrid2D();
    while (idIterator.hasNext()) {
        double sum = 0.;
        id = idIterator.next();
        idArray[j] = id;
        Coordinate coordinate = (Coordinate) pointsToInterpolateId2Coordinates.get(id);
        DirectPosition point = new DirectPosition2D(inInterpolationGrid.getCoordinateReferenceSystem(), coordinate.x, coordinate.y);
        transf.transform(point, gridPoint);
        double[] gridCoord = gridPoint.getCoordinate();
        int x = (int) gridCoord[0];
        int y = (int) gridCoord[1];
        /**
         * StationsSelection is an external class that allows the
         * selection of the stations involved in the study.
         * It is possible to define if to include stations with zero values,
         * station in a define neighborhood or within a max distance from
         * the considered point.
         */
        StationsSelection stations = new StationsSelection();
        stations.idx = coordinate.x;
        stations.idy = coordinate.y;
        stations.inStations = inStations;
        stations.inData = inData;
        stations.doIncludezero = doIncludezero;
        stations.maxdist = maxdist;
        stations.inNumCloserStations = inNumCloserStations;
        stations.fStationsid = fStationsid;
        stations.fStationsZ = fStationsZ;
        stations.execute();
        double[] xStations = stations.xStationInitialSet;
        double[] yStations = stations.yStationInitialSet;
        double[] zStations = stations.zStationInitialSet;
        double[] hStations = stations.hStationInitialSet;
        boolean areAllEquals = stations.areAllEquals;
        int n1 = xStations.length - 1;
        xStations[n1] = coordinate.x;
        yStations[n1] = coordinate.y;
        zStations[n1] = demWR.getSample(x, y, 0);
        double[] hresiduals = hStations;
        if (doDetrended == true) {
            if (zStations[n1] < 0) {
                doDetrended = false;
            } else {
                doDetrended = true;
            }
        }
        if (doDetrended) {
            RegressionLine t = new PolyTrendLine(regressionOrder);
            t.setValues(zStations, hStations);
            double[] regressionParameters = t.getRegressionParameters();
            trend_intercept = regressionParameters[0];
            trend_coefficient = regressionParameters[1];
            hresiduals = t.getResiduals();
        // Regression r = new Regression();
        // 
        // r = new Regression(zStations, hStations);
        // r.polynomial(regressionOrder);
        // 
        // /*If there is a trend for meteorological
        // * variables and elevation and it is statistically significant
        // * then the residuals from this linear trend
        // * are computed for each meteorological stations.
        // */
        // // if (Math.abs(r.getXYcorrCoeff()) > thresholdCorrelation) {
        // 
        // trend_intercept = r.getBestEstimates()[0];
        // trend_coefficient = r.getBestEstimates()[1];
        // hresiduals = r.getResiduals();
        // 
        // // } else {
        // // System.out.println("The trend is not significant");
        // // doDetrended=false;
        // // hresiduals=hStations;
        // 
        // // }
        }
        if (n1 != 0) {
            if (!areAllEquals && n1 > 1) {
                pm.beginTask(msg.message("kriging.working"), pointsToInterpolateId2Coordinates.size());
                double h0 = 0.0;
                /*
                     * calculating the covariance matrix.
                     */
                double[][] covarianceMatrix = covMatrixCalculating(xStations, yStations, zStations, n1);
                double[] knownTerm = knownTermsCalculation(xStations, yStations, zStations, n1);
                /*
                     * solve the linear system, where the result is the weight (moltiplicativeFactor).
                     */
                ColumnVector solution = SimpleLinearSystemSolverFactory.solve(knownTerm, covarianceMatrix, linearSystemSolverType);
                double[] moltiplicativeFactor = solution.copyValues1D();
                for (int k = 0; k < n1; k++) {
                    h0 = h0 + moltiplicativeFactor[k] * hresiduals[k];
                    // sum is computed to check that
                    // the sum of all the weights is 1
                    sum = sum + moltiplicativeFactor[k];
                }
                double trend = (doDetrended) ? zStations[n1] * trend_coefficient + trend_intercept : 0;
                h0 = h0 + trend;
                if (zStations[n1] < 0) {
                    result[j] = HMConstants.doubleNovalue;
                } else {
                    result[j] = h0;
                }
                j++;
                if (Math.abs(sum - 1) >= TOLL) {
                    throw new ModelsRuntimeException("Error in the coffeicients calculation", this.getClass().getSimpleName());
                }
                pm.worked(1);
            } else if (n1 == 1 || areAllEquals) {
                double tmp = hresiduals[0];
                pm.message(msg.message("kriging.setequalsvalue"));
                pm.beginTask(msg.message("kriging.working"), pointsToInterpolateId2Coordinates.size());
                if (zStations[n1] < 0) {
                    result[j] = HMConstants.doubleNovalue;
                } else {
                    result[j] = tmp;
                }
                j++;
                n1 = 0;
                pm.worked(1);
            }
            pm.done();
        } else {
            pm.errorMessage("No value for this time step");
            j = 0;
            double[] value = inData.values().iterator().next();
            if (zStations[n1] < 0) {
                result[j] = HMConstants.doubleNovalue;
            } else {
                result[j] = value[0];
            }
            j++;
        }
    }
    storeResult(result, pointsToInterpolateId2Coordinates);
}
Also used : DirectPosition(org.opengis.geometry.DirectPosition) RegressionLine(org.hortonmachine.gears.utils.math.regressions.RegressionLine) MathTransform(org.opengis.referencing.operation.MathTransform) DirectPosition2D(org.geotools.geometry.DirectPosition2D) ColumnVector(org.hortonmachine.gears.utils.math.matrixes.ColumnVector) PolyTrendLine(org.hortonmachine.gears.utils.math.regressions.PolyTrendLine) StationsSelection(org.hortonmachine.hmachine.modules.statistics.kriging.utils.StationsSelection) Coordinate(org.locationtech.jts.geom.Coordinate) ModelsRuntimeException(org.hortonmachine.gears.libs.exceptions.ModelsRuntimeException) Execute(oms3.annotations.Execute)

Example 2 with StationsSelection

use of org.hortonmachine.hmachine.modules.statistics.kriging.utils.StationsSelection in project hortonmachine by TheHortonMachine.

the class OmsKrigingVectorMode method executeKriging.

/**
 * Executing ordinary kriging.
 * <p>
 * <li>Verify if the parameters are correct.
 * <li>Calculating the matrix of the covariance (a).
 * <li>For each point to interpolated, evalutate the know term vector (b)
 * and solve the system (a x)=b where x is the weight.
 * </p>
 *
 * @throws Exception the exception
 */
@Execute
public void executeKriging() throws Exception {
    verifyInput();
    LinkedHashMap<Integer, Coordinate> pointsToInterpolateId2Coordinates = null;
    pointsToInterpolateId2Coordinates = getCoordinate(0, inInterpolate, fInterpolateid);
    Set<Integer> pointsToInterpolateIdSet = pointsToInterpolateId2Coordinates.keySet();
    Iterator<Integer> idIterator = pointsToInterpolateIdSet.iterator();
    int j = 0;
    double[] result = new double[pointsToInterpolateId2Coordinates.size()];
    int[] idArray = new int[pointsToInterpolateId2Coordinates.size()];
    while (idIterator.hasNext()) {
        double sum = 0.;
        id = idIterator.next();
        idArray[j] = id;
        Coordinate coordinate = (Coordinate) pointsToInterpolateId2Coordinates.get(id);
        /**
         * StationsSelection is an external class that allows the
         * selection of the stations involved in the study.
         * It is possible to define if to include stations with zero values,
         * station in a define neighborhood or within a max distance from
         * the considered point.
         */
        StationsSelection stations = new StationsSelection();
        stations.idx = coordinate.x;
        stations.idy = coordinate.y;
        stations.inStations = inStations;
        stations.inData = inData;
        stations.doIncludezero = doIncludezero;
        stations.maxdist = maxdist;
        stations.inNumCloserStations = inNumCloserStations;
        stations.fStationsid = fStationsid;
        stations.fStationsZ = fStationsZ;
        stations.execute();
        double[] xStations = stations.xStationInitialSet;
        double[] yStations = stations.yStationInitialSet;
        double[] zStations = stations.zStationInitialSet;
        double[] hStations = stations.hStationInitialSet;
        boolean areAllEquals = stations.areAllEquals;
        int n1 = xStations.length - 1;
        xStations[n1] = coordinate.x;
        yStations[n1] = coordinate.y;
        zStations[n1] = coordinate.getZ();
        double[] hresiduals = hStations;
        if (doDetrended) {
            RegressionLine t = new PolyTrendLine(regressionOrder);
            t.setValues(zStations, hStations);
            double[] regressionParameters = t.getRegressionParameters();
            trend_intercept = regressionParameters[0];
            trend_coefficient = regressionParameters[1];
            hresiduals = t.getResiduals();
        // Regression r = new Regression();
        // 
        // r = new Regression(zStations, hStations);
        // r.polynomial(regressionOrder);
        // 
        // /*If there is a trend for meteorological
        // * variables and elevation and it is statistically significant
        // * then the residuals from this linear trend
        // * are computed for each meteorological stations.
        // */
        // 
        // trend_intercept = r.getBestEstimates()[0];
        // trend_coefficient = r.getBestEstimates()[1];
        // hresiduals = r.getResiduals();
        }
        if (n1 != 0) {
            if (!areAllEquals && n1 > 1) {
                pm.beginTask(msg.message("kriging.working"), pointsToInterpolateId2Coordinates.size());
                double h0 = 0.0;
                /*
                     * calculating the covariance matrix.
                     */
                double[][] covarianceMatrix = covMatrixCalculating(xStations, yStations, zStations, n1);
                double[] knownTerm = knownTermsCalculation(xStations, yStations, zStations, n1);
                /*
                     * solve the linear system, where the result is the weight (moltiplicativeFactor).
                     */
                ColumnVector solution = SimpleLinearSystemSolverFactory.solve(knownTerm, covarianceMatrix, linearSystemSolverType);
                double[] moltiplicativeFactor = solution.copyValues1D();
                for (int k = 0; k < n1; k++) {
                    h0 = h0 + moltiplicativeFactor[k] * hresiduals[k];
                    // sum is computed to check that
                    // the sum of all the weights is 1
                    sum = sum + moltiplicativeFactor[k];
                }
                double trend = (doDetrended) ? coordinate.getZ() * trend_coefficient + trend_intercept : 0;
                h0 = h0 + trend;
                result[j] = h0;
                j++;
                if (Math.abs(sum - 1) >= TOLL) {
                    throw new ModelsRuntimeException("Error in the coffeicients calculation", this.getClass().getSimpleName());
                }
                pm.worked(1);
            } else if (n1 == 1 || areAllEquals) {
                double tmp = hresiduals[0];
                pm.message(msg.message("kriging.setequalsvalue"));
                pm.beginTask(msg.message("kriging.working"), pointsToInterpolateId2Coordinates.size());
                result[j] = tmp;
                j++;
                n1 = 0;
                pm.worked(1);
            }
            pm.done();
        } else {
            pm.errorMessage("No value for this time step");
            double[] value = inData.values().iterator().next();
            result[j] = value[0];
            j++;
        }
    }
    storeResult(result, idArray);
}
Also used : RegressionLine(org.hortonmachine.gears.utils.math.regressions.RegressionLine) ColumnVector(org.hortonmachine.gears.utils.math.matrixes.ColumnVector) PolyTrendLine(org.hortonmachine.gears.utils.math.regressions.PolyTrendLine) StationsSelection(org.hortonmachine.hmachine.modules.statistics.kriging.utils.StationsSelection) Coordinate(org.locationtech.jts.geom.Coordinate) ModelsRuntimeException(org.hortonmachine.gears.libs.exceptions.ModelsRuntimeException) Execute(oms3.annotations.Execute)

Example 3 with StationsSelection

use of org.hortonmachine.hmachine.modules.statistics.kriging.utils.StationsSelection in project hortonmachine by TheHortonMachine.

the class ExperimentalVariogram method process.

/**
 * Process.
 *
 * @throws Exception the exception
 */
@Execute
public void process() throws Exception {
    StationsSelection stations = new StationsSelection();
    stations.inStations = inStations;
    stations.inData = inData;
    stations.doIncludezero = doIncludezero;
    stations.fStationsid = fStationsid;
    stations.execute();
    differents = stations.n1;
    double[] xStations = stations.xStationInitialSet;
    double[] yStations = stations.yStationInitialSet;
    double[] hStations = stations.hStationInitialSet;
    int[] idStations = stations.idStationInitialSet;
    // number of different stations
    if (differents > 2) {
        double[][] outResult = processAlgorithm(xStations, yStations, hStations, Cutoffinput);
        storeResult(outResult);
    } else {
        System.out.println("Only 1 data >0 or All the data are equal. Variogram is not running");
    }
}
Also used : StationsSelection(org.hortonmachine.hmachine.modules.statistics.kriging.utils.StationsSelection) Execute(oms3.annotations.Execute)

Example 4 with StationsSelection

use of org.hortonmachine.hmachine.modules.statistics.kriging.utils.StationsSelection in project hortonmachine by TheHortonMachine.

the class OmsKrigingCheckMode method executeKriging.

/**
 * Executing ordinary kriging.
 * <p>
 * <li>Verify if the parameters are correct.
 * <li>Calculating the matrix of the covariance (a).
 * <li>For each point to interpolated, evalutate the know term vector (b)
 * and solve the system (a x)=b where x is the weight.
 * </p>
 *
 * @throws Exception the exception
 */
@Execute
public void executeKriging() throws Exception {
    verifyInput();
    LinkedHashMap<Integer, Coordinate> pointsToInterpolateId2Coordinates = null;
    pointsToInterpolateId2Coordinates = getCoordinate(0, inStations, fStationsid);
    Set<Integer> pointsToInterpolateIdSet = pointsToInterpolateId2Coordinates.keySet();
    Iterator<Integer> idIterator = pointsToInterpolateIdSet.iterator();
    int j = 0;
    double[] result = new double[pointsToInterpolateId2Coordinates.size()];
    int[] idArray = new int[pointsToInterpolateId2Coordinates.size()];
    while (idIterator.hasNext()) {
        double sum = 0.;
        id = idIterator.next();
        idArray[j] = id;
        Coordinate coordinate = (Coordinate) pointsToInterpolateId2Coordinates.get(id);
        /**
         * StationsSelection is an external class that allows the
         * selection of the stations involved in the study.
         * It is possible to define if to include stations with zero values,
         * station in a define neighborhood or within a max distance from
         * the considered point.
         */
        StationsSelection stations = new StationsSelection();
        stations.idx = coordinate.x;
        stations.idy = coordinate.y;
        stations.idOut = id;
        stations.inStations = inStations;
        stations.inData = inData;
        stations.doIncludezero = doIncludezero;
        stations.maxdist = maxdist;
        stations.inNumCloserStations = inNumCloserStations;
        stations.fStationsid = fStationsid;
        stations.fStationsZ = fStationsZ;
        stations.execute();
        double[] xStations = stations.xStationInitialSet;
        double[] yStations = stations.yStationInitialSet;
        double[] zStations = stations.zStationInitialSet;
        double[] hStations = stations.hStationInitialSet;
        boolean areAllEquals = stations.areAllEquals;
        int n1 = xStations.length - 1;
        xStations[n1] = coordinate.x;
        yStations[n1] = coordinate.y;
        zStations[n1] = coordinate.z;
        double[] hresiduals = hStations;
        if (doDetrended) {
            RegressionLine t = new PolyTrendLine(regressionOrder);
            t.setValues(zStations, hStations);
            double[] regressionParameters = t.getRegressionParameters();
            trend_intercept = regressionParameters[0];
            trend_coefficient = regressionParameters[1];
            hresiduals = t.getResiduals();
        // Regression r = new Regression();
        // 
        // r = new Regression(zStations, hStations);
        // r.polynomial(regressionOrder);
        // 
        // /*If there is a trend for meteorological
        // * variables and elevation and it is statistically significant
        // * then the residuals from this linear trend
        // * are computed for each meteorological stations.
        // */
        // 
        // trend_intercept = r.getBestEstimates()[0];
        // trend_coefficient = r.getBestEstimates()[1];
        // hresiduals = r.getResiduals();
        }
        if (n1 != 0) {
            if (!areAllEquals && n1 > 1) {
                pm.beginTask(msg.message("kriging.working"), pointsToInterpolateId2Coordinates.size());
                double h0 = 0.0;
                /*
                     * calculating the covariance matrix.
                     */
                double[][] covarianceMatrix = covMatrixCalculating(xStations, yStations, zStations, n1);
                double[] knownTerm = knownTermsCalculation(xStations, yStations, zStations, n1);
                /*
                     * solve the linear system, where the result is the weight (moltiplicativeFactor).
                     */
                ColumnVector solution = SimpleLinearSystemSolverFactory.solve(knownTerm, covarianceMatrix, linearSystemSolverType);
                double[] moltiplicativeFactor = solution.copyValues1D();
                for (int k = 0; k < n1; k++) {
                    h0 = h0 + moltiplicativeFactor[k] * hresiduals[k];
                    // sum is computed to check that
                    // the sum of all the weights is 1
                    sum = sum + moltiplicativeFactor[k];
                }
                double trend = (doDetrended) ? coordinate.z * trend_coefficient + trend_intercept : 0;
                h0 = h0 + trend;
                result[j] = h0;
                j++;
                if (Math.abs(sum - 1) >= TOLL) {
                    throw new ModelsRuntimeException("Error in the coffeicients calculation", this.getClass().getSimpleName());
                }
                pm.worked(1);
            } else if (n1 == 1 || areAllEquals) {
                double tmp = hresiduals[0];
                pm.message(msg.message("kriging.setequalsvalue"));
                pm.beginTask(msg.message("kriging.working"), pointsToInterpolateId2Coordinates.size());
                result[j] = tmp;
                j++;
                n1 = 0;
                pm.worked(1);
            }
            pm.done();
        } else {
            pm.errorMessage("No value for this time step");
            j = 0;
            double[] value = inData.values().iterator().next();
            result[j] = value[0];
            j++;
        }
    }
    storeResult(result, idArray);
}
Also used : RegressionLine(org.hortonmachine.gears.utils.math.regressions.RegressionLine) ColumnVector(org.hortonmachine.gears.utils.math.matrixes.ColumnVector) PolyTrendLine(org.hortonmachine.gears.utils.math.regressions.PolyTrendLine) StationsSelection(org.hortonmachine.hmachine.modules.statistics.kriging.utils.StationsSelection) Coordinate(org.locationtech.jts.geom.Coordinate) ModelsRuntimeException(org.hortonmachine.gears.libs.exceptions.ModelsRuntimeException) Execute(oms3.annotations.Execute)

Aggregations

Execute (oms3.annotations.Execute)4 StationsSelection (org.hortonmachine.hmachine.modules.statistics.kriging.utils.StationsSelection)4 ModelsRuntimeException (org.hortonmachine.gears.libs.exceptions.ModelsRuntimeException)3 ColumnVector (org.hortonmachine.gears.utils.math.matrixes.ColumnVector)3 PolyTrendLine (org.hortonmachine.gears.utils.math.regressions.PolyTrendLine)3 RegressionLine (org.hortonmachine.gears.utils.math.regressions.RegressionLine)3 Coordinate (org.locationtech.jts.geom.Coordinate)3 DirectPosition2D (org.geotools.geometry.DirectPosition2D)1 DirectPosition (org.opengis.geometry.DirectPosition)1 MathTransform (org.opengis.referencing.operation.MathTransform)1