gernotbelger@8854: /* Copyright (C) 2017 by Bundesanstalt für Gewässerkunde gernotbelger@8877: * Software engineering by gernotbelger@8877: * Björnsen Beratende Ingenieure GmbH gernotbelger@8854: * Dr. Schumacher Ingenieurbüro für Wasser und Umwelt gernotbelger@8854: * gernotbelger@8854: * This file is Free Software under the GNU AGPL (>=v3) gernotbelger@8854: * and comes with ABSOLUTELY NO WARRANTY! Check out the gernotbelger@8854: * documentation coming with Dive4Elements River for details. gernotbelger@8854: */ gernotbelger@8854: package org.dive4elements.river.artifacts.sinfo.flowdepth; gernotbelger@8854: gernotbelger@8883: import java.util.ArrayList; gernotbelger@8891: import java.util.Calendar; gernotbelger@8854: import java.util.Collection; gernotbelger@8883: import java.util.Collections; gernotbelger@8891: import java.util.Date; gernotbelger@8854: import java.util.List; gernotbelger@8854: gernotbelger@8894: import org.apache.commons.lang.math.DoubleRange; gernotbelger@8883: import org.apache.commons.math.FunctionEvaluationException; gernotbelger@8883: import org.apache.commons.math.analysis.UnivariateRealFunction; mschaefer@8898: import org.apache.log4j.Logger; gernotbelger@8879: import org.dive4elements.artifacts.ArtifactDatabase; gernotbelger@8854: import org.dive4elements.artifacts.CallContext; gernotbelger@8854: import org.dive4elements.river.artifacts.BedHeightsArtifact; gernotbelger@8854: import org.dive4elements.river.artifacts.model.Calculation; gernotbelger@8854: import org.dive4elements.river.artifacts.model.CalculationResult; mschaefer@8898: import org.dive4elements.river.artifacts.model.DateRange; gernotbelger@8854: import org.dive4elements.river.artifacts.model.LocationProvider; mschaefer@8901: import org.dive4elements.river.artifacts.model.QKms; gernotbelger@8854: import org.dive4elements.river.artifacts.model.WKms; gernotbelger@8854: import org.dive4elements.river.artifacts.resources.Resources; gernotbelger@8854: import org.dive4elements.river.artifacts.sinfo.SINFOArtifact; gernotbelger@8854: import org.dive4elements.river.artifacts.sinfo.flowdepth.FlowDepthAccess.DifferencesPair; gernotbelger@8894: import org.dive4elements.river.artifacts.sinfo.util.BedHeightInfo; gernotbelger@8894: import org.dive4elements.river.artifacts.sinfo.util.RiverInfo; gernotbelger@8894: import org.dive4elements.river.artifacts.sinfo.util.WstInfo; gernotbelger@8882: import org.dive4elements.river.artifacts.states.WaterlevelData; gernotbelger@8882: import org.dive4elements.river.artifacts.states.WaterlevelFetcher; gernotbelger@8854: import org.dive4elements.river.model.BedHeight; gernotbelger@8883: import org.dive4elements.river.model.BedHeightValue; gernotbelger@8854: import org.dive4elements.river.model.Gauge; gernotbelger@8854: import org.dive4elements.river.model.River; gernotbelger@8883: import org.dive4elements.river.utils.DoubleUtil; gernotbelger@8854: import org.dive4elements.river.utils.GaugeIndex; gernotbelger@8854: import org.dive4elements.river.utils.RiverUtils; gernotbelger@8854: gernotbelger@8854: class FlowDepthCalculation { gernotbelger@8854: mschaefer@8898: private static Logger log = Logger.getLogger(FlowDepthCalculation.class); mschaefer@8898: gernotbelger@8891: private static final int VALID_BED_MEASUREMENT_YEARS = 20; gernotbelger@8891: gernotbelger@8854: private static final String CSV_NOT_IN_GAUGE_RANGE = "export.waterlevel.csv.not.in.gauge.range"; gernotbelger@8854: gernotbelger@8877: private final CallContext context; gernotbelger@8854: gernotbelger@8882: public FlowDepthCalculation(final CallContext context) { gernotbelger@8877: this.context = context; gernotbelger@8854: } gernotbelger@8854: gernotbelger@8877: public CalculationResult calculate(final SINFOArtifact sinfo) { gernotbelger@8854: gernotbelger@8879: /* gernotbelger@8879: * find the user of this artifact, sadly this is not part of the calling context, so instead we determine the gernotbelger@8879: * owner oft the artifact gernotbelger@8879: */ gernotbelger@8879: final ArtifactDatabase database = this.context.getDatabase(); gernotbelger@8879: final String user = database.findArtifactUser(sinfo.identifier()); gernotbelger@8854: gernotbelger@8877: /* access input data */ gernotbelger@8877: final FlowDepthAccess access = new FlowDepthAccess(sinfo); gernotbelger@8877: final River river = access.getRiver(); gernotbelger@8894: final RiverInfo riverInfo = new RiverInfo(river); gernotbelger@8854: gernotbelger@8877: final Collection diffPairs = access.getDifferencePairs(); gernotbelger@8877: gernotbelger@8877: final double from = access.getFrom(); gernotbelger@8877: final double to = access.getTo(); gernotbelger@8894: final DoubleRange calcRange = new DoubleRange(from, to); gernotbelger@8877: gernotbelger@8877: final boolean useTkh = access.isUseTransportBodies(); gernotbelger@8877: gernotbelger@8877: /* calculate results for each diff pair */ gernotbelger@8877: final Calculation problems = new Calculation(); gernotbelger@8877: gernotbelger@8877: final List gauges = river.determineGauges(from, to); gernotbelger@8877: final GaugeIndex gaugeIndex = new GaugeIndex(gauges); gernotbelger@8877: gernotbelger@8882: final String calcModeLabel = Resources.getMsg(this.context.getMeta(), sinfo.getCalculationMode().name()); gernotbelger@8877: gernotbelger@8894: final FlowDepthCalculationResults results = new FlowDepthCalculationResults(calcModeLabel, user, riverInfo, calcRange, useTkh); gernotbelger@8877: gernotbelger@8877: for (final DifferencesPair diffPair : diffPairs) { mschaefer@8898: final FlowDepthCalculationResult result = calculateResult(river, calcRange, diffPair, problems, gaugeIndex, useTkh); gernotbelger@8882: if (result != null) gernotbelger@8877: results.addResult(result); gernotbelger@8877: } gernotbelger@8877: gernotbelger@8882: return new CalculationResult(results, problems); gernotbelger@8877: } gernotbelger@8877: mschaefer@8898: /** mschaefer@8898: * Calculates one W-MSH differences pair. mschaefer@8898: */ gernotbelger@8894: private FlowDepthCalculationResult calculateResult(final River river, final DoubleRange calcRange, final DifferencesPair diffPair, mschaefer@8898: final Calculation problems, final GaugeIndex gaugeIndex, final boolean useTkh) { gernotbelger@8877: gernotbelger@8877: /* access real input data from database */ gernotbelger@8877: final String soundingId = diffPair.getSoundingId(); gernotbelger@8877: final String wstId = diffPair.getWstId(); gernotbelger@8877: gernotbelger@8894: final BedHeight bedHeight = loadBedHeight(soundingId); gernotbelger@8882: if (bedHeight == null) { gernotbelger@8884: final String message = Resources.format(this.context.getMeta(), "Failed to access sounding with id '{0}'", soundingId); gernotbelger@8877: problems.addProblem(message); gernotbelger@8877: return null; gernotbelger@8877: } gernotbelger@8877: gernotbelger@8882: /* REMARK: fetch ALL wst kms, because we want to determine the original reference gauge */ gernotbelger@8884: final WaterlevelData waterlevel = new WaterlevelFetcher().findWaterlevel(this.context, wstId, Double.NaN, Double.NaN); gernotbelger@8882: if (waterlevel == null) { gernotbelger@8884: final String message = Resources.format(this.context.getMeta(), "Failed to access waterlevel with id '{0}'", wstId); gernotbelger@8877: problems.addProblem(message); gernotbelger@8877: return null; gernotbelger@8877: } gernotbelger@8882: final WKms wstKms = waterlevel.getWkms(); gernotbelger@8882: gernotbelger@8883: final String wspLabel = wstKms.getName(); gernotbelger@8883: final String soundingLabel = bedHeight.getDescription(); gernotbelger@8883: final String label = String.format("%s - %s", wspLabel, soundingLabel); gernotbelger@8877: gernotbelger@8883: checkYearDifference(label, waterlevel, bedHeight, problems); gernotbelger@8894: checkWaterlevelDiscretisation(wstKms, calcRange, problems); mschaefer@8901: // TODO: prüfen, ob sohlhöhen die calcRange abdecken/überschneiden gernotbelger@8882: gernotbelger@8882: /* re-determine the reference gauge, in the same way as the WaterlevelArtifact would do it */ gernotbelger@8884: final String notinrange = Resources.getMsg(this.context.getMeta(), CSV_NOT_IN_GAUGE_RANGE, CSV_NOT_IN_GAUGE_RANGE); gernotbelger@8882: gernotbelger@8882: final Gauge refGauge = waterlevel.findReferenceGauge(river); gernotbelger@8882: final String refGaugeName = refGauge == null ? notinrange : refGauge.getName(); gernotbelger@8877: gernotbelger@8877: final BedHeightInfo sounding = BedHeightInfo.from(bedHeight); gernotbelger@8883: final int wspYear = waterlevel.getYear(); gernotbelger@8882: final WstInfo wstInfo = new WstInfo(wspLabel, wspYear, refGaugeName); gernotbelger@8877: gernotbelger@8877: final FlowDepthCalculationResult resultData = new FlowDepthCalculationResult(label, wstInfo, sounding); gernotbelger@8877: mschaefer@8901: boolean doCalcTkh = useTkh; mschaefer@8901: if (doCalcTkh && !(wstKms instanceof QKms)) { gernotbelger@8884: final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingQ", null, label); gernotbelger@8877: problems.addProblem(message); mschaefer@8901: doCalcTkh = false; gernotbelger@8877: } gernotbelger@8877: mschaefer@8898: BedQualityD50KmValueFinder bedMeasurementsFinder = null; mschaefer@8911: if (doCalcTkh) { mschaefer@8901: bedMeasurementsFinder = loadBedMeasurements(river, calcRange, sounding.getYear().intValue()); mschaefer@8911: if (bedMeasurementsFinder == null) { mschaefer@8911: final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingSoilKind", null, label); mschaefer@8911: problems.addProblem(message); mschaefer@8911: doCalcTkh = false; mschaefer@8911: } mschaefer@8911: } gernotbelger@8877: gernotbelger@8877: final String bedHeightLabel = bedHeight.getDescription(); gernotbelger@8877: final String wstLabel = wstKms.getName(); gernotbelger@8877: gernotbelger@8884: final UnivariateRealFunction wstInterpolator = DoubleUtil.getLinearInterpolator(wstKms.allKms(), wstKms.allWs()); mschaefer@8898: UnivariateRealFunction qInterpolator = null; mschaefer@8898: DoubleRange qRange = null; mschaefer@8901: if (doCalcTkh) { mschaefer@8901: qInterpolator = DoubleUtil.getLinearInterpolator(((QKms) wstKms).allKms(), ((QKms) wstKms).allQs()); mschaefer@8911: if (qInterpolator != null) mschaefer@8911: qRange = new DoubleRange( ((QKms) wstKms).allQs().min(), ((QKms) wstKms).allQs().max()); mschaefer@8911: else { mschaefer@8911: final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingQ", null, label); mschaefer@8911: problems.addProblem(message); mschaefer@8911: doCalcTkh = false; mschaefer@8911: } mschaefer@8898: } mschaefer@8898: gernotbelger@8883: // FIXME: sort by station first, but in what direction? mschaefer@8898: // FIXME: using river.getKmUp()? gernotbelger@8883: final List values = bedHeight.getValues(); gernotbelger@8877: gernotbelger@8883: final List sortedValues = new ArrayList<>(values); gernotbelger@8883: Collections.sort(sortedValues, new BedHeightStationComparator()); gernotbelger@8882: mschaefer@8898: // FIXME: wie wird ggf. interpoliert? prüfung ob werte vorhanden? mschaefer@8898: /* SoilKind lastKind = SoilKind.mobil; */ mschaefer@8898: SoilKindKmValueFinder soilKindFinder = null; mschaefer@8901: if (doCalcTkh) { mschaefer@8898: soilKindFinder = new SoilKindKmValueFinder(); mschaefer@8911: if (!soilKindFinder.loadValues(river, calcRange)) { mschaefer@8911: doCalcTkh = false; mschaefer@8911: final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingSoilKind", null, label); mschaefer@8911: problems.addProblem(message); mschaefer@8911: } mschaefer@8898: } mschaefer@8898: mschaefer@8898: FlowVelocityModelKmValueFinder flowVelocitiesFinder = null; mschaefer@8901: if (doCalcTkh) { mschaefer@8898: flowVelocitiesFinder = new FlowVelocityModelKmValueFinder(); mschaefer@8911: if (!flowVelocitiesFinder.loadValues(river, calcRange, qRange)) { mschaefer@8911: doCalcTkh = false; mschaefer@8911: final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingVelocity", null, label); mschaefer@8911: problems.addProblem(message); mschaefer@8911: } mschaefer@8898: } gernotbelger@8886: gernotbelger@8883: for (final BedHeightValue bedHeightValue : sortedValues) { gernotbelger@8877: gernotbelger@8883: final Double station = bedHeightValue.getStation(); gernotbelger@8883: if (station == null || station.isNaN()) gernotbelger@8883: continue; gernotbelger@8883: gernotbelger@8883: final Double meanBedHeightDbl = bedHeightValue.getHeight(); gernotbelger@8883: if (meanBedHeightDbl == null || meanBedHeightDbl.isNaN()) gernotbelger@8883: continue; gernotbelger@8883: gernotbelger@8883: final double km = station; gernotbelger@8883: final double meanBedHeight = meanBedHeightDbl; gernotbelger@8883: gernotbelger@8894: if (!calcRange.containsDouble(km)) gernotbelger@8894: continue; gernotbelger@8894: gernotbelger@8883: try { gernotbelger@8883: // FIXME: check out of range gernotbelger@8883: final double wst = wstInterpolator.value(km); gernotbelger@8883: gernotbelger@8883: final double flowDepth = wst - meanBedHeight; gernotbelger@8883: gernotbelger@8883: // FIXME: piecewise constant interpolation? gernotbelger@8883: // final double discharge = wstKms instanceof QKms ? ((QKms) wstKms).getQ(i) : Double.NaN; mschaefer@8898: double discharge; mschaefer@8898: if (qInterpolator != null) mschaefer@8898: discharge = qInterpolator.value(km); mschaefer@8898: else mschaefer@8898: discharge = Double.NaN; gernotbelger@8883: mschaefer@8911: // Calculate tkh mschaefer@8898: double tkh = 0; mschaefer@8901: if (doCalcTkh) { mschaefer@8911: double d50 = Double.NaN; mschaefer@8911: try { mschaefer@8911: d50 = bedMeasurementsFinder.findD50(km); mschaefer@8911: } catch (Exception e) { mschaefer@8898: final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingD50", null, label); mschaefer@8898: problems.addProblem(km, message); mschaefer@8901: //FIXME: cumulate problems to one message? mschaefer@8898: } mschaefer@8901: if (!Double.isNaN(d50)) { mschaefer@8901: if (flowVelocitiesFinder.findKmQValues(km, discharge)) { mschaefer@8901: tkh = calculateTkh(wst - meanBedHeight, flowVelocitiesFinder.getFindVmainFound(), d50, flowVelocitiesFinder.getFindTauFound()); mschaefer@8911: if (!Double.isNaN(tkh) && (tkh < 0)) mschaefer@8911: tkh = 0; mschaefer@8901: /* log.debug(String.format("calculateTkh km %.3f q %.0f w %.2f mbh %.2f vm %.1f tau %.1f d50(mm) %.1f tkh(cm) %.1f", mschaefer@8901: km, discharge, wst, meanBedHeight, flowVelocitiesFinder.getFindVmainFound(), flowVelocitiesFinder.getFindTauFound(), d50*1000, tkh)); */ mschaefer@8901: } mschaefer@8901: else { mschaefer@8901: final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingQ", null, label); mschaefer@8901: problems.addProblem(km, message); mschaefer@8901: //FIXME: cumulate problems to one message? mschaefer@8901: } mschaefer@8898: } mschaefer@8898: else mschaefer@8898: tkh = Double.NaN; mschaefer@8898: } mschaefer@8898: mschaefer@8898: // Soil kind mschaefer@8901: SoilKind kind = SoilKind.mobil; mschaefer@8901: if (doCalcTkh) { mschaefer@8898: try { mschaefer@8898: kind = soilKindFinder.findSoilKind(km); mschaefer@8898: } catch (Exception e) { mschaefer@8898: final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingSoilKind", null, label); mschaefer@8898: problems.addProblem(km, message); mschaefer@8901: //FIXME: cumulate problems to one message? mschaefer@8898: } mschaefer@8898: } mschaefer@8898: mschaefer@8898: /* // REMARK: bissl spielerei zum testen damit die sohlart nicht zu schnell wechselt mschaefer@8898: SoilKind kind; gernotbelger@8895: final boolean changeKind = Math.random() > 0.95; gernotbelger@8886: if (changeKind) { gernotbelger@8886: switch (lastKind) { gernotbelger@8886: case starr: gernotbelger@8886: kind = SoilKind.mobil; gernotbelger@8886: break; gernotbelger@8886: case mobil: gernotbelger@8886: default: gernotbelger@8886: kind = SoilKind.starr; gernotbelger@8886: break; gernotbelger@8886: } gernotbelger@8886: } else gernotbelger@8886: kind = lastKind; gernotbelger@8886: lastKind = kind; mschaefer@8898: */ mschaefer@8898: mschaefer@8898: //tkh = 100 + 10 * (Math.random() - 0.5); gernotbelger@8886: gernotbelger@8886: final double flowDepthTkh; gernotbelger@8886: final double tkhUp; gernotbelger@8886: final double tkhDown; gernotbelger@8886: switch (kind) { gernotbelger@8886: case starr: gernotbelger@8886: flowDepthTkh = wst - (meanBedHeight + tkh / 100); gernotbelger@8886: tkhUp = tkh; gernotbelger@8886: tkhDown = 0; gernotbelger@8886: break; gernotbelger@8886: gernotbelger@8886: case mobil: gernotbelger@8886: default: gernotbelger@8886: flowDepthTkh = wst - (meanBedHeight + tkh / 200); gernotbelger@8886: tkhUp = tkh / 2; gernotbelger@8886: tkhDown = -tkh / 2; gernotbelger@8886: break; gernotbelger@8886: } gernotbelger@8886: gernotbelger@8883: // REMARK: access the location once only during calculation gernotbelger@8883: final String location = LocationProvider.getLocation(river.getName(), km); gernotbelger@8883: gernotbelger@8883: // REMARK: access the gauge once only during calculation gernotbelger@8883: final Gauge gauge = findGauge(waterlevel, refGauge, gaugeIndex, km); gernotbelger@8883: gernotbelger@8883: final String gaugeLabel = gauge == null ? notinrange : gauge.getName(); gernotbelger@8883: mschaefer@8898: resultData.addRow(km, flowDepth, flowDepthTkh, kind, tkh, tkhUp, tkhDown, wst, discharge, wstLabel, mschaefer@8898: gaugeLabel, meanBedHeight, bedHeightLabel, location); gernotbelger@8883: } gernotbelger@8883: catch (final FunctionEvaluationException e) { gernotbelger@8883: /* should only happen if out of range */ gernotbelger@8883: e.printStackTrace(); gernotbelger@8883: /* simply ignore */ gernotbelger@8883: } gernotbelger@8883: gernotbelger@8877: } gernotbelger@8877: gernotbelger@8877: return resultData; gernotbelger@8877: } gernotbelger@8877: gernotbelger@8883: /** gernotbelger@8891: * Sohlbeschaffenheit (D50 Korndurchmesser aus Seddb) gernotbelger@8891: * Abhängig von Peiljahr gernotbelger@8891: */ mschaefer@8898: private BedQualityD50KmValueFinder loadBedMeasurements(final River river, final DoubleRange kmRange, final int soundingYear) { gernotbelger@8891: gernotbelger@8891: /* construct valid measurement time range */ gernotbelger@8891: final Calendar cal = Calendar.getInstance(); gernotbelger@8891: cal.clear(); gernotbelger@8891: gernotbelger@8891: cal.set(soundingYear - VALID_BED_MEASUREMENT_YEARS, 0, 1); gernotbelger@8891: final Date startTime = cal.getTime(); gernotbelger@8891: gernotbelger@8891: cal.set(soundingYear + VALID_BED_MEASUREMENT_YEARS, 11, 31); gernotbelger@8891: final Date endTime = cal.getTime(); gernotbelger@8891: mschaefer@8898: final BedQualityD50KmValueFinder finder = new BedQualityD50KmValueFinder(); mschaefer@8898: if (finder.loadValues(river, kmRange, new DateRange(startTime, endTime))) mschaefer@8898: return finder; mschaefer@8898: else mschaefer@8898: return null; gernotbelger@8891: } gernotbelger@8891: gernotbelger@8891: /** gernotbelger@8883: * Checks the year difference between waterlevels and sounding, and issues a warning if too big. gernotbelger@8883: * gernotbelger@8883: * Zeitraum Zeitliche Differenz [a] gernotbelger@8883: * X ≥ 1998 ± 3 gernotbelger@8883: * 1958 ≤ X < 1998 ± 6 gernotbelger@8883: * 1918 ≤ X < 1958 ± 12 gernotbelger@8883: * X < 1918 ± 25 gernotbelger@8883: */ gernotbelger@8884: private void checkYearDifference(final String label, final WaterlevelData waterlevel, final BedHeight sounding, final Calculation problems) { gernotbelger@8883: gernotbelger@8883: final Integer soundingYear = sounding.getYear(); gernotbelger@8883: if (soundingYear == null) gernotbelger@8883: return; gernotbelger@8883: gernotbelger@8883: final int wstYear = waterlevel.getYear(); gernotbelger@8883: if (wstYear < 0) gernotbelger@8883: return; gernotbelger@8883: gernotbelger@8883: final int maxDifference = getMaxDifferenceYears(soundingYear); gernotbelger@8883: gernotbelger@8883: final int difference = Math.abs(soundingYear - wstYear); gernotbelger@8883: if (difference > maxDifference) { gernotbelger@8894: final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.year_difference", null, label, wstYear, gernotbelger@8894: soundingYear); gernotbelger@8883: problems.addProblem(message); gernotbelger@8883: } gernotbelger@8883: } gernotbelger@8883: gernotbelger@8883: private int getMaxDifferenceYears(final int year) { gernotbelger@8883: gernotbelger@8883: if (year < 1918) gernotbelger@8883: return 25; gernotbelger@8883: gernotbelger@8883: if (1918 <= year && year < 1958) gernotbelger@8883: return 12; gernotbelger@8883: gernotbelger@8883: if (1958 <= year && year < 1998) gernotbelger@8883: return 6; gernotbelger@8883: gernotbelger@8883: /* >= 1998 */ gernotbelger@8883: return 3; gernotbelger@8883: } gernotbelger@8883: gernotbelger@8884: private Gauge findGauge(final WaterlevelData waterlevel, final Gauge refGauge, final GaugeIndex gaugeIndex, final double km) { gernotbelger@8882: gernotbelger@8882: // REMARK: using same logic as in WaterlevelExporter here gernotbelger@8882: gernotbelger@8882: final boolean showAllGauges = waterlevel.isShowAllGauges(); gernotbelger@8882: gernotbelger@8882: if (showAllGauges) gernotbelger@8882: return gaugeIndex.findGauge(km); gernotbelger@8882: gernotbelger@8882: if (refGauge.getRange().contains(km)) gernotbelger@8882: return refGauge; gernotbelger@8882: gernotbelger@8882: return null; gernotbelger@8882: } gernotbelger@8882: gernotbelger@8882: /* Checks if the discretisation of the waterlevel exceeds 1000m */ gernotbelger@8894: gernotbelger@8894: private void checkWaterlevelDiscretisation(final WKms wstKms, final DoubleRange calcRange, final Calculation problems) { gernotbelger@8894: gernotbelger@8882: final int size = wstKms.size(); gernotbelger@8882: for (int i = 0; i < size - 2; i++) { gernotbelger@8882: final double kmPrev = wstKms.getKm(i); gernotbelger@8882: final double kmNext = wstKms.getKm(i + 1); gernotbelger@8882: gernotbelger@8894: /* only check if we are within the calculation range */ gernotbelger@8894: if (calcRange.overlapsRange(new DoubleRange(kmPrev, kmNext))) { gernotbelger@8894: if (Math.abs(kmPrev - kmNext) > 1) { gernotbelger@8894: final String label = wstKms.getName(); gernotbelger@8882: gernotbelger@8894: final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.waterlevel_discretisation", null, label); gernotbelger@8894: problems.addProblem(kmPrev, message); gernotbelger@8894: } gernotbelger@8882: } gernotbelger@8882: } gernotbelger@8882: } gernotbelger@8882: gernotbelger@8894: private BedHeight loadBedHeight(final String soundingId) { gernotbelger@8877: gernotbelger@8883: // REMARK: absolutely unbelievable.... mschaefer@8898: // The way how bed-heights (and other data too) is accessed is different for nearly every calculation-type gernotbelger@8882: // throughout flys. gernotbelger@8877: // The knowledge on how to parse the datacage-ids is spread through the complete code-base... gernotbelger@8877: gernotbelger@8882: // We use here the way on how bed-heights are accessed by the BedDifferenceAccess/BedDifferenceCalculation, but gernotbelger@8882: // this is plain random gernotbelger@8877: final String[] parts = soundingId.split(";"); gernotbelger@8877: gernotbelger@8877: final BedHeightsArtifact artifact = (BedHeightsArtifact) RiverUtils.getArtifact(parts[0], this.context); gernotbelger@8877: gernotbelger@8877: final Integer bedheightId = artifact.getDataAsInteger("height_id"); gernotbelger@8883: // REMARK: this only works with type 'single'; unclear on how to distinguish from epoch data (or whatever the gernotbelger@8882: // other type means) gernotbelger@8877: // Luckily, the requirement is to only access 'single' data here. gernotbelger@8877: // final String bedheightType = artifact.getDataAsString("type"); gernotbelger@8877: gernotbelger@8883: // REMARK: BedDifferences uses this, but we also need the metadata of the BedHeight gernotbelger@8883: // REMARK: second absolutely awful thing: BedHeight is a hibernate binding class, accessing the database via gernotbelger@8882: // hibernate stuff gernotbelger@8877: // BedHeightFactory uses its own (direct) way of accessing the data, with its own implemented data classes. gernotbelger@8882: // return BedHeightFactory.getHeight(bedheightType, bedheightId, from, to); gernotbelger@8877: gernotbelger@8877: return BedHeight.getBedHeightById(bedheightId); gernotbelger@8877: } mschaefer@8898: mschaefer@8898: /** mschaefer@8898: * Calculates a transport body height mschaefer@8898: * @param h flow depth in m mschaefer@8898: * @param vm flow velocity in m mschaefer@8898: * @param d50 grain diameter D50 in m (!) mschaefer@8898: * @param tau shear stress in N/m^2 mschaefer@8898: * @return transport body height in cm (!) mschaefer@8898: */ mschaefer@8898: private double calculateTkh(double h, double vm, double d50, double tau) { mschaefer@8898: final double PHYS_G = 9.81; mschaefer@8898: final double PHYS_SPECGRAV_S = 2.6; mschaefer@8898: final double PHYS_VELOCCOEFF_N = 6; mschaefer@8898: final double PHYS_FORMCOEFF_ALPHA = 0.7; mschaefer@8898: final double PHYS_VISCOSITY_NUE = 1.3e-6; mschaefer@8898: final double PHYS_GRAIN_DENSITY_RHOS = 2603; mschaefer@8898: final double PHYS_WATER_DENSITY_RHO = 999.97; mschaefer@8898: mschaefer@8898: final double froude = vm / Math.sqrt(PHYS_G * h); mschaefer@8898: final double partReynolds = Math.sqrt((PHYS_SPECGRAV_S - 1) * PHYS_G * d50) / PHYS_VISCOSITY_NUE * d50; mschaefer@8898: final double critShields = 0.22 * Math.pow(partReynolds, -0.6) + 0.06 * Math.pow(10, 7.7 * Math.pow(partReynolds, -0.6)); mschaefer@8898: final double critTau = critShields * (PHYS_GRAIN_DENSITY_RHOS - PHYS_WATER_DENSITY_RHO) * PHYS_G * d50; mschaefer@8898: return 100 * h * (1 - Math.pow(froude, 2)) / (2 * PHYS_VELOCCOEFF_N * PHYS_FORMCOEFF_ALPHA) * (1 - critTau / tau); mschaefer@8898: } gernotbelger@8884: }