view artifacts/src/main/java/org/dive4elements/river/artifacts/model/minfo/SedimentLoadDataCalculation.java @ 8172:f2bbe09e516e

Fix fraction respectively facet names and descriptions.
author "Tom Gottfried <tom@intevation.de>"
date Mon, 01 Sep 2014 12:34:39 +0200
parents 363b82ecf29f
children d2673ca68e70
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 int [] AVERAGE_FLYS = {
        SedimentLoadData.GF_AVERAGE
    };

    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 for bfg parts
        // Grain fraction names are alignt to the grain_fractions table
        GrainFraction.make("total",              TOTAL_LOAD_FLYS),
        GrainFraction.make("bed_load",           BED_LOAD_FLYS),
        GrainFraction.make("bed_load_susp_sand", BED_LOAD_SUSP_SAND_FLYS),
        GrainFraction.make("coarse",             COARSE_FLYS),
        GrainFraction.make("fine_middle",        FINE_MIDDLE_FLYS),
        GrainFraction.make("sand",               SAND_FLYS) ,
        GrainFraction.make("susp_sand",          SUSP_SAND_FLYS),
        GrainFraction.make("susp_sand_bed",      SUSP_SAND_BED_FLYS),
        GrainFraction.make("suspended_sediment", SUSP_SEDIMENT_FLYS),
        GrainFraction.make("average",            AVERAGE_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();

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

        for (int year = min; year <= max; ++year) {
            Value.Filter filter = new And(notEpochs)
                .add(new TimeRangeIntersects(year));
            String period = Integer.toString(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, period);

                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]);

            String period = Integer.toString(epoch[0]) + " - " +
                Integer.toString(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, period);
            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]);

            String period = Integer.toString(epoch[0]) + " - " +
                Integer.toString(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, period);
                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