GeoPandasによる国土数値情報行政区域データからの全国市区町村界データの作成

# 概要

全国市区町村界データや都道府県界データのGeometryを使いたくなることがあるが、思いのほか入手が大変だったりするのでメモ。

例えばESRIジャパンの提供する[全国市区町村界データ](https://www.esrij.com/products/japan-shp/)は使いやすいものの[ArcGIS](https://www.esrij.com/products/arcgis/)以外では使えないという使用規約になっている。

そこで、今回は[国土数値情報（行政区域データ）](https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N03-v2_3.html)をソースにしてGeoPandasの機能によりgeometry情報を統合して市区町村界データを作る方法を調べた。

# 実行

### 【注意】出典の明記（国土数値情報）

- 国土数値情報について、[利用規約](https://nlftp.mlit.go.jp/ksj/other/agreement.html)を確認、遵守する必要がある
    - 商用利用含め割と自由に利用出来るが、出典は明記しないといけない
    - 編集・加工も出来るが、その旨を書かないといけない

## 実行環境（概要）

- Linux(docker使用、ベースイメージ:Ubuntu18.04）
- jupyter lab上で確認
- geopandas==0.8.1

## データ取得と読み込み



- [国土数値情報](https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N03-v2_3.html)から「全国, 世界測地系」のデータをダウンロードして適当なパス上で解凍・展開
  - 例えば"N03-190101_GML.zip"などのファイル名
  - H31年のデータで400MB近くあって結構重い（年々増えてる?）ので注意

```python
import geopandas as gpd

gdf = gpd.read_file("/path/to/N03-19_190101.shp") # 展開したパスのshpファイルを指定
display(gdf) # 中身をざっと確認

	N03_001	N03_002	N03_003	N03_004	N03_007	geometry
0	北海道	オホーツク総合振興局	None	北見市	01208	POLYGON ((144.08144 44.12506, 144.08143 44.125...
1	北海道	オホーツク総合振興局	None	北見市	01208	POLYGON ((143.78333 44.18453, 143.78281 44.183...
2	北海道	オホーツク総合振興局	None	網走市	01211	POLYGON ((144.29495 44.00826, 144.29501 44.008...
3	北海道	オホーツク総合振興局	None	網走市	01211	POLYGON ((144.29602 44.01118, 144.29572 44.010...
4	北海道	オホーツク総合振興局	None	網走市	01211	POLYGON ((144.29266 44.01945, 144.29295 44.018...
...	...	...	...	...	...	...
117575	沖縄県	None	八重山郡	与那国町	47382	POLYGON ((122.98986 24.47191, 122.98999 24.471...
117576	沖縄県	None	八重山郡	与那国町	47382	POLYGON ((122.96591 24.47171, 122.96587 24.471...
117577	沖縄県	None	八重山郡	与那国町	47382	POLYGON ((122.99042 24.47248, 122.99036 24.472...
117578	沖縄県	None	八重山郡	与那国町	47382	POLYGON ((122.99804 24.47512, 122.99820 24.475...
117579	沖縄県	None	None	所属未定地	None	POLYGON ((127.84164 26.43531, 127.84167 26.435...

117580 rows × 6 columns
```

## dissolveによる統合

- https://geopandas.org/aggregation_with_dissolve.html を参考にdissolve機能で市区町村毎に統合していく
    - 指定した1つのカラムの単位でGeometry情報をまとめることが出来る
- 市区町村名のカラムは`N03_004`だが、「北区」「中央区」など被りそうな名前が多いので`N03_007`の行政区域コードを用いることにする
    - https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N03-v2_2.html で仕様（下にスクロールしていって「属性情報」の部分）を確認すると、行政区域コードは「都道府県名 + 市区町村名」に対応らしいので、これでおそらくOKのはず（万一違ってたらすみません）

```python
gdf_dissolved = gdf.dropna(subset=['N03_007']).dissolve(by='N03_007')
display(gdf_dissolved)

	geometry	N03_001	N03_002	N03_003	N03_004
N03_007					
01101	POLYGON ((141.34233 43.06682, 141.34285 43.066...	北海道	石狩振興局	札幌市	中央区
01102	POLYGON ((141.40839 43.18395, 141.40833 43.183...	北海道	石狩振興局	札幌市	北区
01103	POLYGON ((141.44707 43.15616, 141.44694 43.155...	北海道	石狩振興局	札幌市	東区
01104	POLYGON ((141.46244 43.10010, 141.46265 43.100...	北海道	石狩振興局	札幌市	白石区
01105	POLYGON ((141.38463 43.04670, 141.38479 43.046...	北海道	石狩振興局	札幌市	豊平区
...	...	...	...	...	...
47361	MULTIPOLYGON (((126.82034 26.26934, 126.82043 ...	沖縄県	None	島尻郡	久米島町
47362	MULTIPOLYGON (((127.76445 26.12312, 127.76447 ...	沖縄県	None	島尻郡	八重瀬町
47375	MULTIPOLYGON (((124.72860 24.65694, 124.72852 ...	沖縄県	None	宮古郡	多良間村
47381	MULTIPOLYGON (((123.75465 24.06330, 123.75463 ...	沖縄県	None	八重山郡	竹富町
47382	MULTIPOLYGON (((123.01358 24.43628, 123.01352 ...	沖縄県	None	八重山郡	与那国町

1902 rows × 5 columns
```
それらしい数に集約されている

補足：

- 取り扱いが面倒なのでここでは行政区域コードが`null`のものは落としている
- 元データも結構重いため、上記の処理時間は結構かかるかもしれない（筆者は2012年モデルのLet's note使用で3分くらい）


ここでは市区町村界について行ったが、都道府県界であれば`gdf.dissolve(by="N03_001")`のような感じで作れる
→筆者の環境では12minかかりました。。。

## 保存

- https://geopandas.org/io.html#writing-spatial-data を参考
- この際、encodingを`utf-8`にすることで（おそらく）ストレスを減らすことが出来る
    - encoding指定無しで保存したshpファイルに関して、オプション無しで読み込んだら文字化けが起こりました。。。


shpファイル:

```python
gdf_dissolved.to_file("/path/to/<ファイル名>.shp", encoding='utf-8')
```

geojson:

```python
# driverの指定が必要
gdf_dissolved.to_file("/path/to/<ファイル名>.geojson", driver="GeoJSON", encoding='utf-8')
```

## Simplify (オプション)

上記で作った市区町村界データに関して、そのままだと点数が多い？のかESRIジャパンの提供する市区町村界データ（zip圧縮で6MBちょっと）に比べるとデータサイズが桁違いに大きい
そのため、必要に応じて間引く（[simplify](https://geopandas.org/geometric_manipulations.html#GeoSeries.simplify)）ことによってデータサイズを落とす

**境界付近が少しガタガタになったりする**のでそのあたりは注意


- 使用する`GeoSeries`の[simplifyメソッド](https://geopandas.org/geometric_manipulations.html#GeoSeries.simplify)について、引数`tolerance: 数値（必須）`と`preserve_topology: bool(オプション、defaultではTrue)`を指定する
    - 説明が少なくて引数の意味が分かりづらいが、たぶん[shapelyのsimplify](https://shapely.readthedocs.io/en/latest/manual.html#object.simplify)のような感じ


例えば、

```python
gdf_simple = gdf_dissolved.copy()
gdf_simple['geometry'] = gdf_simple['geometry'].simplify(tolerance=0.00085)
```

のようにして作った`gdf_simple`をshpファイルとして保存すると12MBくらいになる。
単にdissolvedしただけのものは200MB以上あったので、結構軽く出来る。
（なおtoleranceの数値の調整は勝手が分からず結構大変だった。。。）

上記の設定で作ったものについて地図上で可視化すると、

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/680173/b4ebf6b5-e6cb-f199-78af-7fd6b049f9db.png)

のような感じ。simplifyの値次第なところもあるが、ぱっと見はそれらしい感じ。

拡大すると、

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/680173/6b81f99d-9590-a4ea-c012-5ea432471127.png)

simplifyを実行した分、境界付近が少し怪しくなることが分かる
→その時々の目的に応じてtolerance, preserve_topologyを適切に調整する必要がありそう

# 参考

手法などの参考：

- https://note.com/kinari_iro/n/nfee9bc97b6d7
- http://xnissy.hatenablog.com/entry/20160205/1454666764
- https://hayatoiijima.jimdofree.com/2017/11/14/%E9%83%BD%E9%81%93%E5%BA%9C%E7%9C%8C%E5%A2%83%E3%81%AE%E3%81%BF%E3%81%AEshp%E3%83%95%E3%82%A1%E3%82%A4%E3%83%AB%E3%81%AE%E4%BD%9C%E3%82%8A%E6%96%B9/

無償で使えるGISデータ：

- https://www.esrij.com/gis-guide/other-dataformat/free-gis-data/