PostGIS 爆管分析之找出上游閥門

環境：

Win10

ArcMap10.4（用於數據處理）

postgresql9.4

postgis2.2.3

pgRouting2.3（postgresql插件）

說明：

繼上一篇文章作了爆管分析找出周圍全部影響閥門後（參見：https://www.cnblogs.com/giser-s/p/11662932.html），發如今業務上使用有侷限性，由於一般爆管之後我並不須要關閉全部周圍閥門，而是隻要關閉上游閥門便可。

下面的方法，是在查找到周圍全部閥門的基礎上繼續的，在周圍閥門中找出與他相接的上游閥門

步驟：

 一、在PostGIS中，創建方向模型，模擬氣體流向。原理就是新增length或者reverse_cost（關鍵字），與source/target方向一致的，則length爲值，reverse_cost爲設一個大點的值9999999；如source/target方向不一致，則給reverse_cost賦值，length爲設一個大點的值9999999。

這裏在查詢時，pgr_dijkstraCost會識別關鍵詞cost和reverse_cost

 

 

 雙向查詢時不用關注reverse_cost（注意cost字段必定要有，或能夠以下例子將length命名cost的別名）

SELECT * FROM pgr_dijkstra( 'SELECT gid AS id, source, target, length as cost  FROM zy', 15139, 13670, directed := false);

單向查詢時須要加上reverse_cost一塊兒查（注意cost字段必定要有，或能夠以下例子將length命名cost的別名）

SELECT * FROM pgr_dijkstra( 'SELECT gid AS id, source, target, length as cost,reverse_cost as reverse_cost  FROM zy', 15139, 13670, directed := true);

 

 

 

 

 

#正向結果：

#逆向結果：

 

 

 

 

二、這裏查詢到爆點周圍全部影響閥門後，須要進行記錄，這裏採用拼接成字符串記錄在數組中（這裏略複雜，不一樣於上一博文直接返回輸出了）

格式：[{v_uptap_gid||','||cast(v_uptap_geom as text)||','||up_temprow.source]}]

說明：[{閥門gid，轉成字符串的geom，管段的頭source}]

raise notice '%' , up_temprow.source;
          --記錄閥門所在管段source
          IF (v_cost @> ARRAY[up_temprow.source::integer]) THEN
          ELSE
            SELECT array_append(v_cost,up_temprow.source) into v_cost;
            SELECT array_append(v_cost,up_temprow.target) into v_cost;
          END IF;
          IF (res_tap_pipe is not null) THEN
            select res_tap_pipe || ARRAY[v_uptap_gid||','||cast(v_uptap_geom as text)||','||up_temprow.source] into res_tap_pipe;
          ELSE
            select ARRAY[v_uptap_gid||','||cast(v_uptap_geom as text)||','||up_temprow.source] into res_tap_pipe;
          END IF;

3、v_cost是周圍全部閥門，利用pgRouting的pgr_dijkstraCost函數，逆向找出與爆點管段相接的閥門。

    --查找上游閥門
    FOREACH m_cost IN ARRAY v_cost 
    LOOP
      SELECT count(*) FROM pgr_dijkstraCost('select gid as id, source, target, length as cost, reverse_cost from zy',m_cost, ARRAY[v_startSource,v_startTarget], true) where agg_cost >= 9999999 into m_cost_value;
      raise notice '%' , 'm_cost---'||cast(m_cost as text);
      ----若是沒有消耗大於9999999的（閾值），則認爲是上游閥門
      IF(m_cost_value = 0) THEN          
        FOREACH m_tap_pipe IN ARRAY res_tap_pipe 
        LOOP
          raise notice '%' , 'm_cost---'||cast(m_cost as text)  ;
          raise notice '%' , 'm_tap_pipe---'||cast(m_tap_pipe as text)  ;
          IF (split_part(m_tap_pipe, ',', 3)::integer = m_cost) THEN 
            --閥門id，閥門圖形要素，閥門類型（上游/下游）
            return query
            select split_part(m_tap_pipe, ',', 1)::integer as res_uptap_gid,split_part(m_tap_pipe, ',', 2)::geometry  as res_uptap_geom ,split_part(m_tap_pipe, ',', 3)::integer as res_source;
          END IF;
        END LOOP;
      END IF;
    END LOOP;

