PostGIS 爆管分析之根据爆点找出所有影响阀门

环境：

Win10

ArcMap10.4（用于数据处理）

postgresql9.4

postgis2.2.3

pgRouting2.3（postgresql插件）

说明：

做爆管分析的第一步，需要先将数据做拓扑处理（方法见博文《PostGIS 结合Openlayers以及Geoserver实现最短路径分析》，共三篇：https://www.cnblogs.com/giser-s/p/11599562.html）

以下在构建拓扑数据成功的基础上继续（保证gid、source、target字段）

步骤：

 1、用户设定的爆管点（startx，starty），会与真实管网位置有差距，这里设置15米容差：假设爆管点距离15米内，最近的管网为爆点管段

execute 'select geom, source, target, ST_StartPoint(geom) as startpoint,ST_EndPoint(geom) as endpoint from ' ||tbl|| 
                            ' where ST_DWithin(geom,ST_Geometryfromtext(''point('|| startx ||' ' || starty ||')'',3857),15)
                            order by ST_Distance(geom,ST_GeometryFromText(''point('|| startx ||' '|| starty ||')'',3857))  limit 1' 
                            into v_startLine, v_startSource ,v_startTarget, v_statpoint ,v_endpoint; 

 2、开始查找，从爆管开始

      --循环第一次，将爆管放入查询条件，查询爆管上是否有阀门
      IF(v_up_idx = 0) THEN
        SELECT array_append(v_up_where, v_startSource) into v_up_where;
      ELSE
        --v_up_where = null;
      END IF;  
   
      --循环开始
      FOR up_temprow IN 
        select zy1.gid,zy1.source,zy1.target from zy zy1 where source = any(v_up_where) or target = any(v_up_where) 
      LOOP

      --查询管网上的点
      select t.gid,t.geom from fm t where t.gid  in (
        select a.gid from fm a,(select c.* from zy c where c.gid = up_temprow.gid) b where ST_intersects(a.geom,b.geom) 
      ) into v_uptap_gid, v_uptap_geom;  

3、如爆管没有阀门，则继续循环下一层级，往下找与爆管相接的管段；如果有阀门，则返回阀门，从下一次循环中剔除（不再找与他相接的管段）

        --如果没查找到阀门，则继续往下查
        IF(v_uptap_gid is null) then
          --source去重，判断如果数组中已有，则不添加
          IF (v_up_where @> ARRAY[up_temprow.source::integer] OR v_all_where @> ARRAY[up_temprow.source::integer]) THEN
          ELSE
              SELECT array_append(v_up_where,up_temprow.source) into v_up_where;
              SELECT array_append(v_all_where,up_temprow.source) into v_all_where;
          END IF;
          --target去重，判断如果数组中已有，则不添加
          IF (v_up_where @> ARRAY[up_temprow.target::integer] OR v_all_where @> ARRAY[up_temprow.target::integer]) THEN
          ELSE
              SELECT array_append(v_up_where,up_temprow.target) into v_up_where;
              SELECT array_append(v_all_where,up_temprow.target) into v_all_where;
          END IF;
        ELSE

4、如果有阀门，则返回阀门gid和geom

          --执行返回结果
          return query
          select v_uptap_gid as res_uptap_gid,v_uptap_geom as res_uptap_geom ;

5、附上全部存储过程

-- Function: test_getpoint5(character varying, double precision, double precision)
-- DROP FUNCTION test_getpoint5(character varying, double precision, double precision);
CREATE OR REPLACE FUNCTION test_getpoint5(
    IN tbl character varying,
    IN startx double precision,
    IN starty double precision)
  RETURNS TABLE(v_gid integer, v_res geometry) AS
$BODY$  

declare  
    v_startLine geometry;--离起点最近的线 
    v_startTarget integer;--距离起点最近线的终点 
    v_startSource integer; 
    v_statpoint geometry;--在v_startLine上距离起点最近的点  
    v_endpoint geometry;--在v_endLine上距离终点最近的点  
    v_up_source integer;--游标，记录是否有记录
    v_up_idx integer;--记录遍历到多少层级
    v_uptap_gid integer;--上游阀门gid
    v_uptap_geom geometry;--上游阀门要素
    v_all_where integer[];--记录所有查询过的管段
    v_up_where integer[];--where条件，将遍历到阀门的管段gid排除
    up_temprow record ;
    test integer;
