用Postgis算最短路径（在任意位置选择起点终点）

1.废话几句

1.前言

阅读本文需要知道什么是shapfile，什么是路径分析，什么是GIS。相比Arcgis的路径分析功能，本文介绍的方法稍微复杂，需要注意的细节更多，但却是完全免费的。PostGis+QGIS+Geoserver开源Gis三剑客用起来真的很舒服。

2.软件环境

pgRouting.png

Postgresql+Postgis+pgRouting+QGIS。

2.数据准备

1.将道路Shp数据导入Postgis

打开PostGIS Shapefile Import/Export Manager，导入线路数据，你的数据只要是单线段类型即可，在导入过程中注意填写SRID、字符编码、并勾选上单体图形，如果导入过程中报错无法生成单体图形，则需要将图形处理成单体图形，可以借助ArcMap的FeatureToLine功能。

image.png

2.数据表修改

现在为刚刚导入的道路数据添加路径分析所需的字段，其中source指该路段的源头（道路节点），target指该路段的终点（道路节点）：

ALTER TABLE daolu_simp
ADD COLUMN source integer,
ADD COLUMN target integer,
ADD COLUMN cost_len double precision,
ADD COLUMN cost_time double precision,
ADD COLUMN rcost_len double precision,
ADD COLUMN rcost_time double precision,
ADD COLUMN x1 double precision,
ADD COLUMN y1 double precision,
ADD COLUMN x2 double precision,
ADD COLUMN y2 double precision,
ADD COLUMN to_cost double precision,
ADD COLUMN rule text,
ADD COLUMN isolated integer;

接下来我们把新增的内容填上，这里主要填道路往返的路程花费和起点终点坐标，时间花费和其他暂时不填：

UPDATE daolu_simp SET x1 = st_x(st_startpoint(geom)),
y1 = st_y(st_startpoint(geom)),
x2 = st_x(st_endpoint(geom)),
y2 = st_y(st_endpoint(geom)),
cost_len = st_length_spheroid(geom, 'SPHEROID["WGS84",6378137,298.25728]'),
rcost_len = st_length_spheroid(geom, 'SPHEROID["WGS84",6378137,298.25728]')

这里有一份可直接用的路网数据（不需要从上一步导入）：

CREATE TABLE edge_table (
id BIGSERIAL,
dir character varying,
source BIGINT,
target BIGINT,
cost_lenFLOAT,
rcost_len FLOAT,
capacity BIGINT,
reverse_capacity BIGINT,
category_id INTEGER,
reverse_category_id INTEGER,
x1 FLOAT,
y1 FLOAT,
x2 FLOAT,
y2 FLOAT,
the_geom geometry
);

INSERT INTO edge_table (
category_id, reverse_category_id,
cost_len, rcost_len,
capacity, reverse_capacity,
x1, y1,
x2, y2) VALUES
(3, 1, 1, 1, 80, 130, 2, 0, 2, 1),
(3, 2, -1, 1, -1, 100, 2, 1, 3, 1),
(2, 1, -1, 1, -1, 130, 3, 1, 4, 1),
(2, 4, 1, 1, 100, 50, 2, 1, 2, 2),
(1, 4, 1, -1, 130, -1, 3, 1, 3, 2),
(4, 2, 1, 1, 50, 100, 0, 2, 1, 2),
(4, 1, 1, 1, 50, 130, 1, 2, 2, 2),
(2, 1, 1, 1, 100, 130, 2, 2, 3, 2),
(1, 3, 1, 1, 130, 80, 3, 2, 4, 2),
(1, 4, 1, 1, 130, 50, 2, 2, 2, 3),
(1, 2, 1, -1, 130, -1, 3, 2, 3, 3),
(2, 3, 1, -1, 100, -1, 2, 3, 3, 3),
(2, 4, 1, -1, 100, -1, 3, 3, 4, 3),
(3, 1, 1, 1, 80, 130, 2, 3, 2, 4),
(3, 4, 1, 1, 80, 50, 4, 2, 4, 3),
(3, 3, 1, 1, 80, 80, 4, 1, 4, 2),
(1, 2, 1, 1, 130, 100, 0.5, 3.5, 1.999999999999,3.5),
(4, 1, 1, 1, 50, 130, 3.5, 2.3, 3.5,4);

