Using PostGIS to calculate where lines intersect using Geography objects

Question

I'm experimenting with PostGIS Geography objects, but it seems like it's not doing the calculations on a sphere using great circles, but instead doing a projection

If I try to find the intersection of a line from LA to Paris (shown below), and another line that starts slightly north of LA ends up slightly south of Paris, I would expect the intersection to be somewhere in Canada

enter image description here

If I run this query, I get a very different result than what I expect:

select
 st_astext(st_intersection(
 ST_GeographyFromText('LINESTRING(-118.4079 33.9434, 2.5559 49.0083)')
 ,ST_GeographyFromText('LINESTRING(-118.4079 35.9434, 2.5559 45.0083)')));
 st_astext
 st_astext

POINT(-82.92390070238243 38.73039642066327)

It seems like ST_Intersection used straight lines (purple) for the calculation instead of the red line. How can I force PostGIS to work with geography here instead of geometry?

Answer 1

Try densifying your lines using ST_SEGMENTIZE:

select
 st_astext(st_intersection(
 st_segmentize(ST_GeographyFromText('LINESTRING(-118.4079 33.9434, 2.5559 49.0083)'),10)
 ,st_segmentize(ST_GeographyFromText('LINESTRING(-118.4079 35.9434, 2.5559 45.0083)'),10)));

gives:

POINT(-98.59590472401487 50.77296552878112)

• Q-bertsuit: Fantastic, thank you! What exactly does segmentize do?
• Ian Turton: follow the link to the manual page

Answer 2

One should be very cautious when using geography function over large distances, as some of them are just a wrapper over the geometry function, using the "best" possible projection, which is UTM for short distances and bad ones for longer distances (it degrades down to world mercator).

ST_Intersection is one them, and we can see its code contains:

SELECT @[email protected](
 @[email protected]_Transform(
 @[email protected]_Intersection(
 @[email protected]_Transform(@[email protected](1ドル), @extschema@._ST_BestSRID(1,ドル 2ドル)), 
 @[email protected]_Transform(@[email protected](2ドル), @extschema@._ST_BestSRID(1,ドル 2ドル)))
, 4326))

In fact, the intersection point when using geometries in 4326 is at POINT(-78.08663333333332 38.96503333333333), about 5 degrees away from the intersection point when using geographies.

One option, as shown in the other answer, is to add vertices (segmentize) the line, so that even in World Mercator the original line follows more or less the great circle arc. The shorter are the segments, the most accurate is the output. It works fairly well but you have to add, and transform, and intersect many vertices.

SELECT ST_NumPoints(st_segmentize(ST_GeographyFromText('LINESTRING(-118.4079 33.9434, 2.5559 49.0083)'),10)::geometry);
 st_numpoints
--------------
 1048577

Another option is to use the gnomonic projection, in which great circles are represented as straight lines. WARNING: this projection works in a single hemisphere. (EPSG ESRI 102034 in the northern hemisphere, EPSG ESRI 102036 in the southern hemisphere).

WITH src(geo1,geo2) as (
 SELECT ST_GeographyFromText('LINESTRING(-118.4079 33.9434, 2.5559 49.0083)'),
 ST_GeographyFromText('LINESTRING(-118.4079 35.9434, 2.5559 45.0083)'))
SELECT
 st_astext(
 st_transform(
 st_intersection(
 st_transform(geo1::geometry,102034),
 st_transform(geo2::geometry,102034)),4326))
FROM src;

Which leads to POINT(-98.59590472403889 50.77296552878083). The answer can be considered the same as when adding vertices, but the lines only had 2 vertices at all time, against 1 million vertices for each of the segmentized line. (ok, to be fair, we get a point 1cm away when using segment size of 1 000m, which still leads to 16 000 vertices)

• Q-bertsuit: Wow, thanks! Very interesting. Seems like there's a lot to be saved in terms of computation time by using the gnomonic projection
• Ian Turton: but as @JGH says only if you are sure your lines stay in one hemisphere which can't always be guaranteed