基於pgrouting的路徑規劃處理方法

2022-04-18 16:01:10

對於GIS業務來說，路徑規劃是非常基礎的一個業務，一般公司如果處理，都會直接選擇呼叫已經成熟的第三方的介面，比如高德、百度等。當然其實路徑規劃的演演算法非常多，像比較著名的Dijkstra、A*演演算法等。當然本篇文章不是介紹演演算法的，本文作者會根據pgrouting已經整合的Dijkstra演演算法來，結合postgresql資料庫來處理最短路徑。

一、資料處理

路徑規劃的核心是資料，資料是一般的路網資料，但是我們拿到路網資料之後，需要對資料進行處理，由於演演算法的思想是基於有向圖的原理，因此首先需要對資料做topo處理，通過topo我們其實就建立了路網中各條道路的頂點關係，下面是主要命令：

--開啟執行路網topo的外掛
create extension postgis;
create extension postgis_topology;
--資料建立拓撲
ALTER TABLE test_road ADD COLUMN source integer;
ALTER TABLE test_road ADD COLUMN target integer;
SELECT pgr_createTopology('test_road',0.00001, 'geom', 'gid');

其中test_road是將路網資料匯入到postgresql中的表名。

處理完topo之後，基本就夠用了，我們就可以藉助pgrouting自帶的函數，其實有很多，我們選擇pgr_dijkstra

CREATE OR REPLACE FUNCTION public.pgr_dijkstra(
    IN edges_sql text,
    IN start_vid bigint,
    IN end_vid bigint,
    IN directed boolean,
    OUT seq integer,
    OUT path_seq integer,
    OUT node bigint,
    OUT edge bigint,
    OUT cost double precision,
    OUT agg_cost double precision)
  RETURNS SETOF record AS
$BODY$
DECLARE
BEGIN
    RETURN query SELECT *
    FROM _pgr_dijkstra(_pgr_get_statement($1), start_vid, end_vid, directed, false);
  END
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100
  ROWS 1000;
ALTER FUNCTION public.pgr_dijkstra(text, bigint, bigint, boolean)
  OWNER TO postgres;

從函數輸入引數可以看到，我們需要一個查詢sql,一個起始點、一個結束點、以及是否考慮方向，好了瞭解到呼叫函數輸入引數，我們就來寫這個函數。

二、原理分析

一般路徑規劃，基本都是輸入一個起點位置、一個終點位置然後直接規劃，那麼對於我們來說，要想套用上面的函數，必須找出起點位置target ，以及終點位置的source，然後規劃根據找出的這兩個topo點，呼叫上面的函數，來返回自己所需要的結果。

如何根據起始點找到對應的target呢，其實就是找離起點最近線的target，同理終點的source，其實就是找離終點最近線的source,當然將這兩個點規劃規劃好之後，基本就可以了，但是最後還需要將起點到起點最近先的target連線起來，終點到終點最近線的source連線起來，這樣整個路徑規劃就算完成了。

下面我們來看具體的實現儲存過程：

CREATE OR REPLACE FUNCTION public.pgr_shortest_road(
IN startx double precision,
IN starty double precision,
IN endx double precision,
IN endy double precision,
OUT road_name character varying,
OUT v_shpath character varying,
OUT cost double precision)
RETURNS SETOF record AS
$BODY$ 
declare 
v_startLine geometry;--離起點最近的線 
v_endLine geometry;--離終點最近的線 
v_startTarget integer;--距離起點最近線的終點 
v_endSource integer;--距離終點最近線的起點 
v_statpoint geometry;--在v_startLine上距離起點最近的點 
v_endpoint geometry;--在v_endLine上距離終點最近的點 
v_res geometry;--最短路徑分析結果 
v_perStart float;--v_statpoint在v_res上的百分比 
v_perEnd float;--v_endpoint在v_res上的百分比 
v_rec record; 
first_name varchar;
end_name varchar;
first_cost double precision;
end_cost double precision;
begin 
--查詢離起點最近的線 
execute 'select geom,target,name from china_road where 
ST_DWithin(geom,ST_Geometryfromtext(''point('|| startx ||' ' || starty||')''),0.01) 
order by ST_Distance(geom,ST_GeometryFromText(''point('|| startx ||' '|| starty ||')'')) limit 1' 
into v_startLine ,v_startTarget,first_name; 
--查詢離終點最近的線 
execute 'select geom,source,name from china_road
where ST_DWithin(geom,ST_Geometryfromtext(''point('|| endx || ' ' || endy ||')''),0.01) 
order by ST_Distance(geom,ST_GeometryFromText(''point('|| endx ||' ' || endy ||')'')) limit 1' 
into v_endLine,v_endSource,end_name; 
--如果沒找到最近的線，就返回null 
if (v_startLine is null) or (v_endLine is null) then 
return; 
end if ; 
select ST_ClosestPoint(v_startLine, ST_Geometryfromtext('point('|| startx ||' ' || starty ||')')) into v_statpoint; 
select ST_ClosestPoint(v_endLine, ST_GeometryFromText('point('|| endx ||' ' || endy ||')')) into v_endpoint; 