3.道路拓扑计算

拓扑分析函数的签名如下：

varchar pgr_createTopology(text edge_table, double precision tolerance,
text the_geom:='the_geom', text id:='id',
text source:='source',text target:='target',
text rows_where:='true', boolean clean:=false)

手册中对各个参数的说明：

edge_table text 路网表的名字. (也可以包含Schema)
tolerance float8 决定路段不连通的阈值. (投影坐标的单位)>
the_geom text 路网数据中道路的图形数据. 默认是the_geom.
id text 路网的Primary key.默认为id.
source text 路段的起点节点. 默认为 source.
target text 路段的终点节点. 默认为target.
rows_where text Condition to SELECT a subset or rows. Default value is true to indicate all
rows that where source or target have a null value, otherwise the condition is used.
clean boolean 是否清除上一次该路网的拓扑数据.默认false.

此函数对路网进行拓扑分析，会生成一个路网节点表，并对路网表中的source和target字段填充入节点编号。
这里我们运行这一句：

SELECT pgr_createTopology('daolu_simp', 0.0001, 'geom', 'gid');

结果为OK说明成功运行。

image.png

4.拓扑检查

拓扑分析检查函数签名：

varchar pgr_analyzeGraph(text edge_table, double precision tolerance,
text the_geom:='the_geom', text id:='id',
text source:='source',text target:='target',text rows_where:='true')

本例我们这样填写：

select pgr_analyzeGraph('daolu_simp', 0.0001, 'geom', 'gid')

检查结果如下：

image.png

3.路径分析

1.pgr_withPoints函数

pgr_dijkstra是pgRouting官方教程中第一个介绍的最短路径规划方法，pgr_withPoints和 pgr_dijkstra一样也是最短路径规划方法，输出结果也基本一致，但是两者的输入不一样。pgr_dijkstra需要输入的是路网节点编号，也就是说只能从道路的端点算起，而pgr_withPoints可以从道路的中间位置选择起点和终点。
来看看最简单使用方法的函数签名：

r_withPoints(edges_sql, points_sql, start_vid, end_vid)

参数说明：edges_sql是从路网表选择出需要参与计算的路径的Sql字符串，points_sql是从POI表（下文会说明）选择出途径点的Sql字符串，start_vid是起点的节点id，正值为路网节点，负值为POI表数据，end_vid为终点节点，同样正值为路网节点，负值为POI表数据。注意在edges_sql中需要指定cost和reverse_cost（去的花费和来的花费），函数通过这两个值的正负情况来获知道路的通过性（单行、双向、不通）

这里我们举一个例子：（一个起点一个终点）

SELECT * FROM pgr_withPoints(
'SELECT id, source, target, cost, reverse_cost FROM edge_table ORDER BY id',
'SELECT pid, edge_id, fraction, side from pointsOfInterest',
-1, -3);

seq | path_seq | node | edge | cost | agg_cost
-----+----------+------+------+------+----------
1 | 1 | -1 | 1 | 0.6 | 0
2 | 2 | 2 | 4 | 1 | 0.6
3 | 3 | 5 | 10 | 1 | 1.6
4 | 4 | 10 | 12 | 0.6 | 2.6
5 | 5 | -3 | -1 | 0 | 3.2
(5 rows)

输出结果中path_seq为路径（点位）顺序，node为途经节点（正值为路网节点，负值为POI节点），edge为路网的路段，cost为该路段的花费。需要注意的是所有的edge加起来可能会多于实际经过路径（如果从POI节点出发），这时需要将一头一尾（结果为-1的edge排除在外）两个路段单独取出来根据fraction计算实际经过部分，其他路段则为完整实际经过路段。

