view flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/minfo/BedQualityCalculation.java @ 3867:d8c99ff18db8

Refactor MINFO bed quality resultsets and facets. flys-artifacts/trunk@5498 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author Raimund Renkert <raimund.renkert@intevation.de>
date Mon, 17 Sep 2012 14:38:47 +0000
parents 7fa38f8bcd8d
children 44e68dc1d87e
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.getBedMeasurements(
                    river,
                    from,
                    to,
                    dr.getFrom(),
                    dr.getTo());
            QualityMeasurements bedMeasurements =
                QualityMeasurementFactory.getBedMeasurements(
                    river,
                    from,
                    to,
                    dr.getFrom(),
                    dr.getTo());
            BedQualityResult result = new BedQualityResult();
            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