postgis+geoserver最短路径

1 安装软件2 数据预处理3 操作postgres4 操作Geoserver5 计算最短路径6 问题记录

1 安装软件

安装PostgreSQL与 postgis，以及geoserver。

2 数据预处理

对路网的shapefile数据进行处理：

1、打断相交线。 2、拓扑检查。

以上处理均可以在ArcMap中实现，在此不再详细列出操作步骤。注意如果不打断相交线或不进行拓扑检查，很有可能导致后面无法进行最短路径的计算，建议可以先用简单数据进行实验。

本实验所用的数据为自己绘制的简单数据，数据及属性表如下图所示 ：

3 操作postgres

打开pgAdmin，新建名为simple_road的数据库。 在simple_road数据库上右键打开CREATE Script，对simple_road数据库增加如下的拓展： CREATE EXTENSION postgis; CREATE EXTENSION postgis_topology; CREATE EXTENSION fuzzystrmatch; CREATE EXTENSION postgis_tiger_geocoder; CREATE EXTENSION pgRouting; CREATE EXTENSION address_standardizer;

打开postgis shapefile import/export manager工具，进行如下图所示的4个步骤： （1）点击“view connection details”连接至postgres的simple_road数据库。 （2）点击“add file”选择shapefile文件上传，注意设置SRID为4326，Geo_Colunm为geom。 （3） 点击“options”进入设置，设置编码为“GBK”以支持中文字段，勾选三个选项，注意必须要勾选最后一个选项，将MULIPOLYLINE简化为POLYLINE。 （4）点击“import”导入数据后，即可在postgres的simple_road数据库中查看到导入的数据。

继续在CREATE Script中输入代码进行以下的步骤：

（1） 添加路网线段的起止点，存储线段的首尾编号

ALTER TABLE simple_road ADD COLUMN source integer; ALTER TABLE simple _road ADD COLUMN target integer;

（2）添加路网道路权重值

ALTER TABLE simple_road ADD COLUMN length double precision;

（3）为权重赋值，这里将路段的长度赋值给权重值

UPDATE simple_road SET length = ST_Length(geom);

（4） 回程成本设置，则可支持回程

ALTER TABLE simple_road ADD COLUMN reverse_cost double precision; UPDATE simple_road SET reverse_cost = length;

（5） 创建拓扑结构（数据量较大时，这一步骤需要等待较长）

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

（6） 为source和target字段创建索引

CREATE INDEX source_idx ON simple_road ("source"); CREATE INDEX target_idx ON simple_road ("target");

（7） 创建起点经度x1、起点纬度y1、终点经度x2、终点纬度y2

ALTER TABLE simple_road ADD COLUMN x1 double precision; ALTER TABLE simple_road ADD COLUMN y1 double precision; ALTER TABLE simple_road ADD COLUMN x2 double precision; ALTER TABLE simple_road ADD COLUMN y2 double precision;

（8） 给x1、y1、x2、y2复制

UPDATE simple_road SET x1 =ST_x(ST_PointN(geom, 1)); UPDATE simple_road SET y1 =ST_y(ST_PointN(geom, 1)); UPDATE simple_road SET x2 =ST_x(ST_PointN(geom, ST_NumPoints(geom))); UPDATE simple_road SET y2 =ST_y(ST_PointN(geom, ST_NumPoints(geom))); 计算任意两点间最短路径的主要思路为：

首先计算距离起点A和终点B最新的路网线段l1和l2，再使用pgr_dijkstra函数计算l1和l2的两线段端点之间的最短路径，也就是l1的起点和终点至l2的起点和终点四种情况的最短路径，求得这四种情况中线段长度最小的情况，然后再在l1求得离起点A最近的点C，以及在l2上离终点B最近的点D，并根据点C、D截取l1和l2，最后将起点A、点C、最短路径、点D、终点B进行拼接，就是起点A到终点B的最短路径。

输入以下代码创建pgr_fromAtoB函数：

