Polygon preparation

Introduction

This markdown file contains all of the code used to prepare the polygons that are used to define the different regions in the Northwest Atlantic. These different regions will then have their SST pixels spatially averaged to create a single time series per region. This is done so that the MHW detection algorithm may then be run on these individual time series to represent a general representation of the SST in those regions, rather than running the algorithm on each pixel individually, which would introduce a host of problems.

# Packages used in this vignette
library(jsonlite, lib.loc = "../R-packages/")
library(tidyverse) # Base suite of functions
library(R.matlab) # For dealing with MATLAB files
library(marmap) # For bathymetry
library(maptools) # Contour tools
library(rgeos) # For intersections

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 kind-hearted man that he is, Benjamin forwarded us the polygons (Richaud et al. 2016; Figure 2) from this work as a MATLAB file so we must first open that up and convert it to an R format for further use. It should be noted that these areas were designed to not encompass depths deeper than 600 m as the investigators were itnerested in characterising the climatolgoies for the shelf and upper slope regions of the north east coast of North America.

# Load the file
NWA_polygons <- R.matlab::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
# vec <- NWA_polygons[1]
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 <- rbind(coords_1, coords_2, coords_3, coords_4, coords_5, coords_6)
colnames(NWA_coords) <- c("region", "lon", "lat")
# saveRDS(NWA_coords, "data/NWA_coords.Rda")

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

# Load polygon coordinates
NWA_coords <- readRDS("data/NWA_coords.Rda")

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

# Quick map
NWA_coords_plot <- ggplot(data = NWA_coords, aes(x = lon, y = lat)) +
  geom_polygon(data = map_base, aes(group = group), show.legend = F) +
  geom_polygon(aes(colour = region, fill = region), size = 1.5, alpha = 0.2) +
  coord_cartesian(xlim = c(min(NWA_coords$lon)-10, max(NWA_coords$lon)+10),
                  ylim = c(min(NWA_coords$lat)-10, max(NWA_coords$lat)+10)) +
  labs(x = NULL, y = NULL, colour = "Region", fill = "Region") +
  theme(legend.position = "bottom")
# ggsave(NWA_coords_plot, filename = "output/NWA_coords_plot.pdf", height = 5, width = 6)

# Visualise
NWA_coords_plot

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.

Before we move on, we’ll do a bit of house keeping to establish a consistent study area for this project based on our polygons. We’ll simply extend the study area by 10 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.

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

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

# Extract the gsl region only
gsl_sub <- NWA_coords[NWA_coords$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)

# Load the corners of the study area
NWA_corners <- readRDS("data/NWA_corners.Rda")

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

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

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

# Plot the new areas to ensure everything worked
NWA_coords_cabot_plot <- ggplot(data = NWA_coords_cabot, aes(x = lon, y = lat)) +
  geom_polygon(data = map_base, aes(group = group), show.legend = F) +
  geom_polygon(aes(colour = region, fill = region), size = 1.5, alpha = 0.2) +
  coord_cartesian(xlim = NWA_corners[1:2],
                  ylim = NWA_corners[3:4]) +
  labs(x = NULL, y = NULL, colour = "Region", fill = "Region") +
  theme(legend.position = "bottom")
# ggsave(NWA_coords_cabot_plot, filename = "output/NWA_coords_cabot_plot.pdf", height = 5, width = 6)

# Visualise
NWA_coords_cabot_plot

Bathymetry polygons

These regions are an excellent start, but because the aim of this research project is to determine the primary drivers of MHWs in the Northwest Atlantic we want to subset the regions not only by their thermal characteristics, but also by their geophysical characteristics, which in this instance is most immediately available as bathymetry. The rationale for this is that surface forcing as a driver of anomalously warm seawater is likely to be more prevalent in shallow waters than in shelf waters or deeper. To this end we will divide the regions into three sub-regions each: 0 – 50 m depth, 51 – 200 m depth, 201+ m depth.

# Donwload NOAA bathy data at highest available resolution
NWA_bathy_hires <- as.xyz(getNOAA.bathy(lon1 = NWA_corners[1]+9, lon2 = NWA_corners[2]-9,
                                        lat1 = NWA_corners[3]+9, lat2 = NWA_corners[4]-9, resolution = 1))
colnames(NWA_bathy_hires) <- c("lon", "lat", "depth")
NWA_bathy_hires <- NWA_bathy_hires %>% 
  filter(depth <= 0) %>% 
  mutate(lon = round(lon, 4),
         lat = round(lat, 4))
# saveRDS(NWA_bathy_hires, file = "data/NWA_bathy_hires.Rda")

