# Using PostGIS to calculate where lines intersect using Geography objects

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

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?

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)
``````
• Fantastic, thank you! What exactly does segmentize do? Commented Feb 14, 2023 at 9:52
• follow the link to the manual page Commented Feb 14, 2023 at 10:00

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)

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