使用PostGIS完成兩點間的河流軌跡及流經長度的計算

基礎准備工作

1.PostGIS 的安裝

在安裝PostGIS前首先必須安裝PostgreSQL，然后再安裝好的Stack Builder中選擇安裝PostGIS組件。具體安裝步驟可參照 PostGIS的安裝與初步使用_不睡覺的怪叔叔的博客-CSDN博客_postgis

2.加載Post GIS擴展

選中指定數據庫，執行加載擴展語句

–添加支持
CREATE EXTENSION postgis;  --添加postgis擴展
CREATE EXTENSION pgrouting;   --添加pgrouting擴展
CREATE EXTENSION postgis_topology;
CREATE EXTENSION fuzzystrmatch;
CREATE EXTENSION postgis_tiger_geocoder;

在做兩點間河流軌跡及流經長度計算過程中，需要加載postgis和pgrouting兩個擴展

可以通過查看加載擴展的版本驗證擴展加載是否成功

–查看postgresql版本
show server_version;

–查看postgis版本
SELECT PostGIS_full_version();

–查看pgrouting版本
select pgr_version();

 

3.河流矢量圖層轉成單線格式

河流包括各種匯入和匯出，為了實現流經流域的計算，河流水系矢量數據需要一個河流一個ID的方式，可以在河流交匯點處將河流進行打段處理。

4.河流矢量數據導入PostgreSQL數據庫

打開位於“開始>所有程序>PostGIS 2.3 bundle for PostgreSQL”之中的PostGIS Shapefile Import/Export Manager。

首先單擊"View connection details"按鈕，打開"PostGIS connection"對話框，輸入用戶名"postgres"及其對應的密碼，設置連接的數據庫，如下圖所示：

 

 連接數據庫之后，單擊"Add file"按鈕，加入***.shp文件，並將其SRID設置為"4326"，如下圖所示。這一步絕對不能省略，否則不能正確導入數據。

5.河流數據拓撲處理

在數據分析過程中，使用到了pgrouting擴展中的 pgr_dijkstra 算法

Dijkstra算法（迪傑斯特拉算法），由荷蘭計算機科學家Edsger Dijkstra於1956年提出。它是一種圖搜索算法，它解決了非負代價邊路徑圖的最短路徑問題，即從起始頂點（start_vid）到結束頂點（end_vid）的最短路徑。此算法可以與有向圖或無向圖一起使用。

函數的簽名摘要：

在實際使用中，需要先明確所有的頂點，並為所有頂點分配唯一的編號，函數的 start_vid 和 end_vid 都是整型數值，函數使用edges_sql參數(sql腳本)篩選出和頂點相鄰的所有邊信息（即河流信息）。

所以，在使用pgr_dijkstra方法前，需要

• 對找到河流的所有頂點信息，並做唯一整型值編號
• 在數據庫中為每條河流設置好起始頂點和結束頂點

 

--篩選出所有頂點信息，st_dump函數主要是將MultiLineString類型 調整成 LineString類型
select  st_astext(st_startpoint((ST_Dump(geom)).geom)) from singleriver
union 
select  st_astext(st_endpoint((ST_Dump(geom)).geom)) from singleriver

 

將查詢結果在Excel中進行整型值編號，再導入到postgresql中的新建表 distinctpoint 中，然后關聯河流數據表，更新河流的開始頂點（source）和結束頂點編號（target）

 

--更新起始頂點編號
update singleriver q
set source=tt.sourcepoint
from singleriver s,
(select gid,p.id as sourcepoint from 
(select gid,st_astext(st_startpoint((ST_Dump(geom)).geom)) as startpoint, st_astext(st_endpoint((ST_Dump(geom)).geom)) as endpoint from singleriver )s
left join distinctpoint p
on s.startpoint=p.point) tt
where q.gid=tt.gid

 

