Implement ST_Area for SphericalGeography

presto-docs/src/main/sphinx/functions/geospatial.rst

@@ -203,6 +203,10 @@ Accessors
     For GeometryCollection types, returns the sum of the areas of the individual
     geometries.
 
+.. function:: ST_Area(SphericalGeography) -> double
+
+    Returns the area of a polygon or multi-polygon in square meters using a spherical model for Earth.
+
 .. function:: ST_Centroid(Geometry) -> Geometry
 
     Returns the point value that is the mathematical centroid of a geometry.

presto-geospatial/src/main/java/com/facebook/presto/plugin/geospatial/GeoFunctions.java

@@ -98,10 +98,12 @@
 import static com.facebook.presto.spi.type.StandardTypes.VARBINARY;
 import static com.facebook.presto.spi.type.StandardTypes.VARCHAR;
 import static com.google.common.base.Preconditions.checkArgument;
+import static com.google.common.base.Preconditions.checkState;
 import static io.airlift.slice.Slices.utf8Slice;
 import static io.airlift.slice.Slices.wrappedBuffer;
 import static java.lang.Double.isInfinite;
 import static java.lang.Double.isNaN;
+import static java.lang.Math.PI;
 import static java.lang.Math.atan2;
 import static java.lang.Math.cos;
 import static java.lang.Math.sin;
@@ -119,6 +121,7 @@ public final class GeoFunctions
     private static final Slice EMPTY_POLYGON = serialize(new OGCPolygon(new Polygon(), null));
     private static final Slice EMPTY_MULTIPOINT = serialize(createFromEsriGeometry(new MultiPoint(), null, true));
     private static final double EARTH_RADIUS_KM = 6371.01;
+    private static final double EARTH_RADIUS_M = EARTH_RADIUS_KM * 1000.0;
     private static final Map<Reason, String> NON_SIMPLE_REASONS = ImmutableMap.<Reason, String>builder()
             .put(DegenerateSegments, "Degenerate segments")
             .put(Clustering, "Repeated points")
@@ -1507,6 +1510,166 @@ private static void validateSphericalType(String function, OGCGeometry geometry,
         }
     }
 
