Examples with Range oms3.annotations.Range used on opensource projects

Example 1 with Range

use of oms3.annotations.Range in project hortonmachine by TheHortonMachine.

the class OmsSigmaFilterPlus method process.

@Execute
public void process() throws Exception {
    checkNull(inGeodata);
    RegionMap regionMap = CoverageUtilities.getRegionParamsFromGridCoverage(inGeodata);
    int cols = regionMap.getCols();
    int rows = regionMap.getRows();
    double[] pixels = CoverageUtilities.renderedImage2DoubleArray(inGeodata.getRenderedImage(), 3);
    /* 
         *  Create a circular kernel of a given radius. Radius = 0.5 includes the 4 neighbors of the
         *  pixel in the center, radius = 1 corresponds to a 3x3 kernel size.
         *  The output is written to class variables kNPoints (number of points inside the kernel) and
         *  lineRadius, which is an array giving the radius of each line. Line length is 2*lineRadius+1.
         */
    if (// this code creates the same sizes as the previous
    pRadius >= 1.5 && pRadius < 1.75)
        // RankFilters
        pRadius = 1.75;
    else if (pRadius >= 2.5 && pRadius < 2.85)
        pRadius = 2.85;
    int r2 = (int) (pRadius * pRadius) + 1;
    kRadius = (int) (Math.sqrt(r2 + 1e-10));
    lineRadius = new int[2 * kRadius + 1];
    lineRadius[kRadius] = kRadius;
    kNPoints = 2 * kRadius + 1;
    for (int y = 1; y <= kRadius; y++) {
        int dx = (int) (Math.sqrt(r2 - y * y + 1e-10));
        lineRadius[kRadius + y] = dx;
        lineRadius[kRadius - y] = dx;
        kNPoints += 4 * dx + 2;
    }
    // process passes
    for (int pass = 0; pass < nPasses; pass++) {
        // min pixels in sigma
        int minPixNumber = (int) (kNPoints * minPixFraction + 0.999999);
        // range
        pixels = doFiltering(pixels, cols, rows, kRadius, lineRadius, pSigmaWidth, minPixNumber, outlierAware);
    }
    WritableRaster outWR = CoverageUtilities.doubleArray2WritableRaster(pixels, cols, rows);
    outGeodata = CoverageUtilities.buildCoverage("sigma", outWR, regionMap, inGeodata.getCoordinateReferenceSystem());
}

Example 2 with Range

use of oms3.annotations.Range in project hortonmachine by TheHortonMachine.

the class LasHeightDistribution method process.

