From 0980ada83c29c5a00daaf53f1c2f8c1c78fc3033 Mon Sep 17 00:00:00 2001 From: adeet1 Date: Wed, 17 Jan 2024 16:30:45 +0000 Subject: [PATCH 1/3] GEOMESA-3308 Fix DWithin filter to compute distance between points across the antimeridian * Translate longitude by +/-360 and compare distances --- .../geomesa/filter/expression/FastDWithin.scala | 16 +++++++++++++--- .../filter/factory/FastFilterFactoryTest.scala | 16 ++++++++++++++++ 2 files changed, 29 insertions(+), 3 deletions(-) diff --git a/geomesa-filter/src/main/scala/org/locationtech/geomesa/filter/expression/FastDWithin.scala b/geomesa-filter/src/main/scala/org/locationtech/geomesa/filter/expression/FastDWithin.scala index a0de19c5151a..f8187749e06c 100644 --- a/geomesa-filter/src/main/scala/org/locationtech/geomesa/filter/expression/FastDWithin.scala +++ b/geomesa-filter/src/main/scala/org/locationtech/geomesa/filter/expression/FastDWithin.scala @@ -41,14 +41,24 @@ object FastDWithin { override def evaluate(obj: AnyRef): Boolean = { val geom = exp1.evaluate(obj).asInstanceOf[Geometry] - if (geom == null || envelope.distance(geom.getEnvelopeInternal) > maxDegrees) { false } else { + + if (geom == null) { + false + } else { val op = new DistanceOp(geometry, geom) op.distance <= minDegrees || { val Array(p1, p2) = op.nearestPoints() val calculator = calculators.get - calculator.setStartingGeographicPoint(p1.x, p1.y) + + calculator.setStartingGeographicPoint(p1.x - 360, p1.y) calculator.setDestinationGeographicPoint(p2.x, p2.y) - calculator.getOrthodromicDistance <= meters + val dist1 = calculator.getOrthodromicDistance + + calculator.setStartingGeographicPoint(p1.x + 360, p1.y) + calculator.setDestinationGeographicPoint(p2.x, p2.y) + val dist2 = calculator.getOrthodromicDistance + + Math.min(dist1, dist2) <= meters } } } diff --git a/geomesa-filter/src/test/scala/org/locationtech/geomesa/filter/factory/FastFilterFactoryTest.scala b/geomesa-filter/src/test/scala/org/locationtech/geomesa/filter/factory/FastFilterFactoryTest.scala index 833ef1913ed9..100c45e4883e 100644 --- a/geomesa-filter/src/test/scala/org/locationtech/geomesa/filter/factory/FastFilterFactoryTest.scala +++ b/geomesa-filter/src/test/scala/org/locationtech/geomesa/filter/factory/FastFilterFactoryTest.scala @@ -120,5 +120,21 @@ class FastFilterFactoryTest extends Specification { FastFilterFactory.toFilter(sft, "s < '100'").evaluate(sf) must beFalse FastFilterFactory.toFilter(sft, "s <= '100'").evaluate(sf) must beFalse } + + "calculate the distance between points horizontally across the antimeridian" >> { + val sft = SimpleFeatureTypes.createType("test", "geom:Point:srid=4326") + val feature = ScalaSimpleFeature.create(sft, "1", "POINT (178 0)") + val ecql = "dwithin('POINT (-179 0)', geom, 2500, kilometers)" + val fast = FastFilterFactory.toFilter(sft, ecql) + fast.evaluate(feature) must beTrue + } + + "calculate the distance between points diagonally across the antimeridian" >> { + val sft = SimpleFeatureTypes.createType("test", "geom:Point:srid=4326") + val feature = ScalaSimpleFeature.create(sft, "1", "POINT (178 3.1415)") + val ecql = "dwithin('POINT (-179 -2.7182)', geom, 2500, kilometers)" + val fast = FastFilterFactory.toFilter(sft, ecql) + fast.evaluate(feature) must beTrue + } } } From db811e3586285e0446f40b1e330449abeb9c4baa Mon Sep 17 00:00:00 2001 From: adeet1 Date: Thu, 18 Jan 2024 15:36:00 +0000 Subject: [PATCH 2/3] Add a unit test to calculate the distance spanning one degree of longitude --- .../geomesa/filter/factory/FastFilterFactoryTest.scala | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/geomesa-filter/src/test/scala/org/locationtech/geomesa/filter/factory/FastFilterFactoryTest.scala b/geomesa-filter/src/test/scala/org/locationtech/geomesa/filter/factory/FastFilterFactoryTest.scala index 100c45e4883e..eb1be2b026e9 100644 --- a/geomesa-filter/src/test/scala/org/locationtech/geomesa/filter/factory/FastFilterFactoryTest.scala +++ b/geomesa-filter/src/test/scala/org/locationtech/geomesa/filter/factory/FastFilterFactoryTest.scala @@ -136,5 +136,14 @@ class FastFilterFactoryTest extends Specification { val fast = FastFilterFactory.toFilter(sft, ecql) fast.evaluate(feature) must beTrue } + + // A degree of longitude is widest at the equator, with a distance of ~69 mi (~111 km) + "calculate the distance spanning one degree of longitude" >> { + val sft = SimpleFeatureTypes.createType("test", "geom:Point:srid=4326") + val feature = ScalaSimpleFeature.create(sft, "1", "POINT (0 0)") + val ecql = "dwithin('POINT (1 0)', geom, 112, kilometers)" + val fast = FastFilterFactory.toFilter(sft, ecql) + fast.evaluate(feature) must beTrue + } } } From ad8ef2390e6913d3f1b2c9c543d2ecee4f0ecfe0 Mon Sep 17 00:00:00 2001 From: adeet1 Date: Thu, 18 Jan 2024 15:46:59 +0000 Subject: [PATCH 3/3] Translate envelope bounds in the same way we do for the geodetic calculations --- .../filter/expression/FastDWithin.scala | 27 +++++++++++++++---- 1 file changed, 22 insertions(+), 5 deletions(-) diff --git a/geomesa-filter/src/main/scala/org/locationtech/geomesa/filter/expression/FastDWithin.scala b/geomesa-filter/src/main/scala/org/locationtech/geomesa/filter/expression/FastDWithin.scala index f8187749e06c..eaced28baf34 100644 --- a/geomesa-filter/src/main/scala/org/locationtech/geomesa/filter/expression/FastDWithin.scala +++ b/geomesa-filter/src/main/scala/org/locationtech/geomesa/filter/expression/FastDWithin.scala @@ -12,7 +12,7 @@ import org.geotools.referencing.GeodeticCalculator import org.locationtech.geomesa.filter.GeometryProcessing import org.locationtech.geomesa.utils.geotools.GeometryUtils import org.locationtech.geomesa.utils.geotools.converters.FastConverter -import org.locationtech.jts.geom.Geometry +import org.locationtech.jts.geom.{Envelope, Geometry} import org.locationtech.jts.operation.distance.DistanceOp import org.opengis.filter.FilterVisitor import org.opengis.filter.MultiValuedFilter.MatchAction @@ -35,14 +35,21 @@ object FastDWithin { */ class DWithinLiteral(exp1: Expression, exp2: Literal, distance: Double, units: String) extends DWithin { private val geometry = FastConverter.evaluate(exp2, classOf[Geometry]) - private val envelope = geometry.getEnvelopeInternal + private val envelope0 = geometry.getEnvelopeInternal private val meters = distance * GeometryProcessing.metersMultiplier(units) private val (minDegrees, maxDegrees) = GeometryUtils.distanceDegrees(geometry, meters) override def evaluate(obj: AnyRef): Boolean = { val geom = exp1.evaluate(obj).asInstanceOf[Geometry] + val dist0 = envelope0.distance(geom.getEnvelopeInternal) - if (geom == null) { + // Translate the longitude +/- 360 + val envelope1 = new Envelope(envelope0.getMinX - 360, envelope0.getMaxX - 360, envelope0.getMinY, envelope0.getMaxY) + val envelope2 = new Envelope(envelope0.getMinX + 360, envelope0.getMaxX + 360, envelope0.getMinY, envelope0.getMaxY) + val dist1 = envelope1.distance(geom.getEnvelopeInternal) + val dist2 = envelope2.distance(geom.getEnvelopeInternal) + + if (geom == null || (dist0 > maxDegrees) && (dist1 > maxDegrees) && (dist2 > maxDegrees)) { false } else { val op = new DistanceOp(geometry, geom) @@ -50,15 +57,25 @@ object FastDWithin { val Array(p1, p2) = op.nearestPoints() val calculator = calculators.get + calculator.setStartingGeographicPoint(p1.x, p1.y) + calculator.setDestinationGeographicPoint(p2.x, p2.y) + val dist0 = calculator.getOrthodromicDistance + if (dist0 <= meters) { + return true + } + + // Translate the longitude +/- 360 calculator.setStartingGeographicPoint(p1.x - 360, p1.y) calculator.setDestinationGeographicPoint(p2.x, p2.y) val dist1 = calculator.getOrthodromicDistance + if (dist1 <= meters) { + return true + } calculator.setStartingGeographicPoint(p1.x + 360, p1.y) calculator.setDestinationGeographicPoint(p2.x, p2.y) val dist2 = calculator.getOrthodromicDistance - - Math.min(dist1, dist2) <= meters + List(dist0, dist1, dist2).min <= meters } } }