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;

結尾：

這里本想將上游閥門和下游閥門分開，但是我們建立的拓撲中並沒有方向，所以改成了查詢出所有的影響閥門。

后續將繼續研究，把方向數據放進去，實現上游閥門、下游閥門、精確找出總閥門的功能