In gulfofmaine/gmRi: Utility Functions for GMRI Research Team

knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE, fig.align = "center", eval = FALSE)

library(ncdf4)
library(sf)
library(raster)
library(stars)
library(tidyverse)
library(patchwork)
library(rnaturalearth)
library(gmRi)

#Geographic boundaries
northeast <- ne_states("united states of america") %>% 
  st_as_sf() %>% 
  filter(region_sub %in% c("New England", "Middle Atlantic", "South Atlantic"))
canada <- ne_states("canada") %>% st_as_sf()

# Plot default
theme_set(theme_bw())

# Potential map theme
gmri_map_theme <- list(
  theme(
    panel.border = element_rect(color = "black", fill = NA),
    plot.background = element_rect(color = "transparent", fill = "transparent"),
    line = element_blank(),
    axis.title.x = element_blank(), # turn off titles
    axis.title.y = element_blank(),
    legend.position = "bottom", 
    legend.title.align = 0.5))

Option 1. Downloading OISST Data Using THREDDS Server

This approach takes an input bounding box for a desired area, and downloads the desired gridded dataset using THREDDS.

NOTE: For most use cases this will be more time consuming than accessing what has already been placed on box. NOT the preferred usage, but still viable.

# Establish Desired Local Destination Path
oisst_path <- "/Users/akemberling/Documents/oisst_local"


# Pull data from Thredds server using bounding box
daily_sst_stack <- env_data_extract(data.set = "OISST", 
                                    dates = c("1982-01-01", "1984-12-31"), 
                                    #box = c(1, 359, -89, 89), 
                                    box = c(-77, -60, 35, 46), 
                                    out.dir = oisst_path, 
                                    mask = NULL)

Loading OISST Files as Raster Stack