@Execute
public void process() throws Exception {
    checkNull(inIndexFile, inVector, inDem);
    // if (outChartsFolder != null) {
    // outChartsFolderFile = new File(outChartsFolder);
    // if (outChartsFolderFile.exists()) {
    // doChart = true;
    // }
    // }
    double percentageOverlap = pOverlapPerc / 100.0;
    File indexFile = new File(inIndexFile);
    SimpleFeatureCollection tilesFC = getVector(inVector);
    List<FeatureMate> tilesMates = FeatureUtilities.featureCollectionToMatesList(tilesFC);
    GridCoverage2D inDemGC = getRaster(inDem);
    CoordinateReferenceSystem crs = inDemGC.getCoordinateReferenceSystem();
    WritableRaster[] finalCoverageWRH = new WritableRaster[1];
    GridCoverage2D outCatsGC = CoverageUtilities.createCoverageFromTemplate(inDemGC, HMConstants.doubleNovalue, finalCoverageWRH);
    WritableRandomIter finalIter = CoverageUtilities.getWritableRandomIterator(finalCoverageWRH[0]);
    try (ALasDataManager dataManager = ALasDataManager.getDataManager(indexFile, inDemGC, pThres, null)) {
        dataManager.open();
        for (int i = 0; i < tilesMates.size(); i++) {
            pm.message("Processing tile: " + i + "/" + tilesMates.size());
            FeatureMate tileMate = tilesMates.get(i);
            String id = tileMate.getAttribute(fId, String.class);
            Geometry tileGeom = tileMate.getGeometry();
            Envelope geomEnvelope = tileGeom.getEnvelopeInternal();
            ReferencedEnvelope refEnvelope = new ReferencedEnvelope(geomEnvelope, crs);
            Envelope2D tileEnvelope = new Envelope2D(refEnvelope);
            WritableRaster[] tmpWrH = new WritableRaster[1];
            GridCoverage2D tmp = CoverageUtilities.createSubCoverageFromTemplate(inDemGC, tileEnvelope, doubleNovalue, tmpWrH);
            RegionMap tileRegionMap = CoverageUtilities.getRegionParamsFromGridCoverage(tmp);
            GridGeometry2D tileGridGeometry = tmp.getGridGeometry();
            List<LasRecord> pointsListForTile = dataManager.getPointsInGeometry(tileGeom, true);
            if (pointsListForTile.size() == 0) {
                pm.errorMessage("No points found in tile: " + id);
                continue;
            }
            if (pointsListForTile.size() < 2) {
                pm.errorMessage("Not enough points found in tile: " + id);
                continue;
            }
            List<double[]> negativeRanges = analyseNegativeLayerRanges(id, pointsListForTile);
            List<GridCoverage2D> rangeCoverages = new ArrayList<GridCoverage2D>();
            for (double[] range : negativeRanges) {
                List<LasRecord> pointsInVerticalRange = ALasDataManager.getPointsInVerticalRange(pointsListForTile, range[0], range[1], true);
                WritableRaster[] wrH = new WritableRaster[1];
                GridCoverage2D tmpCoverage = CoverageUtilities.createSubCoverageFromTemplate(inDemGC, tileEnvelope, doubleNovalue, wrH);
                rangeCoverages.add(tmpCoverage);
                WritableRandomIter tmpIter = CoverageUtilities.getWritableRandomIterator(wrH[0]);
                final DirectPosition2D wp = new DirectPosition2D();
                for (LasRecord lasRecord : pointsInVerticalRange) {
                    wp.setLocation(lasRecord.x, lasRecord.y);
                    GridCoordinates2D gp = tileGridGeometry.worldToGrid(wp);
                    double count = tmpIter.getSampleDouble(gp.x, gp.y, 0);
                    if (isNovalue(count)) {
                        count = 0;
                    }
                    tmpIter.setSample(gp.x, gp.y, 0, count + 1);
                }
                tmpIter.done();
            }
            /*
                 * categorize in non forest/single/double layer
                 */
            /*
                 * 1) check if the multiple layers are double or
                 * single at variable heights. 
                 */
            boolean isDoubleLayered = false;
            if (rangeCoverages.size() > 1) {
                for (int j = 0; j < rangeCoverages.size() - 1; j++) {
                    GridCoverage2D cov1 = rangeCoverages.get(j);
                    GridCoverage2D cov2 = rangeCoverages.get(j + 1);
                    if (overlapForPercentage(cov1, cov2, percentageOverlap)) {
                        isDoubleLayered = true;
                        break;
                    }
                }
            }
            /*
                 * 2) define each pixel of the end map
                 * - 0 = no forest
                 * - 1 = single layer
                 * - 2 = single with variable height
                 * - 3 = double layer
                 */
            GridGeometry2D gridGeometry = outCatsGC.getGridGeometry();
            // RegionMap regionMap = CoverageUtilities.getRegionParamsFromGridCoverage(outCats);
            double[] gridValue = new double[1];
            int cols = tileRegionMap.getCols();
            int rows = tileRegionMap.getRows();
            for (int r = 0; r < rows; r++) {
                for (int c = 0; c < cols; c++) {
                    int value = 0;
                    GridCoordinates2D gridPosition = new GridCoordinates2D(c, r);
                    for (int j = 0; j < rangeCoverages.size(); j++) {
                        GridCoverage2D cov = rangeCoverages.get(j);
                        cov.evaluate(gridPosition, gridValue);
                        if (!isNovalue(gridValue[0])) {
                            value++;
                        }
                    }
                    // set final value in the grid
                    if (value > 1) {
                        // multilayer
                        if (isDoubleLayered) {
                            value = 3;
                        } else {
                            value = 2;
                        }
                    }
                    DirectPosition worldPosition = tileGridGeometry.gridToWorld(gridPosition);
                    GridCoordinates2D worldPositionCats = gridGeometry.worldToGrid(worldPosition);
                    finalIter.setSample(worldPositionCats.x, worldPositionCats.y, 0, value);
                }
            }
        }
    }
    RegionMap regionMap = CoverageUtilities.getRegionParamsFromGridCoverage(outCatsGC);
    int cols = regionMap.getCols();
    int rows = regionMap.getRows();
    for (int r = 0; r < rows; r++) {
        for (int c = 0; c < cols; c++) {
            double value = finalIter.getSampleDouble(c, r, 0);
            if (isNovalue(value)) {
                finalIter.setSample(c, r, 0, 0.0);
            }
        }
    }
    finalIter.done();
    dumpRaster(outCatsGC, outCats);
}

Example 3 with Range

use of oms3.annotations.Range in project hortonmachine by TheHortonMachine.

the class OmsJami method process.