-- Function: pgr_fromAtoB(character varying, double precision, double precision, double precision, double precision) -- DROP FUNCTION pgr_fromAtoB(character varying, double precision, double precision, double precision, double precision); CREATE OR REPLACE FUNCTION pgr_fromAtoB( tbl character varying, startx double precision, starty double precision, endx double precision, endy double precision) RETURNS geometry AS $BODY$ declare startpoint geometry; -- 起点 endpoint geometry; -- 终点 star_line geometry; -- 起点到最近线上点的线段 end_line geometry; -- 终点到最近线上点的线段 v_startLine geometry;--离起点最近的线 v_endLine geometry;--离终点最近的线 v_startTarget integer;--距离起点最近线的终点 v_startSource integer; v_endSource integer;--距离终点最近线的起点 v_endTarget integer; v_statpoint geometry;--在v_startLine上距离起点最近的点 v_endpoint geometry;--在v_endLine上距离终点最近的点 v_res geometry;--最短路径分析结果 v_res_a geometry; v_res_b geometry; v_res_c geometry; v_res_d geometry; v_perStart float; --v_statpoint在v_res上的百分比 v_perEnd float; --v_endpoint在v_res上的百分比 v_shPath_se geometry;--开始到结束 v_shPath_es geometry;--结束到开始 v_shPath geometry;--最终结果 tempnode float; begin --查询离起点最近的线 --4326坐标系 --找起点10000000米范围内的最近线 execute 'select geom, source, target from ' ||tbl|| ' where ST_DWithin(geom,ST_Geometryfromtext(''point('||startx ||' ' || starty||')'',4326), 10000000) order by ST_Distance(geom,ST_GeometryFromText(''point('|| startx ||' '|| starty ||')'',4326)) limit 1' into v_startLine, v_startSource ,v_startTarget; --查询离终点最近的线 --找终点10000000米范围内的最近线 execute 'select geom, source, target from ' ||tbl|| ' where ST_DWithin(geom,ST_Geometryfromtext(''point('|| endx || ' ' || endy ||')'',4326), 10000000) order by ST_Distance(geom,ST_GeometryFromText(''point('|| endx ||' ' || endy ||')'',4326)) limit 1' into v_endLine, v_endSource,v_endTarget; --如果没找到最近的线，就返回null if (v_startLine is null) or (v_endLine is null) then return null; end if ; select ST_ClosestPoint(v_startLine, ST_Geometryfromtext('point('|| startx ||' ' || starty ||')',4326)) into v_statpoint; select ST_ClosestPoint(v_endLine, ST_GeometryFromText('point('|| endx ||' ' || endy ||')',4326)) into v_endpoint; --raise notice '%', v_startSource; -- 12 --raise notice '%', v_endSource; -- 3 --raise notice '%', v_startTarget; -- 1 --raise notice '%', v_endTarget; -- 6 -- 从开始的起点到结束的起点最短路径 execute 'SELECT st_linemerge(st_union(b.geom)) ' || 'FROM pgr_dijkstra( ''SELECT gid as id, source, target, length as cost FROM ' ||tbl|| '''::text,' ||'array['||v_startSource||'] , ' ||'array['||v_endSource||'] , false ) a, ' ||tbl|| ' b WHERE a.edge=b.gid' into v_res ; --从开始的终点到结束的起点最短路径 execute 'SELECT st_linemerge(st_union(b.geom)) ' || 'FROM pgr_dijkstra( ''SELECT gid as id, source, target, length as cost FROM ' ||tbl|| '''::text,' ||'array['||v_startTarget||'] , ' ||'array['||v_endSource||'] , false ) a, ' ||tbl|| ' b WHERE a.edge=b.gid' into v_res_b ; --从开始的起点到结束的终点最短路径 execute 'SELECT st_linemerge(st_union(b.geom)) ' || 'FROM pgr_dijkstra( ''SELECT gid as id, source, target, length as cost FROM ' ||tbl|| '''::text,' ||'array['||v_startSource||'] , ' ||'array['||v_endTarget||'] , false ) a, ' ||tbl|| ' b WHERE a.edge=b.gid' into v_res_c ; --从开始的终点到结束的终点最短路径 execute 'SELECT st_linemerge(st_union(b.geom)) ' || 'FROM pgr_dijkstra( ''SELECT gid as id, source, target, length as cost FROM ' ||tbl|| '''::text,' ||'array['||v_startTarget||'] , ' ||'array['||v_endTarget||'] , false ) a, ' ||tbl|| ' b WHERE a.edge=b.gid' into v_res_d ; if(ST_Length(v_res) > ST_Length(v_res_b)) then v_res = v_res_b; end if; if(ST_Length(v_res) > ST_Length(v_res_c)) then v_res = v_res_c; end if; if(ST_Length(v_res) > ST_Length(v_res_d)) then v_res = v_res_d; end if; --如果找不到最短路径，就返回null --if(v_res is null) then -- return null; --end if; --将v_res,v_startLine,v_endLine进行拼接 select st_linemerge(ST_Union(array[v_res,v_startLine,v_endLine])) into v_res; select st_linelocatepoint(v_res, v_statpoint) into v_perStart; select st_linelocatepoint(v_res, v_endpoint) into v_perEnd; if(v_perStart > v_perEnd) then tempnode = v_perStart; v_perStart = v_perEnd; v_perEnd = tempnode; end if; --截取v_res --拼接线 SELECT st_linesubstring(v_res,v_perStart, v_perEnd) into v_shPath; --拼接起点与终点 select ST_SetSRID( ST_MakePoint(startx , starty), 4326 )into startpoint; select ST_SetSRID( ST_MakePoint(endx , endy), 4326 )into endpoint; select ST_MakeLine( v_statpoint,startpoint) into star_line; select ST_MakeLine( endpoint,v_endpoint) into end_line; select st_linemerge(st_union(array[end_line,v_shPath,star_line])) into v_shPath; return v_shPath; end; $BODY$ LANGUAGE plpgsql VOLATILE STRICT COST 100; ALTER FUNCTION pgr_fromctod(character varying, double precision, double precision, double precision, double precision) OWNER TO postgres;

