Skip to content

GEOMESA-3308 Fix DWithin filter to compute distance between points across the antimeridian #3028

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
adeet1 wants to merge 3 commits into
base: main
Choose a base branch
from
Open

GEOMESA-3308 Fix DWithin filter to compute distance between points across the antimeridian #3028

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
0980ada
GEOMESA-3308 Fix DWithin filter to compute distance between points ac…
adeet1 Jan 17, 2024
db811e3
Add a unit test to calculate the distance spanning one degree of long…
adeet1 Jan 18, 2024
ad8ef23
Translate envelope bounds in the same way we do for the geodetic calc…
adeet1 Jan 18, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
35 changes: 31 additions & 4 deletions geomesa-filter/src/main/scala/org/locationtech/geomesa/filter/expression/FastDWithin.scala
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -35,20 +35,47 @@ 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]
if (geom == null || envelope.distance(geom.getEnvelopeInternal) > maxDegrees) { false } else {
val dist0 = envelope0.distance(geom.getEnvelopeInternal)

// 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)
Copy link
Contributor Author

@adeet1 adeet1 Jan 18, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The variable names geometry and geom are too similar and thus confusing. I think we should change these to something more descriptive.

op.distance <= minDegrees || {
val Array(p1, p2) = op.nearestPoints()
val calculator = calculators.get

calculator.setStartingGeographicPoint(p1.x, p1.y)
calculator.setDestinationGeographicPoint(p2.x, p2.y)
calculator.getOrthodromicDistance <= meters
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
Copy link
Contributor

@elahrvivaz elahrvivaz Jan 18, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

as another optimization here, we could return early if dist1 <= meters and not calculate dist2. but we also need to calculate dist0, which would be the unmodified point. I'm a little surprised none of the existing unit tests failed with this change :( before changing anything, could you add a test for (0 0) to (1 0) (or something similar) and verify that it fails? possibly the calculator is accounting for the anti-meridian wrapping?

Copy link
Contributor Author

@adeet1 adeet1 Jan 18, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I added a unit test (db811e3), but it passes, which I'm not sure is the intended behavior.

However, I think the calculator computes the shortest distance between two points on the Earth's surface, so it's likely that it already accounts for the anti-meridian wrapping.

List(dist0, dist1, dist2).min <= meters
}
}
}
Expand Down
25 changes: 25 additions & 0 deletions ...filter/src/test/scala/org/locationtech/geomesa/filter/factory/FastFilterFactoryTest.scala
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -120,5 +120,30 @@ 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
}

// 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
}
}
}