基礎準備工作
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.函式編寫
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_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 1000
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