view flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/minfo/BedQualityCalculation.java @ 4187:21f4e4b79121

Refactor GaugeDischargeCurveFacet to be able to set a facet name For adding another output of the GaugeDischargeCurveArtifact it is necessary to provide to facet instances with different names. Therefore the GaugeDischargeCurveFacet is extended to set the facet name in the constructor.
author Björn Ricks <bjoern.ricks@intevation.de>
date Fri, 19 Oct 2012 13:25:49 +0200
parents 22cd60315e08
children ca18d8d06218
line wrap: on
line source
package de.intevation.flys.artifacts.model.minfo;

import gnu.trove.TDoubleArrayList;

import java.util.LinkedList;
import java.util.List;
import java.util.Map;

import org.apache.log4j.Logger;

import de.intevation.flys.artifacts.access.BedQualityAccess;
import de.intevation.flys.artifacts.model.Calculation;
import de.intevation.flys.artifacts.model.CalculationResult;
import de.intevation.flys.artifacts.model.DateRange;
import de.intevation.flys.backend.SedDBSessionHolder;


public class BedQualityCalculation extends Calculation {

    private static final Logger logger = Logger
        .getLogger(BedQualityCalculation.class);

    protected String river;
    protected double from;
    protected double to;
    protected List<String> bedDiameter;
    protected List<String> bedloadDiameter;
    protected List<DateRange> ranges;

    public BedQualityCalculation() {
    }

    public CalculationResult calculate(BedQualityAccess access) {
        logger.info("BedQualityCalculation.calculate");

        String river = access.getRiver();
        Double from = access.getFrom();
        Double to = access.getTo();
        List<String> bedDiameter = access.getBedDiameter();
        List<String> bedloadDiameter = access.getBedloadDiameter();
        List<DateRange> ranges = access.getDateRanges();

        if (river == null) {
            // TODO: i18n
            addProblem("minfo.missing.river");
        }

        if (from == null) {
            // TODO: i18n
            addProblem("minfo.missing.from");
        }

        if (to == null) {
            // TODO: i18n
            addProblem("minfo.missing.to");
        }

        if (ranges == null) {
            // TODO: i18n
            addProblem("minfo.missing.periods");
        }

        if (!hasProblems()) {
            this.river = river;
            this.from = from;
            this.to = to;
            this.ranges = ranges;
            this.bedDiameter = bedDiameter;
            this.bedloadDiameter = bedloadDiameter;

            SedDBSessionHolder.acquire();
            try {
                return internalCalculate();
            }
            finally {
                SedDBSessionHolder.release();
            }
        }

        return new CalculationResult();
    }

    protected CalculationResult internalCalculate() {

        List<BedQualityResult> results = new LinkedList<BedQualityResult>();
        // Calculate for all time periods.
        for (DateRange dr : ranges) {
            QualityMeasurements loadMeasurements =
                QualityMeasurementFactory.getBedloadMeasurements(
                    river,
                    from,
                    to,
                    dr.getFrom(),
                    dr.getTo());
            QualityMeasurements bedMeasurements =
                QualityMeasurementFactory.getBedMeasurements(
                    river,
                    from,
                    to,
                    dr.getFrom(),
                    dr.getTo());
            BedQualityResult result = new BedQualityResult();
            result.setDateRange(dr);
            if (bedDiameter != null) {
                result.add(calculateBedParameter(bedMeasurements, dr));
                for (String bd : bedDiameter) {
                    BedDiameterResult bedResult =
                        calculateBed(bedMeasurements, bd, dr);

                    // Avoid adding empty result sets.
                    if (!bedResult.isEmpty()) {
                        result.add(bedResult);
                    }
                }
            }
            if (bedloadDiameter != null) {
                for (String bld : bedloadDiameter) {
                    BedloadDiameterResult loadResult =
                        calculateBedload(loadMeasurements, bld, dr);
                    result.add(loadResult);
                }
            }
            results.add(result);
        }

        return new CalculationResult(
            results.toArray(new BedQualityResult[results.size()]), this);
    }

    private BedParametersResult calculateBedParameter(
        QualityMeasurements qm,
        DateRange dr
    ) {
        List<Double> kms = qm.getKms();
        QualityMeasurements capFiltered = filterCapMeasurements(qm);
        QualityMeasurements subFiltered = filterSubMeasurements(qm);
        TDoubleArrayList location = new TDoubleArrayList();
        TDoubleArrayList porosityCap = new TDoubleArrayList();
        TDoubleArrayList porositySub = new TDoubleArrayList();
        TDoubleArrayList densityCap = new TDoubleArrayList();
        TDoubleArrayList densitySub = new TDoubleArrayList();

        for(double km : kms) {
            double[] pCap = calculatePorosity(capFiltered, km);
            double[] pSub = calculatePorosity(subFiltered, km);
            double[] dCap = calculateDensity(capFiltered, pCap);
            double[] dSub = calculateDensity(subFiltered, pSub);

            double pCapRes = 0d;
            double pSubRes = 0d;
            double dCapRes = 0d;
            double dSubRes = 0d;
            for (int i = 0; i < pCap.length; i++) {
                pCapRes += pCap[i];
                dCapRes += dCap[i];
            }
            for (int i = 0; i < pSub.length; i++) {
                pSubRes += pSub[i];
                dSubRes += dSub[i];
            }
            location.add(km);
            porosityCap.add((pCapRes / pCap.length) * 100 );
            porositySub.add((pSubRes / pSub.length) * 100);
            densityCap.add((dCapRes / dCap.length) / 1000);
            densitySub.add((dSubRes / dSub.length) / 1000);

        }

        return new BedParametersResult(
            location,
            porosityCap,
            porositySub,
            densityCap,
            densitySub);
    }

