FVCOM - Maine Coastal Current

Detailin Maine Coastal Current Information with FVCOM

Published

April 24, 2024

Code 
library(raster)
library(sf) 
library(fvcom) 
library(ncdf4) 
library(tidyverse)
library(gmRi)
library(patchwork)
library(rnaturalearth)

# namespace conflicts
conflicted::conflict_prefer("select", "dplyr")
conflicted::conflict_prefer("filter", "dplyr")

# Paths
proj_path <- cs_path("mills", "Projects/Lobster ECOL")
fvcom_out <- str_c(proj_path, "FVCOM_support/")


maine <- ne_states(country = "United States of America", returnclass = "sf") %>% filter(name == "Maine")
ne_us <- ne_states("United States of America", returnclass = "sf") %>% filter(name %in% c("Maine", "Massachusetts", "Vermont", "New Hampshire", "Connecticut", "New York", "Rhode Islande"))

# Support functions for FVCOM
source(here::here("R/FVCOM_Support.R"))

Processing Maine Coastal Current Information from FVCOM

The Gulf of Maine Coastal Current (GMCC) can be divided into to principal branches (Eastern Maine Coastal Current EMCC & Western Maine Coastal Current WMCC), separated by a junction near Penobscot bay. This junction is important as it is a site where a variable portion of current flow shifts from continuing Southwest, and veers offshore. This shift in current flow has important implications for nutrient and plankton distribution, and larval transport. This is why this section of the current flow is used as an indicator of general circulation patterns within the Gulf of Maine.

Building on Past Work

A previous team member (Matt Dzaugis) was responsible for accessing and storing the data, and we are grateful for his time/effort in doing so.

Matt’s Previous Work on Ecosystem Indicators Code behind the report

Using a principal components analysis of surface current speed for the eastward (u) direction of an area just offshore of Penobscot Bay, the Maine Coastal Current can be decomposed into the Maine Coastal Current Index that captures the connectivity between the EMCC and WMCC. In this analysis, the first principal component, which capture 52.5% of the variability in the dataset, provides an index of connectivity. The second principal component, which captures 15.3% of the variability, is related to vorticity.

Code 
# Loading FVCOM

# Get some file names
gom3_files <- monthly_fvcom_inventory(
  start_yr = 2010,
  end_y = 2010,
  fvcom_folder = cs_path("res", "FVCOM/monthly_means"),
  fvcom_vers = "GOM3"
) 

# Open one
gom3_early <- nc_open(gom3_files[[1]])

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

Maine Coastal Current Region

Code 
# Make a polygon for the area of interest:
mcc_turnoff_coords <- tribble(
  ~"lon", ~"lat",
  -69.34, 43.21, # Bottom Left
  -69.78, 43.62,   # Off Popham
  -67.45, 44.34, # Off Jonesport 
  -67.4, 43.8, # Bottom right
  -69.34, 43.21, # Bottom left again
)

# Make it a polygon
mcc_turnoff_poly <- st_polygon(list(cbind(mcc_turnoff_coords$lon, mcc_turnoff_coords$lat))) %>% 
  st_sfc(crs = 4326) %>% 
  st_as_sf() %>% 
  mutate(area = "Maine Coastal Current Region") 
st_geometry(mcc_turnoff_poly) <- "geometry"


# Map it
ggplot() +
  geom_sf(data = maine) +
  geom_sf(data = mcc_turnoff_poly, fill = "orange", alpha = 0.6) +
  theme_classic() +
  map_theme() +
  labs(title = "Area Coverage of MCC Turnoff FVCOM Input")

Code 
# Trim Mesh to that area

# Project polygon
poly_projected <- st_transform(mcc_turnoff_poly, st_crs(gom3_mesh)) 

# Turn off s2
sf::sf_use_s2(FALSE)

# Trim it
mcc_nodes <- mesh_trim(
  mesh = gom3_mesh, 
  domain = st_as_sf(poly_projected) )

# Map that
ggplot() +
  geom_sf(data = maine) +
  geom_sf(data = mcc_nodes, color = "black", alpha = 0.6) +
  geom_sf(data = mcc_turnoff_poly, fill = "orange", alpha = 0.6) +
  theme_classic() +
  map_theme() +
  coord_sf(xlim = c(-70.2, -66.8), ylim = c(43, 44.75)) +
  labs(title = "Area Coverage of MCC Turnoff FVCOM Input - Mesh trimmed")

