view artifacts/src/main/java/org/dive4elements/river/artifacts/sinfo/tkhcalculation/BedQualityD50KmValueFinder.java @ 8980:b194fa64506a

SINFO - show results themes according to spec, either raw data or floating mean values. Some improvements to error handling and handling of empty results.
author gernotbelger
date Thu, 05 Apr 2018 18:30:34 +0200
parents 45f1ad66560e
children 9b9f5f4ddb80
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.tkhcalculation;

import java.util.Calendar;
import java.util.Date;
import java.util.List;

import org.apache.commons.lang.math.DoubleRange;
import org.apache.commons.math.ArgumentOutsideDomainException;
import org.apache.commons.math.analysis.interpolation.LinearInterpolator;
import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction;
import org.apache.log4j.Logger;
import org.dive4elements.river.artifacts.math.Utils;
import org.dive4elements.river.artifacts.model.Calculation;
import org.dive4elements.river.backend.SedDBSessionHolder;
import org.dive4elements.river.model.River;
import org.hibernate.SQLQuery;
import org.hibernate.Session;
import org.hibernate.type.StandardBasicTypes;

import gnu.trove.TDoubleArrayList;

/**
 * Searchable sorted km array with parallel bed measurements value array and linear interpolation for km and d50 between
 * the array elements.<br />
 * <br />
 * See comment of SQL command on how the values are filtered and aggregated.
 *
 * @author Matthias Schäfer
 *
 */
public class BedQualityD50KmValueFinder {

    /***** FIELDS *****/

    /**
     * Private log to use here.
     */
    private static Logger log = Logger.getLogger(BedQualityD50KmValueFinder.class);

    /**
     * Query selecting all sub layer bed measurements with their d50 for a km range and a time period<br />
     * <br />
     * A km may have bed measurements for multiple dates, multiple distances from the river bank, and multiple depth layers.
     * The query filters by km range, time period and layer (sub layer: below bed to max. 50 cm depth).<br />
     *
     * If PostgreSQL would support a median aggregate function like Oracle does, the aggregation could be placed into this
     * query.
     */
    private static final String SQL_BED_D50_SUBLAYER_MEASUREMENT = "SELECT t.km, t.datum, p.tiefevon, p.tiefebis, a.d50"
            + " FROM sohltest t INNER JOIN station s ON t.stationid = s.stationid" + "    INNER JOIN gewaesser g ON s.gewaesserid = g.gewaesserid"
            + "    INNER JOIN sohlprobe p ON t.sohltestid = p.sohltestid" + "    INNER JOIN siebanalyse a ON p.sohlprobeid = a.sohlprobeid"
            + " WHERE (g.name = :name) AND (s.km BETWEEN :fromkm - 0.0001 AND :tokm + 0.0001)" + "    AND (p.tiefevon > 0.0) AND (p.tiefebis <= 0.5)"
            + "    AND (t.datum BETWEEN :fromdate AND :todate)" + " ORDER BY t.km ASC, a.d50 ASC";

    private Calculation problems;

    /**
     * Real linear interpolator for kms and d50 values (m)
     */
    private final PolynomialSplineFunction interpolator;

    /***** CONSTRUCTORS *****/

    private BedQualityD50KmValueFinder(final Calculation problems, final double[] kms, final double[] values) {
        this.problems = problems;

        // FIXME: check: max distance prüfen? dann D4E-LinearInterpolator verwenden
        this.interpolator = new LinearInterpolator().interpolate(kms, values);
    }

    /***** METHODS *****/

    /**
     * Sohlbeschaffenheit (D50 Korndurchmesser aus Seddb)
     * Abhängig von Peiljahr
     *
     * @param problems
     */
    public static BedQualityD50KmValueFinder loadBedMeasurements(final Calculation problems, final River river, final DoubleRange kmRange,
            final int soundingYear, final int validYears) {

        /* construct valid measurement time range */
        final Calendar cal = Calendar.getInstance();
        cal.clear();

        cal.set(soundingYear - validYears, 0, 1);
        final Date startTime = cal.getTime();

        cal.set(soundingYear + validYears, 11, 31);
        final Date endTime = cal.getTime();

        log.debug(String.format("loadValues km %.3f - %.3f %tF - %tF", kmRange.getMinimumDouble(), kmRange.getMaximumDouble(), startTime, endTime));
        final Session session = SedDBSessionHolder.HOLDER.get();
        final SQLQuery sqlQuery = session.createSQLQuery(SQL_BED_D50_SUBLAYER_MEASUREMENT).addScalar("km", StandardBasicTypes.DOUBLE)
                .addScalar("datum", StandardBasicTypes.DATE).addScalar("tiefevon", StandardBasicTypes.DOUBLE).addScalar("tiefebis", StandardBasicTypes.DOUBLE)
                .addScalar("d50", StandardBasicTypes.DOUBLE);
        final String seddbRiver = river.nameForSeddb();
        sqlQuery.setString("name", seddbRiver);
        sqlQuery.setDouble("fromkm", kmRange.getMinimumDouble());
        sqlQuery.setDouble("tokm", kmRange.getMaximumDouble());
        sqlQuery.setDate("fromdate", startTime);
        sqlQuery.setDate("todate", endTime);

        final List<Object[]> rows = sqlQuery.list();
        final TDoubleArrayList kms = new TDoubleArrayList();
        final TDoubleArrayList values = new TDoubleArrayList();
        final TDoubleArrayList kmd50s = new TDoubleArrayList();

        for (int i = 0; i <= rows.size() - 1; i++) {
            kmd50s.add((double) rows.get(i)[4]);
            if (((i == rows.size() - 1) || !Utils.epsilonEquals((double) rows.get(i)[0], (double) rows.get(i + 1)[0], 0.0001))) {
                final int k = kmd50s.size() / 2;
                values.add(((k + k < kmd50s.size()) ? kmd50s.get(k) : (kmd50s.get(k - 1) + kmd50s.get(k)) / 2) / 1000);
                kms.add((double) rows.get(i)[0]);
                log.debug(String.format("loadValues km %.3f d50(mm) %.1f count %d", kms.get(kms.size() - 1), values.get(values.size() - 1), kmd50s.size()));
                kmd50s.clear();
            }
        }

        if (kms.size() < 2 || values.size() < 2) {
            problems.addProblem("bedqualityd50kmvaluefinder.empty", soundingYear);
            return null;
        }

        try {
            return new BedQualityD50KmValueFinder(problems, kms.toNativeArray(), values.toNativeArray());
        }
        catch (final Exception e) {
            e.printStackTrace();
            problems.addProblem("bedqualityd50kmvaluefinder.error", e.getLocalizedMessage());
            return null;
        }
    }

    /**
     * Returns the d50 value interpolated according to a km
     *
     * @return d50 (m) of the km, or NaN
     */
    public double findD50(final double km) {
        try {
            return this.interpolator.value(km);
        }
        catch (final ArgumentOutsideDomainException e) {
            e.printStackTrace();

            if (this.problems != null) {
                this.problems.addProblem(km, "bedqualityd50kmvaluefinder.missing");
                // Report only once
                this.problems = null;
            }

            return Double.NaN;
        }
    }
}

http://dive4elements.wald.intevation.org