使用PostGIS完成两点间的河流轨迹及流经长度的计算

时间：2022-01-18 作者：HoEn

使用POSTGIS及其扩展pgrouting计算给定两点间的河流流经区域和河流长度。

基础准备工作

域名GIS 的安装

在安装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=域名cepoint
from singleriver s,
(select gid,域名 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 域名tpoint=域名t) tt
where 域名=域名

--插入结束顶点编号
update singleriver q
set target=域名oint
from singleriver s,
(select gid,域名 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 域名oint=域名t) tt
where 域名=域名

至此，河流拓扑数据处理完成

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||\')\'\'),域名) 
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 ||\')\'\'),域名) 
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(域名t);
cost:= 域名;
river_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(域名) as point,COALESCE(域名,\'无名河流\') 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 ||\')\'\'),域名))\', 
v_startTarget, v_endSource , false 
) a, singleriver b 
WHERE 域名 = 域名
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(域名t);
cost:= 域名;
river_name:= 域名;
return next;
end loop; 
end; 
$BODY$
  LANGUAGE plpgsql VOLATILE STRICT
  COST 100
  ROWS 1000