この記事では、QGISとPythonを組み合わせて巡回セールスマン問題(TSP)をLGWAN環境で解決する方法をご紹介します！

「え？インターネット使えないLGWAN環境で？」と思ったそこのあなた！そうなんです、セキュリティが厳しいLGWAN環境では、外部のAPIやサービスが使えず、頭を抱えることも多いですよね…。「GRASS GISのv.salesmanがあるじゃないか」という声も聞こえてきそうですが、GRASS GISを使うには、プラグインのインストール、ワークスペースの作成、データ変換など、意外と面倒な手順が必要なんです😩

そこで！今回は、もっと手軽に、QGISのPythonコンソールだけで実行できるコードでTSPを解く方法をお届けします！これなら、GRASS GISのような複雑な準備は一切不要！QGISの標準機能で完結します！✨

巡回セールスマン問題って何？🤔

そもそもTSPって何？という方のために簡単に説明すると、複数の地点を回るとき、最短距離で効率よく巡回するにはどうしたらいいか？という問題です。訪問先が増えるほどルートの組み合わせが爆発的に増えて、最適なルートを見つけるのが本当に大変なんです🤯

QGISでTSPを解く手順

必要なデータ

巡回する拠点の緯度経度情報（例：エクセルファイル）
道路ネットワークデータ（OpenStreetMapなどから取得し、QGISで読み込める形式に変換）

QGISでの操作

道路ネットワークデータをQGISに読み込みます。
巡回拠点の緯度経度をQGISでポイントデータとして作成します。

Pythonコードの実行

以下のPythonコードをQGISのPythonコンソールにコピー＆ペーストして実行！

# 必要なライブラリをインポート
from qgis.core import (
    QgsProject,
    QgsVectorLayer,
    QgsPointXY,
    QgsFeature,
    QgsGeometry,
    QgsFields,
    QgsField,
    QgsCoordinateReferenceSystem,
    QgsCoordinateTransform,
    QgsDistanceArea
)
from qgis.PyQt.QtCore import QVariant
import networkx as nx

# 1. 開始拠点と巡回拠点のデータ (緯度経度で入力)
#   サンプルデータは愛媛県の観光地を巡るルートです！
start = {'name': '内子座', 'lat': 33.5394, 'lon': 132.6853}
patrol = [
    {'name': '臥龍山荘', 'lat': 33.5031, 'lon': 132.5467},
    {'name': '大洲城', 'lat': 33.5006, 'lon': 132.5461},
    {'name': '宇和島城', 'lat': 33.2231, 'lon': 132.5606},
]

# 座標系変換の設定 (緯度経度を平面直角座標系に変換)
source_crs = QgsCoordinateReferenceSystem("EPSG:4326")  # 緯度経度
target_crs = QgsCoordinateReferenceSystem("EPSG:6672")  # 平面直角座標系
transform_context = QgsProject.instance().transformContext()
coordinate_transform = QgsCoordinateTransform(source_crs, target_crs, transform_context)

# 緯度経度を平面直角座標系のx, yに変換する関数
def transform_coordinates(lat, lon):
    point = QgsPointXY(lon, lat)
    transformed_point = coordinate_transform.transform(point)
    return transformed_point.x(), transformed_point.y()

# 開始拠点と巡回拠点の座標を変換
start['x'], start['y'] = transform_coordinates(start['lat'], start['lon'])
for p in patrol:
    p['x'], p['y'] = transform_coordinates(p['lat'], p['lon'])

# 3. 施設の座標をリストとして取得
locations = [(start['x'], start['y'])] + [(p['x'], p['y']) for p in patrol]

# 4. 道路ネットワークレイヤー名
layer_name = "gis_osm_roads_free_1_6672" # ★★★★★ ご自身の環境に合わせて書き換えてください！★★★★★

# 5. レイヤーを取得
road_layer = QgsProject.instance().mapLayersByName(layer_name)
if not road_layer:
    raise ValueError(f"レイヤー '{layer_name}' が見つかりません。プロジェクトに正しいレイヤーを追加してください。")
road_layer = road_layer[0]

# 6. 距離計算用のクラスを定義
distance_calculator = QgsDistanceArea()
distance_calculator.setSourceCrs(source_crs, transform_context)

# 7. 道路ネットワークをNetworkXグラフに変換
G = nx.Graph()
for feature in road_layer.getFeatures():
    geom = feature.geometry()
    lines = geom.asMultiPolyline() if geom.isMultipart() else [geom.asPolyline()]
    for line in lines:
        for i in range(len(line) - 1):
            point1, point2 = line[i], line[i + 1]
            
            # 距離計算（EPSG:6672）
            p1 = QgsPointXY(point1)
            p2 = QgsPointXY(point2)
            dist = distance_calculator.measureLine(p1, p2)
            G.add_edge(
                (p1.x(), p1.y()),
                (p2.x(), p2.y()),
                weight=dist
            )

# 8. 最寄りノードを探す関数
def find_nearest_node(x, y):
    point = QgsPointXY(x, y)
    nearest_node = min(
        G.nodes,
        key=lambda n: QgsPointXY(n[0], n[1]).distance(point)
    )
    return nearest_node

# 9. 各拠点の最寄りノードを取得
nodes = [find_nearest_node(x, y) for x, y in locations]

# 10. ノード間の距離行列を作成
n = len(nodes)
dist_matrix = [
    [
        nx.shortest_path_length(G, source=nodes[i], target=nodes[j], weight='weight') for j in range(n)
    ]
    for i in range(n)
]

