Handling voyages crossing the antimeridian

[tjansson60]

First of all, thanks for the very nice library. I really enjoy it. I think I have found a bug, which I have detailed below.

I am using the polygon shoreline set GSHHS_c_L1 and I am trying to calculate the distance between the the two points:

p1 = (179, 0)
p2 = (-179, 0)

The true distance between the points is probably something close to the geopy.distance.greater_circle distance of 120 nm as the points are very close to each other with no polygons in the way.

The pyvisgraph distance is 21595 nm and the calculated graph is shown below using folium. It is quite clear that the problem is that pyvisgraph does not take into account the crossing of the antimeridian. Is there any way this could be enabled or is there something I have overlooked?

[TaipanRex]

Yeah this is a problem I have tried to find a solution to. I tried looking at different types of map projections, but you run in to the same problem.

The only idea I have had (but not tried to implement) is to create two (or more) visibility graphs built on maps that are split not just on the 180th meridian, but for example on the 180th and a second on the 90th. I would think you would need to mess around with the coordinates to make that work (I'm not a coordinate system expert at all). You could then do the shortest path on both visibility graphs and choose the shortest one.

Not the most efficient solution, but I have yet to find another possible solution.

[TaipanRex]

Even the GeoJSON specification says they don't have a way to handle it:

3.1.9. Antimeridian Cutting

In representing Features that cross the antimeridian,
interoperability is improved by modifying their geometry. Any
geometry that crosses the antimeridian SHOULD be represented by
cutting it in two such that neither part's representation crosses the
antimeridian.

https://tools.ietf.org/html/rfc7946#section-3.1.9

[tjansson60]

Thanks for the quick reply. I definitely think a option (flag) for potentially slower, but antimeridian safe method would be great.

I have tried to look the code, but I cannot find any references to a specific reference system in it. Another possibility I could imagine would be to shift the coordinates of the polygons and points in latitude by 90 or 180 degrees and recalculate. If the yielded distance were smaller than the original, the distance and path (shifted back) was returned instead. It that way we would no need for a change of reference system?

That being said I definitely also not an expert on projections either, so I do not know if this is feasible? The only problem I could imagine is if the polygons after being shifted then crossed the antimeridian and then caused problems.

Did you imagine to generate two or several visibility graphs up front or in the generation do the calculation in one go?

[TaipanRex]

Did you imagine to generate two or several visibility graphs up front or in the generation do the calculation in one go?

I was thinking of generating two or more visibility graphs up front, run the shortest path on both and choose the shortest one.

Another possibility I could imagine would be to shift the coordinates of the polygons and points in latitude by 90 or 180 degrees

This is a good idea, you do one visibility graph v1 using the unmodified polygons. You then shift all the polygon longitudes by 90 and create a second visibility graph v2. That is done once, then every time you want to find shortest paths, you run shortest path on both v1 and v2 and choose the shortest. For voyage use case, this should handle most antimeridian issues. You could probably also put in some tests to avoid doing the second shortest path, like if the origin and destination "crossing meridian distance" is more than x nm.

[nhoussay]

I am completely new to these issues and this library has helped me a lot.

Should not we duplicate the world side by side in the graph? We can then check the shortest path for both the original coordinates and the coordinates shifted by 180 degrees. Does that make sense?

[TaipanRex]

@nhoussay yeah that is one way we have been thinking.

[robert-bressan]

@TaipanRex, instead of calculating the visibility graph in a projection, you might want to calculate the distances and visibility on the sphere. But, if a projection is desired, AuthaGraph may solve some issues, as it is seamless in any direction. Returning to sphere calculation:

In sphere geometry, two straight lines always cross in two antipodal points. Given two straights AB and CD, such none are antipodal, the cross product P = ±(A×B)×(C×D) returns the direction of those points. Normalize it to find the correct vector P.

To test where these points lie in the shortest or the longest segment, thus, the point is visible or not, a point P over the straight AB is within the shortest segment if (A · B) < (P · A) AND (A · B) < (P · B).
· = dot product.

Then, to apply Lee's algorithm, I would suggest the following steps:

First: rotate the sphere such the "north pole" is coincident with the point you are querying the visibility. A way to acheve this is place an ortonormal linear transform over cartesian coordinates.

Second: find the spherical coordinates (the usual ones: 0 to pi, 0 to 2 pi) where zenith = 0 is the point of interest (the auxiliary north pole). Sort the points by azimuth found (Lee's ray sweep), then, start evaluating the points using the zenith as a marker when to calculate the intersections.

[TaipanRex]

@robert-bressan thanks for your input, that is the first time I have seen a real solution to the problem. I will have to spend some time on trying to understand your math as this is a new area for me, but intuitively using sphere geometry seems sensible. I will also have to try and implement it in a way that users can choose between this and the current projection method, ref. #8. This is computationally more expensive, so a choice of method has to be given. Abstracting the distance function is fine, but as you state, might have to abstract part of Lee's algorithm as well.

[robert-bressan]

(writing from mobile. Sorry for mistakes).
You might consider two aspects:

• the visibility for the map is precalculated. So, the only points you should have this trouble are the starting and the finishing ones.

• you may want to test, before choosing the method, if the route crosses the "border".

• another method is to place copies of the map, putting and assign the vertex number with some bit markers (like v2 = v1+4096). Then in the adjacency list found, apply bit AND to find the real vertex (like v2 = v2 AND 4095). You might want to remove duplicates later.

[robert-bressan]

@TaipanRex , I have studied the problem a little more on the sphere geometry. Using quaternions to rotate might be helpful, but I don't know if it would be fast enough, because of all trig functions involved. In the plane, one may skip most of the trig functions because the sweep sort may be done without the atan2 function.
Steps in sphere geometry, assuming the input in (latitude, longitude).

1. Calculate the cartesian coordinates of all points Cost: 5 trig calculations, 2 multiplications
2. Pick a vertex, and find the quaternion which rotates it to north pole. Cost: 1 division, 1 sum, 2 multiplications, 2 trig calculations
3. Rotate all the points using quaternion rotation Cost: 24 multiplications, 16 sums (per vertex)
4. Convert the rotated points to spherical coordinates (r = 1, phi: 0 to pi, theta: 0 to 2pi) Cost: 2 trig calculations (per vertex)
5. Duplicate the points such a travel by the long edge be feasible Cost: 3 sums (per vertex)
6. Sort the points by the spherical coordinate \theta Cost: n log n
7. Perform Lee's analysis, using the coordinate \phi to measure the visibility. Cost: 4 comparisons in best case, 13 multiplications, 8 sums, 1 square root, 2 trigs and one comparison in worst case, per edge.

In the attachment, you may see some of the calculations:
PyVisGraph (1).pdf