地図からこだわりの道を抽出する方法
この記事は、FOSS4G Advent Calendar 2023 の25日目の記事です
想定する読者
この記事では、オープンデータからこだわりの道を抽出する方法について解説します。
これは、GIS初心者の私が手を動かした時の事を元にしています。
はじめてオープンデータに触れた時に、何から触れば良いのかわからなかったので、同じ境遇の人の参考になればうれしいです。
使用するデータ
地図データには無料のOSM(OpenStreetMap)を使います。
OSM内には車道、歩道、山、川、店舗等が含まれていて、これらを使って目的の道を抽出します。
OSMはポケモンGOやTeslaのカーナビにも使われていたりします。
抽出プロセス
実際に手を動かす前に、ざっくりと抽出プロセスについて解説します。
処理を単純化するために2段構成になっています。
候補を抽出
OSMには膨大なデータが含まれているため、一気にすべての事を行うと大変です。
なので、まずゆるい条件で絞り込み、目的の道の候補を抽出します。
| 絞り込みなし | 絞り込みあり |
|---|---|
 |
 |
この絞り込んだエッジ（点と点の間のライン）が候補となります。
エッジは特定区間の道のような物です。
特徴を評価して絞り込む
さきほどのエッジ使ってさらに絞り込みを行います。
エッジは複数の点（緯度、経度）の集まりからできています。
この点を評価する評価関数を作り、エッジを評価します。
評価値の上位数件がこだわりの道となりなります。
手を動かして実装する
こだわりの道を「らくして運転できる道」と定義し、抽出できるまでをゴールとします。
らくして運転できる道の特徴はこんな感じにしました。
- 一般道
- 2車線以上の道
- 800m以上の区間
- 分岐が少ない道
- 直線的な道
小さい範囲で地図を取得する
まずは「楽して運転できる道」が含まれる5kmの範囲の地図データを取得します。
この地図から「楽して運転できる道」を抽出するまでが第一ゴールです。
地図データ（OSM）の取得にはOSMnxというライブラリを使いました。
このライブラリはOSMのデータ取得、解析、視覚化をするためのライブラリです。
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]')
いい感じに抽出できたら完成です。
完成したコードはこちらです。
応用
これまでの処理をベースにして、評価関数を「頭文字Dで走るような峠道」に変更したものがこちらになります。
osm以外に国土地理院さんの標高モデルも使っています。
赤色が強い所ほど「頭文字Dで走るような峠道」の特徴が強いです。
その他使えそうなデータ
おわりに
素人が素人なりの言葉でまとめてみました。
素人が手を動かして実装したのでおかしい所があるかもしれません。
どこかの誰かの参考になったら幸いです。
Discussion