Survey Program FVCOM Temperature Matching

Approach for Interpolation of Surface + Bottom Temperature Values within FVCOM Mesh

# Packages
{
library(raster)        # netcdf data as raster
library(sf)            # vector spatial mapping/operations
library(fvcom)         # fvcom mesh and variable extractions
library(ncdf4)         # netcdf support
library(tidyverse)     # data wrangling and plotting
library(gmRi)          # color schemes and cloud storage paths
library(patchwork)     # plot arrangement
library(rnaturalearth) # coastlines and state polygons
library(geometry)      # bathycentric coordinates
library(Matrix)        # matrix algebra
library(sysfonts)      # font support
}


# Paths + conflicts
conflicted::conflict_prefer("select", "dplyr")
conflicted::conflict_prefer("filter", "dplyr")
proj_path <- cs_path("mills", "Projects/Lobster ECOL")

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")


# Gom EPU to clip
res_shapes <- cs_path("res", "Shapefiles")
epu_path <- str_c(res_shapes, "EPU/")
gom_poly <- read_sf(str_c(epu_path, "individual_epus/GOM.geojson"))

# Use GMRI style
use_gmri_style_rmd()

FVCOM Temperature Interpolations for Survey Stations Points

To acquire modeled bottom temperature values from the survey programs we need to interpolate values from the FVCOM hindcast data for that point in time using information from the nearest mesh nodes.

A way to do this is using linear interpolation or some other method. The following code will step through performing these steps for the federal bottom trawl survey. We can then compare how the FVCOM hindcast performs relative to the CTD data from the survey.

Load NEFSC 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 State Trawl Data

# # Path(s) to state trawl survey data
asmfc_path <- cs_path("mills", "Projects/Lobster/Trawl_fromASMFC")
# 
# # load the different state survey datasets - ASMFC editions
# load(str_c(asmfc_path, "MassDMF_Lobster_StationsNLengths_1979_2023_230718.Rdata"))
# load(str_c(asmfc_path, "MENH Trawl.Rdata"))
load(str_c(asmfc_path, "RI_LOBSTER_TRAWL_012324.Rdata"))
ri_trawl <- RI_Stations %>% 
  distinct(TrawlIdentity, Season, TowNumb, StratumCode, Year, Month, Day, LAT, LON)

# Maine NH Trawl Survey
menh_path <- cs_path("res", "Maine_NH_Trawl/data_2024")
menh_trawl <- read_csv(str_c(menh_path, "MaineDMR_Trawl_Survey_Tow_Data_2024-05-17.csv")) 

# Mass Trawl
mass_path <- cs_path("res", "MA_Trawl/Pull_20240716/Manipulated")
mass_trawl <- read_csv(str_c(mass_path, "MADMF_SVSTA_SW_2024.csv"))

# Set up dates for later
menh_trawl <- mutate(menh_trawl, date = as.Date(Start_Date))
mass_trawl <- mass_trawl %>% 
  mutate(date = as.Date(
    str_c(Year, 
          str_pad(Month, side = "left", pad = "0", width = 2), 
          str_pad(Day, side = "left", pad = "0", width = 2),
          sep = "-")))
ri_trawl <- ri_trawl %>% 
  mutate(date = as.Date(
    str_c(Year, 
          str_pad(Month, side = "left", pad = "0", width = 2), 
          str_pad(Day, side = "left", pad = "0", width = 2),
          sep = "-")))


# Plot them all
ggplot() + 
  geom_sf(data = new_england) +
  geom_point(data = ri_trawl, aes(LON, LAT, color = "RI Survey"), 
             shape = 3, size = 0.5, alpha = 0.4) +
  geom_point(data = menh_trawl, aes(Start_Longitude, Start_Latitude, color = "ME + NH Survey"), 
             shape = 3, size = 0.5, alpha = 0.4) +
  geom_point(data = mass_trawl, aes(`Start lon`, `Start lat`, color = "MASS Survey"), 
             shape = 3, size = 0.5, alpha = 0.4) +
  coord_sf(xlim = c(-72, -67), ylim = c(41, 44.75)) +
  scale_color_gmri() +
  labs(x = "Longitude", y = "Latitude")

Load VTS Survey Data

We also have point locations for the VTS survey. Load these and get unique date/time information with metadata to rejoin later.

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

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

# Mass
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 FVCOM Inventory

