これは「MIERUNE Advent Calendar 2023」の16日目の記事です。昨日は、@xinmiao1995 さんによる「QGISでデフォルメマップを作る」でした。
はじめに
Pythonには、OSMnxというOpenStreetMapのデータを利用して、道路ネットワークやその他の地物を抽出し、分析・可視化することができるライブラリがあります。
基本的にはOpenStreetMapのデータを使用することが前提のライブラリのため、その他の道路ネットワークデータを使用して、ネットワーク解析をすることは少し厄介です。
そのため、今回はOSMnxも使用しつつ、OSMnxの基盤となっているNetworkXを使用して、OpenStreetMap以外の道路データを使用して、ネットワーク解析する方法を紹介します。
OSMNXとは
OSMnx（OpenStreetMap and NetworkX）は、Pythonで作成されたソフトウェアパッケージであり、地理空間データと街路網を簡単に取得、構築、分析、可視化するためのツールです。このライブラリは特にOpenStreetMap（OSM）からデータを取得するために設計されています。OSMnxはNetworkXと連携して、複雑なネットワーク解析を行うことができます。
NetworkXとは
NetworkXは、複雑なネットワークの構築、操作、研究を行うためのPythonライブラリです。グラフ理論に基づいたデータ構造と分析ツールを提供し、さまざまな種類のネットワークの研究などで利用されます。
使用データの準備
解析に必要なデータをQGISを使用して取得します。
以下の手順で取得したデータは、GeoJSON形式で平面直角座標系にエクスポートして、カレントディレクトリのdataフォルダ内に格納しておくことを前提に説明します。
道路データ
今回は、地理院ベクトルタイルの道路中心データを使用するため、「GSI-VTDownloader」プラグインを使用して、札幌市中心部のデータを入手します。
このデータはマルチラインとして出力されるため、プロセシングの「マルチパートをシングルパートに変換」を使用して、ラインデータに変換しておきます。
データの名称は、road.geojsonとします。
ポイントデータ
ネットワーク解析の起点となるポイントのデータを入手します。
今回は、QuickOSMプラグインを使用して、札幌市中心部にあるコンビニのデータを取得します。
データの名称は、shop.geojsonとします。
パッケージのインストール
データが揃ったので、ここからOSMnxを使用してネットワーク解析をするための準備をしていきます。
Poetrtyを使って環境構築します。
poetry add geopandas
poetry add shapely
poetry add osmnx
poetry add networkx
poetry add scipy
poetry add matplotlib
パッケージのインポート
import geopandas as gpd
from shapely.geometry import Point
import osmnx as ox
import networkx as nx
データのインポート
まずは、道路データからネットワーク解析のための環境を構築する必要があるため、GeoPandasのジオデータフレームとしてインポートします。
road_gdf = gpd.read_file("data/road.geojson")
shop_gdf = gpd.read_file("data/shop.geojson")
グラフの生成
次に、インポートした道路データを使用して、ネットワーク解析で必要となる「グラフ」を作成していきます。
OSMnxにおけるグラフ（Graph）は、主に道路ネットワークをモデル化したもので、このグラフはNetworkXのグラフオブジェクトを基にしており、ノードとエッジで構成されています。
- ノード: 交差点や道路の終端などを表します。ノードは特定の地理的座標に関連付けられています。
- エッジ: 道路の一部分（交差点から交差点までなど）を表します。エッジは、必要に応じて距離、速度制限、道路の種類（例：主要道、裏道など）といった属性を持ちます。
エッジからノードの生成
OSMNXでのグラフの生成には、ノードとエッジが必要となるため、道路のラインデータから所定の形式のエッジとノードを生成します。
# わずかな座標の違いで同じ場所に複数のポイントが作成されるためポイントの座標を丸める
def round_point(point, precision=7):
return Point(round(point.x, precision), round(point.y, precision))
nodes_dict = {}
edges_list = []
id_node = 0
id_edge = 0
for index, row in road_gdf.iterrows():
line = row.geometry
start_point = round_point(Point(line.coords[0]))
end_point = round_point(Point(line.coords[-1]))
if start_point not in nodes_dict:
nodes_dict[start_point] = id_node
id_node += 1
u = nodes_dict[start_point]
if end_point not in nodes_dict:
nodes_dict[end_point] = id_node
id_node += 1
v = nodes_dict[end_point]
edge_data = {"id": id_edge, "u": u, "v": v, "geometry": line}
edges_list.append(edge_data)
id_edge += 1
nodes_gdf = gpd.GeoDataFrame({
"id": list(nodes_dict.values()),
"x": [point.x for point in nodes_dict.keys()],
"y": [point.y for point in nodes_dict.keys()],
"geometry": list(nodes_dict.keys())
})
edges_gdf = gpd.GeoDataFrame(edges_list).set_index(["u", "v", "id"])
edges_gdf['length'] = edges_gdf['geometry'].length
作成したエッジとノードのジオデータフレームを確認してみましょう。
-
edges_gdf
-
nodes_gdf
グラフを生成
では、上記で作成したデータをもとに、OSMnxを使用してグラフデータを作成します。
なお、今回は道路の辺の向きを考慮せずにネットワーク解析をしたいので、G.to_undirectedで無向グラフに変換をしています。
# 札幌市のデータを使用しているため、平面直角座標系12系を指定
G = ox.utils_graph.graph_from_gdfs(nodes_gdf, edges_gdf)
G.graph['crs'] = "EPSG:6680"
# 有向グラフを無向グラフに変換
G_undirected = G.to_undirected()
グラフのプロット
作成されたグラフをプロットしてみましょう。
ox.plot_graph(G_undirected)
エッジとノードをGISデータに出力
この手順は、ネットワーク解析において必ずしも必要な手順ではないですが、グラフのデータをGISデータとして出力するには、以下の手順で行うことが可能です。
エッジ
# エッジをGeojsonで出力
graph_edges_gdf = ox.graph_to_gdfs(G_undirected, nodes=False, edges=True)
graph_edges_gdf.set_crs(epsg=6680, inplace=True)
graph_edges_gdf.to_file("output/edges.geojson", driver="GeoJSON")
ノード
自分の環境ではox.graph_to_gdfsでノードを出力できなかったので、グラフからジオデータフレームを作成して、GISデータとして出力します。
x = []
y = []
nodeid = []
geometry = []
for node, data in list(G_undirected.nodes(data=True)):
if data:
x.append(data["x"])
y.append(data["y"])
nodeid.append(node)
geometry.append(Point(data["x"], data["y"]))
crs = "EPSG:6680"
graph_nodes_gdf = gpd.GeoDataFrame({
"y": y,
"x": x,
"nodeid": nodeid,
"geometry": geometry
},crs=crs)
graph_nodes_gdf.to_file("output/nodes.geojson", driver="GeoJSON")
上記のデータをQGISで表示してみます。
最短経路検索
これで準備ができたので、ネットワーク解析を行なっていきます。
始点と終点のノードを選定
まずは、shop_gdfから適当に視点と終点となるポイントを選定します。
今回は、始点を札幌駅南東付近にあるセブンイレブン、終点を円山公園駅の近くにあるローソンとしてみます。（なお選定に関して特に意味はありません。）
-
始点
-
終点
手順としては、それぞれのosm_idをキーとしてジオメトリを取得し、ox.nearest_nodesで最近傍のノードを取得します。
始点の最近傍のノードを取得
start_geom = shop_gdf[shop_gdf['osm_id'] == "5422615621"]
geometry = start_geom.geometry
coords = (geometry.x, geometry.y)
start_node = ox.nearest_nodes(G_undirected, coords[0], coords[1])
終点の最近傍のノードを取得
end_point = shop_gdf[shop_gdf['osm_id'] == "3390462405"]
end_geom = end_point.geometry
coords = (end_geom.x, end_geom.y)
end_node = ox.distance.nearest_nodes(G_undirected, coords[0], coords[1])
最短経路検索
上記で取得したノードを元に最短経路探索をして、結果をプロットしてみます。結果をみると、最短経路が取得できていることが確認できます。
route = nx.shortest_path(G_undirected, source=start_node[0], target=end_node[0], weight='length')
ox.plot_graph_route(G_undirected, route))
GeoJSONで出力
上記の結果をGISデータとして出力して、QGISで表示してみます。
route_line = ox.utils_graph.graph_to_gdfs(G_undirected.subgraph(route), edges=True, nodes=False)
route_line = route_line.geometry.unary_union
route_gdf = gpd.GeoDataFrame(geometry=[route_line], crs="EPSG:6680")
route_gdf.to_file("output/route.geojson", driver="GeoJSON")
到達圏解析
最後に到達圏解析も行ってみます。
到達圏作成の起点となる場所は、サッポロファクトリーにあるセブンイレブンにしてみます。
始点の最近傍のノードを取得する
start_geom = shop_gdf[shop_gdf['osm_id'] == "4964974720"]
geometry = start_geom.geometry
coords = (geometry.x, geometry.y)
start_node = ox.nearest_nodes(G_undirected, coords[0], coords[1])
到達圏を求める
nx.ego_graphを使用して、到達圏を求めます。
meters = 500 # 到達圏の範囲(m)
subgraph = nx.ego_graph(G_undirected, start_node[0], radius=meters, distance='length') # 単方向
結果のプロット
nc = ['r' if node == start_node else 'b' for node in subgraph.nodes()]
ns = [100 if node == start_node else 8 for node in subgraph.nodes()]
fig, ax = ox.plot_graph(subgraph, node_size=ns, node_color=nc, node_zorder=2)
出典
- 国土地理院ベクトルタイル提供実験
- 地理院タイル