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

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

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

ただ、全くわからない、というわけではなく、これまで[東さんによる地物量の調査](https://qiita.com/higa4/items/a69343a238f571e4bf28)や、[各種統計情報と比較したOSM wikiの比較調査](https://wiki.openstreetmap.org/wiki/Japan/Current_coverage)などを代表として、数々のアプローチが取られてきています。
また、日本では主に瀬戸さんや西村さんを中心として、学術的な立場から、[道路のカバレッジを調査](https://www.gisa-japan.org/conferences/proceedings/2017/papers/C73.pdf)したり、[ゼンリンデータと比較してその状況を把握](https://speakerdeck.com/tosseto/20201107sotmjp?slide=14)し、論文にまとめようとされています。

今回はそうした調査の一環として、国土地理院が配布している基盤地図情報データと比較して、OSMで建物データがどの程度描かれているか、調査してみました。[^1]
[^1]: 注: 本当はもっとちゃんとしたプログラムを組んで半自動で行うのがよいと思うのですが、まずは初回調査ということで、極めて泥臭い方法で行っています。

![スクリーンショット 2021-12-07 17.17.18.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/9928/b162d2f5-0112-27d2-e9a7-225cfcd93f6c.png)

## アプローチ
 下準備
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](https://osm-boundaries.com/)から、対象の地域をダウンロードします。
![スクリーンショット 2021-12-07 12.57.57.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/9928/5f955f62-dd13-3ca6-e63e-1b7789134c3c.png)

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

* Admin level: `4`（[都道府県](https://qiita.com/nyampire/items/423344fa75707dc138af)）を指定します。
* Land/Water: `include water`を指定して、水域を含むようにします。これで、島嶼部についても調査の対象に含まれるようになります。
* Simplify: `Very little`あたりで大丈夫かと思います。対象データの簡素化を行ってファイルサイズを減らす項目です。
![スクリーンショット 2021-12-07 12.58.35.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/9928/82bbe4cc-049f-6966-eda6-56d4e3e0b0a3.png)

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

読み込んだGeoJSONは都道府県単位になっているので、結合して地域としてまとめます。
![スクリーンショット 2021-12-08 15.15.06.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/9928/c64a6fc7-31f8-134e-4418-db664979fe55.png)

QGISのメニューから `ベクタ`->`空間演算ツール`->`融合（dissolve）`で結合し、新しいレイヤとして出力します。
![スクリーンショット 2021-12-08 15.15.00.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/9928/8b4bda94-cb1b-9612-5b64-36a79826de06.png)

なお、結合したレイヤは、Simplifyの影響により、中に空洞ができています。この空洞は、メニューバーから `プロセッシング`->`ベクタ一般`->`ベクタジオメトリ`->`孔（hole）の削除`で削除する事が可能です。
![スクリーンショット 2021-12-08 15.15.22.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/9928/c733704c-7525-4b9c-d127-062bb5c58799.png)
![スクリーンショット 2021-12-08 15.15.42.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/9928/869df6b5-126f-e1bf-28fb-c07e9da3f879.png)

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

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

![スクリーンショット 2021-12-09 11.26.21.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/9928/9ff182e9-e0f8-c1b1-1f32-17f264a6100f.png)

### OpenStreetMap側のデータ調査
#### 下準備で作成した区域データを使って、ohsomeからデータを抽出

ハイデルベルグ大学が開発している[ohsomeというサービスがあります](https://ohsome.org/apps/dashboard/)。ウェブサイトやAPIを使って、OSMに存在しているデータを調査することができます。
[APIドキュメントはこちらです](https://docs.ohsome.org/ohsome-api/v1/index.html)。

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

[^2]: ohsomeのサンプルほぼそのままです。あと、普通はそれぞれのgeojsonをforで回します。

``` 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）、という[指定をしています](https://docs.ohsome.org/ohsome-api/v1/time.html)。`filter`部分の`"building=* and type:way"`は、buildingタグのついたウェイ・オブジェクトすべて、という[指定になっています](https://docs.ohsome.org/ohsome-api/v1/filter.html)。
これらの部分を変更することで、任意の単位やタグを対象とすることができます。

#### 抽出されたデータをSpreadsheetにまとめる
スクリプトが出来上がったら、おもむろにpythonからkickします。

``` json
$ 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にまとめます。[^3]
[^3]: 戻り値をCSVにするオプションがあるのですが、僕の環境ではうまく動きませんでした。

これで、OSM側の調査は完了です。

### 基盤地図情報のデータ調査
#### 基盤地図情報を日本全国分ダウンロード
[基盤地図情報のサイト](https://fgd.gsi.go.jp/download/menu.php)からデータをダウンロードします。これが、地味に難関です。

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

[^4]: ヤバさは例えば、先にログインしておかないとすべて選択しなおしになる、とか、一度にダウンロードする数の上限が非常に低いとか、まだ再チェックできていないのですが極めつけにヤバいな、と思ったのは、全国、を選択して各地方ごとのファイルをダウンロードしたあと、それを結合させると各地域のはじっこのほう（伊豆の東の方とか、鳥取の北半分とか）が対象に含まれていないことがわかる、などでしょうか。 .cab形式で圧縮されていることもあり、mac派のひとは地味にイラッとすることもあるでしょう。[同僚のSmellmanさんも地味に引っかかったようでなによりです](https://smellman.hatenablog.com/entry/2021/12/08/001716)。

また、基盤地図情報は、解凍するとかなり大きな容量になります。
日本全国の解凍後の容量は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
```

何が何やらですが、[基盤地図情報のファイル仕様書](https://fgd.gsi.go.jp/otherdata/spec/FGD_DLFileSpecV4.1.pdf)を確認すると、今回対象となる建物データは`BldA`のファイルであることがわかります。[^5]
[^5]: BldLという、建築物の外周線データもあります。ぶっちゃけ、ほぼ一緒のものです。ちゃんと見てないですが、ラインかポリゴンか、の違いかしら。

ダウンロードして解凍したファイルから対象地域のBldAファイルだけを集めたら、それらをすべておもむろにQGISにドラッグ＆ドロップで投げ入れます。
日本全国だとそれなりの数（4781ファイル）になるので、マシンパワーに余裕のあるときにやりましょう。

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

私はエクスポートする際の出力形式として、今回はGeoPackageを使いました。
全国分のファイル容量は、概ね16.7GB程度になっています。[^6]
[^6]: [FlatGeobuf](https://github.com/flatgeobuf/flatgeobuf)という形式が素晴らしいそうなのですが、私の環境（QGIS 3.18 mac）ではファイル保存時にエラーが出て保存ができませんでした。GeoPackageでもそれなりに速度はでますので、とりあえず今回の用には足りています。

#### 下準備で作成した区域データを用いて、全国分のデータを「切り抜き（clip）」する
全国分のデータから、対象地域のデータのみを切り出します。

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

![スクリーンショット 2021-12-08 14.29.31.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/9928/e7253e07-e253-ba06-2a8b-421005499f95.png)↓
![スクリーンショット 2021-12-08 14.29.50.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/9928/b4384741-4a65-d694-5d5b-7a1cd8c59c10.png)↓
![スクリーンショット 2021-12-08 14.36.27.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/9928/7bc76d73-4f3c-5a55-8924-9321718056d9.png)

#### くり抜かれた地物の数を数える
くり抜き作業の結果、別レイヤとして出力されたデータを、レイヤ一覧で右クリックします。
![スクリーンショット 2021-12-08 14.36.49.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/9928/d303deb5-ea5c-2d8a-c66d-9f91c79312d5.png)

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

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

![スクリーンショット 2021-12-08 14.36.58.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/9928/5bb42c4d-e62e-a00a-3c41-7427a28ed81a.png)

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

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

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

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

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

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

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

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

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