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@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.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: tom@8175: public static final int [] TOTAL_LOAD = { 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: tom@8175: public static final int [] BED_LOAD = { teichmann@8049: SedimentLoadData.GF_COARSE, teichmann@8049: SedimentLoadData.GF_FINE_MIDDLE, teichmann@8049: SedimentLoadData.GF_SAND teichmann@8049: }; teichmann@8049: tom@8175: public static final int [] BED_LOAD_SUSP_SAND = { 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: tom@8175: public static final int [] COARSE = { teichmann@8050: SedimentLoadData.GF_COARSE teichmann@8050: }; teichmann@8050: tom@8175: public static final int [] FINE_MIDDLE = { teichmann@8050: SedimentLoadData.GF_FINE_MIDDLE teichmann@8050: }; teichmann@8050: tom@8175: public static final int [] SAND = { teichmann@8050: SedimentLoadData.GF_SAND teichmann@8050: }; teichmann@8050: tom@8175: public static final int [] SUSP_SAND = { teichmann@8050: SedimentLoadData.GF_SUSP_SAND teichmann@8050: }; teichmann@8050: tom@8175: public static final int [] SUSP_SAND_BED = { teichmann@8050: SedimentLoadData.GF_SUSP_SAND_BED teichmann@8050: }; teichmann@8050: tom@8175: public static final int [] SUSP_SEDIMENT = { teichmann@8050: SedimentLoadData.GF_SUSP_SEDIMENT teichmann@8050: }; teichmann@8050: tom@8185: public static final class LoadSum { teichmann@8055: private String description; teichmann@8055: private int [] grainFractions; teichmann@8055: tom@8185: public LoadSum(String description, int [] grainFractions) { teichmann@8055: this.description = description; teichmann@8055: this.grainFractions = grainFractions; teichmann@8055: } tom@8185: public static final LoadSum make(String description, int [] grainFractions) { tom@8185: return new LoadSum(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: } tom@8193: tom@8193: public int getStationType() { tom@8193: return SedimentLoadData.measurementStationType( tom@8193: SedimentLoadData.grainFractionIndex(this.description)); tom@8193: } tom@8185: } // class LoadSum teichmann@8055: tom@8185: public static final LoadSum [] LOAD_SUMS = { tom@8185: // Names are alignt to the grain_fractions table tom@8185: LoadSum.make("total", TOTAL_LOAD), tom@8185: LoadSum.make("bed_load", BED_LOAD), tom@8185: LoadSum.make("bed_load_susp_sand", BED_LOAD_SUSP_SAND), tom@8185: LoadSum.make("coarse", COARSE), tom@8185: LoadSum.make("fine_middle", FINE_MIDDLE), tom@8185: LoadSum.make("sand", SAND) , tom@8185: LoadSum.make("susp_sand", SUSP_SAND), tom@8185: LoadSum.make("susp_sand_bed", SUSP_SAND_BED), tom@8185: LoadSum.make("suspended_sediment", SUSP_SEDIMENT), 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 tom@8213: return error("minfo.sedimentload.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) { tom@8213: return error("minfo.sedimentload.no.data"); teichmann@8055: } teichmann@8055: teichmann@8062: SedimentLoadDataResult sldr = new SedimentLoadDataResult(); teichmann@8062: 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: tom@8185: for (LoadSum ls: LOAD_SUMS) { tom@8185: teichmann@8055: double [][] result = sum( tom@8193: sld, ls.getGrainFractions(), ls.getStationType(), tom@8213: filter, sum); teichmann@8055: teichmann@8055: if (result[0].length == 0 || DoubleUtil.isNaN(result[1])) { tom@8213: addProblem("sedimentload.missing.fraction." + tom@8213: ls.getDescription(), period); teichmann@8055: continue; teichmann@8055: } teichmann@8066: teichmann@8066: transformT2M3(sd, year, result); teichmann@8066: teichmann@8062: SedimentLoadDataResult.Fraction sldrf = tom@8185: new SedimentLoadDataResult.Fraction(ls.getDescription(), result, period); teichmann@8066: teichmann@8062: sldr.addFraction(sldrf); teichmann@8055: } teichmann@8055: } teichmann@8062: return new CalculationResult(sldr, this); teichmann@8049: } teichmann@8049: teichmann@8053: private CalculationResult calculateEpochs() { teichmann@8057: SedimentLoadData sld = teichmann@8057: SedimentLoadDataFactory.INSTANCE.getSedimentLoadData(river); teichmann@8057: if (sld == null) { tom@8213: return error("minfo.sedimentload.no.data"); teichmann@8057: } teichmann@8057: teichmann@8062: SedimentLoadDataResult sldr = new SedimentLoadDataResult(); teichmann@8062: tom@8176: Sum sum = new Sum(); tom@8176: 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: 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: tom@8185: for (LoadSum ls: LOAD_SUMS) { teichmann@8067: tom@8176: List results = new ArrayList(); tom@8176: tom@8176: for (int year = min; year <= max; ++year) { tom@8176: Value.Filter filter = new And(notEpochs) tom@8176: .add(new TimeRangeIntersects(year)); tom@8176: teichmann@8067: double [][] result = sum( tom@8193: sld, ls.getGrainFractions(), ls.getStationType(), tom@8213: filter, sum); teichmann@8067: teichmann@8067: if (result[0].length == 0 || DoubleUtil.isNaN(result[1])) { teichmann@8067: continue; teichmann@8067: } teichmann@8067: teichmann@8067: transformT2M3(sd, year, result); teichmann@8067: results.add(result); teichmann@8060: } tom@8176: tom@8205: if (results.size() == 0) { tom@8213: addProblem("sedimentload.missing.fraction." + tom@8213: ls.getDescription(), period); tom@8205: continue; tom@8205: } tom@8205: tom@8176: double [][] result = average(results); tom@8176: tom@8205: if (!DoubleUtil.isNaN(result[1])) { tom@8205: SedimentLoadDataResult.Fraction sldrf = tom@8205: new SedimentLoadDataResult.Fraction( tom@8205: ls.getDescription(), result, period); tom@8205: sldr.addFraction(sldrf); tom@8205: } tom@8213: else { tom@8213: addProblem("sedimentload.missing.fraction." + tom@8213: ls.getDescription(), period); tom@8213: } teichmann@8060: } teichmann@8067: teichmann@8060: } 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) { tom@8213: return error("minfo.sedimentload.no.data"); teichmann@8060: } teichmann@8060: teichmann@8062: SedimentLoadDataResult sldr = new SedimentLoadDataResult(); teichmann@8062: teichmann@8068: SedimentDensity sd = getSedimentDensity(); teichmann@8068: 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: tom@8185: for (LoadSum ls: LOAD_SUMS) { teichmann@8057: double [][] result = sum( tom@8193: sld, ls.getGrainFractions(), ls.getStationType(), tom@8213: filter, sum); teichmann@8057: teichmann@8057: if (result[0].length == 0 || DoubleUtil.isNaN(result[1])) { tom@8213: addProblem("sedimentload.missing.fraction." + tom@8213: ls.getDescription(), period); teichmann@8057: continue; teichmann@8057: } tom@8213: teichmann@8068: transformT2M3(sd, year, result); teichmann@8062: SedimentLoadDataResult.Fraction sldrf = tom@8185: new SedimentLoadDataResult.Fraction(ls.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@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, tom@8193: int lsSType, teichmann@8049: Value.Filter filter, tom@8213: Sum sum 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); tom@8193: int sType = station.getType(); tom@8193: teichmann@8049: for (int i = 0; i < grainFractions.length; ++i) { teichmann@8049: int gf = grainFractions[i]; tom@8193: int gfSType = SedimentLoadData.measurementStationType(gf); tom@8193: teichmann@8049: sum.reset(); tom@8193: tom@8193: // Add current single fraction at current station teichmann@8049: station.filterGrainFraction(gf, filter, sum); tom@8193: tom@8193: if (gfSType == Station.UNKNOWN) { tom@8193: log.error("No measurement station type defined for " + tom@8193: "fraction-index" + gf); tom@8193: } tom@8193: else { tom@8193: if (lsSType != Station.BED_LOAD && tom@8193: lsSType != Station.SUSPENDED && tom@8193: (sType == Station.BED_LOAD || tom@8193: sType == Station.SUSPENDED) && tom@8193: sum.getN() == 0) { tom@8193: /* In case the station-type of the load sum is tom@8193: a combined type and we are at non-combined station: tom@8193: we need to add values from previous station of tom@8193: the other type for fractions not given here. */ tom@8193: int otherType = sType == Station.BED_LOAD ? tom@8193: Station.SUSPENDED : Station.BED_LOAD; tom@8193: Station prev = station.prevByType(otherType); 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: 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: 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]; tom@8176: if (Double.isNaN(km)) { 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()) { tom@8176: result[0][i] = entry.getKey(); tom@8176: result[1][i] = entry.getValue().avg(); teichmann@8067: ++i; teichmann@8067: } teichmann@8067: tom@8176: return result; teichmann@8067: } teichmann@8067: teichmann@8049: } teichmann@8049: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf-8 :