Grab the Surface Current Information

Code 
# The code below works, but has been commented out to prevent accidental repeat runs


# # All the fils for gom3
# file_list <- monthly_fvcom_inventory(
#   start_yr = 1978,
#   end_y = 2016,
#   fvcom_folder = cs_path("res", "FVCOM/monthly_means"),
#   fvcom_vers = "GOM3"
# ) 
# 
# 
# # 2 vars for water velocity
# gom3_early$var$u$longname  # Eastward Water Velocity in meters s-1
# gom3_early$var$v$longname  # Northward Water Velocity in meters s-1
# gom3_early$var$ua # Vertically Averaged x-velocity in meters s-1
# gom3_early$var$va # Vertically Averaged y-velocity in meters s-1
# 
# 
# 
# 
# # Getting node elements
# monthly_currents <- map_dfr(
#     # file_list[c(1:5)],
#     file_list,
#     possibly(
#       .f = function(fpath){
#          # Open (lazy-load) the netcdf connection
#         x <- nc_open(fpath)
#         
#         # Time dimension info
#         time_dim <- fvcom_time(x)
#         
#         # Grab surface variables at the elements
#         surface_elems <- get_elem_var(
#           x, # Dataset lazyloaded with ncdf4 from THREDDS 
#           y = 1, # integer, or indices for siglay or siglev (depth indices)
#           # var = c("u", "v"),  # Variables we want
#           var = c("u", "v", "ua", "va"),  # Variables we want
#           elem = unique(mcc_nodes$elem), 
#           time = c(1:time_dim)) %>%  # All time intervals - averaged
#           as.data.frame()
#         
#         # Close connection to the netcdf:
#         nc_close(x)
#         
#         # Join to the mesh
#         surface_mesh <- left_join(mcc_nodes, surface_elems, by = join_by(elem))
#         
#         # Re-grid to a regular-grid raster
#         surface_raster <- raster::stack(
#           sapply(
#             c("u", "v"), 
#             function(f) { 
#               fvcom::rasterize(
#                 x = surface_mesh, 
#                 template = default_template(surface_mesh, res = c(0.05,0.05)),
#                 field = f) }, 
#             simplify = FALSE))
#         
#         # Make it a table to facilitate PCA easier
#         grid_df <-  surface_raster %>% 
#           raster::rasterToPoints() %>% 
#           as.data.frame()
#         return(grid_df)}, 
#       
#       # Spit out NA's if there's trouble
#       otherwise = data.frame(
#         x = NULL,
#         y = NULL,
#         u = NULL,
#         v = NULL)), 
#     .id = "date")


# # Save this out:
# write_csv(monthly_currents, here::here("data/gom3_mcc_turnoff_input_data.csv"))

Checking Failed Files

Some of the files timed out during the download and were empty, this is a targeted chunk of code for re-downloading and re-cropping for those months.

Code 
# The code below works, but has been commented out to prevent accidental repeat runs


# # Need to download these again b/c zero bytes
# missing_months <- c(
#   198109,
#   198303,
#   198502,
#   198610,
#   198705,
#   198707,
#   199212 )



# #### Download Fresh Copies:  ####
#   
# # Assemble year/month structure for file names
#   
# # Build THREDDS Link Structure
# seaplan_hcast <- "http://www.smast.umassd.edu:8080//thredds/fileServer/models/fvcom/NECOFS/Archive/Seaplan_33_Hindcast_v1/monthly_mean/"
# gom3_base <- "gom3_monthly_mean_"
# dest_folder <- cs_path("res", "FVCOM/monthly_means/gom3_mon_means")
# 
# # File names
# fnames <- str_c(gom3_base, missing_months, ".nc")
# 
# Now step through each one and download/save
# purrr::walk(fnames, function(file_name){
# 
#   # Build the download url and out paths for saving
#   url_full <- str_c(seaplan_hcast, file_name)
#   save_full <- str_c(dest_folder, file_name)
# 
#   # Download and save
#   message(str_c("Downloading: ", file_name))
#   # message(str_c("from: ", url_full))
#   # message(str_c("at: ", save_full))
#   download.file(
#     url = url_full,
#     destfile = save_full)
# })





