R leaflet でGFS のデータを可視化

モチベーション

NOAA/GFS が提供する全球の気象データを世界地図上にプロットしたい、かつそれをR+leafletで実現したい、というのが目的である。まずは値を予想しやすい表面の気温のデータをプロットしてみたいと思う。

データのダウンロード

https://www.ncdc.noaa.gov/data-access/model-data/model-datasets/global-forcast-system-gfsから欲しいデータをダウンロードできる。今回の記事では GFS Forecast → GFS(003) → HTTPS → 202103/ → 20210313/ → gfs_3_20210313_0000_000.grb2 をダウンロードして使う。もしくは wget で持ってくる。

$ wget https://www.ncei.noaa.gov/data/global-forecast-system/access/grid-003-1.0-degree/forecast/202103/20210313/gfs_3_20210313_0000_000.grb2

wgrib2 で csv へ変換

wgrib2 のインストール方法については、「GRIB2形式データを処理するプログラムwgrib2（入門編）」を参考にした。wsl2の環境にも問題なくインストールできた。

まずは'TMP'のデータを確認する。

$ wgrib2 gfs_3_20210313_0000_000.grb2 -match 'TMP'
13:423125:d=2021031300:TMP:0.4 mb:anl:
17:642294:d=2021031300:TMP:1 mb:anl:
23:911332:d=2021031300:TMP:2 mb:anl:
29:1181962:d=2021031300:TMP:3 mb:anl:
35:1454074:d=2021031300:TMP:5 mb:anl:
41:1731742:d=2021031300:TMP:7 mb:anl:
47:1986986:d=2021031300:TMP:10 mb:anl:
54:2307057:d=2021031300:TMP:15 mb:anl:
58:2529653:d=2021031300:TMP:20 mb:anl:
65:2865907:d=2021031300:TMP:30 mb:anl:
72:3273588:d=2021031300:TMP:40 mb:anl:
76:3528130:d=2021031300:TMP:50 mb:anl:
89:4072472:d=2021031300:TMP:70 mb:anl:
96:4504379:d=2021031300:TMP:100 mb:anl:
111:5165943:d=2021031300:TMP:150 mb:anl:
126:5808080:d=2021031300:TMP:200 mb:anl:
141:6454240:d=2021031300:TMP:250 mb:anl:
156:7126939:d=2021031300:TMP:300 mb:anl:
171:7847337:d=2021031300:TMP:350 mb:anl:
186:8536624:d=2021031300:TMP:400 mb:anl:
201:9236967:d=2021031300:TMP:450 mb:anl:
215:9841060:d=2021031300:TMP:500 mb:anl:
230:10597059:d=2021031300:TMP:550 mb:anl:
244:11260786:d=2021031300:TMP:600 mb:anl:
258:11924592:d=2021031300:TMP:650 mb:anl:
272:12591263:d=2021031300:TMP:700 mb:anl:
287:13337352:d=2021031300:TMP:750 mb:anl:
301:14007771:d=2021031300:TMP:800 mb:anl:
315:14690381:d=2021031300:TMP:850 mb:anl:
330:15506826:d=2021031300:TMP:900 mb:anl:
344:16229801:d=2021031300:TMP:925 mb:anl:
358:16937091:d=2021031300:TMP:950 mb:anl:
373:17631588:d=2021031300:TMP:975 mb:anl:
386:18193713:d=2021031300:TMP:1000 mb:anl:
404:19083839:d=2021031300:TMP:surface:anl:
415:19412163:d=2021031300:TMP:2 m above ground:anl:
419:19685145:d=2021031300:APTMP:2 m above ground:anl:
445:20877968:d=2021031300:TMP:tropopause:anl:
454:21570612:d=2021031300:TMP:max wind:anl:
463:22273797:d=2021031300:TMP:80 m above ground:anl:
468:22615239:d=2021031300:TMP:100 m above ground:anl:
471:22821663:d=2021031300:TMP:1829 m above mean sea level:anl:
474:23035811:d=2021031300:TMP:2743 m above mean sea level:anl:
477:23251870:d=2021031300:TMP:3658 m above mean sea level:anl:
484:23748674:d=2021031300:TMP:30-0 mb above ground:anl:
497:24611563:d=2021031300:TMP:0.995 sigma level:anl:
510:25265997:d=2021031300:TMP:PV=2e-06 (Km^2/kg/s) surface:anl:
516:25566000:d=2021031300:TMP:PV=-2e-06 (Km^2/kg/s) surface:anl:

id = 404 に TMP:surface のデータがあるのでこれを使う。-csv オプションで csv に変換する。

$ wgrib2 gfs_3_20210313_0000_000.grb2 -match 'TMP:surface' -csv tmp_surface.csv
404:19083839:d=2021031300:TMP:surface:anl:
$ head tmp_surface.csv
"2021-03-13 00:00:00","2021-03-13 00:00:00","TMP","surface",0,-90,232.138
"2021-03-13 00:00:00","2021-03-13 00:00:00","TMP","surface",1,-90,232.138
"2021-03-13 00:00:00","2021-03-13 00:00:00","TMP","surface",2,-90,232.138
"2021-03-13 00:00:00","2021-03-13 00:00:00","TMP","surface",3,-90,232.138
"2021-03-13 00:00:00","2021-03-13 00:00:00","TMP","surface",4,-90,232.138
"2021-03-13 00:00:00","2021-03-13 00:00:00","TMP","surface",5,-90,232.138
"2021-03-13 00:00:00","2021-03-13 00:00:00","TMP","surface",6,-90,232.138
"2021-03-13 00:00:00","2021-03-13 00:00:00","TMP","surface",7,-90,232.138
"2021-03-13 00:00:00","2021-03-13 00:00:00","TMP","surface",8,-90,232.138
"2021-03-13 00:00:00","2021-03-13 00:00:00","TMP","surface",9,-90,232.138

