PipeRapidの情報をQGISに表示させる(2)最長ルートの探索 2022/11/27修正

以前に掲載しておりましたコードの不具合を見直したことにより、探索速度を大幅に改善しました。
なお分水を考慮するか否かを選択するために、コード内のパラメータ変更に加え、QGISのフィールドの設定も必要となりました。手間が増えますがご承知おきください。

また、分水地点で２方向から流入があった場合に、最長パスではない路線を１スパン選択してしまう不具合が解消しきれていませんので、ご注意ください。暫定的な対処方法はソースコード内に注記しましたので、必要に応じてご対応ください。

～以下、操作方法の説明です～

PythonのNetworkXというパッケージを使って、下水管路の最長ルートを探索することができます。

雨水の流量表には最長延長が記載されていますが、PipeRapidではそのルートがどこなのかを知ることが出来ません。（他の縦断ソフトでは知ることが出来るのかは存じ上げません。画面すら見たことが無いので。）

そもそも最長ルートの位置を知りたいという要望もいままで聞いたことが無く、知って何かの役に立つのかすら私自身わかっていないのですが、個人的に以前から知りたかったので、今回スクリプトを作ってみました。

PythonやNetworkXを知る前は、VBとかでツリーを作って・・・みたいな構想をしていましたが、最長ルートを探索するアルゴリズムが分からなくて諦めていました。

やがてPythonに興味を持つことになり、ツリーを操作できるNetworkXというパッケージがあることを知りました。

「最長ルート」などというキーワードではなんの手がかりも得られなかったので、ひたすらNetworkXのドキュメントとgoogle翻訳で手がかりを探しまくりました。

詳しくは割愛しますが、概要としては、(1)選択したスパンに到達するすべての最上流スパンをピックアップする。(2)ピックアップした最上流スパンから選択したスパンまでの延長をそれぞれすべて計算し、最大延長となるルートを最長ルートとする。といった方法です。

あと実は最も苦労したのは最長ルートの探索ではなくて、探索したスパンをQGIS上で選択する方法でした。

QGIS2のフィーチャー(ここでいうスパンに該当)を選択する方法は日本語で詳しく書いてくださっている方がたくさんいらっしゃいますが、私がメインで使っているQGIS3でフィーチャーを選択する方法が日本語のサイトやブログではとうとう見つけられませんでした。結局、これも外国のFAQ掲示板サイトを翻訳しながら見つけました。ただ動作がどうも遅いので、最良の方法かどうかはわかりません。

さて、普段なら自分の書いたスクリプトはソースコードにパクリが含まれていたり、技量不足の露呈を恐れるあまり公開しないのですが、今回は特定のサイトなどからのパクリが無いため恐る恐る掲載します。なおプログラムの組み方とか命名とかの作法は目に余るほど酷いかと思いますがご容赦ください。

では以下にQGISでの動作画面と、最後に見るに堪えないであろうソースコードを載せておきます。~~ご利用の際は、あらかじめQGISのPython環境にnetworkxライブラリをインストールしておいてください。その準備だけで動くと思います。~~（11/27 訂正：QGISのPythonには最初からインストールされてるっぽいので、何もしなくても良いかも知れません）

11/27 追記：QGISにおいて、Linkレイヤに仮想フィールドを追加しておいてください。

フィールド名：UMNo_B　，　計算式：if(right( “UManholeNo” ,2)=’_1′ , substr(“UManholeNo”,1,-2),”UManholeNo”)
フィールド名：LMNo_B　，　計算式：if(right( “LManholeNo” ,2)=’_1′ , substr(“LManholeNo”,1,-2),”LManholeNo”)

# -*- coding: utf-8 -*-

#################################################################
# 【超重要な不具合】
# 分水地点で２方向から流入があった場合に、最長パスではない路線を１スパン
# 選択してしまう不具合がある。現状は、下記の注意※1のパラーメータを1として
# 分水１スパンも抽出し、後で取り除くこと。
#################################################################

from qgis.utils import iface
from qgis.core import *
import networkx as nx

import datetime

# 下流側も追跡する場合は1を入力
endspan_set = 0

# 分水を考慮しない場合は Bunsui変数に 0 を入力し、
# 考慮する場合は 1 を入力する。
# 下記注記の方法によりLinkレイヤに仮想フィールドを設定する。
Bunsui = 0  # 通常は 1 が望ましい。その場合は下記注記により予め仮想フィールドを設定すること。
if Bunsui == 0:
    B_field_U = 'UNodeNo'
    B_field_L = 'LNodeNo'
else:
    B_field_U = 'UMNo_B'
    B_field_L = 'LMNo_B'
    #B_field_L = 'LManholeNo'

