teichmann@5863: /* Copyright (C) 2011, 2012, 2013 by Bundesanstalt für Gewässerkunde teichmann@5863: * Software engineering by Intevation GmbH teichmann@5863: * teichmann@5994: * This file is Free Software under the GNU AGPL (>=v3) teichmann@5863: * and comes with ABSOLUTELY NO WARRANTY! Check out the teichmann@5994: * documentation coming with Dive4Elements River for details. teichmann@5863: */ teichmann@5863: teichmann@5831: package org.dive4elements.river.artifacts.model.minfo; rrenkert@4521: teichmann@7460: import java.io.Serializable; rrenkert@4521: import java.util.ArrayList; rrenkert@4521: import java.util.Collections; teichmann@8064: import java.util.Comparator; teichmann@8063: import java.util.Iterator; rrenkert@4521: import java.util.List; rrenkert@4521: import java.util.Map; rrenkert@4521: import java.util.Set; felix@7200: import java.util.TreeMap; rrenkert@4521: teichmann@7460: import org.apache.log4j.Logger; teichmann@7460: import org.dive4elements.artifacts.common.utils.Config; rrenkert@4521: rrenkert@4521: felix@7200: /** Sediment Densities for multiple years. */ felix@7414: public class SedimentDensity implements Serializable rrenkert@4521: { rrenkert@4521: private static final Logger logger = Logger rrenkert@4521: .getLogger(SedimentDensity.class); rrenkert@4521: teichmann@7460: public static final double DEFAULT_SEDIMNET_DENSITY_FACTOR = 1.9; teichmann@7460: teichmann@7460: public static String SEDIMENT_DENSITY_FACTOR_XPATH = teichmann@7460: "/artifact-database/options/sediment-density-factor/text()"; teichmann@7460: teichmann@7460: public static final double SEDIMNET_DENSITY_FACTOR = teichmann@7460: getSedimentDensityFactor(); teichmann@7460: teichmann@8063: private TreeMap> densities; rrenkert@4521: teichmann@7460: teichmann@7460: /** Figures out the sediment density factor from global config. */ teichmann@7460: private static final double getSedimentDensityFactor() { teichmann@7460: teichmann@7460: double factor = DEFAULT_SEDIMNET_DENSITY_FACTOR; teichmann@7460: teichmann@7460: String factorString = teichmann@7460: Config.getStringXPath(SEDIMENT_DENSITY_FACTOR_XPATH); teichmann@7460: teichmann@7460: if (factorString != null) { teichmann@7460: try { teichmann@7460: factor = Double.parseDouble(factorString.trim()); teichmann@7460: } teichmann@7460: catch (NumberFormatException nfe) { teichmann@7460: logger.error(nfe); teichmann@7460: } teichmann@7460: } teichmann@7460: teichmann@7460: logger.info("Sedmiment density factor: " + factor); teichmann@7460: teichmann@7460: return factor; teichmann@7460: } teichmann@7460: rrenkert@4521: public SedimentDensity() { teichmann@8063: densities = new TreeMap>(); rrenkert@4521: } rrenkert@4521: rrenkert@4521: public Map> getDensities() { rrenkert@4521: return densities; rrenkert@4521: } rrenkert@4521: teichmann@8064: private static final Comparator BY_KM = teichmann@8064: new Comparator() { teichmann@8064: @Override teichmann@8064: public int compare(SedimentDensityValue a, SedimentDensityValue b) { teichmann@8064: double diff = a.getKm() - b.getKm(); teichmann@8064: if (diff < 0.0) return -1; teichmann@8064: if (diff > 0.0) return +1; teichmann@8064: return 0; teichmann@8064: } teichmann@8064: }; teichmann@8064: rrenkert@4521: public void addDensity(double km, double density, int year) { teichmann@7463: teichmann@8065: if (logger.isDebugEnabled()) { teichmann@8065: logger.debug("adding " + year); teichmann@8065: } teichmann@7463: teichmann@7463: Integer key = Integer.valueOf(year); teichmann@7463: teichmann@7463: List list = densities.get(key); teichmann@7463: teichmann@7463: if (list == null) { teichmann@7463: list = new ArrayList(); teichmann@7463: densities.put(key, list); rrenkert@4521: } teichmann@7463: teichmann@8064: // Keep list sorted by km. teichmann@8064: SedimentDensityValue sdv = new SedimentDensityValue(km, density, year); teichmann@8064: int index = Collections.binarySearch(list, sdv, BY_KM); teichmann@8064: teichmann@8064: if (index < 0) { teichmann@8064: // index = -(insertion point) - 1 teichmann@8064: // -(index + 1) = insertion point teichmann@8064: index = -(index + 1); teichmann@8064: } teichmann@8064: teichmann@8064: list.add(index, sdv); rrenkert@4521: } rrenkert@4521: felix@6937: /** felix@6937: * Get the density at year. felix@7200: * Measured densities are valid until the next measurement. felix@7200: * If no measurement was found 1.8 is returned. felix@6937: */ rrenkert@4521: public double getDensity(double km, int year) { teichmann@8063: teichmann@8063: if (densities.isEmpty()) { teichmann@8063: return SEDIMNET_DENSITY_FACTOR; rrenkert@4521: } teichmann@8063: teichmann@8063: if (densities.size() == 1) { teichmann@8063: Map.Entry> entry = densities.firstEntry(); teichmann@8063: return entry.getKey() <= year teichmann@8063: ? getDensityAtKm(entry.getValue(), km) teichmann@8063: : SEDIMNET_DENSITY_FACTOR; teichmann@8063: } teichmann@8064: teichmann@8063: Iterator>> iter = teichmann@8063: densities.entrySet().iterator(); teichmann@8063: teichmann@8063: Map.Entry> last = iter.next(); teichmann@8063: teichmann@8063: while (iter.hasNext()) { teichmann@8063: Map.Entry> current = iter.next(); teichmann@8063: last = current; teichmann@8063: int y1 = last.getKey(); teichmann@8063: int y2 = current.getKey(); teichmann@8063: if (year >= y1 && year < y2) { teichmann@8063: return getDensityAtKm(last.getValue(), km); teichmann@8063: } teichmann@8063: if (year >= y2 && !iter.hasNext()) { teichmann@8063: return getDensityAtKm(current.getValue(), km); rrenkert@4521: } rrenkert@4521: } teichmann@8065: teichmann@7460: return SEDIMNET_DENSITY_FACTOR; rrenkert@4521: } rrenkert@4521: felix@7200: /** Get (sorted) map of km to density of all years. */ felix@7309: public double[][] getAllDensities() felix@7200: { felix@7200: TreeMap map = new TreeMap(); teichmann@8063: // XXX: This looks stupid. teichmann@8063: for (List sdvs: densities.values()) { teichmann@8063: for (SedimentDensityValue sdv: sdvs) { felix@7200: map.put(sdv.getKm(), sdv.getDensity()); felix@7200: } felix@7200: } felix@7200: double[][] points = new double[2][map.keySet().size()]; felix@7200: int i = 0; felix@7200: for (Map.Entry kmDens: map.entrySet()) { felix@7200: points[0][i] = kmDens.getKey(); felix@7309: points[1][i] = kmDens.getValue(); felix@7200: i++; felix@7200: } felix@7200: felix@7200: return points; felix@7200: } felix@7200: felix@7200: /** Get points km,density (sorted by km), for a given year. */ felix@7200: public double[][] getDensities(int year) felix@7200: { teichmann@8063: List list = densities.get(year); teichmann@8063: if (list == null) { teichmann@8063: return new double[2][0]; teichmann@8063: } teichmann@8064: // List is sorted in km. teichmann@8064: double[][] points = new double[2][list.size()]; felix@7200: int i = 0; teichmann@8064: for (SedimentDensityValue sdv: list) { teichmann@8064: points[0][i] = sdv.getKm(); teichmann@8064: points[1][i] = sdv.getDensity(); felix@7200: i++; felix@7200: } felix@7200: felix@7200: return points; felix@7200: } felix@7200: felix@7200: /** Get value at km, interpolated. */ teichmann@8064: private static double getDensityAtKm( rrenkert@4521: List values, rrenkert@4521: double km rrenkert@4521: ) { rrenkert@4521: SedimentDensityValue prev = null; rrenkert@4521: SedimentDensityValue next = null; rrenkert@4521: for (SedimentDensityValue sdv: values) { teichmann@6970: if (Math.abs(sdv.getKm() - km) < 0.00001) { teichmann@6970: return prev.getDensity(); rrenkert@4521: } rrenkert@4521: if (sdv.getKm() > km) { rrenkert@4521: next = sdv; rrenkert@4521: break; rrenkert@4521: } rrenkert@4521: prev = sdv; rrenkert@4521: } teichmann@6970: return spline(prev, next, km); rrenkert@4521: } rrenkert@4521: felix@7200: /** Linearly interpolate between density values. */ teichmann@6970: private static double spline( rrenkert@4521: SedimentDensityValue prev, rrenkert@4521: SedimentDensityValue next, rrenkert@4521: double km rrenkert@4521: ) { teichmann@6970: if (prev == null && next == null) { teichmann@6970: logger.warn("prev and next are null -> NaN"); teichmann@6970: return Double.NaN; teichmann@6970: } teichmann@6970: teichmann@6970: if (prev == null) return next.getDensity(); teichmann@6970: if (next == null) return prev.getDensity(); teichmann@6970: teichmann@6970: // XXX: This is no spline interpolation! rrenkert@4521: double lower = prev.getKm(); rrenkert@4521: double upper = next.getKm(); rrenkert@4521: double upperDensity = next.getDensity(); rrenkert@4521: double lowerDensity = prev.getDensity(); rrenkert@4521: teichmann@6970: double m = (upperDensity - lowerDensity)/(upper - lower); felix@6949: double b = lowerDensity - (m * lower); teichmann@6970: return m * km + b; rrenkert@4521: } rrenkert@4521: felix@6948: felix@6948: /** If multiple values for same year and station are found, felix@6948: * build and store average, dismiss multiple values. */ rrenkert@4521: public void cleanUp() { rrenkert@4521: Set keys = densities.keySet(); felix@6948: // Walk over years rrenkert@4521: for (Integer key : keys) { rrenkert@4521: List list = densities.get(key); rrenkert@4521: if (list.size() == 0) { rrenkert@4521: return; rrenkert@4521: } rrenkert@4521: List cleaned = rrenkert@4521: new ArrayList(); rrenkert@4521: double prevkm = list.get(0).getKm(); rrenkert@4521: int counter = 0; rrenkert@4521: double sum = 0d; rrenkert@4521: for (SedimentDensityValue value : list) { felix@6948: // Apparently we can assume that values are ordered by km. rrenkert@4521: if (value.getKm() == prevkm) { rrenkert@4521: sum += value.getDensity(); rrenkert@4521: counter++; rrenkert@4521: } rrenkert@4521: else { rrenkert@4521: cleaned.add(new SedimentDensityValue( rrenkert@4521: prevkm, rrenkert@4521: sum / counter, rrenkert@4521: value.getYear())); rrenkert@4521: sum = value.getDensity(); rrenkert@4521: counter = 1; rrenkert@4521: } rrenkert@4521: prevkm = value.getKm(); rrenkert@4521: } rrenkert@4521: this.densities.put(key, cleaned); rrenkert@4521: } rrenkert@4521: } rrenkert@4521: } felix@6948: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf-8 :