Study area and regions

Introduction

This markdown file contains all of the code that prepares the polygons used to define the different regions in the Northwest Atlantic. These different regions then have the NOAA OISST pixels found therein spatially averaged to create a single time series per region in the SST preparation vignette. This is done so that the MHW detection algorithm may then be run on these individual time series as a general representation of the SST in those regions, rather than running the algorithm on each pixel individually, which would introduce a host of technical and philosophical problems that I won’t go into here.

# Packages used in this vignette\
library(tidyverse) # Base suite of functions
library(R.matlab) # For dealing with MATLAB files
# library(marmap) # For bathymetry, but not currently used

Coastal region polygons

The first step in this analysis is to broadly define the coastal regions based on previous research into thermally relevant boundaries. We have chosen to use a paper by Richaud et al. (2016) to do this (https://www.sciencedirect.com/science/article/pii/S0278434316303181#f0010). Being the wonderful person that he is, Benjamin forwarded us the polygons [Richaud et al. (2016); Figure 2] from this paper. The only hiccup being that they are a MATLAB file so we must first convert them to an R format. It should be noted that these areas were designed to not encompass depths deeper than 600 m as the investigators were interested in characterising the climatologies for the shelf and upper slope regions of the north east coast of North America. This works for our research purposes as well.

# Load the file
NWA_polygons <- readMat("data/boundaries.mat")

# Remove index list items and attributes
NWA_polygons[grepl("[.]",names(NWA_polygons))] <- NULL
# attributes(NWA_polygons) <- NULL

# Function for neatly converting list items into a dataframe
mat_col <- function(vec){
  df <- as.data.frame(vec)
  df$region <- substr(colnames(df)[1], 2, nchar(colnames(df)[1]))
  colnames(df)[1] <- strtrim(colnames(df)[1], 1)
  df <- df[c(2,1)]
  return(df)
}

# Create multiple smaller data.frames
coords_1 <- cbind(mat_col(NWA_polygons[1]), mat_col(NWA_polygons[2])[2])
coords_2 <- cbind(mat_col(NWA_polygons[3]), mat_col(NWA_polygons[4])[2])
coords_3 <- cbind(mat_col(NWA_polygons[5]), mat_col(NWA_polygons[6])[2])
coords_4 <- cbind(mat_col(NWA_polygons[7]), mat_col(NWA_polygons[8])[2])
coords_5 <- cbind(mat_col(NWA_polygons[9]), mat_col(NWA_polygons[10])[2])
coords_6 <- cbind(mat_col(NWA_polygons[11]), mat_col(NWA_polygons[12])[2])

# Combine them into one full dataframe and save
NWA_coords_base <- rbind(coords_1, coords_2, coords_3, coords_4, coords_5, coords_6)
colnames(NWA_coords_base) <- c("region", "lon", "lat")

With our polygons switched over from MATLAB to R we now want to visualise them to ensure that everything has gone smoothly.

# The base map
map_base <- ggplot2::fortify(maps::map(fill = TRUE, col = "grey80", plot = FALSE)) %>%
  dplyr::rename(lon = long) %>%
  mutate(group = ifelse(lon > 180, group+9999, group),
         lon = ifelse(lon > 180, lon-360, lon)) %>% 
  select(-region, -subregion)
# saveRDS(map_base, "data/map_base.Rda")

# Quick map
ggplot(data = NWA_coords_base, aes(x = lon, y = lat)) +
  geom_polygon(aes(colour = region, fill = region), size = 1.5, alpha = 0.2) +
  geom_polygon(data = map_base, aes(group = group), show.legend = F) +
  coord_cartesian(xlim = c(min(NWA_coords_base$lon)-2, max(NWA_coords_base$lon)+2),
                  ylim = c(min(NWA_coords_base$lat)-2, max(NWA_coords_base$lat)+2)) +
  labs(x = NULL, y = NULL, colour = "Region", fill = "Region") +
  theme(legend.position = "bottom")

The region abbreviations are: “gm” for Gulf of Maine, “gls” for Gulf of St. Lawrence, “ls” for Labrador Shelf, “mab” for Mid-Atlantic Bight, “nfs” for Newfoundland Shelf and “ss” for Scotian Shelf.

# Extract the gsl region only
gsl_sub <- NWA_coords_base[NWA_coords_base$region == "gsl",]

# Add a simple integer column for ease of plotting
gsl_sub$row_count <- 1:nrow(gsl_sub)

ggplot(data = gsl_sub, aes(x = lon, y = lat)) +
  geom_polygon(aes(fill = region)) +
  geom_label(aes(label = row_count)) +
  labs(x = NULL, y = NULL)

# Create smaller gsl polygon
gsl_new <- NWA_coords_base[NWA_coords_base$region == "gsl",] %>% 
  slice(-c(7:9))

# Create new cbs (Cabot Strait) polygon
cbs <- NWA_coords_base[NWA_coords_base$region == "gsl",] %>% 
  slice(6:10) %>% 
  mutate(region = "cbs")

# Attach the new polygons to the original polygons
NWA_coords_cabot <- NWA_coords_base %>% 
  filter(region != "gsl") %>% 
  rbind(., gsl_new, cbs)
# saveRDS(NWA_coords_cabot, "data/NWA_coords_cabot.Rda")

# Plot the new areas to ensure everything worked
ggplot(data = NWA_coords_cabot, aes(x = lon, y = lat)) +
  geom_polygon(aes(colour = region, fill = region), size = 1.5, alpha = 0.2) +
  geom_polygon(data = map_base, aes(group = group), show.legend = F) +
  coord_cartesian(xlim = c(min(NWA_coords_cabot$lon)-2, max(NWA_coords_cabot$lon)+2),
                  ylim = c(min(NWA_coords_cabot$lat)-2, max(NWA_coords_cabot$lat)+2)) +
  labs(x = NULL, y = NULL, colour = "Region", fill = "Region") +
  theme(legend.position = "bottom")

# FIlter out the ls region
NWA_coords <- NWA_coords_cabot %>% 
  filter(region != "ls")

# Save these final study region coordinates
# saveRDS(NWA_coords, "data/NWA_coords.Rda")

# Plot the new areas to ensure everything worked
NWA_study_area <- ggplot(data = NWA_coords, aes(x = lon, y = lat)) +
  geom_polygon(aes(colour = region, fill = region), size = 1.5, alpha = 0.2) +
  geom_polygon(data = map_base, aes(group = group), show.legend = F) +
  coord_cartesian(xlim = c(min(NWA_coords$lon)-2, max(NWA_coords$lon)+2),
                  ylim = c(min(NWA_coords$lat)-2, max(NWA_coords$lat)+0.5),
                  expand = FALSE) +
    scale_x_continuous(breaks = seq(-70, -50, 10),
                     labels = c("70°W", "60°W", "50°W"),
                     position = "top") +
  scale_y_continuous(breaks = c(40, 50),
                     labels = scales::unit_format(suffix = "°N", sep = "")) +
  labs(x = NULL, y = NULL, colour = "Region", fill = "Region") +
  theme_bw() +
  theme(legend.position = "bottom")
# ggsave(NWA_study_area, filename = "output/NWA_study_area.pdf", height = 5, width = 6)
# ggsave(NWA_study_area, filename = "output/NWA_study_area.png", height = 5, width = 6)

# Visualise
NWA_study_area

Study area extent

The final step in this vignette is to establish a consistent study area for this project based on our regions. We’ll simply extend the study area by the nearest 2 whole degrees of longitude and latitude from the furthest edges of the polygons, as seen in the figure above. This will encompass broad synoptic scale variables that may be driving MHWs in our study regions, but should not be so broad as to begin to account for teleconnections, which are currently beyond the scope of this project. The exception to this being that because we want to exclude as much of the Labrador Sea as possible, we will only extend the northern edge of the study area by 0.5 degrees of latitude from the northernmost point of our study regions.

# Set the max/min lon/at values
lon_min <- round(min(NWA_coords$lon)-2)
lon_max <- round(max(NWA_coords$lon)+2)
lat_min <- round(min(NWA_coords$lat)-2)
lat_max <- round(max(NWA_coords$lat)+0.5)

# Combine and save
NWA_corners <- c(lon_min, lon_max, lat_min, lat_max)
# saveRDS(NWA_corners, file = "data/NWA_corners.Rda")

# Pront the values
NWA_corners

[1] -80 -41  32  52

We will now go about creating SST time series for each of the regions and calculate marine heatwaves (MHWs) from them. This work is continued in the SST preparation and MHW detection vignette.

References

Session information

sessionInfo()

R version 3.6.1 (2019-07-05)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.6 LTS

Matrix products: default
BLAS:   /usr/lib/openblas-base/libblas.so.3
LAPACK: /usr/lib/libopenblasp-r0.2.18.so

locale:
 [1] LC_CTYPE=en_CA.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_CA.UTF-8        LC_COLLATE=en_CA.UTF-8    
 [5] LC_MONETARY=en_CA.UTF-8    LC_MESSAGES=en_CA.UTF-8   
 [7] LC_PAPER=en_CA.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_CA.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] bindrcpp_0.2.2  R.matlab_3.6.1  forcats_0.3.0   stringr_1.3.1  
 [5] dplyr_0.7.6     purrr_0.2.5     readr_1.1.1     tidyr_0.8.1    
 [9] tibble_1.4.2    ggplot2_3.0.0   tidyverse_1.2.1

