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

16663 ワード
GIS R wgrib2 leaflet MSM R テキストリンク

はじめに

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

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

方法

手順

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

データ

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

wgrib2

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

ソースコードと解説

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

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

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  #
  )