Mercurial > dive4elements > river

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

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Improved flood duration curve calculation as specified by BfG
author mschaefer
date Mon, 23 Jul 2018 19:20:06 +0200
parents 83ebeb620b5a
children 7100a555607c
comparison
equal deleted inserted replaced
9285:9b16f58c62a7 9286:bcbae86ce7b3
64 } 64 }
65 65
66 /** 66 /**
67 * Calculate the infrastructures flood duration result rows 67 * Calculate the infrastructures flood duration result rows
68 */ 68 */
69 public void execute(final Calculation problems, final String label, final DoubleRange calcRange, 69 public void execute(final Calculation problems, final String label, final DoubleRange calcRange, final RiversideChoiceKey riverside,
70 final RiversideChoiceKey riverside, final WINFOArtifact winfo, final FloodDurationCalculationResults results) { 70 final WINFOArtifact winfo, final FloodDurationCalculationResults results) {
71 71
72 // Find all gauges of the calc range, and create the duration finders 72 // Find all gauges of the calc range, and create the duration finders
73 final Map<Gauge, GaugeDurationValuesFinder> durFinders = new HashMap<>(); 73 final Map<Gauge, GaugeDurationValuesFinder> durFinders = new HashMap<>();
74 Gauge firstGauge = null; 74 Gauge firstGauge = null;
75 for (final Gauge gauge : this.riverInfoProvider.getRiver().determineGauges(calcRange.getMinimumDouble(), calcRange.getMaximumDouble())) { 75 for (final Gauge gauge : this.riverInfoProvider.getRiver().determineGauges(calcRange.getMinimumDouble(), calcRange.getMaximumDouble())) {
89 addRangeStations(allStations, winfo); 89 addRangeStations(allStations, winfo);
90 addGaugeLimits(allStations, durFinders.keySet(), calcRange.getMinimumDouble(), calcRange.getMaximumDouble()); 90 addGaugeLimits(allStations, durFinders.keySet(), calcRange.getMinimumDouble(), calcRange.getMaximumDouble());
91 addInfrastructures(allStations, secondBank, infras); 91 addInfrastructures(allStations, secondBank, infras);
92 final double[] stationsSorted = sortStations(allStations.keySet()); 92 final double[] stationsSorted = sortStations(allStations.keySet());
93 93
94 // Calculate W and Q for all stations and the selected discharge states 94 // Calculate W and Q for all stations and the selected discharge states/waterlevels
95 // TODO Geht das schneller, wenn man WstValueTable statt WINFOArtifact.computeWaterlevelData nutzt? 95 final WQKms[] wqkmsArray = calculateWsts(winfo, stationsSorted, problems);
96 final WQKms[] wqkmsArray = calculateWaterlevels(winfo, stationsSorted, problems); 96
97 97 // Determine discharge state labels of the waterlevels
98 // Determine discharge state labels of the main values 98 updateWstLabels(wqkmsArray, winfo, problems);
99 updateMainValueLabels(wqkmsArray, winfo, problems); 99
100 100 final Map<Gauge, List<Double>> gaugeWstDurations = new HashMap<>();
101 // Load base wst table (river).wst - first run takes long time, then it's cached 101 calcGaugeWstDurations(winfo, new ArrayList<>(durFinders.keySet()), gaugeWstDurations, durFinders);
102
103 // Load base wst table (river).wst
104 // (should be in cache since already used in calculateWaterlevels (winfo.computeWaterlevelData)
102 final WstValueTable wst = WstValueTableFactory.getTable(this.riverInfoProvider.getRiver()); 105 final WstValueTable wst = WstValueTableFactory.getTable(this.riverInfoProvider.getRiver());
103 106
104 // Calculate the durations and create the result rows 107 // Create the result rows, and calculate and add the flood durations etc.
105 for (int i = 0; i <= stationsSorted.length - 1; i++) { 108 for (int i = 0; i <= stationsSorted.length - 1; i++) {
106 final Gauge gauge = this.riverInfoProvider.getGauge(stationsSorted[i], true); 109 final Gauge gauge = this.riverInfoProvider.getGauge(stationsSorted[i], true);
107 final ResultRow row = createRow(stationsSorted[i], gauge, firstGauge, wqkmsArray, durFinders.get(gauge), i); 110 final ResultRow row = createRow(stationsSorted[i], gauge, firstGauge, wqkmsArray, gaugeWstDurations.get(gauge), i);
108 if (allStations.containsKey(stationsSorted[i]) && (allStations.get(stationsSorted[i]) != null)) 111 if (allStations.containsKey(stationsSorted[i]) && (allStations.get(stationsSorted[i]) != null))
109 calculateInfrastructure(row, gauge, allStations.get(stationsSorted[i]), wst, durFinders); 112 calculateInfrastructure(row, gauge, allStations.get(stationsSorted[i]), wst, durFinders);
110 this.rows.add(row); 113 this.rows.add(row);
111 if (secondBank.containsKey(stationsSorted[i])) { 114 if (secondBank.containsKey(stationsSorted[i])) {
112 final ResultRow row2 = ResultRow.create(row); 115 final ResultRow row2 = ResultRow.create(row);
113 calculateInfrastructure(row2, gauge, secondBank.get(stationsSorted[i]), wst, durFinders); 116 calculateInfrastructure(row2, gauge, secondBank.get(stationsSorted[i]), wst, durFinders);
114 this.rows.add(row2); 117 this.rows.add(row2);
115 } 118 }
116 } 119 }
117 120
118 //// Create a finder for Q in the {river}.wst km-w-q table 121 // Get the labels of the selected waterlevels
119 // final WQBaseTableFinder wqFinder = WQBaseTableFinder.loadValues(this.riverInfoProvider.getRiver(), 122 final String[] wstLabels = new String[wqkmsArray.length];
120 //// calcRange.getMinimumDouble(),
121 // calcRange.getMaximumDouble(), problems);
122 //
123 //// Calculate the durations and create the result rows
124 // for (int i = 0; i <= stationsSorted.length - 1; i++) {
125 // final Gauge gauge = this.riverInfoProvider.getGauge(stationsSorted[i], true);
126 // final ResultRow row = createRow(stationsSorted[i], gauge, firstGauge, wqkmsArray, durFinders.get(gauge), i);
127 // if (allStations.containsKey(stationsSorted[i]) && (allStations.get(stationsSorted[i]) != null))
128 // calculateInfrastructure(row, gauge, allStations.get(stationsSorted[i]), wqFinder, durFinders);
129 // this.rows.add(row);
130 // if (secondBank.containsKey(stationsSorted[i])) {
131 // final ResultRow row2 = ResultRow.create(row);
132 // calculateInfrastructure(row2, gauge, secondBank.get(stationsSorted[i]), wqFinder, durFinders);
133 // this.rows.add(row2);
134 // }
135 // }
136
137 final String[] mainValueLabels = new String[wqkmsArray.length];
138 for (int i = 0; i <= wqkmsArray.length - 1; i++) 123 for (int i = 0; i <= wqkmsArray.length - 1; i++)
139 mainValueLabels[i] = wqkmsArray[i].getName(); 124 wstLabels[i] = wqkmsArray[i].getName();
140 results.addResult(new FloodDurationCalculationResult(label, mainValueLabels, this.rows), problems); 125
126 results.addResult(new FloodDurationCalculationResult(label, wstLabels, this.rows), problems);
141 } 127 }
142 128
143 /** 129 /**
144 * Calculates the duration curve for a station 130 * Calculates the duration curve for a station
131 * (other than the version 3.2.1 W-Info Dauerlinie the wst column positions
132 * are taken from the Q values of the gauge's Q-D-table)
145 */ 133 */
146 public WQDay calcWQDays(final Calculation problems, final double station, final WINFOArtifact winfo) { 134 public WQDay calcWQDays(final Calculation problems, final double station, final WINFOArtifact winfo) {
147 135
148 // Same processing as in W-Info DurationCurveState 136 final WstValueTable wst = WstValueTableFactory.getTable(this.riverInfoProvider.getRiver());
149 winfo.addStringData("ld_locations", Double.toString(station)); 137 final Gauge gauge = this.riverInfoProvider.getRiver().determineGaugeByPosition(station);
150 final CalculationResult res = winfo.getDurationCurveData(); 138 final Object[] obj = gauge.fetchDurationCurveData();
151 final WQDay underflow = (WQDay) res.getData(); 139 final int[] udays = (int[]) obj[0];
152 // Transform underflow days into overflow days and re-sort 140 final double[] qs = (double[]) obj[1];
153 final int[] days = new int[underflow.getWs().length]; 141 final int[] odays = new int[udays.length];
154 final double[] ws = new double[days.length]; 142 final double[] oqs = new double[udays.length];
155 final double[] qs = new double[days.length]; 143 final double[] ows = new double[udays.length];
156 for (int i = 0, j = days.length - 1; i <= days.length - 1; i++, j--) { 144 for (int i = 0, j = udays.length - 1; i <= udays.length - 1; i++, j--) {
157 days[j] = 365 - underflow.getDay(i); 145 odays[j] = 365 - udays[i];
158 ws[j] = underflow.getW(i); 146 oqs[j] = qs[i];
159 qs[j] = underflow.getQ(i); 147 final QPosition qpos = wst.getQPosition(gauge.getStation().doubleValue(), qs[i]);
160 } 148 if (qpos != null)
161 return new WQDay(days, ws, qs); 149 ows[j] = wst.interpolateW(station, qpos, problems);
150 else
151 ows[j] = Double.NaN;
152 }
153 return new WQDay(odays, ows, oqs);
162 } 154 }
163 155
164 /** 156 /**
165 * Calculates the data for the Q main value lines in the duration curve chart 157 * Calculates the data for the Q main value lines in the duration curve chart
166 */ 158 */
339 } 331 }
340 332
341 /** 333 /**
342 * Calculates an array of w-q-longitudinal sections for all artifact W/Q options 334 * Calculates an array of w-q-longitudinal sections for all artifact W/Q options
343 */ 335 */
344 private WQKms[] calculateWaterlevels(final WINFOArtifact winfo, final double[] stations, final Calculation problems) { 336 private WQKms[] calculateWsts(final WINFOArtifact winfo, final double[] stations, final Calculation problems) {
345 // REMARK aus TkhCalculation - move to WinfoArtifactWrapper? 337 // First run may take long, further runs are faster since WstValueTable is in cache then
346 // TODO das ist ziemlich langsam - durch den WQBaseTableFinder ersetzen? (vorher W-Optionen in Q umrechnen)
347 // (So funktioniert computeWaterlevelData wohl: 338 // (So funktioniert computeWaterlevelData wohl:
348 // Es sucht die Spalte(n) zum Bezugspegel-Q in der W-Q-Tabelle ({river}.wst in Wst etc.), 339 // Es sucht die Spalte(n) zum Bezugspegel-Q in der W-Q-Tabelle ({river}.wst in Wst etc.),
349 // interpoliert die horizontale Tabellenposition (Q) und dann die vertikale Tabellenposition der station; 340 // interpoliert die horizontale Tabellenposition (Q) und dann die vertikale Tabellenposition der station;
350 // das ergibt das W einer station für einen Abflusszustand; 341 // das ergibt das W einer station für einen Abflusszustand;
351 // bei Vorgabe eines Pegel-W wird vorher anhand der W-Q-Tabelle des Pegels ({gauge}.at in DischargeTable) das Q 342 // bei Vorgabe eines Pegel-W wird vorher anhand der W-Q-Tabelle des Pegels ({gauge}.at in DischargeTable) das Q
352 // interpoliert; 343 // interpoliert;
353 // bei Vorgabe eines W auf freier Strecke wird wohl vorher noch die .wst-Interpolation eingesetzt, um das Q zu bekommen. 344 // bei Vorgabe eines W auf freier Strecke wird wohl vorher noch die .wst-Interpolation eingesetzt, um das Q zu bekommen.
354 345
355 final CalculationResult waterlevelData = winfo.computeWaterlevelData(stations); 346 final CalculationResult wstsData = winfo.computeWaterlevelData(stations);
356 347
357 /* copy all problems */ 348 /* copy all problems */
358 final Calculation winfoProblems = waterlevelData.getReport(); 349 final Calculation winfoProblems = wstsData.getReport();
359 final List<Problem> problems2 = winfoProblems.getProblems(); 350 final List<Problem> problems2 = winfoProblems.getProblems();
360 if (problems2 != null) { 351 if (problems2 != null) {
361 for (final Problem problem : problems2) { 352 for (final Problem problem : problems2) {
362 problems.addProblem(problem); 353 problems.addProblem(problem);
363 } 354 }
364 } 355 }
365 return (WQKms[]) waterlevelData.getData(); 356 return (WQKms[]) wstsData.getData();
366 } 357 }
367 358
368 /** 359 /**
369 * Determines the waterlevel/discharge state labels for the selected Q or W values and sets them in the WQKms array 360 * Determines the waterlevel/discharge state labels for the selected Q or W values and sets them in the WQKms array
370 */ 361 */
371 private void updateMainValueLabels(final WQKms[] wqkmsArray, final WINFOArtifact winfo, final Calculation problems) { 362 private void updateWstLabels(final WQKms[] wqkmsArray, final WINFOArtifact winfo, final Calculation problems) {
372 363
373 for (int i = 0; i <= wqkmsArray.length - 1; i++) 364 for (int i = 0; i <= wqkmsArray.length - 1; i++)
374 wqkmsArray[i].setName(buildWQDescription(wqkmsArray[i], winfo)); 365 wqkmsArray[i].setName(buildWQDescription(wqkmsArray[i], winfo));
375 } 366 }
376 367
377 /** 368 /**
378 * 369 * Builds the description label of a waterlevel
379 * @param wqkms
380 * @param descBuilder
381 * @return
382 */ 370 */
383 private String buildWQDescription(final WQKms wqkms, final WINFOArtifact winfo) { 371 private String buildWQDescription(final WQKms wqkms, final WINFOArtifact winfo) {
372
384 final WaterlevelDescriptionBuilder descBuilder = new WaterlevelDescriptionBuilder(winfo, this.context); 373 final WaterlevelDescriptionBuilder descBuilder = new WaterlevelDescriptionBuilder(winfo, this.context);
374 // TODO Zwischen numerischem Q-Wert und Dauerzahl-Hauptwert (0..364) unterscheiden
385 final String description = descBuilder.getDesc(wqkms); 375 final String description = descBuilder.getDesc(wqkms);
386 if (!description.isEmpty() && Character.isDigit(description.charAt(0))) { 376 if (!description.isEmpty() && Character.isDigit(description.charAt(0))) {
387 if (winfo.isQ()) 377 if (winfo.isQ())
388 return "Q=" + description; 378 return "Q=" + description;
389 else 379 else
392 else 382 else
393 return description; 383 return description;
394 } 384 }
395 385
396 /** 386 /**
387 * Calculates the flood durations of the Qs of the waterlevels/discharge states for a map of gauges
388 */
389 private void calcGaugeWstDurations(final WINFOArtifact winfo, final List<Gauge> gauges, final Map<Gauge, List<Double>> gaugeWstDurations,
390 final Map<Gauge, GaugeDurationValuesFinder> durFinders) {
391
392 final double[] gaugeKms = new double[gauges.size()];
393 for (int i = 0; i <= gauges.size() - 1; i++) {
394 gaugeKms[i] = gauges.get(i).getStation().doubleValue();
395 gaugeWstDurations.put(gauges.get(i), new ArrayList<Double>());
396 }
397 final CalculationResult wstsData = winfo.computeWaterlevelData(gaugeKms);
398 final WQKms[] wsts = (WQKms[]) wstsData.getData();
399 for (int i = 0; i <= gauges.size() - 1; i++) {
400 final GaugeDurationValuesFinder durFinder = durFinders.get(gauges.get(i));
401 for (int j = 0; j <= wsts.length - 1; j++) {
402 final double d = durFinder.getDuration(wsts[j].getQ(i));
403 gaugeWstDurations.get(gauges.get(i)).add(Double.valueOf(underflowDaysToOverflowDays(d)));
404 }
405 }
406 }
407
408 /**
397 * Create a result row for a station and its gauge, and add w-q-values as selected 409 * Create a result row for a station and its gauge, and add w-q-values as selected
398 */ 410 */
399 private ResultRow createRow(final Double station, final Gauge gauge, final Gauge firstGauge, final WQKms[] wqkmsArray, 411 private ResultRow createRow(final Double station, final Gauge gauge, final Gauge firstGauge, final WQKms[] wqkmsArray,
400 final GaugeDurationValuesFinder durationFinder, final int kmIndex) { 412 final List<Double> gaugeDurations, final int kmIndex) {
401 413
402 final ResultRow row = ResultRow.create(); 414 final ResultRow row = ResultRow.create();
403 row.putValue(GeneralResultType.station, station); 415 row.putValue(GeneralResultType.station, station);
404 row.putValue(SInfoResultType.infrastructuretype, null); // is replaced later for an infrastructure 416 row.putValue(SInfoResultType.infrastructuretype, null); // is replaced later for an infrastructure
405 row.putValue(SInfoResultType.floodDuration, Double.NaN); // is replaced later for an infrastructure 417 row.putValue(SInfoResultType.floodDuration, Double.NaN); // is replaced later for an infrastructure
408 row.putValue(SInfoResultType.gaugeLabel, gaugeLabel); 420 row.putValue(SInfoResultType.gaugeLabel, gaugeLabel);
409 421
410 final String location = this.riverInfoProvider.getLocation(station); 422 final String location = this.riverInfoProvider.getLocation(station);
411 row.putValue(SInfoResultType.location, location); 423 row.putValue(SInfoResultType.location, location);
412 424
413 final List<DurationWaterlevel> waterlevels = new ArrayList<>(wqkmsArray.length); 425 final List<DurationWaterlevel> wsts = new ArrayList<>(wqkmsArray.length);
414 426
415 for (final WQKms wqKms : wqkmsArray) { 427 for (int i = 0; i <= wqkmsArray.length - 1; i++) {
416 assert (wqKms.getKm(kmIndex) == station.doubleValue()); 428 assert (wqkmsArray[i].getKm(kmIndex) == station.doubleValue());
417 429
418 final int overflowDays = (int) Math.round(underflowDaysToOverflowDays(durationFinder.getDuration(wqKms.getQ(kmIndex)))); 430 final int overflowDays = (int) Math.round(gaugeDurations.get(i));
419 431
420 final DurationWaterlevel dw = new DurationWaterlevel(wqKms.getW(kmIndex), overflowDays, wqKms.getQ(kmIndex), wqKms.getName()); 432 final DurationWaterlevel dw = new DurationWaterlevel(wqkmsArray[i].getW(kmIndex), overflowDays, wqkmsArray[i].getQ(kmIndex),
421 waterlevels.add(dw); 433 wqkmsArray[i].getName());
422 } 434 wsts.add(dw);
423 row.putValue(SInfoResultType.customMultiRowColWaterlevel, waterlevels); 435 }
436 row.putValue(SInfoResultType.customMultiRowColWaterlevel, wsts);
424 437
425 return row; 438 return row;
426 } 439 }
427
428 /// **
429 // * Calculate the result row fields for one infrastructure - gives wrong duration numbers where Q changes within the
430 /// gauge range
431 // */
432 // private void calculateInfrastructure(final ResultRow row, final Gauge gauge, final InfrastructureValue
433 /// infrastructure, final WQBaseTableFinder wqFinder,
434 // final Map<Gauge, GaugeDurationValuesFinder> durFinders) {
435 //
436 // final double q = wqFinder.getDischarge(infrastructure.getStation(), infrastructure.getHeight());
437 // final double qOut = Double.isInfinite(q) ? Double.NaN : q;
438 // final double dur = underflowDaysToOverflowDays(durFinders.get(gauge).getDuration(q));
439 // row.putValue(SInfoResultType.riverside, infrastructure.getAttributeKey());
440 // row.putValue(SInfoResultType.floodDuration, dur);
441 // row.putValue(SInfoResultType.floodDischarge, qOut);
442 // row.putValue(SInfoResultType.infrastructureHeight, infrastructure.getHeight());
443 // row.putValue(SInfoResultType.infrastructuretype, infrastructure.getInfrastructure().getType().getName());
444 // }
445 440
446 /** 441 /**
447 * Calculate the result row fields for one infrastructure 442 * Calculate the result row fields for one infrastructure
448 */ 443 */
449 private void calculateInfrastructure(final ResultRow row, final Gauge gauge, final InfrastructureValue infrastructure, final WstValueTable wst, 444 private void calculateInfrastructure(final ResultRow row, final Gauge gauge, final InfrastructureValue infrastructure, final WstValueTable wst,
452 // Interpolate the infrastructure height in the wst table to get the corresponding Q 447 // Interpolate the infrastructure height in the wst table to get the corresponding Q
453 final Calculation problems = new Calculation(); 448 final Calculation problems = new Calculation();
454 final double[] qs = wst.findQsForW(infrastructure.getStation().doubleValue(), infrastructure.getHeight().doubleValue(), problems); 449 final double[] qs = wst.findQsForW(infrastructure.getStation().doubleValue(), infrastructure.getHeight().doubleValue(), problems);
455 // TODO Fehlerbehandlung (kein Q gefunden) 450 // TODO Fehlerbehandlung (kein Q gefunden)
456 final double q = (qs.length >= 1) ? qs[0] : Double.NaN; 451 final double q = (qs.length >= 1) ? qs[0] : Double.NaN;
457 final double qOut = Double.isInfinite(q) ? Double.NaN : q; 452 // Determine the relative column position of the Q of the infrastructure height
458 // Determine the relative column position of the Q
459 final QPosition qPos = wst.getQPosition(infrastructure.getStation().doubleValue(), q); 453 final QPosition qPos = wst.getQPosition(infrastructure.getStation().doubleValue(), q);
460 // Get the Q for the found column position for the station of the gauge 454 // Get the Q for the found column position for the station of the gauge
461 final double qGauge = wst.getQ(qPos, gauge.getStation().doubleValue()); 455 final double qGauge = wst.getQ(qPos, gauge.getStation().doubleValue());
462 // Interpolate the Q-D-table of the gauge 456 // Interpolate the Q-D-table of the gauge
463 final double dur = underflowDaysToOverflowDays(durFinders.get(gauge).getDuration(qGauge)); 457 final double dur = underflowDaysToOverflowDays(durFinders.get(gauge).getDuration(qGauge));
464 // Set result row 458 // Set the result row
465 row.putValue(SInfoResultType.riverside, infrastructure.getAttributeKey()); 459 row.putValue(SInfoResultType.riverside, infrastructure.getAttributeKey());
466 row.putValue(SInfoResultType.floodDuration, dur); 460 row.putValue(SInfoResultType.floodDuration, dur);
467 row.putValue(SInfoResultType.floodDischarge, qOut); 461 row.putValue(SInfoResultType.floodDischarge, q);
468 row.putValue(SInfoResultType.infrastructureHeight, infrastructure.getHeight()); 462 row.putValue(SInfoResultType.infrastructureHeight, infrastructure.getHeight());
469 row.putValue(SInfoResultType.infrastructuretype, infrastructure.getInfrastructure().getType().getName()); 463 row.putValue(SInfoResultType.infrastructuretype, infrastructure.getInfrastructure().getType().getName());
470 } 464 }
471 465
472 /** 466 /**
http://dive4elements.wald.intevation.org