    protected BedDiameterResult calculateBed(
        QualityMeasurements qm,
        String diameter,
        DateRange range
    ) {
        List<Double> kms = qm.getKms();
        TDoubleArrayList location = new TDoubleArrayList();
        TDoubleArrayList avDiameterCap = new TDoubleArrayList();
        TDoubleArrayList avDiameterSub = new TDoubleArrayList();
        for (double km : kms) {
            //Filter cap and sub measurements.
            QualityMeasurements capFiltered = filterCapMeasurements(qm);
            QualityMeasurements subFiltered = filterSubMeasurements(qm);

            List<QualityMeasurement> cm = capFiltered.getMeasurements(km);
            List<QualityMeasurement> sm = subFiltered.getMeasurements(km);

            double avCap = calculateAverage(cm, diameter);
            double avSub = calculateAverage(sm, diameter);
            location.add(km);
            avDiameterCap.add(avCap);
            avDiameterSub.add(avSub);
        }
        return new BedDiameterResult(
            diameter,
            avDiameterCap,
            avDiameterSub,
            location);
    }

    private double[] calculateDensity(
        QualityMeasurements capFiltered,
        double[] porosity
    ) {
        double[] density = new double[porosity.length];
        for (int i = 0; i < porosity.length; i++) {
            density[i] = (1 - porosity[i]) * 2650;
        }
        return density;
    }

    private double[] calculatePorosity(
        QualityMeasurements capFiltered,
        double km
    ) {
        List<QualityMeasurement> list = capFiltered.getMeasurements(km);
        double[] results = new double[list.size()];
        int i = 0;
        for (QualityMeasurement qm : list) {
            double deviation = calculateDeviation(qm);
            double p = calculateP(qm);
            double porosity = 0.353 - 0.068 * deviation + 0.146 * p;
            results[i] = porosity;
            i++;
        }

        return results;
    }

    protected BedloadDiameterResult calculateBedload(
        QualityMeasurements qm,
        String diameter,
        DateRange range
    ) {
        List<Double> kms = qm.getKms();
        TDoubleArrayList location = new TDoubleArrayList();
        TDoubleArrayList avDiameter = new TDoubleArrayList();
        for (double km : kms) {
            List<QualityMeasurement> measurements = qm.getMeasurements(km);
            double mid = calculateAverage(measurements, diameter);
            location.add(km);
            avDiameter.add(mid);
        }
        return new BedloadDiameterResult(
            diameter,
            avDiameter,
            location,
            range);
    }

    protected double calculateAverage(
        List<QualityMeasurement> list,
        String diameter
    ) {
        double av = 0;
        for (QualityMeasurement qm : list) {
            av += qm.getDiameter(diameter);
        }
        return av/list.size();
    }

    protected QualityMeasurements filterCapMeasurements(
        QualityMeasurements qms
    ) {
        List<QualityMeasurement> result = new LinkedList<QualityMeasurement>();
        for (QualityMeasurement qm : qms.getMeasurements()) {
            if (qm.getDepth1() == 0d && qm.getDepth2() <= 0.3) {
                result.add(qm);
            }
        }
        return new QualityMeasurements(result);
    }

    protected QualityMeasurements filterSubMeasurements(
        QualityMeasurements qms
    ) {
        List<QualityMeasurement> result = new LinkedList<QualityMeasurement>();
        for (QualityMeasurement qm : qms.getMeasurements()) {
            if (qm.getDepth1() > 0d && qm.getDepth2() <= 0.5) {
                result.add(qm);
            }
        }
        return new QualityMeasurements(result);
    }

    public double calculateDeviation(QualityMeasurement qm) {
        Map<String, Double> dm = qm.getAllDiameter();
        double phiM = 0;
        double[] phis = new double[dm.size()];
        double[] ps = new double[dm.size()];
        int i = 0;
        for (String key : dm.keySet()) {
            double d = dm.get(key);
            double phi = -Math.log(d)/Math.log(2);
            phis[i] = phi;
            double p = calculateWeight(qm, key);
            ps[i] = p;
            phiM += phi * p;
            i++;
        }

        double sig = 0d;
        for (i = 0; i < dm.size(); i++) {
            sig += ps[i] * Math.exp(phis[i] - phiM);
        }
        double deviation = Math.sqrt(sig);
        return deviation;
    }

    protected double calculateP(QualityMeasurement qm) {
        return calculateWeight(qm, "dmin");
    }

    public double calculateWeight(QualityMeasurement qm, String diameter) {
        Map<String, Double> dm = qm.getAllDiameter();
        double value = qm.getDiameter(diameter);

        double sum = 0d;
        for (Double d : dm.values()) {
            sum =+ d.doubleValue();
        }
        double weight = sum/100*value;
        return weight;
    }
}

http://dive4elements.wald.intevation.org