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

Implemented S-INFO Flowdepth TKH calculation
author mschaefer
date Thu, 22 Feb 2018 12:03:31 +0100
parents
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.flowdepth;

import java.util.ArrayList;
import java.util.List;

import org.apache.commons.lang.math.DoubleRange;
import org.apache.log4j.Logger;
import org.hibernate.SQLQuery;
import org.hibernate.Session;
import org.hibernate.type.StandardBasicTypes;

import gnu.trove.TDoubleArrayList;

import org.dive4elements.river.artifacts.math.Linear;
import org.dive4elements.river.artifacts.math.Utils;
import org.dive4elements.river.backend.SessionHolder;
import org.dive4elements.river.model.River;

/**
 * Searchable sorted km array with parallel FlowVelocityKmModelValues array and linear interpolation for km and the model values between the array elements.<br />
 * {@link loadValues} loads all the model values for a given km range of a river.<br />
 * {@link findKmQValues} then searches a km in the values table or the nearest including km interval, resp.
 * The v and tau values for a given discharge are either found directly or also interpolated linearly.<br />
 * 
 * (Created based on a copy of FlowVelocityMeasurementFactory.)
 * 
 * @author Matthias Schäfer
 *
 */
public class FlowVelocityModelKmValueFinder
{
    /***** FIELDS *****/
    
    /**
     * Private log to use here.
     */
    private static Logger log = Logger.getLogger(FlowVelocityModelKmValueFinder.class);
    
    /**
     * Query for a range of stations of a river with all their q, main-v and tau values.<br />
     * (Might be several 10000 rows if many stations and large q range)
     */
    private static final String SQL_SELECT_ALL =
        "SELECT fvmv.station AS station, fvmv.q AS q, fvmv.main_channel AS vmain, fvmv.shear_stress AS tau" 
        + "  FROM (discharge_zone dz INNER JOIN flow_velocity_model fvm ON dz.id = fvm.discharge_zone_id)" 
        + "    INNER JOIN flow_velocity_model_values fvmv ON fvm.id = fvmv.flow_velocity_model_id" 
        + "  WHERE (dz.river_id = :river_id) AND (fvmv.station BETWEEN :kmfrom - 0.0001 AND :kmto + 0.0001)" 
        /* + "  WHERE (dz.river_id = :river_id) AND (fvmv.q BETWEEN :qmin AND :qmax)" */
        + "  ORDER BY fvmv.station ASC, fvmv.q ASC";

    /**
     * Query for a river's max km below a limit with all its q, main-v and tau values.
     */
    private static final String SQL_SELECT_KMLOWER =
        "SELECT fvmv.station AS station, fvmv.q AS q, fvmv.main_channel AS vmain, fvmv.shear_stress AS tau"
        + "  FROM flow_velocity_model_values fvmv"
        + "    INNER JOIN (SELECT MAX(fvmvi.station) AS kmmax"
        + "      FROM(discharge_zone dz INNER JOIN flow_velocity_model fvm ON dz.id = fvm.discharge_zone_id)"
        + "        INNER JOIN flow_velocity_model_values fvmvi ON fvm.id = fvmvi.flow_velocity_model_id"
        + "        WHERE (dz.river_id = :river_id) AND (fvmvi.station < :kmfrom - 0.0001)) finf ON fvmv.station = finf.kmmax" 
        + "  ORDER BY fvmv.q ASC";

    /**
     * Query for a river's min km above a limit with all its q, main-v and tau values.
     */
    private static final String SQL_SELECT_KMUPPER =
        "SELECT fvmv.station AS station, fvmv.q AS q, fvmv.main_channel AS vmain, fvmv.shear_stress AS tau" 
        + "  FROM flow_velocity_model_values fvmv"
        + "    INNER JOIN (SELECT MIN(fvmvi.station) AS kmmin"
        + "      FROM(discharge_zone dz INNER JOIN flow_velocity_model fvm ON dz.id = fvm.discharge_zone_id)" 
        + "        INNER JOIN flow_velocity_model_values fvmvi ON fvm.id = fvmvi.flow_velocity_model_id"
        + "        WHERE (dz.river_id = :river_id) AND (fvmvi.station > :kmto + 0.0001)) fsup ON fvmv.station = fsup.kmmin" 
        + "  ORDER BY fvmv.q ASC";

