PostGIS时空轨迹实战

NSDT工具推荐： Three.js AI纹理开发包 - YOLO合成数据生成器 - GLTF/GLB在线编辑 - 3D模型格式在线转换 - 可编程3D场景编辑器 - REVIT导出3D模型插件 - 3D模型语义搜索引擎 - AI模型在线查看 - Three.js虚拟轴心开发包 - 3D模型在线减面 - STL模型在线切割 - 3D道路快速建模

几个月前，我开始尝试一些涉及时空数据的项目构想。自然地，我想使用 Postgres 和 PostGIS 作为这些项目的主要工具，挑战在于如何将它们整合在一起。在经历了几次错误的开始后，我不得不将这些项目搁置一旁，留待其他优先事项。

在空闲时间，我继续阅读有关该主题的一些相关阅读材料。我发现Anita Graser 的文章为什么应该使用 PostGIS 轨迹，我建议你在继续本文之前阅读它。事实上，最好同时阅读评估 PostGIS 数据库中轨迹的时空数据模型！这些资源中包含大量信息，并带有更多指向其他资源的链接。

本文概述了如何将这些新的 PostGIS 轨迹技巧与我已经可用的 OpenStreetMap 数据一起使用的示例（加载和准备）。通常，轨迹示例假设使用从我们新时代的物联网传感器发送 GPS 点和时间戳收集的数据。此示例改为从数据建模的角度处理轨迹，展示如何使用 pgrouting 合成轨迹数据。数据可视化是共享信息的重要组成部分，QGIS 一直是我最喜欢与 PostGIS 数据一起使用的 GIS 应用程序。

1、路线数据

此示例的路线是使用手动选择的开始 (s_id) 和结束 (e_id) 节点创建的。这些 id 值与表 osm_routing.roads_noded_vertices_pgr 中的 pgrouting 节点数据对齐。为了为时间数据设置阶段，还包括开始时间 (s_time)。此示例为每条路线指定了早上 5:00 的开始时间，但可以轻松调整开始时间以探索影响。

route_id|s_id |e_id  |s_time             |
--------|-----|------|-------------------|
       1|24902| 21887|2020-11-01 05:00:00|
       2|43615| 82508|2020-11-01 05:00:00|
       3|83320| 25649|2020-11-01 05:00:00|
       4|27478| 46247|2020-11-01 05:00:00|
       5| 3487|113816|2020-11-01 05:00:00|
       6|20587|113816|2020-11-01 05:00:00|
       7|30505|113816|2020-11-01 05:00:00|

开始和结束节点（下面的 r.s_id 和 r.e_id）被送入 pgrouting 的 pgr_dijkstra() 路由函数。 osm_routing.roads 表包含由 PgOSM 准备的科罗拉多州丹佛市的 OpenStreetMap 数据。 pgosm.layer_detail 和 pgosm.routable 表用于轻松过滤机动交通（汽车）可用的道路。 osm_routing.roads_noded (edges) 表是用 pg_routing 准备的，它基于基本的 osm_routing.roads 表。

使用以下查询创建路线详细信息并将其保存到新表 (osm_routing.route_details) 中。

CREATE TABLE osm_routing.route_details AS
WITH rte AS (
SELECT * 
    FROM osm_routing.routes
)
SELECT r.route_id, d.*, n.the_geom AS node_geom, e.the_geom AS edge_geom
    FROM rte r
    INNER JOIN pgr_dijkstra(
        'SELECT rn.id, rn.source, rn.target, rn.cost_length AS cost,
                rn.the_geom
            FROM osm_routing.roads_noded rn
            INNER JOIN osm_routing.roads r ON rn.old_id = r.id
            INNER JOIN pgosm.layer_detail ld ON ld.code = r.code
            INNER JOIN pgosm.routable rbl
                ON rbl.layer_detail_id = ld.layer_detail_id
                    AND rbl.route_motor
            ',
        r.s_id, r.e_id, directed := False
        ) d
        ON True
    INNER JOIN osm_routing.roads_noded_vertices_pgr n ON d.node = n.id
    INNER JOIN osm_routing.roads_noded e ON d.edge = e.id
;

pgr_dijkstra() 函数为给定路线的每一步返回一行。 route_id 值链接回 osm_routing.routes。成本使用简单的基于长度的成本，以米为单位。节点和边缘值指向回 pgrouting 使用的各自表中的 id 值。

