🗺️

地図からこだわりの道を抽出する方法

2023/12/25に公開

Python

この記事は、FOSS4G Advent Calendar 2023 の25日目の記事です

想定する読者

この記事では、オープンデータからこだわりの道を抽出する方法について解説します。
これは、GIS初心者の私が手を動かした時の事を元にしています。
はじめてオープンデータに触れた時に、何から触れば良いのかわからなかったので、同じ境遇の人の参考になればうれしいです。

使用するデータ

地図データには無料のOSM(OpenStreetMap)を使います。
OSM内には車道、歩道、山、川、店舗等が含まれていて、これらを使って目的の道を抽出します。
https://openstreetmap.jp/
OSMはポケモンGOやTeslaのカーナビにも使われていたりします。

抽出プロセス

実際に手を動かす前に、ざっくりと抽出プロセスについて解説します。
処理を単純化するために2段構成になっています。

候補を抽出

OSMには膨大なデータが含まれているため、一気にすべての事を行うと大変です。
なので、まずゆるい条件で絞り込み、目的の道の候補を抽出します。

絞り込みなし	絞り込みあり

この絞り込んだエッジ（点と点の間のライン）が候補となります。
エッジは特定区間の道のような物です。
エッジの一覧

特徴を評価して絞り込む

さきほどのエッジ使ってさらに絞り込みを行います。
エッジは複数の点（緯度、経度）の集まりからできています。

この点を評価する評価関数を作り、エッジを評価します。
評価値の上位数件がこだわりの道となりなります。

手を動かして実装する

こだわりの道を「らくして運転できる道」と定義し、抽出できるまでをゴールとします。
らくして運転できる道の特徴はこんな感じにしました。

一般道
2車線以上の道
800m以上の区間
分岐が少ない道
直線的な道

小さい範囲で地図を取得する

まずは「楽して運転できる道」が含まれる5kmの範囲の地図データを取得します。
この地図から「楽して運転できる道」を抽出するまでが第一ゴールです。

地図データ（OSM）の取得にはOSMnxというライブラリを使いました。
このライブラリはOSMのデータ取得、解析、視覚化をするためのライブラリです。
https://osmnx.readthedocs.io/en/stable/

OSMの道には以下のような道路種別と車線の数の情報を含んでいます。

OSMnxはOSMデータを取得する時に条件を指定できるので、抽出条件に「一般道」と「２車線以上の道」を含ませ、ここで確定させます。

import osmnx as ox
import networkx as nx
latitude = 35.336452
longitude = 136.922696
dist = 5000
# キャッシュを使う
ox.config(use_cache=True, log_console=True)
# 5km以内の道路を取得する
graph = ox.graph_from_point(center_point=(latitude,longitude),
                            network_type='drive',
                           simplify=True,
                           retain_all=True,
                           dist=dist,
                           custom_filter='["highway"~"secondary|secondary_link|primary|primary_link|trunk|trunk_link"]["lanes"=2]')
graph2 = ox.graph_from_point(center_point=(latitude,longitude),
                            network_type='drive',
                           simplify=True,
                           retain_all=True,
                           dist=dist,
                           custom_filter='["highway"~"tertiary"]')
# highway:tertiaryは2車線を表すのでlanesの抽出は不要。
# その他の道は2車線を表すタグがついている道のみを抽出する。
# highway: unclassifiedは狭い道なので含めない。
# 複雑なcustom_filter指定方法が分からなかったので、抽出処理を分割して結合してます。
graph = nx.compose(graph, graph2)
# JupyterNotebook上に地図を表示する
map = ox.plot_graph_folium(graph, edge_width=2)
map

800m以下の道を取り除く

短すぎるエッジに価値はないのでここで取り除きます。
Osmnxで取得したデータはgraph形式です。
この形式はエッジ単体での分析には適していません。
これをGeoDataFrame形式に変換すると、エッジ毎の解析がとても行いやすくなります。

# graphからGeoDataFrameに変換
gdf_edges = ox.graph_to_gdfs(graph, nodes=False, edges=True)
# 800m以下の道を取り除く
gdf_edges = gdf_edges[gdf_edges["length"] >= 800]

ここではlengthを指定しましたが、他にも色々なタグがあったりします。

逆方向の道を取り除く

抽出したエッジには、一方通行の道を除き、進行方向と逆方向の2つのエッジが含まれています。
今回は片側だけでよいので取り除きます。

# 逆方向のエッジを削除する
drop_target = []
for index, row in gdf_edges.iterrows():
    if (index[1], index[0], 0) in drop_target:
        continue
    if (index[0], index[1], 0) in drop_target:
        continue
    drop_target.append(index)
gdf_edges = gdf_edges[gdf_edges.index.isin(drop_target)]

分岐が少ない道を評価する関数

さきほど抽出したエッジには、エッジから分岐する道の情報を持っていません。
そのため、エッジを構成する点の座標と一致するノードを探して、分岐数を計算する必要があります。
直前に取得した地図データは、車線数等で絞り込まれているため、道の情報がたりません。
なので、OSmnxですべての道を取得し、その道とエッジを比較して分岐数を計算します。

# 絞り込みなしで地図データを取得する
graph_all = ox.graph_from_point(center_point=(latitude,longitude),
                            network_type='drive',
                            simplify=True,
                            retain_all=True,
                            dist=5000)
# ノードを取り出す。
all_nodes = ox.graph_to_gdfs(graph_all, nodes=True, edges=False)
# エッジの分岐数をカウント。
gdf_edges['branch_cnt'] = 0
for index, row in gdf_edges.iterrows():
  nodes = all_nodes[all_nodes.geometry.intersects(row.geometry)]
  # 分岐数を計算
  gdf_edges.at[index, 'branch_cnt'] = nodes['street_count'].sum() - (len(nodes) * 2)