# Surface and Bottom only, from Dr. Li
fvcom_path <- cs_path("res", "FVCOM/Lobster-ECOL")


# Here are the files we have, loop over them later
fvcom_surfbot_files <- setNames(
  list.files(fvcom_path, full.names = T, pattern = ".nc"),
  str_remove(list.files(fvcom_path, full.names = F, pattern = ".nc"), ".nc"))


# Test File: GOM3 1978
# Load some daily FVCOM that we downloaded and averaged
fvcom_yrx <- nc_open(fvcom_surfbot_files["gom3_1978"])

# Get the mesh itself as a simple feature collection
gom3_mesh <- get_mesh_geometry(fvcom_yrx, what = 'lonlat')

Check Point-Mesh Coverage Overlap

The following map shows the coverage overlap between trawl survey locations and FVCOM GOM3.

# Map everything
ggplot() +
  geom_sf(data = gom3_mesh, alpha = 0.2, linewidth = 0.1, color = "gray30") +
  geom_sf(
    data = st_as_sf(trawl_locations, coords = c("decdeg_beglon", "decdeg_beglat"), crs = 4326),
    aes(color = "NEFSC"),
    shape = 3, size = 0.2, alpha = 0.2) +
  geom_sf(data = vts_trips,
          aes(color = "VTS"),
          shape = 3, size = 0.2, alpha = 0.2) +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  theme(legend.position = "right") +
  scale_fill_gmri() +
  coord_sf(xlim = c(-78, -58), ylim = c(35.5, 46)) +
  theme_bw() + map_theme() +
  labs(
    title = "Coverage overlap of FVCOM and Sample Points",
    fill = "Area")

FVCOM Daily Product Matching

There are a number of projects where biological samples are taken, but environmental conditions might not also be taken simultaneously.

This adds the additional step of needing to index to the correct point in time and pull the correct node values associated with it.

The following code approaches identifying the correct date indices to loop/map through this interpolation process.

Get the Proper Time Index for Each File

# Get the dates within each file
fvcom_dates <- map_dfr(
  fvcom_surfbot_files,
   function(x){
     x_fvcom <- nc_open(x)
     timez <- ncvar_get(x_fvcom, "Times")
     nc_close(x_fvcom)
     return(
       data.frame("fvcom_date" = timez) %>% 
         mutate(time_idx = row_number()))
   }, .id = "fv_file")


# Join them by date to look at their matches
trawl_dates_matched <- trawl_locations %>% 
  mutate(tow_date = as.Date(est_towdate)) %>% 
  left_join(
    mutate(fvcom_dates, fvcom_date = as.Date(fvcom_date)),
    join_by(tow_date == fvcom_date))

Get Proper Element and Node Weights for Points

# Function to add the linear interpolation weights based on node coordinates
triangulation_linear_weighting <- function(pts_sf, fvcom_mesh){

    # Identify the triangles that overlap each point:
    # Use st_join to assign elem, p1, p2, p3 IDs to the points
    pts_assigned <- st_join(
      st_transform(pts_sf, st_crs(fvcom_mesh)),
      gom3_mesh, 
      join = st_within) %>% 
      drop_na(elem)
  

    # Iterate over the rows to add weights:
    pts_weighted <- pts_assigned %>% 
     base::split(., seq_len(nrow(.))) %>%
     purrr::map_dfr(function(pt_assigned){
    
      # Subset the relevant triangle from st_join info
      triangle_match <- fvcom_mesh[pt_assigned$elem,]
      
      # Build matrices for point to interpolate & of surrounding points:
    
      # Matrix for triangle
      # Use the triangles node coordinates from the sf geometries
      # Creates 3x3 matrix: row1 x coords, row 2, y coords, row three rep(1,3)
      node_vertices <- t(st_coordinates(triangle_match[1,])[1:3,1:3])
      
      # Make matrix from the points:
      # creates 3x1 matrix: x, y, 1
      point_coords <- matrix(
        c(st_coordinates(pt_assigned[1,]), 1), 
        nrow = 3)
      
      #### For Linear Interpolation:
      
      # Get inverse of the matrix
      inverse_coordmat <- solve(node_vertices)
      
      # Solve for the weights
      node_wts <- inverse_coordmat %*% point_coords %>%
        t() %>% 
        as.data.frame() %>% 
        setNames(c("p1_wt", "p2_wt", "p3_wt"))
      
      # Return with dataframe
      bind_cols(pt_assigned, node_wts)
    
    
    })
    # End Rowwise
    return(pts_weighted)
}

