Intersecting great circle arcs and polygons

Question

I am trying to find the point of intersection between a great cirle arc and a polygon. I am working in long/lat, where I am given the endpoints of the great circle arc, A and B, and the vertices, C, D, E and F of the polygon in long/lat. The point A is inside the polygon and point B is outside of the polygon. It is obvious by construction that the great circle arc intersects the polygon on the CD side of the polygon. When I compute the intersection of the great circle arc AB, represented by a LINESTRING, with the POLYGON CDEF, I get a LINESTRING from a point on CD to A; when I compute the intersection of the great circle arc AB with the LINESTRING CD, I get a different point on CD. Here are the details:

A:  -122.3 47.4 
B:  -84.4 33.6

C: -107 49
D: -107 31
E: -130 31
F: -130 49

(1) intersect LINESTRING and POLYNOMIAL

select st_astext(st_intersection(st_geographyfromtext('linestring(-84.4 33.6, -122.3 47.4)'), st_geographyfromtext('polygon((-107 49, -107 31, -130 31, 130 49, -107 49))')));

====>
LINESTRING(-107 42.1770201452269, -122.3 47.3999999999286

(2) intersect LINESTRING and LINESTRING

select st_astext(st_intersection(st_geographyfromtext('linestring(-84.4 33.6, -122.3 47.4)'), st_geographyfromtext('LINESTRING(-107 49, -107 31)')));

====>
POINT(-107.040211660479 43.9998296744469)

As you can see, this point is not the same as the endpoing in the output of (1).

What is going on?

I am running PostgreSQL9.3.4 on x86_64-unknown-linux-gnu, compiled by gcc (Ubuntu 4.8.2-16ubuntu6) 4.8.2, 64-bit POSTGIS="2.1.2 r12389" GEOS="3.4.2-CAPI-1.8.2 r3921" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.10.1, released 2013/08/26" LIBXML="2.9.1" LIBJSON="UNKNOWN" RASTER

Finally, a colleague of mine computed the intersection of the two great cirlce arcs directly calling the Java version of GDAL/OGR and got -107 41.8290237496.

Answer 1

The documentation page notes that for ST_Intersection(geography, geography),

For geography this is really a thin wrapper around the geometry implementation. It first determines the best SRID that fits the bounding box of the 2 geography objects (if geography objects are within one half zone UTM but not same UTM will pick one of those) (favoring UTM or Lambert Azimuthal Equal Area (LAEA) north/south pole, and falling back on mercator in worst case scenario) and then intersection in that best fit planar spatial ref and retransforms back to WGS84 geography.

The reason you're getting different results between your line and polygon examples is probably because the line and polygon have different extents, so the heuristic is choosing different planar projections in the two cases in which to perform the calculation.

In any case, the heuristic won't choose a planar projection that will give you an exact result, and that might be what you want. Since you're doing constructive geometry, you might want to build your own heuristic, that chooses a good Gnomic projection to use for your input pairs of geometry.

The Gnomic is terrible for measuring distances or areas, but it's excellent for determining intersection points, since it turns great circles into straight lines. So what you want to do is

• cast your geography to geometry,
• reproject to a "good" gnomic (one with a centre point near the centers of your inputs),
• run the ST_Intersection there in geometry,
• reproject back to geographics,
• cast back to geography, and there you go.