使用PostGIS完成两点间的河流轨迹及流经长度的计算(推荐)
这篇文章主要介绍了使用PostGIS完成两点间的河流轨迹及流经长度的计算,使用POSTGIS及其扩展pgrouting计算给定两点间的河流流经区域和河流长度,需要的朋友可以参考下
目录
基础准备工作
1.PostGIS 的安装
2.加载Post GIS扩展
3.河流矢量图层转成单线格式
4.河流矢量数据导入PostgreSQL数据库
5.河流数据拓扑处理
PG分析处理函数
1.函数编写
2.参数说明
3.内部调用函数说明
4.输出结果验证
基础准备工作
1.PostGIS 的安装
在安装PostGIS前首先必须安装PostgreSQL,然后再安装好的Stack Builder中选择安装PostGIS组件。具体安装步骤可参照PostGIS的安装与初步使用
2.加载Post GIS扩展
选中指定数据库,执行加载扩展语句
1 2 3 4 5 6 | –添加支持 CREATE EXTENSION postgis; --添加postgis扩展 CREATE EXTENSION pgrouting; --添加pgrouting扩展 CREATE EXTENSION postgis_topology; CREATE EXTENSION fuzzystrmatch; CREATE EXTENSION postgis_tiger_geocoder; |
在做两点间河流轨迹及流经长度计算过程中,需要加载postgis和pgrouting两个扩展
可以通过查看加载扩展的版本验证扩展加载是否成功
1 2 3 4 5 6 7 8 | –查看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方法前,需要
对找到河流的所有顶点信息,并做唯一整型值编号
在数据库中为每条河流设置好起始顶点和结束顶点
1 2 3 4 | --筛选出所有顶点信息,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)
1 2 3 4 5 6 7 8 9 | --更新起始顶点编号 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 |
1 2 3 4 5 6 7 8 9 | --插入结束顶点编号 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.函数编写
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 | 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属性)
到此这篇关于使用PostGIS完成两点间的河流轨迹及流经长度的计算的文章就介绍到这了
原文链接:https://www.cnblogs.com/HoEn/p/15816576.html