# Make the trawl locations an sf dataframe
trawl_pts_sf <- trawl_dates_matched %>% 
  st_as_sf(coords = c("decdeg_beglon", "decdeg_beglat"), crs = 4326, remove = F)


# Run for all points:
trawl_pts_weighted <- triangulation_linear_weighting(
  pts_sf = trawl_pts_sf, 
  fvcom_mesh = gom3_mesh) %>% 
  st_drop_geometry()

Apply Interpolations at Proper Time Step

Datewise Interpolation

This step can be looped on each date/year to minimize the amount of fvcom file opening/closing. Within each year we need to identify which timestep to extract data at, and then iterate on them.

#' Interpolate Values from FVCOM Mesh at Timestep
#' 
#' @description Takes a dataframe row containing node and time indices and interpolation weights and returns that contains time index from which to interpolate.
#' 
#' Should contain the following columns:
#' time_idx = integer timestep to use for interpolation
#' p1 = integer index value for node 1 surrounding the interpolation location
#' p2 = integer index value for node 2 surrounding the interpolation location
#' p3 = integer index value for node 3 surrounding the interpolation location
#' p1_wt = linear interpolation weight for p1
#' p2_wt = linear interpolation weight for p1
#' p3_wt = linear interpolation weight for p1
#'
#' @param dated_points_weighted dataframe row containing time index, node index, and node weight information columns
#' @param fvcom_nc FVCOM netcdf file to extract values from
#' @param fvcom_varid = String indicating variable in FVCOM netcdf to interpolate values with
#' @param var_out String to use as variable name in returned dataframe
#'
#' @return
#' @export
#'
#' @examples
interpolate_at_timestep <- function(dated_points_weighted, fvcom_nc, fvcom_varid, var_out){
        
    # Get the values of the variable of interest as vector
    node_vals <- ncvar_get(
      nc = fvcom_nc, 
      varid = fvcom_varid, 
      start = c(1, dated_points_weighted[["time_idx"]]),
      count = c(-1, 1))
      
  # Interpolate using the node numbers and weights
  dated_interpolation <- dated_points_weighted %>% 
    mutate(
      {{var_out}} := node_vals[p1] * p1_wt + node_vals[p2] * p2_wt + node_vals[p3] * p3_wt)
  
  return(dated_interpolation)
}

Iterate on Years to Interpolate All Points

Now that we can pull/interpolate for a specific timestep, we can loop over the yearly files and obtain surface and bottom temperatures. I’m using yearly timesteps to loop because we have yearly files, this way we only need to open each year once, then slice out values at the corresponding timesteps within them.

# Operate over years
trawl_fvcom_temps <- trawl_pts_weighted %>% 
  #filter(year(est_towdate) == 2000) %>% # For testing
  drop_na(time_idx) %>% 
  mutate(year = year(est_towdate)) %>% 
  split(.$year) %>% 
  map_dfr(
    .f = function(samples_year_x){
    
    # Get the file to open
    nc_name <- samples_year_x[["fv_file"]][[1]][[1]]
    
    
    # Open the corresponding Netcdf
    fvcom_yr_x <- nc_open(fvcom_surfbot_files[[nc_name]])
    
    
    # Split the samples by date
    locations_bydate <- samples_year_x %>% 
      base::split(., seq_len(nrow(.))) 
    
    
    # Iterate on those - do bottom temp and surface temp
    dates_interpolated <- locations_bydate %>%
      map(.f = ~interpolate_at_timestep(
        dated_points_weighted = .x,
        fvcom_nc = fvcom_yr_x,
        fvcom_varid = "surface_t", 
        var_out = "surf_temp_c")) %>% 
      map_dfr(.f = ~interpolate_at_timestep(
        dated_points_weighted = .x,
        fvcom_nc = fvcom_yr_x,
        fvcom_varid = "bottom_t", 
        var_out = "bot_temp_c"))
    
    return(dates_interpolated)
    
  })

Check Results

Compare against recorded bottom temperatures:

trawl_fvcom_temps %>% 
  ggplot(aes(bottemp_ctd, bot_temp_c)) + 
  geom_point(alpha = 0.2) +
  geom_abline(slope = 1, intercept = 0, color = "royalblue", linewidth = 1.5) + 
  labs(
    x = "CTD Bottom Temperature",
    y = "FVCOM Bottom Temperature",
    title = "CTD Bottom Temperature vs FVCOM Bottom Temperature")

Show a map of where biases a season looks like:

# Function to convert Fahrenheit to Celsius
# Used to label in F and C
fahrenheit_to_celsius <- function(f) {
  (f - 32) * 5/9
}


# Make Seasonal map
trawl_fvcom_temps %>% 
  mutate(fvcom_bias = bot_temp_c - bottemp_ctd) %>% 
  drop_na(fvcom_bias) %>% 
  st_as_sf(coords = c("decdeg_beglon", "decdeg_beglat"), crs = 4326) %>% 
  ggplot() +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  geom_sf(aes(color = fvcom_bias, fill = fvcom_bias), alpha = 0.5, size = 0.5) +
  geom_abline(slope = 1, intercept = 0, color = "royalblue", linewidth = 1.5) + 
  scale_fill_distiller(
    palette = "RdBu", 
    limits = c(-10, 10), 
    labels = function(x) {
      # Add both Fahrenheit and Celsius labels
      paste0(round(as_fahrenheit(x, data_type = "anomalies")), "°F / ", x, "°C")
    }) +
  scale_color_distiller(
    palette = "RdBu", 
    limits = c(-10, 10), 
    labels = function(x) {
      # Add both Fahrenheit and Celsius labels
      paste0(round(as_fahrenheit(x, data_type = "anomalies")), "°F / ", x, "°C")
    }) +
  facet_wrap(~fct_rev(season)) +
  theme_dark() +
  coord_sf(xlim = c(-78, -65), ylim =c(34, 45)) +
  labs(
    x = "Lon",
    y = "Lat",
    title = "FVCOM Bottom Temperature Biases vs. NEFSC Trawl CTD",
    subtitle = "FVCOM Bias = FVCOM Bottom Temperature - CTD Bottom Temperature")

Interpolate at VTS Locations

Repeat the process for the VTS survey locations. Streamline the code down to a single chunk for brevity, steps remain the same.

# Take the VTS Survey Location and Time information:
# 1. Join them by their date to match them to FVCOM file names and their time index
vts_dates_matched <- vts_trips %>% 
  mutate(Date = as.Date(Date)) %>% 
  left_join(
    mutate(fvcom_dates, fvcom_date = as.Date(fvcom_date)),
    join_by(Date == fvcom_date))

# 2. overlay the points to get the node ids
# Also applies the weights for linear interpolation
# Run for all points:
vts_pts_weighted <- triangulation_linear_weighting(
  pts_sf = vts_dates_matched, 
  fvcom_mesh = gom3_mesh) %>% 
  st_drop_geometry()


# 3. Map over the yearly files, perform interpolation
vts_fvcom_temps <- vts_pts_weighted %>% 
  drop_na(time_idx) %>% 
  mutate(year = year(Date)) %>% 
  split(.$year) %>% 
  map_dfr(
    .f = function(samples_year_x){
    
    # Get the file to open
    nc_name <- samples_year_x[["fv_file"]][[1]][[1]]
    
    
    # Open the corresponding Netcdf
    fvcom_yr_x <- nc_open(fvcom_surfbot_files[[nc_name]])
    
    
    # Split by row to iterate on each point
    locations_bydate <- samples_year_x %>% 
      base::split(., seq_len(nrow(.))) 
    
    # Iterate on those - do bottom temp and surface temp
    dates_interpolated <- locations_bydate %>%
      map(.f = ~interpolate_at_timestep(
        dated_points_weighted = .x,
        fvcom_nc = fvcom_yr_x,
        fvcom_varid = "surface_t", 
        var_out = "surf_temp_c")) %>% 
      map_dfr(.f = ~interpolate_at_timestep(
        dated_points_weighted = .x,
        fvcom_nc = fvcom_yr_x,
        fvcom_varid = "bottom_t", 
        var_out = "bot_temp_c"))
    
    return(dates_interpolated)
    
  })

Interpolate State Survey Locations

# 1. Join them by their date to match them to FVCOM file names and their time index