    /**
     * Query to select all km-q-v-tau of a river that are the q maxima below a q limit
     */
    private static final String SQL_SELECT_QLOWER =
        "SELECT fvmv.station AS station, fvmv.q AS q, fvmv.main_channel AS vmain, fvmv.shear_stress AS tau"
        + "  FROM flow_velocity_model_values fvmv"
        + "    INNER JOIN (SELECT fv2.station, MAX(fv2.q) AS q"
        + "      FROM (discharge_zone dz INNER JOIN flow_velocity_model fvm ON dz.id = fvm.discharge_zone_id)" 
        + "        INNER JOIN flow_velocity_model_values fv2 ON fvm.id = fv2.flow_velocity_model_id" 
        + "      WHERE (dz.river_id = :river_id) AND (fv2.q < :qlim) GROUP BY fv2.station) qx"
        + "    ON (fvmv.station=qx.station) AND (fvmv.q=qx.q)" 
        + "  ORDER BY fvmv.station ASC";
        
    /**
     * Query to select all km-q-v-tau of a river that are the q minima above a q limit
     */
    private static final String SQL_SELECT_QUPPER =
            "SELECT fvmv.station AS station, fvmv.q AS q, fvmv.main_channel AS vmain, fvmv.shear_stress AS tau"
            + "  FROM flow_velocity_model_values fvmv"
            + "    INNER JOIN (SELECT fv2.station, MIN(fv2.q) AS q"
            + "      FROM (discharge_zone dz INNER JOIN flow_velocity_model fvm ON dz.id = fvm.discharge_zone_id)" 
            + "        INNER JOIN flow_velocity_model_values fv2 ON fvm.id = fv2.flow_velocity_model_id" 
            + "      WHERE (dz.river_id = :river_id) AND (fv2.q > :qlim) GROUP BY fv2.station) qx"
            + "    ON (fvmv.station=qx.station) AND (fvmv.q=qx.q)" 
            + "  ORDER BY fvmv.station ASC";
            
    /**
     * Kms of the loaded river range
     */
    private TDoubleArrayList kms;
    
    /**
     * For each km in kms a list of q-v-tau-tupels
     */
    private List<FlowVelocityKmModelValues> values;

    /**
     * Searched km of the last findKmValue
     */
    private double findKm;
    
    /**
     * kms and values index of the interval start found by the last findKmValue
     */
    private int leftIndexFound = -1;

    /**
     * kms and values index of the interval end found by the last findKmValue
     */
    private int rightIndexFound = -1;

    /**
     * Q of the last findKmQValues
     */
    private double findQ;

    /***** METHODS *****/
    
    /**
     * Discharge of the last {@link findKmQValue}
     */
    public double getFindQ() {
        return findQ;
    }
    
    /**
     * Velocity of the last {@link findKmQValues}
     */
    public double getFindVmainFound() {
        if (leftIndexFound < 0)
            return Double.NaN;
        else if (leftIndexFound == rightIndexFound)
            return getLeftValues().getVmainFound();
        else
            return Linear.linear(findKm, getLeftValues().getKm(), getRightValues().getKm(), getLeftValues().getVmainFound(), getRightValues().getVmainFound());
    }
    
    /**
     * Shear stress tau of the last {@link findKmQValues}
     */
    public double getFindTauFound() {
        if (leftIndexFound < 0)
            return Double.NaN;
        else if (leftIndexFound == rightIndexFound)
            return getLeftValues().getTauFound();
        else
            return Linear.linear(findKm, getLeftValues().getKm(), getRightValues().getKm(), getLeftValues().getTauFound(), getRightValues().getTauFound());
    }
    
