Postgresql構建經緯度查詢兩點之間的最短路徑

前言

前段時間遇到了實際的需求，在特定的路網中查詢最短路徑。同時配合 Cesium 進行動態(tài)顯示。

需求

  • 動態(tài)查詢兩點之間的最短路徑（起點固定）；
  • 查詢的路徑高亮顯示；
  • Cesium 對生成的路徑進行小車移動展示。

技術實現(xiàn)路線

  1. 動態(tài)查詢兩點之間的最短路徑 -> Postgresql 中的 PgRouting 實現(xiàn)；
  2. 查詢的路徑高亮顯示 -> Cesium 中 PolylineGlowMaterialProperty 進行高亮顯示；
  3. 路徑進行小車移動展示 ->CZML 進行動態(tài)展示

實現(xiàn)效果

實現(xiàn)效果

Postgresql 最短路徑實現(xiàn)

起初是參考 PgRouting 官網 的做法。但是這種做法是對數(shù)據進行拓撲，生成有向圖（或者無向圖）采用 dijkstra 算法進行最短路徑的生成。這種方法最大的問題就是判斷鼠標點擊的點位于有向圖的位置。相對來說比較麻煩。

官網說明

實現(xiàn)流程：

  • 首先將數(shù)據導入 Postgresql 數(shù)據庫
  • 對數(shù)據進行路網拓撲數(shù)據計算處理，執(zhí)行成功后，執(zhí)行成功后會生產一個 vertices_pgr 的表，里面包含路網相交點的空間數(shù)據
alter table road add column source int;
alter table road add column target int;

create index road_source_idx on road("source");
create index road_target_idx on road("target");

ALTER TABLE road  ADD COLUMN length double precision;  

SELECT pgr_createTopology('road',0.00001, 'geom', 'gid');  

update road set length =st_length(geom);
  • 查詢最短路徑 sql 語句前端傳入的是有向圖中的某兩個端點。
SELECT seq, id1 AS node, id2 AS edge, cost,geom  FROM pgr_dijkstra('  
SELECT gid AS id,              
source::integer,                 
target::integer,                
length::double precision AS cost  
FROM xmpark_road',  
1, 10, false, false) as di  
join xmpark_road pt  
on di.id2 = pt.gid

最方便的還是直接傳入起始點的坐標進行路徑的查詢

感謝 itas109 提供的經緯度查詢的最短路徑的方法。這是 Github 地址 https://github.com/itas109/postgis_navigation。

思路是將創(chuàng)建新的函數(shù)，將鼠標點擊的兩點經緯度傳入獲得最短路徑的返回值。

CREATE OR REPLACE FUNCTION "public"."pgr_fromatob"(IN "tbl" varchar, IN "x1" float8, IN "y1" float8, IN "x2" float8, IN "y2" float8, OUT "seq" int4, OUT "gid" int4, OUT "name" text, OUT "heading" float8, OUT "cost" float8, OUT "geom" "public"."geometry")
  RETURNS SETOF "pg_catalog"."record" AS $BODY$  
DECLARE  
        sql     text;  
        rec     record;  
        source  integer;  
        target  integer;  
        point   integer;  