# Maine + NH
menh_dates_matched <- menh_trawl %>% 
  left_join(
    mutate(fvcom_dates, fvcom_date = as.Date(fvcom_date)),
    join_by(date == fvcom_date)) %>% 
  mutate(
    lon = Start_Longitude,
    lat = Start_Latitude) %>% 
  st_as_sf(coords = c("lon", "lat"), crs = 4326, remove = F)

# Mass
mass_dates_matched <- mass_trawl %>% 
  left_join(
    mutate(fvcom_dates, fvcom_date = as.Date(fvcom_date)),
    join_by(date == fvcom_date)) %>% 
  mutate(
    lon = `Start lon`,
    lat = `Start lat`) %>% 
  st_as_sf(coords = c("lon", "lat"), crs = 4326, remove = F)

# Rhode Island
ri_dates_matched <- ri_trawl %>% 
  left_join(
    mutate(fvcom_dates, fvcom_date = as.Date(fvcom_date)),
    join_by(date == fvcom_date)) %>% 
  mutate(
    lon = LON,
    lat = LAT) %>% 
  drop_na(lon, lat) %>% 
  st_as_sf(coords = c("lon", "lat"), crs = 4326, remove = F)

# 2. overlay the points to get the node ids
# Also applies the weights for linear interpolation

# Maine + NH
menh_pts_weighted <- triangulation_linear_weighting(
  pts_sf = menh_dates_matched, 
  fvcom_mesh = gom3_mesh) %>% 
  st_drop_geometry()

# Mass
mass_pts_weighted <- triangulation_linear_weighting(
  pts_sf = mass_dates_matched, 
  fvcom_mesh = gom3_mesh) %>% 
  st_drop_geometry()

# Rhode Island
ri_pts_weighted <- triangulation_linear_weighting(
  pts_sf = ri_dates_matched, 
  fvcom_mesh = gom3_mesh) %>% 
  st_drop_geometry()

# 3. Map over the yearly files, perform interpolation


# Maine + NH
menh_fvcom_temps <- menh_pts_weighted %>% 
  drop_na(time_idx) %>% 
  mutate(year = year(date)) %>% 
  split(.$year) %>% 
  map_dfr(
    .f = function(samples_year_x){
    
    # Get the file to open
    nc_name <- samples_year_x[["fv_file"]][[1]][[1]]
    
    
    # Open the corresponding Netcdf
    fvcom_yr_x <- nc_open(fvcom_surfbot_files[[nc_name]])
    
    
    # Split by row to iterate on each point
    locations_bydate <- samples_year_x %>% 
      base::split(., seq_len(nrow(.))) 
    
    # Iterate on those - do bottom temp and surface temp
    dates_interpolated <- locations_bydate %>%
      map(.f = ~interpolate_at_timestep(
        dated_points_weighted = .x,
        fvcom_nc = fvcom_yr_x,
        fvcom_varid = "surface_t", 
        var_out = "surf_temp_c")) %>% 
      map_dfr(.f = ~interpolate_at_timestep(
        dated_points_weighted = .x,
        fvcom_nc = fvcom_yr_x,
        fvcom_varid = "bottom_t", 
        var_out = "bot_temp_c"))
    
    return(dates_interpolated)
    
  })


# Mass 
mass_fvcom_temps <- mass_pts_weighted %>% 
  drop_na(time_idx) %>% 
  mutate(year = year(date)) %>% 
  split(.$year) %>% 
  map_dfr(
    .f = function(samples_year_x){
    
    # Get the file to open
    nc_name <- samples_year_x[["fv_file"]][[1]][[1]]
    
    
    # Open the corresponding Netcdf
    fvcom_yr_x <- nc_open(fvcom_surfbot_files[[nc_name]])
    
    
    # Split by row to iterate on each point
    locations_bydate <- samples_year_x %>% 
      base::split(., seq_len(nrow(.))) 
    
    # Iterate on those - do bottom temp and surface temp
    dates_interpolated <- locations_bydate %>%
      map(.f = ~interpolate_at_timestep(
        dated_points_weighted = .x,
        fvcom_nc = fvcom_yr_x,
        fvcom_varid = "surface_t", 
        var_out = "surf_temp_c")) %>% 
      map_dfr(.f = ~interpolate_at_timestep(
        dated_points_weighted = .x,
        fvcom_nc = fvcom_yr_x,
        fvcom_varid = "bottom_t", 
        var_out = "bot_temp_c"))
    
    return(dates_interpolated)
    
  })


