diff --git a/pom.xml b/pom.xml index 742117eb9a43b..4cee753c752a4 100644 --- a/pom.xml +++ b/pom.xml @@ -1145,7 +1145,7 @@ com.esri.geometry esri-geometry-api - 2.2.2 + 2.2.3 com.fasterxml.jackson.core diff --git a/presto-docs/src/main/sphinx/functions/geospatial.rst b/presto-docs/src/main/sphinx/functions/geospatial.rst index 81fbb4ced313d..68fa934250655 100644 --- a/presto-docs/src/main/sphinx/functions/geospatial.rst +++ b/presto-docs/src/main/sphinx/functions/geospatial.rst @@ -186,7 +186,9 @@ Operations .. function:: ST_Buffer(Geometry, distance) -> Geometry Returns the geometry that represents all points whose distance from the - specified geometry is less than or equal to the specified distance. + specified geometry is less than or equal to the specified distance. If the + points of the geometry are extremely close together (``delta < 1e-8``), this + might return an empty geometry. .. function:: ST_Difference(Geometry, Geometry) -> Geometry diff --git a/presto-geospatial/src/main/java/com/facebook/presto/plugin/geospatial/GeoFunctions.java b/presto-geospatial/src/main/java/com/facebook/presto/plugin/geospatial/GeoFunctions.java index df4b74a224018..51d8d5af4654a 100644 --- a/presto-geospatial/src/main/java/com/facebook/presto/plugin/geospatial/GeoFunctions.java +++ b/presto-geospatial/src/main/java/com/facebook/presto/plugin/geospatial/GeoFunctions.java @@ -30,7 +30,6 @@ import com.esri.core.geometry.ogc.OGCGeometry; import com.esri.core.geometry.ogc.OGCGeometryCollection; import com.esri.core.geometry.ogc.OGCLineString; -import com.esri.core.geometry.ogc.OGCMultiPolygon; import com.esri.core.geometry.ogc.OGCPoint; import com.esri.core.geometry.ogc.OGCPolygon; import com.facebook.presto.geospatial.GeometryType; @@ -384,33 +383,7 @@ public static Slice stCentroid(@SqlType(GEOMETRY_TYPE_NAME) Slice input) return serialize(createFromEsriGeometry(new Point(), geometry.getEsriSpatialReference())); } - Point centroid; - try { - switch (geometryType) { - case MULTI_POINT: - centroid = computePointsCentroid((MultiVertexGeometry) geometry.getEsriGeometry()); - break; - case LINE_STRING: - case MULTI_LINE_STRING: - centroid = computeLineCentroid((Polyline) geometry.getEsriGeometry()); - break; - case POLYGON: - centroid = computePolygonCentroid((Polygon) geometry.getEsriGeometry()); - break; - case MULTI_POLYGON: - centroid = computeMultiPolygonCentroid((OGCMultiPolygon) geometry); - break; - default: - throw new PrestoException(INVALID_FUNCTION_ARGUMENT, "Unexpected geometry type: " + geometryType); - } - } - catch (RuntimeException e) { - if (e instanceof PrestoException && ((PrestoException) e).getErrorCode() == INVALID_FUNCTION_ARGUMENT.toErrorCode()) { - throw e; - } - throw new PrestoException(INVALID_FUNCTION_ARGUMENT, format("Cannot compute centroid: %s Use ST_IsValid to confirm that input geometry is valid or compute centroid for a bounding box using ST_Envelope.", e.getMessage()), e); - } - return serialize(createFromEsriGeometry(centroid, geometry.getEsriSpatialReference())); + return serialize(geometry.centroid()); } @Description("Returns the minimum convex geometry that encloses all input geometries") @@ -1360,119 +1333,6 @@ private static void verifySameSpatialReference(OGCGeometry leftGeometry, OGCGeom checkArgument(Objects.equals(leftGeometry.getEsriSpatialReference(), rightGeometry.getEsriSpatialReference()), "Input geometries must have the same spatial reference"); } - // Points centroid is arithmetic mean of the input points - private static Point computePointsCentroid(MultiVertexGeometry multiVertex) - { - double xSum = 0; - double ySum = 0; - for (int i = 0; i < multiVertex.getPointCount(); i++) { - Point point = multiVertex.getPoint(i); - xSum += point.getX(); - ySum += point.getY(); - } - return new Point(xSum / multiVertex.getPointCount(), ySum / multiVertex.getPointCount()); - } - - // Lines centroid is weighted mean of each line segment, weight in terms of line length - private static Point computeLineCentroid(Polyline polyline) - { - double xSum = 0; - double ySum = 0; - double weightSum = 0; - for (int i = 0; i < polyline.getPathCount(); i++) { - Point startPoint = polyline.getPoint(polyline.getPathStart(i)); - Point endPoint = polyline.getPoint(polyline.getPathEnd(i) - 1); - double dx = endPoint.getX() - startPoint.getX(); - double dy = endPoint.getY() - startPoint.getY(); - double length = sqrt(dx * dx + dy * dy); - weightSum += length; - xSum += (startPoint.getX() + endPoint.getX()) * length / 2; - ySum += (startPoint.getY() + endPoint.getY()) * length / 2; - } - return new Point(xSum / weightSum, ySum / weightSum); - } - - // Polygon centroid: area weighted average of centroids in case of holes - private static Point computePolygonCentroid(Polygon polygon) - { - int pathCount = polygon.getPathCount(); - - if (pathCount == 1) { - return getPolygonSansHolesCentroid(polygon); - } - - double xSum = 0; - double ySum = 0; - double areaSum = 0; - - for (int i = 0; i < pathCount; i++) { - int startIndex = polygon.getPathStart(i); - int endIndex = polygon.getPathEnd(i); - - Polygon sansHoles = getSubPolygon(polygon, startIndex, endIndex); - - Point centroid = getPolygonSansHolesCentroid(sansHoles); - double area = sansHoles.calculateArea2D(); - - xSum += centroid.getX() * area; - ySum += centroid.getY() * area; - areaSum += area; - } - - return new Point(xSum / areaSum, ySum / areaSum); - } - - private static Polygon getSubPolygon(Polygon polygon, int startIndex, int endIndex) - { - Polyline boundary = new Polyline(); - boundary.startPath(polygon.getPoint(startIndex)); - for (int i = startIndex + 1; i < endIndex; i++) { - Point current = polygon.getPoint(i); - boundary.lineTo(current); - } - - final Polygon newPolygon = new Polygon(); - newPolygon.add(boundary, false); - return newPolygon; - } - - // Polygon sans holes centroid: - // c[x] = (Sigma(x[i] + x[i + 1]) * (x[i] * y[i + 1] - x[i + 1] * y[i]), for i = 0 to N - 1) / (6 * signedArea) - // c[y] = (Sigma(y[i] + y[i + 1]) * (x[i] * y[i + 1] - x[i + 1] * y[i]), for i = 0 to N - 1) / (6 * signedArea) - private static Point getPolygonSansHolesCentroid(Polygon polygon) - { - int pointCount = polygon.getPointCount(); - double xSum = 0; - double ySum = 0; - double signedArea = 0; - for (int i = 0; i < pointCount; i++) { - Point current = polygon.getPoint(i); - Point next = polygon.getPoint((i + 1) % polygon.getPointCount()); - double ladder = current.getX() * next.getY() - next.getX() * current.getY(); - xSum += (current.getX() + next.getX()) * ladder; - ySum += (current.getY() + next.getY()) * ladder; - signedArea += ladder / 2; - } - return new Point(xSum / (signedArea * 6), ySum / (signedArea * 6)); - } - - // MultiPolygon centroid is weighted mean of each polygon, weight in terms of polygon area - private static Point computeMultiPolygonCentroid(OGCMultiPolygon multiPolygon) - { - double xSum = 0; - double ySum = 0; - double weightSum = 0; - for (int i = 0; i < multiPolygon.numGeometries(); i++) { - Point centroid = computePolygonCentroid((Polygon) multiPolygon.geometryN(i).getEsriGeometry()); - Polygon polygon = (Polygon) multiPolygon.geometryN(i).getEsriGeometry(); - double weight = polygon.calculateArea2D(); - weightSum += weight; - xSum += centroid.getX() * weight; - ySum += centroid.getY() * weight; - } - return new Point(xSum / weightSum, ySum / weightSum); - } - private static boolean envelopes(Slice left, Slice right, EnvelopesPredicate predicate) { Envelope leftEnvelope = deserializeEnvelope(left); diff --git a/presto-geospatial/src/test/java/com/facebook/presto/plugin/geospatial/TestBingTileFunctions.java b/presto-geospatial/src/test/java/com/facebook/presto/plugin/geospatial/TestBingTileFunctions.java index dbb8df50f94bc..977033be8657d 100644 --- a/presto-geospatial/src/test/java/com/facebook/presto/plugin/geospatial/TestBingTileFunctions.java +++ b/presto-geospatial/src/test/java/com/facebook/presto/plugin/geospatial/TestBingTileFunctions.java @@ -382,7 +382,7 @@ public void testBingTileZoomLevel() public void testBingTilePolygon() { assertFunction("ST_AsText(bing_tile_polygon(bing_tile('123030123010121')))", VARCHAR, "POLYGON ((59.996337890625 30.11662158281937, 60.00732421875 30.11662158281937, 60.00732421875 30.12612436422458, 59.996337890625 30.12612436422458, 59.996337890625 30.11662158281937))"); - assertFunction("ST_AsText(ST_Centroid(bing_tile_polygon(bing_tile('123030123010121'))))", VARCHAR, "POINT (60.0018310442288 30.121372968273892)"); + assertFunction("ST_AsText(ST_Centroid(bing_tile_polygon(bing_tile('123030123010121'))))", VARCHAR, "POINT (60.0018310546875 30.121372973521975)"); // Check bottom right corner of a stack of tiles at different zoom levels assertFunction("ST_AsText(apply(bing_tile_polygon(bing_tile(1, 1, 1)), g -> ST_Point(ST_XMax(g), ST_YMin(g))))", VARCHAR, "POINT (180 -85.05112877980659)"); diff --git a/presto-geospatial/src/test/java/com/facebook/presto/plugin/geospatial/TestGeoFunctions.java b/presto-geospatial/src/test/java/com/facebook/presto/plugin/geospatial/TestGeoFunctions.java index 1901faba6c73e..230efc131fb9c 100644 --- a/presto-geospatial/src/test/java/com/facebook/presto/plugin/geospatial/TestGeoFunctions.java +++ b/presto-geospatial/src/test/java/com/facebook/presto/plugin/geospatial/TestGeoFunctions.java @@ -14,9 +14,11 @@ package com.facebook.presto.plugin.geospatial; import com.esri.core.geometry.Point; +import com.esri.core.geometry.ogc.OGCGeometry; import com.esri.core.geometry.ogc.OGCPoint; import com.facebook.presto.geospatial.KdbTreeUtils; import com.facebook.presto.geospatial.Rectangle; +import com.facebook.presto.geospatial.serde.GeometrySerde; import com.facebook.presto.operator.scalar.AbstractTestFunctions; import com.facebook.presto.spi.block.Block; import com.facebook.presto.spi.block.BlockBuilder; @@ -30,6 +32,7 @@ import java.util.stream.Collectors; import static com.facebook.presto.geospatial.KdbTree.buildKdbTree; +import static com.facebook.presto.plugin.geospatial.GeoFunctions.stCentroid; import static com.facebook.presto.plugin.geospatial.GeometryType.GEOMETRY; import static com.facebook.presto.spi.StandardErrorCode.INVALID_FUNCTION_ARGUMENT; import static com.facebook.presto.spi.type.BigintType.BIGINT; @@ -193,6 +196,16 @@ public void testSTBuffer() // infinity() and nan() distance assertFunction("ST_AsText(ST_Buffer(ST_Point(0, 0), infinity()))", VARCHAR, "MULTIPOLYGON EMPTY"); assertInvalidFunction("ST_Buffer(ST_Point(0, 0), nan())", "distance is NaN"); + + // For small polygons, there was a bug in ESRI that throw an NPE. This + // was fixed (https://github.com/Esri/geometry-api-java/pull/243) to + // return an empty geometry instead. Ideally, these would return + // something approximately like `ST_Buffer(ST_Centroid(geometry))`. + assertFunction("ST_IsEmpty(ST_Buffer(ST_Buffer(ST_Point(177.50102959662, 64.726807421691), 0.0000000001), 0.00005))", + BOOLEAN, true); + assertFunction("ST_IsEmpty(ST_Buffer(ST_GeometryFromText(" + + "'POLYGON ((177.0 64.0, 177.0000000001 64.0, 177.0000000001 64.0000000001, 177.0 64.0000000001, 177.0 64.0))'" + + "), 0.01))", BOOLEAN, true); } @Test @@ -209,7 +222,24 @@ public void testSTCentroid() assertCentroid("POLYGON ((0 0, 0 5, 5 5, 5 0, 0 0), (1 1, 1 2, 2 2, 2 1, 1 1))", new Point(2.5416666666666665, 2.5416666666666665)); // invalid geometry - assertInvalidFunction("ST_Centroid(ST_GeometryFromText('MULTIPOLYGON (((4.903234300000006 52.08474289999999, 4.903234265193165 52.084742934806826, 4.903234299999999 52.08474289999999, 4.903234300000006 52.08474289999999)))'))", "Cannot compute centroid: .* Use ST_IsValid to confirm that input geometry is valid or compute centroid for a bounding box using ST_Envelope."); + assertApproximateCentroid("MULTIPOLYGON (((4.903234300000006 52.08474289999999, 4.903234265193165 52.084742934806826, 4.903234299999999 52.08474289999999, 4.903234300000006 52.08474289999999)))", new Point(4.9032343, 52.0847429), 1e-7); + + // Numerical stability tests + assertApproximateCentroid( + "MULTIPOLYGON (((153.492818 -28.13729, 153.492821 -28.137291, 153.492816 -28.137289, 153.492818 -28.13729)))", + new Point(153.49282, -28.13729), 1e-5); + assertApproximateCentroid( + "MULTIPOLYGON (((153.112475 -28.360526, 153.1124759 -28.360527, 153.1124759 -28.360526, 153.112475 -28.360526)))", + new Point(153.112475, -28.360526), 1e-5); + assertApproximateCentroid( + "POLYGON ((4.903234300000006 52.08474289999999, 4.903234265193165 52.084742934806826, 4.903234299999999 52.08474289999999, 4.903234300000006 52.08474289999999))", + new Point(4.9032343, 52.0847429), 1e-6); + assertApproximateCentroid( + "MULTIPOLYGON (((4.903234300000006 52.08474289999999, 4.903234265193165 52.084742934806826, 4.903234299999999 52.08474289999999, 4.903234300000006 52.08474289999999)))", + new Point(4.9032343, 52.0847429), 1e-6); + assertApproximateCentroid( + "POLYGON ((-81.0387349 29.20822, -81.039974 29.210597, -81.0410331 29.2101579, -81.0404758 29.2090879, -81.0404618 29.2090609, -81.040433 29.209005, -81.0404269 29.208993, -81.0404161 29.2089729, -81.0398001 29.20779, -81.0387349 29.20822), (-81.0404229 29.208986, -81.04042 29.2089809, -81.0404269 29.208993, -81.0404229 29.208986))", + new Point(-81.039885, 29.209191), 1e-6); } private void assertCentroid(String wkt, Point centroid) @@ -217,6 +247,14 @@ private void assertCentroid(String wkt, Point centroid) assertFunction(format("ST_AsText(ST_Centroid(ST_GeometryFromText('%s')))", wkt), VARCHAR, new OGCPoint(centroid, null).asText()); } + private void assertApproximateCentroid(String wkt, Point expectedCentroid, double epsilon) + { + OGCPoint actualCentroid = (OGCPoint) GeometrySerde.deserialize( + stCentroid(GeometrySerde.serialize(OGCGeometry.fromText(wkt)))); + assertEquals(actualCentroid.X(), expectedCentroid.getX(), epsilon); + assertEquals(actualCentroid.Y(), expectedCentroid.getY(), epsilon); + } + @Test public void testSTConvexHull() {