距離と分岐数の比率を求めて、この値を分岐が少ない道の評価値とします。

gdf_edges['branch_cnt_to_length_rate'] = 1 - (gdf_edges['branch_cnt'] / gdf_edges['length'])

直線的な道を評価する関数

エッジには形を構成する座標を持っています。
この座標間の角度の合計を求め、エッジが湾曲しているかどうかを判別します。

# 3座標間の角度を計算
def calculate_angle_between_vectors(A, B, C):
    vector_AB = np.array(B) - np.array(A)
    vector_BC = np.array(C) - np.array(B)

    dot_product = np.dot(vector_AB, vector_BC)
    norm_AB = np.linalg.norm(vector_AB)
    norm_BC = np.linalg.norm(vector_BC)

    cosine_theta = dot_product / (norm_AB * norm_BC)
    angle_rad = np.arccos(cosine_theta)
    angle_deg = np.degrees(angle_rad)
    return angle_deg

# 座標間の角度の変化の合計値を求める
gdf_edges['geometory_angle_total'] = gdf_edges['geometry'].apply(
    lambda x: sum([calculate_angle_between_vectors(x.coords[i-1], x.coords[i], x.coords[i+1]) for i in range(1, len(x.coords)-1)])
)

距離と角度の比率を求めて、この値を直線的な道の評価値とします。

# 角度の合計と距離の比率を求める
gdf_edges['geometory_angle_to_length_rate'] = 1 - (gdf_edges['geometory_angle_total'] / gdf_edges['length'])

「らくして運転できる道」の評価値は「分岐数が少ない道」と「直線的な道」の評価値との積で計算できます。

# 楽して運転できる道の評価値
gdf_edges['score'] = gdf_edges['branch_cnt_to_length_rate'] * gdf_edges['geometory_angle_to_length_rate']
# 0~1の範囲に正規化
gdf_edges['score_normalized'] = gdf_edges['score'] / gdf_edges['score'].max()
# 並び替え
gdf_edges = gdf_edges.sort_values('score_normalized', ascending=False)
# ランキング
gdf_edges['rank'] = gdf_edges['score_normalized'].rank(ascending=False)

グラフィカルに確認する

抽出した内容は数値なので、直感的に確認ができません。
地図上で確認するためにFoliumを使い、上位10件を強調して表示させます。

# GeoDataFrameからGraph形式に変換
graph = ox.graph_from_gdfs(gdf_nodes, gdf_edges)
map_osm = ox.plot_graph_folium(graph, edge_width=2)
# 上位10件を強調して表示
map_osm.add_child(
    folium.features.GeoJson(
        gdf_edges.head(10).to_json(),
        style_function=lambda x: {
            'color': "#FF0000",
            'weight': 7,
            'opacity': (1 - (x['properties']['rank'] - 1) * 0.09)
        }
    )
)

さらに直感的に確認するために、Google Mapsやストリートビュー等のリンクを表示させます。

+ # google_map_urlを作成
+ gdf_edges['google_map_url'] = gdf_edges['geometry'].apply(
+     lambda x: f"https://www.google.com/maps/dir/{x.coords[0][1]},{x.coords[0][0]}/'{x.coords[-1][1]},{x.coords[-1][0]}'"
+ )
+ # street_view_urlを作成
+ gdf_edges['street_view_url'] = gdf_edges['geometry'].apply(
+     lambda x: f"https://www.google.com/maps/@{x.coords[0][1]},{x.coords[0][0]},20?layer=c&cbll={x.coords[-1][1]},{x.coords[-1][0]}&cbp=12,0,0,0,0"
+ )
+ # googe_earth_urlを作成
+ def generate_google_earth_url(row):
+     center_index = math.floor(len(row.geometry.coords) / 2) - 1
+     center = row.geometry.coords[center_index]
+     return f"https://earth.google.com/web/search/{center[1]},+{center[0]}"
+ gdf_edges['google_earth_url'] = gdf_edges.apply(generate_google_earth_url, axis=1)

graph = ox.graph_from_gdfs(gdf_nodes, gdf_edges)
# 地図を表示する
map_osm = ox.plot_graph_folium(graph, edge_width=2)

# 候補の上位10件を表示
map_osm.add_child(
    folium.features.GeoJson(
        gdf_edges.head(10).to_json(),
        style_function=lambda x: {
            'color': "#FF0000",
            'weight': 7,
            'opacity': (1 - (x['properties']['rank'] - 1) * 0.09)
        }
    )
+     # エッジにポップアップを表示する
+     .add_child(folium.features.GeoJsonPopup(
+         fields=['geometory_angle_total', 'length', 'branch_cnt', 'score_normalized', 'google_map_url', 'street_view_url', 'google_earth_url'],
+         aliases=['geometory_angle_total', 'length', 'branch_cnt', 'score_normalized', 'google_map_url', 'street_view_url', 'google_earth_url'],
+         localize=True
+     ))
)

これで、抽出と確認が行えるようになりました。

未知のこだわりの道を抽出する

ここまでは、小さい範囲から「楽して運転できる道」の抽出を行いました。
次は、範囲を広げて「未知の楽して運転できる道」を抽出します。

latitude = 35.336452
longitude = 136.922696
+ dist = 10000
- dist = 5000
graph = ox.graph_from_point(center_point=(latitude,longitude),
                            network_type='drive',
                           simplify=True,
                           retain_all=True,
                           dist=dist,
                           custom_filter='["highway"~"secondary|secondary_link|primary|primary_link|trunk|trunk_link"]["lanes"=2]')