PipeRapidの情報をQGISに表示させる(3)選択したスパンに流入する全スパンの探索 2022/11/27修正

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

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

前回投稿した「PipeRapidの情報をQGISに表示させる(2)最長ルートの探索」のスクリプトを応用して、選択スパンに流入する上流側の全スパンを探索することが出来ます。

最長ルートの探索はnetworkx.bellman_ford_pathを使いましたが、上流側の全スパンを探索するにはnetworkx.single_source_shortest_pathを使います。なおグラフも最長ルート探索のときとは逆順で登録する必要があります。ソースコードをご参照ください(説明が手抜きですみません)。

また、関数の詳細については、networkxのドキュメントに記載されています。https://networkx.github.io/documentation/stable/reference/algorithms/shortest_paths.html」

最長ルート探索との違いはソースコードだけです。準備すべきモジュールやQGISの操作に関しては変わりありませんので、冒頭に示した過去記事をご参照ください。

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 -*-

from pkg_resources import working_set
from qgis.utils import iface
from qgis.core import *
import networkx as nx
import datetime

# 分水を考慮しない場合は Bunsui変数に 0 を入力し、
# 考慮する場合は 1 を入力する。
# 下記注記の方法によりLinkレイヤに仮想フィールドを設定する。
Bunsui = 0  # 通常は 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 = []  # 起点スパンリストを初期化

all_in_node = []  # 抽出スパンのノードリスト
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_L], feature[B_field_U], \
        {'fid':feature.id(),'RouteNo': feature['RouteNo'], 'Length': float(feature['Length'])})])
pr_graph.add_edges_from(pr_list)

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

try:
	# networkxの出力結果(辞書型)からkeyを取り出し、List(配列)に変換
    all_in_node = list((nx.single_source_shortest_path(pr_graph, selectfeature[B_field_U])).keys())
except:
    pass

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))}  # 全データを辞書型に変換
# 上記において 左から二番目のpr_listの第二添字を0(LManholeNo) から 1(UManholeNo_B) に入れ替えると、分水スパンを１スパン含めて抽出できる
# これはrapidtest_maxlength.py(最長延長抽出)とは数字が逆であることに注意

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

layer.select(fid_list)

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


##### デバッグ用CSV出力コード #####
'''
import csv
myList = pr_dict.items()
with open('pr_dict.csv','w',encoding='utf-8',newline='') as csvfile:
    writer = csv.writer(csvfile)
    writer.writerows(myList)
'''

print('終了')

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

Follow me!

@GDWING_SW