Pacific-centered_WorldMap.R

###############################################################
# 必要なパッケージの読み込み
###############################################################

pacman::p_load(
  sf,
  rnaturalearth,
  tidyverse,
  ncdf4,
  rvest,
  reshape2)


###############################################################
# NOAAの海面水温データサイトから最新ファイル名と日付を取得
###############################################################

URL <- "https://www.ncei.noaa.gov/data/sea-surface-temperature-optimum-interpolation/v2.1"
URL_base <- paste0(URL, "/access/avhrr/")

lastDir <- read_html(URL_base) |>
  html_nodes("a") |>
  html_attr("href") |>
  tail(1)  
  
lastFile <- read_html(paste0(URL_base, lastDir)) |>
  html_nodes("a") |>
  html_attr("href") |>
  tail(1)
  
YMD <- ymd(substr(strsplit(lastFile , "\\.")[[1]][2], 1, 8))


###############################################################
# 最新ファイルをダウンロードして緯度経度情報を生成
###############################################################

download.file(paste0(URL_base, lastDir, lastFile), 
              destfile = file.path(tempdir(), lastFile), 
              mode = "wb")

nc4data <- nc_open(file.path(tempdir(), lastFile))
sst <- as.vector(ncvar_get(nc4data, "sst")) 
SST <- melt(matrix(sst, ncol = 1440, byrow = T)) |> 
  filter(value != is.na(value)) |> 
  mutate(lon = Var2 / 4 -  0.125,
         lat = Var1 / 4 - 90.125) |> 
  select(lon, lat, sst = value) 


###############################################################
# 世界地図と海面水温データを太平洋中心のEqual Earth図法に変換
###############################################################

world.map0 <- ne_countries(scale = "large", returnclass = "sf")
eq.earth <- "+proj=eqearth +lon_0=-180 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"

world.map <- world.map0 |> 
  st_break_antimeridian(lon_0 = 180) |> 
  st_transform(crs = eq.earth)

SST_sf <- st_as_sf(SST, coords = c("lon","lat")) |> 
  st_set_crs(st_crs(world.map0)) |> 
  st_transform(crs = eq.earth)


###############################################################
# 水温分布と世界地図をプロット
###############################################################

ggplot() +
  geom_sf(data = SST_sf, aes(color = sst), size = 0.1) +
  geom_sf(data = world.map) +
  scale_color_gradient(high = "#e9f0fe", low = "#071f56") +
  theme_void() +
  labs(color = "SST(°C)",
       title = paste0("Gridded and weighted NOAA OISST v2.1 estimates, ",
                      gsub("/0", "/", gsub("-", "/", YMD))),
       caption = URL) +
  theme(plot.margin = unit(c(1, 1, 1, 1), "lines"),
        plot.title = element_text(size = rel(1.1),
                                  hjust = 0.5),
        plot.caption = element_text(size = rel(0.9)),
        legend.title = element_text(size = rel(0.8)),
        legend.position = "bottom",
        legend.key.width = unit(3,"line"),
        legend.key.height = unit(0.5,"line"))