changeset 474:ab29e4ff2fda

Added area interpolation needed for "Horizontalschnitt" gnv-artifacts/trunk@540 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Thu, 14 Jan 2010 10:34:05 +0000
parents a6a33ef35809
children c0504976e606
files gnv-artifacts/ChangeLog gnv-artifacts/src/main/java/de/intevation/gnv/artifacts/context/GNVArtifactContext.java gnv-artifacts/src/main/java/de/intevation/gnv/artifacts/context/GNVArtifactContextFactory.java gnv-artifacts/src/main/java/de/intevation/gnv/math/AreaInterpolation.java gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation2D.java gnv-artifacts/src/main/java/de/intevation/gnv/math/L1Comparator.java gnv-artifacts/src/main/java/de/intevation/gnv/state/profile/horizontalcrosssection/HorizontalCrossSectionMeshOutputState.java gnv-artifacts/src/main/java/de/intevation/gnv/state/profile/verticalcrosssection/VerticalCrossSectionOutputState.java gnv-artifacts/src/test/java/de/intevation/gnv/artifacts/util/ShapeFileWriterTestCase.java
diffstat 9 files changed, 274 insertions(+), 33 deletions(-) [+]
line wrap: on
line diff
--- a/gnv-artifacts/ChangeLog	Wed Jan 13 23:10:56 2010 +0000
+++ b/gnv-artifacts/ChangeLog	Thu Jan 14 10:34:05 2010 +0000
@@ -1,3 +1,28 @@
+2010-01-13	Sascha L. Teichmann	<sascha.teichmann@intevation.de>
+
+	* src/test/java/de/intevation/gnv/artifacts/util/ShapeFileWriterTestCase.java,
+	  src/main/java/de/intevation/gnv/state/profile/verticalcrosssection/VerticalCrossSectionOutputState.java:
+	  Removed needless imports.
+
+	* src/main/java/de/intevation/gnv/math/AreaInterpolation.java: New. Interpolates
+	  area for a given bounding box, taking gaps and DEM into account.
+	  Not very fast. Use bilinear interpolation to match the "Profilschnitt".
+	  Possible TODOs: 
+	  - speed up by assuming the grid is not sparse.
+	  - use higher interpolation methods.
+
+	* src/main/java/de/intevation/gnv/artifacts/context/GNVArtifactContext.java,
+	  src/main/java/de/intevation/gnv/artifacts/context/GNVArtifactContextFactory.java:
+	  Added configuration for ground interpolation.
+
+	* src/main/java/de/intevation/gnv/state/profile/horizontalcrosssection/HorizontalCrossSectionMeshOutputState.java:
+	  Add helper functions to access configuration.
+
+	* src/main/java/de/intevation/gnv/math/Interpolation2D.java: Simplified Code.
+
+	* src/main/java/de/intevation/gnv/math/L1Comparator.java: add setReference()
+	  method.
+
 2010-01-13  Ingo Weinzierl <ingo_weinzierl@web.de>
 
 	* src/main/java/de/intevation/gnv/artifacts/GNVArtifactBase.java: 'advance'
--- a/gnv-artifacts/src/main/java/de/intevation/gnv/artifacts/context/GNVArtifactContext.java	Wed Jan 13 23:10:56 2010 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/artifacts/context/GNVArtifactContext.java	Thu Jan 14 10:34:05 2010 +0000
@@ -46,7 +46,13 @@
         DEFAULT_HORIZONTAL_CROSS_SECTION_PROFILE_SAMPLES = Integer.valueOf(250);
     
     public static final Integer 
-        DEFAULT_HORIZONTAL_CROSS_SECTION_SAMPLES = Integer.valueOf(250);
+        DEFAULT_HORIZONTAL_CROSS_SECTION_SAMPLES = Integer.valueOf(1024);
+
+    public static final String HORIZONTAL_CROSS_SECTION_GROUND_INTERPOLATION_KEY =
+        "gnv.horizontal.cross.section.ground.interpolation";
+
+    public static final String DEFAULT_HORIZONTAL_CROSS_SECTION_GROUND_INTERPOLATION =
+        "bilinear";
 
     public static final String PALETTES_KEY =
         "gnv.color.palettes";