    /**
     * Whether the discharge has been interpolated in the last {@link findKmQValues} 
     */
    public boolean getFindIsQInterpolated() {
        return (getLeftValues() != null) && (getLeftValues().getIsInterpolated() || getRightValues().getIsInterpolated());
    }
    
    /**
     * Queries a range of a river's kms with all their q-v-tau values.
     * @return Whether the load has been successful
     */
    @SuppressWarnings("unchecked")
    public boolean loadValues(River river, DoubleRange kmRange, DoubleRange qRange) {
        // DB session
        log.debug(String.format("loadValues km %.3f - %.3f / q %.1f - %.1f", kmRange.getMinimumDouble(), kmRange.getMaximumDouble(), qRange.getMinimumDouble(), qRange.getMaximumDouble()));
        kms = new TDoubleArrayList();
        values = new ArrayList<FlowVelocityKmModelValues>();
        boolean isDemoValuesCorrection = river.getName().equalsIgnoreCase("beispielfluss");
        final Session session = SessionHolder.HOLDER.get();
        
        // Select km infimum
        SQLQuery sqlQuery = session.createSQLQuery(SQL_SELECT_KMLOWER)
            .addScalar("station", StandardBasicTypes.DOUBLE)
            .addScalar("q", StandardBasicTypes.DOUBLE)
            .addScalar("vmain", StandardBasicTypes.DOUBLE)
            .addScalar("tau", StandardBasicTypes.DOUBLE);
        sqlQuery.setParameter("river_id", river.getId());
        sqlQuery.setParameter("kmfrom", kmRange.getMinimumDouble());
        addKms(sqlQuery.list(), isDemoValuesCorrection);

        // Select km range
        sqlQuery = session.createSQLQuery(SQL_SELECT_ALL)
            .addScalar("station", StandardBasicTypes.DOUBLE)
            .addScalar("q", StandardBasicTypes.DOUBLE)
            .addScalar("vmain", StandardBasicTypes.DOUBLE)
            .addScalar("tau", StandardBasicTypes.DOUBLE);
        sqlQuery.setParameter("river_id", river.getId());
        sqlQuery.setParameter("kmfrom", kmRange.getMinimumDouble());
        sqlQuery.setParameter("kmto", kmRange.getMaximumDouble());
        //sqlQuery.setParameter("qmin", qRange.getMinimumDouble());
        //sqlQuery.setParameter("qmax", qRange.getMaximumDouble());
        int kmcount = kms.size();
        int rowcount = addKms(sqlQuery.list(), isDemoValuesCorrection);
        kmcount = kms.size() - kmcount;
        
        // Select km supremum 
        sqlQuery = session.createSQLQuery(SQL_SELECT_KMUPPER)
            .addScalar("station", StandardBasicTypes.DOUBLE)
            .addScalar("q", StandardBasicTypes.DOUBLE)
            .addScalar("vmain", StandardBasicTypes.DOUBLE)
            .addScalar("tau", StandardBasicTypes.DOUBLE);
        sqlQuery.setParameter("river_id", river.getId());
        sqlQuery.setParameter("kmto", kmRange.getMaximumDouble());
        int supcnt = addKms(sqlQuery.list(), isDemoValuesCorrection);

        // Add copy of last km for search of max km value
        if ((supcnt == 0) && (values.size() >= 1)) {
            kms.add(kms.getQuick(kms.size()) + 0.0001);
            values.add(new FlowVelocityKmModelValues(kms.getQuick(kms.size()-1), values.get(values.size()-1)));
        }
        
        // log.debug
        if (values.size() - 1 >= 0) {
            log.debug(String.format("loadValues %d: km %.3f - %d values", 0, values.get(0).getKm(), values.get(0).size()));
            if (values.size() - 1 >= 1) {
                log.debug(String.format("loadValues %d: km %.3f - %d values", 1, values.get(1).getKm(), values.get(1).size()));
                if (values.size() - 1 >= 2)
                    log.debug("loadValues ...");
                if (values.size() - 2 >= 3)
                    log.debug(String.format("loadValues %d: km %.3f - %d values", values.size()-2, values.get(values.size()-2).getKm(), values.get(values.size()-2).size()));
                if (values.size() - 1 >= 3)
                    log.debug(String.format("loadValues %d: km %.3f - %d values", values.size()-1, values.get(values.size()-1).getKm(), values.get(values.size()-1).size()));
            }
        }
        log.debug(String.format("loadValues %d kms, %d values loaded", kmcount, rowcount));
        return (kms.size() >= 1);
    }
    
