plot_PacificSST.R

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


ymin = -60
ymax =  60

box = c(xmin = -180, ymin = ymin, xmax = 180, ymax = ymax)
sf_use_s2(FALSE)
crs = "+proj=longlat +lon_0=-180 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +no_defs"
world_map <- ne_countries(scale = "large", returnclass = "sf") |> 
  st_break_antimeridian(lon_0 = 180) |>  
  st_transform(crs = crs) |> 
  st_crop(box)

## ggplot(world_map) + geom_sf()

###############################################################

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)
download.file(paste0(URL_base, lastDir, lastFile), 
              destfile = file.path(tempdir(), lastFile), mode = "wb")

YMD <- ymd(substr(strsplit(lastFile , "\\.")[[1]][2], 1, 8))

###############################################################

nc4data <- nc_open(file.path(tempdir(), lastFile))
sst <- as.vector(ncvar_get(nc4data, "sst")) 
sst_sf <- 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) |> 
  filter(lat >= ymin) |> 
  filter(lat <= ymax) |> 
  mutate(lon = lon - 180) |> 
  rowwise() |> 
  mutate(point = list(st_point(c(lon, lat)))) |> 
  st_as_sf() |> 
  st_set_crs(st_crs(world_map))

Ave <- sst_sf |> 
  as_tibble() |> 
  filter(lon >= -90) |> 
  filter(lon <= 120) |> 
  summarize(ave = mean(sst)) |> 
  round(digits = 1) |> 
  pull(ave)

gAve <- round(mean(sst_sf$sst), digits = 1)


ggplot() +
  geom_sf(data = sst_sf, aes(color = sst), size = 2) +
  geom_sf(data = world_map, linewidth = 0.1, 
          color = "gray60", fill = "#dcffca") +
  scale_color_gradient2(low = muted("blue"), mid = "white", high = "red", 
                        midpoint = Ave, breaks = seq(0, 30, 5)) +
  coord_sf(xlim = c(-90, 120), ylim = c(ymin, ymax)) +
  theme_void() +
  labs(color = "SST(°C)",
       title = paste0("Gridded and weighted NOAA OISST v2.1 estimates, ",
                      gsub("/0", "/", gsub("-", "/", YMD))),
       subtitle = paste0("Mean Temperature: ", Ave, "°C", 
                         " (Global Mean: ", gAve, "°C)", ";   60°S ~ 60°N"),
       caption = URL) +
  theme(plot.margin = unit(c(1, 1, 1, 1), "lines"),
        plot.title = element_text(size = rel(1.1)),
        plot.subtitle = element_text(size = rel(0.9)),
        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"))