Mercurial > dive4elements > river

comparison artifacts/src/main/java/org/dive4elements/river/artifacts/sinfo/flowdepth/FlowDepthCalculation.java @ 8898:89f3c5462a16

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Implemented S-INFO Flowdepth TKH calculation
author mschaefer
date Thu, 22 Feb 2018 12:03:31 +0100
parents f431aec10d2c
children d32c22fc686c
comparison
equal deleted inserted replaced
8891:f431aec10d2c 8898:89f3c5462a16
14 import java.util.Collection; 14 import java.util.Collection;
15 import java.util.Collections; 15 import java.util.Collections;
16 import java.util.Date; 16 import java.util.Date;
17 import java.util.List; 17 import java.util.List;
18 18
19 import org.apache.commons.lang.math.DoubleRange;
19 import org.apache.commons.math.FunctionEvaluationException; 20 import org.apache.commons.math.FunctionEvaluationException;
20 import org.apache.commons.math.analysis.UnivariateRealFunction; 21 import org.apache.commons.math.analysis.UnivariateRealFunction;
22 import org.apache.log4j.Logger;
21 import org.dive4elements.artifacts.ArtifactDatabase; 23 import org.dive4elements.artifacts.ArtifactDatabase;
22 import org.dive4elements.artifacts.CallContext; 24 import org.dive4elements.artifacts.CallContext;
23 import org.dive4elements.river.artifacts.BedHeightsArtifact; 25 import org.dive4elements.river.artifacts.BedHeightsArtifact;
24 import org.dive4elements.river.artifacts.model.Calculation; 26 import org.dive4elements.river.artifacts.model.Calculation;
25 import org.dive4elements.river.artifacts.model.CalculationResult; 27 import org.dive4elements.river.artifacts.model.CalculationResult;
28 import org.dive4elements.river.artifacts.model.DateRange;
26 import org.dive4elements.river.artifacts.model.LocationProvider; 29 import org.dive4elements.river.artifacts.model.LocationProvider;
27 import org.dive4elements.river.artifacts.model.QKms; 30 import org.dive4elements.river.artifacts.model.WQKms;
28 import org.dive4elements.river.artifacts.model.WKms; 31 import org.dive4elements.river.artifacts.model.WKms;
29 import org.dive4elements.river.artifacts.model.minfo.QualityMeasurementFactory;
30 import org.dive4elements.river.artifacts.model.minfo.QualityMeasurements;
31 import org.dive4elements.river.artifacts.resources.Resources; 32 import org.dive4elements.river.artifacts.resources.Resources;
32 import org.dive4elements.river.artifacts.sinfo.SINFOArtifact; 33 import org.dive4elements.river.artifacts.sinfo.SINFOArtifact;
33 import org.dive4elements.river.artifacts.sinfo.flowdepth.FlowDepthAccess.DifferencesPair; 34 import org.dive4elements.river.artifacts.sinfo.flowdepth.FlowDepthAccess.DifferencesPair;
34 import org.dive4elements.river.artifacts.states.WaterlevelData; 35 import org.dive4elements.river.artifacts.states.WaterlevelData;
35 import org.dive4elements.river.artifacts.states.WaterlevelFetcher; 36 import org.dive4elements.river.artifacts.states.WaterlevelFetcher;
41 import org.dive4elements.river.utils.GaugeIndex; 42 import org.dive4elements.river.utils.GaugeIndex;
42 import org.dive4elements.river.utils.RiverUtils; 43 import org.dive4elements.river.utils.RiverUtils;
43 44
44 class FlowDepthCalculation { 45 class FlowDepthCalculation {
45 46
47 private static Logger log = Logger.getLogger(FlowDepthCalculation.class);
48
46 private static final int VALID_BED_MEASUREMENT_YEARS = 20; 49 private static final int VALID_BED_MEASUREMENT_YEARS = 20;
47 50
48 private static final String CSV_NOT_IN_GAUGE_RANGE = "export.waterlevel.csv.not.in.gauge.range"; 51 private static final String CSV_NOT_IN_GAUGE_RANGE = "export.waterlevel.csv.not.in.gauge.range";
49 52
50 private final CallContext context; 53 private final CallContext context;
68 71
69 final Collection<DifferencesPair> diffPairs = access.getDifferencePairs(); 72 final Collection<DifferencesPair> diffPairs = access.getDifferencePairs();
70 73
71 final double from = access.getFrom(); 74 final double from = access.getFrom();
72 final double to = access.getTo(); 75 final double to = access.getTo();
76 final DoubleRange calcRange = new DoubleRange(from, to);
73 77
74 final boolean useTkh = access.isUseTransportBodies(); 78 final boolean useTkh = access.isUseTransportBodies();
75 79
76 /* calculate results for each diff pair */ 80 /* calculate results for each diff pair */
77 final Calculation problems = new Calculation(); 81 final Calculation problems = new Calculation();
82 final String calcModeLabel = Resources.getMsg(this.context.getMeta(), sinfo.getCalculationMode().name()); 86 final String calcModeLabel = Resources.getMsg(this.context.getMeta(), sinfo.getCalculationMode().name());
83 87
84 final FlowDepthCalculationResults results = new FlowDepthCalculationResults(calcModeLabel, user, river, from, to, useTkh); 88 final FlowDepthCalculationResults results = new FlowDepthCalculationResults(calcModeLabel, user, river, from, to, useTkh);
85 89
86 for (final DifferencesPair diffPair : diffPairs) { 90 for (final DifferencesPair diffPair : diffPairs) {
87 final FlowDepthCalculationResult result = calculateResult(river, from, to, diffPair, problems, gaugeIndex); 91 final FlowDepthCalculationResult result = calculateResult(river, calcRange, diffPair, problems, gaugeIndex, useTkh);
88 if (result != null) 92 if (result != null)
89 results.addResult(result); 93 results.addResult(result);
90 } 94 }
91 95
92 return new CalculationResult(results, problems); 96 return new CalculationResult(results, problems);
93 } 97 }
94 98
95 private FlowDepthCalculationResult calculateResult(final River river, final double from, final double to, final DifferencesPair diffPair, 99 /**
96 final Calculation problems, final GaugeIndex gaugeIndex) { 100 * Calculates one W-MSH differences pair.
101 */
102 private FlowDepthCalculationResult calculateResult(final River river, final DoubleRange calcRange, final DifferencesPair diffPair,
103 final Calculation problems, final GaugeIndex gaugeIndex, final boolean useTkh) {
97 104
98 /* access real input data from database */ 105 /* access real input data from database */
99 final String soundingId = diffPair.getSoundingId(); 106 final String soundingId = diffPair.getSoundingId();
100 final String wstId = diffPair.getWstId(); 107 final String wstId = diffPair.getWstId();
101 108
102 final BedHeight bedHeight = loadBedHeight(soundingId, from, to); 109 final BedHeight bedHeight = loadBedHeight(soundingId);
103 if (bedHeight == null) { 110 if (bedHeight == null) {
104 final String message = Resources.format(this.context.getMeta(), "Failed to access sounding with id '{0}'", soundingId); 111 final String message = Resources.format(this.context.getMeta(), "Failed to access sounding with id '{0}'", soundingId);
105 problems.addProblem(message); 112 problems.addProblem(message);
106 return null; 113 return null;
107 } 114 }
134 141
135 final FlowDepthCalculationResult resultData = new FlowDepthCalculationResult(label, wstInfo, sounding); 142 final FlowDepthCalculationResult resultData = new FlowDepthCalculationResult(label, wstInfo, sounding);
136 143
137 // FIXME: nur prüfen/beschaffen wenn TKH Berechnung aktiv 144 // FIXME: nur prüfen/beschaffen wenn TKH Berechnung aktiv
138 /* Abflusswerte vorhanden? */ 145 /* Abflusswerte vorhanden? */
139 if (!(wstKms instanceof QKms)) { 146 if (!(wstKms instanceof WQKms)) {
140 final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingQ", null, label); 147 final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingQ", null, label);
141 problems.addProblem(message); 148 problems.addProblem(message);
142 // TODO: keine Berechnung TKH 149 // TODO: keine Berechnung TKH
143 } 150 }
144 151
145 final QualityMeasurements bedMeasurements = getBedMeasurements(river, from, to, sounding.getYear()); 152 BedQualityD50KmValueFinder bedMeasurementsFinder = null;
146 // FIXME: prüfung ob (genug) werte vorhanden sind? was sind genau die kriterien? falls nein, problemhinzufügen und keine 153 if (useTkh)
154 bedMeasurementsFinder = loadBedMeasurements(river, calcRange, sounding.getYear());
155 // FIXME: prüfung ob (genug) werte vorhanden sind? was sind genau die kriterien? falls nein, problem hinzufügen und keine
147 // berechnung tkh 156 // berechnung tkh
148 // FIXME: wie wird ggf. interpoliert? --> absprache?
149 // FIXME: mir direkt aufgefallen, die Beispieldatenbank liefert Werte zum gleichen km und zeitpunkt, die messwerte sind
150 // aber unterschiedlich....???
151 // FIXME: die eigentlichen daten extrahieren, ggf. wenn esswerte zum gleichen datum voriliegen. das neueste nehmen? oder
152 // das näheste zum Peiljahr?
153
154 // FIXME Art der Gewässersohle (starr/mobil)
155 // FIXME: wie wird ggf. interpoliert? prüfung ob werte vorhanden?
156 157
157 final String bedHeightLabel = bedHeight.getDescription(); 158 final String bedHeightLabel = bedHeight.getDescription();
158 final String wstLabel = wstKms.getName(); 159 final String wstLabel = wstKms.getName();
159 160
160 final UnivariateRealFunction wstInterpolator = DoubleUtil.getLinearInterpolator(wstKms.allKms(), wstKms.allWs()); 161 final UnivariateRealFunction wstInterpolator = DoubleUtil.getLinearInterpolator(wstKms.allKms(), wstKms.allWs());
161 162 UnivariateRealFunction qInterpolator = null;
163 DoubleRange qRange = null;
164 if (useTkh && (wstKms instanceof WQKms)) {
165 qInterpolator = DoubleUtil.getLinearInterpolator(((WQKms) wstKms).allKms(), ((WQKms) wstKms).allQs());
166 qRange = new DoubleRange( ((WQKms) wstKms).allQs().min(), ((WQKms) wstKms).allQs().max());
167 }
168
162 // FIXME: sort by station first, but in what direction? 169 // FIXME: sort by station first, but in what direction?
170 // FIXME: using river.getKmUp()?
163 final List<BedHeightValue> values = bedHeight.getValues(); 171 final List<BedHeightValue> values = bedHeight.getValues();
164 172
165 final List<BedHeightValue> sortedValues = new ArrayList<>(values); 173 final List<BedHeightValue> sortedValues = new ArrayList<>(values);
166 Collections.sort(sortedValues, new BedHeightStationComparator()); 174 Collections.sort(sortedValues, new BedHeightStationComparator());
167 175
168 SoilKind lastKind = SoilKind.mobil; 176 // FIXME: wie wird ggf. interpoliert? prüfung ob werte vorhanden?
177 /* SoilKind lastKind = SoilKind.mobil; */
178 SoilKindKmValueFinder soilKindFinder = null;
179 if (useTkh) {
180 soilKindFinder = new SoilKindKmValueFinder();
181 soilKindFinder.loadValues(river, calcRange);
182 }
183
184 FlowVelocityModelKmValueFinder flowVelocitiesFinder = null;
185 if (useTkh) {
186 flowVelocitiesFinder = new FlowVelocityModelKmValueFinder();
187 flowVelocitiesFinder.loadValues(river, calcRange, qRange);
188 }
169 189
170 for (final BedHeightValue bedHeightValue : sortedValues) { 190 for (final BedHeightValue bedHeightValue : sortedValues) {
171 191
172 final Double station = bedHeightValue.getStation(); 192 final Double station = bedHeightValue.getStation();
173 if (station == null || station.isNaN()) 193 if (station == null || station.isNaN())
178 continue; 198 continue;
179 199
180 final double km = station; 200 final double km = station;
181 final double meanBedHeight = meanBedHeightDbl; 201 final double meanBedHeight = meanBedHeightDbl;
182 202
203 if (!calcRange.containsDouble(km))
204 continue;
205
183 try { 206 try {
184 // FIXME: check out of range 207 // FIXME: check out of range
185 final double wst = wstInterpolator.value(km); 208 final double wst = wstInterpolator.value(km);
186 209
187 final double flowDepth = wst - meanBedHeight; 210 final double flowDepth = wst - meanBedHeight;
188 211
189 // FIXME: piecewise constant interpolation? 212 // FIXME: piecewise constant interpolation?
190 // final double discharge = wstKms instanceof QKms ? ((QKms) wstKms).getQ(i) : Double.NaN; 213 // final double discharge = wstKms instanceof QKms ? ((QKms) wstKms).getQ(i) : Double.NaN;
191 final double discharge = Double.NaN; 214 double discharge;
215 if (qInterpolator != null)
216 discharge = qInterpolator.value(km);
217 else
218 discharge = Double.NaN;
192 219
193 // FIXME: calculate tkh 220 // FIXME: calculate tkh
194 221 double tkh = 0;
195 // REMARK: bissl spielerei zum testen damit die sohlart nicht zu schnell wechselt 222 if (useTkh) {
223 double d50 = 0;
224 try {
225 d50 = bedMeasurementsFinder.findD50(km);
226 } catch (Exception e) {
227 final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingD50", null, label);
228 problems.addProblem(km, message);
229 //FIXME: cumulate problems
230 }
231 if (flowVelocitiesFinder.findKmQValues(km, discharge)) {
232 tkh = calculateTkh(wst - meanBedHeight, flowVelocitiesFinder.getFindVmainFound(), d50, flowVelocitiesFinder.getFindTauFound());
233 log.debug(String.format("calculateTkh km %.3f q %.0f w %.2f mbh %.2f vm %.1f tau %.1f d50(mm) %.1f tkh(cm) %.1f",
234 km, discharge, wst, meanBedHeight, flowVelocitiesFinder.getFindVmainFound(), flowVelocitiesFinder.getFindTauFound(), d50*1000, tkh));
235 }
236 else
237 tkh = Double.NaN;
238 }
239
240 // Soil kind
241 SoilKind kind = SoilKind.starr;
242 if (useTkh) {
243 try {
244 kind = soilKindFinder.findSoilKind(km);
245 } catch (Exception e) {
246 final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingSoilKind", null, label);
247 problems.addProblem(km, message);
248 //FIXME: cumulate problems
249 }
250 }
251
252 /* // REMARK: bissl spielerei zum testen damit die sohlart nicht zu schnell wechselt
253 SoilKind kind;
196 final boolean changeKind = Math.random() > 0.95; 254 final boolean changeKind = Math.random() > 0.95;
197 SoilKind kind;
198 if (changeKind) { 255 if (changeKind) {
199 switch (lastKind) { 256 switch (lastKind) {
200 case starr: 257 case starr:
201 kind = SoilKind.mobil; 258 kind = SoilKind.mobil;
202 break; 259 break;
203
204 case mobil: 260 case mobil:
205 default: 261 default:
206 kind = SoilKind.starr; 262 kind = SoilKind.starr;
207 break; 263 break;
208
209 } 264 }
210 } else 265 } else
211 kind = lastKind; 266 kind = lastKind;
212 lastKind = kind; 267 lastKind = kind;
213 268 */
214 final double tkh = 100 + 10 * (Math.random() - 0.5); 269
270 //tkh = 100 + 10 * (Math.random() - 0.5);
215 271
216 final double flowDepthTkh; 272 final double flowDepthTkh;
217 final double tkhUp; 273 final double tkhUp;
218 final double tkhDown; 274 final double tkhDown;
219 switch (kind) { 275 switch (kind) {
238 // REMARK: access the gauge once only during calculation 294 // REMARK: access the gauge once only during calculation
239 final Gauge gauge = findGauge(waterlevel, refGauge, gaugeIndex, km); 295 final Gauge gauge = findGauge(waterlevel, refGauge, gaugeIndex, km);
240 296
241 final String gaugeLabel = gauge == null ? notinrange : gauge.getName(); 297 final String gaugeLabel = gauge == null ? notinrange : gauge.getName();
242 298
243 resultData.addRow(km, flowDepth, flowDepthTkh, kind, tkh, tkhUp, tkhDown, wst, discharge, wstLabel, gaugeLabel, meanBedHeight, bedHeightLabel, 299 resultData.addRow(km, flowDepth, flowDepthTkh, kind, tkh, tkhUp, tkhDown, wst, discharge, wstLabel,
244 location); 300 gaugeLabel, meanBedHeight, bedHeightLabel, location);
245 } 301 }
246 catch (final FunctionEvaluationException e) { 302 catch (final FunctionEvaluationException e) {
247 /* should only happen if out of range */ 303 /* should only happen if out of range */
248 e.printStackTrace(); 304 e.printStackTrace();
249 /* simply ignore */ 305 /* simply ignore */
256 312
257 /** 313 /**
258 * Sohlbeschaffenheit (D50 Korndurchmesser aus Seddb) 314 * Sohlbeschaffenheit (D50 Korndurchmesser aus Seddb)
259 * Abhängig von Peiljahr 315 * Abhängig von Peiljahr
260 */ 316 */
261 private QualityMeasurements getBedMeasurements(final River river, final double from, final double to, final int soundingYear) { 317 private BedQualityD50KmValueFinder loadBedMeasurements(final River river, final DoubleRange kmRange, final int soundingYear) {
262 318
263 /* construct valid measurement time range */ 319 /* construct valid measurement time range */
264 final Calendar cal = Calendar.getInstance(); 320 final Calendar cal = Calendar.getInstance();
265 cal.clear(); 321 cal.clear();
266 322
268 final Date startTime = cal.getTime(); 324 final Date startTime = cal.getTime();
269 325
270 cal.set(soundingYear + VALID_BED_MEASUREMENT_YEARS, 11, 31); 326 cal.set(soundingYear + VALID_BED_MEASUREMENT_YEARS, 11, 31);
271 final Date endTime = cal.getTime(); 327 final Date endTime = cal.getTime();
272 328
273 return QualityMeasurementFactory.getBedMeasurements(river.getName(), from, to, startTime, endTime); 329 final BedQualityD50KmValueFinder finder = new BedQualityD50KmValueFinder();
330 if (finder.loadValues(river, kmRange, new DateRange(startTime, endTime)))
331 return finder;
332 else
333 return null;
274 } 334 }
275 335
276 /** 336 /**
277 * Checks the year difference between waterlevels and sounding, and issues a warning if too big. 337 * Checks the year difference between waterlevels and sounding, and issues a warning if too big.
278 * 338 *
346 problems.addProblem(kmPrev, message); 406 problems.addProblem(kmPrev, message);
347 } 407 }
348 } 408 }
349 } 409 }
350 410
351 private BedHeight loadBedHeight(final String soundingId, final double from, final double to) { 411 private BedHeight loadBedHeight(final String soundingId) {
352 412
353 // REMARK: absolutely unbelievable.... 413 // REMARK: absolutely unbelievable....
354 // The way how bed-heights (and other data too) is accessed is different for nearly ever calculation-type 414 // The way how bed-heights (and other data too) is accessed is different for nearly every calculation-type
355 // throughout flys. 415 // throughout flys.
356 // The knowledge on how to parse the datacage-ids is spread through the complete code-base... 416 // The knowledge on how to parse the datacage-ids is spread through the complete code-base...
357 417
358 // We use here the way on how bed-heights are accessed by the BedDifferenceAccess/BedDifferenceCalculation, but 418 // We use here the way on how bed-heights are accessed by the BedDifferenceAccess/BedDifferenceCalculation, but
359 // this is plain random 419 // this is plain random
373 // BedHeightFactory uses its own (direct) way of accessing the data, with its own implemented data classes. 433 // BedHeightFactory uses its own (direct) way of accessing the data, with its own implemented data classes.
374 // return BedHeightFactory.getHeight(bedheightType, bedheightId, from, to); 434 // return BedHeightFactory.getHeight(bedheightType, bedheightId, from, to);
375 435
376 return BedHeight.getBedHeightById(bedheightId); 436 return BedHeight.getBedHeightById(bedheightId);
377 } 437 }
438
439 /**
440 * Calculates a transport body height
441 * @param h flow depth in m
442 * @param vm flow velocity in m
443 * @param d50 grain diameter D50 in m (!)
444 * @param tau shear stress in N/m^2
445 * @return transport body height in cm (!)
446 */
447 private double calculateTkh(double h, double vm, double d50, double tau) {
448 final double PHYS_G = 9.81;
449 final double PHYS_SPECGRAV_S = 2.6;
450 final double PHYS_VELOCCOEFF_N = 6;
451 final double PHYS_FORMCOEFF_ALPHA = 0.7;
452 final double PHYS_VISCOSITY_NUE = 1.3e-6;
453 final double PHYS_GRAIN_DENSITY_RHOS = 2603;
454 final double PHYS_WATER_DENSITY_RHO = 999.97;
455
456 final double froude = vm / Math.sqrt(PHYS_G * h);
457 final double partReynolds = Math.sqrt((PHYS_SPECGRAV_S - 1) * PHYS_G * d50) / PHYS_VISCOSITY_NUE * d50;
458 final double critShields = 0.22 * Math.pow(partReynolds, -0.6) + 0.06 * Math.pow(10, 7.7 * Math.pow(partReynolds, -0.6));
459 final double critTau = critShields * (PHYS_GRAIN_DENSITY_RHOS - PHYS_WATER_DENSITY_RHO) * PHYS_G * d50;
460 return 100 * h * (1 - Math.pow(froude, 2)) / (2 * PHYS_VELOCCOEFF_N * PHYS_FORMCOEFF_ALPHA) * (1 - critTau / tau);
461 }
378 } 462 }
http://dive4elements.wald.intevation.org