@Execute
public void process() throws Exception {
    pm.message("processing " + tCurrent + " " + pType);
    checkNull(inAltimetry, inAreas, inMeteo, inStations);
    // check fascie num
    int fascieNum = 0;
    for (EIAreas area : inAreas) {
        int idAltimetricBand = area.altimetricBandId;
        if (idAltimetricBand > fascieNum) {
            fascieNum = idAltimetricBand;
        }
    }
    // increment by one. Fascie start from 0, so 0-5 are 6 fascie
    fascieNum++;
    currentTimestamp = formatter.parseDateTime(tCurrent);
    outInterpolatedBand = new HashMap<Integer, double[]>();
    outInterpolated = new HashMap<Integer, double[]>();
    /*
         * get stations
         */
    pm.message("Read stations data.");
    stationCoordinates = new ArrayList<Coordinate>();
    stationFeatures = new ArrayList<SimpleFeature>();
    FeatureIterator<SimpleFeature> featureIterator = inStations.features();
    while (featureIterator.hasNext()) {
        SimpleFeature feature = featureIterator.next();
        Coordinate stationCoord = ((Geometry) feature.getDefaultGeometry()).getCoordinate();
        stationCoordinates.add(stationCoord);
        stationFeatures.add(feature);
    }
    featureIterator.close();
    pm.done();
    /*
         * get basins and create of every basin a buffered geometry at 10,
         * 20 and 50 km, which will be used to find the nearest stations
         * around.
         */
    basinBaricenterCoordinates = new ArrayList<Coordinate>();
    basinFeatures = new ArrayList<SimpleFeature>();
    featureIterator = inInterpolate.features();
    pm.beginTask("Read basins data.", inInterpolate.size());
    while (featureIterator.hasNext()) {
        pm.worked(1);
        SimpleFeature feature = featureIterator.next();
        Coordinate baricenterCoord = ((Geometry) feature.getDefaultGeometry()).getCentroid().getCoordinate();
        basinBaricenterCoordinates.add(baricenterCoord);
        basinFeatures.add(feature);
    }
    featureIterator.close();
    pm.done();
    statElev = new double[stationCoordinates.size()];
    statId = new double[stationCoordinates.size()];
    stationId2CoordinateMap = new HashMap<Integer, Coordinate>();
    extractFromStationFeatures();
    stationid2StationindexMap = new HashMap<Integer, Integer>();
    for (int i = 0; i < statId.length; i++) {
        stationid2StationindexMap.put((int) statId[i], i);
    }
    /*
         * get the basin's id attribute index
         */
    SimpleFeature tmpfeature = basinFeatures.get(0);
    SimpleFeatureType featureType = tmpfeature.getFeatureType();
    basinIdFieldIndex = featureType.indexOf(fBasinid);
    if (basinIdFieldIndex == -1) {
        throw new IllegalArgumentException("The field of the basin id couldn't be found in the supplied basin data.");
    }
    basinid2BasinindexMap = new HashMap<Integer, Integer>();
    basinindex2basinidMap = new HashMap<Integer, Integer>();
    for (int i = 0; i < basinBaricenterCoordinates.size(); i++) {
        int basinid = ((Number) basinFeatures.get(i).getAttribute(basinIdFieldIndex)).intValue();
        basinid2BasinindexMap.put(basinid, i);
        basinindex2basinidMap.put(i, basinid);
    }
    calculateAreas(fascieNum);
    pm.message("Creating the band's elevation for every basin matrix.");
    /*
         * create the altimetric bands matrix
         */
    int basinNum = basinBaricenterCoordinates.size();
    bandsBasins = new double[fascieNum][basinNum];
    for (int i = 0; i < inAltimetry.size(); i++) {
        EIAltimetry eiAltimetry = inAltimetry.get(i);
        int idbasin = eiAltimetry.basinId;
        int idfascia = eiAltimetry.altimetricBandId;
        double elevationInBandBaricenter = eiAltimetry.elevationValue;
        Integer index = basinid2BasinindexMap.get(idbasin);
        if (index != null)
            bandsBasins[idfascia][index] = elevationInBandBaricenter;
    // TODO make it range aware
    // double range = altimetryScalarSet.get(i + 3);
    // bandsBasins[idfascia +
    // 1][basinid2BasinindexMap.get(idbasin)] = baricenter
    // + range / 2.0;
    }
    double[] stationBinsArrays = new double[pBins + 1];
    double maxStatElev = statElev[statElev.length - 1];
    double minStatElev = statElev[0];
    double interval = (maxStatElev - minStatElev) / (double) pBins;
    for (int i = 0; i < stationBinsArrays.length; i++) {
        stationBinsArrays[i] = minStatElev + i * interval;
    }
    /*
         * find all stations inside a elevation band for every basin
         */
    pm.beginTask("Finding all stations inside a elevation band for every basin.", stationBinsArrays.length - 1);
    bin2StationsListMap = new HashMap<Integer, List<Integer>>();
    for (int i = 0; i < stationBinsArrays.length - 1; i++) {
        List<Integer> stationsIds = new ArrayList<Integer>();
        for (int j = 0; j < statId.length; j++) {
            double id = statId[j];
            double elev = statElev[j];
            if (elev >= stationBinsArrays[i] && elev < stationBinsArrays[i + 1]) {
                stationsIds.add((int) id);
            }
        }
        bin2StationsListMap.put(i, stationsIds);
        pm.worked(1);
    }
    pm.done();
    /*
         * get values for current timestep and order them with the stations ids
         */
    double[] statValues = new double[stationCoordinates.size()];
    for (int i = 0; i < statValues.length; i++) {
        statValues[i] = doubleNovalue;
    }
    Set<Integer> stationIdSet = inMeteo.keySet();
    for (Integer stationId : stationIdSet) {
        int id = stationId;
        double[] value = inMeteo.get(id);
        Integer index = stationid2StationindexMap.get((int) id);
        if (index == null)
            continue;
        statValues[index] = value[0];
    }
    // number of active stations for every basin
    int[] activeStationsPerBasin = new int[basinBaricenterCoordinates.size()];
    int[][] stazBacMatrix = createStationBasinsMatrix(statValues, activeStationsPerBasin);
    int[][] stations = new int[stazBacMatrix.length][stazBacMatrix[0].length];
    int contStations = 0;
    // riempimento vettori/matrici
    for (int i = 0; i < stazBacMatrix[0].length; i++) {
        // indice bacino
        contStations = 0;
        for (int j = 0; j < stazBacMatrix.length; j++) {
            // indice stazione
            if (stazBacMatrix[j][i] == 1) {
                stations[contStations][i] = j;
                contStations += 1;
            }
        }
    }
    int bandsNum = bandsBasins.length;
    if (pType == DTDAY || pType == DTMONTH) {
        /*
             * calculate the DT month and day for each station
             */
        // System.out.println("Calculating the dayly and monthly Dt for each station...");
        rangeT(statValues);
    }
    pm.beginTask("Interpolating over bands and basins...", basinBaricenterCoordinates.size());
    for (int i = 0; i < basinBaricenterCoordinates.size(); i++) {
        pm.worked(1);
        // interpolated value for every band
        double[] interpolatedMeteoForBand = new double[bandsNum];
        double interpolatedMeteoForBasin = 0;
        int cont = 0;
        double h;
        int[] jj_av;
        // trova le stazioni che forniscono dati
        // costruisco un nuovo
        jj_av = new int[activeStationsPerBasin[i]];
        // stazioni del bacino in studio
        for (int j = 0; j < activeStationsPerBasin[i]; j++) {
            if (pType != DTDAY && pType != DTMONTH) {
                if (!isNovalue(statValues[stations[j][i]])) {
                    // registro le stazioni
                    jj_av[cont] = stations[j][i];
                    // attive
                    cont += 1;
                }
            } else {
                // NODATA
                if (!isNovalue(minTempPerStation[stations[j][i]]) && isNovalue(maxTempPerStation[stations[j][i]])) {
                    // jj conterrà le stazioni che hanno dati di escursione
                    // termica
                    // giornaliera
                    // registro le stazioni
                    jj_av[cont] = stations[j][i];
                    // attive
                    cont += 1;
                }
            }
        }
        // sopravviva
        if (cont == 0) {
            if (pType == TEMPERATURE) {
                // caso dei dati di temperatura
                pm.errorMessage("ERRORE: PER IL BACINO " + i + " NON SONO DISPONIBILI DATI DI TEMPERATURA, PER QUESTO BACINO STAND-BY");
                for (int f = 0; f < bandsNum; f++) {
                    // per tutte le fasce
                    // altimetriche metto il
                    // dato a -100
                    interpolatedMeteoForBand[f] = doubleNovalue;
                }
                interpolatedMeteoForBasin = doubleNovalue;
            } else if (pType == PRESSURE) {
                // caso dei dati di pressione
                pm.message("  -> Per il bacino " + i + " non sono disponibili dati di pressione, uso valori di default");
                for (int f = 0; f < bandsNum; f++) {
                    // per tutte le fasce
                    // altimetriche considero
                    // un'adiabatica
                    interpolatedMeteoForBand[f] = 1013.25 * Math.exp(-(bandsBasins[f][i]) * 0.00013);
                    interpolatedMeteoForBasin = interpolatedMeteoForBasin + interpolatedMeteoForBand[f] * basinAreasPerFascias[i][f] / basinAreas[i];
                }
            } else if (pType == HUMIDITY) {
                // caso dei dati di umidità
                pm.message("  -> Per il bacino " + i + " non sono disponibili dati di umidita', uso valori di default");
                for (int f = 0; f < bandsNum; f++) {
                    // per tutte le fasce
                    // altimetriche metto NODATA
                    interpolatedMeteoForBand[f] = defaultRh;
                }
                interpolatedMeteoForBasin = defaultRh;
            } else if (pType == WIND) {
                // caso dei dati di velocità del vento
                pm.message("  -> Per il bacino " + i + " non sono disponibili dati di velocita' del vento, uso valori di default");
                for (int f = 0; f < bandsNum; f++) {
                    // per tutte le fasce
                    // altimetriche metto NODATA
                    interpolatedMeteoForBand[f] = defaultW;
                }
                interpolatedMeteoForBasin = defaultW;
            } else if (pType == DTDAY) {
                // caso dei dati di escursione termica
                // giornaliera
                pm.message("  -> Per il bacino " + i + " non sono disponibili dati di escursione termica giornaliera', uso valori di default");
                for (int f = 0; f < bandsNum; f++) {
                    // per tutte le fasce
                    // altimetriche del bacino
                    // assegno all'escursione termica giornaliera il dato
                    // DTd
                    // messo nel file dei parametri
                    interpolatedMeteoForBand[f] = defaultDtday;
                }
                interpolatedMeteoForBasin = defaultDtday;
            } else if (pType == DTMONTH) {
                // caso dei dati di escursione termica
                // mensile
                pm.message("  -> Per il bacino " + i + " non sono disponibili dati di escursione termica mensile', uso valori di default");
                for (int f = 0; f < bandsNum; f++) {
                    /*
                         *  per tutte le fasce
                         * altimetriche del bacino
                         */
                    // assegno all'escursione termica media mensile il
                    // datoDTm
                    // messo nel file dei parametri
                    interpolatedMeteoForBand[f] = defaultDtmonth;
                }
                interpolatedMeteoForBasin = defaultDtmonth;
            }
        } else if (cont == 1) {
            // standard per T e P, valori costanti per RH e V
            for (int f = 0; f < bandsNum; f++) {
                // altimetriche
                if (pType == TEMPERATURE) {
                    // trasformo la temp in K e calcolo T
                    // con
                    // l'adiabatica semplice
                    interpolatedMeteoForBand[f] = (statValues[jj_av[0]] + tk) * Math.exp(-(bandsBasins[f][i] - statElev[jj_av[0]]) * GAMMA / (statValues[jj_av[0]] + tk)) - tk;
                } else if (pType == PRESSURE) {
                    // calcolo P con il gradiente
                    // adiabatico
                    interpolatedMeteoForBand[f] = statValues[jj_av[0]] * Math.exp(-(bandsBasins[f][i] - statElev[jj_av[0]]) * 0.00013);
                } else if (pType == DTDAY) {
                    // se ho una sola stazione assegno il valore della
                    // stazione a tutto il
                    // bacino
                    // altimetriche del bacino assegno il valore di
                    // escursione massima
                    // giornaliera
                    interpolatedMeteoForBand[f] = maxTempPerStation[jj_av[0]] - minTempPerStation[jj_av[0]];
                    if ((maxTempPerStation[jj_av[0]] - minTempPerStation[jj_av[0]]) <= 0) {
                        interpolatedMeteoForBand[f] = defaultDtday;
                    }
                } else if (pType == DTMONTH) {
                    // se ho una sola stazione assegno il valore della
                    // stazione a tutto il
                    // bacino
                    // altimetriche del bacino assegno il valore di
                    // escursione massima mensile
                    interpolatedMeteoForBand[f] = DTmonth[jj_av[0]];
                } else {
                    // RH e V sono costanti al variare delle fasce
                    // altimetriche
                    interpolatedMeteoForBand[f] = statValues[jj_av[0]];
                }
                interpolatedMeteoForBasin = interpolatedMeteoForBasin + interpolatedMeteoForBand[f] * basinAreasPerFascias[i][f] / basinAreas[i];
            }
        } else {
            // caso 2. ci sono almeno 2 stazioni (a quote inferiori alla
            // stazioni piu' bassa considero atmosfera standard come a quote
            // superiori alla staz. piu' alta, in mezzo calcolo LAPSE RATE)
            // alloca L (vettore di dimensioni numero di stazioni attive-1)
            double[] lapseRate = new double[cont - 1];
            for (int j = 0; j < cont - 1; j++) {
                // le stazioni sono in
                // ordine di
                // quota
                // L[j] e' il lapse rate tra la stazione j e j+1, puo'
                // essere
                // calcolato dai dati per j che va da 1 a n-1, dove n e' il
                // numero di stazioni (cont)
                lapseRate[j] = (statValues[jj_av[j]] - statValues[jj_av[j + 1]]) / (statElev[jj_av[j + 1]] - statElev[jj_av[j]]);
            }
            for (int f = 0; f < bandsNum; f++) {
                // più bassa
                if (bandsBasins[f][i] <= statElev[jj_av[0]]) {
                    if (pType == TEMPERATURE) {
                        // T
                        interpolatedMeteoForBand[f] = statValues[jj_av[0]] - GAMMA * (bandsBasins[f][i] - statElev[jj_av[0]]);
                    } else if (pType == PRESSURE) {
                        // P
                        interpolatedMeteoForBand[f] = statValues[jj_av[0]] - (statValues[jj_av[0]] * 0.00013) * (bandsBasins[f][i] - statElev[jj_av[0]]);
                    } else if (pType == DTDAY) {
                        interpolatedMeteoForBand[f] = maxTempPerStation[jj_av[0]] - minTempPerStation[jj_av[0]];
                        if ((maxTempPerStation[jj_av[0]] - minTempPerStation[jj_av[0]]) <= 0) {
                            interpolatedMeteoForBand[f] = defaultDtday;
                        }
                    } else if (pType == DTMONTH) {
                        interpolatedMeteoForBand[f] = DTmonth[jj_av[0]];
                    } else {
                        // RH e V
                        interpolatedMeteoForBand[f] = statValues[jj_av[0]];
                    }
                // per le fasce altimetriche con quote piu' alte della
                // quota
                // della stazione piu' alta prendo i dati della stazione
                // più alta
                } else if (bandsBasins[f][i] >= statElev[jj_av[cont - 1]]) {
                    if (pType == TEMPERATURE) {
                        // T
                        interpolatedMeteoForBand[f] = statValues[jj_av[cont - 1]] - GAMMA * (bandsBasins[f][i] - statElev[jj_av[cont - 1]]);
                    } else if (pType == PRESSURE) {
                        // P
                        interpolatedMeteoForBand[f] = statValues[jj_av[cont - 1]] - (statValues[jj_av[cont - 1]] * 0.00013) * (bandsBasins[f][i] - statElev[jj_av[cont - 1]]);
                    } else if (pType == DTDAY) {
                        interpolatedMeteoForBand[f] = maxTempPerStation[jj_av[cont - 1]] - minTempPerStation[jj_av[cont - 1]];
                        if ((maxTempPerStation[jj_av[0]] - minTempPerStation[jj_av[0]]) <= 0) {
                            interpolatedMeteoForBand[f] = defaultDtday;
                        }
                    } else if (pType == DTMONTH) {
                        interpolatedMeteoForBand[f] = DTmonth[jj_av[cont - 1]];
                    } else {
                        // RH e V
                        interpolatedMeteoForBand[f] = statValues[jj_av[cont - 1]];
                    }
                } else {
                    int k = cont - 1;
                    if (pType == DTDAY) {
                        // la max delle stazioni
                        do {
                            k -= 1;
                            h = statElev[jj_av[k]];
                        } while (bandsBasins[f][i] <= h);
                        // for (int j = 0; j < cont; j++) {
                        // if (f ==0 && i == 100) {
                        // System.out.println(j + " "+ statElev[jj_av[j]]);
                        // }
                        // }
                        // interpolatedMeteoForBand[f] =
                        // ((maxTempPerStation[jj_av[k]] -
                        // minTempPerStation[jj_av[k]])
                        // * (statElev[jj_av[k + 1]] - bandsBasins[f][i]) +
                        // (maxTempPerStation[jj_av[k + 1]] -
                        // minTempPerStation[jj_av[k + 1]])
                        // * (bandsBasins[f][i] - statElev[jj_av[k]]))
                        // / (statElev[jj_av[k + 1]] - statElev[jj_av[k]]);
                        interpolatedMeteoForBand[f] = ((maxTempPerStation[jj_av[k + 1]] - minTempPerStation[jj_av[k + 1]]) - (maxTempPerStation[jj_av[k]] - minTempPerStation[jj_av[k]])) * (bandsBasins[f][i] - statElev[jj_av[k]]) / (statElev[jj_av[k + 1]] - statElev[jj_av[k]]) + (maxTempPerStation[jj_av[k]] - minTempPerStation[jj_av[k]]);
                        if (interpolatedMeteoForBand[f] <= 0) {
                            interpolatedMeteoForBand[f] = defaultDtday;
                        }
                    } else if (pType == DTMONTH) {
                        // la max delle stazioni
                        do {
                            k -= 1;
                            h = statElev[jj_av[k]];
                        } while (bandsBasins[f][i] <= h);
                        interpolatedMeteoForBand[f] = (DTmonth[jj_av[k]] * (statElev[jj_av[k + 1]] - bandsBasins[f][i]) + DTmonth[jj_av[k + 1]] * (bandsBasins[f][i] - statElev[jj_av[k]])) / (statElev[jj_av[k + 1]] - statElev[jj_av[k]]);
                    } else {
                        do {
                            k -= 1;
                            h = statElev[jj_av[k]];
                        } while (bandsBasins[f][i] <= h);
                        interpolatedMeteoForBand[f] = statValues[jj_av[k]] - lapseRate[k] * (bandsBasins[f][i] - statElev[jj_av[k]]);
                    }
                }
                interpolatedMeteoForBasin = interpolatedMeteoForBasin + interpolatedMeteoForBand[f] * basinAreasPerFascias[i][f] / basinAreas[i];
            }
            // controllo su RH>100 e v=0
            if (pType == HUMIDITY) {
                // RH
                double MAX_HUMIDITY = 100;
                double MIN_HUMIDITY = 5;
                for (int f = 0; f < bandsNum; f++) {
                    if (interpolatedMeteoForBand[f] > MAX_HUMIDITY)
                        interpolatedMeteoForBand[f] = MAX_HUMIDITY;
                    if (interpolatedMeteoForBand[f] < MIN_HUMIDITY)
                        interpolatedMeteoForBand[f] = MIN_HUMIDITY;
                }
                if (interpolatedMeteoForBasin > MAX_HUMIDITY)
                    interpolatedMeteoForBasin = MAX_HUMIDITY;
                if (interpolatedMeteoForBasin < MIN_HUMIDITY)
                    interpolatedMeteoForBasin = MIN_HUMIDITY;
            } else if (pType == WIND) {
                // V
                double MIN_WIND = 0.01;
                for (int f = 0; f < bandsNum; f++) {
                    if (interpolatedMeteoForBand[f] < MIN_WIND)
                        interpolatedMeteoForBand[f] = MIN_WIND;
                }
                if (interpolatedMeteoForBasin < MIN_WIND)
                    interpolatedMeteoForBasin = MIN_WIND;
            }
        }
        int basinid = ((Number) basinFeatures.get(i).getAttribute(basinIdFieldIndex)).intValue();
        outInterpolatedBand.put(basinid, interpolatedMeteoForBand);
        outInterpolated.put(basinid, new double[] { interpolatedMeteoForBasin });
    }
    pm.done();
}

