Survey_Locations_GLORYSTemp_extractions

Approach for Extracting Surface and Bottom GLORYS Temperature Values for NEFSC and VTS Surveys

# Loading/Separating Lobster-ECOL Spatial Areas

#### Libraries
library(gmRi)
library(here)
library(sf)
library(tidyverse)
library(rnaturalearth)
library(scales)
library(heatwaveR)
library(sysfonts)


# Degree celsius
deg_c <- "\u00b0C"

# Paths + conflicts
conflicted::conflict_prefer("select", "dplyr")
conflicted::conflict_prefer("filter", "dplyr")
proj_path <- cs_path("mills", "Projects/Lobster ECOL")
glorys_path <- "~/Documents/Repositories/glorys_northeast/GlORYS_surfbot_temps/"
glorys_base_fname <- "CMEMS_Northeast_TempSal_SurfaceBottomTemps_"

source(here::here("R/FVCOM_Support.R"))
theme_set(theme_bw() + map_theme())

# Shapefiles
new_england <- ne_states("united states of america", returnclass = "sf") %>% 
  filter(postal %in% c("VT", "ME", "RI", "MA", "CT", "NH", "NY", "MD", "VA", "NJ", "DE", "NC", "PA", "WV"))
canada <- ne_states("canada", returnclass = "sf")

# Use GMRI style
use_gmri_style_rmd()

Load Trawl Data

# Load the trawl data
trawl_path <- cs_path("res", "nmfs_trawl/SURVDAT_current")
trawl_raw <- read_rds(str_c(trawl_path, "survdat_lw.rds"))
# tidy it a little
trawl_dat <- gmRi::gmri_survdat_prep(
  survdat = trawl_raw$survdat, 
  survdat_source = "most recent", 
  box_location = "cloudstorage")

# Get distinct time/date/tow/locations
trawl_locations <- trawl_dat %>% 
  distinct(cruise6, station, stratum, tow, svvessel, est_towdate, season, decdeg_beglat, decdeg_beglon, bottemp_ctd = bottemp)

Load VTS Data

# Path to resources
vts_path <- cs_path("mills", "Projects/Lobster/VTS_fromASMFC")

# Maine
load(str_c(vts_path, "VTS_data.Rdata"))

# Massachusetts
load(str_c(vts_path, "VTS_MA_Proc_240201 all data for standardized index.Rdata"))


# Need Trips (for date) and Trawls, and Sites I think
vts_trips <- inner_join(
  bind_rows(Trips, mutate(Trips_MA, Fisher = as.character(Fisher))),
  bind_rows(Trawls, Trawls_MA), 
  join_by(TripId)) %>% 
  distinct(TripId, TrawlId, SiteId, Date, Longitude, Latitude) %>% 
  st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, remove = F) 

Load GLORYs Surface and Bottom Temperatures

GLORYS data is processed to contain only surface and bottom temperatures in the github.com/adamkemberling/glorys_northeast repository. These surface+bottom files are saved as annual files containing daily values for surface and bottom layers, as well as the depth at which the bottom layer was taken. This data can now be loaded with raster and point locations extracted from it.

# Load the regional glorys bottom temperature data I prepared
glorys_years <- c(1993:2019)

# Load bottom temperature variable as one big stack
glorys_btemps <- map(glorys_years, function(.x){
  raster::stack(str_c(glorys_path, glorys_base_fname, .x, ".nc"), 
    varname = "bottom_temp")}) %>% 
  raster::stack()

# Load surface too
glorys_stemps <- map(glorys_years, function(.x){
  raster::stack(str_c(glorys_path, glorys_base_fname, .x, ".nc"), 
    varname = "surface_temp")}) %>% 
  raster::stack()


# Plot one day
raster::plot(glorys_btemps[["X1995.05.11"]], main = "GLORYS BT, 1995-05-11")

raster::plot(glorys_stemps[["X1995.05.11"]], main = "GLORYS SST, 1995-05-11")

# Get Dates for date matching
# GLORYS date layers
glorys_dates <- names(glorys_btemps) %>% 
  str_sub(2, -1) %>% 
  str_replace_all("[.]", "-")

Date Matching

# This matchup process is likely 100x faster as a join
glorys_index_key <- data.frame("glorys_date" = glorys_dates) %>% mutate(
  glorys_date = as.Date(glorys_date), 
  stack_index = row_number())


# Can just join this, then split and map over the stack_index column
trawl_date_matches <- left_join(
  x = mutate(trawl_locations, est_towdate = as.Date(est_towdate)), 
  y = glorys_index_key, 
  by = join_by("est_towdate" == "glorys_date"))
