Menu
▾
▴
[Geotools-commits] svn - r38813 - in trunk/modules/unsupported/process-raster/src: main/java/org/geotools/process/raster main/java/org/geotools/process/raster/surface test/java/org/geotools/process/raster test/java/org/geotools/process/raster/surface
Author: mdavis Date: 2012-06-18 11:28:10 -0700 (Mon, 18 Jun 2012) New Revision: 38813 Added: trunk/modules/unsupported/process-raster/src/main/java/org/geotools/process/raster/surface/ trunk/modules/unsupported/process-raster/src/main/java/org/geotools/process/raster/surface/BBOXExpandingFilterVisitor.java trunk/modules/unsupported/process-raster/src/main/java/org/geotools/process/raster/surface/BarnesSurfaceInterpolator.java trunk/modules/unsupported/process-raster/src/main/java/org/geotools/process/raster/surface/BarnesSurfaceProcess.java trunk/modules/unsupported/process-raster/src/main/java/org/geotools/process/raster/surface/BilinearInterpolator.java trunk/modules/unsupported/process-raster/src/main/java/org/geotools/process/raster/surface/GridTransform.java trunk/modules/unsupported/process-raster/src/test/java/org/geotools/process/raster/surface/ trunk/modules/unsupported/process-raster/src/test/java/org/geotools/process/raster/surface/BarnesSurfaceProcessTest.java trunk/modules/unsupported/process-raster/src/test/java/org/geotools/process/raster/surface/BilinearInterpolatorTest.java Log: Added BarnesSurfaceProcess and supporting classes Added: trunk/modules/unsupported/process-raster/src/main/java/org/geotools/process/raster/surface/BBOXExpandingFilterVisitor.java =================================================================== --- trunk/modules/unsupported/process-raster/src/main/java/org/geotools/process/raster/surface/BBOXExpandingFilterVisitor.java (rev 0) +++ trunk/modules/unsupported/process-raster/src/main/java/org/geotools/process/raster/surface/BBOXExpandingFilterVisitor.java 2012-06-18 18:28:10 UTC (rev 38813) @@ -0,0 +1,81 @@ +/* + * GeoTools - The Open Source Java GIS Toolkit + * http://geotools.org + * + * (C) 2011, Open Source Geospatial Foundation (OSGeo) + * (C) 2008-2011 TOPP - www.openplans.org. + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; + * version 2.1 of the License. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + */ +package org.geotools.process.raster.surface; + +import org.geotools.filter.visitor.DuplicatingFilterVisitor; +import org.opengis.filter.expression.Expression; +import org.opengis.filter.spatial.BBOX; + +/** + * A {@link DuplicatingFilterVisitor} which expands the {@link BBOX} of the filter + * by given distances for each box edge. + * + * @author Martin Davis - OpenGeo + * + */ +class BBOXExpandingFilterVisitor extends DuplicatingFilterVisitor { + private double expandMinX; + + private double expandMaxX; + + private double expandMinY; + + private double expandMaxY; + + /** + * Creates a new expanding filter. + * + * @param expandMinX the distance to expand the box X dimension leftwards + * @param expandMaxX the distance to expand the box X dimension rightwards + * @param expandMinY the distance to expand the box Y dimension downwards + * @param expandMaxY the distance to expand the box Y dimension upwards + */ + public BBOXExpandingFilterVisitor(double expandMinX, double expandMaxX, double expandMinY, + double expandMaxY) { + this.expandMinX = expandMinX; + this.expandMaxX = expandMaxX; + this.expandMinY = expandMinY; + this.expandMaxY = expandMaxX; + } + + /** + * Expands the BBOX in the Filter. + * + */ + @SuppressWarnings("deprecation") + @Override + public Object visit(BBOX filter, Object extraData) { + // no need to change the property name + Expression propertyName = filter.getExpression1(); + + /** + * Using the deprecated methods since they are too useful... + */ + double minx = filter.getMinX(); + double miny = filter.getMinY(); + double maxx = filter.getMaxX(); + double maxy = filter.getMaxY(); + String srs = filter.getSRS(); + + return getFactory(extraData).bbox(propertyName, + minx - expandMinX, miny - expandMaxX, + maxx + expandMinY, maxy + expandMaxY, + srs); + } + +} \ No newline at end of file Added: trunk/modules/unsupported/process-raster/src/main/java/org/geotools/process/raster/surface/BarnesSurfaceInterpolator.java =================================================================== --- trunk/modules/unsupported/process-raster/src/main/java/org/geotools/process/raster/surface/BarnesSurfaceInterpolator.java (rev 0) +++ trunk/modules/unsupported/process-raster/src/main/java/org/geotools/process/raster/surface/BarnesSurfaceInterpolator.java 2012-06-18 18:28:10 UTC (rev 38813) @@ -0,0 +1,472 @@ +/* + * GeoTools - The Open Source Java GIS Toolkit + * http://geotools.org + * + * (C) 2011, Open Source Geospatial Foundation (OSGeo) + * (C) 2008-2011 TOPP - www.openplans.org. + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; + * version 2.1 of the License. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + */ +package org.geotools.process.raster.surface; + +import com.vividsolutions.jts.geom.Coordinate; +import com.vividsolutions.jts.geom.Envelope; + +/** + * Interpolates a surface across a regular grid from an irregular set of data points + * using the Barnes Surface Interpolation technique. + * <p> + * Barnes Surface Interpolation is a surface estimating method + * commonly used as an interpolation technique for meteorological datasets. + * The algorithm operates on a regular grid of cells covering a specified extent in the input data space. + * It computes an initial pass to produce an averaged (smoothed) value for each cell in the grid, + * based on the cell's proximity to the points in the input observations. + * Subsequent refinement passes may be performed to improve the surface estimate + * to better approximate the observed values. + * <ul> + * <li>The initial pass produces an averaged (smoothed) value for each grid cell + * using a summation of exponential (Gaussian) decay functions around each observation point. + * <li>Subsequent refinement passes compute an error surface in the same way, using the deltas + * between the previous estimated surface and the observations. + * The error surface is added to the previous estimated surface to refine the estimate + * (by reducing the delta between the estimate and the observations). + * </ul> + * <p> + * For the first pass, the estimated value at each grid cell is: + * + * <pre> + * E<sub>g</sub> = sum(w<sub>i</sub> * o<sub>i</sub>) / sum(w<sub>i</sub>) + * </pre> + * + * where + * <ul> + * <li><code>E<sub>g</sub></code> is the estimated surface value at the grid cell + * <li><code>w<sub>i</sub></code> is the weight value for the i'th observation point (see below for definition) + * <li><code>o<sub>i</sub></code> is the value of the i'th observation point + * </ul> + * <p> + * The weight (decay) function used is: + * + * <pre> + * w<sub>i</sub> = exp(-d<sub>i</sub><sup>2</sup> / L<sup>2</sup>c ) + * </pre> + * + * where: + * <ul> + * <li><code>w<sub>i</sub></code> is the <b>weight</b> of the i'th observation point value + * <li><code>d<sub>i</sub></code> is the <b>distance</b> from the grid cell being estimated to the i'th observation point + * <li><code>L</code> is the <b>length scale</b>, which is determined by the observation spacing + * and the natural scale of the phenomena being measured. + * The length scale is in the units of the coordinate system of the data points. + * It will likely need to be empirically estimated. + * <li><code>c</code> is the <b>convergence factor</b>, which controls how much refinement takes place during each refinement step. + * In the first pass the convergence is automatically set to 1. + * For subsequent passes a value in the range 0.2 - 0.3 is usually effective. + * </ul> + * During refinement passes the value at each grid cell is re-estimated as: + * + * <pre> + * E<sub>g</sub>' = E<sub>g</sub> + sum( w<sub>i</sub> * (o<sub>i</sub> - E<sub>i</sub>) ) / sum( w<sub>i</sub> ) + * </pre> + * + * To optimize performance for large input datasets, it is only necessary to provide + * the data points which affect the surface interpolation within the + * specified output extent. In order to avoid "edge effects", the provided data points + * should be taken from an area somewhat larger than the output extent. + * The extent of the data area depends on the length scale, convergence factor, and data spacing in a complex way. + * A reasonable heuristic for determining the size of the query extent is to expand the output extent by a value of 2L. + * <p> + * Since the visual quality and accuracy of the computed surface is lower further from valid observations, + * the algorithm allows limiting the extent of the + * computed cells. This is done by using the concept of <b>supported grid cells</b>. + * Grid cells are supported by the + * input observations if they are within a specified distance of a specified number of observation points. + * Grid cells which are not supported are not + * computed and are output as NO_DATA values. + * <p> + * <b>References</b> + * <ol> + * <li>Barnes, S. L (1964). "A technique for maximizing details in numerical weather-map analysis". <i>Journal of Applied Meterology</i> 3 (4): 396 - 409 + * </ol> + * + * @author Martin Davis - OpenGeo + * + */ +public class BarnesSurfaceInterpolator { + + /** + * The default grid cell value used to indicate no data was computed for that cell + */ + public static final float DEFAULT_NO_DATA_VALUE = -999; + + private static final double INTERNAL_NO_DATA = Double.NaN; + + // =========== Input parameters + /** + * These parameters control which grid points are considered to be supported, i.e. have enough nearby observation points to be reasonably + * estimated. + * + * A grid point is supported if it has: + * + * count(obs within maxObservationDistance) >= minObservationCount + * + * Using these parameters is optional, but recommended, since estimating grid points which are far from any observations can produce unrealistic + * surfaces. + */ + private int minObservationCount = 2; + + private double maxObservationDistance = 0.0; + + private double convergenceFactor = 0.3; + + private double lengthScale = 0.0; + + private int passCount = 1; + + private Coordinate[] inputObs; + + // ============= Internal parameters (could be exposed) + private float noDataValue = DEFAULT_NO_DATA_VALUE; + + // ============= Computed parameters + // private double effectiveRadius; + + /** + * Indicates whether estimated grid points are filtered based on distance from observations + */ + private boolean useObservationMask; + + // ============ Working data + private float[] estimatedObs; + + /** + * Creates a Barnes Interpolator over a specified dataset of observation values. The observation data is provided as an array of + * {@link Coordinate} values, where the X,Y ordinates are the observation location, and the Z ordinate contains the observation value. + * + * @param data the observed data values + */ + public BarnesSurfaceInterpolator(Coordinate[] observationData) { + this.inputObs = observationData; + } + + /** + * Sets the number of passes performed during Barnes interpolation. + * + * @param passCount the number of estimation passes to perform (1 or more) + */ + public void setPassCount(int passCount) { + if (passCount < 1) + return; + this.passCount = passCount; + } + + /** + * Sets the length scale for the interpolation weighting function. The length scale is determined from the distance between the observation + * points, as well as the scale of the phenomena which is being measured. + * <p> + * + * + * @param lengthScale + */ + public void setLengthScale(double lengthScale) { + this.lengthScale = lengthScale; + } + + /** + * Sets the convergence factor used during refinement passes. The value should be in the range [0,1]. Empirically, values between 0.2 - 0.3 are + * most effective. Smaller values tend to make the interpolated surface too "jittery". Larger values produce less refinement effect. + * + * @param convergenceFactor the factor determining how much to refine the surface estimate + */ + public void setConvergenceFactor(double convergenceFactor) { + this.convergenceFactor = convergenceFactor; + } + + /** + * Sets the maximum distance from an observation for a grid point to be supported by that observation. Empirically determined; a reasonable + * starting point is between 1.5 and 2 times the Length scale. If the value is 0 (which is the default), all grid points are considered to be + * supported, and will thus be computed. + * + * @param maxObsDistance the maximum distance from an observation for a supported grid point + */ + public void setMaxObservationDistance(double maxObsDistance) { + this.maxObservationDistance = maxObsDistance; + } + + /** + * Sets the minimum number of in-range observations which are required for a grid point to be supported. The default is 2. + * + * @param minObsCount the minimum in-range observation count for supported grid points + */ + public void setMinObservationCount(int minObsCount) { + this.minObservationCount = minObsCount; + } + + /** + * Sets the NO_DATA value used to indicate that a grid cell was not computed. This value should be distinct from any potential data value. + * + * @param noDataValue the value to use to represent NO_DATA. + */ + public void setNoData(float noDataValue) { + this.noDataValue = noDataValue; + } + + /** + * Computes the estimated values for a regular grid of cells. The area covered by the grid is specified by an {@link Envelope}. The size of the + * grid is specified by the cell count for the grid width (X) and height (Y). + * + * @param srcEnv the area covered by the grid + * @param xSize the width of the grid + * @param ySize the height of the grid + * + * @return the computed grid of estimated data values (in row-major order) + */ + public float[][] computeSurface(Envelope srcEnv, int xSize, int ySize) { + // not currently used + // effectiveRadius = effectiveRadius(minimumWeight, influenceRadius); + + useObservationMask = minObservationCount > 0 && maxObservationDistance > 0.0; + + float[][] grid = new float[xSize][ySize]; + GridTransform trans = new GridTransform(srcEnv, xSize, ySize); + + estimateGrid(grid, trans); + + if (passCount > 1) { + /** + * First refinement pass requires observation points to be estimated as well + */ + estimatedObs = computeEstimatedObservations(); + refineGrid(grid, trans); + + /** + * For subsequent refinement passes, refine observations then recompute + */ + for (int i = 3; i <= passCount; i++) { + refineEstimatedObservations(estimatedObs); + refineGrid(grid, trans); + } + } + return grid; + } + + private float[] computeEstimatedObservations() { + float[] estimate = new float[inputObs.length]; + for (int i = 0; i < inputObs.length; i++) { + Coordinate dp = inputObs[i]; + float est = (float) estimatedValue(dp.x, dp.y); + if (! Float.isNaN(est)) + estimate[i] = est; + else + estimate[i] = (float) inputObs[i].z; + } + return estimate; + } + + private float[] refineEstimatedObservations(float[] currEst) { + float[] estimate = new float[inputObs.length]; + for (int i = 0; i < inputObs.length; i++) { + Coordinate dp = inputObs[i]; + float del = (float) refinedDelta(dp.x, dp.y, convergenceFactor); + if (! Float.isNaN(del)) + estimate[i] = (float) currEst[i] + del; + else + estimate[i] = (float) inputObs[i].z; + } + return estimate; + } + + /** + * Computes an initial estimate of the interpolated surface. + * + * @param grid the grid matrix buffer to use + * @param trans the transform mapping from data space to the grid + */ + private void estimateGrid(float[][] grid, GridTransform trans) { + for (int i = 0; i < grid.length; i++) { + for (int j = 0; j < grid[0].length; j++) { + double x = trans.x(i); + double y = trans.y(j); + + grid[i][j] = (float) noDataValue; + if (useObservationMask && !isSupportedGridPt(x, y)) + continue; + + float est = (float) estimatedValue(x, y); + if (!Float.isNaN(est)) + grid[i][j] = est; + } + } + } + + /** + * Computes a refined estimate for the interpolated surface. + * + * @param grid the grid matrix buffer to use + * @param trans the transform mapping from data space to the grid + */ + private void refineGrid(float[][] grid, GridTransform trans) { + for (int i = 0; i < grid.length; i++) { + for (int j = 0; j < grid[0].length; j++) { + double x = trans.x(i); + double y = trans.y(j); + + // skip NO_DATA values + if (grid[i][j] == noDataValue) + continue; + + float del = (float) refinedDelta(x, y, convergenceFactor); + /* + // DEBUGGING + if (del < 0) { + float d = (float) refinedDelta(x, y, convergenceFactor); + } + */ + if (! Float.isNaN(del)) + grid[i][j] = grid[i][j] + del; + } + } + } + + private boolean isSupportedGridPt(double x, double y) { + int count = 0; + for (int i = 0; i < inputObs.length; i++) { + double dist = distance(x, y, inputObs[i]); + if (dist <= maxObservationDistance) + count++; + } + return count >= minObservationCount; + } + + private double distance(double x, double y, Coordinate p) { + double dx = x - p.x; + double dy = y - p.y; + return Math.sqrt(dx * dx + dy * dy); + } + + /** + * Computes the initial estimate for a grid point. + * + * @param x the x ordinate of the grid point location + * @param y the y ordinate of the grid point location + * @return the estimated value, or INTERNAL_NO_DATA if the grid cell is not supported + */ + private double estimatedValue(double x, double y) { + Coordinate p = new Coordinate(x, y); + + double sumWgtVal = 0; + double sumWgt = 0; + int dataCount = 0; + for (int i = 0; i < inputObs.length; i++) { + double wgt = weight(p, inputObs[i], lengthScale); + /** + * Skip observation if unusable due to too great a distance + */ + if (Double.isNaN(wgt)) + continue; + + sumWgtVal += wgt * inputObs[i].z; + sumWgt += wgt; + dataCount++; + } + /** + * If grid point is not supported, return NO_DATA + */ + if (dataCount < minObservationCount) + return INTERNAL_NO_DATA; + return sumWgtVal / sumWgt; + } + + /** + * Computes a refinement delta, which is added to a grid point estimated value + * to refine the estimate. + * + * @param x the x ordinate of the grid point location + * @param y the y ordinate of the grid point location + * @param convergenceFactor the convergence factor + * @return the refinement delta value, or INTERNAL_NO_DATA if the grid cell is not supported + */ + private double refinedDelta(double x, double y, double convergenceFactor) { + Coordinate p = new Coordinate(x, y); + + double sumWgtVal = 0; + double sumWgt = 0; + int dataCount = 0; + for (int i = 0; i < inputObs.length; i++) { + double wgt = weight(p, inputObs[i], lengthScale, convergenceFactor); + /** + * Check if observation is unusable (e.g. due to too great a distance) + */ + if (Double.isNaN(wgt)) + continue; + + sumWgtVal += wgt * (inputObs[i].z - estimatedObs[i]); + sumWgt += wgt; + dataCount++; + } + /** + * If grid point is not supported, return NO_DATA + */ + if (dataCount < minObservationCount) + return INTERNAL_NO_DATA; + return sumWgtVal / sumWgt; + } + + private double weight(Coordinate dataPt, Coordinate gridPt, double lengthScale) { + return weight(gridPt, dataPt, lengthScale, 1.0); + } + + private double weight(Coordinate dataPt, Coordinate gridPt, double lengthScale, + double convergenceFactor) { + double dist = dataPt.distance(gridPt); + return weight(dist, lengthScale, convergenceFactor); + } + + private double weight(double dist, double lengthScale, double convergenceFactor) { + /** + * MD - using an effective radius is problematic. + * + * The effective radius grows as a log function of the cutoff weight, so even for very small cutoff weight values, the effective radius is + * only a few times the size of the influence radius. + * + * Also, dropping observation terms from the estimate results in very drastic (discontinuous) changes at distances around the effective + * radius. (Probably because beyond that distance there are very few terms (maybe only 2) contributing to the estimate, so there is no + * smoothing effect from incorporating many estimates) + * + * So - don't use effectiveRadius. + * + * Or, maybe it's ok as long as a observation mask is used as well, since the effect only occurs at large distances from observation points? + */ + /* + * if (dist > effectiveRadius) return INTERNAL_NO_DATA; // + */ + double dr = dist / lengthScale; + double w = Math.exp(-(dr * dr / convergenceFactor)); + // if (dist > cutoffRadius) System.out.println(w); + return w; + } + + /** + * Computes effective radius which is determined by the specified cutoff weight and the radius of the decay function. + * + * @param cutoffWeight + * @param radius + * @return + */ + private double effectiveRadius(double cutoffWeight, double radius) { + double cutoffFactor = Math.sqrt(-Math.log(cutoffWeight)); + double effRadius = radius * cutoffFactor; + double w = weight(effRadius, radius, 1.0); + System.out.println(cutoffWeight + " " + w); + return effRadius; + } + +} Added: trunk/modules/unsupported/process-raster/src/main/java/org/geotools/process/raster/surface/BarnesSurfaceProcess.java =================================================================== --- trunk/modules/unsupported/process-raster/src/main/java/org/geotools/process/raster/surface/BarnesSurfaceProcess.java (rev 0) +++ trunk/modules/unsupported/process-raster/src/main/java/org/geotools/process/raster/surface/BarnesSurfaceProcess.java 2012-06-18 18:28:10 UTC (rev 38813) @@ -0,0 +1,439 @@ +/* + * GeoTools - The Open Source Java GIS Toolkit + * http://geotools.org + * + * (C) 2011, Open Source Geospatial Foundation (OSGeo) + * (C) 2008-2011 TOPP - www.openplans.org. + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; + * version 2.1 of the License. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + */ +package org.geotools.process.raster.surface; + +import java.util.ArrayList; +import java.util.List; + +import javax.measure.unit.NonSI; +import javax.measure.unit.SI; +import javax.measure.unit.Unit; + +import org.geotools.coverage.CoverageFactoryFinder; +import org.geotools.coverage.grid.GridCoverage2D; +import org.geotools.coverage.grid.GridCoverageFactory; +import org.geotools.data.Query; +import org.geotools.data.simple.SimpleFeatureCollection; +import org.geotools.data.simple.SimpleFeatureIterator; +import org.geotools.factory.Hints; +import org.geotools.filter.text.cql2.CQLException; +import org.geotools.filter.text.ecql.ECQL; +import org.geotools.geometry.jts.ReferencedEnvelope; +import org.geotools.process.ProcessException; +import org.geotools.process.factory.DescribeParameter; +import org.geotools.process.factory.DescribeProcess; +import org.geotools.process.factory.DescribeResult; +import org.geotools.process.gs.GSProcess; +import org.geotools.referencing.CRS; +import org.opengis.coverage.grid.GridCoverage; +import org.opengis.coverage.grid.GridGeometry; +import org.opengis.feature.simple.SimpleFeature; +import org.opengis.filter.Filter; +import org.opengis.filter.expression.Expression; +import org.opengis.referencing.FactoryException; +import org.opengis.referencing.crs.CoordinateReferenceSystem; +import org.opengis.referencing.operation.MathTransform; +import org.opengis.util.ProgressListener; + +import com.vividsolutions.jts.geom.Coordinate; +import com.vividsolutions.jts.geom.CoordinateArrays; +import com.vividsolutions.jts.geom.Envelope; +import com.vividsolutions.jts.geom.Geometry; + +/** + * A Process that uses a {@link BarnesSurfaceInterpolator} to compute an interpolated surface + * over a set of irregular data points as a {@link GridCoverage}. + * <p> + * The implementation allows limiting the radius of influence of observations, in order to + * prevent extrapolation into unsupported areas, and to increase performance (by reducing + * the number of observations considered). + * <p> + * To improve performance, the surface grid can be computed at a lower resolution than the requested output image. + * The grid is upsampled to match the required image size. + * Upsampling uses Bilinear Interpolation to maintain visual quality. + * This gives a large improvement in performance, with minimal impact + * on visual quality for small cell sizes (for instance, 10 pixels or less). + * + * To ensure that the computed surface is stable + * (i.e. does not display obvious edge artifacts during zooming, panning and tiling), + * the data query extent should be expanded to be larger than the specified output extent. + * This includes "nearby" points which may affect the value of the surface. + * The expansion distance depends on the + * length scale, convergence factor, and data spacing + * in a complex way, so must be manually determined. + * It does NOT depend on the output window extent. + * (A good heuristic is to set it to expand by at least the size of the length scale.) + * + * To prevent excessive CPU consumption, the process allows limiting the number of data points + * to process. If the limit is exceeded the output is computed consuming and using only the + * maximum number of points specified. + * + * <h3>Parameters</h3> + * <i>M = mandatory, O = optional</i> + * <p> + * <ul> + * <li><b>data</b> (M) - the FeatureCollection containing the point observations + * <li><b>valueAttr</b> (M)- the feature type attribute containing the observed surface value + * <li><b>dataLimit</b> (O)- the maximum number of input points to process + * <li><b>scale</b> (M) - the Length Scale for the interpolation. In units of the input data CRS. + * <li><b>convergence</b> (O) - the convergence factor for refinement. Between 0 and 1 (values below 0.4 are safest). (Default = 0.3) + * <li><b>passes</b> (O) - the number of passes to compute. 1 or greater. (Default = 2) + * <li><b>minObservations</b> (O) - The minimum number of observations required to support a grid cell. (Default = 2) + * <li><b>maxObservationDistance</b> (O) - The maximum distance to an observation for it to support a grid cell. 0 means all observations are used. In units of the input data CRS. (Default = 0) + * <li><b>noDataValue</b> (O) - The NO_DATA value to use for unsupported grid cells in the output coverage. (Default = -999) + * <li><b>pixelsPerCell</b> (O) - The pixels-per-cell value determines the resolution of the computed grid. + * Larger values improve performance, but may degrade appearance. (Default = 1) + * <li><b>queryBuffer</b> (O) - The distance to expand the query envelope by. Larger values provide a more stable surface. In units of the input data CRS. (Default = 0) + * <li><b>outputBBOX</b> (M) - The georeferenced bounding box of the output area + * <li><b>outputWidth</b> (M) - The width of the output raster + * <li><b>outputHeight</b> (M) - The height of the output raster + * </ul> + * The output of the process is a {@linkplain GridCoverage2D} with a single band, + * with cell values in the same domain as the input observation field specified by <code>valueAttr</code>. + * <p> + * Computation of the surface takes places in the CRS of the output. + * If the data CRS is geodetic and the output CRS is planar, or vice-versa, + * the input points are transformed into the output CRS. + * A simple technique is used to convert the surface distance parameters + * <code>scale</code> and <code>maxObservationDistance</code> into the output CRS units. + * + * <h3>Using the process as a Rendering Transformation</h3> + * + * This process can be used as a RenderingTransformation, since it + * implements the <tt>invertQuery(... Query, GridGeometry)</tt> method. + * <p> + * When used as an Rendering Transformation the process rewrites data query to expand the query BBOX. + * This includes "nearby" data points to make the + * computed surface stable under panning and zooming. + * To support this the <code>queryBuffer</code> parameter should be specified to expand + * the query extent appropriately. + * <p> + * The output raster parameters can be determined from the request extents, using the + * following SLD environment variables: + * <p> + * <ul> + * <li><b>outputBBOX</b> - env var = <tt>wms_bbox</tt> + * <li><b>outputWidth</b> - env var = <tt>wms_width</tt> + * <li><b>outputHeight</b> - env var = <tt>wms_height</tt> + * </ul> + * + * <p> + * @author Martin Davis - OpenGeo + * + */ + +@DescribeProcess(title = "BarnesSurface", description = "Uses Barnes Analysis to compute an interpolated surface over a set of irregular data points aa a GridCoverage.") +public class BarnesSurfaceProcess implements GSProcess { + + // no process state is defined, since RenderingTransformation processes must be stateless + + @DescribeResult(name = "result", description = "The interpolated surface as a raster") + public GridCoverage2D execute( + + // process data + @DescribeParameter(name = "data", description = "Features containing the point observations to be interpolated") SimpleFeatureCollection obsFeatures, + @DescribeParameter(name = "valueAttr", description = "Featuretype attribute containing the observed surface value") String valueAttr, + @DescribeParameter(name = "dataLimit", description = "Limits the number of input features processed", min=0, max=1) Integer argDataLimit, + + // process parameters + @DescribeParameter(name = "scale", description = "Length scale to use for the interpolation", min=1, max=1) Double argScale, + @DescribeParameter(name = "convergence", description = "Convergence factor for the interpolation (default: 0.3)", min=0, max=1) Double argConvergence, + @DescribeParameter(name = "passes", description = "Number of passes to compute (default: 2)", min=0, max=1) Integer argPasses, + @DescribeParameter(name = "minObservations", description = "Minimum number of observations required to support a grid cell (default: 2)", min=0, max=1) Integer argMinObsCount, + @DescribeParameter(name = "maxObservationDistance", description = "Maximum distance to a supporting observation (default: 0)", min=0, max=1) Double argMaxObsDistance, + @DescribeParameter(name = "noDataValue", description = "Value to use for NO_DATA cells (default: -999)", min=0, max=1) Double argNoDataValue, + @DescribeParameter(name = "pixelsPerCell", description = "Number of pixels per grid cell (default = 1)", min=0, max=1) Integer argPixelsPerCell, + + // query modification parameters + @DescribeParameter(name = "queryBuffer", description = "Distance by which to expand the query window (default: 0)", min=0, max=1) Double argQueryBuffer, + + // output image parameters + @DescribeParameter(name = "outputBBOX", description = "Georeferenced bounding box of the output") ReferencedEnvelope outputEnv, + @DescribeParameter(name = "outputWidth", description = "Width of the output raster") Integer outputWidth, + @DescribeParameter(name = "outputHeight", description = "Height of the output raster") Integer outputHeight, + + ProgressListener monitor) throws ProcessException { + + /**--------------------------------------------- + * Check that process arguments are valid + * --------------------------------------------- + */ + if (valueAttr == null || valueAttr.length() <= 0) { + throw new IllegalArgumentException("Value attribute must be specified"); + } + + /**--------------------------------------------- + * Set up required information from process arguments. + * --------------------------------------------- + */ + int dataLimit = 0; + if (argDataLimit != null) dataLimit = argDataLimit; + + double lengthScale = argScale; + double convergenceFactor = argConvergence != null ? argConvergence : 0.3; + int passes = argPasses != null ? argPasses : 2; + int minObsCount = argMinObsCount != null ? argMinObsCount : 2; + double maxObsDistance = argMaxObsDistance != null ? argMaxObsDistance : 0.0; + float noDataValue = (float) (argNoDataValue != null ? argNoDataValue : -999); + int pixelsPerCell = 1; + if (argPixelsPerCell != null && argPixelsPerCell > 1) { + pixelsPerCell = argPixelsPerCell; + } + int gridWidth = outputWidth; + int gridHeight = outputHeight; + if (pixelsPerCell > 1) { + gridWidth = outputWidth / pixelsPerCell; + gridHeight = outputHeight / pixelsPerCell; + } + + CoordinateReferenceSystem srcCRS = obsFeatures.getSchema().getCoordinateReferenceSystem(); + CoordinateReferenceSystem dstCRS = outputEnv.getCoordinateReferenceSystem(); + MathTransform trans = null; + try { + trans = CRS.findMathTransform(srcCRS, dstCRS); + } catch (FactoryException e) { + throw new ProcessException(e); + } + /**--------------------------------------------- + * Convert distance parameters to units of the destination CRS. + * --------------------------------------------- + */ + double distanceConversionFactor = distanceConversionFactor(srcCRS, dstCRS); + double dstLengthScale = lengthScale * distanceConversionFactor; + double dstMaxObsDistance = maxObsDistance * distanceConversionFactor; + + /**--------------------------------------------- + * Extract the input observation points + * --------------------------------------------- + */ + Coordinate[] pts = null; + try { + pts = extractPoints(obsFeatures, valueAttr, trans, dataLimit); + } catch (CQLException e) { + throw new ProcessException(e); + } + + /**--------------------------------------------- + * Do the processing + * --------------------------------------------- + */ + //Stopwatch sw = new Stopwatch(); + // interpolate the surface at the specified resolution + float[][] barnesGrid = createBarnesGrid(pts, dstLengthScale, convergenceFactor, passes, minObsCount, dstMaxObsDistance, noDataValue, outputEnv, gridWidth, gridHeight); + + // flip now, since grid size may be smaller + barnesGrid = flipXY(barnesGrid); + + // upsample to output resolution if necessary + float[][] outGrid = barnesGrid; + if (pixelsPerCell > 1) + outGrid = upsample(barnesGrid, noDataValue, outputWidth, outputHeight); + + // convert to the GridCoverage2D required for output + GridCoverageFactory gcf = CoverageFactoryFinder.getGridCoverageFactory(null); + GridCoverage2D gridCov = gcf.create("Process Results", outGrid, outputEnv); + + //System.out.println("************** Barnes Surface computed in " + sw.getTimeString()); + + return gridCov; + } + + /* + * An approximate value for the length of a degree at the equator in meters. + * This doesn't have to be precise, since it is only used to convert + * values which are themselves rough approximations. + */ + private static final double METRES_PER_DEGREE = 111320; + + private static double distanceConversionFactor(CoordinateReferenceSystem srcCRS,CoordinateReferenceSystem dstCRS) + { + Unit<?> srcUnit = srcCRS.getCoordinateSystem().getAxis(0).getUnit(); + Unit<?> dstUnit = dstCRS.getCoordinateSystem().getAxis(0).getUnit(); + if (srcUnit == dstUnit) { + return 1; + } + else if (srcUnit == NonSI.DEGREE_ANGLE && dstUnit == SI.METER) { + return METRES_PER_DEGREE; + } + else if (srcUnit == SI.METER && dstUnit == NonSI.DEGREE_ANGLE) { + return 1.0 / METRES_PER_DEGREE; + } + throw new IllegalStateException("Unable to convert distances from " + srcUnit + " to " + dstUnit); + } + + /** + * Flips an XY matrix along the X=Y axis, and inverts the Y axis. + * Used to convert from "map orientation" into the "image orientation" + * used by GridCoverageFactory. + * The surface interpolation is done on an XY grid, with Y=0 being the bottom of the space. + * GridCoverages are stored in an image format, in a YX grid with 0 being the top. + * + * @param grid the grid to flip + * @return the flipped grid + */ + private static float[][] flipXY(float[][] grid) + { + int xsize = grid.length; + int ysize = grid[0].length; + + float[][] grid2 = new float[ysize][xsize]; + for (int ix = 0; ix < xsize; ix++) { + for (int iy = 0; iy < ysize; iy++) { + int iy2 = ysize - iy - 1; + grid2[iy2][ix] = grid[ix][iy]; + } + } + return grid2; + } + + private float[][] createBarnesGrid(Coordinate[] pts, + double lengthScale, + double convergenceFactor, + int passes, + int minObservationCount, + double maxObservationDistance, + float noDataValue, + Envelope destEnv, + int width, int height) + { + BarnesSurfaceInterpolator barnesInterp = new BarnesSurfaceInterpolator(pts); + barnesInterp.setLengthScale(lengthScale); + barnesInterp.setConvergenceFactor(convergenceFactor); + barnesInterp.setPassCount(passes); + barnesInterp.setMinObservationCount(minObservationCount); + barnesInterp.setMaxObservationDistance(maxObservationDistance); + barnesInterp.setNoData(noDataValue); + + float[][] grid = barnesInterp.computeSurface(destEnv, width, height); + + return grid; + } + + private float[][] upsample(float[][] grid, + float noDataValue, + int width, + int height) { + BilinearInterpolator bi = new BilinearInterpolator(grid, noDataValue); + float[][] outGrid = bi.interpolate(width, height, true); + return outGrid; + } + + /** + * Given a target query and a target grid geometry + * returns the query to be used to read the input data of the process involved in rendering. In + * this process this method is used to: + * <ul> + * <li>determine the extent & CRS of the output grid + * <li>expand the query envelope to ensure stable surface generation + * <li>modify the query hints to ensure point features are returned + * </ul> + * Note that in order to pass validation, all parameters named here must also appear + * in the parameter list of the <tt>execute</tt> method, + * even if they are not used there. + * + * @param argQueryBuffer the distance by which to expand the query window + * @param targetQuery the query used against the data source + * @param targetGridGeometry the grid geometry of the destination image + * @return The transformed query + */ + public Query invertQuery( + @DescribeParameter(name = "queryBuffer", description = "The distance by which to expand the query window", min=0, max=1) Double argQueryBuffer, + Query targetQuery, GridGeometry targetGridGeometry) + throws ProcessException { + + // default is no expansion + double queryBuffer = 0; + if (argQueryBuffer != null) { + queryBuffer = argQueryBuffer; + } + + targetQuery.setFilter(expandBBox(targetQuery.getFilter(), queryBuffer)); + + // clear properties to force all attributes to be read + // (required because the SLD processor cannot see the value attribute specified in the transformation) + // TODO: set the properties to read only the specified value attribute + targetQuery.setProperties(null); + + // set the decimation hint to ensure points are read + Hints hints = targetQuery.getHints(); + hints.put(Hints.GEOMETRY_DISTANCE, 0.0); + + return targetQuery; + } + + private Filter expandBBox(Filter filter, double distance) { + return (Filter) filter.accept( + new BBOXExpandingFilterVisitor(distance, distance, distance, distance), null); + } + + public static Coordinate[] extractPoints(SimpleFeatureCollection obsPoints, String attrName, MathTransform trans, int dataLimit) throws CQLException + { + Expression attrExpr = ECQL.toExpression(attrName); + List<Coordinate> ptList = new ArrayList<Coordinate>(); + SimpleFeatureIterator obsIt = obsPoints.features(); + + double[] srcPt = new double[2]; + double[] dstPt = new double[2]; + + int i = 0; + try { + while (obsIt.hasNext()) { + SimpleFeature feature = obsIt.next(); + + double val = 0; + + try { + + if (dataLimit > 0 && i >= dataLimit) { + //TODO: log this situation + break; + } + i++; + // get the observation value from the attribute (if non-null) + Object valObj = attrExpr.evaluate(feature); + if (valObj != null) { + // System.out.println(evaluate); + Number valNum = (Number) valObj; + val = valNum.doubleValue(); + + // get the point location from the geometry + Geometry geom = (Geometry) feature.getDefaultGeometry(); + Coordinate p = geom.getCoordinate(); + srcPt[0] = p.x; + srcPt[1] = p.y; + trans.transform(srcPt, 0, dstPt, 0, 1); + Coordinate pobs = new Coordinate(dstPt[0], dstPt[1], val); + //Coordinate pobs = new Coordinate(p.x, p.y, val); + ptList.add(pobs); + } + } + catch (Exception e) { + // just carry on for now (debugging) + //throw new ProcessException("Expression " + attrExpr + " failed to evaluate to a numeric value", e); + } + } + } + finally { + obsIt.close(); + } + + Coordinate[] pts = CoordinateArrays.toCoordinateArray(ptList); + return pts; + } + +} Added: trunk/modules/unsupported/process-raster/src/main/java/org/geotools/process/raster/surface/BilinearInterpolator.java =================================================================== --- trunk/modules/unsupported/process-raster/src/main/java/org/geotools/process/raster/surface/BilinearInterpolator.java (rev 0) +++ trunk/modules/unsupported/process-raster/src/main/java/org/geotools/process/raster/surface/BilinearInterpolator.java 2012-06-18 18:28:10 UTC (rev 38813) @@ -0,0 +1,203 @@ +/* + * GeoTools - The Open Source Java GIS Toolkit + * http://geotools.org + * + * (C) 2011, Open Source Geospatial Foundation (OSGeo) + * (C) 2008-2011 TOPP - www.openplans.org. + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; + * version 2.1 of the License. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + */ +package org.geotools.process.raster.surface; + +/** + * Interpolates a grid to a grid of different dimensions using bilinear interpolation. + * <p> + * NO_DATA cell values are supported in the source grid. + * There are two ways of handling the boundary between NO_DATA cells and data cells: + * <ul> + * <li><b>Truncate</b> - If any source cell is NO_DATA, the dest value is NO_DATA. + * This is simple and fast, but does make the data boundaries look a bit ragged. + * <li><b>Smooth</b> - If only one source cell is NO_DATA, the value is interpolated + * using the 3 valid values, across one-half of the interpolated cells. + * This smoothes off the boundary. + * If 2 or more source cells are NO_DATA, then Truncation is used. + * </ul> + * <p> + * Reference: http://en.wikipedia.org/wiki/Bilinear_interpolation. + * + * @author Martin Davis - OpenGeo + * + */ +public class BilinearInterpolator { + + private static final float NULL_NO_DATA = Float.NaN; + private final float[][] src; + private float noDataValue = NULL_NO_DATA; + + /** + * Creates a new interpolator on a given source grid. + * + * @param src the source grid + */ + public BilinearInterpolator(final float[][] src) + { + this(src, NULL_NO_DATA); + } + + /** + * Creates a new interpolator on a given source grid. + * + * @param src the source grid + * @param noDataValue the NO_DATA value (if none use Float.NaN) + */ + public BilinearInterpolator(final float[][] src, final float noDataValue) + { + this.src = src; + this.noDataValue = noDataValue; + } + + /** + * Interpolates the source grid to a new grid of different dimensions. + * + * @param width the width of the destination grid + * @param height the height of the destination grid + * @param smoothBoundary true if boundary smoothing should be performed + * @return the interpolated grid + */ + public float[][] interpolate(final int width, final int height, boolean smoothBoundary) + { + int srcWidth = src.length; + int srcHeight = src[0].length; + + float[][] dest = new float[width][height]; + + float xRatio = ((float) srcWidth - 1) / width ; + float yRatio = ((float) srcHeight - 1) / height ; + + for (int i = 0; i < width; i++) { + for (int j = 0; j < height; j++) { + float x = i * xRatio; + float y = j * yRatio; + int ix = (int) x; + int iy = (int) y; + float xfrac = x - ix; + float yfrac = y - iy; + + float val; + + if (ix < srcWidth - 1 && iy < srcHeight - 1) { + // interpolate if src cell is in grid + float v00 = src[ix][iy]; + float v10 = src[ix + 1][iy]; + float v01 = src[ix][iy + 1]; + float v11 = src[ix + 1][iy + 1]; + if (v00 == noDataValue + || v10 == noDataValue + || v01 == noDataValue + || v11 == noDataValue) { + // handle src cell with NO_DATA value(s) + if (smoothBoundary) { + val = interpolateBoundaryCell(xfrac, yfrac, v00, v10, v01, v11); + } + else { + val = noDataValue; + } + } + else { + // All src cell corners have values + // Compute bilinear interpolation over the src cell + val = ( v00*(1-xfrac)*(1-yfrac) + v10*(xfrac)*(1-yfrac) + + v01*(yfrac)*(1-xfrac) + v11*(xfrac*yfrac) + ) ; + } + } + else { + // dest index at edge of grid + val = src[ix][iy]; + } + dest[i][j] = val; + } + } + return dest; + } + + /** + * Interpolates across an edge grid cell, which has 1 or more NO_DATA values. + * Grid cells with 2 or or NO_DATA values still receive the value NO_DATA. + * Otherwise, the cell is interpolated across the triangle defined by the + * 3 valid corner values. + * This produces a much smoother edge appearance, containing 45-degree lines + * instead of a jagged stairstep boundary. + * + * @param xfrac the fractional x location of the interpolation point within the cell + * @param yfrac the fractional y location of the interpolation point within the cell + * @param v00 the lower left value + * @param v10 the lower right value + * @param v01 the upper left value + * @param v11 the upper right value + * @return the interpolated value + */ + private float interpolateBoundaryCell(float xfrac, float yfrac, float v00, float v10, float v01, float v11) + { + // count noData values + int count = 0; + if (v00 == noDataValue) count++; + if (v10 == noDataValue) count++; + if (v01 == noDataValue) count++; + if (v11 == noDataValue) count++; + + /** + * Cells with >= 2 noData ==> noData + */ + if (count > 1) return noDataValue; + + /** + * Now only one cell has noData value. + * Compute interpolation over cell, with vertex layout normalized to put noData in NE. + * This is done by flipping the cell across the X or Y axis, or both + * (and transforming the point offsets likewise) + */ + if (v00 == noDataValue) return interpolateBoundaryCellNorm(1.0f - yfrac, 1.0f - xfrac, v11, v10, v01); + if (v11 == noDataValue) return interpolateBoundaryCellNorm(xfrac, yfrac, v00, v10, v01); + if (v10 == noDataValue) return interpolateBoundaryCellNorm(xfrac, 1.0f - yfrac, v01, v11, v00); + if (v01 == noDataValue) return interpolateBoundaryCellNorm(1.0f - xfrac, yfrac, v10, v00, v11); + + // should never reach here + return noDataValue; + } + + /** + * Computes an interpolated value across a grid cell which has a single + * NO_DATA value in the NE corner (<tt>v11</tt>), and valid data values + * in the three other corners. + * The interpolated value is computed using linear interpolation across + * the 3D triangle defined by the three valid corners. + * + * @param xfrac the fractional x location of the interpolation point within the cell + * @param yfrac the fractional y location of the interpolation point within the cell + * @param v00 the lower left value + * @param v10 the lower right value + * @param v01 the upper left value + * @return the interpolated value + */ + private float interpolateBoundaryCellNorm(float xfrac, float yfrac, float v00, float v10, float v01) + { + // if point is in NE triangle, value is NO_DATA + if (xfrac + yfrac > 1) return noDataValue; + + // interpolate across plane defined by SW triangle and values + float dx = v10 - v00; + float dy = v01 - v00; + float v = v00 + (xfrac * dx) + (yfrac * dy); + return v; + } + +} Added: trunk/modules/unsupported/process-raster/src/main/java/org/geotools/process/raster/surface/GridTransform.java =================================================================== --- trunk/modules/unsupported/process-raster/src/main/java/org/geotools/process/raster/surface/GridTransform.java (rev 0) +++ trunk/modules/unsupported/process-raster/src/main/java/org/geotools/process/raster/surface/GridTransform.java 2012-06-18 18:28:10 UTC (rev 38813) @@ -0,0 +1,115 @@ +/* + * GeoTools - The Open Source Java GIS Toolkit + * http://geotools.org + * + * (C) 2011, Open Source Geospatial Foundation (OSGeo) + * (C) 2008-2011 TOPP - www.openplans.org. + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; + * version 2.1 of the License. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + */ +package org.geotools.process.raster.surface; + +import com.vividsolutions.jts.geom.Envelope; + +/** + * An affine transformation between two parallel + * coordinate systems, one defined by an {@link Envelope} + * and one defined by a discrete zero-based grid + * representing the same area as the envelope. + * The transformation incorporates an isotropic scaling and a translation. + * + * @author Martin Davis - OpenGeo + * + */ +class GridTransform { + + private Envelope env; + + private int xSize; + + private int ySize; + + private double dx; + + private double dy; + + /** + * Creates a new transform. + * + * @param env the envelope defining one coordinate system + * @param xSize the number of cells along the X axis of the grid + * @param ySize the number of cells along the Y axis of the grid + */ + public GridTransform(Envelope env, int xSize, int ySize) { + this.env = env; + this.xSize = xSize; + this.ySize = ySize; + dx = env.getWidth() / (xSize - 1); + dy = env.getHeight() / (ySize - 1); + } + + /** + * Computes the X ordinate of the i'th grid column. + * @param i the index of a grid column + * @return the X ordinate of the column + */ + public double x(int i) { + if (i >= xSize - 1) + return env.getMaxX(); + return env.getMinX() + i * dx; + } + + /** + * Computes the Y ordinate of the i'th grid row. + * @param j the index of a grid row + * @return the Y ordinate of the row + */ + public double y(int j) { + if (j >= ySize - 1) + return env.getMaxY(); + return env.getMinY() + j * dy; + } + + /** + * Computes the column index of an X ordinate. + * @param x the X ordinate + * @return the column index + */ + public int i(double x) { + if (x > env.getMaxX()) + return xSize; + if (x < env.getMinX()) + return -1; + int i = (int) ((x - env.getMinX()) / dx); + // have already check x is in bounds, so ensure returning a valid value + if (i >= xSize) + i = xSize - 1; + return i; + } + + /** + * Computes the column index of an Y ordinate. + * @param y the Y ordinate + * @return the column index + */ + public int j(double y) { + if (y > env.getMaxY()) + return ySize; + if (y < env.getMinY()) + return -1; + int j = (int) ((y - env.getMinY()) / dy); + // have already check x is in bounds, so ensure returning a valid value + if (j >= ySize) + j = ySize - 1; + return j; + } + +} \ No newline at end of file Added: trunk/modules/unsupported/process-raster/src/test/java/org/geotools/process/raster/surface/BarnesSurfaceProcessTest.java =================================================================== --- trunk/modules/unsupported/process-raster/src/test/java/org/geotools/process/raster/surface/BarnesSurfaceProcessTest.java (rev 0) +++ trunk/modules/unsupported/process-raster/src/test/java/org/geotools/process/raster/surface/BarnesSurfaceProcessTest.java 2012-06-18 18:28:10 UTC (rev 38813) @@ -0,0 +1,128 @@ +/* + * GeoTools - The Open Source Java GIS Toolkit + * http://geotools.org + * + * (C) 2002-2011, Open Source Geospatial Foundation (OSGeo) + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; + * version 2.1 of the License. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + */ +package org.geotools.process.raster.surface; + +import static org.junit.Assert.assertTrue; + +import java.awt.geom.Point2D; + +import org.geotools.coverage.grid.GridCoverage2D; +import org.geotools.data.simple.SimpleFeatureCollection; +import org.geotools.feature.FeatureCollections; +import org.geotools.feature.simple.SimpleFeatureBuilder; +import org.geotools.feature.simple.SimpleFeatureTypeBuilder; +import org.geotools.geometry.jts.ReferencedEnvelope; +import org.geotools.referencing.crs.DefaultGeographicCRS; +import org.junit.Test; +import org.opengis.feature.simple.SimpleFeatureType; +import org.opengis.util.ProgressListener; + +import com.vividsolutions.jts.geom.Coordinate; +import com.vividsolutions.jts.ge... [truncated message content] |