とりあえず抜き出しはできてるよう。

raster, rgdal パッケージのインストール

Rのraster, rgdal というパッケージを使って地図データへのプロットを作っていく。ので、入れてない人はインストールから。

> install.packages("raster")

rgdal を導入するためには gdal が入っていることが必要。Rのエラーメッセージを確認しながら、もろもろのコマンドを apt を使ってインストールする。

$ sudo apt-get install gdal-bin proj-bin libgdal-dev libproj-dev

これを入れた後に

> install.packages("rgdal")

でインストールできる。

地図プロットの作成！

最初に作成した tmp_suraface.csv を raster データとして再構成し、地図上にaddRasterImage でプロットする。分解能が1度のデータという事は既知とする。

library(dplyr)
library(leaflet)
library(raster)

gfs_csv <- read.csv("./tmp_surface.csv", header = FALSE)
gfs_coords <- cbind(gfs_csv$V5, gfs_csv$V6)

cell_size <- 1.0
lon_min <- min(gfs_csv$V5)
lon_max <- max(gfs_csv$V5)
lat_min <- min(gfs_csv$V6)
lat_max <- max(gfs_csv$V6)
ncols <- ((lon_max - lon_min)/cell_size)+1
nrows <- ((lat_max - lat_min)/cell_size)+1 

grid_gfs <- raster(nrows=nrows, ncols=ncols, xmn=lon_min, xmx=lon_max, ymn=lat_min, ymx=lat_max, res=cell_size, crs="+proj=longlat +datum=WGS84")
gfs_rs <- rasterize(gfs_coords, grid_gfs, gfs_csv$V7, fun=mean)

# カラーマップの設定
gfs_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(gfs_rs)),  # カラーパレットの範囲
    na.color = "transparent" # NAのときの色

gfs_map <-
  leaflet() %>%
  # addTiles() %>%
  setView(lng=135,lat=35,zoom=2) %>%
  addProviderTiles(providers$Esri.WorldGrayCanvas) %>%
  addRasterImage(
    x = gfs_rs,
    colors = gfs_pal,
    opacity = 0.7
  ) %>%
  addLegend(
    position = 'topright',
    pal = gfs_pal,
    values = values(gfs_rs),
    opacity = 0.7,
    title = "TMP [K]"

print(gfs_map)

どうやら気温のプロットは作れたよう。ただし二つほど問題点がまだ残っている。

• 上記コードを走らせたときRに7つほどwarning が出る。
• 日付変更線あたりにマップの切れ目ができてしまう。

この二つを何とか解決したい。

マップの切れ目の解消

経度を見ると、最小値が -179 度になっている。

> lon_min
[1] -179

-180 度= 180 度なのでデータの欠損はないのだが、切れ目が出てしまっているのは、-180 度までデータが無いことが問題だと推測する。ので 180 度から -180 度のデータを作り、データフレームを足し合わせてraster データを再構成することで、この問題を解決する。上記のコードをちょっと書き換える。

library(dplyr)
library(leaflet)
library(raster)

gfs_csv <- read.csv("./tmp_suraface.csv", header = FALSE)

# csv に-180 度のデータを加える。
lon_180 <- gfs_csv[(gfs_csv$V5 == 180),]
lon_180$V5 <- lon_180$V5 * -1
gfs_csv_new <- rbind(gfs_csv, lon_180)

gfs_coords <- cbind(gfs_csv_new$V5, gfs_csv_new$V6)
cell_size <- 1.0
lon_min <- min(gfs_csv_new$V5)
lon_max <- max(gfs_csv_new$V5)
lat_min <- min(gfs_csv_new$V6)
lat_max <- max(gfs_csv_new$V6)
ncols <- ((lon_max - lon_min)/cell_size)+1
nrows <- ((lat_max - lat_min)/cell_size)+1 

grid_gfs <- raster(nrows=nrows, ncols=ncols, xmn=lon_min, xmx=lon_max, ymn=lat_min, ymx=lat_max, res=cell_size, crs="+proj=longlat +datum=WGS84")
gfs_rs <- rasterize(gfs_coords, grid_gfs, gfs_csv_new$V7, fun=mean)

gfs_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(gfs_rs)),  # カラーパレットの範囲
    na.color = "transparent" # NAのときの色
  )

gfs_map <-
  leaflet() %>%
  # addTiles() %>%
  setView(lng=135,lat=35,zoom=2) %>%
  addProviderTiles(providers$Esri.WorldGrayCanvas) %>%
  addRasterImage(
    x = gfs_rs,
    colors = gfs_pal,
    opacity = 0.7
  ) %>%
  addLegend(
    position = 'topright',
    pal = gfs_pal,
    values = values(gfs_rs),
    opacity = 0.7,
    title = "TMP [K]"
  )

きれいに線が無くなった。

おわりに

warning の問題は解決していないが、とりあえずの可視化は行えた。今はcsv を経由する原始的な方法でプロットを作っているので、他にもっと賢い方法があれば是非教えて頂きたい。

また、これは何時間もにわたってデータのある予報値の一部だけなので、全予報値をムービーで見せる or スライドバーを作ってデータを選択する、ができるようにしたい。