Whether accessing data downloaded via THREDDS or those stored on box, any single netcdf can be loaded into the environment using raster::stack("file_location.nc)

# Load data into R as Raster stack
daily_sst_stack <- raster::stack(str_c(oisst_path, "/", "OISST.grd"))

Calculating Seasonal Averages from THREDDS Download

NOTE: Seasonal average function built to work off THREDDS download structure. Function does not work with oisst_window_load output.

From these raster stacks it is possible to get average temperature values for seasonal windows of time. For this we can use the oisst_period_means() function.

For this function you need to provide a table indicating the labels for the time windows and the start and end dates

#Set up the Breaks you want
season_breaks <- data.frame(
  "breaks" = c("SPRING", "FALL"),
  "start_date" = c(as.Date("2018-03-01"), 
                   as.Date("2018-09-01")),
  "end_date" = c(as.Date("2020-05-31"), 
                 as.Date("2020-11-30"))
)


#Calculate Means
fall_spring_means <- oisst_period_means(stack_in = daily_sst_stack, 
                                        projection_crs = 4326, 
                                        time_res_df = season_breaks)

Option 2. Access OISST Data from Box

Preferred Method

As a commonly used resource for high resolution sea surface temperature data it makes sense for us to have it on-hand rather than re-download from THREDDS as needed. For this reason there are numerous functions for grabbing data for specific areas and dates from our shared box drive.

This function also works for accessing pre-processed SST anomalies from the 1982-2011 climatology.

# OISST is on box within the NSF OKN Demo Data Folder
box_paths <- research_access_paths(os.use = "unix")
okn_path <- box_paths$okn


# Next we build a time/space extent that we want data for
data_window <- data.frame(
  lon = c(-72, -65), 
  lat = c(35,44), 
  time = as.Date(c("2018-08-01", "2020-12-31")))

# Load the desired data into a raster stack
# anomalies = FALSE returns raw sst
daily_sst_stack <- oisst_window_load(okn_path, data_window, anomalies = FALSE)

4. Extracting Daily Values with Point Locations

Another common need is to match temperature values in space and time with specific point locations. This is achieved with the

# Create spatial points object from survdat station data
res_path <- box_paths$res
load(paste0(res_path, "NMFS_trawl/Survdat_Nye_Aug 2020.RData"))
station_data <- survdat %>% filter(EST_YEAR >= 2018)
rm(survdat)

# Convert station data to sf object for spatial extraction
trawl_sf <- st_as_sf(station_data, coords = c("DECDEG_BEGLON", "DECDEG_BEGLAT"), 
                     crs = 4326)

# # Reproject if necessary
# project_utm <- 26919 #NAD1983 / UTM zone 19N got Maine
# trawl_sf_proj <- trawl_sf %>% st_transform(crs = project_utm)
# project_utm <- st_crs(trawl_sf_proj)


# Format Date Column to match with Raster Stack
trawl_sf <- trawl_sf %>% 
  mutate(DATE = lubridate::ymd_hms(str_c(EST_YEAR, "-", EST_MONTH, "-", EST_DAY, " 12:00:00")))

5. Single Layer Point Extractions

To get values from a single raster layer, regardless of date, you simply use rasters extract function.

#Single Year Test = Fall 1982

#Test Raster 
test_ras <- march_2020


#Test points = all stations all years
test_points <- bind_cols( lon = station_data$DECDEG_BEGLON, lat = station_data$DECDEG_BEGLAT)

# Extract temperature as new column
test_points$sst <- raster::extract(test_ras, test_points)

Mapping Extraction with sf

test_points <- st_as_sf(test_points, coords = c("lon", "lat"), crs = 4326)


#Test plot
ggplot() +
  geom_sf(data = test_points, aes(color = sst)) +
  geom_sf(data = northeast) +
  geom_sf(data = canada) +
  gmri_map_theme +
  scale_color_distiller(palette = "RdBu", na.value = "NA") +
  guides(fill = guide_colorbar(title = "SST - Celsius")) +
  coord_sf(xlim = c(-72, -65), ylim = c(39.5, 45), expand = FALSE) +
  labs(x = "", y = "", caption = "Point Extraction - March 2020 Mean SST")

6. Extracting All Years/Seasons with Raster Brick

#Extracting all years from the brick
test_points <- dplyr::select(station_data, 
                             lon = DECDEG_BEGLON, 
                             lat = DECDEG_BEGLAT, 
                             year = EST_YEAR)



#Extract specific to the sample year (or whatever the raster resolution is)
#Season is generated at random here because imported data does not contain that information
full_extraction <- test_points %>% 
  mutate(season = sample(c("SPRING", "FALL"), replace = T, size = 1),
         raster_res = str_c(season, year, sep = ".")) %>% 
  split(.$raster_res) %>% 
  imap_dfr(function(df, resolution){

    # Build Brick layer ID
    stack_id <- resolution

    # Build Raster Layer Index Number
    layer_index <- which(names(fall_spring_means) == stack_id)

    # Return NA's if there's no layer that matches
    if(length(layer_index) != 0){
        df$sst <- raster::extract(fall_spring_means[[layer_index]], df[, c("lon", "lat")])
      return(df) } else {
        df$sst <- rep(NA, nrow(df))
        return(df)
      }

  })

Mapping it with sf

full_extraction <- st_as_sf(full_extraction, coords = c("lon", "lat"), crs = 4326)

full_extraction %>% 
  filter(year %in% 1982:1984) %>% 
  ggplot() +
    geom_sf(aes(color = sst)) +
    geom_sf(data = northeast) +
    geom_sf(data = canada) +
    scale_color_distiller(palette = "RdBu", na.value = "NA") +
    guides(fill = guide_colorbar(title = "SST - Celsius")) +
    coord_sf(xlim = c(-72, -65), ylim = c(39.5, 45), expand = FALSE) +
    labs(x = "", y = "", caption = "Seasonal Survey Stations Extracting Matching SST Means") +
    facet_grid(year ~ season)