PostGIS航线分析

postgis-users 邮件列表中的一位社区成员最近提出了一个问题：

我有一张高程点表，我想计算出穿过这些点的航线的高程剖面图。如何？

这个问题很好地展示了我最喜欢的一些空间工具，包括索引、球体上的点对点距离查询和最近邻查询。原问题的作者很好地分享了他的高程数据，所以我可以给你看一些图片和一种方法。

1、每2KM高程点

数据集为欧洲1.5M高程点，每2km一个点。当绘制在地图上时，你可以看到它非常密集。

表结构很简单：

 Column   |       Type       | Collation | Nullable | Default
-----------+------------------+-----------+----------+---------
 latlng    | geography        |           |          |
 lat       | double precision |           |          |
 lng       | double precision |           |          |
 elevation | real             |           |          |

我们所需要的只是 latlng 列上的空间索引，以使我们的最近邻搜索更快。

CREATE INDEX elevation_latlng_x ON elevation USING GIST (latlng)

2、航线

航线不同于陆上航线，因为它们不仅“像乌鸦飞过一样”，而且“像乌鸦飞过一个球体”。在更长的距离和更高的纬度上，效果更明显，但大圆路线与点对点路线并不完全相同。

例如，伦敦和华沙之间的直线可以这样构造：

SELECT ST_MakeLine(ST_MakePoint(-0.061,51.482),  -- London
                   ST_MakePoint(21.252,52.150))  -- Warsaw

但是直线（平面上的最短路径）不是我们想要的，我们想要一个大圆（球体上的最短路径）。因此，我们将线转换为地理类型，它使用大圆在顶点之间进行插值，然后对线进行 ST_Segmentize，以在端点之间添加新顶点。

SELECT  ST_Segmentize(geography(
          ST_MakeLine(
            ST_MakePoint(-0.061,51.482),
            ST_MakePoint(21.252,52.150))
        ), 4000)

红线是大圆，蓝线是原来的笛卡尔线。

3、最近点

对于最终查询，我们将构建马德里和瓜达拉哈拉之间的飞行路径，这两个城市距离足够近，我们可以看到实际的高程点。我们将沿着该路径每 4KM 放置一个顶点，并找到离每个顶点最近的高程点。

要找到最近点，我们使用 <-> “索引有序距离”运算符。当放在查询的 ORDER BY 子句中时 <-> 不仅计算左右操作数之间的距离，它还返回从可能的最小距离开始的距离。

因此，查询离马德里最近的高程点 (-3.6996,40.4009) 是：

SELECT
  e.elevation,
  ST_AsText(e.latlng, 3) AS coords,
  geography(ST_MakePoint(-3.6996,40.4009)) <-> e.latlng AS dist
FROM elevation e
ORDER BY dist
LIMIT 1

返回结果如下（在 2ms 内，因为有索引辅助）：

elevation |        coords        |        dist
-----------+----------------------+--------------------
  650.8142 | POINT(-3.707 40.409) | 1045.4185642910293

4、最近的飞行路线

索引辅助最近邻运算符 <-> 是一个奇怪的运算符，因为它需要运算符一侧的文字值。它不能接受两个变量参数，那么我们如何将具有多个点的飞行路径连接到高程网格？

答案是使用 LATERAL JOIN 来控制迭代并将飞行路径点（一次一个）提供给最近的邻居计算。

以块为单位，查询首先构建飞行路径。

WITH path AS (
    SELECT ST_Segmentize(geography(
        ST_MakeLine(
            ST_MakePoint(-3.6996,40.4009),  -- Madrid
            ST_MakePoint(-3.1675,40.6582))  -- Guadalajara
        ), 4000) AS geog
)

然后它使用 ST_DumpPoints 将路径分解为每个顶点的一条记录：

points AS (
    SELECT (ST_DumpPoints(geog::geometry)).geom FROM path
)

最后，它对高程点进行横向连接，为每个路径点找到最近的高程。

SELECT
  d.elevation,
  round(d.dist) AS distance,
  ST_AsText(d.latlng, 3) AS elevation_pt,
  ST_AsText(points.geom, 3) AS path_pt
FROM points
CROSS JOIN LATERAL (
  SELECT
    e.elevation, e.latlng,
    points.geom::geography <-> e.latlng AS dist
  FROM elevation e
  ORDER BY dist
  LIMIT 1
) d

最终结果大约需要 14ms 运行，并找到飞行路线上每个点的高度。

 elevation | distance |     elevation_pt     |       path_pt
-----------+----------+----------------------+----------------------
  650.8142 |     1045 | POINT(-3.707 40.409) | POINT(-3.7 40.401)
  630.8268 |      563 | POINT(-3.66 40.415)  | POINT(-3.666 40.417)
  701.6845 |      728 | POINT(-3.641 40.436) | POINT(-3.633 40.433)
 632.26556 |      882 | POINT(-3.595 40.443) | POINT(-3.6 40.449)
   611.845 |      733 | POINT(-3.575 40.464) | POINT(-3.567 40.465)
  607.7102 |      701 | POINT(-3.533 40.488) | POINT(-3.534 40.482)
  644.5825 |     1057 | POINT(-3.509 40.491) | POINT(-3.501 40.498)
 637.69604 |      166 | POINT(-3.467 40.515) | POINT(-3.467 40.514)
  666.1685 |     1335 | POINT(-3.447 40.536) | POINT(-3.434 40.53)
   673.436 |      373 | POINT(-3.401 40.543) | POINT(-3.401 40.546)
  717.8294 |      938 | POINT(-3.358 40.567) | POINT(-3.368 40.562)
 696.11383 |      910 | POINT(-3.335 40.57)  | POINT(-3.334 40.578)
 709.38855 |      723 | POINT(-3.292 40.594) | POINT(-3.301 40.594)
  716.1025 |      697 | POINT(-3.273 40.615) | POINT(-3.268 40.61)
 655.39514 |      862 | POINT(-3.226 40.621) | POINT(-3.234 40.626)
 659.81757 |      501 | POINT(-3.207 40.642) | POINT(-3.201 40.642)
   654.524 |      943 | POINT(-3.164 40.666) | POINT(-3.168 40.658)

5、结束语

总结一下——我推荐的计算高程剖面和飞行路线的方法是：

创建空间索引以加快最近邻搜索
使用地理类型在线之间查询，因此最短路径在球体上而不是直线上
使用索引有序距离运算符查找两点之间的最短距离
使用最近邻运算符与高程点的横向连接。

原文链接：Elevation Profiles and Flightlines with PostGIS

BimAnt翻译整理，转载请标明出处