平面交差の考慮と計算コストの節約のために、OpenStreetMapから東北地方の道路リンクを交差点基準で分解する
/
平面交差の考慮と計算コストの節約のために、OpenStreetMapから東北地方の道路リンクを交差点基準で分解します。
# ファイル構成
project_dir
├── /data
│ ├── /R2_boundary/04_miyagi/04_miyagi.shp
│ ├── /osm/tohoku/gis_osm_roads_free_1.shp
│ └── (省略)
└── merge_bound.ipynb <- 実行用ノートブック
# ライブラリのインストール
!pip install keplergl
!pip install geopandas
!pip install sodapy
!jupyter labextension install @jupyter-widgets/jupyterlab-manager keplergl-jupyter
!pip install geopy
!pip install simpledbf
# OSMの読み込み
import pandas as pd
from simpledbf import Dbf5
import geopandas as gpd
from keplergl import KeplerGl
import base64
from geopy.distance import geodesic
import hashlib
# # DBFの読み込み
# dbf = Dbf5('"./data/osm/tohoku/gis_osm_roads_free_1.dbf')
# df = dbf.to_dataframe()
# shpファイル読込
gdf = gpd.read_file("./data/osm/tohoku/gis_osm_roads_free_1.shp")
gdf.head()
| osm_id | code | fclass | name | ref | oneway | maxspeed | layer | bridge | tunnel | geometry |
|---|---|---|---|---|---|---|---|---|---|---|
| 20264905 | 5111 | motorway | 東北自動車道 | E4 | F | 100 | 0 | F | F | LINESTRING (141.09415 38.80402, 141.09454 38.8... |
| 22762217 | 5121 | unclassified | None | None | B | 0 | 0 | F | F | LINESTRING (140.88676 38.30354, 140.88661 38.3... |
| 22762218 | 5122 | residential | None | None | B | 0 | 0 | F | F | LINESTRING (140.88614 38.30239, 140.88638 38.3... |
# 道路種別の選定
OSM fclass (opens new window)を参考に、道路種別の選定します。
- footway：歩道
- living_street：生活道路
- motorway:自専道
- motorway_link:自専道(ランプ)
- pedestrian:歩道
- primary:主要地方道
- residential:居住地域内道路
- secondary:都道府県道
- service:私道
- steps:階段
- tertiary:市町村道
- tertiary_link:市町村道(ランプ)
- track:農道・林道,
- track_grade1:1級_農道・林道,
- track_grade4:4級_農道・林道,
- track_grade5:5級_農道・林道,
- trunk:国道
- trunk_link:国道(ランプ)
# 道路酒別を抽出
gdf2 = gdf[gdf["fclass"].isin(["motorway", "motorway_link",
"primary", "secondary",
"tertiary", "tertiary_link",
"living_street",
"trunk", "trunk_link"])]
# 0~1000行を抽出(短時間で確認したいので)
gdf2 = gdf2[0:1000].reset_index(drop=True)
gdf2.head(2)
# Keplerで確認
# データの読み込み
map1 = KeplerGl(height=400)
map1.add_data(data=gdf2, name='data1')
# htmlで表記
orig_html = str(map1._repr_html_(),'utf-8')
b64d_html = base64.b64encode(orig_html.encode('utf-8')).decode('utf-8')
framed_html = f'<iframe src="data:text/html;base64,{b64d_html}" style="width:95%; height: 600px">'
import IPython
IPython.display.HTML(framed_html)
# 座標の抽出とカウント
各座標が各リンクで何回出現するかをカウントします。
# LINESTRINGを整理
geometry = gdf2["geometry"].astype(str).apply(lambda x: x.replace("LINESTRING (", "").replace(")", ""))
# 座標リストに整形
lonlats = geometry.apply(lambda x: x.split(", "))
# osm_id, fclass, init_lon, init_lat, term_lon, term_lat, length
latlon_list = []
for i,lonlat in enumerate(lonlats):
# print(gdf2["osm_id"][i], gdf2["fclass"][i])
for j in range(len(lonlat)):
point = lonlat[j].split(" ")
# リンクの格納
latlon_list.append([point[1], point[0], 1])
# データフレームの作成
df_latlons = pd.DataFrame(latlon_list, columns=["lat","lon", "count"])
# 緯度経度ごとにカウント
df_latlons = df_latlons.groupby(["lat","lon"]).count().reset_index()
df_latlons.head()
| lat | lon | count |
|---|---|---|
| 36.8940682 | 140.7881092 | 1 |
| 36.8941365 | 140.7879983 | 1 |
| 36.894186 | 140.7878967 | 1 |
| 36.8942832 | 140.7876019 | 1 |
| 36.8946206 | 140.7863557 | 1 |
# 平面交差基準でリンクを分解
下記のような分割で、平面交差基準でリンクを分解します。
data ─┬── ノードが2つ ── そのまま格納
└── ノードが2つ以上 ─┬─ 中間点で交差ノードがない ── そのまま格納
└── 中間点で交差ノードがある ── 交差ノードで分解
# 交差ノードの抽出
# 交差ノードの抽出とリスト化
cross_latlons = df_latlons[df_latlons["count"]>1].reset_index(drop=True)
cross_list = [f"{cross_latlons['lon'][i]} {cross_latlons['lat'][i]}" for i in range(len(cross_latlons))]
# LINESTRINGの分解
# LINESTRINGを整理
geometry = gdf2["geometry"].astype(str).apply(lambda x: x.replace("LINESTRING (", "").replace(")", ""))
# 座標リストに整形
lonlats = geometry.apply(lambda x: x.split(", "))
# リンクの分解
def append_link_list(geometry):
# 始点終点の座標
_prev = geometry[0].split(" ")
_next = geometry[-1].split(" ")
# linestringの整形
points = str(geometry).replace("'","").replace("[","").replace("]","")
# links_listに追加
links_list.append([gdf["osm_id"][i], gdf["fclass"][i],
_prev[1], _prev[0], _next[1], _next[0], 1,
f"LINESTRING ({points})"])
# リンクの分解
# osm_id, fclass, init_lon, init_lat, term_lon, term_lat,
# length , raw, geometry
links_list = []
f = 0
for i,lonlat in enumerate(lonlats):
f += 1
# 2点間リンク or 中間点で交差箇所があるかないか
if len(set(lonlat[1:-1])&set(cross_list))==0:
_prev = lonlat[0].split(" ")
_next = lonlat[-1].split(" ")
links_list.append([gdf2["osm_id"][i], gdf2["fclass"][i],
_prev[1], _prev[0], _next[1], _next[0], 0,
gdf2["geometry"].astype(str)[i]])
# 中間点で交差箇所がある
else:
# geometrys = []
geometry = [lonlat[0]]
for point in lonlat[1:-1]:
# 逐次追加
geometry.append(point)
# 交差地点か否か
if point in cross_list:
# geometrys.append(geometry) # 追加
append_link_list(geometry) # 追加
# 交差箇所から再スタート
geometry = [point]
# 終点
geometry.append(lonlat[-1])
append_link_list(geometry) # 追加
# geometrys.append(geometry)
# データフレームの作成
osm_links = pd.DataFrame(links_list, columns=["osm_id", "fclass", "init_lat", "init_lon",
"term_lat", "term_lon", "raw", "geometry"])
osm_links.head(2)
| osm_id | fclass | init_lat | init_lon | term_lat | term_lon | raw | geometry |
|---|---|---|---|---|---|---|---|
| 20264905 | motorway | 38.8040228 | 141.0941549 | 38.8378983 | 141.0989473 | 0 | LINESTRING (141.0941549 38.8040228, 141.094536... |
| 22787759 | primary | 38.2722395 | 140.8690382 | 38.2688936 | 140.8700114 | 0 | LINESTRING (140.8690382 38.2722395, 140.869062... |
# リンク長の算出
# リンク長の算出
def link_length(latlons):
dis = 0
for i in range(len(latlons)-1):
_prev = latlons[i].split(" ")
_next = latlons[i+1].split(" ")
# 逐次加算
dis += geodesic((_prev[1], _prev[0]), (_next[1], _next[0])).km
return dis
def make_lengths(osm_links):
# LINESTRINGを整理
geometry = osm_links["geometry"].astype(str).apply(lambda x: x.replace("LINESTRING (", "").replace(")", ""))
# 座標リストに整形
latlons = geometry.apply(lambda x: x.split(", "))
# リンク長の算出
lengths = latlons.apply(lambda x: link_length(x))
return lengths
# リンク長の算出
osm_links["length"] = make_lengths(osm_links)
osm_links.head(2)
| osm_id | fclass | init_lat | init_lon | term_lat | term_lon | raw | geometry | length | | - | - | - | - | - | - | - | - | | 20264905 | motorway | 38.8040228 | 141.0941549 | 38.8378983 | 141.0989473 | 0 | LINESTRING (141.0941549 38.8040228, 141.094536... | 3.944538 | | 22787759 | primary | 38.2722395 | 140.8690382 | 38.2688936 | 140.8700114 | 0 | LINESTRING (140.8690382 38.2722395, 140.869062... | 0.381072 |
# リンク情報の整理
OSMには交通流の上り下りがないので、上り下りの要素を追記する。
追記方法は下記に示す。
- 上り下り(順流・逆流)のIDを付与
- 始点終点ノードIDとリンクIDを付与
# 上り下り(順流・逆流)のIDを付与
# 新規DF(順流)
osm_links_0 = osm_links.copy()
# 方向のID
osm_links_0["vec"] = "0"
# 新規DF(逆流)
osm_links_1 = osm_links.copy()
# 始点の緯度経度
osm_links_1["init_lat"] = osm_links["term_lat"]
osm_links_1["init_lon"] = osm_links["term_lon"]
# 終点の緯度経度
osm_links_1["term_lat"] = osm_links["init_lat"]
osm_links_1["term_lon"] = osm_links["init_lon"]
# 方向のID
osm_links_1["vec"] = "1"
# 縦に結合
osm_links_vec = pd.concat([osm_links_0, osm_links_1], axis=0)
# 始点終点ノードIDとリンクIDを付与
緯度経度からMD5ハッシュ値に変換して、始点終点ノードを示すIDと、osm_idの代わりとなるリンクIDを付与します。
リンクIDは、違う道路が重なるように配置されていることもあるので、osm_id,geometry,vecごとに付与します。
# # 始点のノードID(MD5ハッシュ値)
dat_init = osm_links_vec["init_lon"].astype(str) + osm_links_vec["init_lat"].astype(str)
osm_links_vec["init_node"] = dat_init.apply(lambda x: hashlib.md5(x.encode()).hexdigest())
# 終点ノードのノードID(MD5ハッシュ値)
dat_term = osm_links_vec["term_lon"].astype(str) + osm_links_vec["term_lat"].astype(str)
osm_links_vec["term_node"] = dat_term.apply(lambda x: hashlib.md5(x.encode()).hexdigest())
# リンクID(MD5ハッシュ値)
dat_link = osm_links_vec["osm_id"] + osm_links_vec["geometry"] + osm_links_vec["vec"]
osm_links_vec["link_id"] = dat_link.apply(lambda x: hashlib.md5(x.encode()).hexdigest())
# 確認
print(len(osm_links_vec), len(set(osm_links_vec["link_id"].tolist()))) #2300 2300
osm_links_vec.head(2)
| osm_id | fclass | init_lat | init_lon | term_lat | term_lon | raw | geometry | length | vec | init_node | term_node | link_id |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 20264905 | motorway | 38.8040228 | 141.0941549 | 38.8378983 | 141.0989473 | 0 | LINESTRING (141.0941549 38.8040228, 141.094536... | 3.944538 | 0 | 8240949eec7aa3e5b5a8e1f51360a602 | 91299405649127a0303c3ef39b5974db | e3571ca226da1b9fb9a7ff2bfe77c562 |
| 22787759 | primary | 38.2722395 | 140.8690382 | 38.2688936 | 140.8700114 | 0 | LINESTRING (140.8690382 38.2722395, 140.869062... | 0.381072 | 0 | fa83c41314133d91cda8a3baa10ec944 | 06996e1cc8fbea04d0036bdbed2fca3b | 8c86d9213130e00cb864a3a05b67e9fc |
# CSVの出力
# CSVの出力
osm_links_vec.to_csv("./data/sample_network.csv", index=False)
# QGISで投影
# リンクの分解
福島県四ツ倉駅周辺において、平面交差を基準に分解したOpenStreetMapを投影した画像とGoogleMapの画像を下記に示す。
上画像のように、平面交差を基準にリンクを分解できている。
# ノードIDの整合性
一応、ハッシュで付与したIDの整合も確認します。OKだった。
前々回の記事で作成
前回の記事で作成
本記事で作成
# まとめ
本記事では、平面交差の考慮と計算コストの節約のために、OpenStreetMapから東北地方の道路リンクを交差点基準で分解しました。
# 参考サイト
【PythonでGIS】GeoPandasまとめ (opens new window)
Kepler.glをJupyterNotebook上で扱ってみた (opens new window)
緯度経度から距離を算出するPythonのライブラリ ― GeoPy (opens new window)
【Pythonプログラミング】データのハッシュ化と実行速度検証 (opens new window)
