Menu
▾
▴
[Geotools-commits] svn - r38815 - in trunk/modules/unsupported/process-feature/src: main/java/org/geotools/process/feature/gs test/java/org/geotools/process/feature/gs
|
From: <svn...@os...> - 2012-06-18 22:28:56
|
Author: mdavis Date: 2012-06-18 15:28:48 -0700 (Mon, 18 Jun 2012) New Revision: 38815 Added: trunk/modules/unsupported/process-feature/src/main/java/org/geotools/process/feature/gs/PointStackerProcess.java trunk/modules/unsupported/process-feature/src/test/java/org/geotools/process/feature/gs/PointStackerProcessTest.java Log: Added PointStacker Rendering Transformation process Added: trunk/modules/unsupported/process-feature/src/main/java/org/geotools/process/feature/gs/PointStackerProcess.java =================================================================== --- trunk/modules/unsupported/process-feature/src/main/java/org/geotools/process/feature/gs/PointStackerProcess.java (rev 0) +++ trunk/modules/unsupported/process-feature/src/main/java/org/geotools/process/feature/gs/PointStackerProcess.java 2012-06-18 22:28:48 UTC (rev 38815) @@ -0,0 +1,384 @@ +/* + * GeoTools - The Open Source Java GIS Toolkit + * http://geotools.org + * + * (C) 2011, Open Source Geospatial Foundation (OSGeo) + * (C) 2001-2007 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.feature.gs; + +import java.util.Collection; +import java.util.HashMap; +import java.util.HashSet; +import java.util.Map; +import java.util.Set; + +import org.geotools.data.collection.ListFeatureCollection; +import org.geotools.data.simple.SimpleFeatureCollection; +import org.geotools.data.simple.SimpleFeatureIterator; +import org.geotools.feature.simple.SimpleFeatureBuilder; +import org.geotools.feature.simple.SimpleFeatureTypeBuilder; +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.feature.simple.SimpleFeature; +import org.opengis.feature.simple.SimpleFeatureType; +import org.opengis.referencing.FactoryException; +import org.opengis.referencing.crs.CoordinateReferenceSystem; +import org.opengis.referencing.operation.MathTransform; +import org.opengis.referencing.operation.TransformException; +import org.opengis.util.ProgressListener; + +import com.vividsolutions.jts.geom.Coordinate; +import com.vividsolutions.jts.geom.Geometry; +import com.vividsolutions.jts.geom.GeometryFactory; +import com.vividsolutions.jts.geom.Point; +import com.vividsolutions.jts.geom.impl.PackedCoordinateSequenceFactory; + +/** + * A Rendering Transformation process which aggregates features into a set of + * visually non-conflicting point features. + * The created points have attributes which provide the total number of points + * aggregated, as well as the number of unique point locations. + * <p> + * This is sometimes called "point clustering". + * The term stacking is used instead, since clustering has multiple + * meanings in geospatial processing - it is also used to + * mean identifying groups defined by point proximity. + * <p> + * The stacking is defined by specifying a grid to aggregate to. + * The grid cell size is specified in pixels relative to the requested output image size. + * This makes it more intuitive to pick an appropriate grid size, + * and ensures that the aggregation works at all zoom levels. + * <p> + * The output is a FeatureCollection containing the following attributes: + * <ul> + * <li><code>geom</code> - the point representing the cluster + * <li><code>count</code> - the total number of points in the cluster + * <li><code>countunique</code> - the number of unique point locations in the cluster + * </ul> + * Note that as required by the Rendering Transformation API, the output + * has the CRS of the input data. + * + * @author mdavis + * + */ +@DescribeProcess(title = "PointStacker", description = "Aggregates a collection of points into a set of stacked points.") +public class PointStackerProcess implements GSProcess { + + public static final String ATTR_GEOM = "geom"; + public static final String ATTR_COUNT = "count"; + public static final String ATTR_COUNT_UNIQUE = "countunique"; + + //TODO: add ability to pick index point selection strategy + //TODO: add ability to set attribute name containing value to be aggregated + //TODO: add ability to specify aggregation method (COUNT, SUM, AVG) + //TODO: ultimately could allow aggregating multiple input attributes, with different methods for each + //TODO: allow including attributes from input data (eg for use with points that are not aggregated) + //TODO: expand query window to avoid edge effects? + + // no process state is defined, since RenderingTransformation processes must be stateless + + @DescribeResult(name = "result", description = "The collection of stacked points") + public SimpleFeatureCollection execute( + + // process data + @DescribeParameter(name = "data", description = "Features containing the data points") SimpleFeatureCollection data, + + // process parameters + @DescribeParameter(name = "cellSize", description = "Cell size for gridding, in pixels") Integer cellSize, + + // output image parameters + @DescribeParameter(name = "outputBBOX", description = "Georeferenced bounding box of the output image") ReferencedEnvelope outputEnv, + @DescribeParameter(name = "outputWidth", description = "Width of the output image, in pixels") Integer outputWidth, + @DescribeParameter(name = "outputHeight", description = "Height of the output image, in pixels") Integer outputHeight, + + ProgressListener monitor) throws ProcessException, TransformException { + + CoordinateReferenceSystem srcCRS = data.getSchema().getCoordinateReferenceSystem(); + CoordinateReferenceSystem dstCRS = outputEnv.getCoordinateReferenceSystem(); + MathTransform crsTransform = null; + MathTransform invTransform = null; + try { + crsTransform = CRS.findMathTransform(srcCRS, dstCRS); + invTransform = crsTransform.inverse(); + } catch (FactoryException e) { + throw new ProcessException(e); + } + + // TODO: allow output CRS to be different to data CRS + // assume same CRS for now... + double cellSizeSrc = cellSize * outputEnv.getWidth() / outputWidth; + + Collection<StackedPoint> stackedPts = stackPoints(data, crsTransform, cellSizeSrc, + outputEnv.getMinX(), outputEnv.getMinY()); + + SimpleFeatureType schema = createType(srcCRS); + SimpleFeatureCollection result = new ListFeatureCollection(schema); + SimpleFeatureBuilder fb = new SimpleFeatureBuilder(schema); + + GeometryFactory factory = new GeometryFactory(new PackedCoordinateSequenceFactory()); + + double[] srcPt = new double[2]; + double[] dstPt = new double[2]; + + + for (StackedPoint sp : stackedPts) { + // create feature for stacked point + Coordinate pt = sp.getLocation(); + + // transform back to src CRS, since RT rendering expects the output to be in the same CRS + srcPt[0] = pt.x; + srcPt[1] = pt.y; + invTransform.transform(srcPt, 0, dstPt, 0, 1); + Coordinate psrc = new Coordinate(dstPt[0], dstPt[1]); + + Geometry point = factory.createPoint(psrc); + fb.add(point); + fb.add(sp.getCount()); + fb.add(sp.getCountUnique()); + + result.add(fb.buildFeature(null)); + } + return result; + } + + /** + * Computes the stacked points for the given data collection. + * All geometry types are handled - for non-point geometries, the centroid is used. + * + * @param data + * @param cellSize + * @param minX + * @param minY + * @return + * @throws TransformException + */ + private Collection<StackedPoint> stackPoints(SimpleFeatureCollection data, + MathTransform crsTransform, + double cellSize, double minX, double minY) throws TransformException { + SimpleFeatureIterator featureIt = data.features(); + + Map<Coordinate, StackedPoint> stackedPts = new HashMap<Coordinate, StackedPoint>(); + + double[] srcPt = new double[2]; + double[] dstPt = new double[2]; + + Coordinate indexPt = new Coordinate(); + try { + while (featureIt.hasNext()) { + SimpleFeature feature = featureIt.next(); + // get the point location from the geometry + Geometry geom = (Geometry) feature.getDefaultGeometry(); + Coordinate p = getRepresentativePoint(geom); + + // reproject data point to output CRS, if required + srcPt[0] = p.x; + srcPt[1] = p.y; + crsTransform.transform(srcPt, 0, dstPt, 0, 1); + Coordinate pout = new Coordinate(dstPt[0], dstPt[1]); + + indexPt.x = pout.x; + indexPt.y = pout.y; + gridIndex(indexPt, cellSize); + + StackedPoint stkPt = stackedPts.get(indexPt); + if (stkPt == null) { + + /** + * Note that the + */ + double centreX = indexPt.x * cellSize + cellSize / 2; + double centreY = indexPt.y * cellSize + cellSize / 2; + + stkPt = new StackedPoint(indexPt, new Coordinate(centreX, centreY)); + stackedPts.put(stkPt.getKey(), stkPt); + } + stkPt.add(pout); + } + + } finally { + featureIt.close(); + } + return stackedPts.values(); + } + + /** + * Gets a point to represent the Geometry. + * If the Geometry is a point, this is returned. + * Otherwise, the centroid is used. + * + * @param g the geometry to find a point for + * @return a point representing the Geometry + */ + private static Coordinate getRepresentativePoint(Geometry g) + { + if (g.getNumPoints() == 1) + return g.getCoordinate(); + return g.getCentroid().getCoordinate(); + } + + /** + * Computes the grid index for a point for the grid determined by the cellsize. + * + * @param griddedPt the point to grid, and also holds the output value + * @param cellSize the grid cell size + */ + private void gridIndex(Coordinate griddedPt, double cellSize) { + + // TODO: is there any situation where this could result in too much loss of precision? + /** + * The grid is based at the origin of the entire data space, + * not just the query window. + * This makes gridding stable during panning. + * + * This should not lose too much precision for any reasonable coordinate system and map size. + * The worst case is a CRS with small ordinate values, and a large cell size. + * The worst case tested is a map in degrees, zoomed out to show about twice the globe - works fine. + */ + // Use longs to avoid possible overflow issues (e.g. for a very small cell size) + long ix = (long) ((griddedPt.x) / cellSize); + long iy = (long) ((griddedPt.y) / cellSize); + + griddedPt.x = ix; + griddedPt.y = iy; + } + + private SimpleFeatureType createType(CoordinateReferenceSystem crs) { + SimpleFeatureTypeBuilder tb = new SimpleFeatureTypeBuilder(); + tb.add(ATTR_GEOM, Point.class, crs); + tb.add(ATTR_COUNT, Integer.class); + tb.add(ATTR_COUNT_UNIQUE, Integer.class); + tb.setName("stackedPoint"); + SimpleFeatureType sfType = tb.buildFeatureType(); + return sfType; + } + + private static class StackedPoint { + private Coordinate key; + + private Coordinate centerPt; + + private Coordinate location = null; + + private int count = 0; + + private Set<Coordinate> uniquePts; + + /** + * Creates a new stacked point grid cell. + * The center point of the cell is supplied + * so that it may be used as or influence the + * location of the final display point + * + * @param key a key for the grid cell (using integer ordinates to avoid precision issues) + * @param centerPt the center point of the grid cell + */ + public StackedPoint(Coordinate key, Coordinate centerPt) { + this.key = new Coordinate(key); + this.centerPt = centerPt; + } + + public Coordinate getKey() { + return key; + } + + public Coordinate getLocation() { + return location; + } + + public int getCount() { + return count; + } + + public int getCountUnique() { + if (uniquePts == null) + return 1; + return uniquePts.size(); + } + + public void add(Coordinate pt) { + count++; + /** + * Only create set if this is the second point seen + * (and assum the first pt is in location) + */ + if (uniquePts == null) { + uniquePts = new HashSet<Coordinate>(); + } + uniquePts.add(pt); + + pickNearestLocation(pt); + //pickCenterLocation(pt); + } + + /** + * Picks the location as the point + * which is nearest to the center of the cell. + * In addition, the nearest location is averaged with the cell center. + * This gives the best chance of avoiding conflicts. + * + * @param pt + */ + private void pickNearestLocation(Coordinate pt) { + // strategy - pick most central point + if (location == null) { + location = average(centerPt, pt); + return; + } + if (pt.distance(centerPt) < location.distance(centerPt)) { + location = average(centerPt, pt); + } + } + + /** + * Picks the location as the centre point of the cell. + * This does not give a good visualization - the gridding is very obvious + * + * @param pt + */ + private void pickCenterLocation(Coordinate pt) { + // strategy - pick first point + if (location == null) { + location = new Coordinate(pt); + return; + } + location = centerPt; + } + + /** + * Picks the first location encountered as the cell location. + * This is sub-optimal, since if the first point is near the cell + * boundary it is likely to collide with neighbouring points. + * + * @param pt + */ + private void pickFirstLocation(Coordinate pt) { + // strategy - pick first point + if (location == null) { + location = new Coordinate(pt); + } + } + + private static Coordinate average(Coordinate p1, Coordinate p2) + { + double x = (p1.x + p2.x) / 2; + double y = (p1.y + p2.y) / 2; + return new Coordinate(x, y); + } + } +} Added: trunk/modules/unsupported/process-feature/src/test/java/org/geotools/process/feature/gs/PointStackerProcessTest.java =================================================================== --- trunk/modules/unsupported/process-feature/src/test/java/org/geotools/process/feature/gs/PointStackerProcessTest.java (rev 0) +++ trunk/modules/unsupported/process-feature/src/test/java/org/geotools/process/feature/gs/PointStackerProcessTest.java 2012-06-18 22:28:48 UTC (rev 38815) @@ -0,0 +1,198 @@ +/* + * 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.feature.gs; + +import static junit.framework.Assert.*; + +import org.geotools.data.simple.SimpleFeatureCollection; +import org.geotools.data.simple.SimpleFeatureIterator; +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.process.ProcessException; +import org.geotools.referencing.CRS; +import org.geotools.referencing.crs.DefaultGeographicCRS; +import org.junit.Test; +import org.opengis.feature.simple.SimpleFeature; +import org.opengis.feature.simple.SimpleFeatureType; +import org.opengis.referencing.FactoryException; +import org.opengis.referencing.NoSuchAuthorityCodeException; +import org.opengis.referencing.crs.CoordinateReferenceSystem; +import org.opengis.referencing.operation.TransformException; +import org.opengis.util.ProgressListener; + +import com.vividsolutions.jts.geom.Coordinate; +import com.vividsolutions.jts.geom.Geometry; +import com.vividsolutions.jts.geom.GeometryFactory; +import com.vividsolutions.jts.geom.MultiPoint; +import com.vividsolutions.jts.geom.Point; +import com.vividsolutions.jts.geom.impl.PackedCoordinateSequenceFactory; + +/** + * Unit test for PointStackerProcess. + * + * @author Martin Davis, OpenGeo + * + */ +public class PointStackerProcessTest { + @Test + public void testSimple() throws ProcessException, TransformException { + ReferencedEnvelope bounds = new ReferencedEnvelope(0, 10, 0, 10, DefaultGeographicCRS.WGS84); + + // Simple dataset with some coincident points + Coordinate[] data = new Coordinate[] { new Coordinate(4, 4), new Coordinate(4.1, 4.1), + new Coordinate(4.1, 4.1), new Coordinate(8, 8) }; + + + SimpleFeatureCollection fc = createPoints(data, bounds); + ProgressListener monitor = null; + + PointStackerProcess psp = new PointStackerProcess(); + SimpleFeatureCollection result = psp.execute(fc, 100, // cellSize + bounds, // outputBBOX + 1000, // outputWidth + 1000, // outputHeight + monitor); + + checkSchemaCorrect(result.getSchema()); + assertEquals(2, result.size()); + checkResultPoint(result, new Coordinate(4, 4), 3, 2); + checkResultPoint(result, new Coordinate(8, 8), 1, 1); + } + + /** + * Tests point stacking when output CRS is different to data CRS. + * The result data should be reprojected. + * + * @throws NoSuchAuthorityCodeException + * @throws FactoryException + * @throws TransformException + * @throws ProcessException + */ + @Test + public void testReprojected() throws NoSuchAuthorityCodeException, FactoryException, ProcessException, TransformException { + + ReferencedEnvelope inBounds = new ReferencedEnvelope(0, 10, 0, 10, DefaultGeographicCRS.WGS84); + + // Dataset with some points located in appropriate area + // points are close enough to create a single cluster + Coordinate[] data = new Coordinate[] { new Coordinate(-121.813201, 48.777343), new Coordinate(-121.813, 48.777) }; + + + SimpleFeatureCollection fc = createPoints(data, inBounds); + ProgressListener monitor = null; + + // Google Mercator BBOX for northern Washington State (roughly) + CoordinateReferenceSystem webMerc = CRS.decode("EPSG:3785"); + ReferencedEnvelope outBounds = new ReferencedEnvelope(-1.4045034049133E7, -1.2937920131607E7, 5916835.1504419, 6386464.2521607, webMerc); + + PointStackerProcess psp = new PointStackerProcess(); + SimpleFeatureCollection result = psp.execute(fc, 100, // cellSize + outBounds, // outputBBOX + 1810, // outputWidth + 768, // outputHeight + monitor); + + checkSchemaCorrect(result.getSchema()); + assertEquals(1, result.size()); + assertEquals(inBounds.getCoordinateReferenceSystem(), result.getBounds().getCoordinateReferenceSystem()); + checkResultPoint(result, new Coordinate(-121.813201, 48.777343), 2, 2); + } + + /** + * Check that a result set contains a stacked point in the right cell with expected attribute + * values. Because it's not known in advance what the actual location of a stacked point will + * be, a nearest-point strategy is used. + * + * @param result + * @param coordinate + * @param i + * @param j + */ + private void checkResultPoint(SimpleFeatureCollection result, Coordinate testPt, + int expectedCount, int expectedCountUnique) { + /** + * Find closest point to loc pt, then check that the attributes match + */ + double minDist = Double.MAX_VALUE; + int count = -1; + int countunique = -1; + + // find nearest result to testPt + for (SimpleFeatureIterator it = result.features(); it.hasNext();) { + SimpleFeature f = it.next(); + Coordinate outPt = ((Point) f.getDefaultGeometry()).getCoordinate(); + double dist = outPt.distance(testPt); + if (dist < minDist) { + minDist = dist; + count = (Integer) f.getAttribute(PointStackerProcess.ATTR_COUNT); + countunique = (Integer) f.getAttribute(PointStackerProcess.ATTR_COUNT_UNIQUE); + } + } + assertEquals(expectedCount, count); + assertEquals(expectedCountUnique, countunique); + } + + private void checkResultBounds(SimpleFeatureCollection result, ReferencedEnvelope bounds) { + boolean isInBounds = true; + for (SimpleFeatureIterator it = result.features(); it.hasNext();) { + SimpleFeature f = it.next(); + Coordinate outPt = ((Point) f.getDefaultGeometry()).getCoordinate(); + if (! bounds.contains(outPt)) { + isInBounds = false; + System.out.println("Found point out of bounds: " + f.getDefaultGeometry()); + } + } + assertTrue(isInBounds); + } + + private void checkSchemaCorrect(SimpleFeatureType ft) { + assertEquals(3, ft.getAttributeCount()); + assertEquals(Point.class, ft.getGeometryDescriptor().getType().getBinding()); + assertEquals(Integer.class, ft.getDescriptor(PointStackerProcess.ATTR_COUNT).getType() + .getBinding()); + assertEquals(Integer.class, ft.getDescriptor(PointStackerProcess.ATTR_COUNT_UNIQUE) + .getType().getBinding()); + + } + + private SimpleFeatureCollection createPoints(Coordinate[] pts, ReferencedEnvelope bounds) { + + SimpleFeatureTypeBuilder tb = new SimpleFeatureTypeBuilder(); + tb.setName("data"); + tb.setCRS(bounds.getCoordinateReferenceSystem()); + tb.add("shape", MultiPoint.class); + tb.add("value", Double.class); + + SimpleFeatureType type = tb.buildFeatureType(); + SimpleFeatureBuilder fb = new SimpleFeatureBuilder(type); + SimpleFeatureCollection fc = FeatureCollections.newCollection(); + + GeometryFactory factory = new GeometryFactory(new PackedCoordinateSequenceFactory()); + + for (Coordinate p : pts) { + Geometry point = factory.createPoint(p); + fb.add(point); + fb.add(p.z); + fc.add(fb.buildFeature(null)); + } + + return fc; + } + +} |