route_id|seq|path_seq|node |edge |cost              |agg_cost          |
--------|---|--------|-----|-----|------------------|------------------|
       1|  1|       1|24902|27505| 20.42259891997424|               0.0|
       1|  2|       2|27462|27506| 70.51812482510054| 20.42259891997425|
       1|  3|       3|27463|27507|21.724521060420614| 90.94072374507479|
       1|  4|       4|27464|27508|23.621809729839022| 112.6652448054954|
       1|  5|       5|27465|27509| 9.095234265004999|136.28705453533442|

每条路线的节点（蓝色圆圈）和边（红线）如下所示。我们可以通过代码 (ST_Intersects()) 轻松确认，并在地图上直观地确认这些路线使用简单的空间检查相交。如果也考虑时间，这些路线会相交吗？

2、轨迹之路

向几何添加时间的一种简单方法是为开始时间和结束时间添加两个 TIMESTAMPTZ 列。不幸的是，这种简单的方法很快就会崩溃，稍后的查询将说明这一点。幸运的是，PostGIS 扩展通过内置的轨迹支持（再次）来拯救。

什么是轨迹（Trajectory）？它至少从两个点开始，其中包括一个表示时间的度量 (M)。这些点可以使用 ST_MakePointM() 创建。包括 M 在内的点可以用 ST_MakeLine() 制作成一条线，现在你就有了一条轨迹。

为形成轨迹而添加的关键细节是 M，即测量（Measurement）。许多轨迹用例使用时间作为衡量标准。不幸的是，无法使用 TIMESTAMPTZ 格式，但在 Postgres 中使用纪元（epoch）是一种非常简单的解决方法。

3、分配速度

为了获得路线每一步的 M 值，我们需要每条边的行进速度。为 edge_mph 添加一列。

ALTER TABLE osm_routing.route_details ADD edge_mph NUMERIC;

相当多的连接才能希尔岛速度限制。将来，此步骤很可能会在数据处理管道中更早地完成。

UPDATE osm_routing.route_details rd_update
    SET edge_mph = rbl.max_speed
    FROM osm_routing.route_details rd
    INNER JOIN osm_routing.routes rte ON rd.route_id = rte.route_id
    INNER JOIN osm_routing.roads_noded rn ON rd.edge = rn.id
    INNER JOIN osm_routing.roads r ON rn.old_id = r.id
    INNER JOIN pgosm.layer_detail ld ON ld.code = r.code
    INNER JOIN pgosm.routable rbl ON rbl.layer_detail_id = ld.layer_detail_id
    WHERE rd_update.route_id = rd.route_id
        AND rd_update.seq = rd.seq
        AND rd_update.path_seq = rd.path_seq
;

4、计算每条边的时间

添加一列来存储每个节点的时间。

ALTER TABLE osm_routing.route_details ADD node_ts TIMESTAMPTZ;

要获得每个 node_ts 的值，我们需要根据边的速度以秒为单位计算成本。我们的距离以米为单位，我们的速度以 MPH 为单位，我们需要秒。功能的完美用例！

以下函数接受两个参数，length_m 和 miles_per_hour，并返回行进边缘所需的时间（以秒为单位）。

CREATE OR REPLACE FUNCTION osm_routing.mph_length_m_to_cost_seconds(
    length_m NUMERIC, miles_per_hour NUMERIC)
 RETURNS NUMERIC
 LANGUAGE sql
 SECURITY DEFINER
 SET search_path TO 'pg_temp'
AS $function$
/*
 * miles per hour * 0.44704  = meters per second
 */
    SELECT ($1 / ($2 * .44704))

$function$
;

这使得根据边缘的长度和速度计算遍历边缘的时间变得容易。以每小时 60 英里的速度行驶 1000 米需要多少秒？

SELECT osm_routing.mph_length_m_to_cost_seconds(1000, 60);

mph_length_m_to_cost_seconds|
----------------------------|
         37.2822715342400382|

以 20 mph 的速度跑 1000 米怎么样？

SELECT osm_routing.mph_length_m_to_cost_seconds(1000, 20);

mph_length_m_to_cost_seconds|
----------------------------|
        111.8468146027201145|