--計算距離起點最近線上的點在該線中的位置
select ST_Line_Locate_Point(st_linemerge(v_startLine), v_statpoint) into v_perStart;

select ST_Line_Locate_Point(st_linemerge(v_endLine), v_endpoint) into v_perEnd;

select ST_Distance_Sphere(v_statpoint,ST_PointN(ST_GeometryN(v_startLine,1), ST_NumPoints(ST_GeometryN(v_startLine,1)))) into first_cost;

select ST_Distance_Sphere(ST_PointN(ST_GeometryN(v_endLine,1),1),v_endpoint) into end_cost; 

if (ST_Intersects(st_geomfromtext('point('|| startx ||' '|| starty ||') '), v_startLine) and ST_Intersects(st_geomfromtext('point('|| endx ||' '|| endy ||') '), v_startLine)) then 
select ST_Distance_Sphere(v_statpoint, v_endpoint) into first_cost;

select ST_Line_Locate_Point(st_linemerge(v_startLine), v_endpoint) into v_perEnd;
for v_rec in 
select ST_Line_SubString(st_linemerge(v_startLine), v_perStart,v_perEnd) as point,COALESCE(end_name,'無名路') as name,end_cost as cost loop
v_shPath:= ST_AsGeoJSON(v_rec.point);
cost:= v_rec.cost;
road_name:= v_rec.name;
return next;
end loop;
return;
end if;
--最短路徑 
for v_rec in 
(select ST_Line_SubString(st_linemerge(v_startLine),v_perStart,1) as point,COALESCE(first_name,'無名路') as name,first_cost as cost
union all
SELECT st_linemerge(b.geom) as point,COALESCE(b.name,'無名路') as name,b.length as cost
FROM pgr_dijkstra(
'SELECT gid as id, source, target, length as cost FROM china_road
where st_intersects(geom,st_buffer(st_linefromtext(''linestring('||startx||' ' || starty ||','|| endx ||' ' || endy ||')''),0.05))', 
v_startTarget, v_endSource , false 
) a, china_road b 
WHERE a.edge = b.gid
union all
select ST_Line_SubString(st_linemerge(v_endLine),0,v_perEnd) as point,COALESCE(end_name,'無名路') as name,end_cost as cost)
loop
v_shPath:= ST_AsGeoJSON(v_rec.point);
cost:= v_rec.cost;
road_name:= v_rec.name;
return next;
end loop; 
end; 
$BODY$
LANGUAGE plpgsql VOLATILE STRICT;

上面這種實現，是將所有查詢道路返回一個集合，然後使用者端來將各個線路進行合併，這種方式對最終效率影響比較大，所以一般會在函數中將道路何合併為一條道路，我們可以使用postgis的st_union函數來處理，小編經過長時間的試驗，在保證效率和準確性的情況下，對上面的儲存過程做了很多優化，最終得出瞭如下：

CREATE OR REPLACE FUNCTION public.pgr_shortest_road(
    startx double precision,
    starty double precision,
    endx double precision,
    endy double precision)
  RETURNS geometry AS
$BODY$   
declare  
    v_startLine geometry;--離起點最近的線  
    v_endLine geometry;--離終點最近的線 
    v_perStart float;--v_statpoint在v_res上的百分比  
    v_perEnd float;--v_endpoint在v_res上的百分比  
    v_shpath geometry;
    distance double precision;
    bufferInstance double precision;
    bufferArray double precision[];  
