使用PostGIS完成兩點間的河流軌跡及流經(jīng)長度的計算(推薦)
基礎準備工作
1.PostGIS 的安裝
在安裝PostGIS前首先必須安裝PostgreSQL,然后再安裝好的Stack Builder中選擇安裝PostGIS組件。具體安裝步驟可參照PostGIS的安裝與初步使用
2.加載Post GIS擴展
選中指定數(shù)據(jù)庫,執(zhí)行加載擴展語句
–添加支持 CREATE EXTENSION postgis; --添加postgis擴展 CREATE EXTENSION pgrouting; --添加pgrouting擴展 CREATE EXTENSION postgis_topology; CREATE EXTENSION fuzzystrmatch; CREATE EXTENSION postgis_tiger_geocoder;
在做兩點間河流軌跡及流經(jīng)長度計算過程中,需要加載postgis和pgrouting兩個擴展
可以通過查看加載擴展的版本驗證擴展加載是否成功
–查看postgresql版本 show server_version; –查看postgis版本 SELECT PostGIS_full_version(); –查看pgrouting版本 select pgr_version();
3.河流矢量圖層轉(zhuǎn)成單線格式
河流包括各種匯入和匯出,為了實現(xiàn)流經(jīng)流域的計算,河流水系矢量數(shù)據(jù)需要一個河流一個ID的方式,可以在河流交匯點處將河流進行打段處理。

4.河流矢量數(shù)據(jù)導入PostgreSQL數(shù)據(jù)庫
打開位于“開始>所有程序>PostGIS 2.3 bundle for PostgreSQL”之中的PostGIS Shapefile Import/Export Manager。
首先單擊"View connection details"按鈕,打開"PostGIS connection"對話框,輸入用戶名"postgres"及其對應的密碼,設置連接的數(shù)據(jù)庫,如下圖所示:

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

5.河流數(shù)據(jù)拓撲處理
在數(shù)據(jù)分析過程中,使用到了pgrouting擴展中的 pgr_dijkstra 算法
Dijkstra算法(迪杰斯特拉算法),由荷蘭計算機科學家Edsger Dijkstra于1956年提出。它是一種圖搜索算法,它解決了非負代價邊路徑圖的最短路徑問題,即從起始頂點(start_vid)到結(jié)束頂點(end_vid)的最短路徑。此算法可以與有向圖或無向圖一起使用。
函數(shù)的簽名摘要:

在實際使用中,需要先明確所有的頂點,并為所有頂點分配唯一的編號,函數(shù)的 start_vid 和 end_vid 都是整型數(shù)值,函數(shù)使用edges_sql參數(shù)(sql腳本)篩選出和頂點相鄰的所有邊信息(即河流信息)。
所以,在使用pgr_dijkstra方法前,需要
- 對找到河流的所有頂點信息,并做唯一整型值編號
- 在數(shù)據(jù)庫中為每條河流設置好起始頂點和結(jié)束頂點
--篩選出所有頂點信息,st_dump函數(shù)主要是將MultiLineString類型 調(diào)整成 LineString類型 select st_astext(st_startpoint((ST_Dump(geom)).geom)) from singleriver union select st_astext(st_endpoint((ST_Dump(geom)).geom)) from singleriver
將查詢結(jié)果在Excel中進行整型值編號,再導入到postgresql中的新建表distinctpoint 中,然后關聯(lián)河流數(shù)據(jù)表,更新河流的開始頂點(source)和結(jié)束頂點編號(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
--插入結(jié)束頂點編號 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
至此,河流拓撲數(shù)據(jù)處理完成
PG分析處理函數(shù)
1.函數(shù)編寫
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;--最短路徑分析結(jié)果
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_geomfromtext('point('|| startx ||' '|| starty ||') '), v_startLine) and ST_Intersects(st_geomfromtext('point('|| endx ||' '|| endy ||') '), v_startLine)) then
select st_distancesphere(v_statpoint, v_endpoint) into first_cost;
select st_linelocatepoint(st_linemerge(v_startLine), v_endpoint) into v_perEnd;
for v_rec in
select st_linesubstring(st_linemerge(v_startLine), v_perStart,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;
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 10002.參數(shù)說明
輸入?yún)?shù):開始點和結(jié)束點的經(jīng)緯度坐標
輸出結(jié)果:river_name:河流名稱;v_shppath:流經(jīng)的河流路徑; cost:河流流經(jīng)長度
3.內(nèi)部調(diào)用函數(shù)說明
函數(shù)調(diào)用過程,根據(jù)postgis不同版本,函數(shù)名稱可能會有偏差,有版本展示形式為st_linesubstring ,有版本展示形式為st_line_substring
4.輸出結(jié)果驗證
為了驗證河流輸出結(jié)果是否正確,流經(jīng)河流路徑是否連通,可以通過在線geojson地圖(geojson.io)呈現(xiàn)出來驗證。
在右側(cè)json-features-geometry 中填充函數(shù)輸出的v_shppath參數(shù)內(nèi)容(按照行單獨輸入,可以輸入多個,注意需要增加json屬性)

