changeset 433:828df3ddb758

Added interpolation capabilities along z axis to XYColumns. gnv-artifacts/trunk@481 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Wed, 23 Dec 2009 12:20:25 +0000
parents 6a70e8883307
children 0eed5749fd63
files gnv-artifacts/ChangeLog gnv-artifacts/src/main/java/de/intevation/gnv/math/AttributedXYColumns.java gnv-artifacts/src/main/java/de/intevation/gnv/math/ConstantFunction.java gnv-artifacts/src/main/java/de/intevation/gnv/math/HeightValue.java gnv-artifacts/src/main/java/de/intevation/gnv/math/LinearFunction.java gnv-artifacts/src/main/java/de/intevation/gnv/math/XYColumn.java gnv-artifacts/src/main/java/de/intevation/gnv/math/XYDepth.java
diffstat 7 files changed, 209 insertions(+), 10 deletions(-) [+]
line wrap: on
line diff
--- a/gnv-artifacts/ChangeLog	Wed Dec 23 09:45:40 2009 +0000
+++ b/gnv-artifacts/ChangeLog	Wed Dec 23 12:20:25 2009 +0000
@@ -1,7 +1,43 @@
+2009-12-23	Sascha L. Teichmann	<sascha.teichmann@intevation.de>
+
+	* src/main/java/de/intevation/gnv/math/XYColumn.java: Added a method
+	  prepare() which generates an continues interpolator for the
+	  given z values. if only one value is given a constant function
+	  is assumed. If the larger z-value is below zero the
+	  next lower value is supplemented at zero. Symmetrically 
+	  if the lowest z-value is above the depth at the given point
+	  the lowest value is repeated at depth. This should guarantee
+	  that the gradient is vansihing towards the surface and the
+	  bottom of the ocean.
+
+	  If after the supplementation there are less than three points
+	  a linear interpolation is performed. If there are more than three
+	  points a higher degree interpolation is used instead. This defaults
+	  to a cubic spline interpolation. Overwrite the getInterpolator()
+	  function to replace this behavior.
+	  
+	* src/main/java/de/intevation/gnv/math/ConstantFunction.java: New.
+	  Constant function used in interpolation.
+
+	* src/main/java/de/intevation/gnv/math/LinearFunction.java: Added
+	  an inner class Univariate which fits into the interpolation
+	  framework.
+
+	* src/main/java/de/intevation/gnv/math/HeightValue.java: Sort
+	  z-Values in descending order because we are below zero.
+
+	* src/main/java/de/intevation/gnv/math/XYDepth.java: New. Interface
+	  to figure out the depth (negative values below surface) for
+	  a given coordinate. TODO: Implement this by query the DEM grid
+	  of the ocean.
+
+	* src/main/java/de/intevation/gnv/math/AttributedXYColumns.java: Added
+	  authors.
+
 2009-12-23	Sascha L. Teichmann	<sascha.teichmann@intevation.de>
 
 	* src/main/java/de/intevation/gnv/state/profile/verticalcrosssection/VerticalCrossSectionOutputState.java:
-	  When preprocessing database data only dissamble WKT points if we
+	  When preprocessing database data only dissemble WKT points if we
 	  have to.
 	  Read z values as double value now.
 	  Commented out CSV export because it takes the database data
--- a/gnv-artifacts/src/main/java/de/intevation/gnv/math/AttributedXYColumns.java	Wed Dec 23 09:45:40 2009 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/AttributedXYColumns.java	Wed Dec 23 12:20:25 2009 +0000
@@ -6,6 +6,10 @@
 
 import java.io.Serializable;
 
