view artifacts/src/main/java/org/dive4elements/river/artifacts/sinfo/flood_duration/FloodDurationCalculator.java @ 9269:83ebeb620b5a

Station specific main value annotations in S-Info flood duration curve, corrected infrastructure flood duration calculation
author mschaefer
date Thu, 19 Jul 2018 08:07:03 +0200
parents 465347d12990
children bcbae86ce7b3
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.sinfo.flood_duration;

import java.util.ArrayList;
import java.util.Collection;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
import java.util.Set;

import org.apache.commons.lang.math.DoubleRange;
import org.dive4elements.artifacts.CallContext;
import org.dive4elements.river.artifacts.WINFOArtifact;
import org.dive4elements.river.artifacts.access.ComputationRangeAccess;
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.WQDay;
import org.dive4elements.river.artifacts.model.WQKms;
import org.dive4elements.river.artifacts.model.WstValueTable;
import org.dive4elements.river.artifacts.model.WstValueTable.QPosition;
import org.dive4elements.river.artifacts.model.WstValueTableFactory;
import org.dive4elements.river.artifacts.resources.Resources;
import org.dive4elements.river.artifacts.sinfo.common.GaugeDurationValuesFinder;
import org.dive4elements.river.artifacts.sinfo.common.RiverInfoProvider;
import org.dive4elements.river.artifacts.sinfo.common.SInfoResultType;
import org.dive4elements.river.artifacts.sinfo.flood_duration.RiversideRadioChoice.RiversideChoiceKey;
import org.dive4elements.river.exports.WaterlevelDescriptionBuilder;
import org.dive4elements.river.jfree.StickyAxisAnnotation;
import org.dive4elements.river.jfree.StickyAxisAnnotation.SimpleAxis;
import org.dive4elements.river.model.Attribute.AttributeKey;
import org.dive4elements.river.model.Gauge;
import org.dive4elements.river.model.sinfo.InfrastructureValue;

import gnu.trove.TDoubleArrayList;

/**
 * Calculation of the result rows of the flood duration of the infrastructures in a river km range
 * and selected main value durations
 *
 * @author Matthias Schäfer
 */
final class FloodDurationCalculator {

    private final List<ResultRow> rows = new ArrayList<>();

    private final RiverInfoProvider riverInfoProvider;

    private final CallContext context;

    public FloodDurationCalculator(final CallContext context, final RiverInfoProvider riverInfoProvider) {
        this.context = context;
        this.riverInfoProvider = riverInfoProvider;
    }