loaded via a namespace (and not attached):
 [1] tidyselect_0.2.4  haven_1.1.2       lattice_0.20-35  
 [4] colorspace_1.3-2  htmltools_0.3.6   yaml_2.2.0       
 [7] rlang_0.2.2       R.oo_1.22.0       pillar_1.3.0     
[10] glue_1.3.0        withr_2.1.2       R.utils_2.7.0    
[13] modelr_0.1.2      readxl_1.1.0      bindr_0.1.1      
[16] plyr_1.8.4        munsell_0.5.0     gtable_0.2.0     
[19] workflowr_1.1.1   cellranger_1.1.0  rvest_0.3.2      
[22] R.methodsS3_1.7.1 evaluate_0.11     labeling_0.3     
[25] knitr_1.20        broom_0.5.0       Rcpp_0.12.18     
[28] scales_1.0.0      backports_1.1.2   jsonlite_1.5     
[31] hms_0.4.2         digest_0.6.16     stringi_1.2.4    
[34] grid_3.6.1        rprojroot_1.3-2   cli_1.0.0        
[37] tools_3.6.1       maps_3.3.0        magrittr_1.5     
[40] lazyeval_0.2.1    crayon_1.3.4      whisker_0.3-2    
[43] pkgconfig_2.0.2   xml2_1.2.0        lubridate_1.7.4  
[46] assertthat_0.2.0  rmarkdown_1.10    httr_1.3.1       
[49] rstudioapi_0.7    R6_2.2.2          nlme_3.1-137     
[52] git2r_0.23.0      compiler_3.6.1