Mercurial > dive4elements > river
view artifacts/src/main/java/org/dive4elements/river/artifacts/bundu/bezugswst/BezugswstCalculation.java @ 9588:c57caff9b00b
Punkt 10.6 CSV-Ausgabe Abflusszeitreihenlänge
author | gernotbelger |
---|---|
date | Thu, 10 Jan 2019 11:56:39 +0100 |
parents | b9c87bbff6a4 |
children | 17414e70746e |
line wrap: on
line source
/** Copyright (C) 2017 by Bundesanstalt für Gewässerkunde * Software engineering by * Björnsen Beratende Ingenieure GmbH * Dr. Schumacher Ingenieurbüro für Wasser und Umwelt * * This file is Free Software under the GNU AGPL (>=v3) * and comes with ABSOLUTELY NO WARRANTY! Check out the * documentation coming with Dive4Elements River for details. */ package org.dive4elements.river.artifacts.bundu.bezugswst; import java.util.ArrayList; import java.util.Calendar; import java.util.List; import org.dive4elements.artifacts.CallContext; import org.dive4elements.river.artifacts.access.FixRealizingAccess; import org.dive4elements.river.artifacts.bundu.BUNDUArtifact; import org.dive4elements.river.artifacts.bundu.BunduResultType; import org.dive4elements.river.artifacts.common.AbstractResultType; import org.dive4elements.river.artifacts.common.GeneralResultType; import org.dive4elements.river.artifacts.common.ResultRow; import org.dive4elements.river.artifacts.model.Calculation; import org.dive4elements.river.artifacts.model.Calculation.Problem; import org.dive4elements.river.artifacts.model.CalculationResult; import org.dive4elements.river.artifacts.model.WQKms; import org.dive4elements.river.artifacts.model.fixings.FixRealizingCalculation; import org.dive4elements.river.artifacts.model.fixings.FixRealizingResult; import org.dive4elements.river.artifacts.model.river.RiverInfoProvider; import org.dive4elements.river.artifacts.resources.Resources; import org.dive4elements.river.artifacts.services.DynamicMainValuesTimeRangeDeterminationService; import org.dive4elements.river.artifacts.services.DynamicMainValuesTimeRangeDeterminationService.GaugeInfoResult; import org.dive4elements.river.artifacts.sinfo.tkhstate.BedHeightsFinder; import org.dive4elements.river.artifacts.sinfo.tkhstate.WinfoArtifactWrapper; import org.dive4elements.river.artifacts.sinfo.util.CalculationUtils; import org.dive4elements.river.artifacts.sinfo.util.RiverInfo; import org.dive4elements.river.artifacts.sinfo.util.WstInfo; import org.dive4elements.river.artifacts.states.WaterlevelData; import org.dive4elements.river.exports.WaterlevelDescriptionBuilder; import org.dive4elements.river.model.BedHeightValueType; import org.dive4elements.river.model.River; import org.dive4elements.river.utils.DoubleUtil; class BezugswstCalculation { // private static Logger log = Logger.getLogger(BezugswstCalculation.class); /** * Additional depth (m) to compute the excavation volume */ private final static double EXCAVATION_DEPTH = 0.2; // REMARK Sollte von außen einstellbar sein /** * Excavation costs (euro) per cubic meter */ private final static double EXPENSE_PER_CBM = 12.0; // REMARK Sollte von außen einstellbar sein private final CallContext context; private final List<ResultRow> rows; private Double missKmFrom; private Double missKmTo; public BezugswstCalculation(final CallContext context) { this.context = context; this.rows = new ArrayList<>(); } /** * Calculates the result rows of a bundu bzws workflow */ public CalculationResult calculate(final BUNDUArtifact bunduartifact) { // Get input data final String user = CalculationUtils.findArtifactUser(this.context, bunduartifact); final BunduAccess access = new BunduAccess(bunduartifact); final River river = access.getRiver(); final RiverInfo riverInfo = new RiverInfo(river); final String calcModeLabel = Resources.getMsg(this.context.getMeta(), "bundu_bezugswst"); final boolean preprocessing = access.getPreprocessing(); final int startYear = access.getStartYear(); final int endYear = access.getBezugsJahr(); final Integer ud = access.getUd(); this.missKmFrom = access.getMissingVolFrom(); this.missKmTo = access.getMissingVolTo(); final GaugeInfoResult gi = DynamicMainValuesTimeRangeDeterminationService .getCommonTimeRangeForGauges(river.determineGauges(access.getLowerKm(), access.getUpperKm()), startYear, endYear, this.context.getMeta()); final int globalAdjustedEndYear = gi.getGlobalEndYear(); final int globalAdjustedStartYear = gi.getGlobalStartYear(); final BezugswstCalculationResults results = new BezugswstCalculationResults(calcModeLabel, user, riverInfo, access.getRange(), access.isCalculateMissingVolume()); final Calculation problems = new Calculation(); // Calculate the wspl for the selected river range as in fixa awspl bunduartifact.addStringData("wq_isq", "true"); final WinfoArtifactWrapper winfo = new WinfoArtifactWrapper(bunduartifact); final RiverInfoProvider riverInfoProvider = RiverInfoProvider.forRange(this.context, river, access.getRange()); final FixRealizingResult fixResult = calculateWspl(bunduartifact, problems); if (fixResult == null) return new CalculationResult(results, problems); final WQKms wqkms = fixResult.getWQKms()[0]; // We have no wst year as the wst is created by a calculation; we do not need it though final int wspYear = -1; // Remark: showAllGauges true for Fixierungsanalyse, false for WInfo, so true here as well final boolean showAllGauges = true; final WaterlevelData waterlevel = new WaterlevelData(wqkms, wspYear, showAllGauges, true); final RiverInfoProvider riverInfoProvider2 = riverInfoProvider.forWaterlevel(waterlevel); final WstInfo wstInfo = new WstInfo(wqkms.getName(), 0, riverInfoProvider2.getReferenceGauge(), true); // Fetch the bed levels of the selected sounding final Integer bedHeightId = access.getBedHeightID(); final BedHeightsFinder bedHeightsFinder = (bedHeightId != null) ? BedHeightsFinder.forId(problems, bedHeightId, access.getRange()) : BedHeightsFinder.NullFinder(); // Fetch the river channel data final ChannelFinder channelFinder = ChannelFinder.loadValues(problems, river, access.getBezugsJahr()); if (channelFinder == null) return new CalculationResult(results, problems); // Compute the result rows for (int i = 0; i <= wqkms.size() - 1; i++) { this.rows.add(createRow(wqkms.getKm(i), wqkms.getW(i), wqkms.getQ(i), bedHeightsFinder, channelFinder, riverInfoProvider2, wstInfo)); } // Compute the missing volumes if (access.isCalculateMissingVolume()) { if ((bedHeightsFinder == null) || bedHeightsFinder.isNull()) return new CalculationResult(results, problems); computeMissingVolumes(problems); final BedQualityCalculator bqCalculator = computeDensities(problems, bunduartifact, access, river); computeMissingMasses(problems, bqCalculator); } // Add the result to the results collection final WaterlevelDescriptionBuilder descBuilder = new WaterlevelDescriptionBuilder(winfo, this.context); final String qtext = descBuilder.getMetadataQ(); final BezugswstMainCalculationResult result = new BezugswstMainCalculationResult("bundu-bzws", this.rows, bedHeightsFinder.getInfo(), wstInfo, access.getFunction(), preprocessing, globalAdjustedStartYear, globalAdjustedEndYear, ud, qtext, wqkms, this.missKmFrom, this.missKmTo); results.addResult(result, problems); // Create the missing volume results if (access.getMissingVolFrom() != null) { final String title1 = Resources.getMsg(this.context.getMeta(), "bundu.export.csv.title.bezugswst.result1"); final BezugswstMissVolCalculationResult1 r1 = new BezugswstMissVolCalculationResult1(title1, this.rows); results.addResult(r1, null); final String title2 = Resources.getMsg(this.context.getMeta(), "bundu.export.csv.title.bezugswst.result2"); final List<ResultRow> rows2 = copyMissRows(); final BezugswstMissVolCalculationResult2 r2 = new BezugswstMissVolCalculationResult2(title2, rows2); results.addResult(r2, null); final String title3 = Resources.getMsg(this.context.getMeta(), "bundu.export.csv.title.bezugswst.result3"); final List<ResultRow> totalRows = new ArrayList<>(); totalRows.add(createTotalsRow(problems)); final BezugswstMissVolCalculationResult3 r3 = new BezugswstMissVolCalculationResult3(title3, totalRows); results.addResult(r3, null); } return new CalculationResult(results, problems); } /** * Calculates a w-q-longitudinal section for a river range and Q specified in an artifact */ private FixRealizingResult calculateWspl(final BUNDUArtifact bundu, final Calculation problems) { final FixRealizingAccess access = new FixRealizingAccess(bundu); final FixRealizingCalculation calc = new FixRealizingCalculation(access); final CalculationResult res = calc.calculate(this.context.getMeta()); final FixRealizingResult fixRes = (FixRealizingResult) res.getData(); final List<Problem> problems2 = res.getReport().getProblems(); if (problems2 != null) { for (final Problem problem : problems2) { problems.addProblem(problem); } } return fixRes; } /** * Create a result row for a station */ private ResultRow createRow(final double station, final double w, final double q, final BedHeightsFinder bedHeightsFinder, final ChannelFinder channelFinder, final RiverInfoProvider riverInfoProv, final WstInfo wstInfo) { // Set W and Q final ResultRow row = ResultRow.create(); row.putValue(GeneralResultType.station, station); row.putValue(BunduResultType.bezugswst, w); row.putValue(GeneralResultType.dischargeQwithUnit, q); row.putValue(GeneralResultType.waterlevelLabel, wstInfo.getLabel()); row.putValue(GeneralResultType.gaugeLabel, riverInfoProv.findGauge(station)); row.putValue(GeneralResultType.location, riverInfoProv.getLocation(station)); if (bedHeightsFinder.getInfo() != null) row.putValue(BunduResultType.sounding, bedHeightsFinder.getInfo().getDescription()); else row.putValue(BunduResultType.sounding, ""); // Set bed and channel bottom height final double msh = bedHeightsFinder.getMeanBedHeight(station); row.putValue(BunduResultType.heightMeanBed, msh); if (!Double.isNaN(w) && !Double.isNaN(msh)) row.putValue(BunduResultType.flowdepthMeanBed, w - msh); else row.putValue(BunduResultType.flowdepthMeanBed, Double.NaN); final double channelDepth = channelFinder.getDepth(station); row.putValue(BunduResultType.channelDepth, channelDepth); double channelHeight; if (!Double.isNaN(w) && !Double.isNaN(channelDepth)) channelHeight = w - channelDepth; else channelHeight = Double.NaN; row.putValue(BunduResultType.channelLowerEdge, channelHeight); final double channelWidth = channelFinder.getWidth(station); row.putValue(BunduResultType.channelWidth, channelWidth); if (!Double.isNaN(channelHeight)) { if (msh > channelHeight + 0.001) row.putValue(BunduResultType.missDepthMeanBed, msh - channelHeight); else row.putValue(BunduResultType.missDepthMeanBed, 0.0); } // Set field heights and missing heights final List<Double> fieldHeights = new ArrayList<>(); final List<Double> fieldDepths = new ArrayList<>(); final List<Double> fieldMissDepths = new ArrayList<>(); final List<Double> fieldMissWidths = new ArrayList<>(); final List<Double> fieldNulls = new ArrayList<>(); int missFieldCnt = 0; for (int i = BedHeightValueType.FIELD_FIRST_INDEX; i <= BedHeightValueType.FIELD_LAST_INDEX; i++) { final double h = bedHeightsFinder.getFieldHeight(station, i); fieldHeights.add(Double.valueOf(h)); fieldDepths.add(Double.valueOf(w - h)); if (h > channelHeight + 0.001) { missFieldCnt++; fieldMissDepths.add(Double.valueOf(h - channelHeight)); fieldMissWidths.add(Double.valueOf(channelWidth / BedHeightValueType.FIELD_LAST_INDEX)); } else { fieldMissDepths.add(Double.valueOf(0.0)); fieldMissWidths.add(Double.valueOf(0.0)); } fieldNulls.add(Double.NaN); } if (isKmInMissingVolumeRange(station)) { row.putValue(BunduResultType.missDepthFields, fieldMissDepths); row.putValue(BunduResultType.missWidthFields, fieldMissWidths); row.putValue(BunduResultType.hasMissingDepth, (missFieldCnt >= 1)); } else { row.putValue(BunduResultType.missDepthFields, fieldNulls); row.putValue(BunduResultType.missWidthFields, fieldNulls); row.putValue(BunduResultType.hasMissingDepth, null); } row.putValue(BunduResultType.missVolumeFields, fieldNulls); row.putValue(BunduResultType.missMassFields, fieldNulls); row.putValue(BunduResultType.bedHeightFields, fieldHeights); row.putValue(BunduResultType.depthFields, fieldDepths); // Preset the missing volume fields with NaN row.putValue(BunduResultType.excavationCosts, Double.NaN); row.putValue(BunduResultType.excavationVolume, Double.NaN); row.putValue(BunduResultType.missVolumeMeanBed, Double.NaN); row.putValue(BunduResultType.missMassMeanBed, Double.NaN); row.putValue(BunduResultType.missVolumeTotal, Double.NaN); row.putValue(BunduResultType.missMassTotal, Double.NaN); row.putValue(BunduResultType.density, Double.NaN); row.putValue(BunduResultType.missStationRangeFrom, Double.NaN); row.putValue(BunduResultType.missStationRangeTo, Double.NaN); return row; } /** * Computes the missing volumes in a km range */ private void computeMissingVolumes(final Calculation problems) { // Search start km int first = -1; for (int j = 0; j <= this.rows.size() - 1; j++) { if (isKmInMissingVolumeRange(this.rows.get(j).getDoubleValue(GeneralResultType.station))) { first = j; break; } } if (first < 0) return; // Calculate all kms in missing volume calc range double km; int last = this.rows.size() - 1; for (int i = first; i <= this.rows.size() - 1; i++) { km = this.rows.get(i).getDoubleValue(GeneralResultType.station); if (!isKmInMissingVolumeRange(km)) break; if (km > this.missKmTo.doubleValue() - 0.0001) last = i; final List<Double> areas = new ArrayList<>(); final List<Double> volumes = new ArrayList<>(); double vTotal = 0.0; double vExcav = 0.0; for (int j = BedHeightValueType.FIELD_FIRST_INDEX; j <= BedHeightValueType.FIELD_LAST_INDEX; j++) { if (getFieldValue(i, BunduResultType.missDepthFields, j) > 0.0001) { computeMissingVolume(volumes, areas, i, first, last, j); vTotal += volumes.get(j - 1); vExcav += volumes.get(j - 1) + areas.get(j - 1) * EXCAVATION_DEPTH; } else { volumes.add(Double.valueOf(0.0)); areas.add(Double.valueOf(0.0)); } } final double[] meanBedVolumeArea = computeMeanBedMissingAreaAndVolume(i, first, last); this.rows.get(i).putValue(BunduResultType.missVolumeMeanBed, meanBedVolumeArea[0]); this.rows.get(i).putValue(BunduResultType.missAreaMeanBed, meanBedVolumeArea[1]); this.rows.get(i).putValue(BunduResultType.missVolumeFields, volumes); this.rows.get(i).putValue(BunduResultType.missAreaFields, areas); this.rows.get(i).putValue(BunduResultType.missVolumeTotal, vTotal); this.rows.get(i).putValue(BunduResultType.excavationVolume, vExcav); this.rows.get(i).putValue(BunduResultType.excavationCosts, vExcav * EXPENSE_PER_CBM); } } /** * Computes the missing volume of a field of a km row */ private void computeMissingVolume(final List<Double> volumes, final List<Double> areas, final int current, final int first, final int last, final int field) { final double areaCurr = missingArea(current, first, last, field); final double areaPrev = missingArea(current - 1, first, last, field); final double areaNext = missingArea(current + 1, first, last, field); final double kmCurr = missingKm(current); final double kmPrev = missingKm(current - 1); final double kmNext = missingKm(current + 1); final double area1 = Double.isNaN(kmPrev) ? 0.0 : 0.5 * (areaCurr + areaPrev); final double area2 = Double.isNaN(kmNext) ? 0.0 : 0.5 * (areaCurr + areaNext); final double volume = Double.valueOf((Math.abs(kmCurr - kmPrev) * 500 * area1) + (Math.abs(kmNext - kmCurr) * 500 * area2)); volumes.add(volume); if (!Double.isNaN(volume)) areas.add(Double.valueOf(area1 + area2)); else areas.add(Double.NaN); } /** * Gets the missing area of a field and a row if in range, otherwise 0.0 */ private double missingArea(final int rowIndex, final int first, final int last, final int fieldIndex) { if ((first <= rowIndex) && (rowIndex <= last)) return getFieldValue(rowIndex, BunduResultType.missDepthFields, fieldIndex) * getFieldValue(rowIndex, BunduResultType.missWidthFields, fieldIndex); else return 0.0; } /** * Computes the missing area and volume of the mean bed level of a km row */ private double[] computeMeanBedMissingAreaAndVolume(final int current, final int first, final int last) { final double areaCurr = meanBedMissingArea(current, first, last); if (areaCurr < 0.0001) return new double[] { 0.0, 0.0 }; final double areaPrev = meanBedMissingArea(current - 1, first, last); final double areaNext = meanBedMissingArea(current + 1, first, last); final double kmCurr = missingKm(current); final double kmPrev = missingKm(current - 1); final double kmNext = missingKm(current + 1); final double area1 = Double.isNaN(kmPrev) ? 0.0 : 0.5 * (areaCurr + areaPrev); final double area2 = Double.isNaN(kmNext) ? 0.0 : 0.5 * (areaCurr + areaNext); final double volume = Double.valueOf((Math.abs(kmCurr - kmPrev) * 500 * area1) + (Math.abs(kmNext - kmCurr) * 500 * area2)); final double area = Double.isNaN(volume) ? Double.NaN : Double.valueOf(area1 + area2); return new double[] { volume, area }; } /** * Gets the missing area of the mean bed level and a row if in range, otherwise 0.0 */ private double meanBedMissingArea(final int rowIndex, final int first, final int last) { if ((first <= rowIndex) && (rowIndex <= last)) { final double dh = this.rows.get(rowIndex).getDoubleValue(BunduResultType.channelDepth) - this.rows.get(rowIndex).getDoubleValue(BunduResultType.flowdepthMeanBed); if (dh > 0.0) return dh * this.rows.get(rowIndex).getDoubleValue(BunduResultType.channelWidth); return 0.0; } return 0.0; } /** * Gets the km of a row if within range, otherwise NaN */ private double missingKm(final int rowIndex) { if ((0 <= rowIndex) && (rowIndex <= this.rows.size() - 1)) return this.rows.get(rowIndex).getDoubleValue(GeneralResultType.station); return Double.NaN; } /** * Create a density calculator and compute the densities of the missing volume km range and time period according to the * reference year */ private BedQualityCalculator computeDensities(final Calculation problems, final BUNDUArtifact bunduartifact, final BunduAccess access, final River river) { final BedQualityCalculator bqCalculator = new BedQualityCalculator(this.context, bunduartifact); // REMARK 10km tolerance at start and end to enable interpolation there final double[] kms = DoubleUtil.explode(access.getMissingVolFrom().doubleValue() - 10.0, access.getMissingVolTo().doubleValue() + 10.0, access.getStep().doubleValue() / 1000); final Calendar endDay = Calendar.getInstance(); endDay.set(access.getBezugsJahr().intValue(), 11, 31); final Calendar startDay = Calendar.getInstance(); // TODO Spezialregelung für den Rhein (bis 1999, 2000 bis 2009, ab 2010) startDay.set(endDay.get(Calendar.YEAR) - 20, 0, 1); bqCalculator.execute(problems, river, kms, startDay.getTime(), endDay.getTime()); return bqCalculator; } /** * Computes the missing masses */ private void computeMissingMasses(final Calculation problems, final BedQualityCalculator densityFinder) { for (final ResultRow row : this.rows) { @SuppressWarnings("unchecked") final List<Double> volumes = (List<Double>) row.getValue(BunduResultType.missVolumeFields); if ((volumes == null) || Double.isNaN(volumes.get(0))) continue; final double density = getDensity(row.getDoubleValue(GeneralResultType.station), densityFinder); final List<Double> masses = new ArrayList<>(); double kmTotal = 0.0; for (int j = BedHeightValueType.FIELD_FIRST_INDEX; j <= BedHeightValueType.FIELD_LAST_INDEX; j++) { final double m = volumes.get(j - 1) * density; masses.add(m); if (!Double.isNaN(m)) kmTotal += m; } row.putValue(BunduResultType.density, density); row.putValue(BunduResultType.missMassFields, masses); row.putValue(BunduResultType.missMassTotal, kmTotal); row.putValue(BunduResultType.missMassMeanBed, row.getDoubleValue(BunduResultType.missVolumeMeanBed) * density); } } /** * Gets a value of one of the field list types of a row * * @param rowIndex * @param type * @param fieldIndex * 1-based field index */ private double getFieldValue(final int rowIndex, final AbstractResultType type, final int fieldIndex) { @SuppressWarnings("unchecked") final List<Double> values = (List<Double>) this.rows.get(rowIndex).getValue(type); return values.get(fieldIndex - 1); } /** * Computes the volume and mass total of all rows with missing volumes */ private ResultRow createTotalsRow(final Calculation problems) { // Search start km double vTotal = 0.0; double mTotal = 0.0; double eTotal = 0.0; double cTotal = 0.0; for (final ResultRow row : this.rows) { final double volume = row.getDoubleValue(BunduResultType.missVolumeTotal); final double mass = row.getDoubleValue(BunduResultType.missMassTotal); final double excavation = row.getDoubleValue(BunduResultType.excavationVolume); final double costs = row.getDoubleValue(BunduResultType.excavationCosts); if (!Double.isNaN(volume) && !Double.isNaN(mass)) { vTotal += volume; mTotal += mass; eTotal += excavation; cTotal += costs; } } final ResultRow sumRow = ResultRow.create(); sumRow.putValue(BunduResultType.missStationRangeFrom, Double.valueOf(this.missKmFrom)); sumRow.putValue(BunduResultType.missStationRangeTo, Double.valueOf(this.missKmTo)); sumRow.putValue(BunduResultType.missVolumeTotal, vTotal); sumRow.putValue(BunduResultType.missMassTotal, mTotal); sumRow.putValue(BunduResultType.excavationVolumeTotal, eTotal); sumRow.putValue(BunduResultType.excavationCostsTotal, cTotal); return sumRow; } /** * Copies the rows of the missing volume calculation range into a new list */ private List<ResultRow> copyMissRows() { final List<ResultRow> missRows = new ArrayList<>(); for (final ResultRow row : this.rows) { final double km = row.getDoubleValue(GeneralResultType.station); if (isKmInMissingVolumeRange(km)) missRows.add(row); } return missRows; } /** * Gets the density of a km from the densities calculation */ private double getDensity(final double km, final BedQualityCalculator densityFinder) { return densityFinder.getDensity(km); } /** * Checks whether a km lies in the missing volume calculation range */ private boolean isKmInMissingVolumeRange(final double km) { if ((this.missKmFrom == null) || (this.missKmTo == null)) return false; return (this.missKmFrom.doubleValue() - 0.0001 < km) && (km < this.missKmTo.doubleValue() + 0.0001); } }