teichmann@8049: /* Copyright (C) 2014 by Bundesanstalt für Gewässerkunde teichmann@8049: * Software engineering by Intevation GmbH teichmann@8049: * teichmann@8049: * This file is Free Software under the GNU AGPL (>=v3) teichmann@8049: * and comes with ABSOLUTELY NO WARRANTY! Check out the teichmann@8049: * documentation coming with Dive4Elements River for details. teichmann@8049: */ teichmann@8049: package org.dive4elements.river.artifacts.model.minfo; teichmann@8049: teichmann@8067: import java.util.ArrayList; teichmann@8049: import java.util.List; teichmann@8067: import java.util.Map; teichmann@8049: import java.util.Set; teichmann@8067: import java.util.TreeMap; teichmann@8055: import java.util.TreeSet; teichmann@8049: teichmann@8049: import org.dive4elements.river.artifacts.access.SedimentLoadAccess; teichmann@8049: import org.dive4elements.river.artifacts.model.Calculation; teichmann@8049: import org.dive4elements.river.artifacts.model.CalculationResult; teichmann@8053: import org.dive4elements.river.artifacts.model.RiverFactory; teichmann@8049: import org.apache.log4j.Logger; teichmann@8049: import org.dive4elements.river.artifacts.model.minfo.SedimentLoadData.Value; teichmann@8049: import org.dive4elements.river.artifacts.model.minfo.SedimentLoadData.Station; teichmann@8055: import org.dive4elements.river.artifacts.model.minfo.SedimentLoadDataValueFilter.And; teichmann@8055: import org.dive4elements.river.artifacts.model.minfo.SedimentLoadDataValueFilter.IsEpoch; teichmann@8060: import org.dive4elements.river.artifacts.model.minfo.SedimentLoadDataValueFilter.IsOfficial; teichmann@8055: import org.dive4elements.river.artifacts.model.minfo.SedimentLoadDataValueFilter.Not; teichmann@8055: import org.dive4elements.river.artifacts.model.minfo.SedimentLoadDataValueFilter.TimeRangeIntersects; teichmann@8053: import org.dive4elements.river.model.River; teichmann@8053: import org.dive4elements.river.utils.DoubleUtil; teichmann@8049: teichmann@8049: public class SedimentLoadDataCalculation teichmann@8049: extends Calculation teichmann@8049: { teichmann@8049: private static final Logger log = Logger teichmann@8049: .getLogger(SedimentLoadDataCalculation.class); teichmann@8049: teichmann@8049: public static final int [] TOTAL_LOAD_FLYS = { teichmann@8049: SedimentLoadData.GF_COARSE, teichmann@8049: SedimentLoadData.GF_FINE_MIDDLE, teichmann@8049: SedimentLoadData.GF_SAND, teichmann@8049: SedimentLoadData.GF_SUSP_SEDIMENT teichmann@8049: }; teichmann@8049: teichmann@8049: public static final int [] BED_LOAD_FLYS = { teichmann@8049: SedimentLoadData.GF_COARSE, teichmann@8049: SedimentLoadData.GF_FINE_MIDDLE, teichmann@8049: SedimentLoadData.GF_SAND teichmann@8049: }; teichmann@8049: teichmann@8049: public static final int [] BED_LOAD_SUSP_SAND_FLYS = { teichmann@8049: SedimentLoadData.GF_COARSE, teichmann@8049: SedimentLoadData.GF_FINE_MIDDLE, teichmann@8049: SedimentLoadData.GF_SAND, teichmann@8049: SedimentLoadData.GF_SUSP_SAND teichmann@8049: }; teichmann@8049: teichmann@8050: public static final int [] TOTAL_LOAD_BFG = { teichmann@8050: SedimentLoadData.GF_TOTAL teichmann@8050: }; teichmann@8050: teichmann@8050: public static final int [] BED_LOAD_BFG = { teichmann@8050: SedimentLoadData.GF_BED_LOAD teichmann@8050: }; teichmann@8050: teichmann@8050: public static final int [] SUSPENDED_LOAD_BFG = { teichmann@8050: SedimentLoadData.GF_SUSPENDED_LOAD teichmann@8050: }; teichmann@8050: teichmann@8050: public static final int [] COARSE_FLYS = { teichmann@8050: SedimentLoadData.GF_COARSE teichmann@8050: }; teichmann@8050: teichmann@8050: public static final int [] FINE_MIDDLE_FLYS = { teichmann@8050: SedimentLoadData.GF_FINE_MIDDLE teichmann@8050: }; teichmann@8050: teichmann@8050: public static final int [] SAND_FLYS = { teichmann@8050: SedimentLoadData.GF_SAND teichmann@8050: }; teichmann@8050: teichmann@8050: public static final int [] SUSP_SAND_FLYS = { teichmann@8050: SedimentLoadData.GF_SUSP_SAND teichmann@8050: }; teichmann@8050: teichmann@8050: public static final int [] SUSP_SAND_BED_FLYS = { teichmann@8050: SedimentLoadData.GF_SUSP_SAND_BED teichmann@8050: }; teichmann@8050: teichmann@8050: public static final int [] SUSP_SEDIMENT_FLYS = { teichmann@8050: SedimentLoadData.GF_SUSP_SEDIMENT teichmann@8050: }; teichmann@8050: andre@8130: public static final int [] AVERAGE_FLYS = { andre@8130: SedimentLoadData.GF_AVERAGE andre@8130: }; andre@8130: teichmann@8055: public static final class GrainFraction { teichmann@8055: private String description; teichmann@8055: private int [] grainFractions; teichmann@8055: teichmann@8055: public GrainFraction(String description, int [] grainFractions) { teichmann@8055: this.description = description; teichmann@8055: this.grainFractions = grainFractions; teichmann@8055: } teichmann@8055: public static final GrainFraction make(String description, int [] grainFractions) { teichmann@8055: return new GrainFraction(description, grainFractions); teichmann@8055: } teichmann@8055: teichmann@8055: public String getDescription() { teichmann@8055: return description; teichmann@8055: } teichmann@8055: teichmann@8055: public int [] getGrainFractions() { teichmann@8055: return grainFractions; teichmann@8055: } teichmann@8055: } // class GrainFraction teichmann@8055: teichmann@8055: public static final GrainFraction [] GRAIN_FRACTIONS = { andre@8130: // TODO: i18n for bfg parts andre@8130: // Grain fraction names are alignt to the grain_fractions table andre@8130: GrainFraction.make("total", TOTAL_LOAD_FLYS), andre@8130: GrainFraction.make("bed_load", BED_LOAD_FLYS), andre@8130: GrainFraction.make("susp_sand_bed", BED_LOAD_SUSP_SAND_FLYS), andre@8130: GrainFraction.make("total.bfg", TOTAL_LOAD_BFG), andre@8130: GrainFraction.make("bed_load.bfg", BED_LOAD_BFG), andre@8130: GrainFraction.make("suspended_load.bfg", SUSPENDED_LOAD_BFG), andre@8130: GrainFraction.make("coarse", COARSE_FLYS), andre@8130: GrainFraction.make("fine_middle", FINE_MIDDLE_FLYS), andre@8130: GrainFraction.make("sand", SAND_FLYS) , andre@8130: GrainFraction.make("susp_sand", SUSP_SAND_FLYS), andre@8130: GrainFraction.make("susp_sand_bed", SUSP_SAND_BED_FLYS), andre@8130: GrainFraction.make("suspended_sediment", SUSP_SEDIMENT_FLYS), andre@8130: GrainFraction.make("average", AVERAGE_FLYS), teichmann@8055: }; teichmann@8055: teichmann@8060: public static class Sum implements Value.Visitor { teichmann@8049: teichmann@8060: protected int n; teichmann@8060: protected double sum; teichmann@8049: teichmann@8049: public Sum() { teichmann@8049: } teichmann@8049: teichmann@8049: public double getSum() { teichmann@8060: return sum; teichmann@8049: } teichmann@8049: teichmann@8049: public int getN() { teichmann@8049: return n; teichmann@8049: } teichmann@8049: teichmann@8049: public void reset() { teichmann@8049: n = 0; teichmann@8049: sum = 0.0; teichmann@8049: } teichmann@8049: teichmann@8049: @Override teichmann@8049: public void visit(Value value) { teichmann@8049: sum += value.getValue(); teichmann@8049: ++n; teichmann@8049: } teichmann@8060: } // class Sum teichmann@8060: teichmann@8052: private String river; teichmann@8052: private String yearEpoch; teichmann@8052: private String unit; teichmann@8055: private int [][] epochs; teichmann@8055: private int [] years; teichmann@8052: private double from; teichmann@8052: private double to; teichmann@8052: teichmann@8052: teichmann@8049: public SedimentLoadDataCalculation() { teichmann@8049: } teichmann@8049: teichmann@8049: public CalculationResult calculate(SedimentLoadAccess access) { teichmann@8049: log.info("SedimentLoadDataCalculation.calculate"); teichmann@8052: teichmann@8052: String river = access.getRiverName(); teichmann@8052: String yearEpoch = access.getYearEpoch(); teichmann@8052: String unit = access.getUnit(); teichmann@8052: teichmann@8055: int [] years = null; teichmann@8055: int [][] epochs = null; teichmann@8052: teichmann@8057: double from = access.getLowerKM(); teichmann@8057: double to = access.getUpperKM(); teichmann@8052: teichmann@8052: if (yearEpoch.equals("year")) { teichmann@8055: years = access.getPeriod(); teichmann@8052: } teichmann@8052: else if (yearEpoch.equals("epoch") || yearEpoch.equals("off_epoch")) { teichmann@8055: epochs = access.getEpochs(); teichmann@8052: } teichmann@8052: else { teichmann@8052: addProblem("minfo.missing.year_epoch"); teichmann@8052: } teichmann@8052: teichmann@8052: if (river == null) { teichmann@8052: // TODO: i18n teichmann@8052: addProblem("minfo.missing.river"); teichmann@8052: } teichmann@8052: teichmann@8055: if (years == null && epochs == null) { teichmann@8052: addProblem("minfo.missing.time"); teichmann@8052: } teichmann@8052: teichmann@8052: if (!hasProblems()) { teichmann@8052: this.river = river; teichmann@8052: this.yearEpoch = yearEpoch; teichmann@8052: this.unit = unit; teichmann@8055: this.years = years; teichmann@8055: this.epochs = epochs; teichmann@8052: this.from = from; teichmann@8052: this.to = to; teichmann@8052: return internalCalculate(); teichmann@8052: } teichmann@8052: teichmann@8055: return error(null); teichmann@8055: } teichmann@8055: teichmann@8055: private CalculationResult error(String msg) { teichmann@8055: if (msg != null) addProblem(msg); teichmann@8053: return new CalculationResult(this); teichmann@8052: } teichmann@8052: teichmann@8052: private CalculationResult internalCalculate() { teichmann@8053: if ("year".equals(yearEpoch)) return calculateYears(); teichmann@8053: if ("epoch".equals(yearEpoch)) return calculateEpochs(); teichmann@8053: if ("off_epoch".equals(yearEpoch)) return calculateOffEpochs(); teichmann@8053: teichmann@8053: // TODO: i18n teichmann@8055: return error("minfo.sediment.load.unknown.calc.mode"); teichmann@8053: } teichmann@8053: teichmann@8053: private CalculationResult calculateYears() { teichmann@8055: SedimentLoadData sld = teichmann@8055: SedimentLoadDataFactory.INSTANCE.getSedimentLoadData(river); teichmann@8055: if (sld == null) { teichmann@8055: // TODO: i18n teichmann@8055: return error("minfo.sediment.load.no.data"); teichmann@8055: } teichmann@8055: teichmann@8062: SedimentLoadDataResult sldr = new SedimentLoadDataResult(); teichmann@8062: teichmann@8057: boolean isKmUp = isKmUp(); teichmann@8055: Set missingFractions = new TreeSet(); teichmann@8055: teichmann@8055: Not notEpochs = new Not(IsEpoch.INSTANCE); teichmann@8055: teichmann@8055: Sum sum = new Sum(); teichmann@8055: teichmann@8066: SedimentDensity sd = getSedimentDensity(); teichmann@8066: tom@8145: int min = Math.min(years[0], years[1]); tom@8145: int max = Math.max(years[0], years[1]); tom@8145: tom@8145: for (int year = min; year <= max; ++year) { teichmann@8068: Value.Filter filter = new And(notEpochs) teichmann@8055: .add(new TimeRangeIntersects(year)); andre@8131: String period = Integer.toString(year); teichmann@8055: teichmann@8055: for (GrainFraction gf: GRAIN_FRACTIONS) { teichmann@8055: double [][] result = sum( teichmann@8055: sld, gf.getGrainFractions(), filter, sum, isKmUp, teichmann@8055: missingFractions); teichmann@8055: teichmann@8055: if (result[0].length == 0 || DoubleUtil.isNaN(result[1])) { teichmann@8055: // TODO: resolve i18n teichmann@8055: addProblem("minfo.sediment.load.no.fractions", teichmann@8055: gf.getDescription()); teichmann@8055: continue; teichmann@8055: } teichmann@8066: teichmann@8066: transformT2M3(sd, year, result); teichmann@8066: teichmann@8062: SedimentLoadDataResult.Fraction sldrf = andre@8131: new SedimentLoadDataResult.Fraction(gf.getDescription(), result, period); teichmann@8066: teichmann@8062: sldr.addFraction(sldrf); teichmann@8055: } tom@8145: tom@8145: // Do not give single year twice tom@8145: if (min == max) break; teichmann@8055: } teichmann@8055: // TODO: Generate messages for missing fractions. teichmann@8062: return new CalculationResult(sldr, this); teichmann@8049: } teichmann@8049: teichmann@8066: teichmann@8053: private CalculationResult calculateEpochs() { teichmann@8057: SedimentLoadData sld = teichmann@8057: SedimentLoadDataFactory.INSTANCE.getSedimentLoadData(river); teichmann@8057: if (sld == null) { teichmann@8057: // TODO: i18n teichmann@8057: return error("minfo.sediment.load.no.data"); teichmann@8057: } teichmann@8057: teichmann@8062: SedimentLoadDataResult sldr = new SedimentLoadDataResult(); teichmann@8062: teichmann@8057: boolean isKmUp = isKmUp(); teichmann@8057: Set missingFractions = new TreeSet(); teichmann@8057: teichmann@8067: SedimentDensity sd = getSedimentDensity(); teichmann@8067: teichmann@8060: // They are not epochs, they are single years! teichmann@8060: Not notEpochs = new Not(IsEpoch.INSTANCE); teichmann@8060: teichmann@8057: for (int [] epoch: epochs) { teichmann@8067: List results = new ArrayList(); teichmann@8060: teichmann@8067: int min = Math.min(epoch[0], epoch[1]); teichmann@8067: int max = Math.max(epoch[0], epoch[1]); teichmann@8067: andre@8131: String period = Integer.toString(epoch[0]) + " - " + andre@8131: Integer.toString(epoch[1]); teichmann@8067: for (int year = min; year <= max; ++year) { teichmann@8068: Value.Filter filter = new And(notEpochs) teichmann@8067: .add(new TimeRangeIntersects(year)); teichmann@8067: teichmann@8067: Sum sum = new Sum(); teichmann@8067: teichmann@8067: for (GrainFraction gf: GRAIN_FRACTIONS) { teichmann@8067: double [][] result = sum( teichmann@8067: sld, gf.getGrainFractions(), filter, sum, isKmUp, teichmann@8067: missingFractions); teichmann@8067: teichmann@8067: if (result[0].length == 0 || DoubleUtil.isNaN(result[1])) { teichmann@8067: // TODO: resolve i18n teichmann@8067: addProblem("minfo.sediment.load.no.fractions", teichmann@8067: gf.getDescription()); teichmann@8067: continue; teichmann@8067: } teichmann@8067: teichmann@8067: transformT2M3(sd, year, result); teichmann@8067: results.add(result); teichmann@8060: } teichmann@8060: } teichmann@8067: teichmann@8067: double [][] result = average(results); teichmann@8067: // TODO: Optionally transform units. teichmann@8067: SedimentLoadDataResult.Fraction sldrf = andre@8131: new SedimentLoadDataResult.Fraction("TODO: nice description", result, period); teichmann@8067: sldr.addFraction(sldrf); teichmann@8060: } teichmann@8060: // TODO: Generate messages for missing fractions. teichmann@8062: return new CalculationResult(sldr, this); teichmann@8060: } teichmann@8060: teichmann@8060: private CalculationResult calculateOffEpochs() { teichmann@8060: SedimentLoadData sld = teichmann@8060: SedimentLoadDataFactory.INSTANCE.getSedimentLoadData(river); teichmann@8060: if (sld == null) { teichmann@8060: // TODO: i18n teichmann@8060: return error("minfo.sediment.load.no.data"); teichmann@8060: } teichmann@8060: teichmann@8062: SedimentLoadDataResult sldr = new SedimentLoadDataResult(); teichmann@8062: teichmann@8068: SedimentDensity sd = getSedimentDensity(); teichmann@8068: teichmann@8060: boolean isKmUp = isKmUp(); teichmann@8060: Set missingFractions = new TreeSet(); teichmann@8060: teichmann@8060: for (int [] epoch: epochs) { teichmann@8068: Value.Filter filter = new And(IsOfficial.INSTANCE) teichmann@8060: .add(new TimeRangeIntersects(epoch[0], epoch[1])); teichmann@8060: teichmann@8068: int year = Math.min(epoch[0], epoch[1]); teichmann@8068: andre@8131: String period = Integer.toString(epoch[0]) + " - " + andre@8131: Integer.toString(epoch[1]); andre@8131: teichmann@8060: Sum sum = new Sum(); teichmann@8057: teichmann@8057: for (GrainFraction gf: GRAIN_FRACTIONS) { teichmann@8057: double [][] result = sum( teichmann@8057: sld, gf.getGrainFractions(), filter, sum, isKmUp, teichmann@8057: missingFractions); teichmann@8057: teichmann@8057: if (result[0].length == 0 || DoubleUtil.isNaN(result[1])) { teichmann@8057: // TODO: resolve i18n teichmann@8057: addProblem("minfo.sediment.load.no.fractions", teichmann@8057: gf.getDescription()); teichmann@8057: continue; teichmann@8057: } teichmann@8068: transformT2M3(sd, year, result); teichmann@8062: SedimentLoadDataResult.Fraction sldrf = andre@8131: new SedimentLoadDataResult.Fraction(gf.getDescription(), result, period); teichmann@8062: sldr.addFraction(sldrf); teichmann@8057: } teichmann@8057: } teichmann@8057: // TODO: Generate messages for missing fractions. teichmann@8062: return new CalculationResult(sldr, this); teichmann@8053: } teichmann@8053: teichmann@8053: /** Figure out flow direction of river. */ teichmann@8068: private final boolean isKmUp() { teichmann@8053: River r = RiverFactory.getRiver(river); teichmann@8053: if (r == null) { teichmann@8053: addProblem("minfo.missing.river"); teichmann@8053: return true; teichmann@8049: } teichmann@8053: return r.getKmUp(); teichmann@8049: } teichmann@8049: teichmann@8066: private final boolean inM3() { teichmann@8066: return unit.equals("m3_per_a"); teichmann@8066: } teichmann@8066: teichmann@8066: private SedimentDensity getSedimentDensity() { teichmann@8066: return inM3() teichmann@8066: ? SedimentDensityFactory.getSedimentDensity(river, from, to) teichmann@8066: : null; teichmann@8066: } teichmann@8066: teichmann@8066: private static void transformT2M3(SedimentDensity sd, int year, double [][] data) { teichmann@8066: if (sd == null) { teichmann@8066: return; teichmann@8066: } teichmann@8066: double [] kms = data[0]; teichmann@8066: double [] values = data[1]; teichmann@8066: for (int i = 0; i < kms.length; ++i) { teichmann@8066: if (Double.isNaN(kms[i]) || Double.isNaN(kms[i])) { teichmann@8066: continue; teichmann@8066: } teichmann@8066: double density = sd.getDensity(kms[i], year); teichmann@8066: values[i] /= density; teichmann@8066: } teichmann@8066: } teichmann@8066: teichmann@8049: public double[][] sum( teichmann@8049: SedimentLoadData sld, teichmann@8049: int [] grainFractions, teichmann@8049: Value.Filter filter, teichmann@8049: Sum sum, teichmann@8049: boolean isKMUp, teichmann@8049: Set missingFractions teichmann@8049: ) { teichmann@8049: List stations = sld.findStations(from, to); teichmann@8049: teichmann@8049: double [] values = new double[grainFractions.length]; teichmann@8049: teichmann@8049: double [][] result = new double[2][stations.size()]; teichmann@8049: teichmann@8049: for (int j = 0, S = stations.size(); j < S; ++j) { teichmann@8049: Station station = stations.get(j); teichmann@8049: for (int i = 0; i < grainFractions.length; ++i) { teichmann@8049: int gf = grainFractions[i]; teichmann@8049: sum.reset(); teichmann@8049: station.filterGrainFraction(gf, filter, sum); teichmann@8049: if (sum.getN() == 0) { // No values found teichmann@8049: int msType = SedimentLoadData.measurementStationType(gf); teichmann@8049: // Station of right fraction type already? No: take previous. teichmann@8049: if (!station.isType(msType)) { teichmann@8049: Station prev = station.prevByType(msType, isKMUp); teichmann@8049: if (prev != null) { teichmann@8049: prev.filterGrainFraction(gf, filter, sum); teichmann@8049: } teichmann@8049: } teichmann@8049: } teichmann@8049: teichmann@8049: if (sum.getN() == 0) { teichmann@8049: missingFractions.add(gf); teichmann@8049: values[i] = Double.NaN; teichmann@8049: } else { teichmann@8049: values[i] = sum.getSum(); teichmann@8049: } teichmann@8049: } teichmann@8049: result[0][j] = station.getStation(); teichmann@8053: result[1][j] = DoubleUtil.sum(values); teichmann@8049: } teichmann@8049: teichmann@8049: // TODO: Handle 'virtual' measument stations 'from' and 'to'. teichmann@8049: teichmann@8049: return result; teichmann@8049: } teichmann@8067: teichmann@8067: private static final class XSum { teichmann@8067: private double sum; teichmann@8067: private int n; teichmann@8067: public XSum() { teichmann@8067: } teichmann@8067: public void add(double v) { teichmann@8067: sum += v; teichmann@8067: ++n; teichmann@8067: } teichmann@8067: public double avg() { teichmann@8067: return sum/n; teichmann@8067: } teichmann@8067: } teichmann@8067: teichmann@8067: private static double [][] average(List data) { teichmann@8067: teichmann@8067: TreeMap map = new TreeMap(); teichmann@8067: teichmann@8067: for (double [][] pair: data) { teichmann@8067: double [] kms = pair[0]; teichmann@8067: double [] vs = pair[1]; teichmann@8067: for (int i = 0; i < kms.length; ++i) { teichmann@8067: double km = kms[i]; teichmann@8067: double v = vs[i]; teichmann@8067: if (Double.isNaN(km) || Double.isNaN(v)) { teichmann@8067: continue; teichmann@8067: } teichmann@8067: XSum xsum = map.get(km); teichmann@8067: if (xsum == null) { teichmann@8067: map.put(km, xsum = new XSum()); teichmann@8067: } teichmann@8067: xsum.add(v); teichmann@8067: } teichmann@8067: } teichmann@8067: teichmann@8067: double [][] result = new double[2][map.size()]; teichmann@8067: int i = 0; teichmann@8067: for (Map.Entry entry: map.entrySet()) { teichmann@8067: result[i][0] = entry.getKey(); teichmann@8067: result[i][1] = entry.getValue().avg(); teichmann@8067: ++i; teichmann@8067: } teichmann@8067: teichmann@8067: return null; teichmann@8067: } teichmann@8067: teichmann@8049: } teichmann@8049: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf-8 :