Geodetic coordinates: How to project a point on a line on the same meridian

Question

Given a geodetic segment, e.g., from Brussels to Moscow, and a geodetic point, e.g., Berlin, which can be expressed in PostGIS as follows

  select geography 'Linestring(4.35 50.85,37.617222 55.755833)', geography 'Point(13.405 52.52)';

I need to find the point on the segment that is exactly at the same meridian of the given point, i.e., their bearing is 0 degrees.

Any idea how to obtain this ?

For example, using PostGIS, I create a segment from Berlin to the North by taking Berlin's longitude with the maximum of the two latitudes of the initial segment and compute the intersection of the two segments as follows

select st_astext(st_intersection(geography 'Linestring(4.35 50.85,37.617222 55.755833)',
  geography 'Linestring(13.405 52.52,13.405 55.755833)'))
                  st_astext
----------------------------------------------
 POINT(13.407059592483968 53.163047143541235)

As can be seen the longitude is not exactly the same and the bearing would not be 0.

with test(inter) as (
  select st_intersection(geography 'Linestring(4.35 50.85,37.617222 55.755833)',
    geography 'Linestring(13.405 52.52,13.405 55.755833)') )
SELECT degrees(bearing(geography 'Point(13.405 52.52)', inter)) from test;
       degrees
---------------------
 0.11002408958628832

where bearing is computed using the traditional formula (e.g., https://gist.github.com/jeromer/2005586)

Answer 1

I find the postgis intersection result in you example surprising (and probably wrong). Doing the inverse of your query, it seems to me that there is a numerical precision issue:

select st_intersects(geography 'POINT(13.407059592484 53.1630471435412)',  geography 'Linestring(13.405 52.52,13.405 55.755833)')
--FALSE

select st_intersects(geography 'POINT(13.405 53.1630471435412)',  geography 'Linestring(13.405 52.52,13.405 55.755833)')
-- TRUE 

Answer 2

Berlin is located North of the line between Bruxelles and Moscow. So your two segments don't intersect each other. Instead of considering the line from Berlin to the highest latitude of the first segment you should directly consider a line going from the North pole to the South Pole, which gives you get the correct answer :

select st_astext(
  st_intersection(
    geography 'Linestring(4.35 50.85,37.617222 55.755833)',
    geography 'Linestring(13.405 89.999999,13.405 -89.999999)')) 
-- POINT(13.405 52.2417230551345)

Answer 3

This strangeness is caused by how intersections are computed when dealing with geographies.

The st_intersection doc says:

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.

Using your coordinates and peeking at the code, we can see the that a LAEA projection is used. If you were to extend the south-north line, a world mercator projection would be used - and the result would be better!

select _ST_BestSRID(geography 'Linestring(4.35 50.85,37.617222 55.755833)',
  geography 'Linestring(13.405 52.52,13.405 55.755833)');
   _st_bestsrid
--------------
       999247
(1 row)


select _ST_BestSRID(geography 'Linestring(4.35 50.85,37.617222 55.755833)',
  geography 'Linestring(13.405 42.52,13.405 65.755833)');
 _st_bestsrid
--------------
       999000
(1 row)


select st_astext(st_intersection(geography 'Linestring(4.35 50.85,37.617222 55.755833)',
  geography 'Linestring(13.405 42.52,13.405 65.755833)'));
           st_astext
--------------------------------
 POINT(13.405 52.2417230551345)
(1 row)  
  

But even though, re-projection over the entire lines length is occuring, which can affect the output. The solution is therefore to subdivide the inputs into many small segments, so the reprojected lines are more aligned with the original ones

select st_astext(st_intersection(  ST_Segmentize(geography 'Linestring(4.35 50.85,37.617222 55.755833)',1000),
  ST_Segmentize(geography 'Linestring(13.405 52.52,13.405 55.755833)',1000)));
                st_astext
------------------------------------------
 POINT(13.4050000115689 53.2633912803655)
(1 row)