# 11. 貪欲法でルートを計算
def greedy_tsp(start_node, dist_matrix):
    visited = [False] * n
    route = [start_node]
    current_node = start_node
    visited[current_node] = True
    total_distance = 0

    for _ in range(n - 1):
        next_node = None
        min_distance = float('inf')
        for i in range(n):
            if not visited[i] and dist_matrix[current_node][i] < min_distance:
                min_distance = dist_matrix[current_node][i]
                next_node = i
        route.append(next_node)
        visited[next_node] = True
        total_distance += min_distance
        current_node = next_node

    # 最後にスタート地点に戻る
    route.append(start_node)
    total_distance += dist_matrix[current_node][start_node]
    return route, total_distance

# 貪欲法でルートを計算
start_node = 0  
best_route, total_distance = greedy_tsp(start_node, dist_matrix)

# 結果を表示
location_names = [start['name']] + [p['name'] for p in patrol]
print("\n貪欲法によるルート:")
for i in range(len(best_route) - 1):
    start_index = best_route[i]
    end_index = best_route[i + 1]
    distance_km = dist_matrix[start_index][end_index] / 1000  # 距離をkmに変換
    print(f"{location_names[start_index]} -> {location_names[end_index]} ({distance_km:.2f} km)")

# 総距離を表示
total_distance_km = total_distance / 1000
print(f"総距離: {total_distance_km:.2f} km")

# 12. 結果をQGISに描画
route_layer = QgsVectorLayer("LineString?crs=EPSG:6672", "Greedy TSP Route", "memory")
route_layer_data = route_layer.dataProvider()

fields = QgsFields()
fields.append(QgsField("ID", QVariant.Int))
route_layer_data.addAttributes(fields)
route_layer.updateFields()

for i in range(len(best_route) - 1):
    start_node = nodes[best_route[i]]
    end_node = nodes[best_route[i + 1]]
    path = nx.shortest_path(G, source=start_node, target=end_node, weight='weight')
    coords = [QgsPointXY(n[0], n[1]) for n in path]
    feature = QgsFeature()
    feature.setGeometry(QgsGeometry.fromPolylineXY(coords))
    feature.setAttributes([i])
    route_layer_data.addFeature(feature)

QgsProject.instance().addMapLayer(route_layer)

コード解説

コードは長いですが、一つずつ見ていけば大丈夫！👌大きく分けて4つの処理を行っています。

データの準備: 巡回拠点と道路ネットワークデータを準備します。
- サンプルデータでは、愛媛県の観光地(内子座、臥龍山荘、大洲城、宇和島城)を巡るケースを想定しています。皆さんの業務に合わせて変更してくださいね！
- 道路ネットワークは、QGISで読み込める形式（例：シェープファイル）で用意しましょう。
前処理:
- 入力された緯度経度を、平面直角座標系（ここではEPSG:6672）に変換します。QGISで正確な距離を計算するには、この変換が必要なんです！もし他のエリアで試される場合はその地域にあったEPSGを選んでください。
- NetworkXライブラリを使って、道路ネットワークをグラフ構造に変換します。交差点を「ノード」、道路を「エッジ」として、ネットワーク構造を表現するイメージです。
経路探索: 今回は貪欲法というアルゴリズムを使って、最短経路を探します。「今いる地点から一番近い拠点へ移動する」を繰り返す、シンプルで高速な方法です！💨
結果の出力: 最適化された巡回ルートをQGISのレイヤーに表示します。地図上で結果を確認できるので、とっても便利！🗺️
貪欲法の限界と、その先へ…
貪欲法は、シンプルで高速な反面、「局所最適解」、つまり必ずしも全体として最も良い解が得られるわけではない、という弱点があります。拠点数が多い場合や、道路ネットワークが複雑な場合は、最適な解とのズレが大きくなる可能性も。
「もっと完璧な解が欲しい！」という方は、遺伝的アルゴリズムや焼きなまし法などの「メタヒューリスティクス」と呼ばれる高度な手法に挑戦するのもアリです！これらの手法は、貪欲法よりも計算時間は長くなりますが、より良い解に辿り着ける可能性があります。

活用例：あなたの業務を、もっとスマートに！

福祉・介護: 訪問介護のルート最適化で、移動時間を短縮し、利用者さん一人ひとりにより長く向き合う時間を創出！💖
災害対策:
- 避難所巡回ルートの事前計画: 災害発生時、限られた時間内に、どの避難所を回るべきか、効率的なルートを迅速に導き出すことは非常に重要です。QGISとPythonで、最適な避難所巡回ルートを事前に計画し、シミュレーションすることで、迅速かつ確実な避難所運営と被災者支援を実現！🏃‍♂️💨
- 物資搬送ルートの事前調査: 災害発生時、道路状況は刻々と変化します。被災状況や道路の寸断情報をQGISにリアルタイムに反映することで、物資搬送の「最適ルート」を素早く見つけ出し、被災者へ迅速かつ確実に支援物資を届けることが可能に！🚚💨
ゴミ収集: 収集ルートを最適化し、コスト削減と環境負荷軽減に貢献！🌱
観光: 観光スポットを効率よく巡るルートを提案し、観光客の満足度UP！✨
営業・配達: 配送ルートの最適化で、ガソリン代節約＆CO2排出量削減！🌍

まとめ

QGISとPythonの強力タッグは、LGWAN環境でもオフラインでTSPを解くための強力な武器になります！今回紹介した貪欲法は、処理が速く、簡単に実装できるのが大きなメリット。ただし、完璧な解が得られない場合があることも頭の片隅に置いておいてくださいね。
さあ、あなたもQGIS×Pythonで、スマートな自治体業務を実現しましょう！

最後に

この記事が、LGWAN環境で奮闘する皆さんの業務効率化に少しでも貢献できれば嬉しいです！
🚀 地図で未来を切り拓きましょう！🚀