view artifacts/src/main/java/org/dive4elements/river/artifacts/model/minfo/SedimentDensity.java @ 7300:83bb52fa0c32

(issue1529) Be more tolerant in the fitting. The invalid value warning is removed because invalid data is expected there when datapoints are not valid for this KM
author Andre Heinecke <aheinecke@intevation.de>
date Fri, 11 Oct 2013 18:40:33 +0200
parents fe32a7f9655e
children a6ceb4b333c3
line wrap: on
line source
/* Copyright (C) 2011, 2012, 2013 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.Collections;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
import java.util.Set;
import java.util.TreeMap;

import org.apache.log4j.Logger;


/** Sediment Densities for multiple years. */
public class SedimentDensity
{
    private static final Logger logger = Logger
        .getLogger(SedimentDensity.class);

    private Map<Integer, List<SedimentDensityValue>> densities;
    private List<Integer> years;

    public SedimentDensity() {
        this.densities = new HashMap<Integer, List<SedimentDensityValue>>();
        this.years = new ArrayList<Integer>();
    }

    public Map<Integer, List<SedimentDensityValue>> getDensities() {
        return densities;
    }

    public void setDensities(Map<Integer, List<SedimentDensityValue>> densities) {
        this.densities = densities;
    }

    public void addDensity(double km, double density, int year) {
        logger.debug("adding " + year);
        if (this.densities.containsKey(year)) {
            List<SedimentDensityValue> list = this.densities.get(year);
            list.add(new SedimentDensityValue(km, density, year));
        }
        else {
            List<SedimentDensityValue> list =
                new ArrayList<SedimentDensityValue>();
            list.add(new SedimentDensityValue(km, density, year));
            densities.put(year, list);
        }
        if (!this.years.contains(new Integer(year))) {
            logger.debug("new year");
            years.add(new Integer(year));
        }
    }

    public List<Integer> getYears() {
        return years;
    }

    public void setYears(List<Integer> years) {
        this.years = years;
    }

    /**
     * Get the density at year.
     * Measured densities are valid until the next measurement.
     * If no measurement was found 1.8 is returned.
     */
    public double getDensity(double km, int year) {
        Collections.sort(this.years);
        if (this.years.size() == 1 && years.get(0) <= year) {
            logger.debug("get density from year " + year + " at km " + km);
            return getDensityAtKm(densities.get(years.get(0)), km);
        }
        else if (this.years.size() > 1) {
            for (int i = 0, I = years.size()-1; i < I; i++) {
                int y1 = years.get(i);
                int y2 = years.get(i + 1);
                if (year >= y1 && year < y2) {
                    return getDensityAtKm(densities.get(y1), km);
                }
                else if (year >= y2 && i == years.size() -1) {
                    return getDensityAtKm(densities.get(y2), km);
                }
            }
        }
        return 1.8d;
    }

    /** Get (sorted) map of km to density of all years. */
    protected double[][] getAllDensities()
    {
        TreeMap<Double, Double> map = new TreeMap<Double,Double>();
        for (int year: years) {
            for (SedimentDensityValue sdv: densities.get(year)) {
                map.put(sdv.getKm(), sdv.getDensity());
            }
        }
        double[][] points = new double[2][map.keySet().size()];
        int i = 0;
        for (Map.Entry<Double, Double> kmDens: map.entrySet()) {
            points[0][i] = kmDens.getKey();
            points[2][i] = kmDens.getValue();
            i++;
        }

        return points;
    }

    /** Get points  km,density (sorted by km), for a given year. */
    public double[][] getDensities(int year)
    {
        TreeMap<Double, Double> map = new TreeMap<Double,Double>();
        for (SedimentDensityValue sdv: densities.get(year)) {
            map.put(sdv.getKm(), sdv.getDensity());
        }
        double[][] points = new double[2][map.keySet().size()];
        int i = 0;
        for (Map.Entry<Double, Double> kmDens: map.entrySet()) {
            points[0][i] = kmDens.getKey();
            points[1][i] = kmDens.getValue();
            i++;
        }

        return points;
    }

    /** Get value at km, interpolated. */
    private double getDensityAtKm(
        List<SedimentDensityValue> values,
        double km
    ) {
        SedimentDensityValue prev = null;
        SedimentDensityValue next = null;
        for (SedimentDensityValue sdv: values) {
            logger.debug("year: " + sdv.getYear());
            if (Math.abs(sdv.getKm() - km) < 0.00001) {
                return prev.getDensity();
            }
            if (sdv.getKm() > km) {
                next = sdv;
                break;
            }
            prev = sdv;
        }
        return spline(prev, next, km);
    }

    /** Linearly interpolate between density values. */
    private static double spline(
        SedimentDensityValue prev,
        SedimentDensityValue next,
        double km
    ) {
        if (prev == null && next == null) {
            logger.warn("prev and next are null -> NaN");
            return Double.NaN;
        }

        if (prev == null) return next.getDensity();
        if (next == null) return prev.getDensity();

        // XXX: This is no spline interpolation!
        double lower = prev.getKm();
        double upper = next.getKm();
        double upperDensity = next.getDensity();
        double lowerDensity = prev.getDensity();

        double m = (upperDensity - lowerDensity)/(upper - lower);
        double b = lowerDensity - (m * lower);
        return m * km + b;
    }


    /** If multiple values for same year and station are found,
     * build and store average, dismiss multiple values. */
    public void cleanUp() {
        Set<Integer> keys = densities.keySet();
        // Walk over years
        for (Integer key : keys) {
            List<SedimentDensityValue> list = densities.get(key);
            if (list.size() == 0) {
                return;
            }
            List<SedimentDensityValue> cleaned =
                new ArrayList<SedimentDensityValue>();
            double prevkm = list.get(0).getKm();
            int counter = 0;
            double sum = 0d;
            for (SedimentDensityValue value : list) {
                // Apparently we can assume that values are ordered by km.
                if (value.getKm() == prevkm) {
                    sum += value.getDensity();
                    counter++;
                }
                else {
                    cleaned.add(new SedimentDensityValue(
                        prevkm,
                        sum / counter,
                        value.getYear()));
                    sum = value.getDensity();
                    counter = 1;
                }
                prevkm = value.getKm();
            }
            this.densities.put(key, cleaned);
        }
    }
}
// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf-8 :

http://dive4elements.wald.intevation.org