begin  
    execute 'select geom,
        case china_road.direction
        when ''3'' then
        source
        else
        target
        end 
        from china_road where  
        ST_DWithin(geom,ST_Geometryfromtext(''point('|| startx ||' ' || starty||')'',4326),0.05) 
        AND width::double precision >= '||roadWidth||'
        order by ST_Distance(geom,ST_GeometryFromText(''point('|| startx ||' '|| starty ||')'',4326))  limit 1'  
    into v_startLine; 
     
    execute 'select geom,
        case china_road.direction
        when ''3'' then
        target
        else
        source
        end 
        from china_road
        where ST_DWithin(geom,ST_Geometryfromtext(''point('|| endx || ' ' || endy ||')'',4326),0.05)  
        AND width::double precision >= '||roadWidth||'
        order by ST_Distance(geom,ST_GeometryFromText(''point('|| endx ||' ' || endy ||')'',4326))  limit 1'  
    into v_endLine;  
    
    if (v_startLine is null) or (v_endLine is null) then  
        return null;
    end if; 
         
    if (ST_equals(v_startLine,v_endLine)) then 
        select ST_LineLocatePoint(st_linemerge(v_startLine), ST_Geometryfromtext('point('|| startx ||' ' || starty ||')',4326)) into v_perStart;
        select ST_LineLocatePoint(st_linemerge(v_endLine), ST_Geometryfromtext('point('|| endx ||' ' || endy ||')',4326)) into v_perEnd;
        select ST_LineSubstring(st_linemerge(v_startLine),v_perStart,v_perEnd) into v_shPath;
        return v_shPath;
    end if;
    
    select ST_DistanceSphere(st_geomfromtext('point('|| startx ||' ' || starty ||')',4326),st_geomfromtext('point('|| endx ||' ' || endy ||')',4326)) into distance;
    if ((distance / 1000) > 50) then
        bufferArray := ARRAY[0.1,0.2,0.3,0.5,0.8];    
    else
        bufferArray := ARRAY[0.02,0.05,0.08,0.1];
    end if;
    
    forEACH bufferInstance IN ARRAY bufferArray
    LOOP
    select _pgr_shortest_road(startx,starty,endx,endy,bufferInstance) into v_shPath;
    if (v_shPath is not null) then    
        return v_shPath;
    end if;
    end loop;
    
    end;  
    $BODY$
  LANGUAGE plpgsql VOLATILE STRICT
  COST 100;
ALTER FUNCTION public.pgr_shortest_road(double precision, double precision, double precision, double precision )
  OWNER TO postgres;

DROP FUNCTION public._pgr_shortest_road(double precision, double precision, double precision, double precision, double precision);

上面的函數，其實對於大部分情況下的操作，基本可以滿足了。

三、效率優化

其實在資料查詢方面，我們使用的是起點和終點之間的線性緩衝來提高效率，如下：

SELECT gid as id, source, target, cost,rev_cost as reverse_cost FROM china_road
      where geom && st_buffer(st_linefromtext(''linestring('||startx||' ' || starty ||','|| endx ||' ' || endy ||')'',4326),'||bufferDistance||')

當然這在大部分情況下，依舊是不錯的，然後在有些情況下，並不能起到很好的作用，因為如果起點和終點之間道路偏移較大（比如直線上的山脈較多的時候，路就會比較繞），這個時候，可能會增大緩衝距離，而增加緩衝距離就會導致，部分割區域的查詢量增大，繼而影響效率，因此其實我們可以考慮使用mapid這個引數，這個引數從哪來呢，一般我們拿到的路網資料都會這個欄位，我們只需要生成一個區域表，而這個區域表就倆個欄位，一個是mapid，一個是這個mapid的polygon範圍，這樣子，上面的查詢條件，就可以換成如下：　　

SELECT gid as id, source, target, cost,rev_cost as reverse_cost FROM china_road
      where mapid in (select mapid from maps where geom && st_buffer(st_linefromtext(''linestring('||startx||' ' || starty ||','|| endx ||' ' || endy ||')''),'||bufferDistance||'))

這樣就可以在很大程度上提高效率。

四、資料bug處理

其實有時候我們拿到的路網資料，並不是非常的準確，或者說是錄入的有瑕疵，我自己遇到的就是生成的topo資料，本來一條路的target應該和它相鄰路的source的點重合，然後實際卻是不一樣，這就導致最終規劃處的有問題，因此，簡單寫了一個處理這種問題的函數　

CREATE OR REPLACE FUNCTION public.modity_road_data()
RETURNS void AS
$BODY$   
declare  
n integer;
begin  
    for n IN (select distinct(source) from china_road )  loop
        update china_road
        set geom = st_multi(st_addpoint(ST_geometryN(geom,1),
        (select st_pointn(ST_geometryN(geom,1),1) from china_road where source = n
        limit 1),
        st_numpoints(ST_geometryN(geom,1))))
        where target = n;
    end loop;
    end;  
    $BODY$
  LANGUAGE plpgsql VOLATILE STRICT
  COST 100;
ALTER FUNCTION public.modity_road_data()
OWNER TO postgres;

五、後續規劃

上面的函數已在百萬資料中做過驗證，後續還會驗證千萬級別的路網資料，當然這種級別，肯定要在策略上做一些調整了，比如最近測試的全國路網中，先規劃起點至起點最近的高速入口，在規劃終點至終點最近的高速出口，然後再高速路網上規劃高速入口到高速出口的路徑，這樣發現效率提升不少，當然，這裡面還有很多邏輯和業務，等所有東西都驗證完畢，會再出一版，千萬級別路徑規劃的文章。

到此這篇關於基於pgrouting的路徑規劃處理的文章就介紹到這了,更多相關pgrouting的路徑規劃內容請搜尋it145.com以前的文章或繼續瀏覽下面的相關文章希望大家以後多多支援it145.com！