changeset 6970:7be97faf5848

flys/issue1235: Same kicks against a few inconsistencies and bugs in the calculation of sediment loads. I believe it _do_not_ delivers the right results.
author Sascha L. Teichmann <teichmann@intevation.de>
date Thu, 05 Sep 2013 17:15:04 +0200
parents c137f5028591
children 164e2f2c9bea
files artifacts/src/main/java/org/dive4elements/river/artifacts/model/minfo/SedimentDensity.java artifacts/src/main/java/org/dive4elements/river/artifacts/model/minfo/SedimentLoad.java artifacts/src/main/java/org/dive4elements/river/artifacts/model/minfo/SedimentLoadCalculation.java artifacts/src/main/java/org/dive4elements/river/artifacts/model/minfo/SedimentLoadFactory.java artifacts/src/main/java/org/dive4elements/river/artifacts/model/minfo/SedimentLoadFraction.java artifacts/src/main/java/org/dive4elements/river/exports/minfo/SedimentLoadExporter.java
diffstat 6 files changed, 122 insertions(+), 161 deletions(-) [+]
line wrap: on
line diff
--- a/artifacts/src/main/java/org/dive4elements/river/artifacts/model/minfo/SedimentDensity.java	Thu Sep 05 15:47:24 2013 +0200
+++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/minfo/SedimentDensity.java	Thu Sep 05 17:15:04 2013 +0200
@@ -95,15 +95,12 @@
         List<SedimentDensityValue> values,
         double km
     ) {
-        boolean found = false;
         SedimentDensityValue prev = null;
         SedimentDensityValue next = null;
         for (SedimentDensityValue sdv: values) {
             logger.debug("year: " + sdv.getYear());
-            if (sdv.getKm() == km) {
-                prev = sdv;
-                found = true;
-                break;
+            if (Math.abs(sdv.getKm() - km) < 0.00001) {
+                return prev.getDensity();
             }
             if (sdv.getKm() > km) {
                 next = sdv;
@@ -111,27 +108,31 @@
             }
             prev = sdv;
         }
-        if (found) {
-            return prev.getDensity();
-        }
-        else {
-            return spline(prev, next, km);
-        }
+        return spline(prev, next, km);
     }
 
-    private double spline(
+    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 m = (upperDensity - lowerDensity)/(upper - lower);
         double b = lowerDensity - (m * lower);
-        return (m * km) + b;
+        return m * km + b;
     }
 
 
--- a/artifacts/src/main/java/org/dive4elements/river/artifacts/model/minfo/SedimentLoad.java	Thu Sep 05 15:47:24 2013 +0200
+++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/minfo/SedimentLoad.java	Thu Sep 05 17:15:04 2013 +0200
@@ -10,10 +10,13 @@
 
 import java.util.Date;
 import java.util.HashMap;
+import java.util.Map;
 import java.util.Set;
+import java.util.TreeMap;
 
 import org.dive4elements.river.artifacts.model.NamedObjectImpl;
 import org.dive4elements.river.artifacts.model.Range;
+import org.dive4elements.river.utils.EpsilonComparator;
 
 
 /** Gives access to Fractions (at kms). */
@@ -26,10 +29,10 @@
     protected boolean isEpoch;
     protected String unit;
 
-    protected HashMap<Double, SedimentLoadFraction> kms;
+    protected Map<Double, SedimentLoadFraction> kms;
 
     public SedimentLoad() {
-        kms = new HashMap<Double, SedimentLoadFraction>();
+        kms = new TreeMap<Double, SedimentLoadFraction>(EpsilonComparator.CMP);
     }
 
     public SedimentLoad(
@@ -88,87 +91,49 @@
     }
 
     public SedimentLoadFraction getFraction(double km) {
-        if (kms.get(km) == null) {
-            return new SedimentLoadFraction();
+        SedimentLoadFraction f = kms.get(km);
+        if (f == null) {
+            f = new SedimentLoadFraction();
+            kms.put(km, f);
         }
-        return kms.get(km);
+        return f;
     }
 
     public void setCoarse(double km, double coarse, Range range) {
-        if (kms.containsKey(km)) {
-            kms.get(km).setCoarse(coarse);
-        }
-        else {
-            SedimentLoadFraction f = new SedimentLoadFraction();
-            f.setCoarse(coarse);
-            f.setCoarseRange(range);
-            kms.put(km, f);
-        }
+        SedimentLoadFraction f = getFraction(km);
+        f.setCoarse(coarse);
+        f.setCoarseRange(range);
     }
 
     public void setFineMiddle(double km, double fine_middle, Range range) {
-        if (kms.containsKey(km)) {
-            kms.get(km).setFineMiddle(fine_middle);
-            kms.get(km).setFineMiddleRange(range);
-        }
-        else {
-            SedimentLoadFraction f = new SedimentLoadFraction();
-            f.setFineMiddle(fine_middle);
-            f.setFineMiddleRange(range);
-            kms.put(km, f);
-        }
+        SedimentLoadFraction f = getFraction(km);
+        f.setFineMiddle(fine_middle);
+        f.setFineMiddleRange(range);
     }
 
+
     public void setSand(double km, double sand, Range range) {
-        if (kms.containsKey(km)) {
-            kms.get(km).setSand(sand);
-            kms.get(km).setSandRange(range);
-        }
-        else {
-            SedimentLoadFraction f = new SedimentLoadFraction();
-            f.setSand(sand);
-            f.setSandRange(range);
-            kms.put(km, f);
-        }
+        SedimentLoadFraction f = getFraction(km);
+        f.setSand(sand);
+        f.setSandRange(range);
     }
 
     public void setSuspSand(double km, double susp_sand, Range range) {
-        if (kms.containsKey(km)) {
-            kms.get(km).setSuspSand(susp_sand);
-            kms.get(km).setSuspSandRange(range);
-        }
-        else {
-            SedimentLoadFraction f = new SedimentLoadFraction();
-            f.setSuspSand(susp_sand);
-            f.setSuspSandRange(range);
-            kms.put(km, f);
-        }
+        SedimentLoadFraction f = getFraction(km);
+        f.setSuspSand(susp_sand);
+        f.setSuspSandRange(range);
     }
 
     public void setSuspSandBed(double km, double susp_sand_bed, Range range) {
-        if (kms.containsKey(km)) {
-            kms.get(km).setSuspSandBed(susp_sand_bed);
-            kms.get(km).setSuspSandBedRange(range);
-        }
-        else {
-            SedimentLoadFraction f = new SedimentLoadFraction();
-            f.setSuspSandBed(susp_sand_bed);
-            f.setSuspSandBedRange(range);
-            kms.put(km, f);
-        }
+        SedimentLoadFraction f = getFraction(km);
+        f.setSuspSandBed(susp_sand_bed);
+        f.setSuspSandBedRange(range);
     }
 
     public void setSuspSediment(double km, double susp_sediment, Range range) {
-        if (kms.containsKey(km)) {
-            kms.get(km).setSuspSediment(susp_sediment);
-            kms.get(km).setSuspSedimentRange(range);
-        }
-        else {
-            SedimentLoadFraction f = new SedimentLoadFraction();
-            f.setSuspSediment(susp_sediment);
-            f.setSuspSedimentRange(range);
-            kms.put(km, f);
-        }
+        SedimentLoadFraction f = getFraction(km);
+        f.setSuspSediment(susp_sediment);
+        f.setSuspSedimentRange(range);
     }
 
     public void setLoadTotal(double km, double total) {
@@ -176,42 +141,21 @@
     }
 
     public void setLoadTotal(double km, double total, Range range) {
-        if (kms.containsKey(km)) {
-            kms.get(km).setLoadTotal(total);
-            kms.get(km).setLoadTotalRange(range);
-        }
-        else {
-            SedimentLoadFraction f = new SedimentLoadFraction();
-            f.setLoadTotal(total);
-            f.setLoadTotalRange(range);
-            kms.put(km, f);
-        }
+        SedimentLoadFraction f = getFraction(km);
+        f.setLoadTotal(total);
+        f.setLoadTotalRange(range);
     }
 
     public void setTotal(double km, double total, Range range) {
-        if (kms.containsKey(km)) {
-            kms.get(km).setTotal(total);
-            kms.get(km).setTotalRange(range);
-        }
-        else {
-            SedimentLoadFraction f = new SedimentLoadFraction();
-            f.setTotal(total);
-            f.setTotalRange(range);
-            kms.put(km, f);
-        }
+        SedimentLoadFraction f = getFraction(km);
+        f.setTotal(total);
+        f.setTotalRange(range);
     }
 
     public void setUnknown(double km, double unknown, Range range) {
-        if (kms.containsKey(km)) {
-            kms.get(km).setUnknown(unknown);
-            kms.get(km).setUnknownRange(range);
-        }
-        else {
-            SedimentLoadFraction f = new SedimentLoadFraction();
-            f.setUnknown(unknown);
-            f.setUnknownRange(range);
-            kms.put(km, f);
-        }
+        SedimentLoadFraction f = getFraction(km);
+        f.setUnknown(unknown);
+        f.setUnknownRange(range);
     }
 
     public String getUnit() {
--- a/artifacts/src/main/java/org/dive4elements/river/artifacts/model/minfo/SedimentLoadCalculation.java	Thu Sep 05 15:47:24 2013 +0200
+++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/minfo/SedimentLoadCalculation.java	Thu Sep 05 17:15:04 2013 +0200
@@ -338,35 +338,6 @@
     }
 
 
-    /** Returns true if all fraction values except SuspSediment are unset. */
-    private boolean hasOnlySuspValues(SedimentLoadFraction fraction) {
-        return (fraction.getSuspSediment() != 0d &&
-            fraction.getCoarse() == 0d &&
-            fraction.getFineMiddle() == 0d &&
-            fraction.getSand() == 0d &&
-            fraction.getSuspSand() == 0d);
-    }
-
-
-    /** Returns true if all fraction values except SuspSediment are set. */
-    private boolean hasButSuspValues(SedimentLoadFraction fraction) {
-        return (fraction.getSuspSediment() == 0d &&
-            fraction.getCoarse() != 0d &&
-            fraction.getFineMiddle() != 0d &&
-            fraction.getSand() != 0d &&
-            fraction.getSuspSand() != 0d);
-    }
-
-
-    /** Returns true if all fraction needed for total calculation are set. */
-    private boolean complete(SedimentLoadFraction fraction) {
-        return (fraction.getCoarse() != 0d &&
-                fraction.getFineMiddle() != 0d &&
-                fraction.getSand() != 0d &&
-                fraction.getSuspSand() != 0d &&
-                fraction.getSuspSediment() != 0d);
-    }
-
 
     /**
      * Set total values in load.
@@ -389,12 +360,10 @@
         Range lastSuspRange = null;
         double lastSuspValue = 0d;
 
-        TreeSet<Double> kms = new TreeSet<Double>(load.getKms());
-
-        for (double km: kms) {
+        for (double km: load.getKms()) { // kms are already sorted!
             logger.debug ("Trying to add at km " + km);
             SedimentLoadFraction fraction = load.getFraction(km);
-            if (complete(fraction)) {
+            if (fraction.isComplete()) {
                 double total = fraction.getCoarse() +
                     fraction.getFineMiddle() +
                     fraction.getSand() +
@@ -430,7 +399,7 @@
                     }
                 }
             }
-            else if (hasOnlySuspValues(fraction) && lastOtherRange != null) {
+            else if (fraction.hasOnlySuspValues() && lastOtherRange != null) {
                 // Split stuff.
                 Range suspSedimentRange = fraction.getSuspSedimentRange();
                 // if intersects with last other range, cool! merge and add!
@@ -455,7 +424,8 @@
                         lastOtherRange.setStart(suspSedimentRange.getEnd());
                         lastSuspRange = null;
                     }
-                    if (Math.abs(suspSedimentRange.getEnd() - lastOtherRange.getEnd()) < 0.1d) {
+                    if (lastOtherRange != null
+                    && Math.abs(suspSedimentRange.getEnd() - lastOtherRange.getEnd()) < 0.1d) {
                         lastOtherRange = null;
                         lastSuspRange = null;
                     }
@@ -467,7 +437,7 @@
                     lastOtherRange = null;
                 }
             }
-            else if (hasButSuspValues(fraction) && lastSuspRange != null) {
+            else if (fraction.hasButSuspValues() && lastSuspRange != null) {
                 // If intersects with last suspsed range, merge and add
                 double total = fraction.getCoarse() +
                     fraction.getFineMiddle() +
@@ -492,7 +462,9 @@
                         lastSuspRange = null;
                         lastOtherValue = total - lastSuspValue;
                     }
-                    if (lastSuspRange != null && Math.abs(lastSuspRange.getEnd() - lastOtherRange.getEnd()) < 0.1d) {
+                    if (lastSuspRange != null
+                    &&  lastOtherRange != null
+                    && Math.abs(lastSuspRange.getEnd() - lastOtherRange.getEnd()) < 0.1d) {
                         lastOtherRange = null;
                         lastSuspRange = null;
                     }
@@ -508,7 +480,7 @@
             else {
                 // Some values are missing or no intersection with former values.
                 // Stay as we are.
-                if (hasButSuspValues(fraction)) {
+                if (fraction.hasButSuspValues()) {
                     double total = fraction.getCoarse() +
                         fraction.getFineMiddle() +
                         fraction.getSand() +
@@ -517,7 +489,7 @@
                     lastOtherValue = total;
                     lastSuspRange = null;
                 }
-                else if (hasOnlySuspValues(fraction)) {
+                else if (fraction.hasOnlySuspValues()) {
                     lastSuspRange = fraction.getSuspSedimentRange();
                     lastSuspValue = fraction.getSuspSediment();
                     lastOtherRange = null;
--- a/artifacts/src/main/java/org/dive4elements/river/artifacts/model/minfo/SedimentLoadFactory.java	Thu Sep 05 15:47:24 2013 +0200
+++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/minfo/SedimentLoadFactory.java	Thu Sep 05 17:15:04 2013 +0200
@@ -13,6 +13,7 @@
 import java.util.Calendar;
 import java.util.Date;
 import java.util.List;
+import java.util.Map;
 import java.util.TreeMap;
 
 import net.sf.ehcache.Cache;
@@ -508,15 +509,24 @@
         TreeMap<Double, MeasurementStation> stations,
         double km
     ) {
-        MeasurementStation station = stations.floorEntry(km).getValue();
-        if (station == null || !station.getRange().contains(km)) {
+        Map.Entry<Double, MeasurementStation> entry = stations.floorEntry(km);
+        if (entry == null) {
+            return null;
+        }
+        MeasurementStation station = entry.getValue();
+        if (station == null
+        || station.getRange() == null
+        || !station.getRange().contains(km)) {
             return null;
         }
 
         double endKm;
 
-        if (stations.ceilingEntry(km + 0.1d) != null) {
-            MeasurementStation nextStation = stations.ceilingEntry(km + 0.1d).getValue();
+        Map.Entry<Double, MeasurementStation> ceilingEntry =
+            stations.ceilingEntry(km + 0.1d);
+
+        if (ceilingEntry != null) {
+            MeasurementStation nextStation = ceilingEntry.getValue();
             endKm = nextStation.getRange().getA().doubleValue();
         }
         else {
--- a/artifacts/src/main/java/org/dive4elements/river/artifacts/model/minfo/SedimentLoadFraction.java	Thu Sep 05 15:47:24 2013 +0200
+++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/minfo/SedimentLoadFraction.java	Thu Sep 05 17:15:04 2013 +0200
@@ -25,15 +25,15 @@
     double total;
     double unknown;
     /** Values are valid within this km range. */
-    Range sandRange = null;
-    Range fineMiddleRange = null;
-    Range coarseRange = null;
-    Range suspSandRange = null;
-    Range suspSandBedRange = null;
-    Range suspSedimentRange = null;
-    Range loadTotalRange = null;
-    Range totalRange = null;
-    Range unknownRange = null;
+    Range sandRange;
+    Range fineMiddleRange;
+    Range coarseRange;
+    Range suspSandRange;
+    Range suspSandBedRange;
+    Range suspSedimentRange;
+    Range loadTotalRange;
+    Range totalRange;
+    Range unknownRange;
 
     public SedimentLoadFraction() {
         sand = 0d;
@@ -189,5 +189,35 @@
     public void setUnknownRange(Range unknownRange) {
         this.unknownRange = unknownRange;
     }
+
+    /** Returns true if all fraction values except SuspSediment are unset. */
+    public boolean hasOnlySuspValues() {
+        return
+            getSuspSediment() != 0d &&
+            getCoarse() == 0d &&
+            getFineMiddle() == 0d &&
+            getSand() == 0d &&
+            getSuspSand() == 0d;
+    }
+
+    /** Returns true if all fraction values except SuspSediment are set. */
+    public boolean hasButSuspValues() {
+        return
+            getSuspSediment() == 0d &&
+            getCoarse() != 0d &&
+            getFineMiddle() != 0d &&
+            getSand() != 0d &&
+            getSuspSand() != 0d;
+    }
+
+    /** Returns true if all fraction needed for total calculation are set. */
+    public boolean isComplete() {
+        return
+            getCoarse() != 0d &&
+            getFineMiddle() != 0d &&
+            getSand() != 0d &&
+            getSuspSand() != 0d &&
+            getSuspSediment() != 0d;
+    }
 }
 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf-8 :
--- a/artifacts/src/main/java/org/dive4elements/river/exports/minfo/SedimentLoadExporter.java	Thu Sep 05 15:47:24 2013 +0200
+++ b/artifacts/src/main/java/org/dive4elements/river/exports/minfo/SedimentLoadExporter.java	Thu Sep 05 17:15:04 2013 +0200
@@ -94,8 +94,12 @@
     protected void writeCSVData(CSVWriter writer) throws IOException {
         writeCSVHeader(writer);
 
+        //@felix: Please fix this!!!!
+        //boolean asSingleYear = toYear == 0;
+        boolean asSingleYear = true;
+
         for (SedimentLoadResult result: results) {
-            String years = (toYear == 0)
+            String years = asSingleYear
                     ? result.getStartYear()+ " "
                     : result.getStartYear() + "-" + result.getEndYear();
             SedimentLoad load = result.getLoad();

http://dive4elements.wald.intevation.org