大地坐标：如何在同一子午线上的线上投影点

Question

给定一个大地线段，例如从布鲁塞尔到莫斯科，以及一个大地线点，例如柏林，在 PostGIS 中可以如下表示

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

我需要找到线段上与给定点完全相同的子午线的点，即它们的方位角为 0 度。

知道如何获得这个吗？

例如，使用 PostGIS，我通过将柏林的经度与初始段的两个纬度中的最大值取值来创建从柏林到北方的段，并按如下方式计算这两个段的交点

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)

可以看出，经度不完全相同，方位也不为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

使用传统公式计算方位（例如，https://gist.github.com/jeromer/2005586）

Answer 1

我发现您的示例中的 postgis 交集结果令人惊讶（并且可能是错误的）。与您的查询相反，在我看来，存在数值精度问题：

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

这种奇怪是由在处理地理时如何计算交叉点造成的。

st_intersection 文档说：

地理：对于地理来说，这实际上是几何实现的薄包装。它首先确定适合 2 个地理对象的边界框的最佳 SRID（如果地理对象在 UTM 的一半区域内，但不同的 UTM 将选择其中一个）（偏爱 UTM 或 Lambert Azimuthal Equal Area (LAEA) 北/南极，并在最坏的情况下退回到墨卡托），然后在最适合的平面空间参考中相交并重新转换回 WGS84 地理。

使用您的坐标并查看代码，我们可以看到使用了 LAEA 投影。如果你要延长南北线，将使用世界墨卡托投影 - 结果会更好！

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)  
  

但即便如此，整个行长度上的重新投影仍在发生，这可能会影响输出。因此解决方案是将输入细分为许多小段，因此重新投影的线与原始线更对齐

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)

Answer 3

柏林位于布鲁塞尔和莫斯科之间的分界线以北。所以你的两个部分不会相互交叉。与其考虑从柏林到第一段最高纬度的线，不如直接考虑从北极到南极的线，这样就可以得到正确答案：

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)