はじめに
Rのleafletパッケージを使ってオープンソースで使えるWeb地図上に，気象庁のメソ数値予報モデルGPV（以下，MSM-GPV）から抽出した初期時刻の10m高度風速値をプロットする方法をメモします．
leafletとは、JavaScriptのオープンソースライブラリである“leaflet.js”をRでも利用できるようにしたパッケージです。これはhtmlwidgetsパッケージにより実現されています。JavaScriptを使わなくてもRだけで利用可能ということで、非常に注目を集めているパッケージです。
方法
手順
- MSM-GPVから必要な初期時刻の10m高度風速値を抽出し，①NetCDF形式のファイルで一時保存する．
- ①をラスタオブジェクトとして入力し，②leafletのマップに表示・保存する．
- ②をpng形式とhtml形式にそれぞれ出力する．
データ
本記事では，気象業務支援センターが配信している2017年12月10日12時（UTC）のサンプルデータを利用させていただきました．
wgrib2
MSM-GPVはGRIB2に準拠したバイナリデータなので，wgrib2というライブラリを利用します．インストール方法は次の記事を参考にしてください．
- Windowsでwgrib2を使いたい <- Windows向け
- 今更ですが、気象データを解析してみた - Qiita <- Mac向け
- wgrib2ツールのメモ - Qiita <- Linux向け
ソースコードと解説
ライブラリと入出力ファイル
まずは，使用するライブラリと入出力ファイルを定義します．
library(tidyverse)
library(ncdf4)
library(raster)
library(leaflet)
ifile_grib <- "Z__C_RJTD_20171210120000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin"
ifile_date <- str_sub(string = basename(ifile_grib), start = 11, end = 24)
ofile_png <- str_c("MSM_", ifile_date, "_U10.png")
ofile_html <- str_c("MSM_", ifile_date, "_U10.html")
msm_netcdf <- tempfile(fileext = ".nc") # 一時ファイル
msm_netcdfはNetCDF形式の中間ファイルです．このファイルは作業ディレクトリに保存しないため，tempfile関数を使って一時ファイルとして扱います．
一時ファイルとは一時的なデータを保持しておくファイルで、通常さくっと作られて用が済めば削除されるものです。Rでも自分で準備して利用することが可能で、ちょっとデータなどを預けておいたりするのに重宝します。
wgirb2
Rでwgrib2のようなLinuxのコマンドライブラリを使用するときには，system関数を実行します．system関数の引数に実行したいwgrib2のコマンドとオプションをpaste関数（あるいはstr_c関数）で記述します．本稿では，wgrib2のコマンドオプションとして，-matchで初期時刻の10 m高度における風速のU, V成分を，-netcdfでNetCDF形式の出力をそれぞれ指定しています．
system(
paste(
"wgrib2", msm_grib,
"-match '(:UGRD:10 m above ground:anl:|:VGRD:10 m above ground:anl:)' ",
"-netcdf", msm_netcdf
)
)
ラスタ
leafletは，GIS（Geophysical Information System; 地理情報システム）の１つであり，地理空間データのデータ構造モデルである「ベクタモデル」と「ラスタモデル」のデータを地図上に表示させることができます．このうち「ラスタモデル」とは，細かく区切られたマス目であるグリッド（grid）を用いて連続面を表現したものを意味し，単に「ラスタ」とも呼ばれます．
プログラムの中では，MSM-GPVから抽出したNetCDF形式の気象データを等緯度経度座標上のラスタとして扱います．RでNetCDF形式のファイルからラスタとしてデータを入力するためには，raster関数を使用し，引数のvarnameで変数を指定します．
# raster ------------------------------------------------------------------
msm_u_rs <- raster(x = msm_netcdf, varname = "UGRD_10maboveground")
msm_v_rs <- raster(x = msm_netcdf, varname = "VGRD_10maboveground")
msm_ws_rs <- sqrt(msm_u_rs ^ 2 + msm_v_rs ^ 2)
# names(msm_ws_rs) <- "Wind_Speed_10m_above_ground"
leaflet
ラスタオブジェクトをleafletのWeb地図上で表示・保存します．下記のURLでは，leafletのさまざまな使い方を紹介しているので参考にしてみてください．
# leaflet -----------------------------------------------------------------
# カラーパレット
msm_ws_pal <-
colorNumeric(
palette = c(
"#000180", "#0007A7", "#0019CA", "#0035E6", "#0058F8",
"#017FFF", "#07A7F8", "#19CAE6", "#35E6CA", "#58F8A7",
"#80FF80", "#A7F858", "#CAE635", "#E6CA19", "#F8A707",
"#FF8001", "#F85800", "#E63500", "#CA1900", "#A70700"
),
domain = range(values(msm_ws_rs)), # カラーパレットの範囲
na.color = "transparent" # NAのときの色
)
msm_map <-
leaflet() %>%
# addTiles() %>%
addProviderTiles(providers$Esri.WorldImagery) %>%
addRasterImage(
x = msm_ws_rs, # ラスタ
colors = msm_ws_pal, # カラーパレット
opacity = 0.7 # 透過度
) %>%
addLegend(
position = 'topright', # 凡例の位置
pal = msm_ws_pal, # カラーパレット
values = values(msm_ws_rs), # 値
opacity = 0.7, # 凡例の透過度
title = "WS [m/s]" # 凡例のタイトル
)
マップの保存
最後に，leafletのオプジェクトをmapview::mapshot関数で出力します．fileでpng形式のファイルを，urlでhtml形式のファイルをそれぞれ指定します．remove_url = TRUEとすると，html形式のファイルが出力されなくなるのでご注意ください．
mapview::mapshot(
x = msm_map, # 保存するオブジェクト
file = ofile_png, # 出力ファイル（png）
url = ofile_html, # 出力ファイル（html）
remove_url = FALSE #
)
