# HG changeset patch # User Sascha L. Teichmann # Date 1264974615 0 # Node ID b248531fa20ba83b18bb016820ccac04bec218c5 # Parent 3c939c95c477e1916b3ce652f0b7f563734840be Added experimental support for extrapolation in "Horizontalschnitte" gnv-artifacts/trunk@648 c6561f87-3c4e-4783-a992-168aeb5c3f6f diff -r 3c939c95c477 -r b248531fa20b gnv-artifacts/ChangeLog --- a/gnv-artifacts/ChangeLog Fri Jan 29 13:28:28 2010 +0000 +++ b/gnv-artifacts/ChangeLog Sun Jan 31 21:50:15 2010 +0000 @@ -1,3 +1,37 @@ +2010-01-31 Sascha L. Teichmann + + * src/main/java/de/intevation/gnv/artifacts/context/GNVArtifactContext.java, + src/main/java/de/intevation/gnv/artifacts/context/GNVArtifactContextFactory.java, + doc/conf/conf.xml: Added configuration for extrapolation in "Horizontalschnitte". + + Use gnv/horizontal-cross-section/extrapolation/@rounds with + integer rounds > 0 to turn extrapolation on (default: 0). + Rounds is a number of successive point extrapolations which means that + the grid is successively filled with missing points based on prior rounds. + The larger 'rounds' get more gaps are filled synthetic generated points. + + Set this to 2 to get good results for the model data FIS. + + * src/main/java/de/intevation/gnv/state/profile/horizontalcrosssection/HorizontalCrossSectionMeshOutputState.java: + Foward configuration to area interpolation. + + * src/main/java/de/intevation/gnv/math/AreaInterpolation.java: + Foward configuration to GridCell. + + * src/main/java/de/intevation/gnv/math/Point2d.java: Added method to extrapolate + point along a line spanned by two points. + Calculate Inverse Distance Weighting (IDW) for a given set of points on + z components. Added method to check if set of points are near a given + point. + + * src/main/java/de/intevation/gnv/math/GridCell.java: Before building the + i/j cells the grid is filled with synthetic generated points. The + position is estimated from the neighboring points. The parameter values + are calculated by IDW. Some care is taken to avoid invalid grid topologies. + + TODO: Implement this for the "Profillschnitt" too to keep the inner + symmetry. + 2010-01-29 Tim Englich * src/test/java/de/intevation/gnv/artifacts/MeshVerticalProfileTestCase.java (testArtifact), @@ -252,7 +286,7 @@ * src/main/java/de/intevation/gnv/chart/VerticalCrossSectionChart.java: Handle the seabad polygon color. -2010-01-24 Sascha L. Teichmann fill-color +2010-01-24 Sascha L. Teichmann * contrib/palette2qgis.xsl: Cosmetic cleanups. diff -r 3c939c95c477 -r b248531fa20b gnv-artifacts/doc/conf/conf.xml --- a/gnv-artifacts/doc/conf/conf.xml Fri Jan 29 13:28:28 2010 +0000 +++ b/gnv-artifacts/doc/conf/conf.xml Sun Jan 31 21:50:15 2010 +0000 @@ -421,6 +421,7 @@ +

