R + leaflet で GFS のデータでコントア（等圧線図）を作る

# はじめに
GFS（全球の気象予報）のデータに海抜0ｍでの気圧データが存在する。それを天気図のようなコントア（等圧線図）にして、leaflet 上に重ねたいというのがこの記事のモチベーションである。GFSデータの詳細については[Global Forecast System (GFS) のデータ概要](https://qiita.com/ryamamoto0406/items/6c52f7eb2615e2c0d479)を参照。

# データの前処理
[R leaflet でGFS のデータを可視化](https://qiita.com/ryamamoto0406/items/34f69e60ed9e33597f49)と同じ方法でcsvに読み込み、前処理、ラスターデータの作成を行う。

まず wgrib2 でデータの抽出。

```bash:shell
$ wget https://nomads.ncep.noaa.gov/pub/data/nccf/com/gfs/prod/gfs.20210416/00/atmos/gfs.t00z.pgrb2.1p00.f000 -O gfs.t00z.pgrb2.1p00.f000.grib2
--中略--
$ wgrib2 gfs.t00z.pgrb2.1p00.f000.grib2 -match "MSLET" -csv mslet.csv
557:33722120:d=2021041600:MSLET:mean sea level:anl:
```

csv の処理、ラスターとして取り込む。

```R:R
#パッケージの読み込み
library(leaflet)
library(raster)

#ファイル読込
gfs_pres_path <- "./data/mslet.csv"
gfs_pres_csv <- read.csv(gfs_pres_path, header = FALSE)

#180度の処理
lon_pres_180 <- gfs_pres_csv[(gfs_pres_csv$V5 == 180),]
lon_pres_180$V5 <- lon_pres_180$V5 * -1
gfs_pres_csv <- rbind(gfs_pres_csv, lon_pres_180)

#cellの構築
gfs_coords <- cbind(gfs_pres_csv$V5, gfs_pres_csv$V6)

cell_size <- 1.0
lon_min <- min(gfs_pres_csv$V5)
lon_max <- max(gfs_pres_csv$V5)
lat_min <- min(gfs_pres_csv$V6)
lat_max <- max(gfs_pres_csv$V6)

lon_ncols <- ((lon_max - lon_min)/cell_size)+1
lat_nrows <- ((lat_max - lat_min)/cell_size)+1

#raster データの構成
grid_gfs <- raster(nrows=lat_nrows, ncols=lon_ncols, xmn=lon_min, xmx=lon_max, ymn=lat_min, ymx=lat_max, res=cell_size, crs="+proj=longlat +ellps=WGS84 +no_defs +datum=WGS84")
gfs_pres_rs <- rasterize(gfs_coords, grid_gfs, gfs_pres_csv$V7, fun=mean)

#カラーパレットの準備
gfs_pres_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_pres_rs)),  # カラーパレットの範囲
    na.color = "transparent"
  )
```

これでラスターデータの準備ができた。

# コントアデータの作成
コントアのプロットを作るだけなら、Rだと`contour`コマンドでできる。

```R:R
contour(gfs_pres_rs)
```
![contour.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1077734/7c5dc354-437e-bdb8-6498-238d58dcba84.png)

これを leaflet 上に表示したい。ので、`rasterToContour`コマンドを使用する。これでラスターデータから等圧線のデータを作ることができる。

```R:R
#等圧線データに変換。
gfs_pres_con <- rasterToContour(gfs_pres_rs)
```

これでコントアのデータの準備ができた。

# コントアの表示
以上、ラスターデータ、コントアを leaflet 上に重ねて表示する。

```R:R
#パッケージの読み込み
library(mapview)

#leaflet の作成
map <- leaflet() %>% addTiles() %>%
  setView(lng = 135, lat = 35, zoom = 4) %>%
  addRasterImage(x = gfs_pres_rs, colors = gfs_pres_pal, opacity = 0.7) %>% #ラスターデータ（カラーマップ）の表示
  addLegend(position = 'topleft', pal = gfs_pres_pal, values = values(gfs_pres_rs), opacity = 0.7, title = "Pressure [Pa]") %>% #凡例
  addPolylines(data = gfs_pres_con, fillOpacity=5, fillColor = "transparent", opacity=10, weight=1) #コントアデータの表示
print(map)

#map の保存
mapview::mapshot(map, file = "./gfs_pres_contour.png")
```
![gfs_pres_contour.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1077734/defc77e8-de68-f222-387e-abbd00113629.png)

ラスターデータに応じたコントアを表示することができた。

# `addPolylines`のバグ？
`addPolylines`にはいくつか引数がある。その中でレイヤーの追加、削除に使う layerId なる引数を指定すると、コントア表示がおかしくなる。

```R:R
#表示に問題が無い場合（layerId 設定無し）
map_np <- leaflet() %>%
  addPolylines(data = gfs_pres_con, fillOpacity=5, fillColor = "transparent", opacity=10, weight=1)
mapview::mapshot(map_np, file = "./contour_noproblem.png")
```
![contour_noproblem.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1077734/a607dbbd-cf2e-78ce-c078-1f96fd699458.png)

```R:R
#表示に問題が有る場合（layerId 設定有り）
map_p <- leaflet() %>%
  addPolylines(data = gfs_pres_con, fillOpacity=5, fillColor = "transparent", opacity=10, weight=1, layerId = "test")
mapview::mapshot(map_p, file = "./contour_problem.png")
```
![contour_problem.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1077734/2beb59b0-b788-0372-b332-2fd1166fd982.png)

一つの layer しか表示できてないようだが、原因が良く分からない。これも対処方が分かっていないので、もし分かる方がいればコメントして頂けると有難いです。

# おわりに
`addPolylines` でコントアを表示することはできたが、数値入りのコントア表示になっていない。`contour`コマンドで表示された図では数値が入っていたので、できれば数値も含めたコントア表示を行いたい。#もし知ってる方が居たらコメント下さい！