2.POI表

POI表在路径分析功能中的角色是，保存用户选择的途径点（起点和终点）并进行必要的处理，表结构sql：

CREATE TABLE pointsOfInterest(
pid BIGSERIAL,
x FLOAT,
y FLOAT,
edge_id BIGINT,
side CHAR,
fraction FLOAT,
the_geom geometry,
newPoint geometry
);

需要注意三个字段：edge_id表示该POI在路网中距离最近（或者其他条件）的路径编号，fraction表示该路径上距离POI最近的位置，the_geom表示POI的图形，newPoint表示POI在路网中对应的位置（fraction表示的位置）。

image.png
通过下面的sql来填写POI表，其中106.72545, 26.56554（sql中出现了4对）经纬度为POI的原始位置，ST_DWithin函数用于选择距离POI最近的路径边，参数1000表示在1000米内搜索；ST_LineLocatePoint函数的作用是计算该路径边上距离POI最近的位置，得到的结果更新表中的fracion字段；ST_LineInterpolatePoint函数的作用则是根据fraction将路径上实际的这个点位获取出来，得到的结果更新表中newPoint字段。

INSERT INTO pointsOfInterest(x, y, edge_id, side, fraction,the_geom,newPoint) 
(
SELECT 106.72545, 26.56554,gid,'b',ST_LineLocatePoint(geom, concat('SRID=4326;POINT(106.72545 26.56554)')::geometry) as fraction,st_makePoint(106.72545,26.56554),ST_LineInterpolatePoint(geom, ST_LineLocatePoint(geom, concat('SRID=4326;POINT(106.72545 26.56554)')::geometry))
FROM daolu_simp
WHERE ST_DWithin(geom, concat('SRID=4326;POINT(106.72545 26.56554)')::geometry, 1000)
order by geom <-> concat('SRID=4326;POINT(106.72545 26.56554)')::geometry
limit 1) 

POI点导入后可以看到有POI点的原始点位（绿色）和对应在路径上的点位（紫色）

image.png

3.开始路径分析

为了方便，我们使用QGIS调试路径分析效果，打开QGIS，连接PostGis数据库，打开数据

image.png
image.png
再次导入一个线路图层，对该图层设置Sql过滤，用于展示规划的线路
image.png
image.png
路径查询Sql如下所示，采用的展示方案是将途径路段的首尾连接起来，核心函数pgr_withPoints的用法上文已经介绍，注意在使用st_makeline生成线段的时候需要对节点按照seq或seq_path进行排序。

select st_makeline(point_geom)
from
(SELECT seq, node ,newpoint as point_geom FROM pgr_withPoints(
'SELECT gid as id, source, target, cost_len as cost, rcost_len as reverse_cost FROM daolu_simp ORDER BY gid',
'SELECT pid as pid, edge_id, fraction, side from pointsOfInterest',
-1, -2) as wp
join pointsofinterest poi
on -wp.node= poi.pid 
union all
SELECT seq, node ,the_geom as point_geom FROM pgr_withPoints(
'SELECT gid as id, source, target, cost_len as cost, rcost_len as reverse_cost FROM daolu_simp ORDER BY gid',
'SELECT pid as pid, edge_id, fraction, side from pointsOfInterest',
-1, -2) as wp
join daolu_simp_vertices_pgr et
on wp.node= et.id

order by seq)as a

计算结果（绿色线条）：

image.png

4.最后

到此为止，我们已经能够计算出路径了（实际上完整的路径还需要另写Sql获取），至于是不是最短路径，需要保证路网数据质量，例如路网连通性、路径花费的合理性等，那些都是数据层面的问题，在程序方面我们已经实现了最基础的版本。若路网存在转弯限制、单行道限制（本文中提到了相应字段，但是未使用）等，还需要我们进一步探讨。