PipeRapidの情報をQGISに表示させる(2)最長ルートの探索
PythonのNetworkXというパッケージを使って、下水管路の最長ルートを探索することができます。
雨水の流量表には最長延長が記載されていますが、PipeRapidではそのルートがどこなのかを知ることが出来ません。（他の縦断ソフトでは知ることが出来るのかは存じ上げません。画面すら見たことが無いので。）
そもそも最長ルートの位置を知りたいという要望もいままで聞いたことが無く、知って何かの役に立つのかすら私自身わかっていないのですが、個人的に以前から知りたかったので、今回スクリプトを作ってみました。
PythonやNetworkXを知る前は、VBとかでツリーを作って・・・みたいな構想をしていましたが、最長ルートを探索するアルゴリズムが分からなくて諦めていました。
やがてPythonに興味を持つことになり、ツリーを操作できるNetworkXというパッケージがあることを知りました。
下水道特有の「最長ルート」などというキーワードではなんの手がかりも得られなかったので、ひたすらNetworkXのドキュメントとgoogle翻訳で手がかりを探しまくりました。
詳しくは割愛しますが、概要としては、(1)選択したスパンに到達するすべての最上流スパンをピックアップする。(2)ピックアップした最上流スパンから選択したスパンまでの延長をそれぞれすべて計算し、最大延長となるルートを最長ルートとする。といった方法です。
あと実は最も苦労したのは最長ルートの探索ではなくて、探索したスパンをQGIS上で選択する方法でした。
QGIS2のフィーチャー(ここでいうスパンに該当)を選択する方法は日本語で詳しく書いてくださっている方がたくさんいらっしゃいますが、私がメインで使っているQGIS3でフィーチャーを選択する方法が日本語のサイトやブログではとうとう見つけられませんでした。結局、これも外国のFAQ掲示板サイトを翻訳しながら見つけました。ただ動作がどうも遅いので、最良の方法かどうかはわかりません。
さて、普段なら自分の書いたスクリプトはソースコードにパクリが含まれていたり、技量不足の露呈を恐れるあまり公開しないのですが、今回は特定のサイトなどからのパクリが無いため恐る恐る掲載します。なおプログラムの組み方とか命名とかの作法は目に余るほど酷いかと思いますがご容赦ください。
では以下にQGISでの動作画面と、最後に見るに堪えないであろうソースコードを載せておきます。ご利用の際は、あらかじめQGISのPython環境にnetworkxライブラリをインストールしておいてください。その準備だけで動くと思います。
# -*- coding: utf-8 -*-
from qgis.utils import iface
from qgis.core import *
import networkx as nx
# 下流側も追跡する場合は1を入力
endspan_set = 1
pr_graph = nx.DiGraph()
pr_list = [] # 全スパンリストを初期化
pr_top = [] # 起点スパンリストを初期化
# 選択したレイヤのフィーチャーを取得(Linkレイヤを事前選択すること)
layer = iface.activeLayer()
# 全部のフィーチャーを取得
features = []
features = [f for f in layer.getFeatures()]
# 選択したフィーチャーを取得(複数選択した場合は、１本のみ取得)
selectfeatures = []
selectfeatures = [sf for sf in layer.selectedFeatures()]
try:
selectfeature = selectfeatures[0]
print(selectfeature['URouteNo'], 'を選択しました。')
# 取得した全フィーチャーのスパン番号とスパン長を有向グラフに登録
for feature in features:
pr_list.extend([(feature['UNodeNo'], feature['LNodeNo'], \
{'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['UNodeNo'], 'Length')
if targetlength > maxlength:
maxspan = pr_span
maxlength = targetlength
except:
pass
# 選択スパンの延長を加算する
maxlength += float(selectfeature['Length'])
maxroute = nx.bellman_ford_path \
(pr_graph, maxspan, selectfeature['UNodeNo'], 'Length')
for pr_num in maxroute:
# QGISフィーチャーからの検索
exp = QgsFeatureRequest().setFilterExpression("UNodeNo ILIKE \'" + pr_num + '\'')
serchfeatures = layer.getFeatures(exp)
# 検索されたフィーチャーを選択状態にする
for feature in serchfeatures:
layer.selectByIds([feature.id()],QgsVectorLayer.AddToSelection)
print(feature['URouteNo'] + ' ' + feature['Length'])
print('最長延長=' + str(maxlength) + 'm')
# 下流側も追跡する場合の処理
if endspan_set == 1:
# 最下流までのスパン延長リスト
underspanfeatures = nx.single_source_bellman_ford_path_length \
(pr_graph, selectfeature['UNodeNo'], 'Length')
# 最下流までのスパン延長を探索し、上流側最長延長に加算
endspan= max(underspanfeatures, key=underspanfeatures.get)
endroutelength = float(underspanfeatures[endspan])
for pr_span in pr_graph:
try:
allroute = nx.single_source_shortest_path(pr_graph, selectfeature['UNodeNo'])
except:
pass
for pr_num in allroute:
# QGISフィーチャーからの検索
exp = QgsFeatureRequest().setFilterExpression \
("UNodeNo ILIKE \'" + pr_num + '\'')
serchfeatures = layer.getFeatures(exp)
# 検索されたフィーチャーを選択状態にする
feature = None
for feature in serchfeatures:
layer.selectByIds([feature.id()],QgsVectorLayer.AddToSelection)
print(feature['URouteNo'] + ' ' + feature['Length'])
# 選択フィーチャーの長さを減じる(２重形状を無くす)
endroutelength -= float(selectfeature['Length'])
print('下流までの延長=' + str(endroutelength) + 'm')
print('合計延長=' + str(maxlength + endroutelength) + 'm')
except:
#selectfeature = None
print('選択したフィーチャーはありませんでした。')
