OpenStreetMapを使って練馬区で最寄り駅から最も遠い町丁目を求めてみた

はじめに

本記事ではOpen Street Mapを使って、練馬区で最寄り駅から最も遠い町丁目を求めます。
筆者が練馬区民のため対象を練馬区としていますが、コードの一部を変えるだけであらゆる市区町村の最寄り駅から最も遠い町丁目を求めることができます。

今回のコードはヘルシンキ大学が公開しているGISに関する教材を元にしています。
筆者はそちらの教材を勉強した後のアウトプットでこの記事を書いています。
GISについて勉強してみたい方はぜひそちらを使ってみてください。そちらの教材についての質問があればQiitaのコメントで議論させていただけると、筆者も勉強になり大変嬉しいです。

※本記事はGoogle Colab上で使用することを想定しています。ローカル環境でのコード実行時には事前に必要なライブラリのインストールが必要になります。

OpenStreetMapとは

OpenStreetMapのプロジェクトロゴ

OpenStreetMap(OSM)は、道路地図などの地理情報データを誰でも利用できるよう、フリーの地理情報データを作成することを目的としたプロジェクトです。誰もが自由に参加して、自由に編集でき、自由に利用する事ができる、いわば地図のWikipediaのようなものです。

PythonではOSMnxと呼ばれるライブラリでOSMを利用することができ、、道路情報の分析、視覚化を行うことができます。OSMnxでは、レストラン、学校など様々なポイントデータの取得から、２地点の徒歩、自動車等での最短経路を算出するアルゴリズムも実装されています。

練馬区の町丁目ごとに重心を求める

まずは必要なライブラリをインポートします。

!pip install geopandas
!pip install osmnx
import osmnx as ox
import geopandas as gpd
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import os
from shapely.geometry import Point, LineString, Polygon
from shapely.ops import nearest_points
from pyproj import CRS

まずは東京都の町丁目ごとの境界データと人口データをe-Statから取得します。
e-Statは日本の統計が閲覧できる政府統計のポータルサイトです。各省庁が公表する統計データがまとまっており、検索したり地図上に表示できたりします。

境界データには町丁目ごとの名称(丸の内一丁目など)が含まれていないため、人口データも取得してマージすることで、境界データに名称データを加えることにします。

サイトからデータを手動でダウンロードしてくることも可能ですが、今回はPythonライブラリのdownloadを使ってデータをダウンロードします。

# 東京都の行政区域ごとの境界データをダウンロードする
data_fp = 'datasetのフォルダのパスを貼り付けてください'
download.download(
"https://www.e-stat.go.jp/gis/statmap-search/data?dlserveyId=A002005212020&code=13&coordSys=1&format=shape&downloadType=5&datum=2000",
   data_fp,
   kind="zip",
   replace=True,
)

# 東京都の行政区域ごとの人口データをダウンロードする
download.download(
   "https://www.e-stat.go.jp/gis/statmap-search/data?statsId=T001081&code=13&downloadType=2",
   data_fp,
   kind="zip",
   replace=True,
)

ダウンロードした東京都の境界データをGeopandasで読み込みます。

fp_tokyo =  os.path.join(data_fp, 'r2ka13.shp')
geo_tokyo = gpd.read_file(fp_tokyo)

東京都の境界データから、練馬区のデータを抽出します。

geo_tokyo = geo_tokyo.loc[:, ['geometry', 'CITY_NAME', 'KEY_CODE']]
geo_tokyo['KEY_CODE'] = geo_tokyo['KEY_CODE'].astype(str)
geo_nerima = geo_tokyo[geo_tokyo["CITY_NAME"]=="練馬区"]

geo_nerima.head()

練馬区の町丁目ごとの重心を求めます。
Geopandasには.centroidメソッドが実装されており、ポリゴンデータの重心を求めることができます。
重心の単位をメートルで求めるために、一旦epsgを6691に変換してから.centroidを適用して、再度epsgを4326に変換しておきます。

geo_nerima = geo_nerima.to_crs(CRS.from_epsg(6691))
geo_nerima["centroid"] = geo_nerima.centroid
geo_nerima = geo_nerima.set_geometry("centroid")
geo_nerima = geo_nerima.to_crs(CRS.from_epsg(4326))
geo_nerima = geo_nerima.set_geometry("geometry")
geo_nerima = geo_nerima.to_crs(CRS.from_epsg(4326))
geo_nerima.head()

次に人口データをPandasで読み込みます。

fp2 =  os.path.join(data_fp, 'tblT001081C13.txt')
data = pd.read_csv(fp2, encoding="CP932", dtype={'KEY_CODE': str, 'HYOSYO': str})
data.head()