下面的 edge_times CTE 使用 mph_length_m_to_cost_seconds() 函数来计算每个线段的遍历时间。 time_to_node CTE 使用 SUM() 作为窗口函数来提供每个节点的总时间成本。最后的 UPDATE 查询将 node_ts 设置为路由的开始时间（来自 osm_routing.routes）加上前一个窗口函数的适当秒数。

WITH edge_times AS (
SELECT r.route_id, rd.seq, rd.path_seq, 
        r.s_time, rd.cost, rd.agg_cost, rd.edge_mph,
        CASE WHEN (rd.edge_mph > 0 AND agg_cost > 0)
            THEN osm_routing.mph_length_m_to_cost_seconds(rd.cost::NUMERIC, rd.edge_mph)
            ELSE NULL
            END AS edge_time_s
    FROM osm_routing.routes r
    INNER JOIN osm_routing.route_details rd ON True
    WHERE r.route_id = rd.route_id
), time_to_node AS (
SELECT SUM(et.edge_time_s)
            OVER (PARTITION BY et.route_id ORDER BY et.seq, et.path_seq)
            AS time_to_node_s,
        *
    FROM edge_times et
)
UPDATE osm_routing.route_details rd
    SET node_ts = ttn.s_time + (ttn.time_to_node_s || ' seconds')::INTERVAL
    FROM time_to_node ttn
    WHERE rd.route_id = ttn.route_id
        AND rd.seq = ttn.seq
        AND rd.path_seq = ttn.path_seq
        AND rd.node_ts IS NULL
;

上面的 UPDATE 查询没有设置路由的初始节点。以下更新修复了每条路线中的第一条记录等于 s_time。

UPDATE osm_routing.route_details rd
    SET node_ts = r.s_time
    FROM osm_routing.routes r
    WHERE rd.route_id = r.route_id
        AND rd.node_ts IS NULL
        AND rd.seq = 1
;

osm_routing.route_details 表现在包含包含路线的每个段的速度限制和时间戳。

route_id|seq|path_seq|edge_mph|node_ts            |
--------|---|--------|--------|-------------------|
       1|  1|       1|   45.00|2020-11-01 05:00:00|
       1|  2|       2|   45.00|2020-11-01 05:00:03|
       1|  3|       3|   45.00|2020-11-01 05:00:04|
       1|  4|       4|   45.00|2020-11-01 05:00:05|
       1|  5|       5|   45.00|2020-11-01 05:00:06|

5、最后，轨迹！

现在要转换为轨迹，我正在使用通过 CTE 创建的物化视图。第一部分（节点）提取每个节点所需的详细信息。第二步 (traj_nodes) 使用 M 值的时间戳的纪元创建点。第三步 (build_line) 使用 LEAD() 窗口函数将路线节点串在一起形成轨迹线。最后一部分计算一些额外的方便列，例如用于简单时间过滤的 time_range。

CREATE MATERIALIZED VIEW osm_routing.route_trajectory_segment AS
WITH nodes AS (
SELECT route_id, seq, path_seq, node, edge, cost, agg_cost,
        edge_mph, node_ts,
        ST_X(ST_Transform(node_geom, 4326)) AS long,
        ST_Y(ST_Transform(node_geom, 4326)) AS lat
    FROM osm_routing.route_details
), traj_nodes AS (
SELECT route_id, seq, path_seq, node_ts,
        ST_SetSRID(ST_MakePointM(long, lat, EXTRACT(epoch FROM node_ts)), 4326) AS traj_point
    FROM nodes
), build_line AS (
SELECT traj_nodes.route_id, traj_nodes.seq, traj_nodes.path_seq,
        node_ts AS time_start,
        LEAD(node_ts) OVER (PARTITION BY route_id ORDER BY seq, path_seq)
            AS time_end,
        ST_Transform(
            ST_MakeLine(traj_nodes.traj_point,
                LEAD(traj_nodes.traj_point)
                    OVER (PARTITION BY traj_nodes.route_id
                            ORDER BY traj_nodes.seq, traj_nodes.path_seq)
                        )
            , 2773)
                AS traj
    FROM traj_nodes
)
SELECT *,
        TSTZRANGE(time_start, time_end) AS time_range,
        ST_Length(traj) AS route_length,
        time_end - time_start AS duration
    FROM build_line
