Polygon preparation

Cabot Strait

It was decided that because we are interested in the geography of the regions, and not just their temperature regimes, the Cabot Strait needed to be defined apart from the Gulf of St. Lawrence region. To do this we will simply snip the “gsl” polygon into two pieces at its narrowest point.

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

It appears from the crude figure above that we should pinch the polygon off into two separate shapes at row 6 and 10.

# 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

Everything is looking good, but we may want to divide the Gulf of Maine (gm) into two polygons as well. This would make the Bay of Fundy in it’s own region. For now however we will move on to the next step, which is dividing the current polygons by bathymetry.

Downloading data

Before we can divide up our polygons by bathymetry, we must first download said bathymetry data. We will use here the NOAA data as it is (to my knowledge) the most convenient bathymetry data available directly through R.

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

With our high res bathymetry data downloaded we will now go about finding the 50 and 200 metre isobaths. These lines will be used to define the edges of our sub-regions. The naming conventions we will use will be

# 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

The resolution on the figure above is fantastic, but perhaps a bit more than what we need as the desire is to create functional polygons for dividing our regions. The level of detail present in the figure above will likely lead to issues later on. Let’s look at the lowres bathymetry data and see if that works better.

# 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

The above plot shows the NOAA bathymetry data downloaded with a resolution of 6. I tried these data at all resolutions from 2 – 8 and visually this looks to be the best balance between removing unnecessary convolutions without loosing the broader features.

Isobath contour based polygons

With a decent resolution picked out for our bathymetry data we now nee to create polygons from the desired isobaths. These will then be used as custom bounding boxes with which we will filter and combine the SST pixels for the next vignette.

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

With the isobath contours now convert to a spatial class object we will need to do the same for our study area polygons so they can play nice with one another.

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

After going through all of this shape file nonsense it dawned on me that I could much more easily just have used the bathymetry data already present in the NAPA model and then clumped the SST pixels based on which of the depth categories they fall within.

We will now go about creating SST time series for each of the depth sub-regions. This work is continued in the SST preparation vignette.