Example 4 with Range

use of oms3.annotations.Range in project hortonmachine by TheHortonMachine.

the class ComponentAccess method rangeCheck.

public static void rangeCheck(Object comp, boolean in, boolean out) throws Exception {
    ComponentAccess cp = new ComponentAccess(comp);
    Collection<Access> acc = new ArrayList<Access>();
    if (in) {
        acc.addAll(cp.inputs());
    }
    if (out) {
        acc.addAll(cp.outputs());
    }
    for (Access a : acc) {
        String name = a.getField().getName();
        Object val = a.getFieldValue();
        Range range = a.getField().getAnnotation(Range.class);
        if (range != null) {
            if (val instanceof Number) {
                double v = ((Number) val).doubleValue();
                if (!Annotations.inRange(range, v)) {
                    throw new ComponentException(name + " not in range " + v);
                }
            } else if (val instanceof double[]) {
                double[] v = (double[]) val;
                for (int i = 0; i < v.length; i++) {
                    if (!Annotations.inRange(range, v[i])) {
                        throw new ComponentException(name + " not in range " + v[i]);
                    }
                }
            }
        }
    }
}

Example 5 with Range

use of oms3.annotations.Range in project hortonmachine by TheHortonMachine.

the class OmsEnergyBalance method process.

@Execute
public void process() throws Exception {
    outPnet = new HashMap<Integer, double[]>();
    outPrain = new HashMap<Integer, double[]>();
    outPsnow = new HashMap<Integer, double[]>();
    outSwe = new HashMap<Integer, double[]>();
    outNetradiation = new HashMap<Integer, double[]>();
    outNetshortradiation = new HashMap<Integer, double[]>();
    if (safePoint == null)
        safePoint = new SafePoint();
    // retrieve number of bands
    num_EI = 0;
    for (EIEnergy energy : inEnergy) {
        int j = energy.energeticBandId;
        if (j > num_EI) {
            num_EI = j;
        }
    }
    num_EI++;
    num_ES = 0;
    for (EIAreas area : inAreas) {
        int j = area.altimetricBandId;
        if (j > num_ES) {
            num_ES = j;
        }
    }
    num_ES++;
    if (basinid2BasinindexMap == null) {
        // get basin features from feature link
        basinsFeatures = new ArrayList<SimpleFeature>();
        FeatureIterator<SimpleFeature> featureIterator = inBasins.features();
        basinNum = inBasins.size();
        SimpleFeatureType featureType = inBasins.getSchema();
        int basinIdFieldIndex = featureType.indexOf(fBasinid);
        if (basinIdFieldIndex == -1) {
            throw new IllegalArgumentException("The field of the basin id couldn't be found in the supplied basin data.");
        }
        if (fBasinlandcover != null) {
            usoFieldIndex = featureType.indexOf(fBasinlandcover);
            if (usoFieldIndex == -1) {
                throw new IllegalArgumentException("The field of the soil type (usofield) couldn't be found in the supplied basin data.");
            }
        }
        basinid2BasinindexMap = new HashMap<Integer, Integer>();
        basinindex2BasinidMap = new HashMap<Integer, Integer>();
        pm.beginTask("Read basins data.", inBasins.size());
        int index = 0;
        Abasin = new double[basinNum];
        while (featureIterator.hasNext()) {
            pm.worked(1);
            SimpleFeature feature = featureIterator.next();
            basinsFeatures.add(feature);
            basinid2BasinindexMap.put(((Number) feature.getAttribute(basinIdFieldIndex)).intValue(), index);
            basinindex2BasinidMap.put(index, ((Number) feature.getAttribute(basinIdFieldIndex)).intValue());
            Geometry basinGeometry = (Geometry) feature.getDefaultGeometry();
            // area in km2 as the input
            Abasin[index] = basinGeometry.getArea() / 1000000.0;
            // area for energetic and
            // altimetric bands
            index++;
            // read land cover if requested
            if (usoFieldIndex != -1) {
                if (usoList == null) {
                    usoList = new ArrayList<Integer>();
                }
                int uso = ((Number) feature.getAttribute(usoFieldIndex)).intValue();
                usoList.add(uso);
            }
        }
        featureIterator.close();
        pm.done();
    }
    // get rain from scalar link
    double[] rain = new double[basinNum];
    Set<Integer> basinIdSet = inRain.keySet();
    pm.beginTask("Read rain data.", basinIdSet.size());
    for (Integer basinId : basinIdSet) {
        pm.worked(1);
        Integer index = basinid2BasinindexMap.get(basinId);
        if (index == null) {
            continue;
        }
        double[] value = inRain.get(basinId);
        if (!isNovalue(value[0])) {
            rain[index] = value[0];
        } else {
            rain[index] = 0.0;
        }
    }
    pm.done();
    // 0,1,2,3,4,5,5,4,3,2,1,0 ones at the beginning of the simulation
    if (EI == null) {
        EI = new double[12][num_EI][basinNum];
        pm.beginTask("Read energy index data.", inEnergy.size());
        for (EIEnergy energy : inEnergy) {
            pm.worked(1);
            int tempId = energy.basinId;
            Integer index = basinid2BasinindexMap.get(tempId);
            if (index == null) {
                // basinid2BasinindexMap.remove(tempId);
                continue;
            }
            int j = energy.energeticBandId;
            int k = energy.virtualMonth;
            int kInverse = 11 - k;
            EI[k][j][index] = energy.energyValue;
            EI[kInverse][j][index] = energy.energyValue;
        }
        pm.done();
    }
    // beginning of the simulation
    if (A == null) {
        A = new double[num_ES][num_EI][basinNum];
        pm.beginTask("Read area per heigth and band data.", inAreas.size());
        HashMap<Integer, HashMap<Integer, HashMap<Integer, Double>>> idbasinMap = new HashMap<Integer, HashMap<Integer, HashMap<Integer, Double>>>();
        for (EIAreas area : inAreas) {
            int idBasin = area.basinId;
            HashMap<Integer, HashMap<Integer, Double>> idfasceMap = idbasinMap.get(idBasin);
            if (idfasceMap == null) {
                idfasceMap = new HashMap<Integer, HashMap<Integer, Double>>();
                idbasinMap.put(idBasin, idfasceMap);
            }
            int idAltimetricBand = area.altimetricBandId;
            HashMap<Integer, Double> idbandeMap = idfasceMap.get(idAltimetricBand);
            if (idbandeMap == null) {
                idbandeMap = new HashMap<Integer, Double>();
                idfasceMap.put(idAltimetricBand, idbandeMap);
            }
            int idEnergeticBand = area.energyBandId;
            double areaValue = area.areaValue;
            idbandeMap.put(idEnergeticBand, areaValue);
            pm.worked(1);
        }
        pm.done();
        for (int i = 0; i < basinNum; i = i + 1) {
            Integer index = basinindex2BasinidMap.get(i);
            if (index == null) {
                basinid2BasinindexMap.remove(i);
                continue;
            }
            HashMap<Integer, HashMap<Integer, Double>> fasceMap = idbasinMap.get(index);
            for (int j = 0; j < num_ES; j++) {
                HashMap<Integer, Double> bandeMap = fasceMap.get(j);
                for (int k = 0; k < num_EI; k++) {
                    A[j][k][i] = bandeMap.get(k);
                }
            }
        }
    }
    // get T (temperatures per basin per band) from scalar input link at each time step
    double[][] T = null;
    if (inTemp != null) {
        T = new double[basinNum][num_ES];
        pm.beginTask("Read temperature data.", inTemp.size());
        Set<Integer> basinIdsSet = inTemp.keySet();
        for (Integer basinId : basinIdsSet) {
            pm.worked(1);
            Integer index = basinid2BasinindexMap.get(basinId);
            if (index == null) {
                // data for a basin that is not considered, ignore it
                continue;
            }
            double[] values = inTemp.get(basinId);
            T[index] = values;
        }
        pm.done();
    }
    // get V (wind speed per basin per band) from scalar link at each time step
    double[][] V = null;
    if (inWind != null) {
        V = new double[basinNum][num_ES];
        pm.beginTask("Read wind speed data.", inWind.size());
        Set<Integer> basinIdsSet = inWind.keySet();
        for (Integer basinId : basinIdsSet) {
            pm.worked(1);
            Integer index = basinid2BasinindexMap.get(basinId);
            if (index == null) {
                // data for a basin that is not considered, ignore it
                continue;
            }
            double[] values = inWind.get(basinId);
            V[index] = values;
        }
    }
    // get P (pressure per basin per band) from scalar link at each time step
    double[][] P = null;
    if (inPressure != null) {
        P = new double[basinNum][num_ES];
        pm.beginTask("Read pressure data.", inPressure.size());
        Set<Integer> basinIdsSet = inPressure.keySet();
        for (Integer basinId : basinIdsSet) {
            pm.worked(1);
            Integer index = basinid2BasinindexMap.get(basinId);
            if (index == null) {
                // data for a basin that is not considered, ignore it
                continue;
            }
            double[] values = inPressure.get(basinId);
            P[index] = values;
        }
        pm.done();
    }
    // get RH (relative humidity per basin per band) from scalar link at each time step
    double[][] RH = null;
    if (inRh != null) {
        RH = new double[basinNum][num_ES];
        pm.beginTask("Read humidity data.", inRh.size());
        Set<Integer> basinIdsSet = inRh.keySet();
        for (Integer basinId : basinIdsSet) {
            pm.worked(1);
            Integer index = basinid2BasinindexMap.get(basinId);
            if (index == null) {
                // data for a basin that is not considered, ignore it
                continue;
            }
            double[] values = inRh.get(basinId);
            RH[index] = values;
        }
        pm.done();
    }
    // get dtday (daily temperature range per basin per band) from scalar link at each time
    // step
    double[][] DTd = null;
    if (inDtday != null) {
        DTd = new double[basinNum][num_ES];
        pm.beginTask("Read dtday data.", inDtday.size());
        Set<Integer> basinIdsSet = inDtday.keySet();
        for (Integer basinId : basinIdsSet) {
            pm.worked(1);
            Integer index = basinid2BasinindexMap.get(basinId);
            if (index == null) {
                // data for a basin that is not considered, ignore it
                continue;
            }
            double[] values = inDtday.get(basinId);
            DTd[index] = values;
        }
        pm.done();
    }
    // get dtmonth (monthly temperature range per basin per band) from scalar link at each
    // time step
    double[][] DTm = null;
    if (inDtmonth != null) {
        DTm = new double[basinNum][num_ES];
        pm.beginTask("Read dtday data.", inDtmonth.size());
        Set<Integer> basinIdsSet = inDtmonth.keySet();
        for (Integer basinId : basinIdsSet) {
            pm.worked(1);
            Integer index = basinid2BasinindexMap.get(basinId);
            if (index == null) {
                // data for a basin that is not considered, ignore it
                continue;
            }
            double[] values = inDtmonth.get(basinId);
            DTm[index] = values;
        }
        pm.done();
    }
    /*
         * set the current time: day, month and hour
         */
    DateTime currentDatetime = formatter.parseDateTime(tCurrent);
    int currentMonth = currentDatetime.getMonthOfYear();
    int currentDay = currentDatetime.getDayOfMonth();
    int currentMinute = currentDatetime.getMinuteOfDay();
    double hour = currentMinute / 60.0;
    System.out.println("ora: " + hour);
    if (averageTemperature == null) {
        averageTemperature = new double[2 * basinNum];
    } else {
        Arrays.fill(averageTemperature, 0.0);
    }
    /*
         * these have to be taken from initial values 
         */
    if (safePoint.SWE == null) {
        if (pInitsafepoint != null && new File(pInitsafepoint).exists()) {
            safePoint = getSafePointData();
        } else {
            safePoint.SWE = new double[num_ES][num_EI][basinNum];
            if (pInitswe == -9999.0) {
                pInitswe = 0.0;
            }
            for (int i = 0; i < basinNum; i++) {
                double sweTmp = pInitswe;
                if (usoList != null) {
                    int usoTmp = usoList.get(i);
                    if (usoTmp == pGlacierid) {
                        sweTmp = GLACIER_SWE;
                    }
                }
                for (int k = 0; k < num_ES; k++) {
                    for (int j = 0; j < num_EI; j++) {
                        safePoint.SWE[j][k][i] = sweTmp;
                    }
                }
            }
            safePoint.U = new double[num_ES][num_EI][basinNum];
            safePoint.SnAge = new double[num_ES][num_EI][basinNum];
            safePoint.Ts = new double[num_ES][num_EI][basinNum];
        }
    }
    // this has to be taken from a file, scalarreader
    // TODO add the input canopyLink for the canopy height for each altimetric band
    /*
         * if there is no canopy input matrix for the model create an empty canopy matrix for each elevation band and for each basin
         */
    double[][] canopy = new double[num_ES][basinNum];
    for (int i = 0; i < canopy.length; i++) {
        for (int j = 0; j < canopy[0].length; j++) {
            canopy[i][j] = pCanopyh;
        }
    }
    checkParametersAndRunEnergyBalance(rain, T, V, P, RH, currentMonth, currentDay, hour, Abasin, A, EI, DTd, DTm, canopy);
}