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!

Saxophone

Each simple polygon dissects the sphere into two parts. Which one is 'inside' and which is 'outside'...? It's quite simple on a plane – the outside is the unbound part, but on a sphere both parts are bounded. And they can even be equal! Consider a 'quadrilateral polygon' whose all vertices are on the equator, with coordinates in degrees: (0, 0), (0, 90), (0,180), (0, -90)... – Wideopen

To my surprise it seems I answered (a near-duplicate of) this question a few years ago #3067595 – Neurasthenia

In my case I will have areas that shouldn't span over half of the world (it will be more or less local regions rather than global ones), so we may assume that 'inside' is smaller than 'outside'. Even without that assumption I don't think that it matters, if we can solve the problem for 'inside' we can do it for the 'outside' as well. – Saxophone

@HighPerformanceMark - in that solution you suggest to pick a point that is pi/2 away from the target, but how can I be sure that this point is outside of the polygon? – Saxophone

Have a look at "Locating a Point on a Spherical Surface Relative to a Spherical Polygon of Arbitrary Shape". – Jeannajeanne

To be honest, I've completely forgotten how that works ! Obviously (ha ha) it must work somehow or the answer wouldn't have been accepted ... – Neurasthenia

@Jean-BaptisteYunès - thanks for the article, it clears all my doubts I had and I was able to migrate the solution to Java! – Saxophone

@Saxophone can you please show your solution in the answer and accept it, for strangers like me? – Recur

@Recur - I was doing it commercially for my employee, so I don't possess rights to publish it - sorry :( but I can say that it was pretty straightforward to map it one-to-one to Java code – Saxophone

s2geometry.io might fit – Deflocculate

Does this answer your question? Is a point inside or outside a polygon which is on the surface of a globe – Semivowel

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.

Batfish

You are correct in pointing out that the antipodes are a problem with this approach - a polygon that contains the antipode but not the point itself would look like it contains the point. Checking the winding of the projected polygon would help in some cases, but would not be foolproof. There are also some cases where a spherical polygon would project to something that is not a polygon in 2d - consider a polygon that is shaped like the stiches of a baseball. One piece of the baseball's leather is "inside" the other piece is "outside". But when projected to 2d, it's just a "C" shape. – Discontinue

Answer 2

The biggest issue is to define what we mean by 'inside the polygon'.

On a sphere, every polygon (as long as the lines are not intersecting) defines two regions of the sphere. Both regions are equally qualified to be called the inside of the polygon.

Consider a simple, 1-meter on a side, yellow square around the south pole.

You can think of the yellow area to be the inside of the square OR you can think of the square enclosing everything north of each line (the rest of the earth).

So, technically, any point on the sphere 'validly' inside the polygon.

The only way to disambiguate is to select which side of the polygon you want. For example, define the interior to always be the area to the right of each edge.

Jovita