I’ve been working on a library for NASA that validates polygons describing areas on the Earth. The polygons are metadata for Earth Science data, such as images from an orbiting satellite. A polygon is represented by a list of points in counter clockwise order defining an outer boundary. The polygon can also have holes which are represented by additional lists of points. One component of polygon validation is checking that all the holes are within the outer boundary.
Valid Polygon
Invalid Polygon
Imperfect Spheres
The Earth is not a perfect sphere. It’s squashed on the poles which makes it closer to a shape called an oblate spheroid. Using a spherical representation of the Earth can result in up to 0.5{f915c993dcd57782d8bbdb9f9291c5d4b3b1a7c844f782feb26c4f553a22f754} errors when calculating distances. I’m only checking polygon validity so this amount of potential error is acceptable. Validating polygons that can appear anywhere on the Earth isn’t trivial. They can cover a pole or cross the antimeridian (The vertical line at 180/-180 degrees longitude). They can be small or very large, concave or convex. While it’s easy for a person to look at a polygon on a map and see that it’s valid it’s much more difficult for a computer. The computer only gets a list of numbers. The program must make sense of the numbers and apply a set of rules to ensure it’s valid.
Checking If a Point Is In a Polygon
Luckily for us this is well worn territory. There are many great resources for spherical math online. To make sure that our “holes” are inside the polygon’s outer boundary we have to check if each point of the hole is inside the outer ring. GeoSpatialMethods.org lists one way for determing if a point is inside a polygon. You first have to know a point external to the polygon. You then count how many times a line from your point to the external point crosses the edges of the polygon. If it crosses an odd number of times it’s in the polygon. If it’s 0 or an even number it’s outside the polygon.
Finding an External Point
Finding a point outside of the polygon is a little tricky. Many resources online assume a human will provide a point external to the polygon. The Earth Science data I’m dealing with doesn’t define an external point though. GeoSpatialMethods.org lists a few ways of guessing the external point like checking the longitude and latitude ranges. That site doesn’t provide the full solution admitting “this is clearly just a guess”. I wasn’t able to find a definitive answer in the Google searching I did of the right way to find an external point to a spherical polygon.
I’ve come up with a method to definitively determine what is inside or outside a spherical polygon. Spherical math is older than Christopher Columbus so I’m sure I’m not the first person to discover this. I couldn’t find it anywhere online so I thought I should at least document it for others. My method involves summing the differences of the courses as you follow the points around the polygon. Let’s look at an example.
This is a picture of a simple polygon with the points 0,1, 4,1, 6,5, 2,5, 0,1 (lon1, lat1, lon2, lat2…). Imagine you’re in an airplane and you’re going to fly around the polygon in the order of the points exactly along the great circle arcs. When flying away from a point you look at your compass heading and note the number. When arriving at a point you note the current compass heading. You continue to note compass headings at all the subsequent points as you’re arriving and leaving. I’ve put the compass headings above in green. There’s a slight difference in the compass headings when leaving a point versus arriving at a point due to the fact the way great circle arcs curve. (I’ll explain more on this later). If you take the difference of each course from the one before it and then sum all the differences you’ll see it adds up to 360. Courses can be calculated using the formula here on the Aviation Formulary. The table below shows all of the compass headings with deltas and sums.
Course Delta Table for Counter Clockwise Polygon
The 360 degrees indicates that we’ve made a complete turn in the counter clockwise direction. If we reverse the order of the points then the course deltas sum up to -360. The -360 degrees indicates a full turn in the clockwise direction.
Course Delta Table for Clockwise Polygon
Let’s see what happens with a polygon around the North Pole. Below is an image showing a polygon around the North Pole with the points -135,85, -45,85, 45,85, 135,85, -135,85. The points are spaced evenly around the pole so that the initial and ending courses match on each arc. There’s about 90 degrees of difference between each initial and ending course. The course or compass heading tells which direction you’re facing in relation to the North Pole. When you start out flying from point 1 you could see the North Pole off to your left. By the time you’ve gotten to point 2 the North Pole is now almost behind you. Even though you’ve been flying in a straight line your compass heading was changing the whole time. It does this anywhere on the Earth except along the equator. The effect is much more noticeable near the poles.
Polygon Around North Pole
Here’s the table of angle deltas. You’ll notice that the sum in now 0.
Course Delta Table for Polygon Around North Pole
The sum is 0 because the course delta sum is really measuring how many times the North Pole has looped around you. When flying around the pole the course changes constantly but the North Pole always stays generally on the left side (or the right side if flying in a clockwise direction around the pole).
Bringing It Together
We’ve seen that a normal counter-clockwise polygon’s course delta sum is 360, a clockwise polygon’s sum is -360, and a polygon around the North Pole sums to 0. It will also be 0 if the polygon is around the South Pole. We can use this information to determine exactly what a polygon contains. If the course delta sum is 360 we know that the polygon does not contain either pole. If the course delta sum is -360 it basically covers everything “outside” of itself. It covers both poles in other words. If the course delta sum is 0 we know that it contains either the North Pole or the South Pole but not both. We can use the same course delta sum approach with the longitude values of the points of the polygon to determine if it wraps counter clockwise around the North or South Poles.
This approach tells us, in a definitive way, what the polygon’s relationship to the poles are. It works for convex and concave polygons with any number of points. I’ve tried it with polygons that spiral around a pole but don’t actually cover it and it’s worked perfectly every time. Once you know the polygon’s relationship to the poles determining an external point is easier. If the polygon covers the North Pole then South Pole is a convenient external point. If the polygon doesn’t cover either pole then both the North and South Poles can be used as external points.
It’s more difficult if the polygon covers both poles. I didn’t have to support this for my project but it should be possible given what you know. Since the polygon covers both poles you know the reverse of the polygon (points in opposite order) cover neither pole. It should be easier to find a point inside the reversed polygon though I’m unsure of a way to accomplish this beyond simply making educated guesses and testing for point containment. It’s also possible to cut the polygon that crosses both poles in half across the equator and deal with the two separately. More work could be done here to find external points for polygons that cover both poles. I probably won’t continue to look for a solution here since it’s not required by our customer. For all other cases this approach has worked perfectly.
Footnote
I didn’t want to get into this in the main article but it’s important to know that a single external point isn’t sufficient. Most spherical algorithms fail when dealing with antipodal points (points that are on complete opposite sides of the earth). You need two external points in case the point being tested is antipodal to the first external point you picked.
Update: I’ve posted a follow up here.
Update: Try out the tool I used to create the screen shots here.