    /**
     * 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) {

        // Find all gauges of the calc range, and create the duration finders
        final Map<Gauge, GaugeDurationValuesFinder> durFinders = new HashMap<>();
        Gauge firstGauge = null;
        for (final Gauge gauge : this.riverInfoProvider.getRiver().determineGauges(calcRange.getMinimumDouble(), calcRange.getMaximumDouble())) {
            durFinders.put(gauge, GaugeDurationValuesFinder.loadValues(gauge, problems));
            if (firstGauge == null)
                firstGauge = gauge;
        }

        // Find all infrastructures within the calc range
        final AttributeKey bankKey = riverside.getAttributeKey();
        final List<InfrastructureValue> infras = InfrastructureValue.getValues(this.riverInfoProvider.getRiver(), calcRange.getMinimumDouble(),
                calcRange.getMaximumDouble(), bankKey);

        // Merge all stations (range/step, borders of gauge ranges, infrastructures)
        final Map<Double, InfrastructureValue> allStations = new HashMap<>();
        final Map<Double, InfrastructureValue> secondBank = new HashMap<>(); // any second infrastructure in case of both-banks-option
        addRangeStations(allStations, winfo);
        addGaugeLimits(allStations, durFinders.keySet(), calcRange.getMinimumDouble(), calcRange.getMaximumDouble());
        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);

        // Determine discharge state labels of the main values
        updateMainValueLabels(wqkmsArray, winfo, problems);

        // Load base wst table (river).wst - first run takes long time, then it's cached
        final WstValueTable wst = WstValueTableFactory.getTable(this.riverInfoProvider.getRiver());

        // 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]), wst, durFinders);
            this.rows.add(row);
            if (secondBank.containsKey(stationsSorted[i])) {
                final ResultRow row2 = ResultRow.create(row);
                calculateInfrastructure(row2, gauge, secondBank.get(stationsSorted[i]), wst, durFinders);
                this.rows.add(row2);
            }
        }

        //// 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);
        // }
        // }

        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);
    }

    /**
     * Calculates the duration curve for a station
     */
    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);
        }
        return new WQDay(days, ws, qs);
    }

    /**
     * Calculates the data for the Q main value lines in the duration curve chart
     */
    public List<StickyAxisAnnotation> calcMainValueQAnnotations(final Calculation problems, final double station, final FloodDurationCalculationResult result) {

        // Search the station in the previously calculated result rows and terminate if no infrastructure row found
        double station1 = station;
        if (Double.isNaN(station)) {
            for (final ResultRow row : result.getRows()) {
                station1 = row.getDoubleValue(GeneralResultType.station);
                break;
            }
        }
        final List<ResultRow> stationRows = searchStation(station1, result.getRows());
        if (stationRows.isEmpty() || (stationRows.get(0).getValue(SInfoResultType.infrastructuretype) == null)) {
            return new ArrayList<>();
        }
        final ResultRow row = stationRows.get(0);
        final List<StickyAxisAnnotation> annotations = new ArrayList<>();
        final List<DurationWaterlevel> wqds = (List<DurationWaterlevel>) row.getValue(SInfoResultType.customMultiRowColWaterlevel);
        for (final DurationWaterlevel wqd : wqds) {
            final String label = wqd.getBezeichnung().startsWith("W=") ? "Q(" + wqd.getBezeichnung() + ")" : wqd.getBezeichnung();
            final StickyAxisAnnotation annotation = new StickyAxisAnnotation(label, (float) wqd.getDischarge(), SimpleAxis.Y_AXIS,
                    FloodDurationCurveGenerator.YAXIS.Q.idx);
            annotation.setHitPoint(wqd.getFloodDurDaysPerYear());
            annotations.add(annotation);
        }
        return annotations;
    }

    /**
     * Calculates the data for the W main value lines in the duration curve chart
     */
    public List<StickyAxisAnnotation> calcMainValueWAnnotations(final Calculation problems, final double station, final FloodDurationCalculationResult result) {

        // Search the station in the previously calculated result rows and terminate if no infrastructure row found
        double station1 = station;
        if (Double.isNaN(station)) {
            for (final ResultRow row : result.getRows()) {
                station1 = row.getDoubleValue(GeneralResultType.station);
                break;
            }
        }
        final List<ResultRow> stationRows = searchStation(station1, result.getRows());
        if (stationRows.isEmpty() || (stationRows.get(0).getValue(SInfoResultType.infrastructuretype) == null)) {
            return new ArrayList<>();
        }
        final List<StickyAxisAnnotation> annotations = new ArrayList<>();
        final ResultRow row = stationRows.get(0);
        final List<DurationWaterlevel> wqds = (List<DurationWaterlevel>) row.getValue(SInfoResultType.customMultiRowColWaterlevel);
        for (final DurationWaterlevel wqd : wqds) {
            final String label = !wqd.getBezeichnung().startsWith("W=") ? "W(" + wqd.getBezeichnung() + ")" : wqd.getBezeichnung();
            final StickyAxisAnnotation annotation = new StickyAxisAnnotation(label, (float) wqd.getWaterlevel(), SimpleAxis.Y_AXIS,
                    FloodDurationCurveGenerator.YAXIS.W.idx);
            annotation.setHitPoint(wqd.getFloodDurDaysPerYear());
            annotations.add(annotation);
        }
        return annotations;
    }

    /**
     * Calculate the data for the W and Q lines in the duration curve chart for the infrastructure height and add to result
     * collection
     */
    public List<StickyAxisAnnotation> calcInfrastructureAnnotations(final Calculation problems, final double station,
            final FloodDurationCalculationResult result) {

        // Search the station in the previously calculated result rows and terminate if no infrastructure row found
        double station1 = station;
        if (Double.isNaN(station)) {
            for (final ResultRow row : result.getRows()) {
                station1 = row.getDoubleValue(GeneralResultType.station);
                break;
            }
        }
        final List<ResultRow> stationRows = searchStation(station1, result.getRows());
        if (stationRows.isEmpty() || (stationRows.get(0).getValue(SInfoResultType.infrastructuretype) == null)) {
            return new ArrayList<>();
        }
        // Same way as in MainValueWFacet and ..QFacet
        final List<StickyAxisAnnotation> annotations = new ArrayList<>();
        for (final ResultRow row : stationRows) {
            annotations.add(calcInfrastructureWAnnotation(row));
            annotations.add(calcInfrastructureQAnnotation(row));
        }
        return annotations;
    }

    /**
     * Searches the one or two rows of a station in a result rows collection
     */
    private List<ResultRow> searchStation(final double station, final Collection<ResultRow> rows) {
        final List<ResultRow> found = new ArrayList<>();
        for (final ResultRow row : rows) {
            if (row.getDoubleValue(GeneralResultType.station) > station + 0.0001)
                break;
            else if (row.getDoubleValue(GeneralResultType.station) > station - 0.0001)
                found.add(row);
        }
        return found;
    }

    /**
     * Calculates the Q annotation lines of an infrastructure
     */
    private StickyAxisAnnotation calcInfrastructureQAnnotation(final ResultRow row) {
        final String label = Resources.getMsg(this.context.getMeta(), "sinfo.chart.flood_duration.curve.infrastructure",
                "sinfo.chart.flood_duration.curve.infrastructure",
                SInfoResultType.infrastructuretype.exportValue(this.context, row.getValue(SInfoResultType.infrastructuretype))
                + ", " + SInfoResultType.riverside.exportValue(this.context, row.getValue(SInfoResultType.riverside)));
        final StickyAxisAnnotation annotation = new StickyAxisAnnotation(label, (float) row.getDoubleValue(SInfoResultType.floodDischarge),
                SimpleAxis.Y_AXIS, FloodDurationCurveGenerator.YAXIS.Q.idx);
        annotation.setHitPoint((float) row.getDoubleValue(SInfoResultType.floodDuration));
        return annotation;
    }

    /**
     * Calculates the W annotation lines of an infrastructure
     */
    private StickyAxisAnnotation calcInfrastructureWAnnotation(final ResultRow row) {
        final String label = Resources.getMsg(this.context.getMeta(), "sinfo.chart.flood_duration.curve.infrastructure",
                "sinfo.chart.flood_duration.curve.infrastructure",
                SInfoResultType.infrastructuretype.exportValue(this.context, row.getValue(SInfoResultType.infrastructuretype))
                + ", " + SInfoResultType.riverside.exportValue(this.context, row.getValue(SInfoResultType.riverside)));
        final StickyAxisAnnotation annotation = new StickyAxisAnnotation(label, (float) row.getDoubleValue(SInfoResultType.infrastructureHeight),
                SimpleAxis.Y_AXIS, FloodDurationCurveGenerator.YAXIS.W.idx);
        annotation.setHitPoint((float) row.getDoubleValue(SInfoResultType.floodDuration));
        return annotation;
    }

    /**
     * Adds to a stations map all stations corresponding to the active range and step
     */
    private void addRangeStations(final Map<Double, InfrastructureValue> allStations, final WINFOArtifact winfo) {
        for (final double station : new ComputationRangeAccess(winfo).getKms())
            allStations.put(Double.valueOf(station), null);
    }

    /**
     * Adds to a stations map all range limits of the gauges within the calc range
     */
    private void addGaugeLimits(final Map<Double, InfrastructureValue> allStations, final Set<Gauge> gauges, final double fromKm, final double toKm) {
        for (final Gauge gauge : gauges) {
            final Double kmA = Double.valueOf(gauge.getRange().getA().doubleValue());
            final Double kmB = Double.valueOf(gauge.getRange().getB().doubleValue());
            if (kmA > fromKm - 0.0001)
                allStations.put(kmA, null);
            if (kmB < toKm + 0.0001)
                allStations.put(kmB, null);
        }
    }

    /**
     * Adds to a stations map all (first) infrastructures of a station, and the second, if any, to another map
     */
    private void addInfrastructures(final Map<Double, InfrastructureValue> allStations, final Map<Double, InfrastructureValue> secondBank,
            final List<InfrastructureValue> infrastructures) {
        for (final InfrastructureValue infrastructure : infrastructures) {
            final Double station = infrastructure.getStation();
            if (!allStations.containsKey(station) || !(allStations.get(station) instanceof InfrastructureValue))
                allStations.put(station, infrastructure);
            else
                secondBank.put(station, infrastructure);
        }
    }

    /**
     * Returns a double array with a sorted stations set
     */
    private double[] sortStations(final Set<Double> stations) {
        final TDoubleArrayList sorted = new TDoubleArrayList();
        for (final Double station : stations)
            sorted.add(station.doubleValue());
        sorted.sort();
        return sorted.toNativeArray();
    }

    /**
     * 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)
        // (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;
        // das ergibt das W einer station für einen Abflusszustand;
        // bei Vorgabe eines Pegel-W wird vorher anhand der W-Q-Tabelle des Pegels ({gauge}.at in DischargeTable) das Q
        // 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);

        /* copy all problems */
        final Calculation winfoProblems = waterlevelData.getReport();
        final List<Problem> problems2 = winfoProblems.getProblems();
        if (problems2 != null) {
            for (final Problem problem : problems2) {
                problems.addProblem(problem);
            }
        }
        return (WQKms[]) waterlevelData.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) {

        for (int i = 0; i <= wqkmsArray.length - 1; i++)
            wqkmsArray[i].setName(buildWQDescription(wqkmsArray[i], winfo));
    }

    /**
     *
     * @param wqkms
     * @param descBuilder
     * @return
     */
    private String buildWQDescription(final WQKms wqkms, final WINFOArtifact winfo) {
        final WaterlevelDescriptionBuilder descBuilder = new WaterlevelDescriptionBuilder(winfo, this.context);
        final String description = descBuilder.getDesc(wqkms);
        if (!description.isEmpty() && Character.isDigit(description.charAt(0))) {
            if (winfo.isQ())
                return "Q=" + description;
            else
                return "W=" + description;
        }
        else
            return description;
    }

    /**
     * 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 ResultRow row = ResultRow.create();
        row.putValue(GeneralResultType.station, station);
        row.putValue(SInfoResultType.infrastructuretype, null); // is replaced later for an infrastructure
        row.putValue(SInfoResultType.floodDuration, Double.NaN); // is replaced later for an infrastructure

        final String gaugeLabel = this.riverInfoProvider.findGauge(station, (gauge == firstGauge));
        row.putValue(SInfoResultType.gaugeLabel, gaugeLabel);

        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 int overflowDays = (int) Math.round(underflowDaysToOverflowDays(durationFinder.getDuration(wqKms.getQ(kmIndex))));

            final DurationWaterlevel dw = new DurationWaterlevel(wqKms.getW(kmIndex), overflowDays, wqKms.getQ(kmIndex), wqKms.getName());
            waterlevels.add(dw);
        }
        row.putValue(SInfoResultType.customMultiRowColWaterlevel, waterlevels);

        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
     */
    private void calculateInfrastructure(final ResultRow row, final Gauge gauge, final InfrastructureValue infrastructure, final WstValueTable wst,
            final Map<Gauge, GaugeDurationValuesFinder> durFinders) {

        // Interpolate the infrastructure height in the wst table to get the corresponding Q
        final Calculation problems = new Calculation();
        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
        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
        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());
    }

    /**
     * Translates underflow duration into overflow duration
     */
    private double underflowDaysToOverflowDays(final double underflowDays) {
        return 365 - underflowDays;
    }
}

http://dive4elements.wald.intevation.org