BEGIN  
    -- 查詢距離出發(fā)點最近的道路節(jié)點  
    EXECUTE 'SELECT id::integer FROM road_vertices_pgr   
            ORDER BY the_geom <-> ST_GeometryFromText(''POINT('   
            || x1 || ' ' || y1 || ')'',900913) LIMIT 1' INTO rec;  
    source := rec.id;  

    -- 查詢距離目的地最近的道路節(jié)點  
    EXECUTE 'SELECT id::integer FROM road_vertices_pgr   
            ORDER BY the_geom <-> ST_GeometryFromText(''POINT('   
            || x2 || ' ' || y2 || ')'',900913) LIMIT 1' INTO rec;  
    target := rec.id;  

    -- 最短路徑查詢   
        seq := 0;  
        sql := 'SELECT gid, geom, cost, source, target,   
                ST_Reverse(geom) AS flip_geom FROM ' ||  
                        'pgr_bdAstar(''SELECT gid as id, source::int, target::int, '  
                                        || 'length::float AS cost,x1,y1,x2,y2 FROM '  
                                        || quote_ident(tbl) || ''', '  
                                        || source || ', ' || target   
                                        || ' ,false, false), '  
                                || quote_ident(tbl) || ' WHERE id2 = gid ORDER BY seq';  

    -- Remember start point  
        point := source;  

        FOR rec IN EXECUTE sql  
        LOOP  
        -- Flip geometry (if required)  
        IF ( point != rec.source ) THEN  
            rec.geom := rec.flip_geom;  
            point := rec.source;  
        ELSE  
            point := rec.target;  
        END IF;  

        -- Calculate heading (simplified)  
        EXECUTE 'SELECT degrees( ST_Azimuth(   
                ST_StartPoint(''' || rec.geom::text || '''),  
                ST_EndPoint(''' || rec.geom::text || ''') ) )'   
            INTO heading;  

        -- Return record  
                seq     := seq + 1;  
                gid     := rec.gid;   
                cost    := rec.cost;  
                geom    := rec.geom;  
                RETURN NEXT;  
        END LOOP;  
        RETURN;  
END;  
$BODY$
  LANGUAGE plpgsql VOLATILE STRICT
  COST 100
  ROWS 1000

使用新函數(shù)之前要對數(shù)據進行處理。路網數(shù)據必須在節(jié)點處打斷，同時在 ArcMap 中進行拓撲糾錯。數(shù)據導入后進行如下處理：

ALTER TABLE road ADD COLUMN source integer;     --添加起點id字段
ALTER TABLE road ADD COLUMN target integer;     --添加終點id字段
ALTER TABLE road ADD COLUMN length double precision;    --添加道路權重字段

SELECT pgr_createTopology('road',0.00001, 'geom', 'gid');   --為source和target賦值，并創(chuàng)建拓撲點表road_vertices_pgr

update road set length =st_length(geom);            --為length賦值

CREATE INDEX source_idx ON road("source");      --為source字段創(chuàng)建索引
CREATE INDEX target_idx ON road("target");          --為target字段創(chuàng)建索引

ALTER TABLE road ADD COLUMN x1 double precision;        --創(chuàng)建起點經度x1
ALTER TABLE road ADD COLUMN y1 double precision;        --創(chuàng)建起點緯度y1
ALTER TABLE road ADD COLUMN x2 double precision;        --創(chuàng)建起點經度x2
ALTER TABLE road ADD COLUMN y2 double precision;        --創(chuàng)建起點經度y2

UPDATE road SET x1 =ST_x(ST_PointN(geom, 1));
UPDATE road SET y1 =ST_y(ST_PointN(geom, 1));
UPDATE road SET x2 =ST_x(ST_PointN(geom, ST_NumPoints(geom)));
UPDATE road SET y2 =ST_y(ST_PointN(geom, ST_NumPoints(geom)));

然后執(zhí)行查詢語句:

SELECT  st_asgeojson(st_makeline(route.geom)) FROM (SELECT geom FROM pgr_fromAtoB('road', 118.574693042441, 31.0002595461945,118.575197797, 31.0068716390001)ORDER BY seq) AS route

查詢結果：

查詢結果

接下來就是對查詢的 GeoJSON 數(shù)據轉換為 CZML 數(shù)據在三維場景中進行展示了

?著作權歸作者所有,轉載或內容合作請聯(lián)系作者
【社區(qū)內容提示】社區(qū)部分內容疑似由AI輔助生成，瀏覽時請結合常識與多方信息審慎甄別。
平臺聲明：文章內容（如有圖片或視頻亦包括在內）由作者上傳并發(fā)布，文章內容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務。

友情鏈接更多精彩內容