気象庁の予報区等GISデータをベクトルタイルにした（全国・地方予報区）
===

<iframe width="100%" height="500" src="https://hfu.github.io/area-forecast/#7.11/36.21/138.973" frameborder="0"></iframe>

 _   _ _   _  __     __        _             
| | | | \ | | \ \   / /__  ___| |_ ___  _ __ 
| | | |  \| |  \ \ / / _ \/ __| __/ _ \| '__|
| |_| | |\  |   \ V /  __/ (__| || (_) | |   
 \___/|_| \_|    \_/ \___|\___|\__\___/|_|   
                                             
 _____ _ _        _____           _ _    _ _   
|_   _(_) | ___  |_   _|__   ___ | | | _(_) |_ 
  | | | | |/ _ \   | |/ _ \ / _ \| | |/ / | __|
  | | | | |  __/   | | (_) | (_) | |   <| | |_ 
  |_| |_|_|\___|   |_|\___/ \___/|_|_|\_\_|\__|

       Make the technology the easy part.

[IT dashboard の行政界データをベクトルタイルにした](https://hackmd.io/k_CpPVFPRpiMi16rO4x9Cg)の類似作業として、[予報区等GISデータ](http://www.data.jma.go.jp/developer/gis.html)をベクトルタイルにしたいと思っています。

# 利用条件に関する確認
[気象庁ホームページの利用規約](http://www.jma.go.jp/jma/kishou/info/coment.html)上、例えば次のように記載すれば大丈夫であると考えています。

> 気象庁「全国・地方予報区」を加工して作成

加えて、予報区等GISデータは「平30情使　第401号」です。[測量法第30条の承認を得て作成された測量成果を利用する際、新たな承認申請は必要ありません](http://www.gsi.go.jp/LAW/2930-qa.html#02)から、予報区等GISデータをベクトルタイルにするにあたって測量法第３０条に係る申請を行う必要はないと考えています。

# 新たな課題は何か
[IT dashboard の行政界データをベクトルタイルにした](https://hackmd.io/k_CpPVFPRpiMi16rO4x9Cg)では直面しなかった課題としては、次のようなものがあると認識しています。

1. データ量が多い。Shapefile 形式のデータを ZIP に固めたファイルで数百メガバイトを超えます。OpenStreetMap の地球規模データを80時間で変換できる国連ベクトルタイルのメソドロジーを使って変換できないものではもちろんありませんが、大規模データなりのしっかりとした工程管理をしたほうが良さそうですね。gh-pages の 1GB の壁についても意識する必要がありますね。今回の場合は、おそらく、データの種類ごとに別レポジトリにすれば大丈夫な気がします。
2. フォーマットが「Shapefile 形式のデータを ZIP に固めたもの」である。この場合、Shapefile を GeoJSON Text Sequence にするところまでで１レポジトリを作ることが良さそうです。

他方で、**ちゃんとしたデータ定義書が付属しているのはとても嬉しい**です。さすがですね。

# とりあえず「全国・地方予報区」を加工してベクトルタイルを作成
## 作業場所
作業場所を https://github.com/hfu/area-forecast に作りました。次の方法でコピーをチェックアウトできます。

git clone [e-mail]:hfu/area-forecast

## ダウンロード
まずはデータをダウンロードということで、 `rake download` すると `curl` でデータをダウンロードします。帯域を 180kB くらいに絞った提供をなさっているらしく、ダウンロードには 50 分程度かかりました。469MB ありますね。

### メモ: 「全国・地方予報区」データの意外な手強さ
「全国・地方予報区」データは意外に手強いです。なぜならば、ファイル名にメタ文字が入っているからです。

area-forecast-src hfu$ unzip -t 20190125_AreaForecast_GIS.zip 
Archive:  20190125_AreaForecast_GIS.zip
    testing: SEn\.dbf   OK
    testing: SEn\.shp   OK
    testing: SEn\.shx   OK
No errors detected in compressed data of 20190125_AreaForecast_GIS.zip.

「府県予報区等」データも同様のようです。

$ unzip -t 20190125_AreaForecastLocalM_prefecture_GIS.zip 
Archive:  20190125_AreaForecastLocalM_prefecture_GIS.zip
    testing: {\.dbf         OK
    testing: {\.shp         OK
    testing: {\.shx         OK
No errors detected in compressed data of 20190125_AreaForecastLocalM_prefecture_GIS.zip.

ですが、こういう手強さに簡単に負けるわけにはいきません。Rakefile で次のようにします（バッドノウハウ）。

task :unzip do
  %w{shp dbf shx}.each {|t|
    sh "unzip -p 20190125_AreaForecast_GIS.zip '*.#{t}' > a.#{t}"
  }
end

そうすると、次のようにして a.shp, a.dbf, a.shx が得られます。

area-forecast-src hfu$ rake unzip
unzip -p 20190125_AreaForecast_GIS.zip '*.shp' > a.shp
unzip -p 20190125_AreaForecast_GIS.zip '*.dbf' > a.dbf
unzip -p 20190125_AreaForecast_GIS.zip '*.shx' > a.shx

a.shp は 728MB もある巨大なファイルです。 `.gitignore` に `a.*` と入れておきます。

## ベクトルタイル変換
### 第一回変換
次のステップは Tippecanoe に GeoJSON Text Sequence を流し込んで mbtiles に変換するところですが、気象庁のデータはよく設計されているので、Node.js でスキーマ調整をする必要はなさそうです。コマンドラインで一気に Tippecanoe にパイプすることを試みましょう。

まず、`ogr2ogr -lco RS=YES -f GeoJSONSeq /vsistdout/ a.shp` で問題なく地物ストリームが出てくるところを確認しました。

その次に、まずは気軽に Tippecanoe にかけてみます。ズームレベルの割り当ては、Vector Tile optimizer を使ってデータ・ドリブンに考えてみることにしましょう。

ogr2ogr -lco RS=YES -f GeoJSONSeq /vsistdout/ a.shp | tippecanoe -o a.mbtiles

上記のコマンドを実行すると、tippecanoe がかなり CPU を使って、長い時間をかけて頑張ってくれます。意義のあるベクトルタイル化を試みている、ということだと思います。（一般的に、ポリゴンのベクトルタイル化はインテンシブな問題のようですね。）

### Vector Tile optimizer による診断
[Vector Tile optimizer](https://github.com/ibesora/vt-optimizer)を使って、a.mbtiles に含まれるベクトルタイルのサイズを分析し、適切な maxzoom / minzoom の設定を探し出します。

$ node ~/vt-optimizer/index.js -m a.mbtiles
  ✔ Parsing VT file contents

Vector Tile Info
Zoom levels:  0,
Format:  pbf
Center:  141.954346,39.951858,14
Layers: 

• a

Vector Tile Summary
Zoom level  Tiles   Total level size (KB)  Average tile size (KB)  Max tile size (KB)
----------  ------  ---------------------  ----------------------  ------------------  -
0           1       3.9091796875           3.9091796875            3.9091796875        ✓
1           1       6.2314453125           6.2314453125            6.2314453125        ✓
2           1       11.4677734375          11.4677734375           11.4677734375       ✓
3           3       24.5791015625          8.193033854166666       11.888671875        ✓
4           3       48.5234375             16.174479166666668      26.7705078125       ✓
5           8       106.0888671875         13.2611083984375        47.8837890625       ✓
6           21      253.919921875          12.091424851190476      129.9345703125      ✓
7           43      471.814453125          10.972429142441861      116.501953125       ✓
8           98      814.7666015625         8.313944913903061       81.015625           ✓
9           256     1422.740234375         5.557579040527344       62.234375           ✓
10          743     2424.291015625         3.2628412054172276      55.677734375        ✓
11          2332    4000.3916015625        1.7154337914075901      34.7490234375       ✓
12          7958    6667.5361328125        0.8378406801724679      35.904296875        ✓
13          29023   13172.517578125        0.45386478234934363     12.5693359375       ✓
14          109905  32686.2548828125       0.2974046211074337      9.5341796875        ✓

ちなみに、a.mbtiles のサイズを見てみると、次のようなものでした。

$ ls -lh a.mbtiles 
-rw-r--r--  1 hfu  staff    73M ---- ------ a.mbtiles

さて、vt-optimizer の出力に戻ってみると、z=6..7 あたりで、Max tile size が 100KB を超えているのがあまり嬉しくはありません。実用上、z=6..７を回避できるのであれば、回避した方が良いでしょう。

ただ、実際に地図を見てみると、そう単純に削除してみても実用性を失いそうですから、**とりあえずズームレベルが浅い方は全部変換してみて、データをみながらフィルタを試みる**ことにしましょう。

ズームレベルが深い方については、平均をみると z=10 と z=11 のあたりでカットオフしても良さそうに思えることから、**このデータは、mazxoom=10 にする**ことにします。こうすることで、ファイル数は概ね 2000 以下になるだろうということがわかりますので、 `gh-pages` で管理できる程度のベクトルタイルセットになるであろうということが想像できます。

### 第二回変換
レイヤ名を area-forecast にし、maxzoom を 10 に定め、ファイルシステムに未圧縮ベクトルタイルを出す方向で、第二回変換を行います。

ogr2ogr -f GeoJSONSeq /vsistdout/ a.shp | tippecanoe --force --output-to-directory=zxy --no-tile-compression --name='area-forecast' --attribution='気象庁「全国・地方予報区」を加工して作成' --layer='area-forecast' --maximum-zoom=10 --minimum-zoom=0

### ウェブサイト作成
IT Dashboard のときの成果を適宜再利用しながら、index.html, app.js, style.hjson を用意して、`rake` を実行して、サイトを作成します。必要に応じて `budo` を使えば、 GitHub へのコミット回数を節約しながら調整を行うこともできます。

http://hfu.github.io/area-forecast

Maputnik でみる場合 https://maputnik.github.io/editor?style=https://hfu.github.io/area-forecast/style.json

# 結語
- 次は「一次細分区域等」をやってみます。
- 今回、ズームレベルを10で打ち切っている結果、未圧縮ベクトルタイルで 27MB と軽量なデータになりましたが、その引き換えに、ズームすると当然に省略がみられます。用途によっては、深いところもタイルを作れば良いですが、最大ズームレベルを1増やすと概ねデータ量が５倍になることに留意する必要があります。