view artifacts/src/main/java/org/dive4elements/river/artifacts/sinfo/flood_duration/FloodDurationCalculator.java @ 9631:6ecd1a28017f

Nachtrag Pos. 20: Q theme for height chart added, calculator corrected (rows for km without any infrastructure)
author mschaefer
date Tue, 15 Oct 2019 12:26:13 +0200
parents 07f02019065e
children
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.HashMap;
import java.util.HashSet;
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.model.river.RiverInfoProvider;
import org.dive4elements.river.artifacts.sinfo.common.GaugeDurationValuesFinder;
import org.dive4elements.river.artifacts.sinfo.common.SInfoResultType;
import org.dive4elements.river.artifacts.sinfo.flood_duration.InfrastructureServerClientXChange.Element;
import org.dive4elements.river.exports.WaterlevelDescriptionBuilder;
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 RiverInfoProvider riverInfoProvider2;
    private final Map<Double, List<InfrastructureValue>> stationInfras;

    private final CallContext context;

    public FloodDurationCalculator(final CallContext context, final RiverInfoProvider riverInfoProvider) {
        this.context = context;
        this.riverInfoProvider = riverInfoProvider;
        this.riverInfoProvider2 = null;
        this.stationInfras = new HashMap<>();
    }

    /**
     * Calculate the infrastructures flood duration result rows
     */
    public void execute(final Calculation problems, final String label, final DoubleRange calcRange, final AttributeKey riverside,
            final List<Element> infrastructureChoices, final boolean withWspl, 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<>();
        for (final Gauge gauge : this.riverInfoProvider.getGauges()) {
            durFinders.put(gauge, GaugeDurationValuesFinder.loadValues(gauge, problems));
        }

        // Find all infrastructures within the calc range
        final Set<String> choices = new HashSet<>();
        for (final Element ifch : infrastructureChoices)
            choices.add(ifch.getGroupLabel() + "\t" + ifch.getTypeLabel());
        final List<InfrastructureValue> infras = InfrastructureValue.getValues(this.riverInfoProvider.getRiver(), calcRange.getMinimumDouble(),
                calcRange.getMaximumDouble(), riverside, choices);

        // Merge all stations (range/step, borders of gauge ranges, infrastructures)
        // final Map<Double, InfrastructureValue> allStations = new HashMap<>();
        this.stationInfras.clear();
        // final Map<Double, InfrastructureValue> secondBank = new HashMap<>(); // any second infrastructure in case of
        // both-banks-option
        addRangeStations(winfo);
        addGaugeLimits(durFinders.keySet(), calcRange.getMinimumDouble(), calcRange.getMaximumDouble());
        addInfrastructures(infras);
        final double[] stationsSorted = sortStations(this.stationInfras.keySet());

        // Calculate W and Q for all stations and the selected discharge states/waterlevels
        final WQKms[] wqkmsArray = calculateWsts(winfo, withWspl, stationsSorted, problems);
        // final WaterlevelData waterlevel = new WaterlevelData(wqkmsArray[0], -1, false, true);
        // this.riverInfoProvider2 = this.riverInfoProvider.forWaterlevel(waterlevel);
        this.riverInfoProvider2 = this.riverInfoProvider.forReferenceRange(calcRange, false);
        // this.riverInfoProvider2.cleanupGaugesAtStart(calcRange);

        // Determine discharge state labels of the waterlevels
        updateWstLabels(wqkmsArray, winfo, problems);

        final Map<Gauge, List<Double>> gaugeWstDurations = new HashMap<>();
        if (withWspl)
            calcGaugeWstDurations(winfo, new ArrayList<>(durFinders.keySet()), gaugeWstDurations, durFinders);

        // Load base wst table (river).wst
        // (should be in cache since already used in calculateWaterlevels (winfo.computeWaterlevelData)
        final WstValueTable wst = WstValueTableFactory.getTable(this.riverInfoProvider2.getRiver());

        final Set<FloodDurationInfrastructureChoice> infrastructures = new HashSet<>();

        // Create the result rows, and calculate and add the flood durations etc.
        ResultRow row;
        boolean starting;
        for (int i = 0; i <= stationsSorted.length - 1; i++) {
            final Gauge gauge = this.riverInfoProvider2.getGauge(stationsSorted[i], true);
            row = createRow(stationsSorted[i], wqkmsArray, gaugeWstDurations.get(gauge), i);
            starting = true;
            if (this.stationInfras.containsKey(stationsSorted[i])) {
                for (final InfrastructureValue infra : this.stationInfras.get(stationsSorted[i])) {
                    if (!starting)
                        row = createRow(stationsSorted[i], wqkmsArray, gaugeWstDurations.get(gauge), i);
                    calculateInfrastructure(row, gauge, infra, wst, durFinders, infrastructures);
                    this.rows.add(row);
                    starting = false;
                }
            }
            if (starting)
                this.rows.add(row);
        }

        // Get the labels of the selected waterlevels
        final String[] wstLabels = new String[wqkmsArray.length];
        for (int i = 0; i <= wqkmsArray.length - 1; i++)
            wstLabels[i] = wqkmsArray[i].getName();

        results.addResult(new FloodDurationCalculationResult(label, wstLabels, this.rows, withWspl, infrastructures), problems);
    }

    /**
     * Calculates the duration curve for a station
     * (other than the version 3.2.1 W-Info Dauerlinie the wst column positions
     * are taken from the Q values of the gauge's Q-D-table)
     */
    public WQDay calcWQDays(final Calculation problems, final double station, final WINFOArtifact winfo) {

        final CalculationResult res = winfo.getDurationCurveData();
        final WQDay wqday = (WQDay) res.getData();
        if (wqday == null)
            return null;

        final int[] odays = new int[wqday.size()];
        for (int i = 0; i <= odays.length - 1; i++)
            odays[i] = 365 - wqday.getDay(i); // TODO Eigentlich 365.25, ist aber mit getDay als int sinnlos
        return new WQDay(odays, wqday.getWs(), wqday.getQs());
    }

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

    /**
     * Adds to a stations map all range limits of the gauges within the calc range
     */
    private void addGaugeLimits(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)
                this.stationInfras.put(kmA, new ArrayList<InfrastructureValue>());
            if (kmB < toKm + 0.0001)
                this.stationInfras.put(kmB, new ArrayList<InfrastructureValue>());
        }
    }

    /**
     * Adds all infrastructures of a station to the station map
     */
    private void addInfrastructures(final List<InfrastructureValue> infrastructures) {
        for (final InfrastructureValue infrastructure : infrastructures) {
            final Double station = infrastructure.getStation();
            if (this.stationInfras.containsKey(station))
                this.stationInfras.get(station).add(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[] calculateWsts(final WINFOArtifact winfo, final boolean withWspl, final double[] stations, final Calculation problems) {
        // First run may take long, further runs are faster since WstValueTable is in cache then
        // (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.

        if (!withWspl)
            return new WQKms[] {};

        final CalculationResult wstsData = winfo.computeWaterlevelData(stations);

        /* copy all problems */
        final Calculation winfoProblems = wstsData.getReport();
        final List<Problem> problems2 = winfoProblems.getProblems();
        if (problems2 != null) {
            for (final Problem problem : problems2) {
                problems.addProblem(problem);
            }
        }
        return (WQKms[]) wstsData.getData();
    }

    /**
     * Determines the waterlevel/discharge state labels for the selected Q or W values and sets them in the WQKms array
     */
    private void updateWstLabels(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));
    }

    /**
     * Builds the description label of a waterlevel
     */
    private String buildWQDescription(final WQKms wqkms, final WINFOArtifact winfo) {

        final WaterlevelDescriptionBuilder descBuilder = new WaterlevelDescriptionBuilder(winfo, this.context);
        // TODO Zwischen numerischem Q-Wert und Dauerzahl-Hauptwert (0..364) unterscheiden
        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;
    }

    /**
     * Calculates the flood durations of the Qs of the waterlevels/discharge states for a map of gauges
     */
    private void calcGaugeWstDurations(final WINFOArtifact winfo, final List<Gauge> gauges, final Map<Gauge, List<Double>> gaugeWstDurations,
            final Map<Gauge, GaugeDurationValuesFinder> durFinders) {

        final double[] gaugeKms = new double[gauges.size()];
        for (int i = 0; i <= gauges.size() - 1; i++) {
            gaugeKms[i] = gauges.get(i).getStation().doubleValue();
            gaugeWstDurations.put(gauges.get(i), new ArrayList<Double>());
        }
        final CalculationResult wstsData = winfo.computeWaterlevelData(gaugeKms);
        final WQKms[] wsts = (WQKms[]) wstsData.getData();
        for (int i = 0; i <= gauges.size() - 1; i++) {
            final GaugeDurationValuesFinder durFinder = durFinders.get(gauges.get(i));
            for (int j = 0; j <= wsts.length - 1; j++) {
                final double d = durFinder.getDuration(wsts[j].getQ(i));
                gaugeWstDurations.get(gauges.get(i)).add(Double.valueOf(underflowDaysToOverflowDays(d)));
            }
        }
    }

    /**
     * Create a result row for a station, and add w-q-values as selected
     */
    private ResultRow createRow(final Double station, final WQKms[] wqkmsArray, final List<Double> gaugeDurations, final int kmIndex) {

        final ResultRow row = ResultRow.create();
        row.putValue(GeneralResultType.station, station);
        row.putValue(SInfoResultType.infrastructuregroup, null); // is replaced later for an infrastructure type
        row.putValue(SInfoResultType.infrastructuretype, null); // is replaced later for an infrastructure part
        row.putValue(SInfoResultType.floodDuration, Double.NaN); // is replaced later for an infrastructure

        final String gaugeLabel = this.riverInfoProvider2.findGauge(station);
        row.putValue(GeneralResultType.gaugeLabel, gaugeLabel);

        final String location = this.riverInfoProvider2.getLocation(station);
        row.putValue(GeneralResultType.location, location);

        final List<DurationWaterlevel> wsts = new ArrayList<>(wqkmsArray.length);

        for (int i = 0; i <= wqkmsArray.length - 1; i++) {
            final DurationWaterlevel dw = new DurationWaterlevel(wqkmsArray[i].getW(kmIndex), gaugeDurations.get(i), wqkmsArray[i].getQ(kmIndex),
                    wqkmsArray[i].getName());
            wsts.add(dw);
        }
        row.putValue(SInfoResultType.customMultiRowColWaterlevel, wsts);

        return row;
    }

    /**
     * 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, final Set<FloodDurationInfrastructureChoice> infrastructures) {

        // 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;
        // Set the result row
        row.putValue(SInfoResultType.riverside, infrastructure.getAttributeKey());
        row.putValue(SInfoResultType.floodDischarge, q);
        row.putValue(SInfoResultType.infrastructureHeight, infrastructure.getHeight());
        row.putValue(SInfoResultType.infrastructuregroup, infrastructure.getInfrastructure().getGroup().getName());
        row.putValue(SInfoResultType.infrastructuretype, infrastructure.getInfrastructure().getType().getName());

        // Determine the relative column position of the Q of the infrastructure height
        final QPosition qPos = wst.getQPosition(infrastructure.getStation().doubleValue(), q);
        if (qPos != null) {
            // 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 D in the result row
            row.putValue(SInfoResultType.floodDuration, dur);
        }
        final FloodDurationInfrastructureChoice groupType = new FloodDurationInfrastructureChoice(row);
        infrastructures.add(groupType);
    }

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

http://dive4elements.wald.intevation.org