【R】Web地図上にメソ数値予報モデルGPVの10m高度風速値をプロット

# はじめに

Rの`leaflet`パッケージを使ってオープンソースで使えるWeb地図上に，気象庁のメソ数値予報モデルGPV（以下，MSM-GPV）から抽出した初期時刻の10m高度風速値をプロットする方法をメモします．

- [ぬるぬる動く！Web地図クライアント「Leaflet」を使おう！](https://shimz.me/blog/leaflet-js/4142)
- [leafletではじめるRによる地図プロット](https://kazutan.github.io/JapanR2015/leaflet_d.html)

> leafletとは、JavaScriptのオープンソースライブラリである“leaflet.js”をRでも利用できるようにしたパッケージです。これはhtmlwidgetsパッケージにより実現されています。JavaScriptを使わなくてもRだけで利用可能ということで、非常に注目を集めているパッケージです。

# 方法

## 手順

1. MSM-GPVから必要な初期時刻の10m高度風速値を抽出し，①NetCDF形式のファイルで一時保存する．
1. ①をラスタオブジェクトとして入力し，②leafletのマップに表示・保存する．
1. ②をpng形式とhtml形式にそれぞれ出力する．

## データ

本記事では，気象業務支援センターが配信している2017年12月10日12時（UTC）のサンプルデータを利用させていただきました．

- [実際に配信されたデータによるサンプル提供 - 一般社団法人 気象業務支援センター](http://www.jmbsc.or.jp/jp/online/c-onlineGsample.html)

## wgrib2

MSM-GPVはGRIB2に準拠したバイナリデータなので，`wgrib2`というライブラリを利用します．インストール方法は次の記事を参考にしてください．

- [Windowsでwgrib2を使いたい](https://teratail.com/questions/120794) <- Windows向け
- [今更ですが、気象データを解析してみた - Qiita](https://qiita.com/ysomei/items/12d6622bed030f6ac793) <- Mac向け
- [wgrib2ツールのメモ - Qiita](https://qiita.com/ShingoOikawa/items/77c837374721a2b945ba) <- Linux向け

# ソースコードと解説

## ライブラリと入出力ファイル

まずは，使用するライブラリと入出力ファイルを定義します．

```{r}
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で一時ファイルの利用 -- Qiita](https://qiita.com/kazutan/items/53253c06aaffcbfe3f7a)

> 一時ファイルとは一時的なデータを保持しておくファイルで、通常さくっと作られて用が済めば削除されるものです。Rでも自分で準備して利用することが可能で、ちょっとデータなどを預けておいたりするのに重宝します。

## wgirb2

Rで`wgrib2`のようなLinuxのコマンドライブラリを使用するときには，`system`関数を実行します．`system`関数の引数に実行したいwgrib2のコマンドとオプションを`paste`関数（あるいは`str_c`関数）で記述します．本稿では，`wgrib2`のコマンドオプションとして，`-match`で初期時刻の10 m高度における風速のU, V成分を，`-netcdf`でNetCDF形式の出力をそれぞれ指定しています．

- [How to use wgrib2 command (hysk)](http://www.hysk.sakura.ne.jp/Linux_tips/how2use_wgrib2)

```{r}
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`で変数を指定します．

```{r}
# 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ではじめるRによる地図プロット](https://kazutan.github.io/JapanR2015/leaflet_d.html)

```{r}
# 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形式のファイルが出力されなくなるのでご注意ください．

```{r}
mapview::mapshot(
  x = msm_map,        # 保存するオブジェクト
  file = ofile_png,   # 出力ファイル（png）
  url = ofile_html,   # 出力ファイル（html）
  remove_url = FALSE  #
  )
```

![MSM_20171210120000_U10.png](https://qiita-image-store.s3.amazonaws.com/0/237116/724d2492-66e0-628b-700a-9ed98dbf8f5d.png)