--插入結束頂點編號
update singleriver q
set target=tt.endpoint
from singleriver s,
(select gid,p.id as endpoint from 
(select gid,st_astext(st_startpoint((ST_Dump(geom)).geom)) as startpoint, st_astext(st_endpoint((ST_Dump(geom)).geom)) as endpoint from singleriver )s
left join distinctpoint p
on s.endpoint=p.point) tt
where q.gid=tt.gid

至此，河流拓撲數據處理完成

PG分析處理函數

1.函數編寫

方式一：TABLE輸出

 

CREATE OR REPLACE FUNCTION "public"."pgr_shortest_river"(IN "startx" float8, IN "starty" float8, IN "endx" float8, IN "endy" float8)
  RETURNS TABLE ("river_name" varchar, "v_shpath" varchar, "cost" float8) AS $BODY$ 
declare 
--river_name varchar;
--v_shPath varchar;
--cost float8;
v_startLine geometry;--離起點最近的線 
v_endLine geometry;--離終點最近的線 
v_startTarget integer;--距離起點最近線的終點 
v_endSource integer;--距離終點最近線的起點 
v_statpoint geometry;--在v_startLine上距離起點最近的點 
v_endpoint geometry;--在v_endLine上距離終點最近的點 
v_res geometry;--最短路徑分析結果 
v_perStart float;--v_statpoint在v_res上的百分比 
v_perEnd float;--v_endpoint在v_res上的百分比 
v_rec record; 
first_name varchar;
end_name varchar;
first_cost double precision;
end_cost double precision;
begin 
--查詢離起點最近的線 
execute 'select (st_dump(geom)).geom as geom,target as target,name from singleriver where 
ST_DWithin(geom,ST_Geometryfromtext(''point('|| startx ||' ' || starty||')''),0.01) 
order by ST_Distance(geom,ST_GeometryFromText(''point('|| startx ||' '|| starty ||')'')) limit 1' 
into v_startLine ,v_startTarget,first_name; 
raise notice '河流名稱%',first_name;
--查詢離終點最近的線 
execute 'select (st_dump(geom)).geom as geom,"source" as source,name from singleriver
where ST_DWithin(geom,ST_Geometryfromtext(''point('|| endx || ' ' || endy ||')''),0.01) 
order by ST_Distance(geom,ST_GeometryFromText(''point('|| endx ||' ' || endy ||')'')) limit 1' 
into v_endLine,v_endSource,end_name; 
--如果沒找到最近的線，就返回null 
if (v_startLine is null) or (v_endLine is null) then 
return; 
end if ; 
select ST_ClosestPoint(v_startLine, ST_Geometryfromtext('point('|| startx ||' ' || starty ||')')) into v_statpoint; 
select ST_ClosestPoint(v_endLine, ST_GeometryFromText('point('|| endx ||' ' || endy ||')')) into v_endpoint; 

--計算距離起點最近線上的點在該線中的位置
select st_linelocatepoint(st_linemerge(v_startLine), v_statpoint) into v_perStart;

select st_linelocatepoint(st_linemerge(v_endLine), v_endpoint) into v_perEnd;

select st_distancesphere(v_statpoint,ST_PointN(ST_GeometryN(v_startLine,1), ST_NumPoints(ST_GeometryN(v_startLine,1)))) into first_cost;

select st_distancesphere(ST_PointN(ST_GeometryN(v_endLine,1),1),v_endpoint) into end_cost; 

