diff artifacts/src/main/java/org/dive4elements/river/artifacts/sinfo/flood_duration/FloodDurationCalculator.java @ 9286:bcbae86ce7b3

Improved flood duration curve calculation as specified by BfG
author mschaefer
date Mon, 23 Jul 2018 19:20:06 +0200
parents 83ebeb620b5a
children 7100a555607c
line wrap: on
line diff
--- a/artifacts/src/main/java/org/dive4elements/river/artifacts/sinfo/flood_duration/FloodDurationCalculator.java	Mon Jul 23 19:18:11 2018 +0200
+++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/sinfo/flood_duration/FloodDurationCalculator.java	Mon Jul 23 19:20:06 2018 +0200
@@ -66,8 +66,8 @@
     /**
      * Calculate the infrastructures flood duration result rows
      */
-    public void execute(final Calculation problems, final String label, final DoubleRange calcRange,
-            final RiversideChoiceKey riverside, final WINFOArtifact winfo, final FloodDurationCalculationResults results) {
+    public void execute(final Calculation problems, final String label, final DoubleRange calcRange, final RiversideChoiceKey riverside,
+            final WINFOArtifact winfo, final FloodDurationCalculationResults results) {
 
         // Find all gauges of the calc range, and create the duration finders
         final Map<Gauge, GaugeDurationValuesFinder> durFinders = new HashMap<>();
@@ -91,20 +91,23 @@
         addInfrastructures(allStations, secondBank, infras);
         final double[] stationsSorted = sortStations(allStations.keySet());
 
-        // Calculate W and Q for all stations and the selected discharge states
-        // TODO Geht das schneller, wenn man WstValueTable statt WINFOArtifact.computeWaterlevelData nutzt?
-        final WQKms[] wqkmsArray = calculateWaterlevels(winfo, stationsSorted, problems);
+        // Calculate W and Q for all stations and the selected discharge states/waterlevels
+        final WQKms[] wqkmsArray = calculateWsts(winfo, stationsSorted, problems);
 
-        // Determine discharge state labels of the main values
-        updateMainValueLabels(wqkmsArray, winfo, problems);
+        // Determine discharge state labels of the waterlevels
+        updateWstLabels(wqkmsArray, winfo, problems);
 
-        // Load base wst table (river).wst - first run takes long time, then it's cached
+        final Map<Gauge, List<Double>> gaugeWstDurations = new HashMap<>();
+        calcGaugeWstDurations(winfo, new ArrayList<>(durFinders.keySet()), gaugeWstDurations, durFinders);
+
+        // Load base wst table (river).wst
+        // (should be in cache since already used in calculateWaterlevels (winfo.computeWaterlevelData)
         final WstValueTable wst = WstValueTableFactory.getTable(this.riverInfoProvider.getRiver());
 
-        // Calculate the durations and create the result rows
+        // Create the result rows, and calculate and add the flood durations etc.
         for (int i = 0; i <= stationsSorted.length - 1; i++) {
             final Gauge gauge = this.riverInfoProvider.getGauge(stationsSorted[i], true);
-            final ResultRow row = createRow(stationsSorted[i], gauge, firstGauge, wqkmsArray, durFinders.get(gauge), i);
+            final ResultRow row = createRow(stationsSorted[i], gauge, firstGauge, wqkmsArray, gaugeWstDurations.get(gauge), i);
             if (allStations.containsKey(stationsSorted[i]) && (allStations.get(stationsSorted[i]) != null))
                 calculateInfrastructure(row, gauge, allStations.get(stationsSorted[i]), wst, durFinders);
             this.rows.add(row);
@@ -115,50 +118,39 @@
             }
         }
 
-        //// Create a finder for Q in the {river}.wst km-w-q table
-        // final WQBaseTableFinder wqFinder = WQBaseTableFinder.loadValues(this.riverInfoProvider.getRiver(),
-        //// calcRange.getMinimumDouble(),
-        // calcRange.getMaximumDouble(), problems);
-        //
-        //// Calculate the durations and create the result rows
-        // for (int i = 0; i <= stationsSorted.length - 1; i++) {
-        // final Gauge gauge = this.riverInfoProvider.getGauge(stationsSorted[i], true);
-        // final ResultRow row = createRow(stationsSorted[i], gauge, firstGauge, wqkmsArray, durFinders.get(gauge), i);
-        // if (allStations.containsKey(stationsSorted[i]) && (allStations.get(stationsSorted[i]) != null))
-        // calculateInfrastructure(row, gauge, allStations.get(stationsSorted[i]), wqFinder, durFinders);
-        // this.rows.add(row);
-        // if (secondBank.containsKey(stationsSorted[i])) {
-        // final ResultRow row2 = ResultRow.create(row);
-        // calculateInfrastructure(row2, gauge, secondBank.get(stationsSorted[i]), wqFinder, durFinders);
-        // this.rows.add(row2);
-        // }
-        // }
+        // Get the labels of the selected waterlevels
+        final String[] wstLabels = new String[wqkmsArray.length];
+        for (int i = 0; i <= wqkmsArray.length - 1; i++)
+            wstLabels[i] = wqkmsArray[i].getName();
 
-        final String[] mainValueLabels = new String[wqkmsArray.length];
-        for (int i = 0; i <= wqkmsArray.length - 1; i++)
-            mainValueLabels[i] = wqkmsArray[i].getName();
-        results.addResult(new FloodDurationCalculationResult(label, mainValueLabels, this.rows), problems);
+        results.addResult(new FloodDurationCalculationResult(label, wstLabels, this.rows), problems);
     }
 
     /**
      * Calculates the duration curve for a station
+     * (other than the version 3.2.1 W-Info Dauerlinie the wst column positions
+     * are taken from the Q values of the gauge's Q-D-table)
      */
     public WQDay calcWQDays(final Calculation problems, final double station, final WINFOArtifact winfo) {
 
-        // Same processing as in W-Info DurationCurveState
-        winfo.addStringData("ld_locations", Double.toString(station));
-        final CalculationResult res = winfo.getDurationCurveData();
-        final WQDay underflow = (WQDay) res.getData();
-        // Transform underflow days into overflow days and re-sort
-        final int[] days = new int[underflow.getWs().length];
-        final double[] ws = new double[days.length];
-        final double[] qs = new double[days.length];
-        for (int i = 0, j = days.length - 1; i <= days.length - 1; i++, j--) {
-            days[j] = 365 - underflow.getDay(i);
-            ws[j] = underflow.getW(i);
-            qs[j] = underflow.getQ(i);
+        final WstValueTable wst = WstValueTableFactory.getTable(this.riverInfoProvider.getRiver());
+        final Gauge gauge = this.riverInfoProvider.getRiver().determineGaugeByPosition(station);
+        final Object[] obj = gauge.fetchDurationCurveData();
+        final int[] udays = (int[]) obj[0];
+        final double[] qs = (double[]) obj[1];
+        final int[] odays = new int[udays.length];
+        final double[] oqs = new double[udays.length];
+        final double[] ows = new double[udays.length];
+        for (int i = 0, j = udays.length - 1; i <= udays.length - 1; i++, j--) {
+            odays[j] = 365 - udays[i];
+            oqs[j] = qs[i];
+            final QPosition qpos = wst.getQPosition(gauge.getStation().doubleValue(), qs[i]);
+            if (qpos != null)
+                ows[j] = wst.interpolateW(station, qpos, problems);
+            else
+                ows[j] = Double.NaN;
         }
-        return new WQDay(days, ws, qs);
+        return new WQDay(odays, ows, oqs);
     }
 
     /**
@@ -341,9 +333,8 @@
     /**
      * Calculates an array of w-q-longitudinal sections for all artifact W/Q options
      */
-    private WQKms[] calculateWaterlevels(final WINFOArtifact winfo, final double[] stations, final Calculation problems) {
-        // REMARK aus TkhCalculation - move to WinfoArtifactWrapper?
-        // TODO das ist ziemlich langsam - durch den WQBaseTableFinder ersetzen? (vorher W-Optionen in Q umrechnen)
+    private WQKms[] calculateWsts(final WINFOArtifact winfo, final double[] stations, final Calculation problems) {
+        // First run may take long, further runs are faster since WstValueTable is in cache then
         // (So funktioniert computeWaterlevelData wohl:
         // Es sucht die Spalte(n) zum Bezugspegel-Q in der W-Q-Tabelle ({river}.wst in Wst etc.),
         // interpoliert die horizontale Tabellenposition (Q) und dann die vertikale Tabellenposition der station;
@@ -352,36 +343,35 @@
         // interpoliert;
         // bei Vorgabe eines W auf freier Strecke wird wohl vorher noch die .wst-Interpolation eingesetzt, um das Q zu bekommen.
 
-        final CalculationResult waterlevelData = winfo.computeWaterlevelData(stations);
+        final CalculationResult wstsData = winfo.computeWaterlevelData(stations);
 
         /* copy all problems */
-        final Calculation winfoProblems = waterlevelData.getReport();
+        final Calculation winfoProblems = wstsData.getReport();
         final List<Problem> problems2 = winfoProblems.getProblems();
         if (problems2 != null) {
             for (final Problem problem : problems2) {
                 problems.addProblem(problem);
             }
         }
-        return (WQKms[]) waterlevelData.getData();
+        return (WQKms[]) wstsData.getData();
     }
 
     /**
      * Determines the waterlevel/discharge state labels for the selected Q or W values and sets them in the WQKms array
      */
-    private void updateMainValueLabels(final WQKms[] wqkmsArray, final WINFOArtifact winfo, final Calculation problems) {
+    private void updateWstLabels(final WQKms[] wqkmsArray, final WINFOArtifact winfo, final Calculation problems) {
 
         for (int i = 0; i <= wqkmsArray.length - 1; i++)
             wqkmsArray[i].setName(buildWQDescription(wqkmsArray[i], winfo));
     }
 
     /**
-     *
-     * @param wqkms
-     * @param descBuilder
-     * @return
+     * Builds the description label of a waterlevel
      */
     private String buildWQDescription(final WQKms wqkms, final WINFOArtifact winfo) {
+
         final WaterlevelDescriptionBuilder descBuilder = new WaterlevelDescriptionBuilder(winfo, this.context);
+        // TODO Zwischen numerischem Q-Wert und Dauerzahl-Hauptwert (0..364) unterscheiden
         final String description = descBuilder.getDesc(wqkms);
         if (!description.isEmpty() && Character.isDigit(description.charAt(0))) {
             if (winfo.isQ())
@@ -394,10 +384,32 @@
     }
 
     /**
+     * Calculates the flood durations of the Qs of the waterlevels/discharge states for a map of gauges
+     */
+    private void calcGaugeWstDurations(final WINFOArtifact winfo, final List<Gauge> gauges, final Map<Gauge, List<Double>> gaugeWstDurations,
+            final Map<Gauge, GaugeDurationValuesFinder> durFinders) {
+
+        final double[] gaugeKms = new double[gauges.size()];
+        for (int i = 0; i <= gauges.size() - 1; i++) {
+            gaugeKms[i] = gauges.get(i).getStation().doubleValue();
+            gaugeWstDurations.put(gauges.get(i), new ArrayList<Double>());
+        }
+        final CalculationResult wstsData = winfo.computeWaterlevelData(gaugeKms);
+        final WQKms[] wsts = (WQKms[]) wstsData.getData();
+        for (int i = 0; i <= gauges.size() - 1; i++) {
+            final GaugeDurationValuesFinder durFinder = durFinders.get(gauges.get(i));
+            for (int j = 0; j <= wsts.length - 1; j++) {
+                final double d = durFinder.getDuration(wsts[j].getQ(i));
+                gaugeWstDurations.get(gauges.get(i)).add(Double.valueOf(underflowDaysToOverflowDays(d)));
+            }
+        }
+    }
+
+    /**
      * Create a result row for a station and its gauge, and add w-q-values as selected
      */
     private ResultRow createRow(final Double station, final Gauge gauge, final Gauge firstGauge, final WQKms[] wqkmsArray,
-            final GaugeDurationValuesFinder durationFinder, final int kmIndex) {
+            final List<Double> gaugeDurations, final int kmIndex) {
 
         final ResultRow row = ResultRow.create();
         row.putValue(GeneralResultType.station, station);
@@ -410,39 +422,22 @@
         final String location = this.riverInfoProvider.getLocation(station);
         row.putValue(SInfoResultType.location, location);
 
-        final List<DurationWaterlevel> waterlevels = new ArrayList<>(wqkmsArray.length);
-
-        for (final WQKms wqKms : wqkmsArray) {
-            assert (wqKms.getKm(kmIndex) == station.doubleValue());
+        final List<DurationWaterlevel> wsts = new ArrayList<>(wqkmsArray.length);
 
-            final int overflowDays = (int) Math.round(underflowDaysToOverflowDays(durationFinder.getDuration(wqKms.getQ(kmIndex))));
+        for (int i = 0; i <= wqkmsArray.length - 1; i++) {
+            assert (wqkmsArray[i].getKm(kmIndex) == station.doubleValue());
 
-            final DurationWaterlevel dw = new DurationWaterlevel(wqKms.getW(kmIndex), overflowDays, wqKms.getQ(kmIndex), wqKms.getName());
-            waterlevels.add(dw);
+            final int overflowDays = (int) Math.round(gaugeDurations.get(i));
+
+            final DurationWaterlevel dw = new DurationWaterlevel(wqkmsArray[i].getW(kmIndex), overflowDays, wqkmsArray[i].getQ(kmIndex),
+                    wqkmsArray[i].getName());
+            wsts.add(dw);
         }
-        row.putValue(SInfoResultType.customMultiRowColWaterlevel, waterlevels);
+        row.putValue(SInfoResultType.customMultiRowColWaterlevel, wsts);
 
         return row;
     }
 
-    /// **
-    // * Calculate the result row fields for one infrastructure - gives wrong duration numbers where Q changes within the
-    /// gauge range
-    // */
-    // private void calculateInfrastructure(final ResultRow row, final Gauge gauge, final InfrastructureValue
-    /// infrastructure, final WQBaseTableFinder wqFinder,
-    // final Map<Gauge, GaugeDurationValuesFinder> durFinders) {
-    //
-    // final double q = wqFinder.getDischarge(infrastructure.getStation(), infrastructure.getHeight());
-    // final double qOut = Double.isInfinite(q) ? Double.NaN : q;
-    // final double dur = underflowDaysToOverflowDays(durFinders.get(gauge).getDuration(q));
-    // row.putValue(SInfoResultType.riverside, infrastructure.getAttributeKey());
-    // row.putValue(SInfoResultType.floodDuration, dur);
-    // row.putValue(SInfoResultType.floodDischarge, qOut);
-    // row.putValue(SInfoResultType.infrastructureHeight, infrastructure.getHeight());
-    // row.putValue(SInfoResultType.infrastructuretype, infrastructure.getInfrastructure().getType().getName());
-    // }
-
     /**
      * Calculate the result row fields for one infrastructure
      */
@@ -454,17 +449,16 @@
         final double[] qs = wst.findQsForW(infrastructure.getStation().doubleValue(), infrastructure.getHeight().doubleValue(), problems);
         // TODO Fehlerbehandlung (kein Q gefunden)
         final double q = (qs.length >= 1) ? qs[0] : Double.NaN;
-        final double qOut = Double.isInfinite(q) ? Double.NaN : q;
-        // Determine the relative column position of the Q
+        // Determine the relative column position of the Q of the infrastructure height
         final QPosition qPos = wst.getQPosition(infrastructure.getStation().doubleValue(), q);
         // Get the Q for the found column position for the station of the gauge
         final double qGauge = wst.getQ(qPos, gauge.getStation().doubleValue());
         // Interpolate the Q-D-table of the gauge
         final double dur = underflowDaysToOverflowDays(durFinders.get(gauge).getDuration(qGauge));
-        // Set result row
+        // Set the result row
         row.putValue(SInfoResultType.riverside, infrastructure.getAttributeKey());
         row.putValue(SInfoResultType.floodDuration, dur);
-        row.putValue(SInfoResultType.floodDischarge, qOut);
+        row.putValue(SInfoResultType.floodDischarge, q);
         row.putValue(SInfoResultType.infrastructureHeight, infrastructure.getHeight());
         row.putValue(SInfoResultType.infrastructuretype, infrastructure.getInfrastructure().getType().getName());
     }

http://dive4elements.wald.intevation.org