data = data.loc[:, ["KEY_CODE", "HYOSYO", "NAME"]]
data = data[data['HYOSYO'] == "4"]
data['KEY_CODE'] = data['KEY_CODE'].astype(str)
data.dropna(inplace=True)
data.head()

KEY_CODEと町丁目名を持ったデータフレームが取得できました。

境界データと、人口データから作成した町丁目名のデータをマージします。

geodata_nerima =  geo_nerima.merge(data, on='KEY_CODE')
geodata_nerima = geodata_nerima.loc[:,['geometry','CITY_NAME', 'NAME','centroid']]
geodata_nerima = geodata_nerima.to_crs(CRS.from_epsg(6691))
geodata_nerima.reset_index(drop=True, inplace=True)
geodata_nerima.head()

町丁目名、ポリゴンデータ、重心のポイントデータを持ったデータフレームを取得することができました。

町丁目ごとの重心の最近傍nodeを取得する

Open Street Mapでは道路がnode(節)とway(道)で構成されています。nodeは交差点や曲がり角など道路の変化点のイメージです。ある緯度経度間の最短経路を求める際には、ある緯度経度の最近傍node間の経路を求めることになります。そのため、上で求めた町丁目ごとの重心の最近傍nodeを求めます。

まずは練馬区の道路情報を取得します。
ox.geocode_to_gdfのジオコードを使って練馬区のポリゴンデータを取得し、network_typeを指定してそのポリゴン内の道路情報を取得します。
network_typeには以下のようなものがあります。
・drive – 運転可能な公道（ただし、整備道路を除く）
・drive_service – 整備道路を含む運転可能な道路
・walk – 歩行者が使用できるすべての道路と小道
・bike – サイクリストが使用できるすべての道路と小道
・all – すべての私道でないOSMの道路と小道
・all_private – 私道を含むすべてのOSMの道路と小道

place_name = 'Nerima, Tokyo, Japan'
place_polygon = ox.geocode_to_gdf(place_name)
graph = ox.graph_from_polygon(place_polygon["geometry"].values[0], network_type='walk')
fig, ax = ox.plot_graph(graph)

取得できました。

各町丁目の重心から最も近いnodeを取得します。

geodata_nerima["osmid"] = geodata_nerima["centroid"].apply(lambda row:ox.nearest_nodes(graph, row.x, row.y))
geodata_nerima.head()

練馬区周辺の駅のリストを取得する

練馬区周辺の駅のリストを取得します。
練馬区のみでなく周辺の駅を取得しようとしている理由は、練馬区の端の方の町丁目は、最寄り駅が他の市区町村の駅になると思われるからです。

鉄道駅のデータは国土交通省から取得します。令和3年の鉄道駅のデータをダウンロードして、datasetフォルダに格納しておいて、Geopandasで読み込みます。

data_fp = '/content/drive/MyDrive/Colab Notebooks/eo_study/GIS/FinalAssignment/UrbanIndicators/dataset'
fp = os.path.join(data_fp, "N02-21_GML/utf8/N02-21_Station.geojson")
StationData = gpd.read_file(fp)
StationData.head()

取得したデータのgeometryはLINESTRING型で線データになっています。対象駅から下り方面への次の駅までの軌跡の座標データが入っています。今回は駅のPOINTデータが必要なため、それぞれの線の開始点のデータを取得します。

get_point = lambda row:Point(row.coords[0])
StationData["geometry"] = StationData["geometry"].apply(get_point)
StationData.head()

POINTデータを取得することができました。

必要なカラムを取り出して、カラム名を変更します。

StationData = StationData.loc[:, ["N02_005", "geometry"]]
StationData = StationData.rename(columns={'N02_005': 'StationName'})
StationData.head()

CRSをgeodata_nerimaに合わせて、単位をメートルにします。

StationData = StationData.to_crs(geodata_nerima.crs)

練馬区の周囲3kmの範囲のポリゴンを取得します。

dissolved = geodata_nerima.dissolve(by="CITY_NAME")
dissolved['geometry'] = dissolved['geometry'].apply(lambda row:row.buffer(3000))
geo_nerima = dissolved.loc[['練馬区']]
geo_nerima.plot()

先ほどの練馬区の形状と比べて全体的に丸くなっており、bufferをとった形状になっていることが見て取れます。

先ほど作成した練馬区周辺の市区町村の駅のリスト(StationData)のうち、練馬区周囲3kmのポリゴン(geodata_nerima)と重なっているものを抽出します。

Station_point = gpd.sjoin(StationData, geo_nerima, how="inner", op="within")
Station_point.head()

