GeoServer+PostgreSQL+PostGIS+pgRouting实现最短路径查询
1 数据准备
1.1 从OSM下载区域SHP数据
https://www.openstreetmap.org/export#map=11/30.6731/104.2019
手动画出需要的范围，导出
2 数据检查处理
2.1 编辑道路
将shp文件导入ArcMap，启用编辑
根据天地图矢量和影像图层，调整和补充道路
2.2 打断线段
全选所有道路，用Topology工具中的Planarize Lines打断线条，使线段相交处形成交点
保存编辑。
3 路网数据处理
3.1 原始数据入库
使用PGIS自带工具入库
点击 View connection details 按钮连接数据库
单击 Add File 可以选择需要入库的 Shapefile文件，可以设置数据库模式（Schema）、数据表的名称（Table，默认为Shapefile文件名称）、Geo Column名称（几何列名称，默认为geom，设置为the_geom）、空间参考ID（SRID，默认为0，设置为 4326）。设置完成后，点击Import导入。
4 入库的数据处理
我们已经处理好数据，现在开始入库并在数据库中处理数据入库的路网数据。
4.1 创建空间字段
在导入数据时使用了“the_geom”指定了原始数据的空间字段名称。创建新的空间字段:
--删除shp线图层自带的id字段
alter table huaxi_road drop column id;
--创建新的空间字段geom
select AddGeometryColumn ('huaxi_road','geom',4326,'LINESTRING',2);
--创建索引
create index gidx_huaxi_road_geom on huaxi_road using gist(geom);
4.2 检查无效的几何对象
select gid from huaxi_road where ST_IsValid(the_geom) IS FALSE;
如果有无效的几何对象，执行下面的命令删除，否则跳过
with cte as(
select gid from huaxi_road where ST_IsValid(the_geom) IS FALSE
) delete from huaxi_road where gid=ANY((select array_agg(gid) from cte)::integer[]);
4.3 将线段转换为Line
update huaxi_road set geom=ST_LineMerge(the_geom);
删除原有的空间字段
ALTER TABLE huaxi_road DROP COLUMN the_geom;
4.4 添加路网分析必须的字段
ALTER TABLE huaxi_road
ADD COLUMN source integer, /*当前线段起点连接至上一线段的id*/
ADD COLUMN target integer, /*当前线段终点连接至下一线段的id*/
ADD COLUMN cost double precision, /*正向成本*/
ADD COLUMN cost_time double precision, /*正向成本所需的时间*/
ADD COLUMN rcost double precision, /*反向成本*/
ADD COLUMN rcost_time double precision, /*反向成本所需的时间*/
ADD COLUMN x1 double precision, /*当前线段起点坐标(x)*/
ADD COLUMN y1 double precision, /*当前线段起点坐标(Y)*/
ADD COLUMN x2 double precision, /*当前线段终点坐标(x)*/
ADD COLUMN y2 double precision, /*当前线段终点坐标(y)*/
ADD COLUMN to_cost double precision,
ADD COLUMN rule text,
ADD COLUMN isolated integer;
4.5 更新属性字段
with base as(
select 'SPHEROID["WGS84",6378137,298.25728]'::spheroid as sph
) update huaxi_road set x1 = st_x(st_startpoint(geom)),
y1 = st_y(st_startpoint(geom)),
x2 = st_x(st_endpoint(geom)),
y2 = st_y(st_endpoint(geom)),
cost = ST_LengthSpheroid(geom, f.sph)::integer,
rcost = ST_LengthSpheroid(geom, f.sph)::integer
from base as f;
VACUUM FULL ANALYZE VERBOSE huaxi_road;
4.6 路网创建拓扑
select pgr_createTopology('huaxi_road', 0.000001, the_geom:='geom', id:='gid', source:='source', target:='target');
VACUUM FULL ANALYZE VERBOSE huaxi_road;
VACUUM FULL ANALYZE VERBOSE huaxi_road_vertices_pgr;
至此所有数据准备工作都已经完成了。
4.7 创建最短路径查询函数
CREATE OR REPLACE FUNCTION "public"."pgr_fromatob"("tbl" varchar, "startx" float8, "starty" float8, "endx" float8, "endy" float8)
RETURNS "public"."geometry" 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
--查询离起点最近的线
execute 'select geom, source, target from ' ||tbl||
' where ST_DWithin(geom,ST_Geometryfromtext(''point('|| startx ||' ' || starty||')'',4326),150)
order by ST_Distance(geom,ST_GeometryFromText(''point('|| startx ||' '|| starty ||')'',4326)) limit 1'
into v_startLine, v_startSource ,v_startTarget;
--查询离终点最近的线
execute 'select geom, source, target from ' ||tbl||
' where ST_DWithin(geom,ST_Geometryfromtext(''point('|| endx || ' ' || endy ||')'',4326),150)
order by ST_Distance(geom,ST_GeometryFromText(''point('|| endx ||' ' || endy ||')'',4326)) 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 ;
-- ST_Distance
--从开始的起点到结束的起点最短路径
execute 'SELECT st_linemerge(st_union(b.geom)) ' ||
'FROM pgr_dijkstra(
''SELECT gid as id, source, target, cost FROM ' || tbl ||''','
||v_startSource || ', ' ||v_endSource||' ,false
) a, '
|| tbl || ' b
WHERE a.edge=b.gid ' into v_res ;
--从开始的终点到结束的起点最短路径
execute 'SELECT st_linemerge(st_union(b.geom)) ' ||
'FROM pgr_dijkstra(
''SELECT gid as id, source, target, cost FROM ' || tbl ||''','
||v_startTarget || ', ' ||v_endSource||' ,false
) a, '
|| tbl || ' b
WHERE a.edge=b.gid ' into v_res_b ;
--从开始的起点到结束的终点最短路径
execute 'SELECT st_linemerge(st_union(b.geom)) ' ||
'FROM pgr_dijkstra(
''SELECT gid as id, source, target, cost FROM ' || tbl ||''','
||v_startSource || ', ' ||v_endTarget||' ,false
) a, '
|| tbl || ' b
WHERE a.edge=b.gid ' into v_res_c ;
--从开始的终点到结束的终点最短路径
execute 'SELECT st_linemerge(st_union(b.geom)) ' ||
'FROM pgr_dijkstra(
''SELECT gid as id, source, target, cost FROM ' || tbl ||''','
||v_startTarget || ', ' ||v_endTarget||' ,false
) a, '
|| tbl || ' b
WHERE a.edge=b.gid ' 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进行拼接
select st_linemerge(ST_Union(array[v_res,v_startLine,v_endLine])) into v_res;
select ST_LineLocatePoint(v_res, ST_Geometryfromtext('point('|| startx ||' ' || starty ||')',4326)) into v_perStart;
select ST_LineLocatePoint(v_res, ST_GeometryFromText('point('|| endx ||' ' || endy ||')',4326)) 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_LineSubstring(v_res,v_perStart, v_perEnd) into v_shPath;
return v_shPath;
end;
$BODY$
LANGUAGE plpgsql VOLATILE STRICT
COST 100
5 发布Geoserver服务
新建一个数据源类型为“PostGIS”的数据存储，然后配置PostGIS的数据库信息
在新建图层中新增SQL视图
配置如下：
SELECT * FROM pgr_fromatob('huaxi_road', %x1%, %y1%, %x2%, %y2%)
验证的正则表达式：^-?[\d.]+$
类型：LingString
SRID：4326
保存视图
保存生成图层。
6 请求Geoserver服务
通过ajax请求Geoserver服务，使用Cesium绘制路径
//服务地址
var host = 'localhost:10100';
//存储名称
var store = 'wms1_ws';
//图层标题
// var layer = 'jiangan_route';
var layer = 'huaxi_route4';
// var layer = 'wangjiang_route';
//图层名称
var typename = store + ':' + layer;
var url = 'http://' + host + '/geoserver/' + store + '/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=' + typename + '&maxFeatures=50&outputFormat=application/json';
var params = [];
//起点坐标
var x1 = 103.991127;
var y1 = 30.556150;
//终点坐标
var x2 = 104.001861;
var y2 = 30.561856;
//huaxi
// var x1 = 104.061424;
// var y1 = 30.643593;
// var x2 = 104.068772;
// var y2 = 30.641216;
//wangjiang
// var x1 = 104.075436;
// var y1 = 30.634966;
// var x2 = 104.084869;
// var y2 = 30.6333556;
console.log(url);
//查询参数
// url +='&viewparams=x1:\''+x1+'\';y1:\''+y1+'\';x2:\''+x2+'\';y2:\''+y2+'\'';
url += '&viewparams=x1:' + x1 + ';y1:' + y1 + ';x2:' + x2 + ';y2:' + y2;
var pos = new Array();
$.ajax({
type: 'post',
url: url,
async: false,
success: function (data) {
console.log(data);
if (data.features) {
console.log(data.features[0].geometry.coordinates);
var lines = data.features[0].geometry.coordinates;
for (var i = 0; i < lines.length; i++) {
// console.log(lines[i][0]+" "+lines[i][1]);
pos.push(lines[i][0], lines[i][1]);
}
var entity = viewer.entities.add({
name: '寻路',
polyline: {
positions: Cesium.Cartesian3.fromDegreesArray(pos),
width: 3,
material: Cesium.Color.GOLD,
clampToGround: false,
followSurface: false
}
});
}
}
});
参考
https://www.cnblogs.com/liongis/p/7692477.html
https://blog.****.net/kmblack1/article/details/78840274