OpenStreetMap日本の建物マッピング率 #GIS

「OSMって、地図の作成率として、日本をどのくらいカバーしてるんですかね」

OSMに関連する立場で企業さんとお話するときによくされるこの質問、これは実に困る質問です。
例えば、比較対象になる正解データ（100％の値）が正確にわからないとか、そもそも地物の書き方がOSMと他データで異なるとか、いろいろと課題があります。なので正確には誰もわからない、というのが第一回答になることが多いです。

ただ、全くわからない、というわけではなく、これまで東さんによる地物量の調査や、各種統計情報と比較したOSM wikiの比較調査などを代表として、数々のアプローチが取られてきています。
また、日本では主に瀬戸さんや西村さんを中心として、学術的な立場から、道路のカバレッジを調査したり、ゼンリンデータと比較してその状況を把握し、論文にまとめようとされています。

今回はそうした調査の一環として、国土地理院が配布している基盤地図情報データと比較して、OSMで建物データがどの程度描かれているか、調査してみました。¹

アプローチ

下準備
1. OSM boundariesから、都道府県のデータをダウンロード
2. ダウンロードしたデータから「関東地方」や「東北地方」など、地方単位の区域データを作成

OpenStreetMap側のデータ調査
1. 下準備で作成した区域データを使って、ohsomeからデータを抽出
2. 抽出されたデータをSpreadsheetにまとめる

基盤地図情報のデータ調査
1. 基盤地図情報を日本全国分ダウンロード
2. 基盤地図情報データから、建物の外周線データを抽出
3. 下準備で作成した区域データを用いて、全国分のデータを「切り抜き（clip）」する
4. 切り抜かれた地物の数を数える

結果

あまり先延ばししてもしょうがないので、先に結果を貼ります。
Kiban、が、基盤地図情報として対象地域に存在する建物ポリゴンの数で、日付のデータ部分は、各年の1月1日時点でのOSM側建物ポリゴンデータの数になります。

	Kiban	2011/1	2012/1	2013/1	2014/1	2015/1	2016/1	2017/1	2018/1	2019/1	2020/1	2021/1	2021/1 割合
北海道	3103498	26	1915	66418	146321	182803	279398	730009	1312566	1741489	1931573	2126447	68.52%
東北	6546804	352	46101	78614	201222	318796	442273	493193	544052	708476	830410	1225956	18.73%
関東	17621338	5930	48998	182014	373296	489714	648725	1007392	1262171	1635746	2472758	2909584	16.51%
中部	14283174	1071	7761	31877	141665	375600	777921	1110706	1465381	2348489	2875531	3260070	22.82%
近畿	11040484	599	10174	59956	180768	379259	535143	648984	837945	1431860	1919579	2261386	20.48%
中国	4958737	152	16623	77753	107862	152943	199308	462166	515985	925612	1015024	1072146	21.62%
四国	2993340	18	343	972	4240	24504	117968	129358	137021	262032	280180	322911	10.79%
九州・沖縄	7890299	124	4222	27816	34432	59600	83731	407431	562098	718271	938974	1378629	17.47%
日本全国	68437674	8272	136137	525420	1189806	1983219	3084467	4989239	6637219	9771975	12264029	14557129	21.27%

計算方法

それでは、実際に行った作業の詳細を説明します。

下準備

「関東地方」や「東北地方」など、地方単位の区域データを作成

まずは、関東地方、などのポリゴンを作成します。
OSM boundariesから、対象の地域をダウンロードします。

なお、エクスポートする際のパラメータは以下のように設定します。

Admin level: 4（都道府県）を指定します。
Land/Water: include waterを指定して、水域を含むようにします。これで、島嶼部についても調査の対象に含まれるようになります。
Simplify: Very littleあたりで大丈夫かと思います。対象データの簡素化を行ってファイルサイズを減らす項目です。

API url項目のところにあるDownloadボタンからデータをダウンロードし、QGISに読み込ませます。
GeoJSONなので、ドラッグ＆ドロップで簡単に読み込ませることができるはずです。

読み込んだGeoJSONは都道府県単位になっているので、結合して地域としてまとめます。

QGISのメニューから ベクタ->空間演算ツール->融合（dissolve）で結合し、新しいレイヤとして出力します。

なお、結合したレイヤは、Simplifyの影響により、中に空洞ができています。この空洞は、メニューバーから プロセッシング->ベクタ一般->ベクタジオメトリ->孔（hole）の削除で削除する事が可能です。

結合され、空洞が削除されたgeojsonはこの後で利用しますが、内部の形式の違いか、ohsomeのAPIが認識できません。

一度geojson.ioに読み込ませてExportする、という、おまじないをすると動くようです。

OpenStreetMap側のデータ調査

下準備で作成した区域データを使って、ohsomeからデータを抽出