到此這篇關于使用PostGIS完成兩點間的河流軌跡及流經(jīng)長度的計算的文章就介紹到這了,更多相關PostGIS兩點間的河流軌跡計算內(nèi)容請搜索腳本之家以前的文章或繼續(xù)瀏覽下面的相關文章希望大家以后多多支持腳本之家!
相關文章
Abp.NHibernate連接PostgreSQl數(shù)據(jù)庫的方法
這篇文章主要為大家詳細介紹了Abp.NHibernate連接PostgreSQl數(shù)據(jù)庫的方法,具有一定的參考價值,感興趣的小伙伴們可以參考一下2018-01-01
解決postgresql 數(shù)字轉(zhuǎn)換成字符串前面會多出一個空格的問題
這篇文章主要介紹了解決postgresql 數(shù)字轉(zhuǎn)換成字符串前面會多出一個空格的問題,具有很好的參考價值,希望對大家有所幫助。一起跟隨小編過來看看吧2020-12-12
PostgreSQL pg_ctl start啟動超時實例分析
這篇文章主要給大家介紹了關于PostgreSQL pg_ctl start啟動超時的相關資料,文中通過示例代碼介紹的非常詳細,對大家的學習或者工作具有一定的參考學習價值,需要的朋友們下面隨著小編來一起學習學習吧2019-01-01
PostgreSQL使用MySQL作為外部表(mysql_fdw)
PostgreSQL 提供了一種訪問和操作外部數(shù)據(jù)源的機制,稱為外部數(shù)據(jù)包裝器,本文主要給大家介紹了PostgreSQL使用MySQL作為外部表的方法,感興趣的朋友跟隨小編一起看看吧2022-11-11
Navicat連接postgresql時出現(xiàn)'datlastsysoid?does?not?exist&
這篇文章主要給大家介紹了關于Navicat連接postgresql時出現(xiàn)'datlastsysoid?does?not?exist'報錯問題的完美解決辦法,文中通過圖文介紹的非常詳細,需要的朋友可以參考下2024-02-02
PostgreSQL進行數(shù)據(jù)導入和導出的操作代碼
在數(shù)據(jù)庫管理中,數(shù)據(jù)的導入和導出是非常常見的操作,特別是在 PostgreSQL 中,提供了多種工具和方法來實現(xiàn)數(shù)據(jù)的有效管理,本文將詳細介紹在 PostgreSQL 中如何進行數(shù)據(jù)導入和導出,并給出具體的命令及示例,需要的朋友可以參考下2024-10-10
PostgreSQL 重復數(shù)據(jù)處理的操作方法
這篇文章主要介紹了PostgreSQL 重復數(shù)據(jù)處理的操作方法,本文給大家介紹的非常詳細,對大家的學習或工作具有一定的參考借鑒價值,需要的朋友可以參考下2020-12-12