+/**
+ *  @author Ingo Weinzierl
+ *  @author Sascha L. Teichmann
+ */
 public class AttributedXYColumns
 implements   Serializable
 {
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/ConstantFunction.java	Wed Dec 23 12:20:25 2009 +0000
@@ -0,0 +1,24 @@
+package de.intevation.gnv.math;
+
+import org.apache.commons.math.analysis.UnivariateRealFunction;
+
+/**
+ * @author Sascha L. Teichmann <sascha.teichmann@intevation.de>
+ */
+public class ConstantFunction
+implements   UnivariateRealFunction
+{
+    protected double value;
+    
+    public ConstantFunction() {
+    }
+
+    public ConstantFunction(double value) {
+        this.value = value;
+    }
+
+    public double value(double x) {
+        return value;
+    }
+}
+// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8:
--- a/gnv-artifacts/src/main/java/de/intevation/gnv/math/HeightValue.java	Wed Dec 23 09:45:40 2009 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/HeightValue.java	Wed Dec 23 12:20:25 2009 +0000
@@ -2,16 +2,28 @@
 
 import java.io.Serializable;
 
+import java.util.Comparator;
+
 /**
  * @author Ingo Weinzierl <ingo.weinzierl@intevation.de>
+ * @author Sascha L. Teichmann <sascha.teichmann@intevation.de>
  */
 public class HeightValue
 implements   Serializable
 {
-    private double z;
-    private double v;
-    private int    k;
+    public static final Comparator INV_Z_COMPARATOR = new Comparator() {
+        public int compare(Object a, Object b) {
+            HeightValue ha = (HeightValue)a;
+            HeightValue hb = (HeightValue)b;
+            if (ha.z > hb.z) return -1;
+            if (ha.z < hb.z) return +1;
+            return 0;
+        }
+    };
 
+    public double z;
+    public double v;
+    public int    k;
 
     public HeightValue(double z, double v, int k) {
         this.z = z;
@@ -19,18 +31,16 @@
         this.k = k;
     }
 
-
     public double getZ() {
         return z;
     }
 
-
     public double getV() {
         return v;
     }
 
-
     public double getK() {
         return k;
     }
 }
+// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8:
--- a/gnv-artifacts/src/main/java/de/intevation/gnv/math/LinearFunction.java	Wed Dec 23 09:45:40 2009 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/LinearFunction.java	Wed Dec 23 12:20:25 2009 +0000
@@ -4,6 +4,8 @@
 
 import org.apache.commons.math.FunctionEvaluationException;
 
+import org.apache.commons.math.analysis.UnivariateRealFunction;
+
 /**
  *  @author Sascha L. Teichmann
  */
@@ -12,6 +14,31 @@
 {
     public static final LinearFunction INSTANCE = new LinearFunction();
 
+    public static class Univariate 
+    implements          UnivariateRealFunction
+    {
+        protected double m;
+        protected double b;
+
+        public Univariate() {
+        }
+
+        public Univariate(double x1, double y1, double x2, double y2) {
+            if (y1 == y2) {
+                m = 0d;
+                b = (x1 + x2)*0.5d;
+            }
+            else {
+                m = (x1 - x2)/(y1 - y2);
+                b = y1 - m*x1;
+            }
+        }
+
+        public double value(double x) {
+            return m*x + b;
+        }
+    } // class Univariate
+
     public LinearFunction() {
     }
 
--- a/gnv-artifacts/src/main/java/de/intevation/gnv/math/XYColumn.java	Wed Dec 23 09:45:40 2009 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/XYColumn.java	Wed Dec 23 12:20:25 2009 +0000
@@ -2,6 +2,19 @@
 
 import java.util.ArrayList;
 import java.util.List;
+import java.util.Collections;
+
+import org.apache.commons.math.analysis.interpolation.SplineInterpolator;
+import org.apache.commons.math.analysis.interpolation.UnivariateRealInterpolator;
+
+import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction;
+
+import org.apache.commons.math.analysis.UnivariateRealFunction;
+
+import org.apache.commons.math.MathException;
+import org.apache.commons.math.FunctionEvaluationException;
+
+import org.apache.log4j.Logger;
 
 /**
  * @author Ingo Weinzierl <ingo.weinzierl@intevation.de>
@@ -9,25 +22,98 @@
  */
 public class XYColumn
 extends      Point2d
+implements   UnivariateRealFunction
 {
+    private static Logger log = Logger.getLogger(Interpolation2D.class);
+
     protected List<HeightValue> values;
 
+    protected transient UnivariateRealFunction curve;
+
     public XYColumn() {
+        values = new ArrayList<HeightValue>();
     }
 
     public XYColumn(double x, double y, int i, int j) {
         super(x, y, i, j);
+        values = new ArrayList<HeightValue>();
     }
 
     public void add(HeightValue value) {
-        if (values == null) {
-            values = new ArrayList<HeightValue>();
-        }
         values.add(value);
     }
 
     public List<HeightValue> getValues() {
         return values;
     }
+
+    public double value(double depth) {
+        try {
+            if (curve != null) {
+                return curve.value(depth);
+            }
+        }
+        catch (FunctionEvaluationException fee) {
+            log.error("evaluation failed", fee);
+        }
+
+        return Double.NaN;
+    }
+
+    public void prepare(XYDepth xyDepth) {
+        int N = values.size();
+        if (curve == null && N > 0) {
+            if (N == 1) {
+                // only one value -> constant function
+                curve = new ConstantFunction(values.get(0).v);
+            }
+            else { // more than on value
+                double depth = xyDepth.depth(this);
+                Collections.sort(values, HeightValue.INV_Z_COMPARATOR);
+
+                // if there is no value at 0 repeat first value
+                HeightValue first = values.get(0);
+                if (first.z < 0d) {
+                    values.add(0, new HeightValue(0d, first.z, first.k-1));
+                }
+
+                // if there is no value at depth repeat last value
+                HeightValue last = values.get(N-1);
+                if (last.z > depth) {
+                    values.add(new HeightValue(depth, last.z, last.k+1));
+                }
+                N = values.size();
+                if (N < 3) { // interpolate linear
+                    first = values.get(0);
+                    last  = values.get(N-1);
+                    curve = new LinearFunction.Univariate(
+                        first.z, first.v,
+                        last.z,  last.v);
+                }
+                else { // higher degree interpolation
+                    double [] z = new double[N];
+                    double [] v = new double[N];
+                    for (int i = 0; i < N; ++i) {
+                        HeightValue h = values.get(i);
+                        z[i] = h.z;
+                        v[i] = h.v;
+                    }
+                    try {
+                        curve = getInterpolator().interpolate(z, v);
+                    }
+                    catch (MathException me) {
+                        log.error("interpolation failed", me);
+                    }
+                }
+            }
+        }
+        else {
+            log.error("no points for interpolation");
+        }
+    }
+
+    protected UnivariateRealInterpolator getInterpolator() {
+        return new SplineInterpolator();
+    }
 }
 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8:
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/XYDepth.java	Wed Dec 23 12:20:25 2009 +0000
@@ -0,0 +1,12 @@
+package de.intevation.gnv.math;
+
+import com.vividsolutions.jts.geom.Coordinate;
+
+/**
+ *  @author Sascha L. Teichmann
+ */
+public interface XYDepth
+{
+    double depth(Coordinate coordinate);
+}
+// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8:

http://dive4elements.wald.intevation.org