Skip to content

Commit

Permalink
Improve numerical stability for polygon centroids
Browse files Browse the repository at this point in the history 
Degenerate polygons with near-0 areas are returning centroids
with infinite coordinates.  This detects such cases, and uses
the bounding box centroid instead, which is a good approximation
for small polygons.

This centroid calculations were failing in user queries; I've used
a couple of the failing polygons as test cases.
  • Loading branch information
@jagill
jagill committed Jul 31, 2019
1 parent 3db6534 commit d7d348f
Show file tree
Hide file tree
Showing 2 changed files with 49 additions and 15 deletions.
11 changes: 11 additions & 0 deletions presto-geospatial/src/main/java/com/facebook/presto/plugin/geospatial/GeoFunctions.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,7 @@
import com.esri.core.geometry.ogc.OGCPoint;
import com.esri.core.geometry.ogc.OGCPolygon;
import com.facebook.presto.geospatial.GeometryType;
import com.facebook.presto.geospatial.GeometryUtils;
import com.facebook.presto.geospatial.KdbTree;
import com.facebook.presto.geospatial.Rectangle;
import com.facebook.presto.geospatial.serde.GeometrySerde;
Expand Down Expand Up @@ -104,6 +105,7 @@
import static java.lang.Double.isInfinite;
import static java.lang.Double.isNaN;
import static java.lang.Math.PI;
import static java.lang.Math.abs;
import static java.lang.Math.atan2;
import static java.lang.Math.cos;
import static java.lang.Math.sin;
Expand Down Expand Up @@ -1453,6 +1455,15 @@ private static Point getPolygonSansHolesCentroid(Polygon polygon)
ySum += (current.getY() + next.getY()) * ladder;
signedArea += ladder / 2;
}
if (abs(signedArea) < Double.MIN_NORMAL) {
// We can't trust the calculation, so default to bounding box.
// We've already handled the empty case, so the Envelope is non-null.
OGCGeometry ogcGeom = OGCGeometry.createFromEsriGeometry(polygon, null);
Envelope bounds = GeometryUtils.getEnvelope(ogcGeom);
return new Point(
(bounds.getXMax() + bounds.getXMin()) / 2,
(bounds.getYMax() + bounds.getYMin()) / 2);
}
return new Point(xSum / (signedArea * 6), ySum / (signedArea * 6));
}

Expand Down
53 changes: 38 additions & 15 deletions presto-geospatial/src/test/java/com/facebook/presto/plugin/geospatial/TestGeoFunctions.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,14 +14,17 @@
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;
import com.facebook.presto.spi.type.ArrayType;
import com.google.common.collect.ImmutableList;
import io.airlift.slice.Slice;
import org.testng.annotations.BeforeClass;
import org.testng.annotations.Test;

Expand All @@ -30,6 +33,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 All @@ -44,6 +48,17 @@
public class TestGeoFunctions
extends AbstractTestFunctions
{
private static String makeKdbTreeJson()
{
ImmutableList.Builder<Rectangle> rectangles = ImmutableList.builder();
for (double x = 0; x < 10; x += 1) {
for (double y = 0; y < 5; y += 1) {
rectangles.add(new Rectangle(x, y, x + 1, y + 2));
}
}
return KdbTreeUtils.toJson(buildKdbTree(10, new Rectangle(0, 0, 9, 4), rectangles.build()));
}

@BeforeClass
protected void registerFunctions()
{
Expand Down Expand Up @@ -87,17 +102,6 @@ public void testSpatialPartitions()
assertSpatialPartitions(kdbTreeJson, "MULTIPOINT (2 6, 3 7)", 1.2, ImmutableList.of());
}

private static String makeKdbTreeJson()
{
ImmutableList.Builder<Rectangle> rectangles = ImmutableList.builder();
for (double x = 0; x < 10; x += 1) {
for (double y = 0; y < 5; y += 1) {
rectangles.add(new Rectangle(x, y, x + 1, y + 2));
}
}
return KdbTreeUtils.toJson(buildKdbTree(10, new Rectangle(0, 0, 9, 4), rectangles.build()));
}

private void assertSpatialPartitions(String kdbTreeJson, String wkt, List<Integer> expectedPartitions)
{
assertFunction(format("spatial_partitions(cast('%s' as KdbTree), ST_GeometryFromText('%s'))", kdbTreeJson, wkt), new ArrayType(INTEGER), expectedPartitions);
Expand Down Expand Up @@ -212,11 +216,30 @@ public void testSTCentroid()
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.");
}

@Test
public void testSTCentroidNumericStability()
{
String wkt1 = "MULTIPOLYGON (((153.492818 -28.13729, 153.492821 -28.137291, 153.492816 -28.137289, 153.492818 -28.13729)))";
assertApproximateCentroid(wkt1, new Point(153.49282, -28.13729), 1e-5);

String wkt2 = "MULTIPOLYGON (((153.112475 -28.360526, 153.1124759 -28.360527, 153.1124759 -28.360526, 153.112475 -28.360526)))";
assertApproximateCentroid(wkt2, new Point(153.112475, -28.360526), 1e-5);
}

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 centroid, double epsilon)
{
OGCGeometry geometry = OGCGeometry.fromText(wkt);
Slice inputSlice = GeometrySerde.serialize(geometry);
OGCPoint calculatedCentroid = (OGCPoint) GeometrySerde.deserialize(stCentroid(inputSlice));
assertEquals(calculatedCentroid.X(), centroid.getX(), epsilon);
assertEquals(calculatedCentroid.Y(), centroid.getY(), epsilon);
}

@Test
public void testSTConvexHull()
{
Expand Down Expand Up @@ -1084,8 +1107,8 @@ private void assertMultiPoint(String expectedWkt, String... pointWkts)
format(
"ST_MultiPoint(array[%s])",
Arrays.stream(pointWkts)
.map(wkt -> wkt == null ? "null" : format("ST_GeometryFromText('%s')", wkt))
.collect(Collectors.joining(","))),
.map(wkt -> wkt == null ? "null" : format("ST_GeometryFromText('%s')", wkt))
.collect(Collectors.joining(","))),
GEOMETRY,
expectedWkt);
}
Expand All @@ -1096,8 +1119,8 @@ private void assertInvalidMultiPoint(String errorMessage, String... pointWkts)
format(
"ST_MultiPoint(array[%s])",
Arrays.stream(pointWkts)
.map(wkt -> wkt == null ? "null" : format("ST_GeometryFromText('%s')", wkt))
.collect(Collectors.joining(","))),
.map(wkt -> wkt == null ? "null" : format("ST_GeometryFromText('%s')", wkt))
.collect(Collectors.joining(","))),
INVALID_FUNCTION_ARGUMENT,
format("Invalid input to ST_MultiPoint: %s", errorMessage));
}
Expand Down

0 comments on commit d7d348f

Please sign in to comment.