基于pgrouting的路径分析

一.技术背景，相关技术介绍

      PgRouting是基于开源空间数据库PostGIS用于网络分析的扩展模块，最初它被称作pgDijkstra，因为它只是利用Dijkstra算法实现最短路径搜索，之后慢慢添加了其他的路径分析算法，如A算法，双向A算法，Dijkstra算法，双向Dijkstra算法，tsp货郎担算法等，然后被更名为pgRouting[1]。该扩展库依托PostGIS自身的gist索引，丰富的坐标系与图形类型，强大的几何处理能力，如空间查询，空间处理，线性参考等优势，能保障在较大数据级别下的网络分析效果更快更好。
　　PostGIS早已奠定了最优秀的开源空间数据库地位，在新时代GIS中的应用将会越来越普遍。其实，网络分析算法很多服务端语言如java，C#等虽能实现，但基于真实城市道路数据量较大且查询分析操作步骤复杂与数据库交互频繁，以这类服务端频繁访问数据库导致数据库开销压力较大，分析较慢，故选择PgRouting在数据库内部实现算法，提升分析效率。最后，路径分析不仅仅是最短路径，在实际应用中还有最短耗时，最近距离，道路对车辆类型限制，道路对速度限制等因素，交通事故、市政事故导致的交通障碍点等问题，所有的问题本质其实是对路径分析权重（Weight）的设置问题。

二.数据准备

1.准备路网数据，从OSM下载中国区的数据，然后使用ArcMap对数据进行裁剪，获得自己想要的区域路网数据；

2.通过Postgis自带的shp导入工具导入被裁剪好的路网数据（注意shp路径不要太深，不要还有中文；对于中文乱码，可以设置GBK等转码）；

3.设置导入的SRID为3857进行坐标系转化；

4.勾选Options里面的"Generate simple geometries instead of MULTI geometries",因为路径分析只支持LineString类型，不支持MuliLineString类型，这里介绍一个方法，因为裁剪的过程或者数据自身问题，当勾选这个的时候可能会出现导入数据失败，当出现这个情况的时候，先不勾选此选项导入到数据库中，然后通过"SELECT gid from road_hz_two  where get_txt_count(st_astext(geom),'(')>2;"获得多线数据，手动删除这些数据，然后将数据通过shp2sql导出程shp数据，然后在重复上诉过程，实现支持路网分析的数据的导入；

5.查看刚才导入的表的数据，增加source,target字段，后面分析需要用到；

--在road_two中添加source,target字段
alter TABLE road_two add COLUMN source int;
alter table road_two add COLUMN target int;

6.通过pgrouting提供的pgr_createTopology方法，对道路数据创建拓扑关系

--创建连通性topo
--road_two是表名称，geom是该表的图形字段名称，gid是改变的主键
--一般我们使用shp2pgsql工具会自动创建gid为主键，geom为图形，如果是其他工具，注意对应字段
select pgr_createtopology('road_two', 0.0001, 'geom', 'gid');

7.对连通性字段建立索引，增加查询速度

create index road_source_index on road_two("source");
create index road_target_index on road_two("target");

8.到此为止，数据的准备工作完成；

三.路径分析

相关说明：osm下载的路网数据，里面包含"oneway"道路方向的说明 ，"B"代表双向，"T"代表仅反向，  "F"代表仅正向；

3.1.道路成本权重说明

在算法中分为有向图，无向图，图的长度一般设置为权重，路网分析中，具体到此交通领域，也分为双向通行道路，单向通行道路，交通师父导致的临时交通阻塞无法通行（障碍点），不同等级道路对车辆类型限制，以下是不同条件下如何设置通行成本权重的相关示例：

• 双向通行

update road_two set length=st_length(geom),rev_length=st_length(geom) where oneway='B';
--采用真实地理距离是这样：
update road_two set lenght=st_length(st_transform(geom),4326),true),rev_length=st_length((st_transform(geom),4326),true) where oneway='B';

• 单向通行