#### Get Timeseries for those months  ####
# new_files <- str_c("/Users/akemberling/Library/CloudStorage/Box-Box/RES_Data/FVCOM/monthly_means/gom3_mon_means/gom3_monthly_mean_", missing_months, ".nc")
# new_files <- setNames(new_files, missing_months)
# 
# # Run them through
# missing_currents <- map_dfr(
#     new_files,
#     possibly(
#       .f = function(fpath){
#          # Open (lazy-load) the netcdf connection
#         x <- nc_open(fpath)
#         
#         # Time dimension info
#         time_dim <- fvcom_time(x)
#         
#         # Grab surface variables at the elements
#         surface_elems <- get_elem_var(
#           x, # Dataset lazyloaded with ncdf4 from THREDDS 
#           y = 1, # integer, or indices for siglay or siglev (depth indices)
#           var = c("u", "v"),  # Variables we want
#           elem = unique(mcc_nodes$elem), 
#           time = c(1:time_dim)) %>%  # All time intervals - averaged
#           as.data.frame()
#         
#         # Close connection to the netcdf:
#         nc_close(x)
#         
#         # Join to the mesh
#         surface_mesh <- left_join(mcc_nodes, surface_elems, by = join_by(elem))
#         
#         # Re-grid to a regular-grid raster
#         surface_raster <- raster::stack(
#           sapply(
#             c("u", "v"), 
#             function(f) { 
#               fvcom::rasterize(
#                 x = surface_mesh, 
#                 template = default_template(surface_mesh, res = c(0.05,0.05)),
#                 field = f) }, 
#             simplify = FALSE))
#         
#         # Make it a table to facilitate PCA easier
#         grid_df <-  surface_raster %>% 
#           raster::rasterToPoints() %>% 
#           as.data.frame()
#         return(grid_df)}, 
#       
#       # Spit out NA's if there's trouble
#       otherwise = data.frame(
#         x = NULL,
#         y = NULL,
#         u = NULL,
#         v = NULL)), 
#     .id = "date")


# # Join these to the rest of the months that worked before:
# all_gom3_currents <- bind_rows(mutate(monthly_currents, date = as.character(date)), missing_currents)



# Now re-save
# write_csv(all_gom3_currents, here::here("data/gom3_mcc_turnoff_input_data.csv"))

Visualizing Current Vectors

Code 
# Load it and plot a year (so we can verify it looks reasonable)
monthly_currents <- read_csv(here::here("data/gom3_mcc_turnoff_input_data.csv")) %>% 
  rename(label = date) %>% 
  mutate(year = str_sub(label, 1,4),
         month = as.numeric(str_sub(label, 5, 6)),
         date = as.Date(str_c(year, month, "01", sep = "-")))

# Make ID's from lat/lon pairs
cell_labels <- monthly_currents %>% 
  distinct(x,y) %>% 
  mutate(cell_id = row_number())

# Use them + variables to reshape into a matrix
monthly_currents <- right_join(cell_labels, monthly_currents)


 # Get xy coordinates
 # all_gom3_currents %>% 
 monthly_currents %>% 
   filter(str_sub(label, 1, 4) == 2000) %>% 
   ggplot() + 
    metR::geom_vector(
    aes(x, y = y, dx = u, dy = v), 
    arrow.angle = 30, arrow.type = "open", arrow.length = .25, 
    pivot = 0, 
    preserve.dir = TRUE, 
    direction = "ccw", 
    color = "gray20", 
    alpha = 0.85)+
    scale_fill_distiller(palette = "RdBu") +
   theme_classic() +
  map_theme() +
   facet_wrap(~date)

Performing PCA

Each lat/lon pair can be given a unique ID. We can then pivot wider so that ID-variable matches each get one row per month. Then we can perform the PCA on these variables.

We can do this on all the current variables, or just one: https://github.com/dzaugis/Ecosystem_Indicators/blob/6d21e553614cb06eb7ea02e4546535cf038d7678/Code/MCC_index_report.Rmd#L69-L76

Code 
# Use them + variables to reshape into a matrix
pca_mat <- monthly_currents %>% 
  select(-x, -y, -year, -month) %>% 
  pivot_wider(names_from = cell_id, values_from = c(u, v)) %>% 
  column_to_rownames("date")

# # One row per month - check
# dim(pca_mat)[1] ==length(unique(monthly_currents$label))


# Do PCA
mcc_pca <- prcomp(pca_mat, scale. = TRUE, center = TRUE)

# Pull Principal Components
mcc_pca_pc1 <- data.frame(mcc_pca$x[,1]) %>% 
  rownames_to_column( "date") %>% 
  rename("PC1" = mcc_pca.x...1.)
