Intersection algorithm that correctly handles 180° meridian and poles?

Question

Is there a well known intersection algorithm that correctly handles 180 meridian and poles?

For example, let's say we have a list of latitude and longitude values that represent Antarctica. Let's say we also have a simple polygon that represents an airplane. We want to know if the airplane is over Antarctica.

Generic 2D polygon intersection does not work in this case if you simply use latitude for y and longitude for x, because the flat coordinate system will have edges at -180 and 180 longitude, and -90 and 90 latitude. Antarctica will go off the page on three of the edges.

Answer 1

GIS (as a field) hasn't done too hot when it comes to really grappling with the surface of the globe.

For example, your problem isn't fully defined. Unlike in 2D, where we know the edges of a polygon are composed of straight lines, what are they on the globe? Arcs of great circles, minimizing distance between vertices, are a good choice but not the only one. Straight lines (that thus travel below the surface of the globe) are another choice, for example. (Rhumb lines are a choice as well but potentially a dumb one, especially near the poles...) Here's a great link WRT edge interpretation issues: http://blog.opengeo.org/2010/08/10/shape-of-a-polygon/

Another big issue with the project-to-2d-then-intersect approach, besides the singularities others have mentioned, is that any newly created intersection points (where polygon edges cross) will be out-of-place and dependent on the projection used. (I believe this is where the recommendation for densification comes from: By adding a ton of intermediate vertices, you get a reduced error from projection of the middle of polygon edges.)

Assuming you don't want to make all the compromises and workarounds of projection to 2D and are looking to think through and code something up yourself, I've done a little bit of prototype (and prototype-only!) code for a client dealing with this.