#FT是道路方向与数字化方向一致，那么正向通行成本为道路长度，反向成本为正无穷（以极大值代替）
update road_two set length=st_length(geom),rev_length=99999999999 where oneway='F';
update road_two set length=99999999999,rev_length=st_length(geom) where oneway='T';

• 障碍点

#假设gid=20的道路因事故，修路暂时不能通行
update road_two set lenght=99999999999,rev_length=99999999999 where gid=20;

• 限制通行

假设当前是一辆大货车，通过有限高限重的道路，在为他做规划时，先获取车辆类型，再查询road表中是否有对其限制的因素（以下纯逻辑描述sql）

#假设道路表有字段restrict，该字段是array，记录了不可通行的车辆类型
update road_two set lenght=99999999999,rev_length=99999999999 where 'lorry'=any(restrict);

 3.2设置道路权重

1.本文只涉及单向，双向等简单情况，复杂情况不考虑。

通过以下方式设置道路权重

--FT是道路方向与数字化方向一致，那么正向通行成本为道路长度，反向成本为正无穷（以极大值代替）
UPDATE road_two SET cost = CASE WHEN oneway='B'THEN length --双向
WHEN oneway='T' THEN 999999999 -- 反向
WHEN oneway='F' THEN length -- 正向
ELSE length END;
--------------------- 
UPDATE road_two SET rev_cost = CASE WHEN oneway='B' THEN length --双向
WHEN oneway='T' THEN length -- 反向
WHEN oneway='F' THEN 999999999 -- 正向
ELSE length END;
--------------------- 

2.计算最短路径

• 最短路径，不包含方向

SELECT line.*,pt.geom FROM pgr_dijkstra(
    'SELECT gid AS id,
         source,
         target,
         cost,  rev_cost::double precision AS reverse_cost
        FROM road_two',
    1060, 1661,false,false) line LEFT JOIN road_two pt on line.id2=pt.gid;

--通过id转edge来实现对边表的对应
drop table if exists dijkstra_res_edge_dretrue;
SELECT seq,id1 AS node, id2 AS edge,line."cost",pt.geom  into  dijkstra_res_node_point_true FROM pgr_dijkstra(
    'SELECT gid AS id,
         source,
         target,
         cost, 
         rev_cost::double precision AS reverse_cost
        FROM road_two',
    1,10,true,false) line LEFT JOIN road_two pt on line.id2=pt.gid;

--通过edge直接对应边表
drop table if exists dijkstra_res_node_dretrue;
SELECT line.* ,pt.the_geom as geom  into dijkstra_res_true_one  FROM pgr_dijkstra(
    'SELECT gid AS id,
         source,
         target,
         cost
        FROM road_two',
    1,10,true) line LEFT JOIN road_two_vertices_pgr pt on line.node=pt.id;

--以上两种方式都是包含方向的，但第四个参数的使用还需要查明-todo

• 最短路径，包含方向

SELECT seq,id1 AS node, id2 AS edge,line."cost",pt.geom   FROM pgr_dijkstra(
    'SELECT gid AS id,
         source,
         target,
         cost, 
         rev_cost::double precision AS reverse_cost
        FROM road_two',
    1, 10,true,true) line LEFT JOIN road_two pt on line.id2=pt.gid;

按顺序经过多个点的最短路径：从1->10->20（待完善）

SELECT line.*,pt.geom      FROM pgr_dijkstraVia(
    'SELECT gid AS id,
         source,
         target,
         cost, 
         rev_cost::double precision AS reverse_cost
        FROM road_two order by gid',
    ARRAY[1,10,20],true,false,true) line LEFT JOIN road_two pt on line.edge=pt.gid;

关于pgr_dijkstraVia，pgr_dijkstra的详细用法可以参看pgrouting官网

3.结果查看

可以将查询到的数据生成到一个新表，然后通过udig查看路径分析的结果是否正确

 

参考资料

1.https://www.jianshu.com/p/34c8378c3da9

2.https://blog.csdn.net/redtomemory/article/details/78293447

3.https://workshop.pgrouting.org/2.5/en/chapters/shortest_path.html