;

在轨迹和时间范围列中创建索引以快速查询。

CREATE INDEX gix_osm_routing_route_trajectory_segment
    ON osm_routing.route_trajectory_segment
    USING GIST (traj);
CREATE INDEX ix_osm_routing_route_trajectory_segment_trange
    ON osm_routing.route_trajectory_segment
    USING GIST (time_range);

虽然上面创建的实体化视图包括专用时间戳和时间范围计算，但这些值可以从轨迹数据本身中提取出来。

SELECT route_id, seq,
        TO_TIMESTAMP(ST_M(ST_StartPoint(traj))) AS t_start_traj, 
        TO_TIMESTAMP(ST_M(ST_StartPoint(traj))) AS t_end_traj
    FROM osm_routing.route_trajectory_segment
    ORDER BY route_id, seq, path_seq
    LIMIT 5
;


route_id|seq|t_start_traj       |t_end_traj         |
--------|---|-------------------|-------------------|
       1|  1|2020-11-01 05:00:00|2020-11-01 05:00:00|
       1|  2|2020-11-01 05:00:03|2020-11-01 05:00:03|
       1|  3|2020-11-01 05:00:04|2020-11-01 05:00:04|
       1|  4|2020-11-01 05:00:05|2020-11-01 05:00:05|
       1|  5|2020-11-01 05:00:06|2020-11-01 05:00:06|

6、检查有效性

使用 ST_IsValidTrajectory() 检查你的轨迹是否有无效数据是一个好习惯。我花了不止几次尝试才能正确构建最终轨迹数据。我的大部分问题都是由于我自己在排序和其他逻辑方面的错误。

7、时间和空间的交叉点

此查询是如何使用时间范围数据伪造轨迹的示例。如前所述，除了最基本的示例之外，此方法无法按预期工作。连接使用 a.time_range && b.time_range AND ST_DWithin(a.traj, b.traj, 2) 来过滤时间范围重叠且几何形状彼此相距 2 米的行。使用这种方法返回 45 行，其中大部分是误报。

SELECT a.route_id AS route_id_a, b.route_id AS route_id_b,
        a.seq AS seq_a, b.seq AS seq_b,
        a.time_range, b.time_range,
        a.traj, b.traj AS traj_line_intersect
    FROM osm_routing.route_trajectory_segment a
    INNER JOIN osm_routing.route_trajectory_segment b
        ON a.time_range && b.time_range 
            AND ST_DWithin(a.traj, b.traj, 2)
            AND a.route_id < b.route_id
;

仔细查看一条记录，我们可以看到一条路线的尾部与另一条路线的起点位于同一点 (POINT (-104.99831660000001 39.74985710023273))。不同的路线（05:00:05 与 05:00:11）相隔 6 秒穿过这些节点。沿路线行进的物体之间的距离真的在 2 米以内吗？

SELECT a.route_id AS route_id_a, b.route_id AS route_id_b,
        a.seq AS seq_a, b.seq AS seq_b,
        ST_Transform(ST_StartPoint(a.traj), 4326) AS start_a,
        ST_Transform(ST_EndPoint(b.traj), 4326) AS end_b,
        TO_TIMESTAMP(ST_M(ST_StartPoint(a.traj))) AS start_time_a,
        TO_TIMESTAMP(ST_M(ST_EndPoint(b.traj))) AS end_time_b
    FROM osm_routing.route_trajectory_segment a
    INNER JOIN osm_routing.route_trajectory_segment b
        ON a.time_range && b.time_range 
            AND ST_DWithin(a.traj, b.traj, 2)
            AND a.route_id < b.route_id
            AND (a.route_id = 3 AND b.route_id = 4 AND a.seq = 3 AND b.seq = 6)
;

Name        |Value                                        |
------------|---------------------------------------------|
route_id_a  |3                                            |
route_id_b  |4                                            |
seq_a       |3                                            |
seq_b       |6                                            |
start_a     |POINT (-104.99831660000001 39.74985710023273)|
end_b       |POINT (-104.99831660000001 39.74985710023273)|
start_time_a|2020-11-01 05:00:05                          |
end_time_b  |2020-11-01 05:00:11                          |