'''
※重要：上流追跡に分水を考慮する場合の仮想フィールド設定について

Linkレイヤに「UMNo_B」「LMNo_B」という名前の仮想フィールドを作成し、
そのフィールド計算式にはそれぞれ下記のように入力する。
if((right("UManholeNo" ,2)='_1') or (right("UManholeNo" ,2)='_2'), substr("UManholeNo",1,-2),"UManholeNo")
if((right("LManholeNo" ,2)='_1') or (right("LManholeNo" ,2)='_2'), substr("LManholeNo",1,-2),"LManholeNo")
'''

pr_graph = nx.DiGraph()
pr_list = []  # 全スパンリストを初期化
pr_top = []  # 起点スパンリストを初期化

fid_list = []  # 抽出スパンのIDリスト

# 選択したレイヤのフィーチャーを取得(Linkレイヤを事前選択すること)
layer = iface.activeLayer()

# 全部のフィーチャーを取得
features = []
# features = [f for f in layer.getFeatures()]  # 下の方法で出来るため、この行は廃止した
features = layer.getFeatures()

# 選択したフィーチャーを取得(複数選択した場合は、１本のみ取得)
selectfeatures = []

print(__name__)

#selectfeatures = [sf for sf in layer.selectedFeatures()]  # 下の方法で出来るため、この行は廃止した
selectfeatures = layer.selectedFeatures()
selectfeature = selectfeatures[0]
print(selectfeature['RouteNo'], 'を選択しました。')

# 取得した全フィーチャーのスパン番号とスパン長を有向グラフに登録
time_0 = datetime.datetime.now()
for feature in features:
    # LengthがNULLの場合は0を代入する
    if feature['Length'] == None:
        feature['Length'] = 0
    
    pr_list.extend([(feature[B_field_U], feature[B_field_L], \
        {'fid':feature.id(),'RouteNo': feature['RouteNo'], 'Length': float(feature['Length'])})])
pr_graph.add_edges_from(pr_list)

# 最上流スパンの抽出
for pt in pr_graph:
    indeg = pr_graph.in_degree(pt)
    if indeg == 0:
        pr_top.append(pt)

maxlength = 0
for pr_span in pr_top:
    try:
        targetlength = nx.bellman_ford_path_length \
            (pr_graph, pr_span, selectfeature[B_field_U], 'Length')
        if targetlength > maxlength:
            maxspan = pr_span
            maxlength = targetlength
    except:
        pass

# 選択スパンの延長を加算する
maxlength += float(selectfeature['Length'])

time_1 = datetime.datetime.now()
print ('グラフ登録時間 ' + str(time_1 - time_0))

maxroute = nx.bellman_ford_path \
    (pr_graph, maxspan, selectfeature[B_field_U], 'Length')

time_2 = datetime.datetime.now()
print ('上流追跡時間 ' + str(time_2 - time_1) + ' Total ' + str(time_2 - time_0))

pr_dict = {pr_list[i][2]['fid'] : pr_list[i][0] for i in range(0, len(pr_list))}  # 全データを辞書型に変換
# 注意※1
# 上記において 左から二番目のpr_listの第二添字を1(LManholeNo) から 0(UManholeNo_B) に入れ替えると、分水スパンを１スパン含めて抽出できる
# これはrapidtest_all_in_area.py(全流入スパン抽出)とは数字が逆であることに注意

fid_list.extend(k for k, v in pr_dict.items() if v in maxroute)

layer.select(fid_list)

print('最長延長=' + str(maxlength) + 'm')

endroutelength = 0
# 下流側も追跡する場合の処理
if endspan_set == 1:
    # 最下流までのスパン延長リスト
    lowerspanfeatures = nx.single_source_bellman_ford_path_length \
        (pr_graph, selectfeature[B_field_U], 'Length')
    
    # 最下流までのスパン延長を探索し、上流側最長延長に加算
    endspan= max(lowerspanfeatures, key=lowerspanfeatures.get)
    endroutelength = float(lowerspanfeatures[endspan])
        
    for pr_span in pr_graph:
        try:
            lowerroute = nx.single_source_shortest_path(pr_graph, selectfeature[B_field_U])
        except:
            pass

    pr_dict = {pr_list[i][2]['fid'] : pr_list[i][0] for i in range(0, len(pr_list))}  # 全データを辞書型に変換
    # 上記において 左から二番目のpr_listの第二添字を0(LManholeNo) から 1(UManholeNo_B) に入れ替えると、分水スパンを１スパン含めて抽出できる
    # これはrapidtest_all_in_area.py(全流入スパン抽出)とは数字が逆であることに注意
    
    fid_list.extend(k for k, v in pr_dict.items() if v in lowerroute)
    print('下流までの延長=' + str(endroutelength) + 'm')
    
print('合計延長=' + str(maxlength + endroutelength) + 'm')
layer.select(fid_list)

time_3 = datetime.datetime.now()
print ('地物選択時間 ' + str(time_3 - time_2) + ' Total ' + str(time_3 - time_0))

print('終了')

#except:
    #selectfeature = None
    #print('選択したフィーチャーはありませんでした。')

Follow me!

@GDWING_SW