前言
前段時(shí)間遇到了實(shí)際的需求，在特定的路網(wǎng)中查詢最短路徑邀层。同時(shí)配合 Cesium 進(jìn)行動(dòng)態(tài)顯示棕洋。
需求
- 動(dòng)態(tài)查詢兩點(diǎn)之間的最短路徑（起點(diǎn)固定）；
- 查詢的路徑高亮顯示逗扒；
- Cesium 對(duì)生成的路徑進(jìn)行小車移動(dòng)展示古戴。
技術(shù)實(shí)現(xiàn)路線
- 動(dòng)態(tài)查詢兩點(diǎn)之間的最短路徑 -> Postgresql 中的 PgRouting 實(shí)現(xiàn)；
- 查詢的路徑高亮顯示 -> Cesium 中
PolylineGlowMaterialProperty進(jìn)行高亮顯示矩肩； - 路徑進(jìn)行小車移動(dòng)展示 ->CZML 進(jìn)行動(dòng)態(tài)展示
實(shí)現(xiàn)效果
實(shí)現(xiàn)效果
Postgresql 最短路徑實(shí)現(xiàn)
起初是參考 PgRouting 官網(wǎng) 的做法现恼。但是這種做法是對(duì)數(shù)據(jù)進(jìn)行拓?fù)洌捎邢驁D（或者無(wú)向圖）采用 dijkstra 算法進(jìn)行最短路徑的生成。這種方法最大的問(wèn)題就是判斷鼠標(biāo)點(diǎn)擊的點(diǎn)位于有向圖的位置述暂。相對(duì)來(lái)說(shuō)比較麻煩痹升。
官網(wǎng)說(shuō)明
實(shí)現(xiàn)流程：
- 首先將數(shù)據(jù)導(dǎo)入 Postgresql 數(shù)據(jù)庫(kù)
- 對(duì)數(shù)據(jù)進(jìn)行路網(wǎng)拓?fù)鋽?shù)據(jù)計(jì)算處理，執(zhí)行成功后畦韭，執(zhí)行成功后會(huì)生產(chǎn)一個(gè) vertices_pgr 的表疼蛾，里面包含路網(wǎng)相交點(diǎn)的空間數(shù)據(jù)
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 語(yǔ)句前端傳入的是有向圖中的某兩個(gè)端點(diǎn)。
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
最方便的還是直接傳入起始點(diǎn)的坐標(biāo)進(jìn)行路徑的查詢
感謝 itas109 提供的經(jīng)緯度查詢的最短路徑的方法艺配。這是 Github 地址 https://github.com/itas109/postgis_navigation察郁。
思路是將創(chuàng)建新的函數(shù)，將鼠標(biāo)點(diǎn)擊的兩點(diǎn)經(jīng)緯度傳入獲得最短路徑的返回值转唉。
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ā)點(diǎn)最近的道路節(jié)點(diǎn)
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é)點(diǎn)
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ù)之前要對(duì)數(shù)據(jù)進(jìn)行處理皮钠。路網(wǎng)數(shù)據(jù)必須在節(jié)點(diǎn)處打斷，同時(shí)在 ArcMap 中進(jìn)行拓?fù)浼m錯(cuò)赠法。數(shù)據(jù)導(dǎo)入后進(jìn)行如下處理：
ALTER TABLE road ADD COLUMN source integer; --添加起點(diǎn)id字段
ALTER TABLE road ADD COLUMN target integer; --添加終點(diǎn)id字段
ALTER TABLE road ADD COLUMN length double precision; --添加道路權(quán)重字段
SELECT pgr_createTopology('road',0.00001, 'geom', 'gid'); --為source和target賦值麦轰，并創(chuàng)建拓?fù)潼c(diǎn)表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)建起點(diǎn)經(jīng)度x1
ALTER TABLE road ADD COLUMN y1 double precision; --創(chuàng)建起點(diǎn)緯度y1
ALTER TABLE road ADD COLUMN x2 double precision; --創(chuàng)建起點(diǎn)經(jīng)度x2
ALTER TABLE road ADD COLUMN y2 double precision; --創(chuàng)建起點(diǎn)經(jī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í)行查詢語(yǔ)句:
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
查詢結(jié)果：
查詢結(jié)果
接下來(lái)就是對(duì)查詢的 GeoJSON 數(shù)據(jù)轉(zhuǎn)換為 CZML 數(shù)據(jù)在三維場(chǎng)景中進(jìn)行展示了