ハイデルベルグ大学が開発しているohsomeというサービスがあります。ウェブサイトやAPIを使って、OSMに存在しているデータを調査することができます。
APIドキュメントはこちらです。

「関東地方」など、調査対象となる区域のGeoJSONが完成したら、それを適当なフォルダに配置します。
同じフォルダに、以下の内容でスクリプトを作成します。pythonです。²

import requests
GEOJSON = open("./Region-Kyushu-okinawa.geojson","r")
URL = 'https://api.ohsome.org/v1/elements/count'
data = {"bpolys": GEOJSON, "format": "json", "time": "2011-01-01/2021-01-01/P1Y", "filter": "building=* and type:way"}
response = requests.post(URL, data=data)
print(response.json())

スクリプト内でopenしているファイルは、対象の地域の名前に修正します。
また、timeは、今回、2011年1月1日から2021年1月1日まで（2011-01-01/2021-01-01）、1年間ごとに（P1Y）、という指定をしています。filter部分の"building=* and type:way"は、buildingタグのついたウェイ・オブジェクトすべて、という指定になっています。
これらの部分を変更することで、任意の単位やタグを対象とすることができます。

抽出されたデータをSpreadsheetにまとめる

スクリプトが出来上がったら、おもむろにpythonからkickします。

$ python ./test.py 
{u'attribution': {u'url': u'https://ohsome.org/copyrights', u'text': u'\xa9 OpenStreetMap contributors'}, u'result': [{u'timestamp': u'2011-01-01T00:00:00Z', u'value': 124.0}, {u'timestamp': u'2012-01-01T00:00:00Z', u'value': 4222.0}, {u'timestamp': u'2013-01-01T00:00:00Z', u'value': 27816.0}, {u'timestamp': u'2014-01-01T00:00:00Z', u'value': 34432.0}, {u'timestamp': u'2015-01-01T00:00:00Z', u'value': 59600.0}, {u'timestamp': u'2016-01-01T00:00:00Z', u'value': 83731.0}, {u'timestamp': u'2017-01-01T00:00:00Z', u'value': 407431.0}, {u'timestamp': u'2018-01-01T00:00:00Z', u'value': 562098.0}, {u'timestamp': u'2019-01-01T00:00:00Z', u'value': 718271.0}, {u'timestamp': u'2020-01-01T00:00:00Z', u'value': 938974.0}, {u'timestamp': u'2021-01-01T00:00:00Z', u'value': 1378629.0}], u'apiVersion': u'1.6.1'}

今回は countのエンドポイントにクエリを投げているので、しばらくするとJSON形式で、上記のようにデータが戻ってきます。

JSON形式で戻ってきたデータをSpreadsheetにまとめます。³
これで、OSM側の調査は完了です。

基盤地図情報のデータ調査

基盤地図情報を日本全国分ダウンロード

基盤地図情報のサイトからデータをダウンロードします。これが、地味に難関です。

単純にひたすらファイルをダウンロードするだけなのですが、ここは、圧倒的なUI&UXの悪さとの格闘です。僕も思わず、知人の地理院のひとに愚痴ってしまいました。。。⁴

また、基盤地図情報は、解凍するとかなり大きな容量になります。
日本全国の解凍後の容量は300Gを超えますので、残りディスク容量には注意してください。

基盤地図情報データから、建物の外周線データを抽出

解答した基盤地図情報は、以下のように1ファイル1レイヤの構造になっています。

FG-GML-534064-AdmArea-20210101-0001.xml
FG-GML-534064-AdmBdry-20210101-0001.xml
FG-GML-534064-AdmPt-20210101-0001.xml
FG-GML-534064-BldA-20210101-0001.xml
FG-GML-534064-BldL-20210101-0001.xml
FG-GML-534064-Cntr-20210101-0001.xml
FG-GML-534064-CommBdry-20210101-0001.xml
FG-GML-534064-CommPt-20210101-0001.xml
FG-GML-534064-ElevPt-20210101-0001.xml
FG-GML-534064-GCP-20210101-0001.xml
FG-GML-534064-RailCL-20210101-0001.xml
FG-GML-534064-RdCompt-20210101-0001.xml
FG-GML-534064-RdEdg-20210101-0001.xml
FG-GML-534064-WA-20210101-0001.xml
FG-GML-534064-WL-20210101-0001.xml
FG-GML-534064-WStrA-20210101-0001.xml
FG-GML-534064-WStrL-20210101-0001.xml
fmdid21-0101.xml

何が何やらですが、基盤地図情報のファイル仕様書を確認すると、今回対象となる建物データはBldAのファイルであることがわかります。⁵
ダウンロードして解凍したファイルから対象地域のBldAファイルだけを集めたら、それらをすべておもむろにQGISにドラッグ＆ドロップで投げ入れます。
日本全国だとそれなりの数（4781ファイル）になるので、マシンパワーに余裕のあるときにやりましょう。

