QGISで一番遠いところに行こう

これは MIERUNE AdventCalendar 2024 16日目の記事です! 昨日は@Takayuki_Kawajiriさんによる記事でした。

はじめに

ネットワークの経路探索では最短経路を探すのが最も定番なトピックの一つで、それに向けてQGISのプロセッシングツールとして「Shortest path(point to point)」というツールも用意されています。効率の良い最短経路以外に、より長い距離を網羅するために、最も遠い距離や最も離れているポイントも実は現実世界では必要とされています。例えば、観光客になるべく長いルートで途中の景色を楽しめるような経路設計、なるべく多くの場所を経由する物流ルートや公共交通ルートの設計や、水分分析で下流点から上流側の最も長い谷線を探索することで主要な水流経路を特定するなど、様々な場面で応用されている。

この記事ではQGISでスクリプトの実行によって簡単に始発点から最も遠いポイントへの経路を探す方法を紹介します。

動作環境

• QGIS 3.34
• Mac OS

データ

• 国土地理院 - 数値地図（国土基本情報）- 道路中心線データ
  • https://github.com/gsi-cyberjapan/optimal_bvmap?tab=readme-ov-file#pmtiles%E7%89%88
  • PMTiles地物の入手はtippecanoe-decodeを利用います。本題から離れますので、詳細手順は割愛します。
• サンプルとして利用する札幌駅北口の座標値はQGISから取得しています。

スクリプト

ネットワークの理論計算ができるPythonのパッケージnetworkxを利用します。QGIS3.34の実行環境にデフォルト状態でnetworkxが内包されており、スクリプトで直接importして利用できます。処理の流れですが、ジオメトリの読み込みはgeopandasで行い、networkxの解析を経由し、QGISで描画できるQgsVectorLayerに変換するという形になります。

まずはGeoDataFrameの入出力ができる関数を作成します。

def search_path_to_farthest_node(
    start_point: gpd.GeoDataFrame, road: gpd.GeoDataFrame
) -> gpd.GeoDataFrame:

道路網をグラフへの変換するために、LineStringからネットワークのedgeに変換します。

    G = nx.Graph()
    for line in road.geometry:
        if isinstance(line, (LineString, MultiLineString)):
            lines = [line] if isinstance(line, LineString) else list(line.geoms)
            for segment_line in lines:
                coords = list(segment_line.coords)
                edges = [(coords[i], coords[i + 1]) for i in range(len(coords) - 1)]
                G.add_edges_from(edges, weight=segment_line.length)

始発点はネットワーク上のノードに限らないので、ポイントから最近接ノードを特定する関数を用意し、経路探索の起点とします。

    def find_nearest_node(point, graph):
        nearest_node = None
        min_distance = float("inf")
        for node in graph.nodes:
            distance = Point(node).distance(point)
            if distance < min_distance:
                min_distance = distance
                nearest_node = node
        return nearest_node

最も遠いノードへの経路を探索します。

  if start_node:
      # Dijkstra法を使用し、起点から全ノードへの最短経路長を計算する
      lengths = nx.single_source_dijkstra_path_length(G, start_node)
      
      # 最も遠いノード（farthest_node）を、計算結果（lengths）の中から特定する
      farthest_node = max(lengths, key=lengths.get)
      # edgeの重みとして距離（weight）を考慮し、start_node から farthest_node への最短経路を確定する
      path_to_farthest_node = nx.shortest_path(
          G, source=start_node, target=farthest_node, weight="weight"
      )
			# LineStringに変換する
      path_line = LineString(path_to_farthest_node)

最後にGeoDataFrameに格納して返します。

    return gpd.GeoDataFrame(geometry=[path_line])

ここまで最も遠いノードへの経路を特定する関数を作成しました。

このようなスクリプトで実行します。戻り値のGeoDataFrameをjson形式に変換し、そのままQgsVectorLayerとして作成できます。

from qgis.core import QgsVectorLayer, QgsProject

road = gpd.read_file(
    "道路データの絶対パス"
)
start_point = gpd.GeoDataFrame(geometry=[Point(141.3510449, 43.0675969)])  # 札幌駅
path_to_farthest_node = search_path_to_farthest_node(start_point, road)

path_to_farthest_node_str = path_to_farthest_node.to_json()
path_layer = QgsVectorLayer(path_to_farthest_node_str, "path_to_farthest_node", "ogr")
QgsProject.instance().addMapLayer(path_layer)

処理結果

札幌市中心部の道で、札幌駅北口からの最も遠い場所を計算してみたら、駅前通り〜大通公園〜薄野〜市電通りを経由して中島公園の出口（豊平川沿い）へ向かうルートが生成されました。賑やかな都市景観から豊かな自然まで、札幌中心部を十分楽しめる観光ルートになったのではないでしょうか！

さらに北海道全体（離島を除く）を全て含めてみたら…なんと道南の「北海道道607号石崎松前線」の突き当たりになりました。

最後に

繰り返しになりますが、まとまったソースコードも載せておきます。ぜひQGISで試してみてください〜

import geopandas as gpd
import networkx as nx
from shapely.geometry import LineString, MultiLineString, Point

def search_path_to_farthest_node(
    start_point: gpd.GeoDataFrame, road: gpd.GeoDataFrame
) -> gpd.GeoDataFrame:
    # 道路データをNetworkXのグラフに変換
    G = nx.Graph()
    for line in road.geometry:
        if isinstance(line, (LineString, MultiLineString)):
            lines = [line] if isinstance(line, LineString) else list(line.geoms)
            for segment_line in lines:
                coords = list(segment_line.coords)
                edges = [(coords[i], coords[i + 1]) for i in range(len(coords) - 1)]
                G.add_edges_from(edges, weight=segment_line.length)

    # start-pointに最も近いノードを見つける
    def find_nearest_node(point, graph):
        nearest_node = None
        min_distance = float("inf")
        for node in graph.nodes:
            distance = Point(node).distance(point)
            if distance < min_distance:
                min_distance = distance
                nearest_node = node
        return nearest_node

    start_node = find_nearest_node(start_point.geometry[0], G)

    # 最長経路を探索
    if start_node:
        # NetworkXの関数を使って最長経路を探索
        lengths = nx.single_source_dijkstra_path_length(G, start_node)
        farthest_node = max(lengths, key=lengths.get)
        path_to_farthest_node = nx.shortest_path(
            G, source=start_node, target=farthest_node, weight="weight"
        )

        # 最長経路をGeoJSONに変換
        path_line = LineString(path_to_farthest_node)

    else:
        path_line = LineString()
    return gpd.GeoDataFrame(geometry=[path_line])
    

from qgis.core import QgsVectorLayer, QgsProject

road = gpd.read_file(
     "道路データの絶対パス"
)
start_point = gpd.GeoDataFrame(geometry=[Point(141.3510449, 43.0675969)])  # 札幌駅
longest_path = search_path_to_farthest_node(start_point, road)

longestpath_str = longest_path.to_json()
longestpath_layer = QgsVectorLayer(longestpath_str, "longest_path", "ogr")
QgsProject.instance().addMapLayer(longestpath_layer)

明日は@kntoshiyaさんによる記事です！お楽しみに！！