R/leafletで波浪モデルの波高を可視化

8243 ワード
R wgrib2 leaflet cdo GEFS R テキストリンク

GEFSデータのダウンロード、wgrib2でnetcdfに変換。

wget https://nomads.ncep.noaa.gov/pub/data/nccf/com/gens/prod/gefs.20210401/06/wave/gridded/gefs.wave.t06z.c00.global.0p25.f003.grib2
wgrib2 gefs.wave.t06z.c00.global.0p25.f003.grib2 -netcdf tmp.nc

leafletで可視化。

library(raster)
library(leaflet)

r <- raster(x = "tmp.nc", varname = "HTSGW_surface")
leaflet() %>%
  addTiles() %>%
  addRasterImage(r, opacity = 0.5) %>%
  setView(lng = 0, lat = 0, zoom = 2)

lat&lngの最小最大値の確認。

library(ncdf4)

nc <- nc_open(filename = "tmp.nc")
lat <- ncvar_get(nc = nc, varid = "latitude")
range(lat)
#[1] -90  90
lng <- ncvar_get(nc = nc, varid = "longitude")
range(lng)
#[1]   0.00 359.75
tmp <- ifelse(lng > 180, -360 + lng, lng)
range(tmp)
#[1] -180.00  179.75

leafletで可視化するときは、lngが0~359.75度だと180度以上が表示できない（らしい）。

ということで、cdoで0~359.75度を-180~179.75度に変換。

cdo sellonlatbox,-180,179.75,-90,90 tmp.nc gefs.wave.t06z.c00.global.0p25.f003.nc

再度可視化。

library(raster)
library(leaflet)

r <- raster(x = "gefs.wave.t06z.c00.global.0p25.f003.nc", varname = "HTSGW_surface")
leaflet() %>%
  addTiles() %>%
  addRasterImage(r, opacity = 0.5) %>%
  setView(lng = 0, lat = 0, zoom = 2)

、、、180度のデータが無いから隙間があるけど、、、補間しないといけないんだと思うけど、、、まぁいいや。。。