Postgresql构建经纬度查询两点之间的最短路径
前言
前段时间遇到了实际的需求，在特定的路网中查询最短路径。同时配合 Cesium 进行动态显示。
需求
- 动态查询两点之间的最短路径（起点固定）；
- 查询的路径高亮显示；
- Cesium 对生成的路径进行小车移动展示。
技术实现路线
- 动态查询两点之间的最短路径 -> Postgresql 中的 PgRouting 实现；
- 查询的路径高亮显示 -> Cesium 中
PolylineGlowMaterialProperty进行高亮显示； - 路径进行小车移动展示 ->CZML 进行动态展示
实现效果
Postgresql 最短路径实现
起初是参考 PgRouting 官网 的做法。但是这种做法是对数据进行拓扑，生成有向图（或者无向图）采用 dijkstra 算法进行最短路径的生成。这种方法最大的问题就是判断鼠标点击的点位于有向图的位置。相对来说比较麻烦。
实现流程：
- 首先将数据导入 Postgresql 数据库
- 对数据进行路网拓扑数据计算处理，执行成功后，执行成功后会生产一个 vertices_pgr 的表，里面包含路网相交点的空间数据
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。
思路是将创建新的函数，将鼠标点击的两点经纬度传入获得最短路径的返回值。
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
-- 查询距离出发点最近的道路节点
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;
-- 查询距离目的地最近的道路节点
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
使用新函数之前要对数据进行处理。路网数据必须在节点处打断，同时在 ArcMap 中进行拓扑纠错。数据导入后进行如下处理：
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赋值，并创建拓扑点表road_vertices_pgr
update road set length =st_length(geom); --为length赋值
CREATE INDEX source_idx ON road("source"); --为source字段创建索引
CREATE INDEX target_idx ON road("target"); --为target字段创建索引
ALTER TABLE road ADD COLUMN x1 double precision; --创建起点经度x1
ALTER TABLE road ADD COLUMN y1 double precision; --创建起点纬度y1
ALTER TABLE road ADD COLUMN x2 double precision; --创建起点经度x2
ALTER TABLE road ADD COLUMN y2 double precision; --创建起点经度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)));
然后执行查询语句:
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 数据转换为 CZML 数据在三维场景中进行展示了