要检查由轨迹表示的两个移动物体之间的最近距离，我们可以使用 ST_DistanceCPA。以下查询使用 ST_DistanceCPA(a.traj, b.traj) AS route_cpa_meters 更新。这说明两条路线的最近点竟然相距3.9米，远远超过了2米的门槛！因此，虽然时间重叠，并且几何形状在 2 米以内，但移动的物体在沿着各自的路线移动时不会在 2 米以内。

SELECT a.route_id AS route_id_a, b.route_id AS route_id_b,
        a.seq AS seq_a, b.seq AS seq_b,
        ST_DistanceCPA(a.traj, b.traj) AS route_cpa_meters
    FROM osm_routing.route_trajectory_segment a
    INNER JOIN osm_routing.route_trajectory_segment b
        ON a.time_range && b.time_range 
            AND ST_DWithin(a.traj, b.traj, 2)
            AND a.route_id < b.route_id
            AND a.route_id = 3 AND b.route_id = 4
            AND a.seq = 3 AND b.seq = 6
;

route_id_a|route_id_b|seq_a|seq_b|route_cpa_meters |
----------|----------|-----|-----|-----------------|
         3|         4|    3|    6|3.938249529587584|

8、轨迹特定函数

为了更准确地进行时空过滤，请使用轨迹特定函数 ST_CPAWithin。这用 ST_CPAWithin(a.traj, b.traj, 2) 替换了之前查询中关于时间范围和几何的两部分连接。这只返回 3 条记录，而不是 45 条。这还包括一个检查列，以显示所有三个轨迹段的时间都在 2 米或更短的范围内。

SELECT a.route_id AS route_id_a, b.route_id AS route_id_b,
        a.seq AS seq_a, b.seq AS seq_b,
        ST_DistanceCPA(a.traj, b.traj) AS traj_closest_point_meters,
        a.traj, b.traj AS traj_line_intersect
    FROM osm_routing.route_trajectory_segment a
    INNER JOIN osm_routing.route_trajectory_segment b
        ON ST_CPAWithin(a.traj, b.traj, 2)
            AND a.route_id < b.route_id
;

route_id_a|route_id_b|seq_a|seq_b|traj_closest_point_meters|
----------|----------|-----|-----|-------------------------|
         3|         4|   10|   13|         1.87513147035179|
         3|         4|   10|   14|       1.8751314051078767|
         3|         4|   11|   14|       1.8751314051078767|

9、在 QGIS 中制作动画

查看实际数据总是有帮助的。 QGIS 现在包含一个名为 Temporal Controller 的功能，在 QGIS 3.14 中可用。和更新。我正在使用 QGIS v3.16 当前的功能，效果很好。 Nyall Dawson 有一个关于如何使用此功能的很棒的教程。上面创建的轨迹路线在下面使用 QGIS 的时间控制器创建的动画图像中进行了说明。

从 QGIS 保存动画目前似乎只有一个选项，将一系列 PNG 文件保存到一个目录中。为了将它们拼接成动画 GIF，我使用了转换工具。我无法获得与图像一起导出的累积范围功能，并且获得正确的分辨率和范围设置仍然有点挑剔。也就是说，这是一项新功能，我相信这些细微差别将来会得到改进。

现在，将导出的 PNG 文件转换为动画 GIF。

cd /path/to/qgis/dir
convert -delay 10 -loop 0 *.png postgis-trajectory-qgis-animation.gif

10、结束语

这篇文章展示了如何使用 OpenStreetMap 数据和 pgrouting 开始使用 PostGIS 轨迹数据。既然时空建模的基础知识已经到位，还有很多探索要做。

原文链接：PostGIS Trajectory: Space plus Time

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

PostGIS时空轨迹实战

1、路线数据

2、轨迹之路

3、分配速度

4、计算每条边的时间

5、最后，轨迹！

6、检查有效性

7、时间和空间的交叉点

8、轨迹特定函数

9、在 QGIS 中制作动画

10、结束语

admin

PostGIS最近邻搜索内幕

Geometry vs. Geography

1、路线数据

2、轨迹之路

3、分配速度

4、计算每条边的时间

5、最后，轨迹！

6、检查有效性

7、时间和空间的交叉点

8、轨迹特定函数

9、在 QGIS 中制作动画

10、结束语

PostGIS最近邻搜索内幕

Geometry vs. Geography

You might also like...

You might also like...