2020/10/14，利用Postgis数据库对LineString的Shp文件进行路径规划，利用Geoserver对规划的结果进行发布，在Mapbox中使用WMS服务调用并呈现在地图上。
本文是以经纬度坐标为依据进行计算，参考其他网上已有教程可以选择以节点进行计算，其他链接放置本文底部。
软件版本
- Geoserver：geoserver-2.17.2-bin；
- PostGIS：2.1.7
- PostgreSQL：9.3
- QGIS：3.14.15
操作步骤
- 在
pgAdmin中导入shp文件，本次示例以roadLine文件为例，效果如图：
坐标系如图：
文件字段如图，框选内容字段是之前实验其他博文所添加，在本文中没有作用：
使用pgAdmin连接到数据库，直接导入文件，这里的模式选择为public，表名为gyroad，坐标系设置4326。
注意：红框内容强调，最后一个必须勾选才能生成单个的geometry。
- 数据导入完毕之后会在当前数据库的
public模式生成表，
接下来对数据库使用SQL语句，这里采用QGIS的SQL窗口，pgAdmin下的SQL也可以使用。
--向数据库导入支持PostGIS和pgRouting的函数和基础表（安装后只需要执行一次，以后会一直存在于本数据库中）。
CREATE EXTENSION postgis;
CREATE EXTENSION pgrouting;
CREATE EXTENSION postgis_topology;
CREATE EXTENSION fuzzystrmatch;
CREATE EXTENSION postgis_tiger_geocoder;
CREATE EXTENSION address_standardizer;
-----------------------------------------
--以起始点坐标为参数进行最短路径分析
--以gyroad表作为实例--
-----------------------------------------
ALTER TABLE gyroad ADD COLUMN source integer; --设置起点字段
ALTER TABLE gyroad ADD COLUMN target integer; --设置终点字段
ALTER TABLE gyroad ADD COLUMN length double precision; --增加路线长度字段（根据长度设置计算权重）
UPDATE gyroad SET length = ST_Length(geom); --计算路线长度
select pgr_createTopology('gyroad', 0.0001, 'geom', 'gid'); --创建拓扑模式
--为表增加起始点坐标x,y字段
ALTER TABLE gyroad ADD COLUMN x1 double precision;
ALTER TABLE gyroad ADD COLUMN y1 double precision;
ALTER TABLE gyroad ADD COLUMN x2 double precision;
ALTER TABLE gyroad ADD COLUMN y2 double precision;
--计算起始点坐标
UPDATE gyroad SET x1 =ST_x(ST_PointN(geom, 1));
UPDATE gyroad SET y1 =ST_y(ST_PointN(geom, 1));
UPDATE gyroad SET x2 =ST_x(ST_PointN(geom, ST_NumPoints(geom)));
UPDATE gyroad SET y2 =ST_y(ST_PointN(geom, ST_NumPoints(geom)));
--测试1
--A*算法路径查询
SELECT seq, id1 AS node, id2 AS edge, cost FROM pgr_astar('SELECT gid AS id,
source::integer,target::integer,length::double precision AS cost,
x1, y1, x2, y2 FROM gyroad',30, 60, false,false);
--查询从出发点到目的地的消耗
SELECT seq, id1 AS source, id2 AS target,cost FROM pgr_kdijkstraCost('SELECT gid AS id,
source::integer,target::integer,length::double precision AS cost FROM gyroad',30,
array[60,70,100], false, false);
--pgr_kdijkstraPath函数路径分析
SELECT seq, id1 AS path, id2 AS edge, cost FROM pgr_kdijkstraPath('SELECT gid AS id,
source::integer,target::integer,length::double precision AS cost
FROM gyroad',30, array[60,100], false, false);
--清除function pgr_fromAtoB(var)
--若之前没有设定该函数可不执行一些代码
--DROP FUNCTION pgr_fromAtoB(varchar, double precision, double precision,double precision, double precision);
--基于任意两点之间的最短路径分析
CREATE OR REPLACE FUNCTION pgr_fromAtoB(
IN tbl varchar,--数据库表名
IN x1 double precision,--起点x坐标
IN y1 double precision,--起点y坐标
IN x2 double precision,--终点x坐标
IN y2 double precision,--终点y坐标
OUT seq integer,--道路序号
OUT gid integer,
OUT name text,--道路名 注意：本文实验的数据中不存在name字段 若需要使用可以删除本函数中所有关于name的设定，也可以修改原数据增加name字段。以下函数中删除name的部分以（DELECT）为标记。
OUT heading double precision,
OUT cost double precision,--消耗
OUT geom geometry--道路几何集合
)
RETURNS SETOF record AS
$BODY$
DECLARE
sql text;
rec record;
source integer;
target integer;
point integer;
BEGIN
-- 查询距离出发点最近的道路节点
EXECUTE 'SELECT id::integer FROM '|| quote_ident(tbl) ||'_vertices_pgr
ORDER BY the_geom <-> ST_GeometryFromText(''POINT('
|| x1 || ' ' || y1 || ')'',4326) LIMIT 1' INTO rec;
source := rec.id;
-- 查询距离目的地最近的道路节点
EXECUTE 'SELECT id::integer FROM '|| quote_ident(tbl) ||'_vertices_pgr
ORDER BY the_geom <-> ST_GeometryFromText(''POINT('
|| x2 || ' ' || y2 || ')'',4326) LIMIT 1' INTO rec;
target := rec.id;
-- 最短路径查询
seq := 0;
-- （DELECT）
sql := 'SELECT gid, geom, name, 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;
name := rec.name; -- （DELECT）
cost := rec.cost;
geom := rec.geom;
RETURN NEXT;
END LOOP;
RETURN;
END;
$BODY$
LANGUAGE 'plpgsql' VOLATILE STRICT;
--测试函数
--注意：模式若不在public中，会出现找不到表的错误情况，需要修改以上函数从具体路径找表，本人SQL较生疏，望交流反馈。
SELECT st_astext(ST_MakeLine(route.geom))
FROM(SELECT seq,gid,heading,cost,geom
FROM pgr_fromAtoB('gyroad',119.21182, 34.60405, 119.21695, 34.60855)ORDER BY seq) AS route;
以文本的形式进行查询，但是在Geoserver创建SQL视图的时候需要修改SQL语句，如下：
SELECT seq,gid,heading,cost,geom
FROM pgr_fromAtoB('gyroad', %x1%, %y1%, %x2%, %y2%) ORDER BY seq
若在上方能够成功出现st_astext的内容，可以进行下一步Geoserver的发布。
- 利用
Geoserver进行视图URL的发布。
创建一个新的数据存储，选择PostGIS模式，输入数据库名称、账号密码和模式名称，本文之前创立过数据存储，这里不再展示。
在下一步的新建图层界面，直接发布链接到数据库中的gyroad图层，选择4326坐标系，之后重新回到本页面配置新的SQL视图。
视图名称可以随意设置，SQL语句在上文中提到可以直接使用；从SQL猜想的参数分别是始末经纬度坐标，默认值设置为0，修改验证的正则表达式为^-?[\d.]+$。
刷新属性内容，设置geom为LineString ，SRID保持一致为4326。
保存后设置SRID、经纬度边框等，这里不再赘述。
在图层预览中查看之前发布的视图
效果如图所示，直线的原因是之前设置的参数默认值为0；修改URL地址，在原URL的后面增加viewparams数据，以此为例：&viewparams=x1:119.21114;y1:34.60776;x2:119.21692;y2:34.60626。
点击地图可以获得要素信息
- 在Mapbox中调用Geoserver的WMS服务。本质还是根据Mapbox官方的WMS example进行服务调用。问题在于如何获取始末点的经纬度坐标，这里选择设定
function来完成相关操作。全部代码如下：
<!DOCTYPE html>
<html lang="en">
<head>
<meta charset="UTF-8">
<title>MapBox加载WMS地图服务</title>
<script src='https://api.mapbox.com/mapbox-gl-js/v0.50.0/mapbox-gl.js'></script>
<link href='https://api.mapbox.com/mapbox-gl-js/v0.50.0/mapbox-gl.css' rel='stylesheet' />
<style>
body {
margin: 0;
padding: 0;
}
#map {
position: absolute;
top: 0;
bottom: 0;
width: 100%;
}
/*隐藏logo*/
.mapboxgl-ctrl.mapboxgl-ctrl-attrib {
display: none !important;
}
.mapboxgl-ctrl-logo {
display: none !important;
}
</style>
</head>
<body>
<div id="map"></div>
<script>
mapboxgl.accessToken =
'YOUR TOKEN';
var map = new mapboxgl.Map({
container: 'map',
style: 'mapbox://styles/mapbox/streets-v10',
center: [119.21, 34.604],
zoom: 7
});
function getRouting(x1, y1, x2, y2) {
let url =
`http://localhost:8080/geoserver/postNavigate/wms?service=WMS&version=1.1.0&request=GetMap&layers=postNavigate%3Agyroad_test&bbox={bbox-epsg-3857}&width=330&height=768&srs=EPSG%3A3857&format=image%2Fpng&TRANSPARENT=TRUE&viewparams=x1:${x1};y1:${y1};x2:${x2};y2:${y2}`;
console.log(url);
map.addLayer({
'id': 'wms-test-layer',
'type': 'raster',
'source': {
'type': 'raster',
'tiles': [url],
},
// 以下内容是分析URL地址各个属性。
// 'tiles': [
// //这里可以设置各个图层lw_base_mapbox
// 'http://localhost:8080/geoserver/postNavigate/wms?service=WMS' //wms地图服务地址
// +
// '&version=1.1.0' //wms的版本 网上有教程是1.3.0
// +
// '&request=getmap' //调用的方法名称，获取地图必须是这个方法
// +
// '&BGCOLOR=ff00ff' //生成的图片背景颜色
// //layers 使用ArcGIS进行发布的时候默认图层名称为0，1，2的索引值，在发布地图服务的时候可以勾选
// //“使用地图文档中的图层名称”，在访问的时候图层名称就是地图上图层名称，不然就是0，1，2，3的索引值
// +
// '&layers=postNavigate%3Agyroad_test' //要显示的图层名称，多个图层用,隔开；
// // +
// // '&styles=' //图层显示样式，同样多个样式名称间用，隔开；
// +
// // '&crs=EPSG:3857' +
// '&bbox={bbox-epsg-3857}' //使用map加载的wms的时候，使用这个标识来同步要获取地图的范围的坐标，必须为这个值
// +
// '&width=512' //返回的图片的大小
// +
// '&height=512' //返回的图片的大小
// +
// '&format=image/png' //返回的图片格式
// +
// `viewparams=x1:${x1};y1:${y1};x2:${x2};y2:${y2}` //可以使用字符串模板获取参数
// +
// '&TRANSPARENT=TRUE' //设置背景是不是可以透明，没有数据的地方就进行透明
// +
// '&SLD=http://120.77.2.35:8090/lw_base_mapbox.sld' //客户端可以自定义各个图层显示的样式，该文件定义各个图层显示样式
// //在使用ArcGIS发布服务的时候 ，可以设置“SLD 路径或 URL”来指定样式文件，这样默认显示样式就按照该样式文件进行显示:
// ],
// 'tileSize': 512 //图片显示的大小，最好和上面大小保持一致
// },
// 'paint': {
// "raster-opacity": 1, //图层显示透明度
// //raster-hue-rotate 设置该值以后，显示的颜色就不会是图层样式里面设置的颜色，所以最好不要设置
// //"raster-hue-rotate": 60,//在色轮上旋转色相的角度
// }
});
}
</script>
</body>
</html>
最后效果如图：
参考文档：
Mapbox 添加WMS服务
PostGIS pgrouting路径分析
基于PostGIS+PgRouting的最短路径查询的实现
5613
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
