view artifacts/src/main/java/org/dive4elements/river/artifacts/model/minfo/SedimentLoadCalculation.java @ 6665:b7945db8a43b

issue1413: Only show unknown sediment loads of selected unit type. Therefore, adjusted the factory to take the units name. Unfortunately, names in db do not match values of data items. Thus do manual replacing. In Facet and Calculate, take the chosen unit via access and to the string replacement. In Facet, do not transform data (we assume it comes in unit as labeled in the db), and removed the possibility of m3/a-data of unknown yields in a t/a diagram and vice versa.
author Felix Wolfsteller <felix.wolfsteller@intevation.de>
date Thu, 25 Jul 2013 15:08:13 +0200
parents 1ca0688dddc7
children 0a70a320bfca
line wrap: on
line source
/* Copyright (C) 2011, 2012, 2013 by Bundesanstalt für Gewässerkunde
 * Software engineering by Intevation GmbH
 *
 * 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.model.minfo;

import gnu.trove.TDoubleArrayList;

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

import org.apache.log4j.Logger;

import org.dive4elements.river.artifacts.access.SedimentLoadAccess;
import org.dive4elements.river.artifacts.model.Calculation;
import org.dive4elements.river.artifacts.model.CalculationResult;
import org.dive4elements.river.artifacts.model.Range;


/** Calculate sediment load. */
public class SedimentLoadCalculation
extends Calculation
{

    /** Private logger. */
    private static final Logger logger = Logger
        .getLogger(SedimentLoadCalculation.class);

    protected String river;
    protected String yearEpoch;
    protected double kmUp;
    protected double kmLow;
    protected int[] period;
    protected int[][] epoch;
    protected String unit;

    public SedimentLoadCalculation() {
    }

    public CalculationResult calculate(SedimentLoadAccess access) {
        logger.info("SedimentLoadCalculation.calculate");

        String river = access.getRiver();
        String yearEpoch = access.getYearEpoch();
        String unit = access.getUnit();
        int[] period = null;
        int[][] epoch = null;
        double kmUp = access.getUpperKM();
        double kmLow = access.getLowerKM();
        if (yearEpoch.equals("year")) {
            period = access.getPeriod();
            epoch = null;
        }
        else if (yearEpoch.equals("epoch") || yearEpoch.equals("off_epoch")) {
            epoch = access.getEpochs();
            period = null;
        }
        else {
            addProblem("minfo.missing.year_epoch");
        }

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

        if (period == null && epoch == null) {
            addProblem("minfo.missing.time");
        }

        if (!hasProblems()) {
            this.river = river;
            this.yearEpoch = yearEpoch;
            this.unit = unit;
            this.period = period;
            this.epoch = epoch;
            this.kmUp = kmUp;
            this.kmLow = kmLow;
            return internalCalculate();
        }

        return new CalculationResult();
    }

    private CalculationResult internalCalculate() {
        logger.debug("internalCalulate; mode:" + yearEpoch);
        if (yearEpoch.equals("year")) {
            List<SedimentLoadResult> results =
                new ArrayList<SedimentLoadResult>();
            for (int i = period[0]; i <= period[1]; i++) {
                SedimentLoadResult res = calculateYear(i);
                results.add(res);
            }
            return new CalculationResult(
                results.toArray(new SedimentLoadResult[results.size()]), this);
        }
        else if (yearEpoch.equals("epoch")) {
            List<SedimentLoadResult> results =
                new ArrayList<SedimentLoadResult>();
            for (int i = 0; i < epoch.length; i++) {
                SedimentLoadResult res = calculateEpoch(i);
                results.add(res);
            }
            return new CalculationResult(
                results.toArray(new SedimentLoadResult[results.size()]), this);
        }
        else if (yearEpoch.equals("off_epoch")) {
            List<SedimentLoadResult> results =
                new ArrayList<SedimentLoadResult>();
            for (int i = 0; i < epoch.length; i++) {
                SedimentLoadResult res = calculateOffEpoch(i);
                results.add(res);
            }
            return new CalculationResult(
                results.toArray(new SedimentLoadResult[results.size()]), this);
        }
        return null;
    }

    private SedimentLoadResult calculateEpoch(int i) {
        List<SedimentLoad> epochLoads = new ArrayList<SedimentLoad>();
        for (int j = epoch[i][0]; j < epoch[i][1]; j++) {
            epochLoads.add(SedimentLoadFactory.getLoadWithData(
                this.river,
                this.yearEpoch,
                this.kmLow,
                this.kmUp,
                j,
                j));
        }

        SedimentLoad resLoad = new SedimentLoad();
        TDoubleArrayList kms = new TDoubleArrayList();

        for (SedimentLoad load : epochLoads) {
            for (double km : load.getKms()) {
                if (!kms.contains(km)) {
                    kms.add(km);
                }
            }
        }

        for (int j = 0; j < kms.size(); j++) {
            int cSum = 0;
            int fmSum = 0;
            int sSum = 0;
            int ssSum = 0;
            int ssbSum = 0;
            int sseSum = 0;
            double km = kms.get(j);
            for (SedimentLoad load : epochLoads) {
                SedimentLoadFraction f = load.getFraction(km);
                if (f.getCoarse() > 0d) {
                    double c = resLoad.getFraction(km).getCoarse();
                    resLoad.setCoarse(km, c + f.getCoarse(), f.getCoarseRange());
                    cSum++;
                }
                if (f.getFineMiddle() > 0d) {
                    double fm = resLoad.getFraction(km).getFineMiddle();
                    resLoad.setFineMiddle(km, fm + f.getFineMiddle(), f.getFineMiddleRange());
                    fmSum++;
                }
                if (f.getSand() > 0d) {
                    double s = resLoad.getFraction(km).getSand();
                    resLoad.setSand(km, s + f.getSand(), f.getSandRange());
                    sSum++;
                }
                if (f.getSuspSand() > 0d) {
                    double s = resLoad.getFraction(km).getSuspSand();
                    resLoad.setSuspSand(km, s + f.getSuspSand(), f.getSuspSandRange());
                    ssSum++;
                }
                if (f.getSuspSandBed() > 0d) {
                    double s = resLoad.getFraction(km).getSuspSandBed();
                    resLoad.setSuspSandBed(km, s + f.getSuspSandBed(), f.getSuspSandBedRange());
                    ssbSum++;
                }
                if (f.getSuspSediment() > 0d) {
                    double s = resLoad.getFraction(km).getSuspSediment();
                    resLoad.setSuspSediment(km, s + f.getSuspSediment(), f.getSuspSedimentRange());
                    sseSum++;
                }
            }
            SedimentLoadFraction fr = resLoad.getFraction(km);
            // Prevent divisions by zero, the fraction defaults to 0d.
            if (cSum != 0) {
                resLoad.setCoarse(km, fr.getCoarse()/cSum, fr.getCoarseRange());
            }
            if (fmSum != 0) {
                resLoad.setFineMiddle(km, fr.getFineMiddle()/fmSum,
                    fr.getFineMiddleRange());
            }
            if (sSum != 0) {
                resLoad.setSand(km, fr.getSand()/sSum, fr.getSandRange());
            }
            if (ssSum != 0) {
                resLoad.setSuspSand(km, fr.getSuspSand()/ssSum,
                    fr.getSuspSandRange());
            }
            if (ssbSum != 0) {
                resLoad.setSuspSandBed(km, fr.getSuspSandBed()/ssbSum,
                    fr.getSuspSandBedRange());
            }
            if (sseSum != 0) {
                resLoad.setSuspSediment(km, fr.getSuspSediment()/sseSum, fr.getSuspSedimentRange());
            }
        }
        resLoad.setDescription("");
        resLoad.setEpoch(true);

        SedimentLoadResult result;
        SedimentLoad sl = calculateTotalLoad(resLoad, this.epoch[i][0]);
        if (this.unit.equals("m3_per_a")) {
            SedimentLoad slu = calculateUnit(sl, this.epoch[i][0]);
            result = new SedimentLoadResult(
                this.epoch[i][0],
                this.epoch[i][1],
                slu);
        }
        else {
            result = new SedimentLoadResult(
                this.epoch[i][0],
                this.epoch[i][1],
                sl);
        }

        return result;
    }

    private SedimentLoadResult calculateOffEpoch(int i) {
        SedimentLoad load = SedimentLoadFactory.getLoadWithData(
            this.river,
            this.yearEpoch,
            kmLow,
            kmUp,
            this.epoch[i][0],
            this.epoch[i][1]);
        SedimentLoadResult result;
        SedimentLoad sl = calculateTotalLoad(load, this.epoch[i][0]);
        if (unit.equals("m3_per_a")) {
            SedimentLoad slu = calculateUnit(sl, epoch[i][0]);
            result = new SedimentLoadResult(
                this.epoch[i][0],
                this.epoch[i][1],
                slu);
        }
        else {
            result = new SedimentLoadResult(
                this.epoch[i][0],
                this.epoch[i][1],
                sl);
        }

        return result;
    }

    /** Fetch loads for a single year, calculate total and
     * return the result containing both. */
    private SedimentLoadResult calculateYear(int y) {
        SedimentLoad load = SedimentLoadFactory.getLoadWithData(
            this.river,
            this.yearEpoch,
            this.kmLow,
            this.kmUp,
            y,
            y);

        SedimentLoadResult result;
        SedimentLoad sl = calculateTotalLoad(load, y);
        if (unit.equals("m3_per_a")) {
            SedimentLoad slu = calculateUnit(sl, y);
            result = new SedimentLoadResult(y, 0, slu);
        }
        else {
            result = new SedimentLoadResult(y, 0, sl);
        }
        return result;
    }

    /** Add up the loads of a year. */
    private SedimentLoad calculateTotalLoad(SedimentLoad load, int year) {
        logger.debug("calculateTotalLoad");
        boolean problemThisYear = false;
        if (!load.hasCoarse()) {
            addProblem("missing.fraction.coarse", Integer.toString(year));
            problemThisYear = true;
        }
        if (!load.hasFineMiddle()) {
            addProblem("missing.fraction.fine_middle", Integer.toString(year));
            problemThisYear = true;
        }
        if (!load.hasSand()) {
            addProblem("missing.fraction.sand", Integer.toString(year));
            problemThisYear = true;
        }
        if (!load.hasSuspSand()) {
            addProblem("missing.fraction.susp_sand", Integer.toString(year));
            problemThisYear = true;
        }
        if (!load.hasSuspSediment()) {
            addProblem("missing.fraction.susp_sediment", Integer.toString(year));
            problemThisYear = true;
        }
        if (problemThisYear) {
            logger.warn("Some problem, not calculating total load.");
            return load;
        }
        return partialTotal(load);
    }


    /** Returns true if all fraction values except SuspSediment are unset. */
    private boolean hasOnlySuspValues(SedimentLoadFraction fraction) {
        return (fraction.getSuspSediment() != 0d &&
            fraction.getCoarse() == 0d &&
            fraction.getFineMiddle() == 0d &&
            fraction.getSand() == 0d &&
            fraction.getSuspSand() == 0d);
    }


    /** Returns true if all fraction values except SuspSediment are set. */
    private boolean hasButSuspValues(SedimentLoadFraction fraction) {
        return (fraction.getSuspSediment() == 0d &&
            fraction.getCoarse() != 0d &&
            fraction.getFineMiddle() != 0d &&
            fraction.getSand() != 0d &&
            fraction.getSuspSand() != 0d);
    }


    /** Returns true if all fraction needed for total calculation are set. */
    private boolean complete(SedimentLoadFraction fraction) {
        return (fraction.getCoarse() != 0d &&
                fraction.getFineMiddle() != 0d &&
                fraction.getSand() != 0d &&
                fraction.getSuspSand() != 0d &&
                fraction.getSuspSediment() != 0d);
    }


    /**
     * Set total values in load.
     * Therefore, run over the kms and find ranges where either all
     * or all Geschiebe or just the Schwebstoff fractions are set.
     * Merge these ranges and add (maybe new) respective fractions to
     * load.
     * @param load SedimentLoad to add total values (and ranges) to.
     * @return input param load.
     */
    private SedimentLoad partialTotal(SedimentLoad load) {
        SedimentLoad fairLoad = load;

        Range lastOtherRange = null;
        double lastOtherValue = 0d;

        Range lastSuspRange = null;
        double lastSuspValue = 0d;

        TreeSet<Double> kms = new TreeSet<Double>(load.getKms());

        for (double km: kms) {
            logger.debug ("Trying to add at km " + km);
            SedimentLoadFraction fraction = load.getFraction(km);
            if (complete(fraction)) {
                double total = fraction.getCoarse() +
                    fraction.getFineMiddle() +
                    fraction.getSand() +
                    fraction.getSuspSand() +
                    fraction.getSuspSediment();
                // Easiest case. Add values up and set'em.
                if (fraction.getCoarseRange().equals(
                    fraction.getSuspSedimentRange())) {
                    lastOtherRange = null;
                    lastSuspRange = null;
                    fairLoad.setTotal(km, total, fraction.getCoarseRange());
                }
                else {
                    // Need to split a range.
                    if (fraction.getCoarseRange().getEnd()
                        < fraction.getSuspSedimentRange().getEnd()) {
                        // Schwebstoff is longer.
                        // Adjust and remember schwebstoffs range and value.
                        lastSuspRange = (Range) fraction.getSuspSedimentRange().clone();
                        lastSuspRange.setStart(fraction.getCoarseRange().getEnd());
                        lastSuspValue = fraction.getSuspSediment();
                        lastOtherRange = null;
                        fairLoad.setTotal(km, total, fraction.getCoarseRange());
                    }
                    else {
                        // Geschiebe is longer.
                        // Adjust and remember other values.
                        lastOtherRange = (Range) fraction.getSuspSedimentRange().clone();
                        lastOtherRange.setStart(fraction.getSuspSedimentRange().getEnd());
                        lastOtherValue = (total - fraction.getSuspSediment());
                        lastSuspRange = null;
                        fairLoad.setTotal(km, total, fraction.getSuspSedimentRange());
                    }
                }
            }
            else if (hasOnlySuspValues(fraction) && lastOtherRange != null) {
                // Split stuff.
                Range suspSedimentRange = fraction.getSuspSedimentRange();
                // if intersects with last other range, cool! merge and add!
                if (lastOtherRange.contains(km)) {
                    double maxStart = 0d;
                    double minEnd = 0d;
                    maxStart = Math.max(suspSedimentRange.getStart(),
                        lastOtherRange.getStart());

                    minEnd = Math.min(suspSedimentRange.getEnd(),
                        lastOtherRange.getEnd());
                    double total = lastOtherValue + fraction.getSuspSediment();
                    Range totalRange = new Range(maxStart, minEnd);
                    if (suspSedimentRange.getEnd() > lastOtherRange.getEnd()) {
                        lastSuspRange = (Range) suspSedimentRange.clone();
                        lastSuspRange.setStart(lastOtherRange.getEnd());
                        lastSuspValue = fraction.getSuspSediment();
                        lastOtherRange = null;
                    }
                    else {
                        // Other is "longer".
                        lastOtherRange.setStart(suspSedimentRange.getEnd());
                        lastSuspRange = null;
                    }
                    if (Math.abs(suspSedimentRange.getEnd() - lastOtherRange.getEnd()) < 0.1d) {
                        lastOtherRange = null;
                        lastSuspRange = null;
                    }
                    fairLoad.setTotal(km, total + fraction.getSuspSediment(), totalRange);
                }
                else {
                    lastSuspRange = suspSedimentRange;
                    lastSuspValue = fraction.getSuspSediment();
                    lastOtherRange = null;
                }
            }
            else if (hasButSuspValues(fraction) && lastSuspRange != null) {
                // If intersects with last suspsed range, merge and add
                double total = fraction.getCoarse() +
                    fraction.getFineMiddle() +
                    fraction.getSand() +
                    fraction.getSuspSand() +
                    lastSuspValue;
                double maxStart = Math.max(fraction.getCoarseRange().getStart(),
                    lastSuspRange.getStart());
                if (lastSuspRange.contains(km)) {
                    double minEnd = Math.min(fraction.getCoarseRange().getEnd(),
                        lastSuspRange.getEnd());
                    Range totalRange = new Range(maxStart, minEnd);
                    if (lastSuspRange.getEnd() > fraction.getCoarseRange().getEnd()) {
                        // SuspSed longer.
                        lastSuspRange.setStart(fraction.getCoarseRange().getEnd());
                        lastOtherRange = null;
                    }
                    else {
                        // Other longer
                        lastOtherRange = (Range) fraction.getCoarseRange().clone();
                        lastOtherRange.setStart(lastSuspRange.getEnd());
                        lastSuspRange = null;
                        lastOtherValue = total - lastSuspValue;
                    }
                    if (lastSuspRange != null && Math.abs(lastSuspRange.getEnd() - lastOtherRange.getEnd()) < 0.1d) {
                        lastOtherRange = null;
                        lastSuspRange = null;
                    }
                    fairLoad.setTotal(km, total, totalRange);
                }
                else {
                    // Ranges are disjoint.
                    lastOtherRange = fraction.getCoarseRange();
                    lastOtherValue = total - fraction.getSuspSediment();
                    lastSuspRange = null;
                }
            }
            else {
                // Some values are missing or no intersection with former values.
                // Stay as we are.
                if (hasButSuspValues(fraction)) {
                    double total = fraction.getCoarse() +
                        fraction.getFineMiddle() +
                        fraction.getSand() +
                        fraction.getSuspSand();
                    lastOtherRange = fraction.getCoarseRange();
                    lastOtherValue = total;
                    lastSuspRange = null;
                }
                else if (hasOnlySuspValues(fraction)) {
                    lastSuspRange = fraction.getSuspSedimentRange();
                    lastSuspValue = fraction.getSuspSediment();
                    lastOtherRange = null;
                }
            }
        }
        return fairLoad;
    }


    private SedimentLoad calculateUnit(SedimentLoad load, int year) {
        SedimentDensity density =
            SedimentDensityFactory.getSedimentDensity(river, kmLow, kmUp, year);
        for (double km: load.getKms()) {
            double dens = density.getDensity(km, year);
            SedimentLoadFraction fraction = load.getFraction(km);
            double coarse = fraction.getCoarse();
            double fineMiddle = fraction.getFineMiddle();
            double sand = fraction.getSand();
            double suspSand = fraction.getSuspSand();
            double bedSand = fraction.getSuspSandBed();
            double sediment = fraction.getSuspSediment();
            double total = fraction.getTotal();
            load.setCoarse(km, (coarse * dens), fraction.getCoarseRange());
            load.setFineMiddle(km, (fineMiddle * dens), fraction.getFineMiddleRange());
            load.setSand(km, (sand * dens), fraction.getSandRange());
            load.setSuspSand(km, (suspSand * dens), fraction.getSuspSandRange());
            load.setSuspSandBed(km, (bedSand * dens), fraction.getSuspSandBedRange());
            load.setSuspSediment(km, (sediment * dens), fraction.getSuspSedimentRange());
            load.setTotal(km, (total * dens), fraction.getTotalRange());
        }
        return load;
    }
}
// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf-8 :

http://dive4elements.wald.intevation.org