search

benbovy / spherely

Manipulation and analysis of geometric objects on the sphere.
https://spherely.readthedocs.io/
BSD 3-Clause "New" or "Revised" License
119 stars 8 forks source link

Issue with exact equality and roundtrip of lat/lng coordinates #40

Open jorisvandenbossche opened 3 months ago

jorisvandenbossche commented 3 months ago

Something I ran into while testing https://github.com/benbovy/spherely/pull/20 (xref https://github.com/benbovy/spherely/pull/20#discussion_r1141150124), is that it is not that easy to pass an exact equality test with equals().

Illustration:

# construct a point through centroid to get something that is not exactly Point(1, 0)
>>> point = spherely.centroid(spherely.LineString([(0, 0), (0, 2)]))
>>> point
POINT (0.9999999999999998 0)

# reconstruct this point through its lat/long coordinates
>>> point2 = spherely.Point(spherely.get_y(point), spherely.get_x(point))
>>> point2
POINT (0.9999999999999998 0)

# those two points (with the same WKT repr and equal coordinates in degrees) are not considered equal
>>> spherely.equals(point, point2)
False

What I assume is happening here is that the roundtrip from S2Point to S2LatLng and back is not exact (the conversion from S2Point to LatLng essentially happens in the WKT repr or in get_x/get_y, and the conversion from S2LatLng to S2Point is essentially happening in the Point(..) constructor.

While that might be expected based on knowledge of the S2 internals (how a S2Point is stored, which is not as its lon/lat coords), I suppose this will not be expected in general by users? And it also gives some usability issues (also when testing)

For equals() specifically, it might be possible to have an option there to snap round coordinates while checking equality.

benbovy commented 3 months ago

The implementation of equals ultimately relies on s2geometry's S2BooleanOperation::Equals, which provides a precision option. Unfortunately Precision.SNAPPED is not yet implemented, but I guess it would still be possible to add some tolerance by passing a snap function to the Options constructor (e.g., IntLatLngSnapFunction or IdentitySnapFunction with a given snap_radius). We could do that in Spherely by customizing

https://github.com/benbovy/spherely/blob/38d3bf9df81e86140dacb962b573fc95e44f9316/src/predicates.cpp#L33

In terms of API, I see that shapely also provides equals_exact. I'm not sure how they behave exactly for point geometries, though. Would it make sense in spherely for equals to have some internal minimal tolerance set by default so that the example above works, and provide a separate equals_exact function that accepts a user-defined tolerance?