You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@sedona.apache.org by ji...@apache.org on 2023/08/29 04:12:56 UTC
[sedona] branch master updated: [SEDONA-286] Correctly handle poles and anti-meridian in optimized sphere/spheroid distance join (#989)
This is an automated email from the ASF dual-hosted git repository.
jiayu pushed a commit to branch master
in repository https://gitbox.apache.org/repos/asf/sedona.git
The following commit(s) were added to refs/heads/master by this push:
new 7ed65fe71 [SEDONA-286] Correctly handle poles and anti-meridian in optimized sphere/spheroid distance join (#989)
7ed65fe71 is described below
commit 7ed65fe7197f7c75c7e6c679cf234d6b244a3881
Author: Kristin Cowalcijk <bo...@wherobots.com>
AuthorDate: Tue Aug 29 12:12:50 2023 +0800
[SEDONA-286] Correctly handle poles and anti-meridian in optimized sphere/spheroid distance join (#989)
---
.../org/apache/sedona/common/sphere/Haversine.java | 38 +++
.../org/apache/sedona/common/FunctionsTest.java | 52 ++++
.../common/sphere/HaversineEnvelopeTest.java | 282 +++++++++++++++++++++
docs/api/sql/Optimizer.md | 2 +-
.../scala/org/apache/sedona/sql/UDF/Catalog.scala | 2 +-
.../sql/sedona_sql/expressions/Functions.scala | 2 +-
.../sql/sedona_sql/expressions/st_functions.scala | 4 +-
.../strategy/join/JoinQueryDetector.scala | 12 +-
.../sedona_sql/strategy/join/JoinedGeometry.scala | 38 ++-
.../sedona/sql/BroadcastIndexJoinSuite.scala | 5 -
.../sedona/sql/SphereDistanceJoinSuite.scala | 157 ++++++++++++
.../apache/sedona/sql/predicateJoinTestScala.scala | 5 -
12 files changed, 556 insertions(+), 43 deletions(-)
diff --git a/common/src/main/java/org/apache/sedona/common/sphere/Haversine.java b/common/src/main/java/org/apache/sedona/common/sphere/Haversine.java
index ba10fa966..2a2970fea 100644
--- a/common/src/main/java/org/apache/sedona/common/sphere/Haversine.java
+++ b/common/src/main/java/org/apache/sedona/common/sphere/Haversine.java
@@ -19,6 +19,7 @@
package org.apache.sedona.common.sphere;
import org.locationtech.jts.geom.Coordinate;
+import org.locationtech.jts.geom.Envelope;
import org.locationtech.jts.geom.Geometry;
import static java.lang.Math.atan2;
@@ -62,4 +63,41 @@ public class Haversine
{
return distance(geom1, geom2, 6371008.0);
}
+
+ /**
+ * Expand the given envelope on sphere by the given distance in meter.
+ *
+ * @param envelope the envelope to expand
+ * @param distance in meter
+ * @param sphereRadius radius of the sphere in meter
+ * @return expanded envelope
+ */
+ public static Envelope expandEnvelope(Envelope envelope, double distance, double sphereRadius)
+ {
+ // 10% buffer to get rid of false negatives
+ double scaleFactor = 1.1;
+ double latDeltaRadian = distance / sphereRadius;
+ double latDeltaDegree = Math.toDegrees(latDeltaRadian);
+ double newMinY = envelope.getMinY() - latDeltaDegree * scaleFactor;
+ double newMaxY = envelope.getMaxY() + latDeltaDegree * scaleFactor;
+ if (newMinY <= -90 || newMaxY >= 90) {
+ // The expanded envelope covers the pole, so its longitude range should be expanded to [-180, 180]
+ return new Envelope(-180, 180, Math.max(newMinY, -90), Math.min(newMaxY, 90));
+ }
+
+ double minLatRadian = Math.toRadians(newMinY);
+ double maxLatRadian = Math.toRadians(newMaxY);
+ double lonDeltaRadian = Math.max(Math.abs(distance / (sphereRadius * Math.cos(maxLatRadian))),
+ Math.abs(distance / (sphereRadius * Math.cos(minLatRadian))));
+ double lonDeltaDegree = Math.toDegrees(lonDeltaRadian);
+ double newMinX = envelope.getMinX() - lonDeltaDegree * scaleFactor;
+ double newMaxX = envelope.getMaxX() + lonDeltaDegree * scaleFactor;
+ if (newMinX <= -180 || newMaxX >= 180) {
+ // The expanded envelope crosses the anti-meridian. Expanding its longitude range to [-180, 180] is
+ // the best thing we can do, since we treat lon-lat as planar coordinates when running spatial join.
+ return new Envelope(-180, 180, newMinY, newMaxY);
+ } else {
+ return new Envelope(newMinX, newMaxX, newMinY, newMaxY);
+ }
+ }
}
diff --git a/common/src/test/java/org/apache/sedona/common/FunctionsTest.java b/common/src/test/java/org/apache/sedona/common/FunctionsTest.java
index 946949c54..b58cfbe5e 100644
--- a/common/src/test/java/org/apache/sedona/common/FunctionsTest.java
+++ b/common/src/test/java/org/apache/sedona/common/FunctionsTest.java
@@ -660,6 +660,32 @@ public class FunctionsTest {
p1 = GEOMETRY_FACTORY.createPoint(new Coordinate(113.914603, 22.308919));
p2 = GEOMETRY_FACTORY.createPoint(new Coordinate(-79.630556, 43.677223));
assertEquals(1.2548548944238186E7, Haversine.distance(p1, p2), 0.1);
+
+ // Crossing the anti-meridian
+ p1 = GEOMETRY_FACTORY.createPoint(new Coordinate(179.999, 0));
+ p2 = GEOMETRY_FACTORY.createPoint(new Coordinate(-179.999, 0));
+ assertTrue(Haversine.distance(p1, p2) < 300);
+ assertTrue(Haversine.distance(p2, p1) < 300);
+ p1 = GEOMETRY_FACTORY.createPoint(new Coordinate(179.999, 60));
+ p2 = GEOMETRY_FACTORY.createPoint(new Coordinate(-179.999, 60));
+ assertTrue(Haversine.distance(p1, p2) < 300);
+ assertTrue(Haversine.distance(p2, p1) < 300);
+ p1 = GEOMETRY_FACTORY.createPoint(new Coordinate(179.999, -60));
+ p2 = GEOMETRY_FACTORY.createPoint(new Coordinate(-179.999, -60));
+ assertTrue(Haversine.distance(p1, p2) < 300);
+ assertTrue(Haversine.distance(p2, p1) < 300);
+
+ // Crossing the North Pole
+ p1 = GEOMETRY_FACTORY.createPoint(new Coordinate(-60, 89.999));
+ p2 = GEOMETRY_FACTORY.createPoint(new Coordinate(120, 89.999));
+ assertTrue(Haversine.distance(p1, p2) < 300);
+ assertTrue(Haversine.distance(p2, p1) < 300);
+
+ // Crossing the South Pole
+ p1 = GEOMETRY_FACTORY.createPoint(new Coordinate(-60, -89.999));
+ p2 = GEOMETRY_FACTORY.createPoint(new Coordinate(120, -89.999));
+ assertTrue(Haversine.distance(p1, p2) < 300);
+ assertTrue(Haversine.distance(p2, p1) < 300);
}
@Test
@@ -694,6 +720,32 @@ public class FunctionsTest {
p1 = GEOMETRY_FACTORY.createPoint(new Coordinate(113.914603, 22.308919));
p2 = GEOMETRY_FACTORY.createPoint(new Coordinate(-79.630556, 43.677223));
assertEquals(1.2568775317073349E7, Spheroid.distance(p1, p2), 0.1);
+
+ // Crossing the anti-meridian
+ p1 = GEOMETRY_FACTORY.createPoint(new Coordinate(179.999, 0));
+ p2 = GEOMETRY_FACTORY.createPoint(new Coordinate(-179.999, 0));
+ assertTrue(Spheroid.distance(p1, p2) < 300);
+ assertTrue(Spheroid.distance(p2, p1) < 300);
+ p1 = GEOMETRY_FACTORY.createPoint(new Coordinate(179.999, 60));
+ p2 = GEOMETRY_FACTORY.createPoint(new Coordinate(-179.999, 60));
+ assertTrue(Spheroid.distance(p1, p2) < 300);
+ assertTrue(Spheroid.distance(p2, p1) < 300);
+ p1 = GEOMETRY_FACTORY.createPoint(new Coordinate(179.999, -60));
+ p2 = GEOMETRY_FACTORY.createPoint(new Coordinate(-179.999, -60));
+ assertTrue(Spheroid.distance(p1, p2) < 300);
+ assertTrue(Spheroid.distance(p2, p1) < 300);
+
+ // Crossing the North Pole
+ p1 = GEOMETRY_FACTORY.createPoint(new Coordinate(-60, 89.999));
+ p2 = GEOMETRY_FACTORY.createPoint(new Coordinate(120, 89.999));
+ assertTrue(Spheroid.distance(p1, p2) < 300);
+ assertTrue(Spheroid.distance(p2, p1) < 300);
+
+ // Crossing the South Pole
+ p1 = GEOMETRY_FACTORY.createPoint(new Coordinate(-60, -89.999));
+ p2 = GEOMETRY_FACTORY.createPoint(new Coordinate(120, -89.999));
+ assertTrue(Spheroid.distance(p1, p2) < 300);
+ assertTrue(Spheroid.distance(p2, p1) < 300);
}
@Test
diff --git a/common/src/test/java/org/apache/sedona/common/sphere/HaversineEnvelopeTest.java b/common/src/test/java/org/apache/sedona/common/sphere/HaversineEnvelopeTest.java
new file mode 100644
index 000000000..b087e83a6
--- /dev/null
+++ b/common/src/test/java/org/apache/sedona/common/sphere/HaversineEnvelopeTest.java
@@ -0,0 +1,282 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements. See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership. The ASF licenses this file
+ * to you 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 org.apache.sedona.common.sphere;
+
+import org.junit.Assert;
+import org.junit.Test;
+import org.locationtech.jts.geom.Coordinate;
+import org.locationtech.jts.geom.Envelope;
+import org.locationtech.jts.geom.GeometryFactory;
+import org.locationtech.jts.geom.Point;
+
+public class HaversineEnvelopeTest {
+ private static final int SPHERE_RADIUS = 6357000;
+ private static final GeometryFactory factory = new GeometryFactory();
+
+ @Test
+ public void testExpandPoint() {
+ // The null island
+ Envelope env = new Envelope(0, 0, 0, 0);
+ Envelope expandedEnv = Haversine.expandEnvelope(env, 500_000, SPHERE_RADIUS);
+ validateExpandedEnv(env, expandedEnv, 500_000);
+
+ // North semi-sphere
+ env = new Envelope(60, 60, 40, 40);
+ expandedEnv = Haversine.expandEnvelope(env, 500_000, SPHERE_RADIUS);
+ validateExpandedEnv(env, expandedEnv, 500_000);
+ Assert.assertFalse(coverEntireLongitudeRange(expandedEnv));
+
+ env = new Envelope(-60, -60, 70, 70);
+ expandedEnv = Haversine.expandEnvelope(env, 500_000, SPHERE_RADIUS);
+ validateExpandedEnv(env, expandedEnv, 500_000);
+ Assert.assertFalse(coverEntireLongitudeRange(expandedEnv));
+
+ // South semi-sphere
+ env = new Envelope(60, 60, -40, -40);
+ expandedEnv = Haversine.expandEnvelope(env, 500_000, SPHERE_RADIUS);
+ validateExpandedEnv(env, expandedEnv, 500_000);
+ Assert.assertFalse(coverEntireLongitudeRange(expandedEnv));
+
+ env = new Envelope(-60, -60, -70, -70);
+ expandedEnv = Haversine.expandEnvelope(env, 500_000, SPHERE_RADIUS);
+ validateExpandedEnv(env, expandedEnv, 500_000);
+ Assert.assertFalse(coverEntireLongitudeRange(expandedEnv));
+ }
+
+ @Test
+ public void testExpandPointCrossingAntiMeridian() {
+ // On the anti-meridian
+ Envelope env = new Envelope(-180, -180, 50, 50);
+ Envelope expandedEnv = Haversine.expandEnvelope(env, 500_000, SPHERE_RADIUS);
+ validateExpandedEnv(env, expandedEnv, 500_000);
+ Assert.assertTrue(coverEntireLongitudeRange(expandedEnv));
+
+ env = new Envelope(180, 180, -60, -60);
+ expandedEnv = Haversine.expandEnvelope(env, 500_000, SPHERE_RADIUS);
+ validateExpandedEnv(env, expandedEnv, 500_000);
+ Assert.assertTrue(coverEntireLongitudeRange(expandedEnv));
+
+ // Expand to cross the anti-meridian
+ env = new Envelope(-178, -178, -50, -50);
+ expandedEnv = Haversine.expandEnvelope(env, 500_000, SPHERE_RADIUS);
+ validateExpandedEnv(env, expandedEnv, 500_000);
+ Assert.assertTrue(coverEntireLongitudeRange(expandedEnv));
+
+ env = new Envelope(178, 178, 60, 60);
+ expandedEnv = Haversine.expandEnvelope(env, 500_000, SPHERE_RADIUS);
+ validateExpandedEnv(env, expandedEnv, 500_000);
+ Assert.assertTrue(coverEntireLongitudeRange(expandedEnv));
+ }
+
+ @Test
+ public void testExpandPointCoveringThePoles() {
+ // Point on the pole
+ Envelope env = new Envelope(0, 0, 90, 90);
+ Envelope expandedEnv = Haversine.expandEnvelope(env, 500_000, SPHERE_RADIUS);
+ validateExpandedEnv(env, expandedEnv, 500_000);
+ Assert.assertTrue(coverTheNorthPole(expandedEnv));
+
+ env = new Envelope(-120, -120, -90, -90);
+ expandedEnv = Haversine.expandEnvelope(env, 500_000, SPHERE_RADIUS);
+ validateExpandedEnv(env, expandedEnv, 500_000);
+ Assert.assertTrue(coverTheSouthPole(expandedEnv));
+
+ // Expand to cover the North Pole
+ env = new Envelope(-120, -120, 89, 89);
+ expandedEnv = Haversine.expandEnvelope(env, 500_000, SPHERE_RADIUS);
+ validateExpandedEnv(env, expandedEnv, 500_000);
+ Assert.assertTrue(coverTheNorthPole(expandedEnv));
+
+ // Expand to cover the South Pole
+ env = new Envelope(160, 160, -89, -89);
+ expandedEnv = Haversine.expandEnvelope(env, 500_000, SPHERE_RADIUS);
+ validateExpandedEnv(env, expandedEnv, 500_000);
+ Assert.assertTrue(coverTheSouthPole(expandedEnv));
+ }
+
+ @Test
+ public void testExpandPointCoveringTheGlobe() {
+ Envelope env = new Envelope(10, 10, 20, 20);
+ Envelope expandedEnv = Haversine.expandEnvelope(env, 2 * SPHERE_RADIUS * Math.PI, SPHERE_RADIUS);
+ Assert.assertTrue(coverTheNorthPole(expandedEnv));
+ Assert.assertTrue(coverTheSouthPole(expandedEnv));
+ }
+
+ @Test
+ public void testExpandEnvelope() {
+ Envelope env = new Envelope(-10, 10, -20, 20);
+ Envelope expandedEnv = Haversine.expandEnvelope(env, 500_000, SPHERE_RADIUS);
+ validateExpandedEnv(env, expandedEnv, 500_000);
+
+ env = new Envelope(10, 20, 30, 40);
+ expandedEnv = Haversine.expandEnvelope(env, 500_000, SPHERE_RADIUS);
+ validateExpandedEnv(env, expandedEnv, 500_000);
+
+ env = new Envelope(-20, 10, 60, 70);
+ expandedEnv = Haversine.expandEnvelope(env, 500_000, SPHERE_RADIUS);
+ validateExpandedEnv(env, expandedEnv, 500_000);
+
+ env = new Envelope(10, 20, -40, -30);
+ expandedEnv = Haversine.expandEnvelope(env, 500_000, SPHERE_RADIUS);
+ validateExpandedEnv(env, expandedEnv, 500_000);
+
+ env = new Envelope(10, 20, -40, -30);
+ expandedEnv = Haversine.expandEnvelope(env, 500_000, SPHERE_RADIUS);
+ validateExpandedEnv(env, expandedEnv, 500_000);
+
+ env = new Envelope(-20, 10, -70, -60);
+ expandedEnv = Haversine.expandEnvelope(env, 500_000, SPHERE_RADIUS);
+ validateExpandedEnv(env, expandedEnv, 500_000);
+ }
+
+ @Test
+ public void testExpandEnvelopeCrossingAntiMeridian() {
+ Envelope env = new Envelope(-200, 160, -20, 20);
+ Envelope expandedEnv = Haversine.expandEnvelope(env, 500_000, SPHERE_RADIUS);
+ validateExpandedEnv(env, expandedEnv, 500_000);
+ Assert.assertTrue(coverEntireLongitudeRange(expandedEnv));
+
+ env = new Envelope(160, 200, -20, 20);
+ expandedEnv = Haversine.expandEnvelope(env, 500_000, SPHERE_RADIUS);
+ validateExpandedEnv(env, expandedEnv, 500_000);
+ Assert.assertTrue(coverEntireLongitudeRange(expandedEnv));
+
+ env = new Envelope(-179, -160, -40, 60);
+ expandedEnv = Haversine.expandEnvelope(env, 500_000, SPHERE_RADIUS);
+ validateExpandedEnv(env, expandedEnv, 500_000);
+ Assert.assertTrue(coverEntireLongitudeRange(expandedEnv));
+
+ env = new Envelope(160, 179, -40, 60);
+ expandedEnv = Haversine.expandEnvelope(env, 500_000, SPHERE_RADIUS);
+ validateExpandedEnv(env, expandedEnv, 500_000);
+ Assert.assertTrue(coverEntireLongitudeRange(expandedEnv));
+ }
+
+ @Test
+ public void testExpandEnvelopeCoveringPoles() {
+ Envelope env = new Envelope(40, 50, 89, 90);
+ Envelope expandedEnv = Haversine.expandEnvelope(env, 500_000, SPHERE_RADIUS);
+ validateExpandedEnv(env, expandedEnv, 500_000);
+ Assert.assertTrue(coverTheNorthPole(expandedEnv));
+
+ env = new Envelope(40, 50, -90, -89);
+ expandedEnv = Haversine.expandEnvelope(env, 500_000, SPHERE_RADIUS);
+ validateExpandedEnv(env, expandedEnv, 500_000);
+ Assert.assertTrue(coverTheSouthPole(expandedEnv));
+
+ env = new Envelope(40, 50, 60, 89);
+ expandedEnv = Haversine.expandEnvelope(env, 500_000, SPHERE_RADIUS);
+ validateExpandedEnv(env, expandedEnv, 500_000);
+ Assert.assertTrue(coverTheNorthPole(expandedEnv));
+
+ env = new Envelope(40, 50, -89, -80);
+ expandedEnv = Haversine.expandEnvelope(env, 500_000, SPHERE_RADIUS);
+ validateExpandedEnv(env, expandedEnv, 500_000);
+ Assert.assertTrue(coverTheSouthPole(expandedEnv));
+ }
+
+ @Test
+ public void testExpandEnvelopeCoverTheGlobe() {
+ Envelope env = new Envelope(-178, 178, -89, 89);
+ Envelope expandedEnv = Haversine.expandEnvelope(env, 500_000, SPHERE_RADIUS);
+ validateExpandedEnv(env, expandedEnv, 500_000);
+ Assert.assertTrue(coverTheNorthPole(expandedEnv));
+ Assert.assertTrue(coverTheSouthPole(expandedEnv));
+
+ env = new Envelope(10, 20, -89, 89);
+ expandedEnv = Haversine.expandEnvelope(env, 500_000, SPHERE_RADIUS);
+ validateExpandedEnv(env, expandedEnv, 500_000);
+ Assert.assertTrue(coverTheNorthPole(expandedEnv));
+ Assert.assertTrue(coverTheSouthPole(expandedEnv));
+ }
+
+ private void validateExpandedEnv(Envelope env, Envelope expandedEnv, double distance) {
+ if (env.getMinX() >= -180 && env.getMaxX() <= 180 && env.getMinY() >= -90 && env.getMinY() <= 90) {
+ Assert.assertTrue(expandedEnv.contains(env));
+ }
+
+ // Validate that points on the edge of the expanded envelope are further than the distance
+ if (!coverEntireLongitudeRange(expandedEnv)) {
+ // Upper left corner of original envelope
+ Coordinate original = new Coordinate(env.getMinX(), env.getMaxY());
+ Coordinate expanded = new Coordinate(expandedEnv.getMinX(), env.getMaxY());
+ validateDistance(original, expanded, distance);
+
+ // Left middle of original envelope
+ original = new Coordinate(env.getMinX(), env.centre().y);
+ expanded = new Coordinate(expandedEnv.getMinX(), env.centre().y);
+ validateDistance(original, expanded, distance);
+
+ // Lower left corner of original envelope
+ original = new Coordinate(env.getMinX(), env.getMinY());
+ expanded = new Coordinate(expandedEnv.getMinX(), env.getMinY());
+ validateDistance(original, expanded, distance);
+
+ // Upper right of the original envelope
+ original = new Coordinate(env.getMaxX(), env.getMaxY());
+ expanded = new Coordinate(expandedEnv.getMaxX(), env.getMaxY());
+ validateDistance(original, expanded, distance);
+
+ // Right middle of the original envelope
+ original = new Coordinate(env.getMaxX(), env.centre().y);
+ expanded = new Coordinate(expandedEnv.getMaxX(), env.centre().y);
+ validateDistance(original, expanded, distance);
+
+ // Lower right of the original envelope
+ original = new Coordinate(env.getMaxX(), env.getMinY());
+ expanded = new Coordinate(expandedEnv.getMaxX(), env.getMinY());
+ validateDistance(original, expanded, distance);
+ }
+
+ if (!coverTheNorthPole(expandedEnv)) {
+ // Validate the center of top edge
+ Coordinate original = new Coordinate(env.centre().x, env.getMaxY());
+ Coordinate expanded = new Coordinate(env.centre().x, expandedEnv.getMaxY());
+ validateDistance(original, expanded, distance);
+ }
+
+ if (!coverTheSouthPole(expandedEnv)) {
+ // Validate the center of bottom edge
+ Coordinate original = new Coordinate(env.centre().x, env.getMinY());
+ Coordinate expanded = new Coordinate(env.centre().x, expandedEnv.getMinY());
+ validateDistance(original, expanded, distance);
+ }
+ }
+
+ private boolean coverEntireLongitudeRange(Envelope env) {
+ return env.getMinX() <= -180 && env.getMaxX() >= 180;
+ }
+
+ private boolean coverTheNorthPole(Envelope env) {
+ return coverEntireLongitudeRange(env) && env.getMaxY() >= 90;
+ }
+
+ private boolean coverTheSouthPole(Envelope env) {
+ return coverEntireLongitudeRange(env) && env.getMinY() <= -90;
+ }
+
+ private void validateDistance(Coordinate original, Coordinate expanded, double distance) {
+ Point ptOriginal = factory.createPoint(original);
+ Point ptExpanded = factory.createPoint(expanded);
+ double actualSphereDistance = Haversine.distance(ptOriginal, ptExpanded);
+ Assert.assertTrue(actualSphereDistance >= distance);
+ double actualSpheroidDistance = Spheroid.distance(ptOriginal, ptExpanded);
+ Assert.assertTrue(actualSpheroidDistance >= distance);
+ }
+}
diff --git a/docs/api/sql/Optimizer.md b/docs/api/sql/Optimizer.md
index a95b287ae..d96bcaee8 100644
--- a/docs/api/sql/Optimizer.md
+++ b/docs/api/sql/Optimizer.md
@@ -118,7 +118,7 @@ WHERE ST_DistanceSpheroid(pointdf1.pointshape1,pointdf2.pointshape2) <= 2
```
!!!warning
- If you use `ST_DistanceSpheroid ` or `ST_DistanceSphere` as the predicate, the unit of the distance is meter. Currently, distance join with geodesic distance calculators work best for point data. For non-point data, it only considers their centroids. The distance join algorithm internally uses an approximate distance buffer which might lead to inaccurate results if your data is close to the poles or antimeridian.
+ If you use `ST_DistanceSpheroid ` or `ST_DistanceSphere` as the predicate, the unit of the distance is meter. Currently, distance join with geodesic distance calculators work best for point data. For non-point data, it only considers their centroids.
## Broadcast index join
diff --git a/sql/common/src/main/scala/org/apache/sedona/sql/UDF/Catalog.scala b/sql/common/src/main/scala/org/apache/sedona/sql/UDF/Catalog.scala
index ee104a411..91c5bf50c 100644
--- a/sql/common/src/main/scala/org/apache/sedona/sql/UDF/Catalog.scala
+++ b/sql/common/src/main/scala/org/apache/sedona/sql/UDF/Catalog.scala
@@ -152,7 +152,7 @@ object Catalog {
function[ST_Split](),
function[ST_S2CellIDs](),
function[ST_GeometricMedian](1e-6, 1000, false),
- function[ST_DistanceSphere](6371008.0),
+ function[ST_DistanceSphere](),
function[ST_DistanceSpheroid](),
function[ST_AreaSpheroid](),
function[ST_LengthSpheroid](),
diff --git a/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/Functions.scala b/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/Functions.scala
index f063919fc..0bfedc38d 100644
--- a/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/Functions.scala
+++ b/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/Functions.scala
@@ -971,7 +971,7 @@ case class ST_GeometricMedian(inputExpressions: Seq[Expression])
}
case class ST_DistanceSphere(inputExpressions: Seq[Expression])
- extends InferredExpression(inferrableFunction3(Haversine.distance)) {
+ extends InferredExpression(inferrableFunction2(Haversine.distance), inferrableFunction3(Haversine.distance)) {
protected def withNewChildrenInternal(newChildren: IndexedSeq[Expression]) = {
copy(inputExpressions = newChildren)
diff --git a/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/st_functions.scala b/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/st_functions.scala
index db993d483..b6b6caeae 100644
--- a/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/st_functions.scala
+++ b/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/st_functions.scala
@@ -302,8 +302,8 @@ object st_functions extends DataFrameAPI {
def ST_GeometricMedian(geometry: Column, tolerance: Column, maxIter: Column, failIfNotConverged: Column): Column = wrapExpression[ST_GeometricMedian](geometry, tolerance, maxIter, failIfNotConverged)
def ST_GeometricMedian(geometry: String, tolerance: Double, maxIter: Int, failIfNotConverged: Boolean): Column = wrapExpression[ST_GeometricMedian](geometry, tolerance, maxIter, failIfNotConverged)
- def ST_DistanceSphere(a: Column, b: Column): Column = wrapExpression[ST_DistanceSphere](a, b, 6371008.0)
- def ST_DistanceSphere(a: String, b: String): Column = wrapExpression[ST_DistanceSphere](a, b, 6371008.0)
+ def ST_DistanceSphere(a: Column, b: Column): Column = wrapExpression[ST_DistanceSphere](a, b)
+ def ST_DistanceSphere(a: String, b: String): Column = wrapExpression[ST_DistanceSphere](a, b)
def ST_DistanceSphere(a: Column, b: Column, c: Column): Column = wrapExpression[ST_DistanceSphere](a, b, c)
def ST_DistanceSphere(a: String, b: String, c: Double): Column = wrapExpression[ST_DistanceSphere](a, b, c)
diff --git a/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/strategy/join/JoinQueryDetector.scala b/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/strategy/join/JoinQueryDetector.scala
index 4c65cef08..ca52f6189 100644
--- a/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/strategy/join/JoinQueryDetector.scala
+++ b/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/strategy/join/JoinQueryDetector.scala
@@ -147,17 +147,17 @@ class JoinQueryDetector(sparkSession: SparkSession) extends Strategy {
case Some(And(_, LessThan(ST_Distance(Seq(leftShape, rightShape)), distance))) =>
Some(JoinQueryDetection(left, right, leftShape, rightShape, SpatialPredicate.INTERSECTS, false, condition, Some(distance)))
// ST_DistanceSphere
- case Some(LessThanOrEqual(ST_DistanceSphere(Seq(leftShape, rightShape, radius)), distance)) =>
+ case Some(LessThanOrEqual(ST_DistanceSphere(Seq(leftShape, rightShape)), distance)) =>
Some(JoinQueryDetection(left, right, leftShape, rightShape, SpatialPredicate.INTERSECTS, true, condition, Some(distance)))
- case Some(And(LessThanOrEqual(ST_DistanceSphere(Seq(leftShape, rightShape, radius)), distance), _)) =>
+ case Some(And(LessThanOrEqual(ST_DistanceSphere(Seq(leftShape, rightShape)), distance), _)) =>
Some(JoinQueryDetection(left, right, leftShape, rightShape, SpatialPredicate.INTERSECTS, true, condition, Some(distance)))
- case Some(And(_, LessThanOrEqual(ST_DistanceSphere(Seq(leftShape, rightShape, radius)), distance))) =>
+ case Some(And(_, LessThanOrEqual(ST_DistanceSphere(Seq(leftShape, rightShape)), distance))) =>
Some(JoinQueryDetection(left, right, leftShape, rightShape, SpatialPredicate.INTERSECTS, true, condition, Some(distance)))
- case Some(LessThan(ST_DistanceSphere(Seq(leftShape, rightShape, radius)), distance)) =>
+ case Some(LessThan(ST_DistanceSphere(Seq(leftShape, rightShape)), distance)) =>
Some(JoinQueryDetection(left, right, leftShape, rightShape, SpatialPredicate.INTERSECTS, true, condition, Some(distance)))
- case Some(And(LessThan(ST_DistanceSphere(Seq(leftShape, rightShape, radius)), distance), _)) =>
+ case Some(And(LessThan(ST_DistanceSphere(Seq(leftShape, rightShape)), distance), _)) =>
Some(JoinQueryDetection(left, right, leftShape, rightShape, SpatialPredicate.INTERSECTS, true, condition, Some(distance)))
- case Some(And(_, LessThan(ST_DistanceSphere(Seq(leftShape, rightShape, radius)), distance))) =>
+ case Some(And(_, LessThan(ST_DistanceSphere(Seq(leftShape, rightShape)), distance))) =>
Some(JoinQueryDetection(left, right, leftShape, rightShape, SpatialPredicate.INTERSECTS, true, condition, Some(distance)))
// ST_DistanceSpheroid
case Some(LessThanOrEqual(ST_DistanceSpheroid(Seq(leftShape, rightShape)), distance)) =>
diff --git a/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/strategy/join/JoinedGeometry.scala b/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/strategy/join/JoinedGeometry.scala
index 745abd828..77383c0c0 100644
--- a/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/strategy/join/JoinedGeometry.scala
+++ b/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/strategy/join/JoinedGeometry.scala
@@ -18,6 +18,7 @@
*/
package org.apache.spark.sql.sedona_sql.strategy.join
+import org.apache.sedona.common.sphere.Haversine
import org.locationtech.jts.geom.{Envelope, Geometry}
/**
@@ -33,33 +34,26 @@ object JoinedGeometry {
*/
def geometryToExpandedEnvelope(geom: Geometry, distance: Double, isGeography: Boolean): Geometry = {
val envelope = geom.getEnvelopeInternal.copy()
- // Here we use the polar radius of the spheroid as the radius of the sphere, so that the expanded
- // envelope will work for both spherical and spheroidal distances.
- expandEnvelope(envelope, distance, 6357000.0, isGeography)
- geom.getFactory.toGeometry(envelope)
+ val expandedEnvelope = if (isGeography) {
+ expandEnvelopeForGeography(envelope, distance)
+ } else {
+ val newEnvelope = envelope.copy()
+ newEnvelope.expandBy(distance)
+ newEnvelope
+ }
+ geom.getFactory.toGeometry(expandedEnvelope)
}
/**
* Expand the given envelope by the given distance in meter.
- * For geography, we expand the envelope by the given distance in both longitude and latitude.
*
- * @param envelope the envelope to expand
- * @param distance in meter
- * @param sphereRadius in meter
- * @param isGeography whether the envelope is on a sphere
+ * @param envelope the envelope to expand
+ * @param distance in meter
*/
- private def expandEnvelope(envelope: Envelope, distance: Double, sphereRadius: Double, isGeography: Boolean): Unit = {
- if (isGeography) {
- val scaleFactor = 1.1 // 10% buffer to get rid of false negatives
- val latRadian = Math.toRadians((envelope.getMinY + envelope.getMaxY) / 2.0)
- val latDeltaRadian = distance / sphereRadius;
- val latDeltaDegree = Math.toDegrees(latDeltaRadian)
- val lonDeltaRadian = Math.max(Math.abs(distance / (sphereRadius * Math.cos(latRadian + latDeltaRadian))),
- Math.abs(distance / (sphereRadius * Math.cos(latRadian - latDeltaRadian))))
- val lonDeltaDegree = Math.toDegrees(lonDeltaRadian)
- envelope.expandBy(latDeltaDegree * scaleFactor, lonDeltaDegree * scaleFactor)
- } else {
- envelope.expandBy(distance)
- }
+ private def expandEnvelopeForGeography(envelope: Envelope, distance: Double): Envelope = {
+ // Here we use the polar radius of the spheroid as the radius of the sphere, so that the expanded
+ // envelope will work for both spherical and spheroidal distances.
+ val sphereRadius = 6357000.0
+ Haversine.expandEnvelope(envelope, distance, sphereRadius)
}
}
diff --git a/sql/common/src/test/scala/org/apache/sedona/sql/BroadcastIndexJoinSuite.scala b/sql/common/src/test/scala/org/apache/sedona/sql/BroadcastIndexJoinSuite.scala
index 6971a6a56..489ece9cf 100644
--- a/sql/common/src/test/scala/org/apache/sedona/sql/BroadcastIndexJoinSuite.scala
+++ b/sql/common/src/test/scala/org/apache/sedona/sql/BroadcastIndexJoinSuite.scala
@@ -343,11 +343,6 @@ class BroadcastIndexJoinSuite extends TestBaseScala {
pointDf2.alias("pointDf2"), expr(s"ST_DistanceSphere(pointDf1.pointshape, pointDf2.pointshape) <= $distance"))
assert(distanceJoinDf.queryExecution.sparkPlan.collect { case p: BroadcastIndexJoinExec => p }.size === 1)
assert(distanceJoinDf.count() == expected)
-
- distanceJoinDf = broadcast(pointDf1).alias("pointDf1").join(
- pointDf2.alias("pointDf2"), expr(s"ST_DistanceSphere(pointDf1.pointshape, pointDf2.pointshape, 6371008.0) < $distance"))
- assert(distanceJoinDf.queryExecution.sparkPlan.collect { case p: BroadcastIndexJoinExec => p }.size === 1)
- assert(distanceJoinDf.count() == expected)
})
}
diff --git a/sql/common/src/test/scala/org/apache/sedona/sql/SphereDistanceJoinSuite.scala b/sql/common/src/test/scala/org/apache/sedona/sql/SphereDistanceJoinSuite.scala
new file mode 100644
index 000000000..40e38e5fc
--- /dev/null
+++ b/sql/common/src/test/scala/org/apache/sedona/sql/SphereDistanceJoinSuite.scala
@@ -0,0 +1,157 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements. See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership. The ASF licenses this file
+ * to you 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 org.apache.sedona.sql
+
+import org.apache.sedona.common.sphere.{Haversine, Spheroid}
+import org.apache.spark.sql.DataFrame
+import org.apache.spark.sql.sedona_sql.strategy.join.{BroadcastIndexJoinExec, DistanceJoinExec}
+import org.locationtech.jts.geom.{Coordinate, Geometry, GeometryFactory}
+import org.scalatest.prop.TableDrivenPropertyChecks
+
+import scala.util.Random
+
+class SphereDistanceJoinSuite extends TestBaseScala with TableDrivenPropertyChecks {
+ private val spatialJoinPartitionSideConfKey = "sedona.join.spatitionside"
+ private val spatialJoinPartitionSide = sparkSession.sparkContext.getConf.get(spatialJoinPartitionSideConfKey, "left")
+ private val factory = new GeometryFactory()
+
+ private val testData1: Seq[(Int, Double, Geometry)] = generateTestData()
+ private val testData2: Seq[(Int, Double, Geometry)] = generateTestData()
+
+ override def beforeAll(): Unit = {
+ super.beforeAll()
+ prepareTempViewsForTestData()
+ }
+
+ describe("Sedona-SQL Spatial Join Test") {
+ val joinConditions = Table("join condition",
+ "ST_DistanceSphere(df1.geom, df2.geom) < 2000000",
+ "ST_DistanceSphere(df2.geom, df1.geom) < 2000000",
+ "ST_DistanceSpheroid(df1.geom, df2.geom) < 2000000",
+ "ST_DistanceSpheroid(df2.geom, df1.geom) < 2000000",
+
+ "ST_DistanceSphere(df1.geom, df2.geom) < df1.dist",
+ "ST_DistanceSphere(df2.geom, df1.geom) < df1.dist",
+ "ST_DistanceSpheroid(df1.geom, df2.geom) < df1.dist",
+ "ST_DistanceSpheroid(df2.geom, df1.geom) < df1.dist",
+
+ "ST_DistanceSphere(df1.geom, df2.geom) < df2.dist",
+ "ST_DistanceSphere(df2.geom, df1.geom) < df2.dist",
+ "ST_DistanceSpheroid(df1.geom, df2.geom) < df2.dist",
+ "ST_DistanceSpheroid(df2.geom, df1.geom) < df2.dist"
+ )
+
+ try {
+ forAll(joinConditions) { joinCondition =>
+ val expected = buildExpectedResult(joinCondition)
+ it(s"sphere distance join ON $joinCondition, with left side as dominant side") {
+ sparkSession.sparkContext.getConf.set(spatialJoinPartitionSideConfKey, "left")
+ val result = sparkSession.sql(s"SELECT df1.id, df2.id FROM df1 JOIN df2 ON $joinCondition")
+ verifyResult(expected, result)
+ }
+ it(s"sphere distance join ON $joinCondition, with right side as dominant side") {
+ sparkSession.sparkContext.getConf.set(spatialJoinPartitionSideConfKey, "right")
+ val result = sparkSession.sql(s"SELECT df1.id, df2.id FROM df1 JOIN df2 ON $joinCondition")
+ verifyResult(expected, result)
+ }
+ it(s"sphere distance ON $joinCondition, broadcast df1") {
+ val result = sparkSession.sql(s"SELECT /*+ BROADCAST(df1) */ df1.id, df2.id FROM df1 JOIN df2 ON $joinCondition")
+ verifyResult(expected, result)
+ }
+ it(s"sphere distance ON $joinCondition, broadcast df2") {
+ val result = sparkSession.sql(s"SELECT /*+ BROADCAST(df2) */ df1.id, df2.id FROM df1 JOIN df2 ON $joinCondition")
+ verifyResult(expected, result)
+ }
+ }
+ } finally {
+ sparkSession.sparkContext.getConf.set(spatialJoinPartitionSideConfKey, spatialJoinPartitionSide)
+ }
+ }
+
+ describe("Sphere distance join with custom sphere radius") {
+ it("do not optimize distance join with custom sphere radius") {
+ val df = sparkSession.sql("SELECT df1.id, df2.id FROM df1 JOIN df2 ON ST_DistanceSphere(df1.geom, df2.geom, 100) > 10")
+ assert(!isUsingOptimizedSpatialJoin(df))
+ }
+ }
+
+ private def prepareTempViewsForTestData(): Unit = {
+ import sparkSession.implicits._
+ testData1.toDF("id", "dist", "geom").createOrReplaceTempView("df1")
+ testData2.toDF("id", "dist", "geom").createOrReplaceTempView("df2")
+ }
+
+ private def generateTestData(): Seq[(Int, Double, Geometry)] = {
+ val geometries = (-180 to 180 by 10).flatMap { x =>
+ (-80 to 80 by 10).map { y =>
+ factory.createPoint(new Coordinate(x, y))
+ }
+ } ++ Seq(factory.createPoint(new Coordinate(0, -90)), factory.createPoint(new Coordinate(0, 90)))
+ val rand = new Random()
+ geometries.zipWithIndex.map { case (geom, idx) =>
+ (idx, 110000 + 2000000 * rand.nextDouble, geom)
+ }
+ }
+
+ private def buildExpectedResult(joinCondition: String): Seq[(Int, Int)] = {
+ val evaluate = joinCondition match {
+ case "ST_DistanceSphere(df1.geom, df2.geom) < 2000000" |
+ "ST_DistanceSphere(df2.geom, df1.geom) < 2000000" =>
+ (g1: Geometry, d1: Double, g2: Geometry, d2: Double) => Haversine.distance(g1, g2) < 2000000
+ case "ST_DistanceSphere(df1.geom, df2.geom) < df1.dist" |
+ "ST_DistanceSphere(df2.geom, df1.geom) < df1.dist" =>
+ (g1: Geometry, d1: Double, g2: Geometry, d2: Double) => Haversine.distance(g1, g2) < d1
+ case "ST_DistanceSphere(df1.geom, df2.geom) < df2.dist" |
+ "ST_DistanceSphere(df2.geom, df1.geom) < df2.dist" =>
+ (g1: Geometry, d1: Double, g2: Geometry, d2: Double) => Haversine.distance(g1, g2) < d2
+ case "ST_DistanceSpheroid(df1.geom, df2.geom) < 2000000" |
+ "ST_DistanceSpheroid(df2.geom, df1.geom) < 2000000" =>
+ (g1: Geometry, d1: Double, g2: Geometry, d2: Double) => Spheroid.distance(g1, g2) < 2000000
+ case "ST_DistanceSpheroid(df1.geom, df2.geom) < df1.dist" |
+ "ST_DistanceSpheroid(df2.geom, df1.geom) < df1.dist" =>
+ (g1: Geometry, d1: Double, g2: Geometry, d2: Double) => Spheroid.distance(g1, g2) < d1
+ case "ST_DistanceSpheroid(df1.geom, df2.geom) < df2.dist" |
+ "ST_DistanceSpheroid(df2.geom, df1.geom) < df2.dist" =>
+ (g1: Geometry, d1: Double, g2: Geometry, d2: Double) => Spheroid.distance(g1, g2) < d2
+ }
+ testData1.flatMap { case (id1, dist1, geom1) =>
+ testData2.flatMap { case (id2, dist2, geom2) =>
+ if (evaluate(geom1, dist1, geom2, dist2)) {
+ Some((id1, id2))
+ } else {
+ None
+ }
+ }
+ }
+ }
+
+ private def verifyResult(expected: Seq[(Int, Int)], result: DataFrame): Unit = {
+ isUsingOptimizedSpatialJoin(result)
+ val actual = result.collect().map(row => (row.getInt(0), row.getInt(1))).sorted
+ assert(actual.nonEmpty)
+ assert(actual === expected)
+ }
+
+ private def isUsingOptimizedSpatialJoin(df: DataFrame): Boolean = {
+ df.queryExecution.executedPlan.collect {
+ case _: BroadcastIndexJoinExec |
+ _: DistanceJoinExec => true
+ }.nonEmpty
+ }
+}
diff --git a/sql/common/src/test/scala/org/apache/sedona/sql/predicateJoinTestScala.scala b/sql/common/src/test/scala/org/apache/sedona/sql/predicateJoinTestScala.scala
index 001f91837..0ac021c55 100644
--- a/sql/common/src/test/scala/org/apache/sedona/sql/predicateJoinTestScala.scala
+++ b/sql/common/src/test/scala/org/apache/sedona/sql/predicateJoinTestScala.scala
@@ -444,11 +444,6 @@ class predicateJoinTestScala extends TestBaseScala {
pointDf2.alias("pointDf2"), expr(s"ST_DistanceSphere(pointDf1.pointshape, pointDf2.pointshape) < $distance"))
assert(distanceJoinDf.queryExecution.sparkPlan.collect { case p: DistanceJoinExec => p }.size === 1)
assert(distanceJoinDf.count() == expected)
-
- distanceJoinDf = pointDf1.alias("pointDf1").join(
- pointDf2.alias("pointDf2"), expr(s"ST_DistanceSphere(pointDf1.pointshape, pointDf2.pointshape, 6371008.0) < $distance"))
- assert(distanceJoinDf.queryExecution.sparkPlan.collect { case p: DistanceJoinExec => p }.size === 1)
- assert(distanceJoinDf.count() == expected)
})
}