--查詢的兩個位置在同一條河流上時
if (ST_Intersects(st_buffer(v_statpoint,0.0001), v_startLine) and ST_Intersects(st_buffer(v_endpoint,0.0001), v_startLine)) then 
select st_linelocatepoint(st_linemerge(v_startLine), v_endpoint) into v_perEnd;
for v_rec in 
select st_linesubstring(st_linemerge(v_startLine),  case when v_perStart<v_perEnd then v_perStart else v_perEnd end,case when v_perStart>v_perEnd then v_perStart else v_perEnd end) as point,COALESCE(end_name,'無名河流') as name,st_length ( st_transform ( st_setsrid( st_linesubstring ( v_startLine,  case when v_perStart<v_perEnd then v_perStart else v_perEnd end,case when v_perStart>v_perEnd then v_perStart else v_perEnd end), 4326), 3857 ) ) as cost
loop
v_shPath:= ST_AsGeoJSON(v_rec.point);
cost:= v_rec.cost;
river_name:= v_rec.name;
return next;
end loop;
return;
end if;
--最短路徑 
for v_rec in 
(select st_linesubstring(st_linemerge(v_startLine),v_perStart,1) as point,COALESCE(first_name,'無名河流') as name,first_cost as cost
union all
SELECT st_linemerge(b.geom) as point,COALESCE(b.name,'無名河流') as name,st_length(geom, false) as cost
FROM pgr_dijkstra(
'SELECT gid as id, source, target, st_length(geom, false) as cost FROM singleriver
where st_intersects(geom,st_buffer(st_linefromtext(''linestring('||startx||' ' || starty ||','|| endx ||' ' || endy ||')''),0.05))', 
v_startTarget, v_endSource , false 
) a, singleriver b 
WHERE a.edge = b.gid
union all
select st_linesubstring(st_linemerge(v_endLine),0,v_perEnd) as point,COALESCE(end_name,'無名河流') as name,end_cost as cost)
loop
v_shPath:= ST_AsGeoJSON(v_rec.point);
cost:= v_rec.cost;
river_name:= v_rec.name;
return next;
end loop; 
end; 
$BODY$
  LANGUAGE plpgsql VOLATILE STRICT
  COST 100
  ROWS 1000

 

函數調用方式

select * from pgr_shortest_river(114.34339155124664,27.85120700316823,114.44430856192854,27.839804751500512)

  

方式二： setof record 輸出

 

CREATE OR REPLACE FUNCTION "public"."pgr_shortest_river"(IN "startx" float8, IN "starty" float8, IN "endx" float8, IN "endy" float8, OUT "river_name" varchar, OUT "v_shpath" varchar, OUT "cost" float8)
  RETURNS SETOF "pg_catalog"."record" AS $BODY$ 
declare 
v_startLine geometry;--離起點最近的線 
v_endLine geometry;--離終點最近的線 
v_startTarget integer;--距離起點最近線的終點 
v_endSource integer;--距離終點最近線的起點 
v_statpoint geometry;--在v_startLine上距離起點最近的點 
v_endpoint geometry;--在v_endLine上距離終點最近的點 
v_res geometry;--最短路徑分析結果 
v_perStart float;--v_statpoint在v_res上的百分比 
v_perEnd float;--v_endpoint在v_res上的百分比 
v_rec record; 
first_name varchar;
end_name varchar;
first_cost double precision;
end_cost double precision;
begin 
--查詢離起點最近的線 
execute 'select (st_dump(geom)).geom as geom,target as target,name from singleriver where 
ST_DWithin(geom,ST_Geometryfromtext(''point('|| startx ||' ' || starty||')''),0.01) 
order by ST_Distance(geom,ST_GeometryFromText(''point('|| startx ||' '|| starty ||')'')) limit 1' 
into v_startLine ,v_startTarget,first_name; 
raise notice '起點線段%',v_startLine;
raise notice '起點位置%',v_startTarget;
raise notice '河流名稱%',first_name;
--查詢離終點最近的線 
execute 'select (st_dump(geom)).geom as geom,"source" as source,name from singleriver
where ST_DWithin(geom,ST_Geometryfromtext(''point('|| endx || ' ' || endy ||')''),0.01) 
order by ST_Distance(geom,ST_GeometryFromText(''point('|| endx ||' ' || endy ||')'')) limit 1' 
into v_endLine,v_endSource,end_name; 
--如果沒找到最近的線，就返回null 
if (v_startLine is null) or (v_endLine is null) then 
return; 
end if ; 
select ST_ClosestPoint(v_startLine, ST_Geometryfromtext('point('|| startx ||' ' || starty ||')')) into v_statpoint; 
select ST_ClosestPoint(v_endLine, ST_GeometryFromText('point('|| endx ||' ' || endy ||')')) into v_endpoint; 

