diff --git a/presto-geospatial-toolkit/src/main/java/com/facebook/presto/geospatial/SphericalGeographyUtils.java b/presto-geospatial-toolkit/src/main/java/com/facebook/presto/geospatial/SphericalGeographyUtils.java new file mode 100644 index 0000000000000..c75e5ba3e3003 --- /dev/null +++ b/presto-geospatial-toolkit/src/main/java/com/facebook/presto/geospatial/SphericalGeographyUtils.java @@ -0,0 +1,116 @@ +/* + * 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.geospatial; + +import com.esri.core.geometry.Point; +import com.esri.core.geometry.ogc.OGCGeometry; +import com.facebook.presto.spi.PrestoException; +import com.google.common.base.Joiner; + +import java.util.EnumSet; +import java.util.Set; + +import static com.facebook.presto.geospatial.GeometryType.POINT; +import static com.facebook.presto.spi.StandardErrorCode.INVALID_FUNCTION_ARGUMENT; +import static java.lang.Math.atan2; +import static java.lang.Math.cos; +import static java.lang.Math.sin; +import static java.lang.Math.sqrt; +import static java.lang.Math.toRadians; +import static java.lang.String.format; + +public class SphericalGeographyUtils +{ + public static final double EARTH_RADIUS_KM = 6371.01; + public static final double EARTH_RADIUS_M = EARTH_RADIUS_KM * 1000.0; + private static final float MIN_LATITUDE = -90; + private static final float MAX_LATITUDE = 90; + private static final float MIN_LONGITUDE = -180; + private static final float MAX_LONGITUDE = 180; + private static final Joiner OR_JOINER = Joiner.on(" or "); + private static final Set ALLOWED_SPHERICAL_DISTANCE_TYPES = EnumSet.of(POINT); + + private SphericalGeographyUtils() {} + + public static void checkLatitude(double latitude) + { + if (Double.isNaN(latitude) || Double.isInfinite(latitude) || latitude < MIN_LATITUDE || latitude > MAX_LATITUDE) { + throw new PrestoException(INVALID_FUNCTION_ARGUMENT, "Latitude must be between -90 and 90"); + } + } + + public static void checkLongitude(double longitude) + { + if (Double.isNaN(longitude) || Double.isInfinite(longitude) || longitude < MIN_LONGITUDE || longitude > MAX_LONGITUDE) { + throw new PrestoException(INVALID_FUNCTION_ARGUMENT, "Longitude must be between -180 and 180"); + } + } + + public static Double sphericalDistance(OGCGeometry leftGeometry, OGCGeometry rightGeometry) + { + if (leftGeometry.isEmpty() || rightGeometry.isEmpty()) { + return null; + } + + validateSphericalType("ST_Distance", leftGeometry, ALLOWED_SPHERICAL_DISTANCE_TYPES); + validateSphericalType("ST_Distance", rightGeometry, ALLOWED_SPHERICAL_DISTANCE_TYPES); + Point leftPoint = (Point) leftGeometry.getEsriGeometry(); + Point rightPoint = (Point) rightGeometry.getEsriGeometry(); + + // greatCircleDistance returns distance in KM. + return greatCircleDistance(leftPoint.getY(), leftPoint.getX(), rightPoint.getY(), rightPoint.getX()) * 1000; + } + + /** + * Calculate the distance between two points on Earth. + * + * This assumes a spherical Earth, and uses the Vincenty formula. + * (https://en.wikipedia.org/wiki/Great-circle_distance) + */ + public static double greatCircleDistance( + double latitude1, + double longitude1, + double latitude2, + double longitude2) + { + checkLatitude(latitude1); + checkLongitude(longitude1); + checkLatitude(latitude2); + checkLongitude(longitude2); + + double radianLatitude1 = toRadians(latitude1); + double radianLatitude2 = toRadians(latitude2); + + double sin1 = sin(radianLatitude1); + double cos1 = cos(radianLatitude1); + double sin2 = sin(radianLatitude2); + double cos2 = cos(radianLatitude2); + + double deltaLongitude = toRadians(longitude1) - toRadians(longitude2); + double cosDeltaLongitude = cos(deltaLongitude); + + double t1 = cos2 * sin(deltaLongitude); + double t2 = cos1 * sin2 - sin1 * cos2 * cosDeltaLongitude; + double t3 = sin1 * sin2 + cos1 * cos2 * cosDeltaLongitude; + return atan2(sqrt(t1 * t1 + t2 * t2), t3) * EARTH_RADIUS_KM; + } + + public static void validateSphericalType(String function, OGCGeometry geometry, Set validTypes) + { + GeometryType type = GeometryType.getForEsriGeometryType(geometry.geometryType()); + if (!validTypes.contains(type)) { + throw new PrestoException(INVALID_FUNCTION_ARGUMENT, format("When applied to SphericalGeography inputs, %s only supports %s. Input type is: %s", function, OR_JOINER.join(validTypes), type)); + } + } +} 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 0cac1de123281..7d6a7d0a50d49 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 @@ -17,11 +17,9 @@ import com.esri.core.geometry.GeometryCursor; import com.esri.core.geometry.GeometryException; import com.esri.core.geometry.ListeningGeometryCursor; -import com.esri.core.geometry.MultiPath; import com.esri.core.geometry.NonSimpleResult.Reason; import com.esri.core.geometry.OperatorUnion; import com.esri.core.geometry.Point; -import com.esri.core.geometry.Polygon; import com.esri.core.geometry.Polyline; import com.esri.core.geometry.ogc.OGCConcreteGeometryCollection; import com.esri.core.geometry.ogc.OGCGeometry; @@ -70,7 +68,6 @@ import java.util.Objects; import java.util.Set; -import static com.esri.core.geometry.Geometry.Type; import static com.esri.core.geometry.NonSimpleResult.Reason.Clustering; import static com.esri.core.geometry.NonSimpleResult.Reason.Cracking; import static com.esri.core.geometry.NonSimpleResult.Reason.CrossOver; @@ -101,7 +98,6 @@ import static com.facebook.presto.geospatial.serde.JtsGeometrySerde.serialize; import static com.facebook.presto.plugin.geospatial.GeometryType.GEOMETRY; import static com.facebook.presto.plugin.geospatial.GeometryType.GEOMETRY_TYPE_NAME; -import static com.facebook.presto.plugin.geospatial.SphericalGeographyType.SPHERICAL_GEOGRAPHY_TYPE_NAME; import static com.facebook.presto.spi.StandardErrorCode.INVALID_FUNCTION_ARGUMENT; import static com.facebook.presto.spi.type.StandardTypes.BIGINT; import static com.facebook.presto.spi.type.StandardTypes.BOOLEAN; @@ -111,18 +107,11 @@ 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; -import static java.lang.Math.sqrt; import static java.lang.Math.toIntExact; -import static java.lang.Math.toRadians; import static java.lang.String.format; import static java.util.Arrays.setAll; import static java.util.Objects.requireNonNull; @@ -132,8 +121,6 @@ public final class GeoFunctions { private static final Joiner OR_JOINER = Joiner.on(" or "); private static final Slice EMPTY_POLYGON = serialize(createJtsEmptyPolygon()); - 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 NON_SIMPLE_REASONS = ImmutableMap.builder() .put(DegenerateSegments, "Degenerate segments") .put(Clustering, "Repeated points") @@ -146,14 +133,6 @@ public final class GeoFunctions private static final int NUMBER_OF_DIMENSIONS = 3; private static final Block EMPTY_ARRAY_OF_INTS = IntegerType.INTEGER.createFixedSizeBlockBuilder(0).build(); - private static final float MIN_LATITUDE = -90; - private static final float MAX_LATITUDE = 90; - private static final float MIN_LONGITUDE = -180; - private static final float MAX_LONGITUDE = 180; - - private static final EnumSet GEOMETRY_TYPES_FOR_SPHERICAL_GEOGRAPHY = EnumSet.of( - Type.Point, Type.Polyline, Type.Polygon, Type.MultiPoint); - private GeoFunctions() {} @Description("Returns a Geometry type LineString object from Well-Known Text representation (WKT)") @@ -272,48 +251,6 @@ public static Slice stGeomFromBinary(@SqlType(VARBINARY) Slice input) return EsriGeometrySerde.serialize(geomFromBinary(input)); } - @Description("Converts a Geometry object to a SphericalGeography object") - @ScalarFunction("to_spherical_geography") - @SqlType(SPHERICAL_GEOGRAPHY_TYPE_NAME) - public static Slice toSphericalGeography(@SqlType(GEOMETRY_TYPE_NAME) Slice input) - { - // "every point in input is in range" <=> "the envelope of input is in range" - Envelope envelope = deserializeEnvelope(input); - if (!envelope.isEmpty()) { - checkLatitude(envelope.getYMin()); - checkLatitude(envelope.getYMax()); - checkLongitude(envelope.getXMin()); - checkLongitude(envelope.getXMax()); - } - OGCGeometry geometry = EsriGeometrySerde.deserialize(input); - if (geometry.is3D()) { - throw new PrestoException(INVALID_FUNCTION_ARGUMENT, "Cannot convert 3D geometry to a spherical geography"); - } - - GeometryCursor cursor = geometry.getEsriGeometryCursor(); - while (true) { - com.esri.core.geometry.Geometry subGeometry = cursor.next(); - if (subGeometry == null) { - break; - } - - if (!GEOMETRY_TYPES_FOR_SPHERICAL_GEOGRAPHY.contains(subGeometry.getType())) { - throw new PrestoException(INVALID_FUNCTION_ARGUMENT, "Cannot convert geometry of this type to spherical geography: " + subGeometry.getType()); - } - } - - return input; - } - - @Description("Converts a SphericalGeography object to a Geometry object.") - @ScalarFunction("to_geometry") - @SqlType(GEOMETRY_TYPE_NAME) - public static Slice toGeometry(@SqlType(SPHERICAL_GEOGRAPHY_TYPE_NAME) Slice input) - { - // Every SphericalGeography object is a valid geometry object - return input; - } - @Description("Returns the Well-Known Text (WKT) representation of the geometry") @ScalarFunction("ST_AsText") @SqlType(VARCHAR) @@ -1258,7 +1195,8 @@ public static Block spatialPartitions(@SqlType(KdbTreeType.NAME) Object kdbTree, return spatialPartitions((KdbTree) kdbTree, expandedEnvelope2D); } - private static Block spatialPartitions(KdbTree kdbTree, Rectangle envelope) + // Package visible for SphericalGeoFunctions + /*package*/ static Block spatialPartitions(KdbTree kdbTree, Rectangle envelope) { Map partitions = kdbTree.findIntersectingLeaves(envelope); if (partitions.isEmpty()) { @@ -1289,51 +1227,6 @@ private static Block spatialPartitions(KdbTree kdbTree, Rectangle envelope) return blockBuilder.build(); } - @ScalarFunction - @Description("Calculates the great-circle distance between two points on the Earth's surface in kilometers") - @SqlType(DOUBLE) - public static double greatCircleDistance( - @SqlType(DOUBLE) double latitude1, - @SqlType(DOUBLE) double longitude1, - @SqlType(DOUBLE) double latitude2, - @SqlType(DOUBLE) double longitude2) - { - checkLatitude(latitude1); - checkLongitude(longitude1); - checkLatitude(latitude2); - checkLongitude(longitude2); - - double radianLatitude1 = toRadians(latitude1); - double radianLatitude2 = toRadians(latitude2); - - double sin1 = sin(radianLatitude1); - double cos1 = cos(radianLatitude1); - double sin2 = sin(radianLatitude2); - double cos2 = cos(radianLatitude2); - - double deltaLongitude = toRadians(longitude1) - toRadians(longitude2); - double cosDeltaLongitude = cos(deltaLongitude); - - double t1 = cos2 * sin(deltaLongitude); - double t2 = cos1 * sin2 - sin1 * cos2 * cosDeltaLongitude; - double t3 = sin1 * sin2 + cos1 * cos2 * cosDeltaLongitude; - return atan2(sqrt(t1 * t1 + t2 * t2), t3) * EARTH_RADIUS_KM; - } - - private static void checkLatitude(double latitude) - { - if (Double.isNaN(latitude) || Double.isInfinite(latitude) || latitude < MIN_LATITUDE || latitude > MAX_LATITUDE) { - throw new PrestoException(INVALID_FUNCTION_ARGUMENT, "Latitude must be between -90 and 90"); - } - } - - private static void checkLongitude(double longitude) - { - if (Double.isNaN(longitude) || Double.isInfinite(longitude) || longitude < MIN_LONGITUDE || longitude > MAX_LONGITUDE) { - throw new PrestoException(INVALID_FUNCTION_ARGUMENT, "Longitude must be between -180 and 180"); - } - } - private static OGCGeometry geomFromBinary(Slice input) { requireNonNull(input, "input is null"); @@ -1383,232 +1276,6 @@ private interface EnvelopesPredicate { boolean apply(Envelope left, Envelope right); } - - @SqlNullable - @Description("Returns the great-circle distance in meters between two SphericalGeography points.") - @ScalarFunction("ST_Distance") - @SqlType(DOUBLE) - public static Double stSphericalDistance(@SqlType(SPHERICAL_GEOGRAPHY_TYPE_NAME) Slice left, @SqlType(SPHERICAL_GEOGRAPHY_TYPE_NAME) Slice right) - { - OGCGeometry leftGeometry = EsriGeometrySerde.deserialize(left); - OGCGeometry rightGeometry = EsriGeometrySerde.deserialize(right); - if (leftGeometry.isEmpty() || rightGeometry.isEmpty()) { - return null; - } - - // TODO: support more SphericalGeography types. - validateSphericalType("ST_Distance", leftGeometry, EnumSet.of(POINT)); - validateSphericalType("ST_Distance", rightGeometry, EnumSet.of(POINT)); - Point leftPoint = (Point) leftGeometry.getEsriGeometry(); - Point rightPoint = (Point) rightGeometry.getEsriGeometry(); - - // greatCircleDistance returns distance in KM. - return greatCircleDistance(leftPoint.getY(), leftPoint.getX(), rightPoint.getY(), rightPoint.getX()) * 1000; - } - - private static void validateSphericalType(String function, OGCGeometry geometry, Set validTypes) - { - GeometryType type = GeometryType.getForEsriGeometryType(geometry.geometryType()); - if (!validTypes.contains(type)) { - throw new PrestoException(INVALID_FUNCTION_ARGUMENT, format("When applied to SphericalGeography inputs, %s only supports %s. Input type is: %s", function, OR_JOINER.join(validTypes), type)); - } - } - - @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 = EsriGeometrySerde.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); - } - - @SqlNullable - @Description("Returns the great-circle length in meters of a linestring or multi-linestring on Earth's surface") - @ScalarFunction("ST_Length") - @SqlType(DOUBLE) - public static Double stSphericalLength(@SqlType(SPHERICAL_GEOGRAPHY_TYPE_NAME) Slice input) - { - OGCGeometry geometry = EsriGeometrySerde.deserialize(input); - if (geometry.isEmpty()) { - return null; - } - - validateSphericalType("ST_Length", geometry, EnumSet.of(LINE_STRING, MULTI_LINE_STRING)); - MultiPath lineString = (MultiPath) geometry.getEsriGeometry(); - - double sum = 0; - - // sum up paths on (multi)linestring - for (int path = 0; path < lineString.getPathCount(); path++) { - if (lineString.getPathSize(path) < 2) { - continue; - } - - // sum up distances between adjacent points on this path - int pathStart = lineString.getPathStart(path); - Point prev = lineString.getPoint(pathStart); - for (int i = pathStart + 1; i < lineString.getPathEnd(path); i++) { - Point next = lineString.getPoint(i); - sum += greatCircleDistance(prev.getY(), prev.getX(), next.getY(), next.getX()); - prev = next; - } - } - - return sum * 1000; - } - - 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 getGeometrySlicesFromBlock(Block block) { requireNonNull(block, "block is null"); diff --git a/presto-geospatial/src/main/java/com/facebook/presto/plugin/geospatial/GeoPlugin.java b/presto-geospatial/src/main/java/com/facebook/presto/plugin/geospatial/GeoPlugin.java index b5305edf4a4ba..a44efc500f736 100644 --- a/presto-geospatial/src/main/java/com/facebook/presto/plugin/geospatial/GeoPlugin.java +++ b/presto-geospatial/src/main/java/com/facebook/presto/plugin/geospatial/GeoPlugin.java @@ -50,6 +50,7 @@ public Set> getFunctions() .add(KdbTreeCasts.class) .add(SpatialPartitioningAggregateFunction.class) .add(SpatialPartitioningInternalAggregateFunction.class) + .add(SphericalGeoFunctions.class) .build(); } } diff --git a/presto-geospatial/src/main/java/com/facebook/presto/plugin/geospatial/SphericalGeoFunctions.java b/presto-geospatial/src/main/java/com/facebook/presto/plugin/geospatial/SphericalGeoFunctions.java new file mode 100644 index 0000000000000..5735ef1c79c8c --- /dev/null +++ b/presto-geospatial/src/main/java/com/facebook/presto/plugin/geospatial/SphericalGeoFunctions.java @@ -0,0 +1,366 @@ +/* + * 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 com.esri.core.geometry.Envelope; +import com.esri.core.geometry.Geometry; +import com.esri.core.geometry.Geometry.Type; +import com.esri.core.geometry.GeometryCursor; +import com.esri.core.geometry.MultiPath; +import com.esri.core.geometry.Point; +import com.esri.core.geometry.Polygon; +import com.esri.core.geometry.ogc.OGCGeometry; +import com.facebook.presto.geospatial.KdbTree; +import com.facebook.presto.geospatial.Rectangle; +import com.facebook.presto.geospatial.SphericalGeographyUtils; +import com.facebook.presto.geospatial.serde.EsriGeometrySerde; +import com.facebook.presto.spi.PrestoException; +import com.facebook.presto.spi.block.Block; +import com.facebook.presto.spi.function.Description; +import com.facebook.presto.spi.function.ScalarFunction; +import com.facebook.presto.spi.function.SqlNullable; +import com.facebook.presto.spi.function.SqlType; +import com.facebook.presto.spi.type.KdbTreeType; +import io.airlift.slice.Slice; + +import java.util.EnumSet; + +import static com.facebook.presto.geospatial.GeometryType.LINE_STRING; +import static com.facebook.presto.geospatial.GeometryType.MULTI_LINE_STRING; +import static com.facebook.presto.geospatial.GeometryType.MULTI_POLYGON; +import static com.facebook.presto.geospatial.GeometryType.POLYGON; +import static com.facebook.presto.geospatial.SphericalGeographyUtils.EARTH_RADIUS_M; +import static com.facebook.presto.geospatial.SphericalGeographyUtils.checkLatitude; +import static com.facebook.presto.geospatial.SphericalGeographyUtils.checkLongitude; +import static com.facebook.presto.geospatial.SphericalGeographyUtils.sphericalDistance; +import static com.facebook.presto.geospatial.SphericalGeographyUtils.validateSphericalType; +import static com.facebook.presto.geospatial.serde.EsriGeometrySerde.deserializeEnvelope; +import static com.facebook.presto.plugin.geospatial.GeometryType.GEOMETRY_TYPE_NAME; +import static com.facebook.presto.plugin.geospatial.SphericalGeographyType.SPHERICAL_GEOGRAPHY_TYPE_NAME; +import static com.facebook.presto.spi.StandardErrorCode.INVALID_FUNCTION_ARGUMENT; +import static com.facebook.presto.spi.type.StandardTypes.DOUBLE; +import static com.google.common.base.Preconditions.checkState; +import static java.lang.Double.isInfinite; +import static java.lang.Double.isNaN; +import static java.lang.Math.PI; +import static java.lang.Math.toRadians; + +public final class SphericalGeoFunctions +{ + private static final EnumSet GEOMETRY_TYPES_FOR_SPHERICAL_GEOGRAPHY = EnumSet.of( + Type.Point, Type.Polyline, Type.Polygon, Type.MultiPoint); + + private SphericalGeoFunctions() {} + + @Description("Converts a Geometry object to a SphericalGeography object") + @ScalarFunction("to_spherical_geography") + @SqlType(SPHERICAL_GEOGRAPHY_TYPE_NAME) + public static Slice toSphericalGeography(@SqlType(GEOMETRY_TYPE_NAME) Slice input) + { + // "every point in input is in range" <=> "the envelope of input is in range" + Envelope envelope = deserializeEnvelope(input); + if (!envelope.isEmpty()) { + checkLatitude(envelope.getYMin()); + checkLatitude(envelope.getYMax()); + checkLongitude(envelope.getXMin()); + checkLongitude(envelope.getXMax()); + } + OGCGeometry geometry = EsriGeometrySerde.deserialize(input); + if (geometry.is3D()) { + throw new PrestoException(INVALID_FUNCTION_ARGUMENT, "Cannot convert 3D geometry to a spherical geography"); + } + + GeometryCursor cursor = geometry.getEsriGeometryCursor(); + while (true) { + com.esri.core.geometry.Geometry subGeometry = cursor.next(); + if (subGeometry == null) { + break; + } + + if (!GEOMETRY_TYPES_FOR_SPHERICAL_GEOGRAPHY.contains(subGeometry.getType())) { + throw new PrestoException(INVALID_FUNCTION_ARGUMENT, "Cannot convert geometry of this type to spherical geography: " + subGeometry.getType()); + } + } + + return input; + } + + @Description("Converts a SphericalGeography object to a Geometry object.") + @ScalarFunction("to_geometry") + @SqlType(GEOMETRY_TYPE_NAME) + public static Slice toGeometry(@SqlType(SPHERICAL_GEOGRAPHY_TYPE_NAME) Slice input) + { + // Every SphericalGeography object is a valid geometry object + return input; + } + + @SqlNullable + @Description("Returns the great-circle distance in meters between two SphericalGeography points.") + @ScalarFunction("ST_Distance") + @SqlType(DOUBLE) + public static Double stSphericalDistance(@SqlType(SPHERICAL_GEOGRAPHY_TYPE_NAME) Slice left, @SqlType(SPHERICAL_GEOGRAPHY_TYPE_NAME) Slice right) + { + return sphericalDistance(EsriGeometrySerde.deserialize(left), EsriGeometrySerde.deserialize(right)); + } + + @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 = EsriGeometrySerde.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); + } + + @ScalarFunction + @Description("Calculates the great-circle distance between two points on the Earth's surface in kilometers") + @SqlType(DOUBLE) + public static double greatCircleDistance( + @SqlType(DOUBLE) double latitude1, + @SqlType(DOUBLE) double longitude1, + @SqlType(DOUBLE) double latitude2, + @SqlType(DOUBLE) double longitude2) + { + return SphericalGeographyUtils.greatCircleDistance(latitude1, longitude1, latitude2, longitude2); + } + + @ScalarFunction + @SqlNullable + @Description("Returns an array of spatial partition IDs for a given geometry") + @SqlType("array(int)") + public static Block spatialPartitions(@SqlType(KdbTreeType.NAME) Object kdbTree, @SqlType(SPHERICAL_GEOGRAPHY_TYPE_NAME) Slice geometry) + { + Envelope envelope = deserializeEnvelope(geometry); + if (envelope.isEmpty()) { + // Empty geometry + return null; + } + + return GeoFunctions.spatialPartitions((KdbTree) kdbTree, new Rectangle(envelope.getXMin(), envelope.getYMin(), envelope.getXMax(), envelope.getYMax())); + } + + @ScalarFunction + @SqlNullable + @Description("Returns an array of spatial partition IDs for a geometry representing a set of points within specified distance from the input geometry") + @SqlType("array(int)") + public static Block spatialPartitions(@SqlType(KdbTreeType.NAME) Object kdbTree, @SqlType(SPHERICAL_GEOGRAPHY_TYPE_NAME) Slice geometry, @SqlType(DOUBLE) double distance) + { + if (isNaN(distance)) { + throw new PrestoException(INVALID_FUNCTION_ARGUMENT, "distance is NaN"); + } + + if (isInfinite(distance)) { + throw new PrestoException(INVALID_FUNCTION_ARGUMENT, "distance is infinite"); + } + + if (distance < 0) { + throw new PrestoException(INVALID_FUNCTION_ARGUMENT, "distance is negative"); + } + + Envelope envelope = deserializeEnvelope(geometry); + if (envelope.isEmpty()) { + return null; + } + + Rectangle expandedEnvelope2D = new Rectangle(envelope.getXMin() - distance, envelope.getYMin() - distance, envelope.getXMax() + distance, envelope.getYMax() + distance); + return GeoFunctions.spatialPartitions((KdbTree) kdbTree, expandedEnvelope2D); + } + + @SqlNullable + @Description("Returns the great-circle length in meters of a linestring or multi-linestring on Earth's surface") + @ScalarFunction("ST_Length") + @SqlType(DOUBLE) + public static Double stSphericalLength(@SqlType(SPHERICAL_GEOGRAPHY_TYPE_NAME) Slice input) + { + OGCGeometry geometry = EsriGeometrySerde.deserialize(input); + if (geometry.isEmpty()) { + return null; + } + + validateSphericalType("ST_Length", geometry, EnumSet.of(LINE_STRING, MULTI_LINE_STRING)); + MultiPath lineString = (MultiPath) geometry.getEsriGeometry(); + + double sum = 0; + + // sum up paths on (multi)linestring + for (int path = 0; path < lineString.getPathCount(); path++) { + if (lineString.getPathSize(path) < 2) { + continue; + } + + // sum up distances between adjacent points on this path + int pathStart = lineString.getPathStart(path); + Point prev = lineString.getPoint(pathStart); + for (int i = pathStart + 1; i < lineString.getPathEnd(path); i++) { + Point next = lineString.getPoint(i); + sum += greatCircleDistance(prev.getY(), prev.getX(), next.getY(), next.getX()); + prev = next; + } + } + + return sum * 1000; + } + + 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; + } + } +} diff --git a/presto-geospatial/src/test/java/com/facebook/presto/plugin/geospatial/BenchmarkSTArea.java b/presto-geospatial/src/test/java/com/facebook/presto/plugin/geospatial/BenchmarkSTArea.java index fc8889052ce72..a6eacdee435e8 100644 --- a/presto-geospatial/src/test/java/com/facebook/presto/plugin/geospatial/BenchmarkSTArea.java +++ b/presto-geospatial/src/test/java/com/facebook/presto/plugin/geospatial/BenchmarkSTArea.java @@ -35,9 +35,9 @@ 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.createCirclePolygon; import static com.facebook.presto.plugin.geospatial.GeometryBenchmarkUtils.loadPolygon; +import static com.facebook.presto.plugin.geospatial.SphericalGeoFunctions.toSphericalGeography; import static io.airlift.slice.Slices.utf8Slice; import static org.testng.Assert.assertEquals; @@ -52,13 +52,13 @@ public class BenchmarkSTArea @Benchmark public Object stSphericalArea(BenchmarkData data) { - return GeoFunctions.stSphericalArea(data.geography); + return SphericalGeoFunctions.stSphericalArea(data.geography); } @Benchmark public Object stSphericalArea500k(BenchmarkData data) { - return GeoFunctions.stSphericalArea(data.geography500k); + return SphericalGeoFunctions.stSphericalArea(data.geography500k); } @Benchmark diff --git a/presto-geospatial/src/test/java/com/facebook/presto/plugin/geospatial/TestExtractSpatialInnerJoin.java b/presto-geospatial/src/test/java/com/facebook/presto/plugin/geospatial/TestExtractSpatialInnerJoin.java index d2383f2ee15d9..539d5a9779885 100644 --- a/presto-geospatial/src/test/java/com/facebook/presto/plugin/geospatial/TestExtractSpatialInnerJoin.java +++ b/presto-geospatial/src/test/java/com/facebook/presto/plugin/geospatial/TestExtractSpatialInnerJoin.java @@ -106,30 +106,7 @@ public void testDoesNotFire() p.values(p.variable("b", GEOMETRY))))) .doesNotFire(); - // SphericalGeography operand - assertRuleApplication() - .on(p -> - p.filter( - sqlToRowExpression( - "ST_Distance(a, b) < 5", - ImmutableMap.of("a", SPHERICAL_GEOGRAPHY, "b", SPHERICAL_GEOGRAPHY)), - p.join(INNER, - p.values(p.variable("a", SPHERICAL_GEOGRAPHY)), - p.values(p.variable("b", SPHERICAL_GEOGRAPHY))))) - .doesNotFire(); - - // to_spherical_geography() operand - assertRuleApplication() - .on(p -> - p.filter( - sqlToRowExpression( - "ST_Distance(to_spherical_geography(ST_GeometryFromText(wkt)), point) < 5", - ImmutableMap.of("wkt", VARCHAR, "point", SPHERICAL_GEOGRAPHY)), - p.join(INNER, - p.values(p.variable("wkt", VARCHAR)), - p.values(p.variable("point", SPHERICAL_GEOGRAPHY))))) - .doesNotFire(); - + // We do not support spherical ST_Contains yet assertRuleApplication() .on(p -> p.filter(PlanBuilder.expression("ST_Contains(to_spherical_geography(ST_GeometryFromText(wkt)), point)"), @@ -287,6 +264,46 @@ private void testPointAndRadiusExpressionsInDistanceQuery(String filter, String values(ImmutableMap.of("lat_b", 0, "lng_b", 1, "name_b", 2)))))); } + @Test + public void testSphericalGeographiesDoesFire() + { + assertRuleApplication() + .on(p -> + p.filter( + sqlToRowExpression( + "ST_Distance(a, b) < 5000", + ImmutableMap.of("a", SPHERICAL_GEOGRAPHY, "b", SPHERICAL_GEOGRAPHY)), + p.join(INNER, + p.values(p.variable("a", SPHERICAL_GEOGRAPHY)), + p.values(p.variable("b", SPHERICAL_GEOGRAPHY))))) + .matches( + spatialJoin("ST_Distance(a, b) < radius", + values(ImmutableMap.of("a", 0)), + project(ImmutableMap.of( + "b", expression("b"), + "radius", expression("BIGINT '5000'")), + values(ImmutableMap.of("b", 0))))); + + // to_spherical_geography() operand + assertRuleApplication() + .on(p -> + p.filter( + sqlToRowExpression( + "ST_Distance(to_spherical_geography(ST_GeometryFromText(wkt)), point) < 5", + ImmutableMap.of("wkt", VARCHAR, "point", SPHERICAL_GEOGRAPHY)), + p.join(INNER, + p.values(p.variable("wkt", VARCHAR)), + p.values(p.variable("point", SPHERICAL_GEOGRAPHY))))) + .matches( + spatialJoin("ST_Distance(to_spherical_geography, point) < radius", + project(ImmutableMap.of("to_spherical_geography", expression("to_spherical_geography(ST_GeometryFromText(wkt))")), + values(ImmutableMap.of("wkt", 0))), + project(ImmutableMap.of( + "point", expression("point"), + "radius", expression("INTEGER '5'")), + values(ImmutableMap.of("point", 0))))); + } + @Test public void testConvertToSpatialJoin() { 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 e7068843fff20..7ae0fb00e74a2 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 @@ -833,29 +833,6 @@ public void testSTSymmetricDifference() assertFunction("ST_AsText(ST_SymDifference(ST_GeometryFromText('MULTIPOLYGON (((0 0, 0 2, 2 2, 2 0, 0 0)), ((2 2, 2 4, 4 4, 4 2, 2 2)))'), ST_GeometryFromText('POLYGON ((0 0, 0 3, 3 3, 3 0, 0 0))')))", VARCHAR, "MULTIPOLYGON (((2 0, 2 2, 3 2, 3 0, 2 0)), ((0 2, 0 3, 2 3, 2 2, 0 2)), ((3 2, 3 3, 2 3, 2 4, 4 4, 4 2, 3 2)))"); } - @Test - public void testGreatCircleDistance() - { - assertFunction("great_circle_distance(36.12, -86.67, 33.94, -118.40)", DOUBLE, 2886.448973436703); - assertFunction("great_circle_distance(33.94, -118.40, 36.12, -86.67)", DOUBLE, 2886.448973436703); - assertFunction("great_circle_distance(42.3601, -71.0589, 42.4430, -71.2290)", DOUBLE, 16.73469743457461); - assertFunction("great_circle_distance(36.12, -86.67, 36.12, -86.67)", DOUBLE, 0.0); - - assertInvalidFunction("great_circle_distance(100, 20, 30, 40)", "Latitude must be between -90 and 90"); - assertInvalidFunction("great_circle_distance(10, 20, 300, 40)", "Latitude must be between -90 and 90"); - assertInvalidFunction("great_circle_distance(10, 200, 30, 40)", "Longitude must be between -180 and 180"); - assertInvalidFunction("great_circle_distance(10, 20, 30, 400)", "Longitude must be between -180 and 180"); - - assertInvalidFunction("great_circle_distance(nan(), -86.67, 33.94, -118.40)", "Latitude must be between -90 and 90"); - assertInvalidFunction("great_circle_distance(infinity(), -86.67, 33.94, -118.40)", "Latitude must be between -90 and 90"); - assertInvalidFunction("great_circle_distance(36.12, nan(), 33.94, -118.40)", "Longitude must be between -180 and 180"); - assertInvalidFunction("great_circle_distance(36.12, infinity(), 33.94, -118.40)", "Longitude must be between -180 and 180"); - assertInvalidFunction("great_circle_distance(36.12, -86.67, nan(), -118.40)", "Latitude must be between -90 and 90"); - assertInvalidFunction("great_circle_distance(36.12, -86.67, infinity(), -118.40)", "Latitude must be between -90 and 90"); - assertInvalidFunction("great_circle_distance(36.12, -86.67, 33.94, nan())", "Longitude must be between -180 and 180"); - assertInvalidFunction("great_circle_distance(36.12, -86.67, 33.94, infinity())", "Longitude must be between -180 and 180"); - } - @Test public void testSTInteriorRings() { diff --git a/presto-geospatial/src/test/java/com/facebook/presto/plugin/geospatial/TestSpatialJoinPlanning.java b/presto-geospatial/src/test/java/com/facebook/presto/plugin/geospatial/TestSpatialJoinPlanning.java index a6cc14c5d7b03..125eb55c887e1 100644 --- a/presto-geospatial/src/test/java/com/facebook/presto/plugin/geospatial/TestSpatialJoinPlanning.java +++ b/presto-geospatial/src/test/java/com/facebook/presto/plugin/geospatial/TestSpatialJoinPlanning.java @@ -311,6 +311,56 @@ public void testDistanceQuery() anyTree(values(ImmutableMap.of("b_lng", 0, "b_lat", 1)))))))))); } + @Test + public void testSphericalDistanceQuery() + { + String queryFormat = "SELECT b.name, a.name " + + "FROM (VALUES (2.1, 2.1, 'x')) AS a (lng, lat, name), (VALUES (2.1, 2.1, 'x')) AS b (lng, lat, name) " + + "WHERE ST_Distance(to_spherical_geography(ST_Point(a.lng, a.lat)), to_spherical_geography(ST_Point(b.lng, b.lat))) %s 3.1"; + String queryLessThan = format(queryFormat, "<"); + String queryLessThanEquals = format(queryFormat, "<="); + + // broadcast + assertPlan(queryLessThan, anyTree( + spatialJoin("st_distance(st_point_a, st_point_b) < radius", + project(ImmutableMap.of("st_point_a", expression("to_spherical_geography(ST_Point(cast(a_lng as double), cast(a_lat as double)))")), + anyTree(values(ImmutableMap.of("a_lng", 0, "a_lat", 1)))), + anyTree(project(ImmutableMap.of("st_point_b", expression("to_spherical_geography(ST_Point(cast(b_lng as double), cast(b_lat as double)))"), "radius", expression("3.1e0")), + anyTree(values(ImmutableMap.of("b_lng", 0, "b_lat", 1)))))))); + + assertPlan(queryLessThanEquals, anyTree( + spatialJoin("st_distance(st_point_a, st_point_b) <= radius", + project(ImmutableMap.of("st_point_a", expression("to_spherical_geography(ST_Point(cast(a_lng as double), cast(a_lat as double)))")), + anyTree(values(ImmutableMap.of("a_lng", 0, "a_lat", 1)))), + anyTree(project(ImmutableMap.of("st_point_b", expression("to_spherical_geography(ST_Point(cast(b_lng as double), cast(b_lat as double)))"), "radius", expression("3.1e0")), + anyTree(values(ImmutableMap.of("b_lng", 0, "b_lat", 1)))))))); + + // distributed + assertDistributedPlan(queryLessThan, withSpatialPartitioning("memory.default.kdb_tree"), + anyTree(spatialJoin("st_distance(st_point_a, st_point_b) < radius", Optional.of(KDB_TREE_JSON), + anyTree(unnest( + project(ImmutableMap.of("partitions", expression(format("spatial_partitions(cast('%s' as kdbtree), st_point_a)", KDB_TREE_JSON))), + project(ImmutableMap.of("st_point_a", expression("to_spherical_geography(ST_Point(cast(a_lng as double), cast(a_lat as double)))")), + anyTree(values(ImmutableMap.of("a_lng", 0, "a_lat", 1))))))), + anyTree(unnest( + project(ImmutableMap.of("partitions", expression(format("spatial_partitions(cast('%s' as kdbtree), st_point_b, 3.1e0)", KDB_TREE_JSON)), + "radius", expression("3.1e0")), + project(ImmutableMap.of("st_point_b", expression("to_spherical_geography(ST_Point(cast(b_lng as double), cast(b_lat as double)))")), + anyTree(values(ImmutableMap.of("b_lng", 0, "b_lat", 1)))))))))); + + assertDistributedPlan(queryLessThanEquals, withSpatialPartitioning("memory.default.kdb_tree"), anyTree( + spatialJoin("st_distance(st_point_a, st_point_b) <= radius", Optional.of(KDB_TREE_JSON), + anyTree(unnest( + project(ImmutableMap.of("partitions", expression(format("spatial_partitions(cast('%s' as kdbtree), st_point_a)", KDB_TREE_JSON))), + project(ImmutableMap.of("st_point_a", expression("to_spherical_geography(ST_Point(cast(a_lng as double), cast(a_lat as double)))")), + anyTree(values(ImmutableMap.of("a_lng", 0, "a_lat", 1))))))), + anyTree( + project(ImmutableMap.of("partitions", expression(format("spatial_partitions(cast('%s' as kdbtree), st_point_b, 3.1e0)", KDB_TREE_JSON)), + "radius", expression("3.1e0")), + project(ImmutableMap.of("st_point_b", expression("to_spherical_geography(ST_Point(cast(b_lng as double), cast(b_lat as double)))")), + anyTree(values(ImmutableMap.of("b_lng", 0, "b_lat", 1))))))))); + } + @Test public void testNotContains() { diff --git a/presto-geospatial/src/test/java/com/facebook/presto/plugin/geospatial/TestSpatialJoins.java b/presto-geospatial/src/test/java/com/facebook/presto/plugin/geospatial/TestSpatialJoins.java index 7c40638014ecf..7675268f1e583 100644 --- a/presto-geospatial/src/test/java/com/facebook/presto/plugin/geospatial/TestSpatialJoins.java +++ b/presto-geospatial/src/test/java/com/facebook/presto/plugin/geospatial/TestSpatialJoins.java @@ -385,4 +385,35 @@ public void testRelationshipBroadcastSpatialJoin() testRelationshipSpatialJoin(getSession(), "ST_Overlaps", OVERLAPS_PAIRS); testRelationshipSpatialJoin(getSession(), "ST_Crosses", CROSSES_PAIRS); } + + @Test + public void testSphericalSpatialJoin() + { + double tooSmallRadius = 10; + double sufficientRadius = 111200; + // Ensure that correct spherical geography filter predicate is used for + // the join. If planar geometry is used, the distance will be 1, but if + // spherical geography is used, it will be ~111195. + assertSphericalSpatialJoin("<", sufficientRadius, true); + assertSphericalSpatialJoin("<", tooSmallRadius, false); + assertSphericalSpatialJoin("<=", sufficientRadius, true); + assertSphericalSpatialJoin("<=", tooSmallRadius, false); + } + + private void assertSphericalSpatialJoin(String comparison, double radius, boolean shouldMatch) + { + String expected; + if (shouldMatch) { + expected = "SELECT 'p0', 'p1'"; + } + else { + expected = "SELECT * FROM (VALUES 1) LIMIT 0"; + } + + assertQuery(format("SELECT a.name, b.name FROM " + + "( VALUES (to_spherical_geography(ST_Point(0, 0)), 'p0') ) as a(point, name) " + + "JOIN ( VALUES (to_spherical_geography(ST_Point(0, 1)), 'p1') ) as b(point, name) " + + "ON ST_Distance(a.point, b.point) %s %s", comparison, radius), + expected); + } } diff --git a/presto-geospatial/src/test/java/com/facebook/presto/plugin/geospatial/TestSphericalGeoFunctions.java b/presto-geospatial/src/test/java/com/facebook/presto/plugin/geospatial/TestSphericalGeoFunctions.java index 491057508a207..49ad1d4a2a54f 100644 --- a/presto-geospatial/src/test/java/com/facebook/presto/plugin/geospatial/TestSphericalGeoFunctions.java +++ b/presto-geospatial/src/test/java/com/facebook/presto/plugin/geospatial/TestSphericalGeoFunctions.java @@ -30,6 +30,7 @@ import java.util.stream.Collectors; import static com.facebook.presto.metadata.FunctionExtractor.extractFunctions; +import static com.facebook.presto.plugin.geospatial.SphericalGeoFunctions.toSphericalGeography; import static com.facebook.presto.plugin.geospatial.SphericalGeographyType.SPHERICAL_GEOGRAPHY; import static com.facebook.presto.spi.type.DoubleType.DOUBLE; import static com.facebook.presto.spi.type.VarcharType.VARCHAR; @@ -72,7 +73,7 @@ public void testGetObjectValue() BlockBuilder builder = SPHERICAL_GEOGRAPHY.createBlockBuilder(null, wktList.size()); for (String wkt : wktList) { - SPHERICAL_GEOGRAPHY.writeSlice(builder, GeoFunctions.toSphericalGeography(GeoFunctions.stGeometryFromText(utf8Slice(wkt)))); + SPHERICAL_GEOGRAPHY.writeSlice(builder, toSphericalGeography(GeoFunctions.stGeometryFromText(utf8Slice(wkt)))); } Block block = builder.build(); for (int i = 0; i < wktList.size(); i++) { @@ -133,6 +134,29 @@ private void assertInvalidLatitude(String wkt) assertInvalidFunction(format("to_spherical_geography(ST_GeometryFromText('%s'))", wkt), "Latitude must be between -90 and 90"); } + @Test + public void testGreatCircleDistance() + { + assertFunction("great_circle_distance(36.12, -86.67, 33.94, -118.40)", DOUBLE, 2886.448973436703); + assertFunction("great_circle_distance(33.94, -118.40, 36.12, -86.67)", DOUBLE, 2886.448973436703); + assertFunction("great_circle_distance(42.3601, -71.0589, 42.4430, -71.2290)", DOUBLE, 16.73469743457461); + assertFunction("great_circle_distance(36.12, -86.67, 36.12, -86.67)", DOUBLE, 0.0); + + assertInvalidFunction("great_circle_distance(100, 20, 30, 40)", "Latitude must be between -90 and 90"); + assertInvalidFunction("great_circle_distance(10, 20, 300, 40)", "Latitude must be between -90 and 90"); + assertInvalidFunction("great_circle_distance(10, 200, 30, 40)", "Longitude must be between -180 and 180"); + assertInvalidFunction("great_circle_distance(10, 20, 30, 400)", "Longitude must be between -180 and 180"); + + assertInvalidFunction("great_circle_distance(nan(), -86.67, 33.94, -118.40)", "Latitude must be between -90 and 90"); + assertInvalidFunction("great_circle_distance(infinity(), -86.67, 33.94, -118.40)", "Latitude must be between -90 and 90"); + assertInvalidFunction("great_circle_distance(36.12, nan(), 33.94, -118.40)", "Longitude must be between -180 and 180"); + assertInvalidFunction("great_circle_distance(36.12, infinity(), 33.94, -118.40)", "Longitude must be between -180 and 180"); + assertInvalidFunction("great_circle_distance(36.12, -86.67, nan(), -118.40)", "Latitude must be between -90 and 90"); + assertInvalidFunction("great_circle_distance(36.12, -86.67, infinity(), -118.40)", "Latitude must be between -90 and 90"); + assertInvalidFunction("great_circle_distance(36.12, -86.67, 33.94, nan())", "Longitude must be between -180 and 180"); + assertInvalidFunction("great_circle_distance(36.12, -86.67, 33.94, infinity())", "Longitude must be between -180 and 180"); + } + @Test public void testDistance() { diff --git a/presto-main/src/main/java/com/facebook/presto/sql/planner/LocalExecutionPlanner.java b/presto-main/src/main/java/com/facebook/presto/sql/planner/LocalExecutionPlanner.java index 1e951007a07ae..cd4476acf0c61 100644 --- a/presto-main/src/main/java/com/facebook/presto/sql/planner/LocalExecutionPlanner.java +++ b/presto-main/src/main/java/com/facebook/presto/sql/planner/LocalExecutionPlanner.java @@ -141,6 +141,7 @@ import com.facebook.presto.spi.relation.RowExpression; import com.facebook.presto.spi.relation.VariableReferenceExpression; import com.facebook.presto.spi.type.Type; +import com.facebook.presto.spi.type.TypeSignature; import com.facebook.presto.spiller.PartitioningSpillerFactory; import com.facebook.presto.spiller.SingleStreamSpillerFactory; import com.facebook.presto.spiller.SpillerFactory; @@ -216,6 +217,7 @@ import java.util.Set; import java.util.concurrent.atomic.AtomicInteger; import java.util.function.Function; +import java.util.function.Predicate; import java.util.function.Supplier; import java.util.stream.Collectors; import java.util.stream.IntStream; @@ -232,6 +234,7 @@ import static com.facebook.presto.SystemSessionProperties.isSpillEnabled; import static com.facebook.presto.expressions.LogicalRowExpressions.TRUE_CONSTANT; import static com.facebook.presto.expressions.RowExpressionNodeInliner.replaceExpression; +import static com.facebook.presto.geospatial.SphericalGeographyUtils.sphericalDistance; import static com.facebook.presto.operator.DistinctLimitOperator.DistinctLimitOperatorFactory; import static com.facebook.presto.operator.NestedLoopBuildOperator.NestedLoopBuildOperatorFactory; import static com.facebook.presto.operator.NestedLoopJoinOperator.NestedLoopJoinOperatorFactory; @@ -247,11 +250,13 @@ import static com.facebook.presto.operator.WindowFunctionDefinition.window; import static com.facebook.presto.operator.unnest.UnnestOperator.UnnestOperatorFactory; import static com.facebook.presto.spi.StandardErrorCode.COMPILER_ERROR; +import static com.facebook.presto.spi.StandardErrorCode.GENERIC_INTERNAL_ERROR; import static com.facebook.presto.spi.plan.AggregationNode.Step.FINAL; import static com.facebook.presto.spi.plan.AggregationNode.Step.INTERMEDIATE; import static com.facebook.presto.spi.plan.AggregationNode.Step.PARTIAL; import static com.facebook.presto.spi.relation.ExpressionOptimizer.Level.OPTIMIZED; import static com.facebook.presto.spi.type.BigintType.BIGINT; +import static com.facebook.presto.spi.type.TypeSignature.parseTypeSignature; import static com.facebook.presto.spi.type.TypeUtils.writeNativeValue; import static com.facebook.presto.sql.gen.LambdaBytecodeGenerator.compileLambdaProvider; import static com.facebook.presto.sql.planner.RowExpressionInterpreter.rowExpressionInterpreter; @@ -316,6 +321,8 @@ public class LocalExecutionPlanner private final OrderingCompiler orderingCompiler; private final JsonCodec tableCommitContextCodec; + private static final TypeSignature SPHERICAL_GEOGRAPHY_TYPE_SIGNATURE = parseTypeSignature("SphericalGeography"); + @Inject public LocalExecutionPlanner( Metadata metadata, @@ -1744,7 +1751,24 @@ private Optional removeExpressionFromFilter(RowExpression filter, private SpatialPredicate spatialTest(CallExpression functionCall, boolean probeFirst, Optional comparisonOperator) { - QualifiedFunctionName functionName = metadata.getFunctionManager().getFunctionMetadata(functionCall.getFunctionHandle()).getName(); + FunctionMetadata functionMetadata = metadata.getFunctionManager().getFunctionMetadata(functionCall.getFunctionHandle()); + QualifiedFunctionName functionName = functionMetadata.getName(); + List argumentTypes = functionMetadata.getArgumentTypes(); + Predicate isSpherical = (typeSignature) + -> typeSignature.equals(SPHERICAL_GEOGRAPHY_TYPE_SIGNATURE); + if (argumentTypes.stream().allMatch(isSpherical)) { + return sphericalSpatialTest(functionName, comparisonOperator); + } + else if (argumentTypes.stream().noneMatch(isSpherical)) { + return euclideanSpatialTest(functionName, comparisonOperator, probeFirst); + } + else { + throw new PrestoException(GENERIC_INTERNAL_ERROR, "Mixing spherical and euclidean geometric types"); + } + } + + private SpatialPredicate euclideanSpatialTest(QualifiedFunctionName functionName, Optional comparisonOperator, boolean probeFirst) + { if (functionName.equals(ST_CONTAINS)) { if (probeFirst) { return (buildGeometry, probeGeometry, radius) -> probeGeometry.contains(buildGeometry); @@ -1784,12 +1808,28 @@ else if (comparisonOperator.get() == OperatorType.LESS_THAN_OR_EQUAL) { return (buildGeometry, probeGeometry, radius) -> buildGeometry.distance(probeGeometry) <= radius.getAsDouble(); } else { - throw new UnsupportedOperationException("Unsupported comparison operator: " + comparisonOperator.get()); + throw new UnsupportedOperationException("Unsupported comparison operator: " + comparisonOperator); } } throw new UnsupportedOperationException("Unsupported spatial function: " + functionName); } + private SpatialPredicate sphericalSpatialTest(QualifiedFunctionName functionName, Optional comparisonOperator) + { + if (functionName.equals(ST_DISTANCE)) { + if (comparisonOperator.get() == OperatorType.LESS_THAN) { + return (buildGeometry, probeGeometry, radius) -> sphericalDistance(buildGeometry, probeGeometry) < radius.getAsDouble(); + } + else if (comparisonOperator.get() == OperatorType.LESS_THAN_OR_EQUAL) { + return (buildGeometry, probeGeometry, radius) -> sphericalDistance(buildGeometry, probeGeometry) <= radius.getAsDouble(); + } + else { + throw new UnsupportedOperationException("Unsupported spherical comparison operator: " + comparisonOperator); + } + } + throw new UnsupportedOperationException("Unsupported spherical spatial function: " + functionName); + } + private Set getSymbolReferences(Collection variables) { return variables.stream().map(VariableReferenceExpression::getName).map(SymbolReference::new).collect(toImmutableSet()); diff --git a/presto-main/src/main/java/com/facebook/presto/sql/planner/iterative/rule/ExtractSpatialJoins.java b/presto-main/src/main/java/com/facebook/presto/sql/planner/iterative/rule/ExtractSpatialJoins.java index 6326fb01473a6..c3b76272658ea 100644 --- a/presto-main/src/main/java/com/facebook/presto/sql/planner/iterative/rule/ExtractSpatialJoins.java +++ b/presto-main/src/main/java/com/facebook/presto/sql/planner/iterative/rule/ExtractSpatialJoins.java @@ -386,9 +386,10 @@ private static Result tryCreateSpatialJoin( RowExpression secondArgument = arguments.get(1); Type sphericalGeographyType = metadata.getType(SPHERICAL_GEOGRAPHY_TYPE_SIGNATURE); - // TODO check if we can just call arguments.getType().equals(sphericalGeographyType) as opposed to checking individual args if (firstArgument.getType().equals(sphericalGeographyType) || secondArgument.getType().equals(sphericalGeographyType)) { - return Result.empty(); + if (joinNode.getType() != INNER) { + return Result.empty(); + } } Set firstVariables = extractUnique(firstArgument);