Here's a sketch of the approach. You'll need to know what a vector is as well as the meaning of the dot and cross products. (Caveat: wikipedia is convenient for quick links, but sufficient if it's your first introduction to the topics. A good 3D graphics tutorial will help.)

• Represent a point on the sphere by a unit Cartesian 3D vector. The point on the surface of the earth is where, if you extended the vector into an infinitely-long ray, it would intersect the surface of the earth.
• Represent great circles by a plane through the origin. (In 3D a single unit vector suffices to define a plane through the origin; it is the normal to the plane.) The great circle is the intersection of the entire plane with the surface of the earth.
• You can find the intersection points of two great circles by intersecting their two planes.
• Define an arc of a great circle by two points. The vector that defines the great circle normal is the cross product of the start and end points' vectors.
• To determine if a point that we know is on the great circle lies within the arc, create two planes such that: they are perpendicular to the plane of the great circle, one contains the start point and the other contains the end point, and they are oriented facing each other. Then the point is on the arc iff it lies to the 'inside' of both planes. (To help visualize: you have created pac-man's smile as he chomps down upon the arc. If the test point is between his jaws, then it lies on the arc, as we already know it is upon the great circle.)
• To determine if two arcs cross: find the two points of intersection of their corresponding great circles, then test each point to see if it lies within both arcs.
• A less-ambiguous definition: a polygon is a collection of points, each of which are connected by edges which consist of arcs of great circles. The points are ordered such that if you walk along the surface of the earth along the polygon's edges, the 'inside' of the polygon will lie to your left. Let's leave complex polygons (islands, holes, and self-intersections) out of it for now.
• You can tell whether a point is on the right or left side of a plane via the sign of the dot product of their corresponding vectors. (This is equivalent to, as you travel around the great circle, whether the point on the surface of the sphere lies to your left or your right.)
• An accurate test to determine if a point is inside a polygon: Does it lie on the left side of all the edges?
• Now we have the ability to test if a point is inside a polygon and to determine edge crossings: the ingredients for polygon-polygon intersections! The margin of this comment is too small to write out a full algorithm but the basic steps are to: (a) find all intersections, then (b) walk edges, alternating which polygon you're walking on as you encounter intersection points.
• Once all of the above is working, start thinking about indexing strategies to make it faster, since the point-in-polygon I outlines is O(n) in the number of edges, and polygon intersection O(m * n) in number of edges.

Pheew.

There are some great advantages to this approach: All of the operations above boil down to only multiplications and additions. (After converting data to this representation: e.g. lat/long coordinated require some trig to get a Cartesian XYZ vector.) There are no singularities at the poles or in the coordinate system, and not too many special cases to worry about (example special case: parallel overlapping edges).

Peeking at the code it looks like the Spheres package someone else linked follows some of this approach, although it also looks a bit half baked.

PostGIS may also use a similar approach over its geography datatype, but I haven't looked under the hood. I do know that for spatial indexing, at least, they use an R-tree over 3D cartesian.

(Note: this answer became long enough that I'm probably going to edit up into a blog post... Feedback/comments very welcome!)

Answer 2

I'm actually looking for a general solution that will work for many large polygons.

Most of the answers I have read so far are naturally focusing on finding a suitable projection for your features. This may be counter-intuitive, but here you need to consider where a projection fails to work, not where it works well.

Some projections fail to work only at one point: these are among the so-called "global" projections. A good example is the azimuthal equidistant: the only point you cannot project is the one diametrically opposite to its point of origin. (Actually, it's not strictly true that you can't project that last point; I'm glossing over a technical issue. Let it suffice that the projection experiences severe problems at the opposite pole.)

Thus, if you can find a single point P outside the union of all the polygons you plan to work with, then all you need do is work within a global projection that fails only at this single point. (E.g., use an oblique equidistant azimuthal centered at the point diametrically opposite P.)

(You can find such a point automatically. For example, create a polygon covering the earth, slightly buffer the union of all your features, remove this buffer from the earth polygon, and select any arbitrary point within what is left. The reason for buffering is that it's best to stay away from the singular points where a projection fails, the further the better: this is an issue of computational accuracy.)

There are plenty of global projections with desirable properties, including equidistant, conformal (e.g., stereographic), and equal-area. Thus, choosing a global projection for your analysis is not a severe limitation.

The main limitation is that the most popular software out there, ArcGIS, does not support many oblique versions of global projections! :-(

Answer 3

You can solve this with projections!

Just because the "ends" of the scales are at -180 and +180, doesn't mean you have to consider the map with those ends.

Use EPSG:32761 ( http://spatialreference.org/ref/epsg/32761/ )

Convert your latitude and longitude values into coordinates in the UPS space, and then you can use regular geometry calculations to determine if your airplane is over Antarctica.

More information about the Universal Polar Stereographic coordinate system is here: http://en.wikipedia.org/wiki/Universal_Polar_Stereographic_coordinate_system

If you're just trying to determine an intersection algorithm for gigantic polygons, then treat the earth as a spheroid (or even a sphere) and use spherical geometry to analyze your polygons.

I found a library called Spheres that's GPLed, so you could review its code to see precisely what algorithms they use for arc intersecion and "points within a spherical polygon".

However, when you want to DRAW these polygons on a flat surface (like paper, or your screen) then you will need to pick a projection, or at least a projection style.

Answer 4

One solution to this problem is to use some sort of 3-dimensional coordinate system for the earth instead of a two dimensional coordinate system. If the 2D polygons are given a large virtual amount of depth, such as half way to the center of the earth, then the polygons don't need to be modeled as curves on the surface of the earth. It should be enough for only the verticies to be at the surface.

Answer 5

While there are issues with using Cartesian coordinates, I don't see the rub in the example you've provided. Here's what it looks like to me, with a polygon covering Antarctica bounded by the edges of the coordinate space:

_{(source: geohack.net)}

The intersection between the plane and the continent can safely be computed without performing any complex transformations.

Answer 6

A possible solution is to cut the polygon into sub-polygons at all the edges (180 meridian and poles). The software could even maintain some sort of reference between the original polygon and the sub-polygons.

Answer 7

Take the centroid of the airplane and use it as the point of tangency to for an azimuthal projection. You'd also need to densify the geometries before projecting them for cases like the area of the earth currently seeing daylight.