4 操作Geoserver

添加新的数据存储，类型为postgis数据库：

数据源名称为simple_road： 连接的数据库为simple_road，user和password为连接数据库的用户及密码： 其他选项为默认配置，点击最下方的“保存”按钮。

保存后，选择“配置新的SQL视图”， 输入视图名称为“new”，SQL语句为 select * from pgr_fromAtoB('simple_road',%x1%, %y1%, %x2%, %y2%)

然后点击下方的“从SQL猜想的参数”，输入默认值为0，验证的正则表达式为*^-?[\d.]+$* 点击下方的“刷新”按钮，在新弹出的栏目中，选择类型为“LineString”，SRI为“4326”： 然后计算数据的边框： 最后点击保存即可。

5 计算最短路径

使用如下url即可获取计算得到的最短路径的geojson格式，其中，参数如下表所示：

参数示例含义testgeoserver发布的postgis图层的工作区typeName工具区:视图名viewparamsx1、y1、x2、y2是最短路径分析的起点和终点，用“ ; ”隔开

http://localhost:8080/geoserver/test/wfs?service=WFS&version=1.1.0&request=GetFeature&typeName=test:new&outputformat=JSON&viewparams=x1:111.882740557158;y1:29.8983406704794;x2:113.047556556274;y2:29.0935587074538

6 问题记录

使用postgis shapefile import/export manager工具上传shapefile至数据库时，忘记设定SRID为4326。

网上很多帖子的计算最短路径的函数都存在问题： （1）本实验使用的pgrouting为3.0版本，其中pgr_dijkstra函数的参数都为4个，具体函数参数可以在postgres中查看。但是网上贴子中使用的pgr_dijkstra函数的参数都为5个，因此需要进行修改。可参考如下网址中对该函数的介绍，需要注意返回的结果也发生了变化，不再是id、id2等字段表示节点，而是node字段标识节点，edge字段表示边id。

http://www.mamicode.com/info-detail-2781808.html

（2) 使用时候发现st_linelocatepoint函数参数有时候会报错，报错结果如下图所示，查看使用该函数的传入参数，可以发现传入参数为multilinestring类型，而不是linestring类型。 继续检查可以发现前面在执行pgr_dijkstra函数的时候，就已经返回了multilinestring，这是因为sql语句里面写了后面的group和order导致在使用st_linemerge并未合并线段，因此需要将所有的group和order删除。

-- 改写前 execute 'SELECT st_linemerge(st_union(b.geom)) ' || 'FROM pgr_dijkstra( ''SELECT gid as id, source, target, length as cost FROM ' ||tbl|| '''::text,' ||'array['||v_startSource||'] , ' ||'array['||v_endTarget||'] , false ) a, ' ||tbl|| ' b WHERE a.edge=b.gid GROUP by edge ORDER by edge' into v_res_c ; -- 改写后 execute 'SELECT st_linemerge(st_union(b.geom)) ' || 'FROM pgr_dijkstra( ''SELECT gid as id, source, target, length as cost FROM ' ||tbl|| '''::text,' ||'array['||v_startSource||'] , ' ||'array['||v_endTarget||'] , false ) a, ' ||tbl|| ' b WHERE a.edge=b.gid' into v_res_c ;

在使用pgr_dijkstra函数计算完最短路径后，还需要在路径上截取起点至终点需要的部分，再将截取的结果与起点、终点拼接起来，得到最后最短路径的结果。必须要截取和拼接，否则返回结果是错误的。

pgsql代码难以调试的时候，可以使用如下方法进行变量的输出，再某行代码查看输出结果，以判断错误。

raise notice '%', 变量名

在geoserver发布新视图后，点击从数据计算边框的时候，经常会报错。这时候需要拉到页面底部查看错误，此时错误多维pgsql中的错误，可以根据该错误调整函数代码，重新定义函数。