Examples with EIAreas org.hortonmachine.gears.io.eicalculator.EIAreas used on opensource projects

Example 1 with EIAreas

use of org.hortonmachine.gears.io.eicalculator.EIAreas 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 2 with EIAreas

use of org.hortonmachine.gears.io.eicalculator.EIAreas in project hortonmachine by TheHortonMachine.

the class OmsJami method calculateAreas.

private void calculateAreas(int fascieNum) {
    if (basinAreas == null) {
        pm.beginTask("Calculate areas per basin and altimetric band.", inAreas.size());
        HashMap<Integer, HashMap<Integer, Double>> idbasin2InfoMap = new HashMap<Integer, HashMap<Integer, Double>>();
        HashMap<Integer, Double> idbasin2AreaMap = new HashMap<Integer, Double>();
        for (EIAreas area : inAreas) {
            int idBasin = area.basinId;
            HashMap<Integer, Double> idfasceMap = idbasin2InfoMap.get(idBasin);
            if (idfasceMap == null) {
                idfasceMap = new HashMap<Integer, Double>();
                idbasin2InfoMap.put(idBasin, idfasceMap);
            }
            int idAltimetricBand = area.altimetricBandId;
            double areaValue = area.areaValue;
            // area fascia
            Double areaFascia = idfasceMap.get(idAltimetricBand);
            if (areaFascia == null) {
                idfasceMap.put(idAltimetricBand, areaValue);
            } else {
                Double sum = areaValue + areaFascia;
                idfasceMap.put(idAltimetricBand, sum);
            }
            // total basin area
            Double basinArea = idbasin2AreaMap.get(idBasin);
            if (basinArea == null) {
                idbasin2AreaMap.put(idBasin, areaValue);
            } else {
                Double sum = areaValue + basinArea;
                idbasin2AreaMap.put(idBasin, sum);
            }
            pm.worked(1);
        }
        pm.done();
        basinAreas = new double[idbasin2AreaMap.size()];
        basinAreasPerFascias = new double[idbasin2AreaMap.size()][fascieNum];
        Set<Entry<Integer, Integer>> entrySet = basinid2BasinindexMap.entrySet();
        for (Entry<Integer, Integer> entry : entrySet) {
            Integer basinId = entry.getKey();
            Integer basinIndex = entry.getValue();
            Double area = idbasin2AreaMap.get(basinId);
            basinAreas[basinIndex] = area;
            HashMap<Integer, Double> fascia2AreaMap = idbasin2InfoMap.get(basinId);
            for (int fasciaIndex = 0; fasciaIndex < fascieNum; fasciaIndex++) {
                basinAreasPerFascias[basinIndex][fasciaIndex] = fascia2AreaMap.get(fasciaIndex);
            }
        }
    }
}

Example 3 with EIAreas

use of org.hortonmachine.gears.io.eicalculator.EIAreas 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);
}

Example 4 with EIAreas

use of org.hortonmachine.gears.io.eicalculator.EIAreas in project hortonmachine by TheHortonMachine.

the class OmsEnergyIndexCalculator method process.

@Execute
public void process() throws Exception {
    HashMap<String, Double> regionMap = CoverageUtilities.getRegionParamsFromGridCoverage(inBasins);
    cols = regionMap.get(CoverageUtilities.COLS).intValue();
    rows = regionMap.get(CoverageUtilities.ROWS).intValue();
    dx = regionMap.get(CoverageUtilities.XRES);
    dy = regionMap.get(CoverageUtilities.YRES);
    double n = regionMap.get(CoverageUtilities.NORTH);
    double s = regionMap.get(CoverageUtilities.SOUTH);
    double w = regionMap.get(CoverageUtilities.WEST);
    double e = regionMap.get(CoverageUtilities.EAST);
    double meanX = w + (e - w) / 2.0;
    double meanY = s + (n - s) / 2.0;
    Coordinate tmp = new Coordinate(meanX, meanY);
    MathTransform mathTransform = CRS.findMathTransform(inAspect.getCoordinateReferenceSystem(), DefaultGeographicCRS.WGS84);
    Coordinate newC = JTS.transform(tmp, null, mathTransform);
    avgLatitude = newC.y;
    RenderedImage idbasinImage = inBasins.getRenderedImage();
    idbasinImageIterator = RandomIterFactory.create(idbasinImage, null);
    RenderedImage elevImage = inElev.getRenderedImage();
    elevImageIterator = RandomIterFactory.create(elevImage, null);
    RenderedImage tmpImage = inCurvatures.getRenderedImage();
    curvatureImage = CoverageUtilities.createWritableRaster(tmpImage.getWidth(), tmpImage.getHeight(), null, null, null);
    RandomIter tmpIterator = RandomIterFactory.create(tmpImage, null);
    // TODO check what this is for?!?!?
    for (int i = 0; i < tmpImage.getHeight(); i++) {
        for (int j = 0; j < tmpImage.getWidth(); j++) {
            double value = tmpIterator.getSampleDouble(j, i, 0);
            curvatureImage.setSample(j, i, 0, value);
        }
    }
    RenderedImage aspectImage = inAspect.getRenderedImage();
    aspectImageIterator = RandomIterFactory.create(aspectImage, null);
    RenderedImage slopeImage = inSlope.getRenderedImage();
    slopeImageIterator = RandomIterFactory.create(slopeImage, null);
    avgLatitude *= (PI / 180.0);
    // $NON-NLS-1$
    pm.message(msg.message("eicalculator.preparing_inputs"));
    eibasinNum = prepareInputsOutputs();
    // $NON-NLS-1$
    pm.beginTask(msg.message("eicalculator.computing"), 6);
    for (int m = 0; m < 6; m++) {
        pm.worked(1);
        compute_EI(m + 1);
    }
    pm.done();
    average_EI(10, 6);
    // $NON-NLS-1$
    pm.beginTask(msg.message("eicalculator.calc_areas"), eibasinNum);
    for (int i = 0; i < eibasinNum; i++) {
        pm.worked(1);
        area(i);
    }
    pm.done();
    /*
         * putting the results together
         */
    outAltimetry = new ArrayList<EIAltimetry>();
    outEnergy = new ArrayList<EIEnergy>();
    outArea = new ArrayList<EIAreas>();
    for (int i = 0; i < eibasinNum; i++) {
        int realBasinId = index2idMap.get(i + 1);
        /*
             * ENERGY BANDS
             * 
             * Cycle over the virtual months:
             * 0: 22 DICEMBRE - 20 GENNAIO
             * 1: 21 GENNAIO - 20 FEBBRAIO
             * 2: 21 FEBBRAIO - 22 MARZO
             * 3: 23 MARZO - 22 APRILE
             * 4: 23 APRILE - 22 MAGGIO
             * 5: 23 MAGGIO - 22 GIUGNO
             */
        for (int j = 0; j < 6; j++) {
            for (int k = 0; k < pEi; k++) {
                EIEnergy tmpEi = new EIEnergy();
                // the basin id
                tmpEi.basinId = realBasinId;
                tmpEi.energeticBandId = k;
                tmpEi.virtualMonth = j;
                tmpEi.energyValue = eibasinEI[0][k][i];
                outEnergy.add(tmpEi);
            }
        }
        /*
             * ALTIMETRIC BANDS
             */
        for (int k = 0; k < pEs; k++) {
            EIAltimetry tmpAl = new EIAltimetry();
            tmpAl.basinId = realBasinId;
            tmpAl.altimetricBandId = k;
            tmpAl.elevationValue = eibasinES[k][i];
            tmpAl.bandRange = eibasinESrange[k][i];
            outAltimetry.add(tmpAl);
        }
        /*
             * AREAS
             */
        for (int j = 0; j < pEs; j++) {
            for (int k = 0; k < pEi; k++) {
                EIAreas tmpAr = new EIAreas();
                tmpAr.basinId = realBasinId;
                tmpAr.altimetricBandId = j;
                tmpAr.energyBandId = k;
                tmpAr.areaValue = eibasinA[j][k][i];
                outArea.add(tmpAr);
            }
        }
    }
}

