leafletとleaflet-velocityで風量情報を地図に表示してみる(windows)
はじめに
こんにちは．マエカワです．
突然ですが皆さん，気象情報を地図に表示してみたいですよね？分かります．みなまで言わなくとも，皆さんが気象情報を地図上に表示してwebアプリケーションチックにしてみたいことはお見通しです．
…
ということで今回は，気象情報(風量)を地図上に表示していこうと思います．具体的には，2020年3月1日の地表面から10メートル地点の風量を地図上に表示していきます．
実生活で必要になり試してみてはいるのですが，いろいろなサイトを経由して実装を進めなければならなかったので一回まとめておこうと思います．有識者の方からすると初歩中の初歩の内容になるので，あくまでも自分のような超初心者が「こういう記事があったらいいなぁ」と思っている内容を目指していこうと思います．この記事が「風量を地図に表示したいけど，どうやっていいのか分からない！！」というお悩みをお持ち方の助けになることができれば幸いです．
また，ハッカソン等の「なんでもいいからいい感じのもの作ってみて」系の課題に対して，いい見栄えのプロダクトを完成させることができるとも思っています．地図って災害情報をマップ上に表示するとか，意義のあるものにしやすいんですよね．
そんなこんなでやっていきましょう．どうかお付き合いのほどよろしくお願いいたします．
下準備
データの取得
京都大学が提供している気象庁データがあるので，そちらを使用させていただきましょう．まず，以下のサイトへ飛びましょう．
database.rish.kyoto-u.ac.jp
次に，「数値予報GPV」(下画像の赤枠で囲ってある場所です)をクリックすることでデータが格納されているページに飛ぶことができます．
このページから取得したい日付のディレクトリにアクセスしていきます．今回は2020年3月1日のデータが欲しいので，「2020/03/01」ディレクトリにアクセスしましょう．飛んだ先のページはこんな感じ．
www.data.jma.go.jp
今回欲しいのは地球全体に関する地表面からNメートル地点の風量を対象とするので，「Z_C_RJTD_{yyyymmddhhMMss}_GSM_GPV_{Rjp, Rgl}_FD{高度}_grib2.bin」と命名されたファイルをダウンロードしましょう．もし日本近辺のデータを対象とするならば，RglではなくRjpが名前に含まれるバイナリデータをダウンロードしましょう(glはglobal，jpはjapanの略)．2020年3月1日0時，地表面から10メートルの気象データが欲しいので，ダウンロードするファイルは「Z__C_RJTD_20200301000000_GSM_GPV_Rgl_FD0000_grib2.bin」になります(http://database.rish.kyoto-u.ac.jp/arch/jmadata/data/gpv/original/2020/03/01/Z__C_RJTD_20200301000000_GSM_GPV_Rgl_FD0000_grib2.bin)
ここまでで，データの取得は完了です．
wgrib2の導入
取得した気象データはバイナリで保存されているため，デコードして中身をGRIB2形式で取り出さなければなりません．そのためのツールとしてwgrib2が提供されているので，使わせてもらいましょう．ダウンロードはいたって簡単．マウスをポチポチするだけです．
www.ftp.cpc.ncep.noaa.gov
上のサイトでwindows用のパッケージが配布されていますので，ここからダウンロードしましょう．DLLファイルは使用しているPCによって必要なものが変わってくるようなので，まずは`wgrib2.exe`をダウンロード．実行すると「このDLLファイルが足りねぇ!!」という趣旨のアラートが出てきますので，足りないと言われたDLLファイルをその都度ダウンロードしましょう．
その後，wgrib2.exeと各DLLファイルが入っているフォルダにパスを通し，コマンドプロンプト上でwgrib2コマンドを使用できるようにしておきましょう．
grib2jsonの導入
今回使用するleaflet-velocityというライブラリは風量データをJSON形式で取得し，地図上に表示するものになります．なので，GRIB2形式で取得した風量データをJSON形式に変換しなければなりません．そのためにgrib2jsonを使用します．
github.com
上のリポジトリからクローンしましょう．
git clone "リポジトリのURL"
これでgrib2jsonの導入は完了です．
Mavenの導入
grib2jsonの導入は終わりましたが，これだけでは使うことができません．Mavenでパッケージングする必要があります．Mavenを導入していない場合は下記のサイト経由でインストールしましょう．
maven.apache.org
qiita.com
Mavenのパス通しが完了したら，grib2jsonをパッケージングしましょう．方法は簡単です．クローンしたgrib2jsonファイルに移動して，Mavenのpackageコマンドを打つだけ．
cd grib2json
mvn package
この操作によって，grib2jsonフォルダにパッケージングされたtargetフォルダが作成されます．次にtargetフォルダに入っているgrib2json-x.x.x-SNAPSHOT.tar.gzを解凍します．解凍したフォルダ内のbinフォルダにパスを通すことでgrib2jsonコマンドを使用することができるようになります．ここまでして，ようやくgrib2jsonが使えるようになりました．
leafletの導入
今回地図を表示するにあたってleafletというライブラリを使用していくので，ダウンロードしていきましょう．方法は簡単です．以下のサイトに飛んで，最新版leafletをダウンロードするだけ．
leafletjs.com
ZIP形式でダウンロードされるので，解凍してJSファイルとCSSファイルを取得しましょう．これで完了．
今回の記事ではleafletの使い方に関して書いていませんが，公式サイトのチュートリアルが優秀なのでそちらを参考にしてください．自分でもまとめようと思っているので，完成次第追記しますね．
leafletjs.com
leaflet-velocityの導入
以下のリポジトリからdistの中に入っているJSファイルとCSSファイルをダウンロードしましょう．
github.com
git clone "リポジトリのURL"
これで導入は終わり．簡単です．
ただし，今回はdistフォルダに入っている「.min」が含まれるJSファイルとCSSファイルを使用していきます．自分の環境だと，正規の(.minが入っていない)ファイルが欠損してしまっていたんですよね...．対応していない文字コードで開いたのか何なのかはわかりませんが，分かり次第追記したいと思います．
本実装
では，本実装を進めていきましょう．今回，ローカルで用意したデータをHTMLで表示することになるので，CORSエラーを回避するために動作確認用のサーバを立てて作業していきます．ブラウザ起動時にオプションを付けることでCORSエラーを回避することができるという記事がいっぱいあったのですが，なんだかうまくいかなかったのでApacheを使ってサーバを立てていこうと思います．Apacheの使い方(超初心者向け)は以下のサイトを参考に．
www.adminweb.jp
ファイル構造
今回想定するファイル構造は以下の通りです．今回はmapフォルダ下で作業していきます．
Apache
┗htdocs
┗map
┣ data
┃┗ Z__C_RJTD_20200301000000_GSM_GPV_Rgl_FD0000_grib2.bin
┣ js
┃┣ leaflet.js
┃┣ leaflet-velocity.js (今回は使わない)
┃┗ leaflet-velocity.min.js
┣ css
┃┣ leaflet.css
┃┗ leaflet-velocity.css (今回は使わない)
┃┗ leaflet-velocity.min.css
┗ sampleMap.html
データ処理
GRIB2ファイルの作成
それでは，気象庁のデータをleaflet-velocityで扱えるJSON形式まで変換していきましょう．ここまでの作業を正常に終えていれば，wgrib2コマンドとgrib2jsonコマンドを使用できるはずです．まずはwgrib2コマンドを使用し，バイナリファイルにどのようなデータが入っているのかを見てみましょう．
C:\map>cd data C:\map\data>wgrib2 Z__C_RJTD_20200301000000_GSM_GPV_Rgl_FD0000_grib2.bin 1.1:0:d=2020030100:PRMSL:mean sea level:anl: 1.2:0:d=2020030100:PRES:surface:anl: 1.3:0:d=2020030100:UGRD:10 m above ground:anl: 1.4:0:d=2020030100:VGRD:10 m above ground:anl: 1.5:0:d=2020030100:TMP:2 m above ground:anl: 1.6:0:d=2020030100:RH:2 m above ground:anl: 1.7:0:d=2020030100:LCDC:surface:anl: 1.8:0:d=2020030100:MCDC:surface:anl: 1.9:0:d=2020030100:HCDC:surface:anl: 1.10:0:d=2020030100:TCDC:surface:anl: (中略) 1.94:0:d=2020030100:HGT:20 mb:anl: 1.95:0:d=2020030100:UGRD:20 mb:anl: 1.96:0:d=2020030100:VGRD:20 mb:anl: 1.97:0:d=2020030100:TMP:20 mb:anl: 1.98:0:d=2020030100:VVEL:20 mb:anl: 1.99:0:d=2020030100:HGT:10 mb:anl: 1.100:0:d=2020030100:UGRD:10 mb:anl: 1.101:0:d=2020030100:VGRD:10 mb:anl: 1.102:0:d=2020030100:TMP:10 mb:anl: 1.103:0:d=2020030100:VVEL:10 mb:anl:
こんな感じにデータがずらっと並びます．これらのデータは「:(コロン)」で区切られており，左から
- ID
- 符号
- 日付(yyyymmddhh)
- 気象項目
- 条件
- anl(初期値) または fcst(予測値)
を意味しています．例えば，4行目の「1.2:0:d=2020030100:PRES:surface:anl:」は
- IDが1.2
- 符号は0
- 2020/03/01 00:00のデータ
- データタイプは「気圧」
- 地表面における
- 実値
が格納されているということが分かります．データタイプの意味は，以下のサイトに詳しく書かれているので参考にしてください．
qiita.com
www.hysk.sakura.ne.jp
このあたりの情報は一回自分でまとめようと思うので，記事にでき次第追記します．
今回は地表面から10メートル地点の風量データが欲しいので，「UGRD(東西風)」と「VGRD(南北風)」を格納しているデータを探していきましょう．条件で表示するデータを絞り込む場合は，wgrib2のmatchオプションを使用します．例えば，UGRD(東西風)のデータのみを抽出したい場合は以下のようにします．
C:\map\data>wgrib2 Z__C_RJTD_20200301000000_GSM_GPV_Rgl_FD0000_grib2.bin -match UGRD 1.3:0:d=2020030100:UGRD:10 m above ground:anl: 1.12:0:d=2020030100:UGRD:1000 mb:anl: 1.18:0:d=2020030100:UGRD:925 mb:anl: 1.24:0:d=2020030100:UGRD:850 mb:anl: 1.30:0:d=2020030100:UGRD:700 mb:anl: 1.36:0:d=2020030100:UGRD:600 mb:anl: 1.42:0:d=2020030100:UGRD:500 mb:anl: 1.48:0:d=2020030100:UGRD:400 mb:anl: 1.54:0:d=2020030100:UGRD:300 mb:anl: 1.60:0:d=2020030100:UGRD:250 mb:anl: 1.65:0:d=2020030100:UGRD:200 mb:anl: 1.70:0:d=2020030100:UGRD:150 mb:anl: 1.75:0:d=2020030100:UGRD:100 mb:anl: 1.80:0:d=2020030100:UGRD:70 mb:anl: 1.85:0:d=2020030100:UGRD:50 mb:anl: 1.90:0:d=2020030100:UGRD:30 mb:anl: 1.95:0:d=2020030100:UGRD:20 mb:anl: 1.100:0:d=2020030100:UGRD:10 mb:anl:
このデータの中で，条件が「10 m above ground」のものを探していきましょう．ちょうど2行目に書かれているIDが1.3のデータが該当します．次にVGRD(南北風)についても同様の作業をしていきましょう．
C:\map\data>wgrib2 Z__C_RJTD_20200301000000_GSM_GPV_Rgl_FD0000_grib2.bin -match VGRD 1.4:0:d=2020030100:VGRD:10 m above ground:anl: 1.13:0:d=2020030100:VGRD:1000 mb:anl: 1.19:0:d=2020030100:VGRD:925 mb:anl: 1.25:0:d=2020030100:VGRD:850 mb:anl: 1.31:0:d=2020030100:VGRD:700 mb:anl: 1.37:0:d=2020030100:VGRD:600 mb:anl: 1.43:0:d=2020030100:VGRD:500 mb:anl: 1.49:0:d=2020030100:VGRD:400 mb:anl: 1.55:0:d=2020030100:VGRD:300 mb:anl: 1.61:0:d=2020030100:VGRD:250 mb:anl: 1.66:0:d=2020030100:VGRD:200 mb:anl: 1.71:0:d=2020030100:VGRD:150 mb:anl: 1.76:0:d=2020030100:VGRD:100 mb:anl: 1.81:0:d=2020030100:VGRD:70 mb:anl: 1.86:0:d=2020030100:VGRD:50 mb:anl: 1.91:0:d=2020030100:VGRD:30 mb:anl: 1.96:0:d=2020030100:VGRD:20 mb:anl: 1.101:0:d=2020030100:VGRD:10 mb:anl:
条件に該当するのは2行目，IDが1.4のデータになります．
ということで，今回欲しいデータのIDが1.3と1.4であることが分かったので，これらのデータをGRIB2形式に変換していきましょう．変換にはwgrib2のdオプションでIDを指定した後，gribオプションをつけて保存ファイル名を記述するだけです．つまり，以下のようになります．
C:\map\data>wgrib2 Z__C_RJTD_20200301000000_GSM_GPV_Rgl_FD0000_grib2.bin -d 1.3 -grib u_wind.grib2 1.3:0:d=2020030100:UGRD:10 m above ground:anl: C:\map\data>wgrib2 Z__C_RJTD_20200301000000_GSM_GPV_Rgl_FD0000_grib2.bin -d 1.4 -grib v_wind.grib2 1.4:0:d=2020030100:VGRD:10 m above ground:anl:
これで，UGRDデータはu_wind.grib2，VGRDデータはv_wind.grib2というファイルにGRIB2形式で保存されました．
しかし，これだけではgrib2jsonでJSON形式に変換できません．先ほど作成した二つのファイルを連結する必要があります．なので，
C:\map\data>copy /b u_wind.grib2 + v_wind.grib2 wind_2020030100_10_global.grib2 u_wind.grib2 v_wind.grib2 1 個のファイルをコピーしました。
ここまで来てようやく，grib2jsonで扱うことができるUGRDとVGRDが統合されたwind_2020030100_10_global.grib2を作成することができました！！
JSONファイルの作成
ここまでで，GRIB2形式で保存されているwind_2020030100_10_global.grib2が作成できました．次は，このGRIB2データをleaflet-velocityが対応できるようにJSON形式に変換していきます．と言っても，今度は一行のコマンドで完了してしまいます．
C:\map\data>grib2json --names --data --o wind_2020030100_10_global.json wind_2020030100_10_global.grib2
完了．これで，「2019年3月1日0時における地表面から10メートルの風量データ」が格納されているJSON形式ファイルwind_2020030100_10_global.jsonを作成できました！！
ここまでの作業で，ディレクトリは以下のようになっているはずです．
Apache
┗htdocs
┗map
┣ data
┃┣ Z__C_RJTD_20200301000000_GSM_GPV_Rgl_FD0000_grib2.bin
┃┣ u_wind.grib2 (新しく追加！)
┃┣ v_wind.grib2 (新しく追加！)
┃┣ wind_2020030100_10_global.grib2 (新しく追加！)
┃┗ wind_2020030100_10_global.json (新しく追加！)
┣ js
┃┣ leaflet.js
┃┣ leaflet-velocity.js (今回は使わない)
┃┗ leaflet-velocity.min.js
┣ css
┃┣ leaflet.css
┃┗ leaflet-velocity.css (今回は使わない)
┃┗ leaflet-velocity.min.css
┗ sampleMap.html
コーディング
データの用意もできたので，いよいよ地図を表示させるコードを書いていきましょう．といっても，leaflet-velocity内のL.velocityLayerメソッドを使用するだけです．では，書いていきましょう．
<html> <head> <title>Sample Map | Show Wind Velocity</title> <!-- leaflet の導入 --> <script src="./js/leaflet.js"></script> <link rel="stylesheet" href="./css/leaflet.css"> <!-- leaflet-velocity の導入 --> <script src="./js/leaflet-velocity.min.js"></script> <link rel="stylesheet" href="./css/leaflet-velocity.min.css"> <!-- jQuery の導入 --> <script src="https://code.jquery.com/jquery-2.2.4.min.js"></script> <!-- JavaScript --> <script> let map; //地図の変数 function initMap(){ //地図の初期化 map = L.map("sampleMap", { minZoom: 2, maxZoom: 19 }); map.setView([36.5, 136], 3); //黒地マップのタイルを用意 let blackMap = L.tileLayer( 'https://cartodb-basemaps-{s}.global.ssl.fastly.net/dark_nolabels/{z}/{x}/{y}.png', { attribution: '© <a href="http://www.openstreetmap.org/copyright">OpenStreetMap</a> © <a href="http://cartodb.com/attributions">CartoDB</a>', subdomains: 'abcd', minZoom: 2, maxZoom: 19 } ); //JSONファイルの読み込み $.getJSON("./data/wind_2020030100_10_global.json", function(data){ //風量レイヤの準備 let velocityLayer = L.velocityLayer({ displayValues: true, //ここはtrueにしましょう data: data, //データの設定 minVelocity: 0, //表示する最小風速 maxVelocity: 10, //表示する最大風速 velocityScale: 0.03, //表示の長さ(大きくするほど表示時間が長くなります) colorScale: [ //表示の色設定 "#ffffff", //最低風速の色 "#ffff32", "#ffff64", "#ffff96", "#ffffc8", "#ffff00" //最高風速の色 ], opacity: 1 //透過率(1が最大) }); velocityLayer.addTo(map); }); blackMap.addTo(map); //地図タイルの設定 } </script> </head> <body onload="initMap()"> <div id="sampleMap" style="width: 100%; height: 100%;"></div> </body> </html>
コメントにも大体何をしているのかを書いていますが，補足説明をしていきましょう．
- 4～13行目：ライブラリの読み込み
- 21行目：地図変数の初期化．
- L.mapの第一引数は，HTMLにおけるid属性に対応するので，body内では設定したid属性を指定することで地図を表示することができます(61行目)．
- setViewメソッドで初期化(25行目)しないと，そのあとのaddToメソッドなどを使用できないので注意しましょう．
- 28行目：黒字の地図タイルの用意
- 今回風表示を際立たせて見せたいので，黒ベースの地図を使用することにしました．
- 38行目：JSONファイルを読み込み
- jQuery標準のgetJSONメソッドを使用しています．
- 41行目：風量レイヤの準備．いろいろオプションを付けることができます．詳しいオプションはgithubのREADME.mdに記述されているので，ご参考に．
- displayValuesはtrueにしましょう．表示されません．
- colorScaleでは，色コードを配列で指定することで表示する際の色を変更できます．低風速はインデックスの小さい色で表示されます．おそらく，風速を指定した色コードの数で等分割して，インデックスの小さい色コードほど低風速に割り振られているのだと思います．今回は，高速度ほど黄色で表示されるように調整してあります．
- 57行目：風量レイヤの追加．
- いろいろ試しましたが，どうやらgetJSONメソッド内でvelocityLayerの処理をする必要があるみたいです．このネストから離れてしまうと適用できなくなるので，ここだけは注意．使用でそうなっているのか，コードの書き方が悪いのかがまだわからないので，調査している最中です．
- 60行目：地図タイルの設定
ここまでで用意するものは全て用意できました．やったぜ．残るはhttpd.exeを起動しサーバを立て，ブラウザで「localhost:8080/map/sampleMap.html」にアクセスするだけです．
完成イメージ
以下の画像は，この記事に沿って準備・実装を進めた場合に想定される「localhost:8080/map/sampleMap.html」の表示画面です．多分，こうなる…はず…．なってほしい．
おわりに
お読みいただきありがとうございました．
今回は風量をleaflet-velocityを使って地図上に表示していきました．wgrib2やgrib2jsonの導入方法からデータ処理，コーディングまで包括的に，また誰でもわかるように書けていればいいなあと思います．書き間違いなどがあれば，コメントなどでご指摘してくださるとうれしいです．環境構築が一番面倒くさいんですよね．linuxならnpmでパッケージ配布されていたり，一つのコマンドで環境構築完了みたいになるんですけどね...．windowsは面倒くさい...．今回の記事が困っているどなたかの助けになれれば，とてもうれしいです．扱うデータ次第で下の図のような地図を作ることもできるので楽しいです．
地図アプリケーション，いろいろ応用できるとかなり楽しいので，興味を持たれた方はぜひ作ってみてはどうでしょうか．
追記：異なる空間参照系への対応(2020/03/04)
leaflet-velocityは様々な空間参照系(CRS)にもフレキシブルに対応してくれます．CRSとは，簡単に言ってしまばえば「メルカトル図法」などに代表される図法の違いと理解してとりあえず問題ありません．面積比を正しくしたものだったり，距離の比を正しくしたものだったり，それぞれの特徴によって地図の計上は変化していきます．QGISの公式ページにいろいろ載っているので，参考にしてみてください．
docs.qgis.org
この記事のデモとして使用したのはEPSG3857という，leafletのL.mapのデフォルトとして設定されているCRSです．また，オープンソースで提供されている地図タイルのほとんどはこのEPSG3857形式です．しかし，業務上どうしても違うCRSで表示したいという場合があると思います．その場合もleaflet-velocityは対応してくれます．どうすればいいのかというと，sampleMap.htmlの21行目の変数mapの定義を以下のように書き換えるだけです．
map = L.map("sampleMap", {
crs: L.CRS.EPSG4326,
minZoom: 2,
maxZoom: 19
});
何をしたのかというと，L.mapのcrsオプションをleafletに用意してあるオブジェクトL.CRS.EPSG4326に差し替えました．EPSG4326とは等距投影で表示されるもので，簡単に言えば「距離の縮尺がどの地点でも同じ」地図になります．ちょうど下の画像みたいな地図です．
ということで，異なるCRSで地図を表示ししたいという方向けの追記でした．
追記：デモページを公開しました(2020/03/16)
以下のページで，実際にぬるぬる風情報が動いている様子を見ることができます．簡単なものになりますが，よければ参考にしてみてください．
ymaekawa.php.xdomain.jp
読んでくださった方に感謝を．
マエカワ．