begin 
    --查询离起点最近的线 
    --3857坐标系
    --找起点15米范围内的最近线 
    execute 'select geom, source, target, ST_StartPoint(geom) as startpoint,ST_EndPoint(geom) as endpoint from ' ||tbl|| 
                            ' where ST_DWithin(geom,ST_Geometryfromtext(''point('|| startx ||' ' || starty ||')'',3857),15)
                            order by ST_Distance(geom,ST_GeometryFromText(''point('|| startx ||' '|| starty ||')'',3857))  limit 1' 
                            into v_startLine, v_startSource ,v_startTarget, v_statpoint ,v_endpoint; 
    IF(v_startLine is not null) THEN
    --查找上游阀门
    v_up_idx = 0;
    v_up_source = 1;
    test = 0;
    SELECT array_append(v_up_where, v_startSource) into v_up_where;
    WHILE array_length(v_up_where,1) > 0 
    LOOP
      --游标归零
      v_up_source = 0; 
      --IF(v_up_idx = 0) THEN
        --SELECT array_append(v_up_where, v_startSource) into v_up_where;
        --SELECT array_append(v_up_where, v_startTarget) into v_up_where;
      --ELSE
        --v_up_where = null;
      --END IF;      
      --记录层级
      v_up_idx = v_up_idx + 1;
      --获取当前层级节点      
      FOR up_temprow IN 
        select zy1.gid,zy1.source,zy1.target from zy zy1 where source = any(v_up_where) or target = any(v_up_where) 
      LOOP
        test = test +1;
        --清空需要查的点
        IF(v_up_source = 0) THEN
          v_up_where = null;
        END IF;
        --清空初始执行节点
        v_startSource = 0;
        --标志执行有数据
        v_up_source = 1;
        --查询管网上的点
        select t.gid,t.geom from fm t where t.gid  in (
          select a.gid from fm a,(select c.* from zy c where c.gid = up_temprow.gid) b where ST_intersects(a.geom,b.geom) 
        ) into v_uptap_gid, v_uptap_geom;  
        
        --如果没查找到阀门，则继续往下查
        IF(v_uptap_gid is null) then
          --source去重，判断如果数组中已有，则不添加
          IF (v_up_where @> ARRAY[up_temprow.source::integer] OR v_all_where @> ARRAY[up_temprow.source::integer]) THEN
          ELSE
              SELECT array_append(v_up_where,up_temprow.source) into v_up_where;
              SELECT array_append(v_all_where,up_temprow.source) into v_all_where;
          END IF;
          --target去重，判断如果数组中已有，则不添加
          IF (v_up_where @> ARRAY[up_temprow.target::integer] OR v_all_where @> ARRAY[up_temprow.target::integer]) THEN
          ELSE
              SELECT array_append(v_up_where,up_temprow.target) into v_up_where;
              SELECT array_append(v_all_where,up_temprow.target) into v_all_where;
          END IF;
        ELSE
          raise notice '%' , v_uptap_gid ||'---'||cast(test as text);
          --执行返回结果
          return query
          select v_uptap_gid as res_uptap_gid,v_uptap_geom as res_uptap_geom ;
        END IF;
        --return next;
      END LOOP;
    END LOOP;
    END IF;
end;  

$BODY$
  LANGUAGE plpgsql VOLATILE STRICT
  COST 100
  ROWS 1000;
ALTER FUNCTION test_getpoint5(character varying, double precision, double precision)
  OWNER TO postgres;

结尾：

这里本想将上游阀门和下游阀门分开，但是我们建立的拓扑中并没有方向，所以改成了查询出所有的影响阀门。

后续将继续研究，把方向数据放进去，实现上游阀门、下游阀门、精确找出总阀门的功能