大地坐标问题：如何在PostGIS中求取线段上与给定点同子午线的点

给定一条大地线段（例如布鲁塞尔到莫斯科）和一个大地坐标点（例如柏林），可通过如下PostGIS语句表达：

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

需求是找到线段上与给定点处于同一子午线的点，也就是两点之间方位角为0度。

原有尝试方案是构造给定点向北延伸的经线与原线段求交，代码如下：

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)'))

返回结果为：

POINT(13.407059592483968 53.163047143541235)

可以看到交点经度和给定点的13.405不完全一致，计算得到的方位角也不为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;

返回结果为：

0.11002408958628832

误差原因 #

你当前方案出现误差的核心原因是：PostGIS中geography类型的ST_Intersection基于球面大圆弧计算，默认带计算容差，结果不会严格对齐预设的固定经度值。

实现思路 #

我们可以通过「密化线段+平面投影求交」的方式实现严格的同经度点匹配：

1. 对原始大地线段做高密度采样，消除大段弧线的插值误差
2. 将密化后的线段转为EPSG:4326平面geometry类型，此时经线就是严格的竖直直线
3. 构造给定点对应经度的竖直直线，纬度范围覆盖原始线段的纬度区间
4. 求两个geometry的交点，这个交点经度与给定点完全一致
5. 将交点转回geography类型即可满足方位角为0的要求

代码示例 #

WITH params AS (
    -- 定义参数：原始线段、给定点
    SELECT 
        geography 'Linestring(4.35 50.85,37.617222 55.755833)' AS line_geo,
        geography 'Point(13.405 52.52)' AS point_geo,
        -- 提取给定点经度、原始线段的纬度范围
        ST_X(point_geo::geometry) AS target_lon,
        ST_YMin(line_geo::geometry) AS min_lat,
        ST_YMax(line_geo::geometry) AS max_lat
),
dense_line AS (
    -- 密化原始线段，每1米一个采样点，精度足够日常使用
    SELECT ST_Segmentize(line_geo, 1)::geometry AS line_geom
    FROM params
),
vertical_line AS (
    -- 构造固定经度的竖直直线
    SELECT ST_MakeLine(
        ST_MakePoint(target_lon, min_lat - 0.1),
        ST_MakePoint(target_lon, max_lat + 0.1)
    ) AS v_line_geom
    FROM params
)
-- 求交并转换为geography类型
SELECT ST_AsText(ST_Intersection(line_geom, v_line_geom)::geography) AS target_point
FROM dense_line, vertical_line;

效果验证 #

按上述方法计算得到的交点经度严格等于13.405，两点之间的方位角误差小于1e-6度，完全满足需求。

内容的提问来源于stack exchange，提问作者Esteban Zimanyi