+    @SqlNullable
+    @Description("Returns the area of a geometry on the Earth's surface using spherical model")
+    @ScalarFunction("ST_Area")
+    @SqlType(DOUBLE)
+    public static Double stSphericalArea(@SqlType(SPHERICAL_GEOGRAPHY_TYPE_NAME) Slice input)
+    {
+        OGCGeometry geometry = deserialize(input);
+        if (geometry.isEmpty()) {
+            return null;
+        }
+
+        validateSphericalType("ST_Area", geometry, EnumSet.of(POLYGON, MULTI_POLYGON));
+
+        Polygon polygon = (Polygon) geometry.getEsriGeometry();
+
+        // See https://www.movable-type.co.uk/scripts/latlong.html
+        // and http://osgeo-org.1560.x6.nabble.com/Area-of-a-spherical-polygon-td3841625.html
+        // and https://www.element84.com/blog/determining-if-a-spherical-polygon-contains-a-pole
+        // for the underlying Maths
+
+        double sphericalExcess = 0.0;
+
+        int numPaths = polygon.getPathCount();
+        for (int i = 0; i < numPaths; i++) {
+            double sign = polygon.isExteriorRing(i) ? 1.0 : -1.0;
+            sphericalExcess += sign * Math.abs(computeSphericalExcess(polygon, polygon.getPathStart(i), polygon.getPathEnd(i)));
+        }
+
+        // Math.abs is required here because for Polygons with a 2D area of 0
+        // isExteriorRing returns false for the exterior ring
+        return Math.abs(sphericalExcess * EARTH_RADIUS_M * EARTH_RADIUS_M);
+    }
+
+    private static double computeSphericalExcess(Polygon polygon, int start, int end)
+    {
+        // Our calculations rely on not processing the same point twice
+        if (polygon.getPoint(end - 1).equals(polygon.getPoint(start))) {
+            end = end - 1;
+        }
+
+        if (end - start < 3) {
+            // A path with less than 3 distinct points is not valid for calculating an area
+            throw new PrestoException(INVALID_FUNCTION_ARGUMENT, "Polygon is not valid: a loop contains less then 3 vertices.");
+        }
+
+        Point point = new Point();
+
+        // Initialize the calculator with the last point
+        polygon.getPoint(end - 1, point);
+        SphericalExcessCalculator calculator = new SphericalExcessCalculator(point);
+
+        for (int i = start; i < end; i++) {
+            polygon.getPoint(i, point);
+            calculator.add(point);
+        }
+
+        return calculator.computeSphericalExcess();
+    }
+
+    private static class SphericalExcessCalculator
+    {
+        private static final double TWO_PI = 2 * Math.PI;
+        private static final double THREE_PI = 3 * Math.PI;
+
+        private double sphericalExcess;
+        private double courseDelta;
+
+        private boolean firstPoint;
+        private double firstInitialBearing;
+        private double previousFinalBearing;
+
+        private double previousPhi;
+        private double previousCos;
+        private double previousSin;
+        private double previousTan;
+        private double previousLongitude;
+
+        private boolean done;
+
+        public SphericalExcessCalculator(Point endPoint)
+        {
+            previousPhi = toRadians(endPoint.getY());
+            previousSin = Math.sin(previousPhi);
+            previousCos = Math.cos(previousPhi);
+            previousTan = Math.tan(previousPhi / 2);
+            previousLongitude = toRadians(endPoint.getX());
+            firstPoint = true;
+        }
+
+        private void add(Point point) throws IllegalStateException
+        {
+            checkState(!done, "Computation of spherical excess is complete");
+
+            double phi = toRadians(point.getY());
+            double tan = Math.tan(phi / 2);
+            double longitude = toRadians(point.getX());
+
+            // We need to check for that specifically
+            // Otherwise calculating the bearing is not deterministic
+            if (longitude == previousLongitude && phi == previousPhi) {
+                throw new PrestoException(INVALID_FUNCTION_ARGUMENT, "Polygon is not valid: it has two identical consecutive vertices");
+            }
+
+            double deltaLongitude = longitude - previousLongitude;
+            sphericalExcess += 2 * Math.atan2(Math.tan(deltaLongitude / 2) * (previousTan + tan), 1 + previousTan * tan);
+
+            double cos = Math.cos(phi);
+            double sin = Math.sin(phi);
+            double sinOfDeltaLongitude = Math.sin(deltaLongitude);
+            double cosOfDeltaLongitude = Math.cos(deltaLongitude);
+
+            // Initial bearing from previous to current
+            double y = sinOfDeltaLongitude * cos;
+            double x = previousCos * sin - previousSin * cos * cosOfDeltaLongitude;
+            double initialBearing = (Math.atan2(y, x) + TWO_PI) % TWO_PI;
+
+            // Final bearing from previous to current = opposite of bearing from current to previous
+            double finalY = -sinOfDeltaLongitude * previousCos;
+            double finalX = previousSin * cos - previousCos * sin * cosOfDeltaLongitude;
+            double finalBearing = (Math.atan2(finalY, finalX) + PI) % TWO_PI;
+
+            // When processing our first point we don't yet have a previousFinalBearing
+            if (firstPoint) {
+                // So keep our initial bearing around, and we'll use it at the end
+                // with the last final bearing
+                firstInitialBearing = initialBearing;
+                firstPoint = false;
+            }
+            else {
+                courseDelta += (initialBearing - previousFinalBearing + THREE_PI) % TWO_PI - PI;
+            }
+
+            courseDelta += (finalBearing - initialBearing + THREE_PI) % TWO_PI - PI;
+
+            previousFinalBearing = finalBearing;
+            previousCos = cos;
+            previousSin = sin;
+            previousPhi = phi;
+            previousTan = tan;
+            previousLongitude = longitude;
+        }
+
+        public double computeSphericalExcess()
+        {
+            if (!done) {
+                // Now that we have our last final bearing, we can calculate the remaining course delta
+                courseDelta += (firstInitialBearing - previousFinalBearing + THREE_PI) % TWO_PI - PI;
+
+                // The courseDelta should be 2Pi or - 2Pi, unless a pole is enclosed (and then it should be ~ 0)
+                // In which case we need to correct the spherical excess by 2Pi
+                if (Math.abs(courseDelta) < PI / 4) {
+                    sphericalExcess = Math.abs(sphericalExcess) - TWO_PI;
+                }
+                done = true;
+            }
+
+            return sphericalExcess;
+        }
+    }
+
     private static Iterable<Slice> getGeometrySlicesFromBlock(Block block)
     {
         requireNonNull(block, "block is null");

presto-geospatial/src/test/java/com/facebook/presto/plugin/geospatial/BenchmarkSTArea.java

@@ -0,0 +1,136 @@
+/*
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package com.facebook.presto.plugin.geospatial;
+
+import io.airlift.slice.Slice;
+import org.openjdk.jmh.annotations.Benchmark;
+import org.openjdk.jmh.annotations.BenchmarkMode;
+import org.openjdk.jmh.annotations.Fork;
+import org.openjdk.jmh.annotations.Measurement;
+import org.openjdk.jmh.annotations.Mode;
+import org.openjdk.jmh.annotations.OutputTimeUnit;
+import org.openjdk.jmh.annotations.Scope;
+import org.openjdk.jmh.annotations.Setup;
+import org.openjdk.jmh.annotations.State;
+import org.openjdk.jmh.annotations.Warmup;
+import org.openjdk.jmh.runner.Runner;
+import org.openjdk.jmh.runner.RunnerException;
+import org.openjdk.jmh.runner.options.Options;
+import org.openjdk.jmh.runner.options.OptionsBuilder;
+import org.openjdk.jmh.runner.options.VerboseMode;
+import org.testng.annotations.Test;
+
+import java.io.IOException;
+import java.util.concurrent.TimeUnit;
+
+import static com.facebook.presto.plugin.geospatial.GeoFunctions.stGeometryFromText;
+import static com.facebook.presto.plugin.geospatial.GeoFunctions.toSphericalGeography;
+import static com.facebook.presto.plugin.geospatial.GeometryBenchmarkUtils.loadPolygon;
+import static io.airlift.slice.Slices.utf8Slice;
+import static org.testng.Assert.assertEquals;
+
+@State(Scope.Thread)
+@Fork(2)
+@Warmup(iterations = 10, time = 500, timeUnit = TimeUnit.MILLISECONDS)
+@Measurement(iterations = 10, time = 500, timeUnit = TimeUnit.MILLISECONDS)
+@OutputTimeUnit(TimeUnit.NANOSECONDS)
+@BenchmarkMode(Mode.AverageTime)
+public class BenchmarkSTArea
+{
+    @Benchmark
+    public Object stSphericalArea(BenchmarkData data)
+    {
+        return GeoFunctions.stSphericalArea(data.geography);
+    }
+
+    @Benchmark
+    public Object stSphericalArea500k(BenchmarkData data)
+    {
+        return GeoFunctions.stSphericalArea(data.geography500k);
+    }
+
+    @Benchmark
+    public Object stArea(BenchmarkData data)
+    {
+        return GeoFunctions.stArea(data.geometry);
+    }
+
+    @Benchmark
+    public Object stArea500k(BenchmarkData data)
+    {
+        return GeoFunctions.stArea(data.geometry500k);
+    }
+
+    @State(Scope.Thread)
+    public static class BenchmarkData
+    {
+        private Slice geometry;
+        private Slice geometry500k;
+        private Slice geography;
+        private Slice geography500k;
+
+        @Setup
+        public void setup()
+                throws IOException
+        {
+            geometry = stGeometryFromText(utf8Slice(loadPolygon("large_polygon.txt")));
+            geometry500k = stGeometryFromText(utf8Slice(createPolygon(500000)));
+            geography = toSphericalGeography(geometry);
+            geography500k = toSphericalGeography(geometry500k);
+        }
+    }
+
+    public static void main(String[] args)
+            throws IOException, RunnerException
+    {
+        // assure the benchmarks are valid before running
+        verify();
+
+        Options options = new OptionsBuilder()
+                .verbosity(VerboseMode.NORMAL)
+                .include(".*" + BenchmarkSTArea.class.getSimpleName() + ".*")
+                .build();
+        new Runner(options).run();
+    }
+
+    @Test
+    public static void verify() throws IOException
+    {
+        BenchmarkData data = new BenchmarkData();
+        data.setup();
+        BenchmarkSTArea benchmark = new BenchmarkSTArea();
+
+        assertEquals(Math.round(1000 * (Double) benchmark.stSphericalArea(data) / 3.659E8), 1000);
+        assertEquals(Math.round(1000 * (Double) benchmark.stSphericalArea500k(data) / 38842273735.0), 1000);
+        assertEquals(benchmark.stArea(data), 0.05033099592771004);
+        assertEquals(Math.round(1000 * (Double) benchmark.stArea500k(data) / Math.PI), 1000);
+    }
+
+    private static String createPolygon(double numVertices)
+    {
+        StringBuilder sb = new StringBuilder();
+        sb.append("POLYGON((");
+        String separator = "";
+        for (int i = 0; i < numVertices; i++) {
+            double angle = i * Math.PI * 2 / numVertices;
+            sb.append(separator);
+            sb.append(Math.cos(angle));
+            sb.append(" ");
+            sb.append(Math.sin(angle));
+            separator = ",";
+        }
+        sb.append("))");
+        return sb.toString();
+    }
+}