vts_date_matches <- left_join(
  x = vts_trips, 
  y = glorys_index_key, 
  by = join_by("Date" == "glorys_date"))

Map through Date Matches, Extract Values

# Map through those, extract bottom temperature values
trawl_location_temps <- trawl_date_matches %>% 
  drop_na(stack_index) %>% 
  split(.$stack_index) %>% 
  imap_dfr(function(locations, stack_index){
    stack_index = as.numeric(as.character(stack_index))
  
    # Make an sf point from that row
    locations_sf <- st_as_sf(
      locations, 
      coords = c("decdeg_beglon", "decdeg_beglat"), 
      remove = FALSE, crs = 4326)
    
    # Use extract to pull the value out of that raster stack at the appropriate date
    bot_temp_extract <- raster::extract(x = glorys_btemps[[stack_index]], y = locations_sf)
    surf_temp_extract <- raster::extract(x = glorys_stemps[[stack_index]], y = locations_sf)
    locations_sf <- mutate(
      locations_sf, 
      surf_temp_c = surf_temp_extract,
      bot_temp_c = bot_temp_extract)
    return(locations_sf)
  
})


# Plot bottom temperatures
trawl_location_temps %>% 
  ggplot() +
  geom_sf(data = new_england) +
  geom_sf(aes(color = bot_temp_c), size = 0.5, alpha = 0.4) +
  scale_color_distiller(palette = "RdBu") +
  facet_wrap(~fct_rev(season)) +
  map_theme()  +
  coord_sf(xlim =c(-76, -66), ylim = c(35.5, 45)) +
  labs(x = "Lon",
       y = "Lat",
       title = "GLORYS Bottom Temperature + NEFSC Trawl Locations")

# Plot bottom temperatures
trawl_location_temps %>% 
  ggplot() +
  geom_sf(data = new_england) +
  geom_sf(aes(color = surf_temp_c), size = 0.5, alpha = 0.4) +
  scale_color_distiller(palette = "RdBu") +
  facet_wrap(~fct_rev(season)) +
  map_theme()  +
  coord_sf(xlim =c(-76, -66), ylim = c(35.5, 45)) +
  labs(x = "Lon",
       y = "Lat",
       title = "GLORYS Surface Temperature + NEFSC Trawl Locations")

Repeat for VTS Survey Locations

# Map through those, extract bottom temperature values
vts_location_temps <-  vts_date_matches %>% 
  drop_na(stack_index) %>% 
  split(.$stack_index) %>% 
  imap_dfr(function(locations, stack_index){
    stack_index = as.numeric(as.character(stack_index))
    
  # Make an sf point from that row
    locations_sf <- st_as_sf(
      locations, 
      coords = c("Longitude", "Latitude"), 
      remove = FALSE, crs = 4326)
    
   # Use extract to pull the value out of that raster stack at the appropriate date
   bot_temp_extract <- raster::extract(x = glorys_btemps[[stack_index]], y = locations_sf)
    surf_temp_extract <- raster::extract(x = glorys_stemps[[stack_index]], y = locations_sf)
    locations_sf <- mutate(
      locations_sf, 
      surf_temp_c = surf_temp_extract,
      bot_temp_c = bot_temp_extract)
    return(locations_sf)
  
})


# Plot those on a map
vts_location_temps %>% 
  ggplot() +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  geom_sf(aes(color = bot_temp_c), size = 0.5, alpha = 0.4) +
  scale_color_distiller(palette = "RdBu") +
  map_theme()  +
  coord_sf(xlim =c(-71, -67), ylim = c(41, 45)) +
  labs(x = "Lon",
       y = "Lat",
       color = "Bottom Temperature",
       title = "GLORYS Bottom Temperature + VTS Locations")

# Plot those on a map
vts_location_temps %>% 
  ggplot() +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  geom_sf(aes(color = surf_temp_c), size = 0.5, alpha = 0.4) +
  scale_color_distiller(palette = "RdBu") +
  map_theme()  +
  coord_sf(xlim =c(-71, -67), ylim = c(41, 45)) +
  labs(x = "Lon",
       y = "Lat",
       color = "Surface Temperature",
       title = "GLORYS Surfaace Temperature + VTS Locations")

Export

# # Save NEFSC Trawl and VTS Temperatures
# write_csv(
#   st_drop_geometry(trawl_location_temps),
#   here::here("data/point_location_temperatures", "NEFSCtrawl_location_glorys_temps.csv"))
# write_csv(
#   st_drop_geometry(vts_location_temps),
#   here::here("data/point_location_temperatures", "VTS_location_glorys_temps.csv"))