四、附上所有存儲過程

  

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排除
    v_down_where integer[];--where條件，將遍歷到閥門的管段gid排除
    up_temprow record ;
    --v_cost record;--記錄閥門管段source（用於計算消耗，判斷方向）
    m_cost integer;
    m_cost_value integer;
    temprow record;
    v_cost integer[];
    res_source integer;
    res_tap_pipe text[];
    m_tap_pipe text;
    idx_tap_pipe integer; --遍歷結果遊標
    m_up_cost integer;--上游閥門
    v_up_cost integer[];--上游閥門集合
    res_main_pipe integer[];--總閥門集合
    m_main_pipe 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; 
    raise notice '%' , 'v_startSource---'||cast(v_startSource as text);
    IF(v_startLine is not null) THEN
    --查找上游閥門
    v_up_idx = 0;
    v_up_source = 1;
    --尋找上游閥門
    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; 
      --記錄層級
      --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) 
        --select zy1.gid,zy1.source,zy1.target from zy zy1 where target = any(v_up_where)--找上游
      LOOP
        --清空須要查的點
        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;  
        raise notice '%' , 'UP---'||up_temprow.gid;
        --若是沒查找到閥門，則繼續往下查
        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 '%' , up_temprow.source;
          --記錄閥門所在管段source
          IF (v_cost @> ARRAY[up_temprow.source::integer]) THEN
          ELSE
            SELECT array_append(v_cost,up_temprow.source) into v_cost;
            SELECT array_append(v_cost,up_temprow.target) into v_cost;
          END IF;
          IF (res_tap_pipe is not null) THEN
            select res_tap_pipe || ARRAY[v_uptap_gid||','||cast(v_uptap_geom as text)||','||up_temprow.source] into res_tap_pipe;
          ELSE
            select ARRAY[v_uptap_gid||','||cast(v_uptap_geom as text)||','||up_temprow.source] into res_tap_pipe;
          END IF;
        END IF;
      END LOOP;
    END LOOP;
    --raise notice '%' , v_cost;
    raise notice '%' , 'res_tap_pipe---'||cast(res_tap_pipe as text);
    --return query select * from v_cost;
    raise notice '%' , 'v_cost---'||cast(v_cost as text);

    --查找上游閥門
    FOREACH m_cost IN ARRAY v_cost 
    LOOP
      SELECT count(*) FROM pgr_dijkstraCost('select gid as id, source, target, length as cost, reverse_cost from zy',m_cost, ARRAY[v_startSource,v_startTarget], true) where agg_cost >= 9999999 into m_cost_value;
      raise notice '%' , 'm_cost---'||cast(m_cost as text);
      ----若是沒有消耗大於9999999的（閾值），則認爲是上游閥門
      --IF(m_cost_value = 0) THEN
      --  --判斷上游閥門間是否有上下游關係
      --  SELECT array_append(v_up_cost,m_cost) into v_up_cost;
      --END IF;
      IF(m_cost_value = 0) THEN          
        FOREACH m_tap_pipe IN ARRAY res_tap_pipe 
        LOOP
          raise notice '%' , 'm_cost---'||cast(m_cost as text)  ;
          raise notice '%' , 'm_tap_pipe---'||cast(m_tap_pipe as text)  ;
          IF (split_part(m_tap_pipe, ',', 3)::integer = m_cost) THEN 
            --閥門id，閥門圖形要素，閥門類型（上游/下游）
            return query
            select split_part(m_tap_pipe, ',', 1)::integer as res_uptap_gid,split_part(m_tap_pipe, ',', 2)::geometry  as res_uptap_geom ,split_part(m_tap_pipe, ',', 3)::integer as res_source;
          END IF;
        END LOOP;
      END IF;
    END LOOP;

     raise notice '%' , '上游閥門---'||cast(v_up_cost as text);
  END IF;
end;  

結尾：

此文算是半成品文章，代碼很亂也沒有進行優化，不少沒有用到的變量沒有刪掉，這裏暫且看成思路的記錄。

算法上還有不少須要改進，這裏只是在找到全部周邊閥門的基礎上繼續往下寫的，其實能夠不用分兩塊，直接上來就開始找上游閥門（算法留待後面繼續優化）