Mercurial > dive4elements > gnv-client
view gnv-artifacts/src/main/java/de/intevation/gnv/math/AreaInterpolation.java @ 592:3c939c95c477
Removed accidental checkin of Test-Shapefiles
gnv-artifacts/trunk@647 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author | Tim Englich <tim.englich@intevation.de> |
---|---|
date | Fri, 29 Jan 2010 13:28:28 +0000 |
parents | 44415ae01ddb |
children | b248531fa20b |
line wrap: on
line source
package de.intevation.gnv.math; import com.vividsolutions.jts.geom.Coordinate; import com.vividsolutions.jts.geom.Envelope; import com.vividsolutions.jts.index.strtree.STRtree; import java.awt.Dimension; import java.io.Serializable; import java.util.Arrays; import java.util.List; 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 static final double EPS = 1e-6d; public boolean interpolate( List<? extends Point2d> points, Envelope boundingBox, Dimension samples, XYDepth depth ) { boolean debug = log.isDebugEnabled(); if (points == null || points.isEmpty()) { log.warn("no points to interpolate"); return false; } List<GridCell> cells = GridCell.pointsToGridCells( points, Interpolation2D.relevantArea( boundingBox, points)); if (cells.isEmpty()) { log.warn("no cells to interpolate"); return false; } int W = samples.width; int H = samples.height; double cellWidth = boundingBox.getWidth() / W; double cellHeight = boundingBox.getHeight() / H; STRtree spatialIndex = new STRtree(); for (GridCell cell: cells) { spatialIndex.insert(cell.getEnvelope(), cell); } if (debug) { log.debug("width: " + boundingBox.getWidth()); log.debug("height: " + boundingBox.getHeight()); log.debug("sample width: " + W); log.debug("sample height: " + H); log.debug("cell width: " + cellWidth); log.debug("cell height: " + cellHeight); } Envelope queryBuffer = new Envelope(); Coordinate center = new Coordinate(); GridCell.CellFinder finder = new GridCell.CellFinder(); double [] raster = new double[W*H]; Arrays.fill(raster, Double.NaN); double minX = boundingBox.getMinX(); double minY = boundingBox.getMinY(); long startTime = System.currentTimeMillis(); int pos = 0; for (int j = 0; j < H; ++j) { double y = j*cellHeight + 0.5d*cellHeight + minY; double x = 0.5d*cellWidth + minX; for (int end = pos + W; pos < end; ++pos, x += cellWidth) { queryBuffer.init(x - EPS, x + EPS, y - EPS, y + EPS); center.x = x; center.y = y; finder.prepare(center); spatialIndex.query(queryBuffer, finder); GridCell found = finder.found; if (found == null || depth.depth(center) > 0d) { continue; } double z1 = Interpolation2D.interpolate( found.p1.x, found.p1.z, found.p2.x, found.p2.z, center.x); double z2 = Interpolation2D.interpolate( found.p3.x, found.p3.z, found.p4.x, found.p4.z, center.x); double y1 = Interpolation2D.interpolate( found.p1.x, found.p1.y, found.p2.x, found.p2.y, center.x); double y2 = Interpolation2D.interpolate( found.p3.x, found.p3.y, found.p4.x, found.p4.y, center.x); raster[pos] = Interpolation2D.interpolate( y1, z1, y2, z2, center.y); } } long stopTime = System.currentTimeMillis(); if (debug) { log.debug("interpolation took: " + (stopTime - startTime)/1000f + " secs"); } this.raster = raster; this.width = W; this.height = H; return true; } } // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :