PostGIS 結合Openlayers以及Geoserver實現最短路徑分析（二）

前文講述了怎么用ArcMap制作了測試數據，並導入了PostGIS，接下來我們需要結合PgRouting插件，對入庫的數據再進行一下處理。

0、引入擴展包

postgis

pgrouting

postgis_topology

fuzzystrmatch

1、在pgAdmin中，執行下面的sql語句

 --添加起點字段source
ALTER TABLE zy ADD COLUMN source integer;

--添加終點字段target
ALTER table zy add column target integer;

--添加道路權重值字段
ALTER TABLE zy ADD COLUMN length double precision;

--建立拓撲關系，這里會為source和target字段賦值（表明，閾值，要素字段名，主鍵id）
SELECT pgr_createTopology('zy', 0.00001, 'geom', 'gid');
###這里在創建完拓撲，source和target賦值了以后，可以給上一篇文章中創建的點圖層預留的Id字段賦值了（參看第5步）

--為source和target字段創建索引（加快查詢速度）
CREATE INDEX source_idx ON zy("source");
CREATE INDEX target_idx ON zy("target");

--添加線段端點坐標
ALTER TABLE ZY ADD COLUMN x1 double precision;        --創建起點經度x1
ALTER TABLE ZY ADD COLUMN y1 double precision;        --創建起點緯度y1
ALTER TABLE ZY ADD COLUMN x2 double precision;        --創建起點經度x2
ALTER TABLE ZY ADD COLUMN y2 double precision;        --創建起點經度y2
--給x1、y1、x2、y2賦值
UPDATE ZY SET x1 =ST_x(ST_PointN(geom, 1));    
UPDATE ZY SET y1 =ST_y(ST_PointN(geom, 1));    
UPDATE ZY SET x2 =ST_x(ST_PointN(geom, ST_NumPoints(geom)));    
UPDATE ZY SET y2 =ST_y(ST_PointN(geom, ST_NumPoints(geom)));     

--為length賦值（以長度作為權重）
UPDATE zy SET length =st_length(geom);

--將長度值賦給reverse_cost，作為路線選擇標准（消耗）
ALTER TABLE zy ADD COLUMN reverse_cost double precision;
UPDATE zy SET reverse_cost = st_length(geom);

2、測試數據，如果執行下面的語句有結果，則表示數據預處理成功

--通過起點號、終點號查詢最短路徑
--source為線表起點字段名稱
--target為線表終點字段名稱
--起點終點前后順序無固定要求
--length為長度字段，也可以使用自己的評價體系
--1、9為測試使用起點號\終點號
--zy表名
--id1經過節點號
--id2經過路網線的gid
SELECT seq, id1 AS node, id2 AS edge, cost FROM pgr_dijkstra('
                SELECT  gid AS id,
                         source::integer,
                         target::integer,
                         length::double precision AS cost
                        FROM zy',
                1, 4, false, false);

3、創建最短路徑查詢函數（核心）

-- FUNCTION: public.pgr_shortestpath(character varying, double precision, double precision, double precision, double precision)

-- DROP FUNCTION public.pgr_shortestpath(character varying, double precision, double precision, double precision, double precision);

CREATE OR REPLACE FUNCTION public.pgr_shortestpath(
tbl character varying,
startx double precision,
starty double precision,
endx double precision,
endy double precision)
RETURNS geometry
LANGUAGE 'plpgsql'
COST 100
VOLATILE STRICT
AS $BODY$

declare
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
--查詢離起點最近的線
--3857坐標系
--找起點15米范圍內的最近線
execute 'select geom, source, target from ' ||tbl||
' where ST_DWithin(geom,ST_Geometryfromtext(''point('||startx ||' ' || starty||')'',3857),15)
order by ST_Distance(geom,ST_GeometryFromText(''point('|| startx ||' '|| starty ||')'',3857)) limit 1'
into v_startLine, v_startSource ,v_startTarget;
--查詢離終點最近的線
--找終點15米范圍內的最近線
execute 'select geom, source, target from ' ||tbl||
' where ST_DWithin(geom,ST_Geometryfromtext(''point('|| endx || ' ' || endy ||')'',3857),15)
order by ST_Distance(geom,ST_GeometryFromText(''point('|| endx ||' ' || endy ||')'',3857)) 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 ;

raise notice '%', v_startLine;
raise notice '%', v_endLine;

select ST_ClosestPoint(v_startLine, ST_Geometryfromtext('point('|| startx ||' ' || starty ||')',3857)) into v_statpoint;
select ST_ClosestPoint(v_endLine, ST_GeometryFromText('point('|| endx ||' ' || endy ||')',3857)) into v_endpoint;
-- ST_Distance
--從開始的起點到結束的起點最短路徑
execute 'SELECT st_linemerge(st_union(b.geom)) ' ||
'FROM pgr_kdijkstraPath(
''SELECT gid as id, source, target, length as cost FROM ' || tbl ||''','
||v_startSource|| ', ' ||'array['||v_endSource||'] , false, false
) a, '
||tbl|| ' b
WHERE a.id3=b.gid
GROUP by id1
ORDER by id1' into v_res ;
--從開始的終點到結束的起點最短路徑
execute 'SELECT st_linemerge(st_union(b.geom)) ' ||
'FROM pgr_kdijkstraPath(
''SELECT gid as id, source, target, length as cost FROM ' || tbl ||''','
||v_startTarget|| ', ' ||'array['||v_endSource||'] , false, false
) a, '
||tbl|| ' b
WHERE a.id3=b.gid
GROUP by id1
ORDER by id1' into v_res_b ;
--從開始的起點到結束的終點最短路徑
execute 'SELECT st_linemerge(st_union(b.geom)) ' ||
'FROM pgr_kdijkstraPath(
''SELECT gid as id, source, target, length as cost FROM ' || tbl ||''','
||v_startSource || ', ' ||'array['||v_endTarget||'] , false, false
) a, '
|| tbl || ' b
WHERE a.id3=b.gid
GROUP by id1
ORDER by id1' into v_res_c ;
--從開始的終點到結束的終點最短路徑
execute 'SELECT st_linemerge(st_union(b.geom)) ' ||
'FROM pgr_kdijkstraPath(
''SELECT gid as id, source, target, length as cost FROM ' || tbl ||''','
||v_startTarget || ', ' ||'array['||v_endTarget||'] , false, false
) a, '
|| tbl || ' b
WHERE a.id3=b.gid
GROUP by id1
ORDER by id1' 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進行拼接
raise notice '%', v_res;
select ST_LineMerge(ST_Union(array[v_res,v_startLine,v_endLine])) into v_res;
raise notice '%', v_res;
select ST_Line_Locate_Point(v_res, v_statpoint) into v_perStart;
select ST_Line_Locate_Point(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_Line_SubString(v_res,v_perStart, v_perEnd) into v_shPath;
raise notice '%', v_shPath;
return v_shPath;
end;
$BODY$;

ALTER FUNCTION public.pgr_shortestpath(character varying, double precision, double precision, double precision, double precision)
OWNER TO postgres;

 4、測試

在調用函數做測試的時候，遇到ST_Line_Locate_Point時報錯：“line_locate_point : 1st arg isnt a line”，查資料解釋說是函數第一個參數格式應該是LineString，我檢查了一下我的v_res格式，是MultiLineString所以報錯。（參考了https://www.jianshu.com/p/4b9d22406bce，在此謝謝作者大大，很詳細）

 

知道原因以后怎么解決呢？還是個頭疼的問題。我把查詢出來的MultiLineString放到ArcMap中，圖形化的顯示出來，發現這兩段竟然是不連通的兩端，最后分析應該是數據不連通導致的，兩點間查不到結果，沒有最短路徑。

在ArcMap中利用網絡分析工具(networkAnalyse tool)試了一下也不行，果然是這原因，后來換了兩個連通的點查詢結果就沒問題，所以函數沒有問題，結果控制一下就行。 

5、至此數據預處理結束了，接下來在下一篇中介紹Geoserver中發布數據並調用

6、后續補充

在代碼中提到做ST_Distance時，分別查詢了四次查詢：

--從開始的起點到結束的起點最短路徑 

--從開始的終點到結束的起點最短路徑 

--從開始的起點到結束的終點最短路徑 

--從開始的終點到結束的終點最短路徑

原來在做的時候不是很理解為什么要做這四次查詢，是否多余。后來在做帶路徑的最短路徑查詢時，有點心得。

 

實際查詢的是A到C，A到D，B到C，B到D

如果不做這步的話，取任意點做路徑，萬一取到BC或者BD，其實就變成同一個線段了，沒有意義了