ill-identified diary

所属組織の見解などとは一切関係なく小難しい話しかしません

[R] [OpenStreatMap] 東京の道路データをグラフに要約する

この記事は最終更新日から3年以上が経過しています

概要

ということで, 地理情報関係のことに関心が向いたところで, こんなツイートが流れてきた.

都市の道路をベクトルとして扱い, polar plotする*1ことで, 都市構造の傾向を視覚化している. 本人が言及してるように, 東京がどうなっているのか知りたいので試してみた.

当初は Python でやろうとしたが, 結局 R の osmdata, sf, ggplot2 パッケージなどを使ってやった. 例によって, 地理情報処理システム (GIS) は私の学生時代の専門とも, 過去現在の職歴とも関係なく, 完全に趣味で追っているだけなので注意*2.

Python を使った方法 (未遂)

もとのツイートのリンク先は, Pythonosmnx モジュール開発者による紹介記事である. osmnx は, 地図情報を公開するサービス OpenStreetMap (以下, OSM) のデータを取得し, 地図情報を描画するモジュールである. なので名前はたぶん OpenStreatMapNetworkX の略だと思う. 今回の場合は道路情報をネットワークデータの形式でダウンロードして利用しているが, ポリゴンデータとしてダウンロードすることで, 都市の輪郭を描いたりもできる.

インストールするには, 上記のgithubリンクに書いてあるとおり, rtree, geopandas モジュールをまずインストールする. さらに, 外部プログラム GDAL と libspatialindex に依存しているので, それらもインストールする*3. 最後に pip または condaosmnx をインストールする.

試しに, モスクワ市*4について輪郭と公道を描画してみる. network_type=’drive’ というのが, 「車が通行できる車道 (ただし側道除く)」だけを取り出すというオプションになっている. 道路のネットワークデータはそのままだとデータが膨大になりすぎるので, こうやってフィルタをかけたほうが安全だと思う.

import osmnx as ox
import matplotlib.pyplot as plt
fig, _ = ox.plot_shape(
        ox.gdf_from_place("Moscow, Russia"))
fig.savefig("shape_Moscow.png")
fig, _ = ox.plot_graph(
        ox.graph_from_place("Moscow, Russia", network_type="drive"))
fig.savefig("load_Moscow.png")

画像は以下の通り. 輪郭図のアスペクト比がおかしいが, ox.plot_shape には簡単に補正する機能がなく, 自分で直すのも面倒だったので放置した. クレムリンを中心とした放射線状の道路と, 環状道路, そしてところどころモスクワ川に沿って曲がった道路があるのが見える. 右上の大きな空白は Лосиный Остров 国立公園だろう.

f:id:ill-identified:20180720005817p:plainf:id:ill-identified:20180720005826p:plain

続いて, ペテルブルク市*5.

fig, _ = ox.plot_shape(
        ox.gdf_from_place("St-Peterburg, Russia"))
fig.savefig("shape_Leningrad.png")
fig, _ = ox.plot_graph(
        ox.graph_from_place("St-Peterburg, Russia", network_type="drive")) 
fig.savefig("loads_Leningrad.png")

f:id:ill-identified:20180720010032p:plainf:id:ill-identified:20180720010041p:plain

ペテルブルク市も中心部から放射線状に伸びているが, モスクワより区画整理されているようだ. ネヴァ川やペトロパヴロフスク要塞跡地の輪郭もわかる.

で, 次に東京都の道路網を試してみたのだが, メモリを大量占拠してマシンがクラッシュして終わった. 結局このモジュールはAPI (Overpass API) を叩いてデータを落としているだけなので, クエリがうまく設定できてないのかもしれないが, クエリの使用を細かく調べるには時間がなさすぎる. 別に Python じゃなくてもいいだろうということで, Rに切り替えた. むしろ, グラフィック関係なら R のほうが楽にできそうだ*6.

Rを使った方法

目標であるグラフを書くのに必要な処理は,

  1. OSM から必要なデータを取得

  2. R が処理できる形式に変換

  3. (2) で変換したデータをグラフを描画できるように道路ごとの角度を計算し, 集計

  4. (3) の結果をもとにグラフ描画

という4つのステージに分けられる. このうち, 最後のは ggplot2 で可能だと予想できる. また tidy data の形に変換できれば集計も簡単なので, たぶん (3) もそうむずかしくない. (2) については, R が扱える形式として, sf オブジェクトがある. よって, OSM からデータを取得し sf パッケージに変換する必要がある. 探してみたところ, osmdata パッケージというのがあるのでこれを使ってみる*7. OSM では, 全ての地物に keyvalue のペアという形で属性データが登録されており, これを使って取得する地物を条件抽出することができる.

keyvalue の一覧は, 以下で確認できる.

