diff --git a/docs/src/main/sphinx/functions/geospatial.md b/docs/src/main/sphinx/functions/geospatial.md index 448bc8f2e202..19473464ada1 100644 --- a/docs/src/main/sphinx/functions/geospatial.md +++ b/docs/src/main/sphinx/functions/geospatial.md @@ -174,7 +174,8 @@ Returns the closure of the combinatorial boundary of this geometry. :::{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. +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/plugin/trino-geospatial/src/main/java/io/trino/plugin/geospatial/GeoFunctions.java b/plugin/trino-geospatial/src/main/java/io/trino/plugin/geospatial/GeoFunctions.java index cd36a36e0044..74a7fad720cf 100644 --- a/plugin/trino-geospatial/src/main/java/io/trino/plugin/geospatial/GeoFunctions.java +++ b/plugin/trino-geospatial/src/main/java/io/trino/plugin/geospatial/GeoFunctions.java @@ -33,7 +33,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.google.common.base.Joiner; @@ -433,25 +432,7 @@ public static Slice stCentroid(@SqlType(GEOMETRY_TYPE_NAME) Slice input) return serialize(createFromEsriGeometry(new Point(), geometry.getEsriSpatialReference())); } - Point centroid; - 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 TrinoException(INVALID_FUNCTION_ARGUMENT, "Unexpected geometry type: " + geometryType); - } - return serialize(createFromEsriGeometry(centroid, geometry.getEsriSpatialReference())); + return serialize(geometry.centroid()); } @Description("Returns the minimum convex geometry that encloses all input geometries") @@ -1609,119 +1590,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); - } - - 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/plugin/trino-geospatial/src/test/java/io/trino/plugin/geospatial/TestBingTileFunctions.java b/plugin/trino-geospatial/src/test/java/io/trino/plugin/geospatial/TestBingTileFunctions.java index e38fe7fde072..4044d3610577 100644 --- a/plugin/trino-geospatial/src/test/java/io/trino/plugin/geospatial/TestBingTileFunctions.java +++ b/plugin/trino-geospatial/src/test/java/io/trino/plugin/geospatial/TestBingTileFunctions.java @@ -454,7 +454,7 @@ public void testBingTilePolygon() assertThat(assertions.function("ST_AsText", "ST_Centroid(bing_tile_polygon(bing_tile('123030123010121')))")) .hasType(VARCHAR) - .isEqualTo("POINT (60.0018310442288 30.121372968273892)"); + .isEqualTo("POINT (60.0018310546875 30.121372973521975)"); // Check bottom right corner of a stack of tiles at different zoom levels assertThat(assertions.function("ST_AsText", "apply(bing_tile_polygon(bing_tile(1, 1, 1)), g -> ST_Point(ST_XMax(g), ST_YMin(g)))")) diff --git a/plugin/trino-geospatial/src/test/java/io/trino/plugin/geospatial/TestGeoFunctions.java b/plugin/trino-geospatial/src/test/java/io/trino/plugin/geospatial/TestGeoFunctions.java index 96e7b5765316..bd8620850350 100644 --- a/plugin/trino-geospatial/src/test/java/io/trino/plugin/geospatial/TestGeoFunctions.java +++ b/plugin/trino-geospatial/src/test/java/io/trino/plugin/geospatial/TestGeoFunctions.java @@ -14,10 +14,12 @@ package io.trino.plugin.geospatial; import com.esri.core.geometry.Point; +import com.esri.core.geometry.ogc.OGCGeometry; import com.esri.core.geometry.ogc.OGCPoint; import com.google.common.collect.ImmutableList; import io.trino.geospatial.KdbTreeUtils; import io.trino.geospatial.Rectangle; +import io.trino.geospatial.serde.GeometrySerde; import io.trino.spi.block.Block; import io.trino.spi.block.BlockBuilder; import io.trino.spi.type.ArrayType; @@ -35,6 +37,7 @@ import static com.google.common.collect.ImmutableList.toImmutableList; import static io.trino.geospatial.KdbTree.buildKdbTree; +import static io.trino.plugin.geospatial.GeoFunctions.stCentroid; import static io.trino.plugin.geospatial.GeometryType.GEOMETRY; import static io.trino.spi.type.BooleanType.BOOLEAN; import static io.trino.spi.type.DoubleType.DOUBLE; @@ -42,6 +45,7 @@ import static io.trino.spi.type.VarcharType.VARCHAR; import static io.trino.testing.assertions.TrinoExceptionAssert.assertTrinoExceptionThrownBy; import static org.assertj.core.api.Assertions.assertThat; +import static org.junit.jupiter.api.Assertions.assertEquals; import static org.junit.jupiter.api.TestInstance.Lifecycle.PER_CLASS; import static org.junit.jupiter.api.parallel.ExecutionMode.CONCURRENT; @@ -259,6 +263,19 @@ public void testSTBuffer() assertTrinoExceptionThrownBy(assertions.function("ST_Buffer", "ST_Point(0, 0)", "nan()")::evaluate) .hasMessage("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))`. + assertThat(assertions.function("ST_IsEmpty", "ST_Buffer(ST_Buffer(ST_Point(177.50102959662, 64.726807421691), 0.0000000001), 0.00005)")) + .hasType(BOOLEAN) + .isEqualTo(true); + + assertThat(assertions.function("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)")) + .hasType(BOOLEAN) + .isEqualTo(true); } @Test @@ -299,6 +316,33 @@ public void testSTCentroid() assertThat(assertions.function("ST_AsText", "ST_Centroid(ST_GeometryFromText('POLYGON ((0 0, 0 5, 5 5, 5 0, 0 0), (1 1, 1 2, 2 2, 2 1, 1 1))'))")) .hasType(VARCHAR) .isEqualTo("POINT (2.5416666666666665 2.5416666666666665)"); + + 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 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