Point in polygon on Earth globe

Question

I have a list of coordinates (latitude, longitude) that define a polygon. Its edges are created by connecting two points with the arc that is the shortest path between those points.

My problem is to determine whether another point (let's call it U) lays in or out of the polygon. I've been searching web for hours looking for an algorithm that will be complete and won't have any flaws. Here's what I want my algorithm to support and what to accept (in terms of possible weaknesses):

1. The Earth may be treated as a perfect sphere (from what I've read it results in 0.3% precision loss that I'm fine with).
2. It must correctly handle polygons that cross International Date Line.
3. It must correctly handle polygons that span over the North Pole and South Pole.

I've decided to implement the following approach (as a modification of ray casting algorithm that works for 2D scenario).

1. I want to pick the point S (latitude, longitude) that is outside of the polygon.
2. For each pair of vertices that define a single edge, I want to calculate the great circle (let's call it G).
3. I want to calculate the great circle for pair of points S and U.
4. For each great circle defined in point 2, I want to calculate whether this great circle intersects with G. If so, I'll check if the intersection point lays on the edge of the polygon.
5. I will count how many intersections there are, and based on that (even/odd) I'll decide if point U is inside/outside of the polygon.

I know how to implement the calculations from points 2 to 5, but I don't have a clue how to pick a starting point S. It's not that obvious as on 2D plane, since I can't just pick a point that is to the left of the leftmost point.

Any ideas on how can I pick this point (S) and if my approach makes sense and is optimal?

Thanks for any input!

Answer 1

If your polygons are local, you can just take the plane tangent to the earth sphere at the point B, and then calculate the projection of the polygon vertices on that plane, so that the problem becomes reduced to a 2D one.

This method introduces a small error as you are approximating the spherical arcs with straight lines in the projection. If your polygons are small it would probably be insignificant, otherwise, you can add intermediate points along the arcs when doing the projection.

You should also take into account the polygons on the antipodes of B, but those could be discarded taking into account the polygons orientation, or checking the distance between B and some polygon vertex.

Finally, if you have to query too many points for that, you may like to pick some fixed projection planes (for instance, those forming an octahedron wrapping the sphere) and precalculate the projection of the polygons on then. You could even create some 2d indexing structure as a quadtree for every one in order to speed up the lookup.