--- a/gnv-artifacts/src/main/java/de/intevation/gnv/artifacts/context/GNVArtifactContextFactory.java	Wed Jan 13 23:10:56 2010 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/artifacts/context/GNVArtifactContextFactory.java	Thu Jan 14 10:34:05 2010 +0000
@@ -72,6 +72,9 @@
     public final static String HORIZONTAL_CROSS_SECTION_PROFILE_SAMPLES =
         "/artifact-database/gnv/horizontal-cross-section-profile/samples/@number";
 
+    public final static String HORIZONTAL_CROSS_SECTION_GROUND_INTERPOLATION =
+        "/artifact-database/gnv/horizontal-cross-section/ground/@interpolation";
+
     public final static String HORIZONTAL_CROSS_SECTION_SAMPLES =
         "/artifact-database/gnv/horizontal-cross-section/samples/@number";
     
@@ -323,6 +326,33 @@
 
         configureHorizontalCrossSectionSamples(config, context);
         configureHorizontalCrossSectionResultShapeFilePath(config, context);
+        configureHorizontalCrossSectionGroundInterpolation(config, context);
+    }
+
+    protected void configureHorizontalCrossSectionGroundInterpolation(
+        Document           config,
+        GNVArtifactContext context
+    )
+    {
+        log.info(
+            "configuration of horizontal cross section ground interpolation");
+
+        String interpolation = Config.getStringXPath(
+            config,
+            HORIZONTAL_CROSS_SECTION_GROUND_INTERPOLATION);
+
+        if (interpolation == null 
+        || (interpolation = interpolation.trim()).length() == 0) {
+            interpolation = GNVArtifactContext
+                .DEFAULT_HORIZONTAL_CROSS_SECTION_GROUND_INTERPOLATION;
+        }
+
+        log.info("ground interpolation: " + interpolation);
+
+        context.put(
+            GNVArtifactContext
+                .HORIZONTAL_CROSS_SECTION_GROUND_INTERPOLATION_KEY,
+            interpolation);
     }
 
     protected void configureHorizontalCrossSectionResultShapeFilePath(
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/AreaInterpolation.java	Thu Jan 14 10:34:05 2010 +0000
@@ -0,0 +1,156 @@
+package de.intevation.gnv.math;
+
+import java.io.Serializable;
+
+import java.util.ArrayList;
+import java.util.Arrays;
+import java.util.List;
+import java.util.HashMap;
+import java.util.Collections;
+
+import java.awt.Dimension;
+
+import com.vividsolutions.jts.geom.Coordinate;
+import com.vividsolutions.jts.geom.Envelope;
+
+import com.vividsolutions.jts.index.quadtree.Quadtree;
+
+import org.apache.log4j.Logger;
+
+/**
+ *  @author Sascha L. Teichmann (sascha.teichmann@intevation.de)
+ */
+public class AreaInterpolation
+implements   Serializable
+{
+    private static Logger log = Logger.getLogger(AreaInterpolation.class);
+
+    protected double [] raster;
+
+    protected int width;
+    protected int height;
+
+    public AreaInterpolation() {
+    }
+
+    public int getWidth() {
+        return width;
+    }
+
+    public int getHeight() {
+        return height;
+    }
+
+    public double [] getRaster() {
+        return raster;
+    }
+
+    public boolean interpolate(
+        List<? extends Point2d> points,
+        Envelope                boundingBox,
+        Dimension               samples,
+        XYDepth                 depth
+    ) {
+        if (points == null || points.isEmpty()) {
+            return false;
+        }
+
+        double [] buffer = Interpolation2D.calculateBuffer(points);
+        double dxMax = buffer[0];
+        double dyMax = buffer[1];
+
+        Quadtree spatialIndex = new Quadtree();
+
+        for (int i = points.size()-1; i >= 0; --i) {
+            Point2d p = points.get(i);
+            spatialIndex.insert(p.envelope(), p);
+        }
+
+        int W = samples.width;
+        int H = samples.height;
+
+        double cellWidth  = boundingBox.getWidth()  / W;
+        double cellHeight = boundingBox.getHeight() / H;
+
+        Envelope     queryBuffer = new Envelope();
+        Point2d []   neighbors   = new Point2d[4];
+        Coordinate   center      = new Coordinate();
+        L1Comparator invL1       = new L1Comparator(center);
+
+        double [] raster = new double[W*H];
+        Arrays.fill(raster, Double.NaN);
+
+        int row = 0;
+        for (int j = 0; j < H; ++j, row += W) {
+            double y = j*cellHeight + 0.5d*cellHeight;
+            double x = 0.5d*cellWidth;
+            for (int i = 0; i < W; ++i, x += cellWidth) {
+                queryBuffer.init(
+                    x - dxMax, x + dxMax, 
+                    y - dyMax, y + dyMax);
+
+                List potential = spatialIndex.query(queryBuffer);
+
+                if (potential.isEmpty()) {
+                    continue;
+                }
+
+                center.x = x; center.y = y;
+                Collections.sort(potential, invL1);
+
+                neighbors[0] = neighbors[1] = 
+                neighbors[2] = neighbors[3] = null;
+
+                int mask = 0;
+                // reversed order is essential here!
+                for (int k = potential.size()-1; k >= 0; --k) {
+                    Point2d n = (Point2d)potential.get(k);
+                    int code = n.x > center.x ? 1 : 0;
+                    if (n.y > center.y) code |= 2;
+                    neighbors[code] = n;
+                    mask |= 1 << code;
+                }
+
+                int numNeighbors = Integer.bitCount(mask);
+
+                // only interpolate if we have all four neighbors
+                // we do not have any gaps and we are below surface.
+                if (numNeighbors == 4
+                && !neighbors[0].hasIGap(neighbors[1])
+                && !neighbors[1].hasJGap(neighbors[3])
+                && !neighbors[3].hasIGap(neighbors[2])
+                && !neighbors[2].hasJGap(neighbors[0])
+                && depth.depth(center) <= 0d
+                ) {
+                    double z1 = Interpolation2D.interpolate(
+                        neighbors[0].x, neighbors[0].z,
+                        neighbors[1].x, neighbors[1].z,
+                        center.x);
+                    double z2 = Interpolation2D.interpolate(
+                        neighbors[2].x, neighbors[2].z,
+                        neighbors[3].x, neighbors[3].z,
+                        center.x);
+                    double y1 = Interpolation2D.interpolate(
+                        neighbors[0].x, neighbors[0].y,
+                        neighbors[1].x, neighbors[1].y,
+                        center.x);
+                    double y2 = Interpolation2D.interpolate(
+                        neighbors[2].x, neighbors[2].y,
+                        neighbors[3].x, neighbors[3].y,
+                        center.x);
+                    raster[row+i] = Interpolation2D.interpolate(
+                        y1, z1,
+                        y2, z2,
+                        center.y);
+                }
+            }
+        }
+
+        this.raster = raster;
+        this.width  = W;
+        this.height = H;
+
+        return true;
+    }
+}
+// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
--- a/gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation2D.java	Wed Jan 13 23:10:56 2010 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation2D.java	Thu Jan 14 10:34:05 2010 +0000
@@ -135,6 +135,8 @@
         int missedInterpolations = 0;
         int interpolations = 0;
 
+        L1Comparator invL1 = new L1Comparator(center);
+
         for (double p = from; p <= to; p += dP) {
             if (!linearToMap.locate(p, center)) {
                 continue;
@@ -145,7 +147,6 @@
 
             List potential = spatialIndex.query(queryBuffer);
 
-            L1Comparator invL1 = new L1Comparator(center);
             Collections.sort(potential, invL1);
 
             neighbors[0] = neighbors[1] = 
--- a/gnv-artifacts/src/main/java/de/intevation/gnv/math/L1Comparator.java	Wed Jan 13 23:10:56 2010 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/L1Comparator.java	Thu Jan 14 10:34:05 2010 +0000
@@ -12,10 +12,17 @@
 {
     private Coordinate ref;
 
+    public L1Comparator() {
+    }
+
     public L1Comparator(Coordinate ref) {
         this.ref = ref;
     }
 
+    public void setReference(Coordinate ref) {
+        this.ref = ref;
+    }
+
     public int compare(Object a, Object b) {
         Coordinate pa = (Coordinate)a;
         Coordinate pb = (Coordinate)b;
--- a/gnv-artifacts/src/main/java/de/intevation/gnv/state/profile/horizontalcrosssection/HorizontalCrossSectionMeshOutputState.java	Wed Jan 13 23:10:56 2010 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/state/profile/horizontalcrosssection/HorizontalCrossSectionMeshOutputState.java	Thu Jan 14 10:34:05 2010 +0000
@@ -3,10 +3,12 @@
  */
 package de.intevation.gnv.state.profile.horizontalcrosssection;
 
+import java.io.File;
 import java.io.IOException;
 import java.io.OutputStream;
 import java.io.OutputStreamWriter;
 import java.io.UnsupportedEncodingException;
+
 import java.util.Arrays;
 import java.util.Collection;
 import java.util.Iterator;
@@ -16,42 +18,53 @@
 import org.apache.log4j.Logger;
 
 import org.jfree.chart.ChartTheme;
+
 import org.w3c.dom.Node;
 
 import com.vividsolutions.jts.geom.Coordinate;
 import com.vividsolutions.jts.geom.Envelope;
 import com.vividsolutions.jts.geom.Polygon;
-import com.vividsolutions.jts.io.ParseException;
-import com.vividsolutions.jts.io.WKTReader;
 
 import au.com.bytecode.opencsv.CSVWriter;
+
 import de.intevation.gnv.artifacts.cache.CacheFactory;
 import de.intevation.gnv.artifacts.context.GNVArtifactContext;
+
 import de.intevation.gnv.chart.Chart;
 import de.intevation.gnv.chart.ChartLabels;
+
 import de.intevation.gnv.geobackend.base.Result;
+
 import de.intevation.gnv.geobackend.base.query.QueryExecutor;
 import de.intevation.gnv.geobackend.base.query.QueryExecutorFactory;
 import de.intevation.gnv.geobackend.base.query.exception.QueryException;
+
 import de.intevation.gnv.state.InputData;
+
 import de.intevation.gnv.state.exception.StateException;
+
 import de.intevation.gnv.state.timeseries.TimeSeriesOutputState;
+
 import de.intevation.gnv.statistics.Statistics;
+
 import de.intevation.gnv.utils.StringUtils;
 import de.intevation.gnv.utils.WKTUtils;
 
 import de.intevation.artifactdatabase.Config;
+
 import de.intevation.artifacts.CallContext;
 
+import de.intevation.gnv.geobackend.sde.datasources.RasterObject;
+
 /**
- * @author Tim Englich <tim.englich@intevation.de>
- * 
+ * @author Tim Englich         (tim.englich@intevation.de)
+ * @author Sascha L. Teichmann (sascha.teichmann@intevation.de)
  */
 public class HorizontalCrossSectionMeshOutputState 
-                                      extends TimeSeriesOutputState {
-
+extends      TimeSeriesOutputState
+{
     private static Logger log = Logger
-    .getLogger(HorizontalCrossSectionMeshOutputState.class);
+        .getLogger(HorizontalCrossSectionMeshOutputState.class);
     
     /**
      * The UID of this Class
@@ -95,23 +108,6 @@
         log.info("Chart not in cache yet.");
         
         log.warn("This sort of chart is not implemented yet.");
-        /* TODO Implement a special chart for this sort of charts.
-        chart = new HorizontalProfileChart(
-            chartLables,
-            chartTheme,
-            parameters,
-            measurements,
-            result,
-            dates,
-            locale
-        );
-        chart.generateChart();
-
-        if (CACHE_CHART) {
-            log.info("Put chart into cache.");
-            purifyChart(chart, uuid);
-        }
-        */
         
         InputData meshPolygon = inputData.get("mesh_polygon");
         String meshPolygonWkt = null;
@@ -226,11 +222,12 @@
      *      java.util.Collection)
      */
     @Override
-    protected void createCSV(OutputStream outputStream,
-                             Collection<Result> chartResult)
-                                            throws UnsupportedEncodingException,
-                                                   IOException,
-                                                   StateException {
+    protected void createCSV(
+        OutputStream       outputStream,
+        Collection<Result> chartResult
+    ) throws UnsupportedEncodingException, IOException, StateException 
+    {
+        /*
         if (chartResult != null) {
             try {
                 CSVWriter writer = new CSVWriter(new OutputStreamWriter(
@@ -260,6 +257,7 @@
             throw new StateException(
                     "No Data given for generating an CSV-File.");
         }
+        */
     }
     
     /**
@@ -282,4 +280,24 @@
             : GNVArtifactContext.DEFAULT_HORIZONTAL_CROSS_SECTION_SAMPLES;
     }
 
+    private static File shapefileDirectory(CallContext callContext) {
+        GNVArtifactContext context =
+            (GNVArtifactContext)callContext.globalContext();
+        File dir = (File)context.get(
+            GNVArtifactContext.HORIZONTAL_CROSS_SECTION_RESULT_SHAPEFILE_PATH_KEY);
+        return dir != null
+            ? dir
+            : GNVArtifactContext.DEFAULT_HORIZONTAL_CROSS_SECTION_PROFILE_SHAPEFILE_PATH;
+    }
+
+    private static int getGroundInterpolation(CallContext callContext) {
+        GNVArtifactContext context = 
+            (GNVArtifactContext)callContext.globalContext();
+
+        String interpolation = (String)context.get(
+            GNVArtifactContext.HORIZONTAL_CROSS_SECTION_GROUND_INTERPOLATION_KEY);
+
+        return RasterObject.getInterpolationType(interpolation);
+    }
 }
+// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
--- a/gnv-artifacts/src/main/java/de/intevation/gnv/state/profile/verticalcrosssection/VerticalCrossSectionOutputState.java	Wed Jan 13 23:10:56 2010 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/state/profile/verticalcrosssection/VerticalCrossSectionOutputState.java	Thu Jan 14 10:34:05 2010 +0000
@@ -57,7 +57,6 @@
 import de.intevation.gnv.math.IJKey;
 import de.intevation.gnv.math.LinearMetrics;
 import de.intevation.gnv.math.Interpolation3D;
-import de.intevation.gnv.math.ConstantXYDepth;
 
 import de.intevation.gnv.state.InputData;
 
--- a/gnv-artifacts/src/test/java/de/intevation/gnv/artifacts/util/ShapeFileWriterTestCase.java	Wed Jan 13 23:10:56 2010 +0000
+++ b/gnv-artifacts/src/test/java/de/intevation/gnv/artifacts/util/ShapeFileWriterTestCase.java	Thu Jan 14 10:34:05 2010 +0000
@@ -13,7 +13,6 @@
 import org.apache.log4j.Logger;
 
 import com.vividsolutions.jts.geom.MultiLineString;
-import com.vividsolutions.jts.io.ParseException;
 import com.vividsolutions.jts.io.WKTReader;
 
 import de.intevation.gnv.utils.Pair;

http://dive4elements.wald.intevation.org