# Projection for polygon intersection confined into one hemisphere?

I have polygons where each vertex is a (lon,lat) geographic coordinate. I want to find out if they intersect or contain each other. The polygons are guaranteed to never exceed a hemisphere (so they'll never 'wrap around' beyond a plane that cuts the Earth in half), but they can vary dramatically in scale.

I think one way this might be done is by transforming the polygons using some specific projection, doing the math for intersection in cartesian coordinates and then 'unprojecting' back to geographic coordinates.

Is there a projection I can generate that is fast to apply, even if its at the cost of somewhat not-so-accurate (~50m-100m) results?

I'm not using any database tools or large libs here, just C++ as I need a small, portable set of code. It would be handy if I could come up with the projection on the fly based on the geographic coordinate polygons rather than having to refer to a list of 'proper' SRSs.

• I'd think that what you're looking for is a library capable of doing geometry relationship calculation, not a projection (which you would need to implement, and which wouldn't handle the intersection operation). Jul 10, 2014 at 18:01
• I already have a lib that does that in cartesian space.
– Pris
Jul 10, 2014 at 18:09
• So long as your feature vertex density is appropriate, there shouldn't be much difference between Cartesian comparison in projected and geographic coordinate systems. Jul 10, 2014 at 18:35
• I'm assuming you mean simply map lon and lat to 'x' and 'y' when you say use a geographic coordinate system. How do I determine the appropriate vertex density for a given accuracy? My features are typically large fences or bounds and have only a few vertices, so I'd need to generate additional ones.
– Pris
Jul 10, 2014 at 18:46
• How far apart are these vertices? You wouldn't need great circle densification short of a few nautical miles. In fact, if necessary, great circle densification is your solution, rather than a custom projection with each operation. Jul 10, 2014 at 19:18

## 2 Answers

Yes, as whuber says, the gnomonic projection is the way to go. GeographicLib (written by me) includes an ellipsoidal generalization of this projection in which geodesics are very nearly staight. For example, if you consider two points 4000 km apart, the straight line connecting them in the gnomonic projection deviates from the geodesic by less than 30 m. You can access this projection from C++, Matlab, or via a command line utility.

• Your library looks really helpful and seems like a plug and play solution for my problem.
– Pris
Jul 11, 2014 at 19:43

Do the computations using the Cartesian formulas in a Gnomonic projection. These projections cover a single hemisphere and map all portions of great circles to line segments. Presumably, in the representation of your spherical polygons as sequences of (lon,lat) values it is implicit that each successive pair of vertices is to be traveled along a geodesic. These are great circles in a spherical earth model and close to great circles in an ellipsoidal model. The Cartesian formulas assume that the mapped images of these vertices will be connected by straight line segments, which are exactly the images of the great circle arcs implied in the original representation. Thus essentially no error is made by the projection: relationships of containment and intersection will be preserved.

The formulas for the polar Gnomonic projections on the sphere are relatively simple. (This is one of the earliest projections invented.) For oblique aspects (needed if the polygons cross the Equator) consider rotating all the (lon,lat) coordinates to place the vertex of the projection at one of the poles and then apply the Gnomonic formulas.