    /**
     * Adds the km-q-v-tau values of a query result row to the last km of the list, or a new one resp.
     * @return Number of rows
     */
    private int addKms(List<Object[]> rows, boolean isDemoValuesCorrection) {
        for (Object[] row : rows) {
            if ((kms.size() == 0) || !Utils.epsilonEquals(kms.get(kms.size()-1), (double) row[0], 0.0001)) {
                kms.add((double) row[0]);
                values.add(new FlowVelocityKmModelValues(kms.get(kms.size()-1)));
            }
            if (isDemoValuesCorrection)
                // "Verfremdung" der v-Werte etwas korrigieren (Originalwerte wurden mit Zufallswert zwischen 10 und 20 multipliziert)
                values.get(values.size()-1).addValues((double) row[1], ((double) row[2]) / 10, (double) row[3]);
            else
                values.get(values.size()-1).addValues((double) row[1], (double) row[2], (double) row[3]);
        }
        return rows.size();
    }
    
    /**
     * Searches a km and finds or interpolates the velocity and shear stress values for a discharge<br />
     * The values may be got via {@link getVmainFound} etc.
     * @return Whether values have been found
     */
    public boolean findKmQValues(double km, double q) {
        findQ = q;
        if (!searchKm(km))
            return false;
        if (leftIndexFound == rightIndexFound) {
            // Exact km match
            final double qfound = getLeftValues().findQ(q);
            log.debug(String.format("findKmQValues km %.3f q %.0f = %.0f (%d)", km, q, qfound, leftIndexFound));
            return !Double.isNaN(qfound);
        }
        else {
            final double[] qfound = {getLeftValues().findQ(q), getRightValues().findQ(q)};
            log.debug(String.format("findKmQValues km %.3f q %.0f = %.0f (%d, %.3f) - %.0f (%d, %.3f)", km, q, qfound[0], leftIndexFound,
                    getLeftValues().getKm(), qfound[1], rightIndexFound, getRightValues().getKm()));
            return !Double.isNaN(qfound[0]) && !Double.isNaN(qfound[1]);
        }
    }
    
    /**
     * Searches a km
     * @return Whether the km was within the supported range
     */
    private boolean searchKm(double km) {
        findKm = km;
        leftIndexFound = -1;
        rightIndexFound = -1;
        if ((kms == null) || (kms.size() == 0))
            return false;
        int i = kms.binarySearch(km);
        if (i >= 0) {
            // Exact km match
            leftIndexFound = i;
            rightIndexFound = i;
            return true;
        }
        else {
            // Out of range or within km interval
            if (i < 0)
                i = -i - 1;
            if ((i <= 0) || (i >= kms.size()))
                return false;
            leftIndexFound = i - 1;
            rightIndexFound = i;
            return true;
        }
    }
    
    private FlowVelocityKmModelValues getLeftValues() {
        return values.get(leftIndexFound);
    }
    private FlowVelocityKmModelValues getRightValues() {
        return values.get(rightIndexFound);
    }

}

http://dive4elements.wald.intevation.org