mcc_pca_pc2 <- data.frame(mcc_pca$x[,2]) %>% 
  rownames_to_column("date") %>% 
  rename( "PC2" = mcc_pca.x...2.) 
mcc_pca_both <- left_join(mcc_pca_pc1, mcc_pca_pc2, by= "date") %>% 
  mutate(date = as.Date(date), 
         mon = lubridate::month(date), 
         yr = lubridate::year(date))

# Pull Loadings
mcc_loadings.pc1 <- data.frame(mcc_pca$rotation[,1])  %>% 
  rownames_to_column("loc") %>% 
  rename("PC1" = mcc_pca.rotation...1.)
mcc_loadings.pc2 <- data.frame(mcc_pca$rotation[,2]) %>% 
  rownames_to_column("loc") %>% 
  rename("PC2" = mcc_pca.rotation...2.) 

# reshape loadings
mcc_loadings <- left_join(mcc_loadings.pc1, mcc_loadings.pc2, by= "loc") %>% 
  separate(loc, sep = "_", into = c("var", "cell_id")) %>% 
  filter(var != "label") %>% 
  pivot_longer(cols=c(PC1, PC2), names_to = "PC", values_to = "values") %>%
  mutate(plots = paste(PC, var, sep = "_"))

Plotting Correlations in Space

Code 
# little bonus code for direction and velocity
mcc_corrs <- monthly_currents %>% 
  left_join(mcc_pca_both, by = "date") %>% 
  filter(year > 1980) %>% 
  mutate(dir = ifelse(
    REdaS::rad2deg(atan2(v,u)) < 0, 
    REdaS::rad2deg(atan2(v,u)) + 360, 
    REdaS::rad2deg(atan2(v,u))), 
    vel = sqrt(u^2+v^2)) %>% 
  group_by(x, y) %>% 
  summarise(
    PC1_vel = cor(vel, PC1),
    PC2_vel = cor(vel, PC2),
    PC1_dir = cor(dir, PC1),
    PC2_dir = cor(dir, PC2),
    PC1_u   = cor(u, PC1),
    PC2_u   = cor(u, PC2),
    PC1_v   = cor(v, PC1),
    PC2_v   = cor(v, PC2)) %>% 
  pivot_longer(
    cols = c(-y, -x), 
    names_to = "Cor", 
    values_to = "Correlation")


# Map correlations between Principal Components and u/v
mcc_corrs %>% 
  filter(Cor %in% c("PC1_u", "PC2_u", "PC1_v", "PC2_v")) %>% 
  ggplot() + 
  geom_sf(data = maine) +
  geom_raster(aes(x, y, fill = Correlation)) + 
  scale_fill_distiller(palette = "RdBu") + 
  facet_wrap(~Cor) + 
  theme_classic() +
  map_theme() +
  coord_sf(xlim = c(-70.8, -66.7), ylim = c(43, 44.8), expand = F,
           datum = "+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0")

Vertical Current Velocity?

(It seems like Matt used the vertical movement velocity info in the PCA as well, here is how he grabbed it from the netcdfs, one value each time step)[https://github.com/dzaugis/Ecosystem_Indicators/blob/6d21e553614cb06eb7ea02e4546535cf038d7678/Code/FVCOM_shp_extract.Rmd#L254C1-L279C2]

Code 
# # Need u vert "ua" for this
# monthly_currents %>% 
#  filter(date %in% c(as.Date("2010-03-01"), as.Date("1986-03-01"), as.Date("2011-03-01"))) %>% 
#   mutate(
#     vel = sqrt(u_vert^2 + v_vert^2), 
#     PC = if_else(yr == 2010, "positive PC1", if_else(yr ==1986, "negative PC1", "neutral PC1"))) %>% 
#   ggplot() + geom_sf(data= ne_us, fill = "grey") + 
#   geom_segment(aes(x = lon, y = lat, xend=lon+u, yend=lat+v, color = vel), 
#                arrow = arrow(angle = 30, length = unit(0.05, "inches"), type = "closed")) + 
#   scale_color_viridis_c() + theme_bw() + 
#   coord_sf(datum = "+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0") + facet_wrap(~PC)+ 
#   theme(panel.background = element_blank(), panel.grid = element_blank(), 
#         axis.title = element_blank(), axis.text = element_blank(),
#         axis.ticks = element_blank())