また、ファイルが複数に分かれていると処理に時間がかかります。
なので、こちらもQGISのメニューからベクタ管理ツール->ベクタレイヤのマージで1つのレイヤに統合し、その後、出力されたレイヤを右クリックしてエクスポートすることで、ファイルやDBなどに保存しておくと幸せになれます。

私はエクスポートする際の出力形式として、今回はGeoPackageを使いました。
全国分のファイル容量は、概ね16.7GB程度になっています。⁶

下準備で作成した区域データを用いて、全国分のデータを「切り抜き（clip）」する

全国分のデータから、対象地域のデータのみを切り出します。

全国など広域のデータをQGISに読み込ませ、さらに、調査対象区域のGeoJSONを別レイヤとして読み込ませます。
QGISのメニューからベクタ->空間演算ツール->切り抜く（clip）を実行します。
これで、対象区域の建物ポリゴンだけを切り出して別レイヤにすることができます。

↓
↓

くり抜かれた地物の数を数える

くり抜き作業の結果、別レイヤとして出力されたデータを、レイヤ一覧で右クリックします。

地物の数を表示という選択肢がありますので、クリックすると、そのレイヤ内に存在するポリゴンの数が表示されます。

このポリゴンの数が、対象地域に存在するポリゴン、つまり基盤地図情報側の建物ポリゴンの数になります。

あとは、基盤地図情報側のデータもSpreadsheetにまとめて、完成です。

補遺

なぜ、都道府県ではなく各地方単位？

各都道府県レベルでも算出することは可能なのですが、まずは日本全国を見渡した傾向が掴みたかった、というのが理由のひとつ。それから、47都道府県分スクリプト回したりコピペしたりデータ作成したりするのが面倒だった、というのがもう一つの理由です。

割合、の意味について

今回の調査ではあくまでも、ウェイで描かれた建物ポリゴンの数だけを対象にした比較を行っています。
完全な、オブジェクト単位での1:1対応でデータを比較しているわけではないことに留意してください。

例えば、基盤地図情報上では2つの建物として描かれているものが、OSM上では1つになっていたりする場合もあります。また、OSMではマルチポリゴンとして描かれている建物も、少数ながら存在します。
そのため、地域によっては、描かれている率が100％を超える地域も存在すると思われます。

あくまでも、傾向調査であることに留意ください。

基盤地図情報の精度について

基盤地図情報は、1/2500精度でデータが作成されている地域と、1/25000精度で作成されている地域が存在します。

1/2500精度の地域（都市部中心）はそれなりに期待できるデータ品質なのですが、他方、1/25000精度の地域はかなり信頼性の低い品質となっています。正直、衛星画像からOSMでトレースするデータのほうがよほど信頼性が高いデータを作成できます。
こちらも、あくまでも傾向調査であることに留意ください。

Future Works

それぞれの建物データを誰が描いたかをoshファイルから掌握することで、特定地域で誰が主に活動し、誰の活動が止まってしまったのか、などの傾向を掴みたいと考えています。

これにより、次に日本全体のOSM活動を支援するためにどのようなアクションをとればよいか、検討したいと考えています。

注: 本当はもっとちゃんとしたプログラムを組んで半自動で行うのがよいと思うのですが、まずは初回調査ということで、極めて泥臭い方法で行っています。 ↩
ohsomeのサンプルほぼそのままです。あと、普通はそれぞれのgeojsonをforで回します。 ↩
戻り値をCSVにするオプションがあるのですが、僕の環境ではうまく動きませんでした。 ↩
ヤバさは例えば、先にログインしておかないとすべて選択しなおしになる、とか、一度にダウンロードする数の上限が非常に低いとか、まだ再チェックできていないのですが極めつけにヤバいな、と思ったのは、全国、を選択して各地方ごとのファイルをダウンロードしたあと、それを結合させると各地域のはじっこのほう（伊豆の東の方とか、鳥取の北半分とか）が対象に含まれていないことがわかる、などでしょうか。 .cab形式で圧縮されていることもあり、mac派のひとは地味にイラッとすることもあるでしょう。同僚のSmellmanさんも地味に引っかかったようでなによりです。 ↩
BldLという、建築物の外周線データもあります。ぶっちゃけ、ほぼ一緒のものです。ちゃんと見てないですが、ラインかポリゴンか、の違いかしら。 ↩
FlatGeobufという形式が素晴らしいそうなのですが、私の環境（QGIS 3.18 mac）ではファイル保存時にエラーが出て保存ができませんでした。GeoPackageでもそれなりに速度はでますので、とりあえず今回の用には足りています。 ↩