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

几个月前,我开始尝试一些涉及时空数据的项目构想。 自然地,我想使用 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_detailpgosm.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翻译整理,转载请标明出处