取得した練馬区周辺の駅のリストと練馬区周辺のポリゴンを表示します。

fig, ax = plt.subplots(figsize=(15,8))
geo_nerima.plot(ax=ax)
Station_point.plot(ax=ax, color='black', markersize=20)

ポリゴンと重なっている駅を抽出できたことがわかります。

各駅の最近傍nodeを取得する

上で町丁目ごとの重心の最近傍nodeを取得したのと同様に、各駅の最近傍nodeを取得します。
駅のリストは練馬区の周囲3kmの範囲と重なるものを抽出したため、道路情報も練馬区の周囲3kmの範囲で取得します。

place_name = "Nerima, Tokyo, Japan"
place_polygon = ox.geocode_to_gdf(place_name)
place_polygon = place_polygon.to_crs(epsg=6691)
place_polygon["geometry"] = place_polygon.buffer(3000)
place_polygon = place_polygon.to_crs(epsg=4326)
graph = ox.graph_from_polygon(place_polygon["geometry"].values[0], network_type='walk')
fig, ax = ox.plot_graph(graph)

def get_nearest_nodes(row, graph):
  osmid = ox.nearest_nodes(graph, row.x, row.y)
  return osmid

上と同様に最近傍nodeを取得します。
Station_point["osmid"] = Station_point["geometry"].apply(lambda row:get_nearest_nodes(row, graph))
Station_point.reset_index(drop=True, inplace=True)
Station_point.head()

取得することができました。

各町丁目からの最寄り駅と最短経路を求める

今回の主題である、各町丁目からの最寄り駅と最短経路を求めていきます。
上で作成した練馬区の周囲3kmの道路のnode情報をGeopandasに変換します。
最短経路はメートルで求めたいため、CRSをepsg=6691に変換しておきます。

nodes = ox.graph_to_gdfs(graph, nodes=True, edges=False)
nodes = nodes.to_crs(epsg=6691)
nodes.head()

各町丁目ごとに各駅までの最短経路を取得し、最短経路の距離が最も短い駅を最寄り駅として、"町丁目名"、"最寄り駅名"、"最短経路のgeometry"、"最短経路の距離"をデータフレーム(routes)に格納していきます。
※総当たりで最短経路を求めているので、計算時間が長いので注意してください。

columns = ['Name', 'nearest_Station', 'geometry', 'route_length']
routes = gpd.GeoDataFrame(columns=columns, crs=nodes.crs)

for i in range(len(geodata_nerima)):
  min_length = 1e9
  for j in range(len(Station_point)):
    origin_osmid = geodata_nerima.loc[i]["osmid"]
    dest_osmid = Station_point.loc[j]["osmid"]
    if origin_osmid == dest_osmid:
      continue
    route = ox.shortest_path(graph, origin_osmid, dest_osmid, weight='length')
    if route == None:
      continue

    route_nodes = nodes.loc[route]
    route_line = LineString(list(route_nodes.geometry.values))
    # print(geodata_nerima.loc[i]["NAME"], Station_point.loc[j]["StationName"], route_line.length)
    if min_length > route_line.length:
      min_route = route_line
      min_length = route_line.length
      nearest_Station = Station_point.loc[j]["StationName"]
  routes = routes.append({'Name': geodata_nerima.loc[i]["NAME"], 'nearest_Station': nearest_Station, 'geometry': min_route, 'route_length':min_length}, ignore_index=True)
routes

各町丁目の最寄り駅と最短経路が求められました。

最寄り駅から最も遠い町丁目を求めるため、最短経路の距離(route_length)が最大の行を取得します。

index_max = routes["route_length"].idxmax()
routes.loc[index_max]

下記の結果が得られました。
練馬区で最寄り駅から最も遠い町丁目は「大泉学園町八丁目」で最寄り駅は和光市駅、距離は2.915kmと算出されました。
確認のため、GoogleMapで大泉学園町八丁目から和光市駅までの経路検索したところ、2.9kmでした。かなり近い値が得られているようです。

まとめ

今回はOpen Street Mapを使用して、練馬区で最寄り駅から最も遠い町丁目を求めてみました。
扱うデータを増やして東京23区の”陸の孤島”を見つけるつもりでしたが、ポリゴンデータの取得の手間と計算時間の観点から今回は練馬区に絞って検証してみました。

ポリゴンデータの取得を効率化したり、各町丁目から各駅の最短経路の取得を総当たりでなく近傍の何点かに絞ったり、工夫を凝らして東京23区での検証も実現したいと思っています。（余裕がある時に…

これからも勉強しながら成果を記事にしていこうと思うので、よろしくお願いします。