# Donwload NOAA bathy data at a courser resolution
NWA_bathy_lowres <- as.xyz(getNOAA.bathy(lon1 = NWA_corners[1]+9, lon2 = NWA_corners[2]-9,
                                         lat1 = NWA_corners[3]+9, lat2 = NWA_corners[4]-9, resolution = 6))
colnames(NWA_bathy_lowres) <- c("lon", "lat", "depth")
NWA_bathy_lowres <- NWA_bathy_lowres %>% 
  filter(depth <= 0) %>% 
  mutate(lon = round(lon, 4),
         lat = round(lat, 4))
# saveRDS(NWA_bathy_lowres, file = "data/NWA_bathy_lowres.Rda")

# Load hires bathymetry data
NWA_bathy_hires <- readRDS("data/NWA_bathy_hires.Rda")

# Plot to see what we've got
bathy_hires_plot <- ggplot(NWA_coords_cabot, aes(x = lon, y = lat)) +
  geom_polygon(data = map_base, aes(group = group), show.legend = F) +
  geom_polygon(aes(group = region, fill = region, colour = region), alpha = 0.2) +
  geom_contour(data = NWA_bathy_hires, aes(z = depth),
               breaks = c(-50),  size = c(0.3),  colour = "grey70") +
    geom_contour(data = NWA_bathy_hires, aes(z = depth),
               breaks = c(-200),  size = c(0.3),  colour = "grey30") +
  coord_cartesian(xlim = NWA_corners[1:2],
                  ylim = NWA_corners[3:4]) +
  labs(x = NULL, y = NULL)

# Visualise
bathy_hires_plot

# Load hires bathymetry data
NWA_bathy_lowres <- readRDS("data/NWA_bathy_lowres.Rda")

# Plot to see what we've got
bathy_lowres_plot <- ggplot(NWA_coords_cabot, aes(x = lon, y = lat)) +
  geom_polygon(data = map_base, aes(group = group), show.legend = F) +
  geom_polygon(aes(group = region, fill = region, colour = region), alpha = 0.2) +
  geom_contour(data = NWA_bathy_lowres, aes(z = depth),
               breaks = c(-50),  size = c(0.3),  colour = "grey70") +
    geom_contour(data = NWA_bathy_lowres, aes(z = depth),
               breaks = c(-200),  size = c(0.3),  colour = "grey30") +
  coord_cartesian(xlim = NWA_corners[1:2],
                  ylim = NWA_corners[3:4]) +
  labs(x = NULL, y = NULL)

# Visualise
bathy_lowres_plot

# Spread bathy dataframe to matrix to play nice with shapefiles
NWA_bathy_matrix <- NWA_bathy_lowres %>% 
  mutate(depth = round(depth, -1)) %>% 
  reshape2::acast( lon~lat, value.var = "depth")

# Generate contours
cont <- contourLines(x = as.numeric(row.names(NWA_bathy_matrix)), y = as.numeric(colnames(NWA_bathy_matrix)), 
                     z = NWA_bathy_matrix, levels = c(-50, -200))
cont_lines <- ContourLines2SLDF(cont, proj4string = CRS("+proj=longlat +datum=WGS84")) # create SpatialLines DF

# Have a look
plot(cont_lines)

# First convert the dataframe of different areas to a list
poly_list <- split(NWA_coords_cabot, NWA_coords_cabot$region)

# We then only want lon/lat in the list, not the names
poly_list <- lapply(poly_list, function(x) { x["region"] <- NULL; x })

# Convert coords to polygon type
poly_poly <- sapply(poly_list, Polygon)

# Add id variables back in
poly_poly <- lapply(seq_along(poly_poly), function(i) Polygons(list(poly_poly[[i]]), 
                                            ID = names(poly_list)[i]))

# Create SpatialPolygons object
spatial_poly <- SpatialPolygons(poly_poly, proj4string = CRS("+proj=longlat +datum=WGS84"))

# Clip the contour lines to the study area polygons
spatial_poly_lines <- gIntersection(cont_lines, spatial_poly) # clip the contours to subarea0

# Visualise
plot(spatial_poly)
lines(spatial_poly_lines)

Session information

sessionInfo()

R version 3.6.0 (2019-04-26)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.5 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  rgeos_0.3-28    maptools_0.9-4  sp_1.3-1       
 [5] marmap_1.0.1    R.matlab_3.6.1  forcats_0.3.0   stringr_1.3.1  
 [9] dplyr_0.7.6     purrr_0.2.5     readr_1.1.1     tidyr_0.8.1    
[13] tibble_1.4.2    ggplot2_3.0.0   tidyverse_1.2.1 jsonlite_1.6   

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