diff -r 3c939c95c477 -r b248531fa20b gnv-artifacts/src/main/java/de/intevation/gnv/artifacts/context/GNVArtifactContext.java --- a/gnv-artifacts/src/main/java/de/intevation/gnv/artifacts/context/GNVArtifactContext.java Fri Jan 29 13:28:28 2010 +0000 +++ b/gnv-artifacts/src/main/java/de/intevation/gnv/artifacts/context/GNVArtifactContext.java Sun Jan 31 21:50:15 2010 +0000 @@ -36,6 +36,9 @@ public static final String HORIZONTAL_CROSS_SECTION_SAMPLES_KEY = "gnv.horizontal.cross.section.samples"; + public static final String HORIZONTAL_CROSS_SECTION_EXTRAPOLATION_ROUNDS_KEY = + "gnv.horizontal.cross.section.extrapolation.rounds"; + public static final String HORIZONTAL_CROSS_SECTION_RESULT_SHAPEFILE_PATH_KEY = "gnv.horizontal.cross.section.result.shapefile"; @@ -50,6 +53,9 @@ public static final Integer DEFAULT_HORIZONTAL_CROSS_SECTION_SAMPLES = Integer.valueOf(1024); + public static final Integer + DEFAULT_HORIZONTAL_CROSS_SECTION_EXTRAPOLATION_ROUNDS = Integer.valueOf(0); + public static final String HORIZONTAL_CROSS_SECTION_GROUND_INTERPOLATION_KEY = "gnv.horizontal.cross.section.ground.interpolation"; diff -r 3c939c95c477 -r b248531fa20b gnv-artifacts/src/main/java/de/intevation/gnv/artifacts/context/GNVArtifactContextFactory.java --- a/gnv-artifacts/src/main/java/de/intevation/gnv/artifacts/context/GNVArtifactContextFactory.java Fri Jan 29 13:28:28 2010 +0000 +++ b/gnv-artifacts/src/main/java/de/intevation/gnv/artifacts/context/GNVArtifactContextFactory.java Sun Jan 31 21:50:15 2010 +0000 @@ -1,6 +1,3 @@ -/** - * - */ package de.intevation.gnv.artifacts.context; import java.awt.Color; @@ -81,6 +78,9 @@ public final static String HORIZONTAL_CROSS_SECTION_SAMPLES = "/artifact-database/gnv/horizontal-cross-section/samples/@number"; + public final static String HORIZONTAL_CROSS_SECTION_EXTRAPOLATION_ROUNDS = + "/artifact-database/gnv/horizontal-cross-section/extrapolation/@rounds"; + public final static String HORIZONTAL_CROSS_SECTION_RESULT_SHAPEFILE_PATH = "/artifact-database/gnv/horizontal-cross-section/result-shapefile-directory/@path"; @@ -364,10 +364,47 @@ log.info("configuration of horizontal cross section"); configureHorizontalCrossSectionSamples(config, context); + configureHorizontalCrossSectionExtrapolation(config, context); configureHorizontalCrossSectionResultShapeFilePath(config, context); configureHorizontalCrossSectionGroundInterpolation(config, context); } + protected void configureHorizontalCrossSectionExtrapolation( + Document config, + GNVArtifactContext context + ) + { + log.info( + "configuration of horizontal cross section extrapolation"); + + String extrapolationRoundsValue = Config.getStringXPath( + config, + HORIZONTAL_CROSS_SECTION_EXTRAPOLATION_ROUNDS); + + Integer extrapolationRounds = + GNVArtifactContext.DEFAULT_HORIZONTAL_CROSS_SECTION_EXTRAPOLATION_ROUNDS; + + if (extrapolationRoundsValue != null + && (extrapolationRoundsValue = extrapolationRoundsValue.trim()).length() > 0 + ) { + try { + extrapolationRounds = + Integer.valueOf(extrapolationRoundsValue); + } + catch (NumberFormatException nfe) { + log.error( + "'" + extrapolationRoundsValue + "' is not a valid integer"); + } + } + + log.info("extrapolation rounds: " + extrapolationRounds); + + context.put( + GNVArtifactContext + .HORIZONTAL_CROSS_SECTION_EXTRAPOLATION_ROUNDS_KEY, + extrapolationRounds); + } + protected void configureHorizontalCrossSectionGroundInterpolation( Document config, GNVArtifactContext context diff -r 3c939c95c477 -r b248531fa20b gnv-artifacts/src/main/java/de/intevation/gnv/math/AreaInterpolation.java --- a/gnv-artifacts/src/main/java/de/intevation/gnv/math/AreaInterpolation.java Fri Jan 29 13:28:28 2010 +0000 +++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/AreaInterpolation.java Sun Jan 31 21:50:15 2010 +0000 @@ -49,7 +49,8 @@ List points, Envelope boundingBox, Dimension samples, - XYDepth depth + XYDepth depth, + int extrapolationRounds ) { boolean debug = log.isDebugEnabled(); @@ -62,7 +63,8 @@ points, Interpolation2D.relevantArea( boundingBox, - points)); + points), + extrapolationRounds); if (cells.isEmpty()) { log.warn("no cells to interpolate"); diff -r 3c939c95c477 -r b248531fa20b gnv-artifacts/src/main/java/de/intevation/gnv/math/GridCell.java --- a/gnv-artifacts/src/main/java/de/intevation/gnv/math/GridCell.java Fri Jan 29 13:28:28 2010 +0000 +++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/GridCell.java Sun Jan 31 21:50:15 2010 +0000 @@ -1,5 +1,7 @@ package de.intevation.gnv.math; +import com.vividsolutions.jts.algorithm.CGAlgorithms; + import com.vividsolutions.jts.geom.Coordinate; import com.vividsolutions.jts.geom.Envelope; import com.vividsolutions.jts.geom.Geometry; @@ -73,8 +75,21 @@ } protected void createPolygon() { + Coordinate [] coords = new Coordinate [] { p1, p2, p3, p4, p1 }; + if (!CGAlgorithms.isCCW(coords)) { + for (int i = 0, j = coords.length-1; i < j; ++i, --j) { + Coordinate c = coords[i]; + coords[i] = coords[j]; + coords[j] = c; + } + } LinearRing shell = GEOMETRY_FACTORY - .createLinearRing(new Coordinate [] { p1, p2, p3, p4, p1 }); + .createLinearRing(coords); + + if (!shell.isValid()) { + log.warn("linear ring is not valid"); + } + polygon = GEOMETRY_FACTORY.createPolygon(shell, null); } @@ -96,6 +111,111 @@ List points, Envelope relevantArea ) { + return pointsToGridCells(points, relevantArea, 0); + } + + private static final int NEIGHBORS [][] = { + { -1, -1 }, // 0 + { -1, 0 }, // 1 + { -1, +1 }, // 2 + { 0, +1 }, // 3 + { +1, +1 }, // 4 + { +1, 0 }, // 5 + { +1, -1 }, // 6 + { 0, -1 } // 7 + }; + + public static int extrapolate( + HashMap> rows, + int minI, int maxI, + int minJ, int maxJ, + int rounds, + Envelope relevantArea + ) { + Point2d [] neighbors = new Point2d[NEIGHBORS.length]; + Point2d [] positions = new Point2d[NEIGHBORS.length]; + + int total = 0; + + ArrayList> prio = + new ArrayList>(NEIGHBORS.length); + + for (int i = 0; i < NEIGHBORS.length; ++i) { + prio.add(new ArrayList()); + } + + while (rounds-- > 0) { + for (int i = minI; i <= maxI; ++i) { + for (int j = minJ; j <= maxJ; ++j) { + Point2d p = get(rows, i, j); + if (p != null) { + continue; + } + + int count = 0; + + for (int k = 0; k < neighbors.length; ++k) { + neighbors[k] = positions[k] = null; + int dij [] = NEIGHBORS[k]; + Point2d n1 = get(rows, i+ dij[0], j+ dij[1]); + Point2d n2 = get(rows, i+2*dij[0], j+2*dij[1]); + if (n1 != null && n2 != null) { + ++count; + } + } + + if (count > 0) { + prio.get(count-1).add(new IJKey(i, j)); + } + } // for all columns + } // for all rows + + int N = 0; + + for (int l = NEIGHBORS.length-1; l >= 0; --l) { + ArrayList list = prio.get(l); + for (IJKey ij: list) { + int i = ij.i; + int j = ij.j; + for (int k = 0; k < neighbors.length; ++k) { + neighbors[k] = positions[k] = null; + int dij [] = NEIGHBORS[k]; + Point2d n1 = get(rows, i+ dij[0], j+ dij[1]); + Point2d n2 = get(rows, i+2*dij[0], j+2*dij[1]); + if (n1 != null && n2 != null) { + neighbors[k] = n1; + positions[k] = n1.extrapolate(-1d, n2); + } + } + + Point2d avg = Point2d.average(positions); + + if (avg != null && avg.near(positions) + && (relevantArea == null || relevantArea.contains(avg.x, avg.y))) { + avg.i = i; + avg.j = j; + avg.inverseDistanceWeighting(neighbors); + put(rows, avg, i, j); + } + } + N += list.size(); + list.clear(); + } + + if (N == 0) { + break; + } + total += N; + } // for all rounds + + return total; + } + + public static List pointsToGridCells( + List points, + Envelope relevantArea, + int extrapolationRounds + ) { int minI = Integer.MAX_VALUE; int maxI = Integer.MIN_VALUE; int minJ = Integer.MAX_VALUE; @@ -129,25 +249,36 @@ ArrayList cells = new ArrayList(points.size()); + int extrapolated = extrapolate( + rows, + minI, maxI, + minJ, maxJ, + extrapolationRounds, + relevantArea); + for (int i = minI + 1; i <= maxI; ++i) { HashMap row1 = rows.get(i-1); HashMap row2 = rows.get(i); - if (row1 != null && row2 != null) { - for (int j = minJ + 1; j < maxJ; ++j) { - Point2d p1 = row1.get(j-1); - Point2d p2 = row1.get(j); - Point2d p3 = row2.get(j); - Point2d p4 = row2.get(j-1); - if (p1 != null && p2 != null && p3 != null && p4 != null) { - cells.add(new GridCell(p1, p2, p3, p4)); - } - } // for all columns in row + if (row1 == null || row2 == null) { + continue; + } + + for (int j = minJ + 1; j < maxJ; ++j) { + Point2d p1 = row1.get(j-1); + Point2d p2 = row1.get(j); + Point2d p3 = row2.get(j); + Point2d p4 = row2.get(j-1); + + if (p1 != null && p2 != null && p3 != null && p4 != null) { + cells.add(new GridCell(p1, p2, p3, p4)); + } } } // for all rows if (log.isDebugEnabled()) { log.debug("culled points: " + culled); + log.debug("extrapolated points: " + extrapolated); log.debug("min/max i: " + minI + " / " + maxI); log.debug("min/max j: " + minJ + " / " + maxJ); log.debug("cells found: " + cells.size()); @@ -155,5 +286,26 @@ return cells; } + + private static Point2d get( + HashMap> rows, + int i, int j + ) { + HashMap row = rows.get(i); + return row != null ? row.get(j) : null; + } + + private static void put( + HashMap> rows, + Point2d point, + int i, int j + ) { + Integer I = Integer.valueOf(i); + HashMap row = rows.get(I); + if (row == null) { + rows.put(I, row = new HashMap()); + } + row.put(j, point); + } } // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 : diff -r 3c939c95c477 -r b248531fa20b gnv-artifacts/src/main/java/de/intevation/gnv/math/Point2d.java --- a/gnv-artifacts/src/main/java/de/intevation/gnv/math/Point2d.java Fri Jan 29 13:28:28 2010 +0000 +++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/Point2d.java Sun Jan 31 21:50:15 2010 +0000 @@ -1,16 +1,20 @@ package de.intevation.gnv.math; +import com.vividsolutions.jts.geom.Coordinate; +import com.vividsolutions.jts.geom.Envelope; + import java.util.Comparator; -import com.vividsolutions.jts.geom.Envelope; -import com.vividsolutions.jts.geom.Coordinate; +import org.apache.log4j.Logger; /** - * @author Sascha L. Teichmann + * @author Sascha L. Teichmann (sascha.teichmann@intevation.de) */ public class Point2d extends Coordinate { + private static Logger log = Logger.getLogger(Point2d.class); + public static final double EPSILON = 1e-3d; public static final Comparator X_COMPARATOR = new Comparator() { @@ -59,6 +63,10 @@ public Point2d() { } + public Point2d(double x, double y) { + super(x, y); + } + public Point2d(double x, double y, int i, int j) { super(x, y); this.i = i; @@ -97,5 +105,95 @@ public boolean hasJGap(Point2d other) { return Math.abs(j - other.j) > 1; } + + public Point2d extrapolate(double t, Point2d other) { + if (other == null) { + return null; + } + Point2d p = newPoint(); + p.x = t*(other.x - x) + x; + p.y = t*(other.y - y) + y; + return p; + } + + public Point2d newPoint() { + return new Point2d(0d, 0d); + } + + public Point2d newPoint(double x, double y) { + return new Point2d(x, y); + } + + public void inverseDistanceWeighting(Point2d [] ps) { + + double sum = 0d; + + double [] d = new double[ps.length]; + + for (int i = 0; i < ps.length; ++i) { + Point2d p = ps[i]; + if (p != null) { + double di = distance(p); + if (di < 1e-5d) { z = p.z; return; } + di = 1d/di; + d[i] = di; + sum += di; + } + } + + if (sum == 0d) { + return; + } + + double v = 0d; + + for (int i = 0; i < ps.length; ++i) { + Point2d p = ps[i]; + if (p != null) { + v += p.z*d[i]; + } + } + z = v/sum; + } + + public static Point2d average(Point2d [] ps) { + + Point2d p = null; + int count = 0; + + for (int i = 0; i < ps.length; ++i) { + Point2d t = ps[i]; + if (t != null) { + ++count; + if (p == null) { + p = t.newPoint(t.x, t.y); + } + else { + p.x += t.x; + p.y += t.y; + } + } + } + + if (p != null) { + double s = 1d/count; + p.x *= s; + p.y *= s; + } + + return p; + } + + public boolean near(Point2d [] ps) { + + for (int i = 0; i < ps.length; ++i) { + Point2d p = ps[i]; + if (p != null && distance(p) > EPSILON) { + return false; + } + } + + return true; + } } // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8: diff -r 3c939c95c477 -r b248531fa20b gnv-artifacts/src/main/java/de/intevation/gnv/state/profile/horizontalcrosssection/HorizontalCrossSectionMeshOutputState.java --- a/gnv-artifacts/src/main/java/de/intevation/gnv/state/profile/horizontalcrosssection/HorizontalCrossSectionMeshOutputState.java Fri Jan 29 13:28:28 2010 +0000 +++ b/gnv-artifacts/src/main/java/de/intevation/gnv/state/profile/horizontalcrosssection/HorizontalCrossSectionMeshOutputState.java Sun Jan 31 21:50:15 2010 +0000 @@ -522,12 +522,14 @@ int numSamples = numSamples(callContext); int groundInterpolation = getGroundInterpolation(callContext); + int extrapolationRounds = extrapolationRounds(callContext); if (!interpolation.interpolate( input.getPoints(), boundingBox, new Dimension(numSamples, numSamples), - new QueriedXYDepth(groundInterpolation) + new QueriedXYDepth(groundInterpolation), + extrapolationRounds )) { log.error("interpolation failed"); return null; @@ -636,6 +638,16 @@ : GNVArtifactContext.DEFAULT_HORIZONTAL_CROSS_SECTION_SAMPLES; } + private static int extrapolationRounds(CallContext callContext) { + GNVArtifactContext context = + (GNVArtifactContext)callContext.globalContext(); + Integer extrapolationRounds = (Integer)context.get( + GNVArtifactContext.HORIZONTAL_CROSS_SECTION_EXTRAPOLATION_ROUNDS_KEY); + return extrapolationRounds != null + ? extrapolationRounds.intValue() + : GNVArtifactContext.DEFAULT_HORIZONTAL_CROSS_SECTION_EXTRAPOLATION_ROUNDS; + } + private static File shapefileDirectory(CallContext callContext) { GNVArtifactContext context = (GNVArtifactContext)callContext.globalContext();