--計算距離起點最近線上的點在該線中的位置
select st_linelocatepoint(st_linemerge(v_startLine), v_statpoint) into v_perStart;

select st_linelocatepoint(st_linemerge(v_endLine), v_endpoint) into v_perEnd;

select st_distancesphere(v_statpoint,ST_PointN(ST_GeometryN(v_startLine,1), ST_NumPoints(ST_GeometryN(v_startLine,1)))) into first_cost;

select st_distancesphere(ST_PointN(ST_GeometryN(v_endLine,1),1),v_endpoint) into end_cost; 

--查詢的兩個位置在同一條河流上時
if (ST_Intersects(st_buffer(v_statpoint,0.0001), v_startLine) and ST_Intersects(st_buffer(v_endpoint,0.0001), v_startLine)) then 
select st_linelocatepoint(st_linemerge(v_startLine), v_endpoint) into v_perEnd;
for v_rec in 
select st_linesubstring(st_linemerge(v_startLine),  case when v_perStart<v_perEnd then v_perStart else v_perEnd end,case when v_perStart>v_perEnd then v_perStart else v_perEnd end) as point,COALESCE(end_name,'無名河流') as name,st_length ( st_transform ( st_setsrid( st_linesubstring ( v_startLine,  case when v_perStart<v_perEnd then v_perStart else v_perEnd end,case when v_perStart>v_perEnd then v_perStart else v_perEnd end), 4326), 3857 ) ) as cost
loop
v_shPath:= ST_AsGeoJSON(v_rec.point);
cost:= v_rec.cost;
river_name:= v_rec.name;
return next;
end loop;
return;
end if;
--最短路徑 
for v_rec in 
(select st_linesubstring(st_linemerge(v_startLine),v_perStart,1) as point,COALESCE(first_name,'無名河流') as name,first_cost as cost
union all
SELECT st_linemerge(b.geom) as point,COALESCE(b.name,'無名河流') as name,st_length(geom, false) as cost
FROM pgr_dijkstra(
'SELECT gid as id, source, target, st_length(geom, false) as cost FROM singleriver
where st_intersects(geom,st_buffer(st_linefromtext(''linestring('||startx||' ' || starty ||','|| endx ||' ' || endy ||')''),0.05))', 
v_startTarget, v_endSource , false 
) a, singleriver b 
WHERE a.edge = b.gid
union all
select st_linesubstring(st_linemerge(v_endLine),0,v_perEnd) as point,COALESCE(end_name,'無名河流') as name,end_cost as cost)
loop
v_shPath:= ST_AsGeoJSON(v_rec.point);
cost:= v_rec.cost;
river_name:= v_rec.name;
return next;
end loop; 
end; 
$BODY$
  LANGUAGE plpgsql VOLATILE STRICT
  COST 100
  ROWS 1000

 

 

 

函數調用方式

select pgr_shortest_river(114.34339155124664,27.85120700316823,114.44430856192854,27.839804751500512)

 

2.參數說明

輸入參數：開始點和結束點的經緯度坐標

輸出結果：river_name:河流名稱；v_shppath:流經的河流路徑； cost：河流流經長度

3.內部調用函數說明

函數調用過程，根據postgis不同版本，函數名稱可能會有偏差，有版本展示形式為  st_linesubstring ，有版本展示形式為 st_line_substring

4.輸出結果驗證

為了驗證河流輸出結果是否正確，流經河流路徑是否連通，可以通過在線geojson地圖（geojson.io）呈現出來驗證。

在右側json-features-geometry 中填充函數輸出的v_shppath參數內容（按照行單獨輸入，可以輸入多個，注意需要增加json屬性）

右側json數據：https://files.cnblogs.com/files/HoEn/%E6%B2%B3%E6%B5%81%E6%95%B0%E6%8D%AE%E6%B5%8B%E8%AF%95.json