Example 5 with EIAreas

use of org.hortonmachine.gears.io.eicalculator.EIAreas in project hortonmachine by TheHortonMachine.

the class TestEnergyIndexCalculator method testEnergyIndexCalculator.

/**
 * TODO make this test a bit more serious.
 *
 * @throws Exception
 */
public void testEnergyIndexCalculator() throws Exception {
    HashMap<String, Double> envelopeParams = HMTestMaps.getEnvelopeparams();
    CoordinateReferenceSystem crs = HMTestMaps.getCrs();
    // PrintStreamProgressMonitor pm = new PrintStreamProgressMonitor(System.out, System.out);
    double[][] aspectData = HMTestMaps.aspectDataRadiants;
    GridCoverage2D aspectCoverage = CoverageUtilities.buildCoverage("aspect", aspectData, envelopeParams, crs, true);
    double[][] nablaData = HMTestMaps.nablaData0;
    GridCoverage2D nablaCoverage = CoverageUtilities.buildCoverage("nabla", nablaData, envelopeParams, crs, true);
    double[][] pitData = HMTestMaps.pitData;
    GridCoverage2D pitCoverage = CoverageUtilities.buildCoverage("pit", pitData, envelopeParams, crs, true);
    double[][] slopeData = HMTestMaps.slopeData;
    GridCoverage2D slopeCoverage = CoverageUtilities.buildCoverage("slope", slopeData, envelopeParams, crs, true);
    double[][] subbasinsData = HMTestMaps.basinDataNN0;
    GridCoverage2D subbasinsCoverage = CoverageUtilities.buildCoverage("subbasins", subbasinsData, envelopeParams, crs, true);
    OmsEnergyIndexCalculator eiCalculator = new OmsEnergyIndexCalculator();
    eiCalculator.inAspect = aspectCoverage;
    eiCalculator.inCurvatures = nablaCoverage;
    eiCalculator.inElev = pitCoverage;
    eiCalculator.inSlope = slopeCoverage;
    eiCalculator.inBasins = subbasinsCoverage;
    eiCalculator.pDt = 1;
    eiCalculator.pEi = 2;
    eiCalculator.pEs = 2;
    eiCalculator.pm = new DummyProgressMonitor();
    eiCalculator.process();
    List<EIAltimetry> altimetricValues = eiCalculator.outAltimetry;
    List<EIEnergy> energeticValues = eiCalculator.outEnergy;
    List<EIAreas> areaValues = eiCalculator.outArea;
    EIAltimetry eiAltimetry = altimetricValues.get(0);
    assertEquals(1, eiAltimetry.basinId);
    assertEquals(0, eiAltimetry.altimetricBandId);
    assertEquals(737.5, eiAltimetry.elevationValue);
    assertEquals(75.0, eiAltimetry.bandRange);
    EIEnergy eiEnergy = energeticValues.get(0);
    assertEquals(1, eiEnergy.basinId);
    assertEquals(0, eiEnergy.energeticBandId);
    assertEquals(0, eiEnergy.virtualMonth);
    assertEquals(0.09808943859674346, eiEnergy.energyValue, 0.0001);
    EIAreas eiAreas = areaValues.get(0);
    assertEquals(1, eiAreas.basinId);
    assertEquals(0, eiAreas.altimetricBandId);
    assertEquals(0, eiAreas.energyBandId);
    assertEquals(0.0, eiAreas.areaValue, 0.0001);
}