https://wiki.openstreetmap.org/wiki/Map_Features なお, 道路全般を表すキーは, highway である*8.

osmnxソースコードを確認して, 同じクエリを再現したつもりだ*9.

しかし, osmdata で指定できるクエリの構文の種類は現状では少なく、選択した地物の全ての属性情報がダウンロードされたりと無駄が多く*10*11, やはりメモリが不足した*12. そこで, いきなりRに読み込むのではなく, 一旦ローカルに xml ファイルを保存してから改めて R に読み込むことにした.

ダウンロードしたファイルには東京23区の外まで伸びた道路も含まれているため, 東京都の行政区画の形に切り抜いた. 行政区画も OSM から読み込めるはずだが, なぜか海岸線を読み込めなかったため, どうも不格好な形になってしまう. そこで, 行政区画は国土交通省・国土数値情報から最新のものをダウンロードした.

これらの加工を sf パッケージの機能を利用して行ったが, このパッケージをまともに触ったのも今回が初めてなので, ひょっとしたら間違ってるかもしれない*13.

以下が, OSMで取得できた東京23区内の道路図である. データは道路の種類という属性も持っていた*14ので、種類別に色分けしてみた.

f:id:ill-identified:20180720010202p:plain

ところで, 行政区画 (ポリゴンデータ) に 道路 (線データ) を重ねるのは工夫が必要である. 以下ににかかれている方法を参考にしてほしい.

そして, 最終目標の polar plot は, やはり ggplot2 で簡単に作成できる. polar plot は, 要はヒストグラムを円周上に並べたものなので, geom_histogram()coord_polar() で作成できるはずだ。方向を表すデータに関する扱いは、同じく Tokyo R #71 で発表された、方向統計学に詳しい Kotsubo さん発表の資料 How to predict the wind direction of tomorrow も参考になる*15.

多少デザインを調整して作成したグラフが以下になる. polar plot も, 道路の種類で色分けしている.(2018/7/28: 南北の逆転などを微修正)

f:id:ill-identified:20180728125430p:plain

南側の凹みが気になる. 角度を atan2 関数で求めたので, なにか丸め誤差的なものが悪さしているのかもしれないが, 当初掲載したグラフは誤って coord_polar(start=pi) としていたので、つなぎ目がおかしかったので修正した。東西南北に伸びる道路が多いという傾向はつかめる. また, 一番面積の多そうな深青色は residential つまり生活道路なので, 住宅地が区画整理されていると読み取れる.

さらに, ggplot2 の強みとして, 集計方法を簡単に変更できるというものがある. 今回は, 道路の角度だけでなく長さも計算したので, 道路の長さで重み付けしたヒストグラムも作成してみた. これは, geom_histogram() 関数の aes 属性に, weight として追加することでできる.

f:id:ill-identified:20180728125445p:plain

すると, 今度は東西方向のビンが潰れてしまった. つまり, 東京23区内には, 南北に伸びる道路が多いということか.

プログラムはここでダウンロードできる

github.com


*1:見た目が鶏頭図と似ているが, polar plot と鶏頭図は違う.

*2:この記事を描いたのは月曜日のはずなのに今日まで問題が解決しなかったことからも, 不慣れであることがわかる.

*3:たとえば Ubuntu 系ならば apt install libgdal-dev libspatialindex-dev でインストールできる.

*4:先日の Tokyo.R のFIFA ワールドカップ 2018 絡みの発表にもろに影響されている.

*5:影響されている.

*6:osmnx の機能としてはpolar plot のプログラムは提供されておらず, matplotlib を使って描いている.

*7:sp オブジェクトや xml にも対応している.

*8:OSM はどうやら用語にイギリス英語を使うようにしているので, highway は北米英語でいう高速道路ではなく, 道路全般になる. 参考: https://wiki.openstreetmap.org/wiki/JA:%E9%81%93%E8%B7%AF

*9:ただし、osmnx は json 形式で落とすようになっているのに対し, osmdata は xml 形式でダウンロードする仕様になっている

*10:また, クエリを再現したつもりが検索条件で排除されているはずの属性をもつ地物が含まれてしまったりしている. ただし osmnx のほうも正しくダウンロードできているか検証しているわけではない.

*11:クエリを直接入力することはできる.

*12:データをダウンロードしてRオブジェクトに返還するパッケージにありがちな現象だ.

*13:特に線データを座標に変換する処理

*14:他にも、車線数なんかもあり、こういう情報と組み合わせても面白いかもしれないが、明らかに入力ミスと思われるものも数カ所あった。公共事業で作られているわけではないので、こういった品質の悪さには眼をつぶる必要がある.

*15:coord_polar() は円グラフを作るくらいにしか役に立たない (≒全く役に立たない) と思っていたが, こういう使い方があった.