Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix ST_Centroid and ST_Buffer for tiny geometries #13323

Merged
merged 3 commits into from
Sep 4, 2019
Merged

Fix ST_Centroid and ST_Buffer for tiny geometries #13323

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
b004edf
Upgrade ESRI geometry to v2.2.3
jagill Sep 3, 2019
1e1579f
Test for ST_Buffer NPE fix
jagill Sep 3, 2019
a0fead3
Fix centroids for small/degenerate polygons
jagill Sep 3, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion pom.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1145,7 +1145,7 @@
<dependency>
<groupId>com.esri.geometry</groupId>
<artifactId>esri-geometry-api</artifactId>
<version>2.2.2</version>
<version>2.2.3</version>
<exclusions>
<exclusion>
<groupId>com.fasterxml.jackson.core</groupId>
Expand Down
4 changes: 3 additions & 1 deletion presto-docs/src/main/sphinx/functions/geospatial.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
142 changes: 1 addition & 141 deletions presto-geospatial/src/main/java/com/facebook/presto/plugin/geospatial/GeoFunctions.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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")
Expand Down Expand Up @@ -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);
Expand Down
2 changes: 1 addition & 1 deletion ...geospatial/src/test/java/com/facebook/presto/plugin/geospatial/TestBingTileFunctions.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)");
Expand Down
40 changes: 39 additions & 1 deletion presto-geospatial/src/test/java/com/facebook/presto/plugin/geospatial/TestGeoFunctions.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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;
Expand Down Expand Up @@ -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
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Would you add a note in ST_Buffer's documentation about this behavior?

// 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
Expand All @@ -209,14 +222,39 @@ 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)
{
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()
{
Expand Down