# Rhode Island
ri_fvcom_temps <- ri_pts_weighted %>% 
  drop_na(time_idx) %>% 
  mutate(year = year(date)) %>% 
  split(.$year) %>% 
  map_dfr(
    .f = function(samples_year_x){
    
    # Get the file to open
    nc_name <- samples_year_x[["fv_file"]][[1]][[1]]
    
    
    # Open the corresponding Netcdf
    fvcom_yr_x <- nc_open(fvcom_surfbot_files[[nc_name]])
    
    
    # Split by row to iterate on each point
    locations_bydate <- samples_year_x %>% 
      base::split(., seq_len(nrow(.))) 
    
    # Iterate on those - do bottom temp and surface temp
    dates_interpolated <- locations_bydate %>%
      map(.f = ~interpolate_at_timestep(
        dated_points_weighted = .x,
        fvcom_nc = fvcom_yr_x,
        fvcom_varid = "surface_t", 
        var_out = "surf_temp_c")) %>% 
      map_dfr(.f = ~interpolate_at_timestep(
        dated_points_weighted = .x,
        fvcom_nc = fvcom_yr_x,
        fvcom_varid = "bottom_t", 
        var_out = "bot_temp_c"))
    
    return(dates_interpolated)
    
  })

Map All Interpolated Values

# Make a plot of everything
vts_sf <- vts_fvcom_temps %>% 
  mutate(Month = month(Date)) %>% 
  st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, remove = F)
nefsc_sf <- trawl_fvcom_temps %>% 
  mutate(Month = month(est_towdate)) %>% 
  st_as_sf(coords = c("decdeg_beglon", "decdeg_beglat"), crs = 4326)
menh_sf <- menh_fvcom_temps %>% 
  mutate(Month = month(date)) %>% 
  st_as_sf(coords = c("lon", "lat"), crs = 4326)
mass_sf <- mass_fvcom_temps %>% 
  st_as_sf(coords = c("lon", "lat"), crs = 4326)
ri_sf <- ri_fvcom_temps %>% 
  st_as_sf(coords = c("lon", "lat"), crs = 4326) 

# Build the map
ggplot() +
  geom_sf(
    data = vts_sf,
    aes(color = bot_temp_c, shape = "VTS Survey"), 
    size = 0.4, alpha = 0.9) +
  geom_sf(
    data =  nefsc_sf,
    aes(color = bot_temp_c, shape = "NEFSC Trawl Survey"),
     size = 0.4, alpha = 0.9) +
  geom_sf(
    data =  menh_sf,
    aes(color = bot_temp_c, shape = "MENH Trawl Survey"),
     size = 0.4, alpha = 0.9) +
  geom_sf(
    data =  mass_sf,
    aes(color = bot_temp_c, shape = "Mass Trawl Survey"),
     size = 0.4, alpha = 0.9) +
  geom_sf(
    data =  ri_sf,
    aes(color = bot_temp_c, shape = "RI Trawl Survey"),
     size = 0.4, alpha = 0.9) +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  facet_wrap(~Month, ncol = 3) +
  scale_color_distiller(palette = "RdBu") +
  # coord_sf(xlim = c(-71, -67), ylim = c(41, 45)) +
  coord_sf(xlim = c(-78, -65), ylim =c(34, 45)) +
  map_theme() +
  labs(
    x = "",
    y = "",
    title = "Trawl + VTS FVCOM Interpolations",
    color = "FVCOM\nBottom Temperature\n\u00b0C",
    shape = "Survey")

Export

# # # Save VTS Bottom temperatures
# write_csv(
#   x = trawl_fvcom_temps, 
#   file = here::here("data/point_location_temperatures", "NEFSCtrawl_location_fvcom_temps.csv"))
# write_csv(
#   x = vts_fvcom_temps, 
#   file = here::here("data/point_location_temperatures", "VTS_location_fvcom_temps.csv"))
# write_csv(
#   x = menh_fvcom_temps, 
#   file = here::here("data/point_location_temperatures", "MENH_trawl_fvcom_temps.csv"))
# write_csv(
#   x = mass_fvcom_temps, 
#   file = here::here("data/point_location_temperatures", "MA_trawl_fvcom_temps.csv"))
# write_csv(
#   x = ri_fvcom_temps, 
#   file = here::here("data/point_location_temperatures", "RI_trawl_fvcom_temps.csv"))