view artifacts/src/main/java/org/dive4elements/river/artifacts/model/minfo/SedimentLoadDataCalculation.java @ 8098:09725b65955a

Add new and simplyfied SedimentLoadFacet The SedimentLoadFacet is intended to work with the Measurement stations. It uses the same mechanismn to access the Mesurement station values as the calculation does. SedimentLoadLS values need a different facet that will come soon.
author Andre Heinecke <andre.heinecke@intevation.de>
date Fri, 15 Aug 2014 18:27:19 +0200
parents 640fffb0ad60
children bbad52b073a4
line wrap: on
line source
/* Copyright (C) 2014 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 java.util.ArrayList;
import java.util.List;
import java.util.Map;
import java.util.Set;
import java.util.TreeMap;
import java.util.TreeSet;

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.RiverFactory;
import org.apache.log4j.Logger;
import org.dive4elements.river.artifacts.model.minfo.SedimentLoadData.Value;
import org.dive4elements.river.artifacts.model.minfo.SedimentLoadData.Station;
import org.dive4elements.river.artifacts.model.minfo.SedimentLoadDataValueFilter.And;
import org.dive4elements.river.artifacts.model.minfo.SedimentLoadDataValueFilter.IsEpoch;
import org.dive4elements.river.artifacts.model.minfo.SedimentLoadDataValueFilter.IsOfficial;
import org.dive4elements.river.artifacts.model.minfo.SedimentLoadDataValueFilter.Not;
import org.dive4elements.river.artifacts.model.minfo.SedimentLoadDataValueFilter.TimeRangeIntersects;
import org.dive4elements.river.model.River;
import org.dive4elements.river.utils.DoubleUtil;

public class SedimentLoadDataCalculation
extends      Calculation
{
    private static final Logger log = Logger
        .getLogger(SedimentLoadDataCalculation.class);

    public static final int [] TOTAL_LOAD_FLYS = {
        SedimentLoadData.GF_COARSE,
        SedimentLoadData.GF_FINE_MIDDLE,
        SedimentLoadData.GF_SAND,
        SedimentLoadData.GF_SUSP_SEDIMENT
    };

    public static final int [] BED_LOAD_FLYS = {
        SedimentLoadData.GF_COARSE,
        SedimentLoadData.GF_FINE_MIDDLE,
        SedimentLoadData.GF_SAND
    };

    public static final int [] BED_LOAD_SUSP_SAND_FLYS = {
        SedimentLoadData.GF_COARSE,
        SedimentLoadData.GF_FINE_MIDDLE,
        SedimentLoadData.GF_SAND,
        SedimentLoadData.GF_SUSP_SAND
    };

    public static final int [] TOTAL_LOAD_BFG = {
        SedimentLoadData.GF_TOTAL
    };

    public static final int [] BED_LOAD_BFG = {
        SedimentLoadData.GF_BED_LOAD
    };

    public static final int [] SUSPENDED_LOAD_BFG = {
        SedimentLoadData.GF_SUSPENDED_LOAD
    };

    public static final int [] COARSE_FLYS = {
        SedimentLoadData.GF_COARSE
    };

    public static final int [] FINE_MIDDLE_FLYS = {
        SedimentLoadData.GF_FINE_MIDDLE
    };

    public static final int [] SAND_FLYS = {
        SedimentLoadData.GF_SAND
    };

    public static final int [] SUSP_SAND_FLYS = {
        SedimentLoadData.GF_SUSP_SAND
    };

    public static final int [] SUSP_SAND_BED_FLYS = {
        SedimentLoadData.GF_SUSP_SAND_BED
    };

    public static final int [] SUSP_SEDIMENT_FLYS = {
        SedimentLoadData.GF_SUSP_SEDIMENT
    };

    public static final class GrainFraction {
        private String description;
        private int [] grainFractions;

        public GrainFraction(String description, int [] grainFractions) {
            this.description = description;
            this.grainFractions = grainFractions;
        }
        public static final GrainFraction make(String description, int [] grainFractions) {
            return new GrainFraction(description, grainFractions);
        }

        public String getDescription() {
            return description;
        }

        public int [] getGrainFractions() {
            return grainFractions;
        }
    } // class GrainFraction

    public static final GrainFraction [] GRAIN_FRACTIONS = {
        // TODO: i18n
        GrainFraction.make("minfo.total.load.flys",         TOTAL_LOAD_FLYS),
        GrainFraction.make("minfo.bed.load.flys",           BED_LOAD_FLYS),
        GrainFraction.make("minfo.bed.load.susp.sand.flys", BED_LOAD_SUSP_SAND_FLYS),
        GrainFraction.make("minfo.total.load.bfg",          TOTAL_LOAD_BFG),
        GrainFraction.make("minfo.bed.load.bfg",            BED_LOAD_BFG),
        GrainFraction.make("minfo.suspended.load.bfg",      SUSPENDED_LOAD_BFG),
        GrainFraction.make("minfo.coarse.flys",             COARSE_FLYS),
        GrainFraction.make("minfo.fine.middle.flys",        FINE_MIDDLE_FLYS),
        GrainFraction.make("minfo.sand.flys",               SAND_FLYS) ,
        GrainFraction.make("minfo.susp.sand.flys",          SUSP_SAND_FLYS),
        GrainFraction.make("minfo.susp.sand.bed.flys",      SUSP_SAND_BED_FLYS),
        GrainFraction.make("minfo.susp.sediment.flys",      SUSP_SEDIMENT_FLYS),
    };

    public static class Sum implements Value.Visitor {

        protected int    n;
        protected double sum;

        public Sum() {
        }

        public double getSum() {
            return sum;
        }

        public int getN() {
            return n;
        }

        public void reset() {
            n   = 0;
            sum = 0.0;
        }

        @Override
        public void visit(Value value) {
            sum += value.getValue();
            ++n;
        }
    } // class Sum

    private String   river;
    private String   yearEpoch;
    private String   unit;
    private int [][] epochs;
    private int []   years;
    private double   from;
    private double   to;


    public SedimentLoadDataCalculation() {
    }

    public CalculationResult calculate(SedimentLoadAccess access) {
        log.info("SedimentLoadDataCalculation.calculate");

        String river     = access.getRiverName();
        String yearEpoch = access.getYearEpoch();
        String unit      = access.getUnit();

        int [] years  = null;
        int [][] epochs = null;

        double from = access.getLowerKM();
        double to   = access.getUpperKM();

        if (yearEpoch.equals("year")) {
            years = access.getPeriod();
        }
        else if (yearEpoch.equals("epoch") || yearEpoch.equals("off_epoch")) {
            epochs = access.getEpochs();
        }
        else {
            addProblem("minfo.missing.year_epoch");
        }

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

        if (years == null && epochs == null) {
            addProblem("minfo.missing.time");
        }

        if (!hasProblems()) {
            this.river     = river;
            this.yearEpoch = yearEpoch;
            this.unit      = unit;
            this.years     = years;
            this.epochs    = epochs;
            this.from      = from;
            this.to        = to;
            return internalCalculate();
        }

        return error(null);
    }

    private CalculationResult error(String msg) {
        if (msg != null) addProblem(msg);
        return new CalculationResult(this);
    }

    private CalculationResult internalCalculate() {
        if ("year".equals(yearEpoch))      return calculateYears();
        if ("epoch".equals(yearEpoch))     return calculateEpochs();
        if ("off_epoch".equals(yearEpoch)) return calculateOffEpochs();

        // TODO: i18n
        return error("minfo.sediment.load.unknown.calc.mode");
    }

    private CalculationResult calculateYears() {
        SedimentLoadData sld =
            SedimentLoadDataFactory.INSTANCE.getSedimentLoadData(river);
        if (sld == null) {
            // TODO: i18n
            return error("minfo.sediment.load.no.data");
        }

        SedimentLoadDataResult sldr = new SedimentLoadDataResult();

        boolean isKmUp = isKmUp();
        Set<Integer> missingFractions = new TreeSet<Integer>();

        Not notEpochs = new Not(IsEpoch.INSTANCE);

        Sum sum = new Sum();

        SedimentDensity sd = getSedimentDensity();

        for (int year: years) {
            Value.Filter filter = new And(notEpochs)
                .add(new TimeRangeIntersects(year));

            for (GrainFraction gf: GRAIN_FRACTIONS) {
                double [][] result = sum(
                    sld, gf.getGrainFractions(), filter, sum, isKmUp,
                    missingFractions);

                if (result[0].length == 0 || DoubleUtil.isNaN(result[1])) {
                    // TODO: resolve i18n
                    addProblem("minfo.sediment.load.no.fractions",
                        gf.getDescription());
                    continue;
                }

                transformT2M3(sd, year, result);

                SedimentLoadDataResult.Fraction sldrf =
                    new SedimentLoadDataResult.Fraction(gf.getDescription(), result);

                sldr.addFraction(sldrf);
            }
        }
        // TODO: Generate messages for missing fractions.
        return new CalculationResult(sldr, this);
    }


    private CalculationResult calculateEpochs() {
        SedimentLoadData sld =
            SedimentLoadDataFactory.INSTANCE.getSedimentLoadData(river);
        if (sld == null) {
            // TODO: i18n
            return error("minfo.sediment.load.no.data");
        }

        SedimentLoadDataResult sldr = new SedimentLoadDataResult();

        boolean isKmUp = isKmUp();
        Set<Integer> missingFractions = new TreeSet<Integer>();

        SedimentDensity sd = getSedimentDensity();

        // They are not epochs, they are single years!
        Not notEpochs = new Not(IsEpoch.INSTANCE);

        for (int [] epoch: epochs) {
            List<double [][]> results = new ArrayList<double [][]>();

            int min = Math.min(epoch[0], epoch[1]);
            int max = Math.max(epoch[0], epoch[1]);

            for (int year = min; year <= max; ++year) {
                Value.Filter filter = new And(notEpochs)
                    .add(new TimeRangeIntersects(year));

                Sum sum = new Sum();

                for (GrainFraction gf: GRAIN_FRACTIONS) {
                    double [][] result = sum(
                        sld, gf.getGrainFractions(), filter, sum, isKmUp,
                        missingFractions);

                    if (result[0].length == 0 || DoubleUtil.isNaN(result[1])) {
                        // TODO: resolve i18n
                        addProblem("minfo.sediment.load.no.fractions",
                            gf.getDescription());
                        continue;
                    }

                    transformT2M3(sd, year, result);
                    results.add(result);
                }
            }

            double [][] result = average(results);
            // TODO: Optionally transform units.
            SedimentLoadDataResult.Fraction sldrf =
                new SedimentLoadDataResult.Fraction("TODO: nice description", result);
            sldr.addFraction(sldrf);
        }
        // TODO: Generate messages for missing fractions.
        return new CalculationResult(sldr, this);
    }

    private CalculationResult calculateOffEpochs() {
        SedimentLoadData sld =
            SedimentLoadDataFactory.INSTANCE.getSedimentLoadData(river);
        if (sld == null) {
            // TODO: i18n
            return error("minfo.sediment.load.no.data");
        }

        SedimentLoadDataResult sldr = new SedimentLoadDataResult();

        SedimentDensity sd = getSedimentDensity();

        boolean isKmUp = isKmUp();
        Set<Integer> missingFractions = new TreeSet<Integer>();

        for (int [] epoch: epochs) {
            Value.Filter filter = new And(IsOfficial.INSTANCE)
                .add(new TimeRangeIntersects(epoch[0], epoch[1]));

            int year = Math.min(epoch[0], epoch[1]);

            Sum sum = new Sum();

            for (GrainFraction gf: GRAIN_FRACTIONS) {
                double [][] result = sum(
                    sld, gf.getGrainFractions(), filter, sum, isKmUp,
                    missingFractions);

                if (result[0].length == 0 || DoubleUtil.isNaN(result[1])) {
                    // TODO: resolve i18n
                    addProblem("minfo.sediment.load.no.fractions",
                        gf.getDescription());
                    continue;
                }
                transformT2M3(sd, year, result);
                SedimentLoadDataResult.Fraction sldrf =
                    new SedimentLoadDataResult.Fraction(gf.getDescription(), result);
                sldr.addFraction(sldrf);
            }
        }
        // TODO: Generate messages for missing fractions.
        return new CalculationResult(sldr, this);
    }

    /** Figure out flow direction of river. */
    private final boolean isKmUp() {
        River r = RiverFactory.getRiver(river);
        if (r == null) {
            addProblem("minfo.missing.river");
            return true;
        }
        return r.getKmUp();
    }

    private final boolean inM3() {
        return unit.equals("m3_per_a");
    }

    private SedimentDensity getSedimentDensity() {
        return inM3()
            ? SedimentDensityFactory.getSedimentDensity(river, from, to)
            : null;
    }

    private static void transformT2M3(SedimentDensity sd, int year, double [][] data) {
        if (sd == null) {
            return;
        }
        double [] kms = data[0];
        double [] values = data[1];
        for (int i = 0; i < kms.length; ++i) {
            if (Double.isNaN(kms[i]) || Double.isNaN(kms[i])) {
                continue;
            }
            double density = sd.getDensity(kms[i], year);
            values[i] /= density;
        }
    }

    public double[][] sum(
        SedimentLoadData sld,
        int []           grainFractions,
        Value.Filter     filter,
        Sum              sum,
        boolean          isKMUp,
        Set<Integer>     missingFractions
    ) {
        List<Station> stations = sld.findStations(from, to);

        double [] values = new double[grainFractions.length];

        double [][] result = new double[2][stations.size()];

        for (int j = 0, S = stations.size(); j < S; ++j) {
            Station station = stations.get(j);
            for (int i = 0; i < grainFractions.length; ++i) {
                int gf = grainFractions[i];
                sum.reset();
                station.filterGrainFraction(gf, filter, sum);
                if (sum.getN() == 0) { // No values found
                    int msType = SedimentLoadData.measurementStationType(gf);
                    // Station of right fraction type already? No: take previous.
                    if (!station.isType(msType)) {
                        Station prev = station.prevByType(msType, isKMUp);
                        if (prev != null) {
                            prev.filterGrainFraction(gf, filter, sum);
                        }
                    }
                }

                if (sum.getN() == 0) {
                    missingFractions.add(gf);
                    values[i] = Double.NaN;
                } else {
                    values[i] = sum.getSum();
                }
            }
            result[0][j] = station.getStation();
            result[1][j] = DoubleUtil.sum(values);
        }

        // TODO: Handle 'virtual' measument stations 'from' and 'to'.

        return result;
    }

    private static final class XSum {
        private double sum;
        private int n;
        public XSum() {
        }
        public void add(double v) {
            sum += v;
            ++n;
        }
        public double avg() {
            return sum/n;
        }
    }

    private static double [][] average(List<double [][]> data) {

        TreeMap<Double, XSum> map = new TreeMap<Double, XSum>();

        for (double [][] pair: data) {
            double [] kms = pair[0];
            double [] vs = pair[1];
            for (int i = 0; i < kms.length; ++i) {
                double km = kms[i];
                double v = vs[i];
                if (Double.isNaN(km) || Double.isNaN(v)) {
                    continue;
                }
                XSum xsum = map.get(km);
                if (xsum == null) {
                    map.put(km, xsum = new XSum());
                }
                xsum.add(v);
            }
        }

        double [][] result = new double[2][map.size()];
        int i = 0;
        for (Map.Entry<Double, XSum> entry: map.entrySet()) {
            result[i][0] = entry.getKey();
            result[i][1] = entry.getValue().avg();
            ++i;
        }

        return null;
    }

}
// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf-8 :

http://dive4elements.wald.intevation.org