Calculating Warm Core Ring Metrics

Setting up your environment

Load libraries

## Wrangling
library(ncdf4)
library(dplyr)
library(tidyverse)
## Diagnostics
library(car)          ## Anova, influence plot
library(patchwork)    ## multi-panel figures
library(swfscMisc)    ## create circular spatial object
## Extraction/graphics
library(geosphere)
library(ggplot2)
## From GitHub

Import data

df <- read.csv('ring_data.csv')
df<- df[,-1]
nafo <- read.csv('NAFO_timeshelf_weeks.csv')
sys_state <- read.csv('sys_states.csv')
bathy <- nc_open('SRTM30-NESGRID-BATHY-SRTM30_PLUS.nc')
# Convert Occupancy data to long format
occupancy <- gather(nafo, id, occupancy, z6C:z5Ze, factor_key = TRUE)

Modify data

The following code sets up our dataframe with place holders for metrics that will be calculated and included

df$id <- as.character(df$id)
ringnames <- unique(df$id)
# Initiate columns to add
df$weeks_old <- NA
df$days_old <- NA
df$subarea <- NA # nafo
df$score <- NA # system state info
# Open bathymetry file 
lat <- ncvar_get(bathy, 'latitude')
lon <- ncvar_get(bathy, 'longitude')
isobath <- ncvar_get(bathy, attributes(bathy$var)$names[1])
fillvalue <- ncatt_get(bathy, 'depth', "_fillvalue")
isobath[isobath == fillvalue$value] <- NA 
nc_close(bathy) 
grid <- expand.grid(lon = lon, lat = lat)
isobath[(isobath > -200) | (isobath < -200)] <- NA

Add in information about system state as defined in Paul Rago’s work (2020)

locs1 <- which(df$year == 2011)
df[locs1,12] <- 'A'
locs2 <- which(df$year == 2012)
df[locs2,12] <- 'A'
locs3 <- which(df$year == 2013)
df[locs3,12] <- 'P'
locs4 <- which(df$year == 2014)
df[locs4,12] <- 'A'
locs5 <- which(df$year == 2015)
df[locs5,12] <- 'P'
locs6 <- which(df$year == 2016)
df[locs6,12] <- 'P'
locs7 <- which(df$year == 2017)
df[locs7,12] <- 'G'
locs8 <- which(df$year == 2018)
df[locs8,12] <- 'G'
locs9 <- which(df$year == 2019)
df[locs9,12] <- 'G'
locs10 <- which(df$year == 2020)
df[locs10,12] <- 'G'

Calculate first metric : Ages of rings

This bit of code calculates the age of each ring in terms of weeks old and days old.

These new metrics are then added to the original data frame, ‘df’. From this data frame, the age of demise is calculated by determining the maximum days old for each ring. This information is stored in a new object called ‘demise’. Examining the first few rows of our updated data frame ‘df’:

for(i in 1:length(ringnames)){
  tmp2 <- subset(df, id == ringnames[i])
  for(j in 1:length(tmp2[,1])){
    tmp2[j, 9] <-  (j)
    tmp2[j, 10] <- ((tmp2[j,9]+1) - tmp2[1, 9])*7
    locs <- which(df$id == ringnames[i])
    df$weeks_old[locs] <- tmp2[,9]
    df$days_old[locs] <- tmp2[,10]
  }
 }


# calculate the max days old for each ring, age of demise == max(days_old)

# set up matrix to hold only ring information pertaining to date of demise 
demise <- matrix(nrow = length(ringnames), ncol = ncol(df)) 

for (i in 1:length(ringnames)){
  ring <- subset(df, id == ringnames[i]) # create a subset for each ring
  dtmp <- subset(ring, weeks_old == max(weeks_old)) # isolate last week 
  demise[i,] <- as.matrix(dtmp) # place into holding matrix
}
colnames(demise) <- colnames(df) # name columns
demise <- as.data.frame(demise)

head(df)

##      id year week       date    lon   lat     area   radius weeks_old days_old
## 1 T2010 2011    1 2011-01-03 -62.85 41.98 17359.99 74.33610         1        7
## 2 T2010 2011    2 2011-01-10 -62.85 42.01 18133.01 75.97313         2       14
## 3 T2010 2011    3 2011-01-17 -62.88 42.00 17683.99 75.02659         3       21
## 4 T2010 2011    4 2011-01-24 -62.86 41.99 17830.87 75.33752         4       28
## 5 T2010 2011    5 2011-01-31 -62.91 42.02 18063.17 75.82668         5       35
## 6 T2010 2011    6 2011-02-07 -62.85 42.01 17291.07 74.18840         6       42
##   subarea score
## 1      NA     A
## 2      NA     A
## 3      NA     A
## 4      NA     A
## 5      NA     A
## 6      NA     A

nrow(demise)

## [1] 281

Check that the new data frame has the same number of rows as there are unique rings (281) in the data set: r{unique(df$id) == nrow(demise)}

Focus on rings West of -65 the fishing season

The following lines of code subset the data based on a latitudinal and temporal range that best describes the geographical area and weeks during which fishing typically occurs.

# all rings that exist west of -65
shelf_rings1 <-  subset(df, lon <= -65.00)
# all rings that exist during fishing weeks and are west of -65
shelf_rings_fs1 <- subset(shelf_rings1, week %in% c(18:44))
# (fs == fishing season)

Examine differences in the number of rings across each of the categories we have just created

(i) Total number of rings:

## [1] 281

(ii) Number of rings on shelf:

## [1] 123

(iii) Number of rings on shelf during the fishing season:

## [1] 99

Rings that exist west of -65 represent 35 percent of total rings.

Below is a visual representation of the number of rings West of -65 degrees in the years 2011 through 2020.

The chunk of code summarizes the color of each bar represents the ‘score’ of the system state that year was ascribed via Paul Rago’s (2020) analysis. Note, the extent of P. Rago’s analysis is 2011:2019 (excluding 2020).

## Anova Table (Type III tests)
## 
## Response: number
##              Sum Sq Df F value    Pr(>F)    
## (Intercept) 2640.33  1 45.5979 0.0002644 ***
## score        129.57  2  1.1188 0.3787817    
## Residuals    405.33  7                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Anova Table (Type III tests)
## 
## Response: days_old
##                 Sum Sq  Df F value    Pr(>F)    
## (Intercept)      74504   1 23.4763 2.131e-06 ***
## as.factor(year)  55757   9  1.9521   0.04507 *  
## Residuals       860039 271                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Visualize: age data

The following plots look at the frequency of ages of rings that are near the shelf during the fishing season

Quantify

Run a quick Anova to test differences between the maximum age of ring across system states:

aov1.3 <- aov(days_old ~ score, data = demise_fs)
Anova(aov1.3, type = 'III')

## Anova Table (Type III tests)
## 
## Response: days_old
##             Sum Sq  Df F value    Pr(>F)    
## (Intercept) 187904   1 63.2237 2.839e-13 ***
## score        11896   2  2.0013    0.1385    
## Residuals   487417 164                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The variable score is not significant (p = 0.1385). This suggests is not a significant difference in the age of demise of rings between Good years and Poor years.

Here we see which specific groups means (compared with each other) are different.

Ring Occupancy

The following figures explore ring occupancy for rings existing West of -65 during the fishing season and hatching season

For the following analyses: Fishing season is classified as weeks 18:44, (spanning May through October).

Hatching season is a term used to reference the timing of mass hatching events via the two presumed cohorts that become available to the fishery during the fishing season. ** NOTE: hatching season is a term I created just for this concept **

There are assumed to be 2 cohorts each year. The first cohort has a ‘hatching season’ which peaks from approximately January through March. The following Jan/March subsets represent the paralarvae that are expected to hatch during the winter/early spring months and likely make up the majority of the spring/early summer catch, based on a 6-7 month maturation phase (Hendrickson_ea_04). The second hatching season is proposed to take place in June - July and the mature and spawning squid would then likely make up the majority of the Fall catch. This 6-7 month maturation phase, fairly consistent with other Illex species (Schwarz_Perez_12). Schwarz and colleagues have found Illex argentinus generally spend ~30 days as paralarvae, ~ 130 days as juveniles and 30-60 days as mature/spawning adults.

The following code uses these assumptions about the timing of maturation to explore WCR occupancy during paralarvae and juvenile phases - where rings may serve as transport or food resource aggregation mechanisms.

occ <- subset(occupancy, week %in% c(18:44))
sping_catch_hatch_occ2 <- subset(occupancy, week %in% c(49:53)) # Dec prev year

sping_catch_hatch_occ <- subset(occupancy, week %in% c(1:13)) # J:March
# fall_catch_hatch_occ <- subset(occupancy, week %in% c(10:17)) # 
fall_catch_hatch_occ <- subset(occupancy, week %in% c(23:30)) # 23:30 June/july
  
fs_occ <- subset(occupancy, week %in% c(18:44))
hs_occ <- subset(occupancy, week %in% c(1:30))

hs_occ_ann <- hs_occ %>%
  group_by(year, score) %>%
  summarise(n = sum(total.occ),
            mean_occ = mean(total.occ),
            sd_occ = sd(total.occ),
            max_occ = max(total.occ),
            min_occ = min(total.occ))
hs_occ_wk <- hs_occ %>%
  group_by(year,week, score) %>%
  summarise(n = sum(total.occ),
            mean_occ = mean(total.occ),
            sd_occ = sd(total.occ),
            max_occ = max(total.occ),
            min_occ = min(total.occ))
sc_hatch_occ_annd <- sping_catch_hatch_occ2 %>%
  group_by(year, score) %>%
  summarise(n = sum(total.occ),
            mean_occ = mean(total.occ),
            sd_occ = sd(total.occ),
            max_occ = max(total.occ),
            min_occ = min(total.occ))
sc_hatch_occ_wkd <- sping_catch_hatch_occ2 %>%
  group_by(year,week, score) %>%
  summarise(n = sum(total.occ),
            mean_occ = mean(total.occ),
            sd_occ = sd(total.occ),
            max_occ = max(total.occ),
            min_occ = min(total.occ))
sc_hatch_occ_ann <- sping_catch_hatch_occ %>%
  group_by(year, score) %>%
  summarise(n = sum(total.occ),
            mean_occ = mean(total.occ),
            sd_occ = sd(total.occ),
            max_occ = max(total.occ),
            min_occ = min(total.occ))
sc_hatch_occ_wk <- sping_catch_hatch_occ %>%
  group_by(year,week, score) %>%
  summarise(n = sum(total.occ),
            mean_occ = mean(total.occ),
            sd_occ = sd(total.occ),
            max_occ = max(total.occ),
            min_occ = min(total.occ))
fc_hatch_occ_wk <- fall_catch_hatch_occ %>%
  group_by(year, week, score) %>%
  summarise(n = sum(total.occ),
            mean_occ = mean(total.occ),
            sd_occ = sd(total.occ),
            max_occ = max(total.occ),
            min_occ = min(total.occ))
fc_hatch_occ_ann <- fall_catch_hatch_occ %>%
  group_by(year, score) %>%
  summarise(n = sum(total.occ),
            mean_occ = mean(total.occ),
            sd_occ = sd(total.occ),
            max_occ = max(total.occ),
            min_occ = min(total.occ))
occ_stats_week <- occ %>%
  group_by(year, week, score) %>%
  summarise(n = sum(total.occ),
            mean_occ = mean(total.occ),
            sd_occ = sd(total.occ),
            max_occ = max(total.occ),
            min_occ = min(total.occ))
occ_stats_ann <- occ %>%
  group_by(year, score) %>%
  summarise(n = sum(total.occ),
            mean_occ = mean(total.occ),
            sd_occ = sd(total.occ),
            max_occ = max(total.occ),
            min_occ = min(total.occ))
dt <- occ %>%
  dplyr::group_by(score, total.occ)%>%
  dplyr::tally()%>%
  dplyr::mutate(percent=n/sum(n))

# Generate figures 
year_plot <- ggplot(occ, aes(x = (as.factor(year)), y = total.occ, fill = score)) +
  geom_bar(stat="identity" ) +
  xlab('Year') +
  ylab('Ring days on shelf') +
  labs(title = 'Annual Occupancy During Fishing Season') +
  theme_minimal()

# Annual tally 
tyear = ggplot(data = occ_stats_ann,
       mapping = aes(x = year,
                     y = n)) +
  xlim(c(2011,2020)) +
  geom_line(color = 'black') +
  geom_point(aes(fill = factor(score)), shape = 21,  size = 4) +
  scale_x_continuous(breaks = seq(2011,2020, by = 1)) +
  xlab('Year') +
  ylab('Ring days on shelf') + 
  labs(title = 'Total Annual Ring Occupancy', fill = 'System state') +
  theme_minimal()

myear = ggplot(data = occ_stats_ann,
       mapping = aes(x = year,
                     y = mean_occ)) +
  geom_line(color = 'black') +
  geom_point(aes(fill = factor(score)), shape = 21,  size = 4) +
  scale_x_continuous(breaks = seq(2011,2020, by = 1)) +
  xlab('Year') +
  ylab('Ring days on shelf') + 
  labs(title = 'Mean Annual Ring Occupancy', fill = 'System state') +
  geom_errorbar(aes(ymin = mean_occ - sd_occ, ymax = mean_occ + sd_occ),
                width = .2, position = position_dodge(.9)) + 
  theme_minimal()


g1 <- ggplot(occ_stats_week, aes(x = week, y = max_occ, fill = score)) +
  geom_bar(stat="identity", alpha = 0.7 ) +
  xlab('Week') +
  ylab('Max Occupancy') +
  labs(title = 'Max Occupancy Across Fishing Season') +
  theme_minimal() +
  facet_wrap(~score)
g2 <- ggplot(occ_stats_week, aes(x = week, y = mean_occ, fill = score)) +
  geom_bar(stat="identity", alpha = 0.7 ) +
  xlab('Week') +
  ylab('Mean Occupancy') +
  labs(title = 'Mean Occupancy Across Fishing Season') +
  theme_minimal() +
  facet_wrap(~score)
g3 <- ggplot(sc_hatch_occ_wk, aes(x = week, y = mean_occ, fill = score)) +
  geom_bar(stat="identity", alpha = 0.7 ) +
  xlab('Week') +
  ylab('Mean Occupancy') +
  labs(title = 'Mean Ring Occupancy Jan:March Hatching Season', 
       subtitle = 'Spring catch') +
  theme_minimal() +
  facet_wrap(~score)
g4 = ggplot(data = sc_hatch_occ_ann,
       mapping = aes(x = year,
                     y = mean_occ)) +
  geom_line(color = 'black') +
  geom_point(aes(fill = factor(score)), shape = 21,  size = 4) +
  scale_x_continuous(breaks = seq(2011,2020, by = 1)) +
  xlab('Year') +
  ylab('Ring days on shelf') + 
  labs(title = 'Mean Ring Occupancy Jan:March Hatching Season', 
       subtitle = 'Spring catch',
       fill = 'System state') +
  geom_errorbar(aes(ymin = mean_occ - sd_occ, ymax = mean_occ + sd_occ),
                width = .2, position = position_dodge(.9)) + 
  theme_minimal()
# g5 <- ggplot(sc_hatch_occ_wkd, aes(x = week, y = mean_occ, fill = score)) +
#   geom_bar(stat="identity", alpha = 0.7 ) +
#   xlab('Week') +
#   ylab('Mean Occupancy') +
#   labs(title = 'Mean Occupancy Dec Paralarvae', 
#        subtitle = 'Spring catch',) +
#   theme_minimal() +
#   facet_wrap(~score)
# g6 = ggplot(data = sc_hatch_occ_annd,
#        mapping = aes(x = year,
#                      y = mean_occ)) +
#   geom_line(color = 'black') +
#   geom_point(aes(fill = factor(score)), shape = 21,  size = 4) +
#   scale_x_continuous(breaks = seq(2011,2020, by = 1)) +
#   xlab('Year') +
#   ylab('Ring days on shelf') + 
#   labs(title = 'Mean Occupancy Dec Paralarvae', subtitle = 'Spring catch',
#        fill = 'System state') +
#   geom_errorbar(aes(ymin = mean_occ - sd_occ, ymax = mean_occ + sd_occ),
#                 width = .2, position = position_dodge(.9)) + 
#   theme_minimal()
g7 <- ggplot(fc_hatch_occ_wk, aes(x = week, y = mean_occ, fill = score)) +
  geom_bar(stat="identity", alpha = 0.7 ) +
  xlab('Week') +
  ylab('Mean Occupancy') +
  labs(title = 'Mean Ring Occupancy June:July Hatching Season',
       subtitle = 'Fall catch') +
  theme_minimal() +
  facet_wrap(~score)
g8 = ggplot(data = fc_hatch_occ_ann,
       mapping = aes(x = year,
                     y = mean_occ)) +
  geom_line(color = 'black') +
  geom_point(aes(fill = factor(score)), shape = 21,  size = 4) +
  scale_x_continuous(breaks = seq(2011,2020, by = 1)) +
  xlab('Year') +
  ylab('Ring days on shelf') + 
  labs(title = 'Mean Ring Occupancy June:July Hatching Season', 
       subtitle = 'Fall catch', fill = 'System state') +
  geom_errorbar(aes(ymin = mean_occ - sd_occ, ymax = mean_occ + sd_occ),
                width = .2, position = position_dodge(.9)) + 
  theme_minimal()
g9 <- ggplot(hs_occ_wk, aes(x = week, y = mean_occ, fill = score)) +
  geom_bar(stat="identity", alpha = 0.7 ) +
  xlab('Week') +
  ylab('Mean Occupancy') +
  labs(title = 'Mean Ring Occupancy Jan:July Full Hatching Season',
       subtitle = 'Both cohorts') +
  theme_minimal() +
  facet_wrap(~score)
g10 = ggplot(data = hs_occ_ann,
            mapping = aes(x = year,
                          y = mean_occ)) +
  geom_line(color = 'black') +
  geom_point(aes(fill = factor(score)), shape = 21,  size = 4) +
  scale_x_continuous(breaks = seq(2011,2020, by = 1)) +
  xlab('Year') +
  ylab('Ring days on shelf') + 
  labs(title = 'Mean Ring Occupancy June:July Hatching Season', 
       subtitle = 'Fall catch', fill = 'System state') +
  geom_errorbar(aes(ymin = mean_occ - sd_occ, ymax = mean_occ + sd_occ),
                width = .2, position = position_dodge(.9)) + 
  theme_minimal()

Looking at Ring occupancy during fishing season

# Plot
year_plot

#occ_plot
tyear

myear

#g1
g2

Quick ANOVA to test differences:

# Quantify
aov2.1 <- aov(total.occ ~ score, data = occ)
Anova(aov2.1, type = 'III')

## Anova Table (Type III tests)
## 
## Response: total.occ
##              Sum Sq   Df F value    Pr(>F)    
## (Intercept)  5371.9    1 279.185 < 2.2e-16 ***
## score         551.4    2  14.329 6.951e-07 ***
## Residuals   25918.2 1347                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(aov2.1)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = total.occ ~ score, data = occ)
## 
## $score
##           diff       lwr        upr     p adj
## G-A -0.5493827 -1.225923  0.1271575 0.1375346
## P-A -1.6172840 -2.340536 -0.8940321 0.0000005
## P-G -1.0679012 -1.744441 -0.3913610 0.0006472

# Consider only rings that were on shelf for greater than 0 days
# aov2.1 <- aov(timeshelf ~ score, data = occ[which(occ$timeshelf>0),])
# Anova(aov2.1, type = 'III')
# TukeyHSD(aov2.1)

Results suggest that there is a significant difference in ring occupancy between Good and Poor years (p < 0.0001). Specifically, there is a significant difference in ring occupancy in Good vs Poor years.

Looking at Ring occupancy during prelarvae and juvenile time periods

# Plot 
g9/g10

g3/g4

#g5/g6

g7/g8

Digging deeper

It is clear from the above analyses that ‘Good’ years (years with higher than average catch) are not associated with a higher frequency of shorter lived WCRs. Furthermore, there was not a significant difference in number of rings per year (as expected re: Gangopadhyay_ea_20).

The next question to consider is whether there is a difference in the origin (location of birth) or demise (location of demise) of rings located near the shelf during the fishing season.

Silva_ea_20 found that rings born in Zone 2 and demised in zone 1 (Z2:Z1) had the highest survival probability. Additionally, WCRs last seen or demised in Zone 1 had a higher survival probability, with the exception of those also born in Zone 1 (Z1:Z1).

The following figures summarize analyses to determine how many rings were born and demised in each of the four zones described by Silva_ea_20 and assess how that plays into year differences. Specifically, does the location of origin (or demise) display a signal of ‘system state’?

Zone of Origin

The code below isolates the first week a ring was observed and uses the associated latitude and longitude during that week as the reference point to calculate its zone of origin. Geographic zones follow four 5◦ zones between 75 and 55◦ as defined by Silva and colleagues (2020).

# Calculate Origin
origin <- df %>% 
  dplyr::filter(weeks_old == 1) # create new dataframe which isolates 
                                # the first week of observation for each ring
z4 <-origin %>%
  dplyr::filter(lon > -60.000, lon < -55.000) # subset 'origin' df by zone
z3 <-origin %>%
  dplyr::filter(lon > -65.000, lon < -60.000)
z2 <-origin %>%
  dplyr::filter(lon > -70.000, lon < -65.000)
z1 <-origin %>%
  dplyr::filter(lon > -75.000, lon < -70.000)
# Number of rings from each zone should sum to 281 (total rings)..
test = sum(nrow(z4)+ nrow(z3)+ nrow(z2) + nrow(z1))

# Identify where rings that hit the shelf during the fishing season originated,
# keeping an eye out for those that were born in zones 2 or 1 in particular. 
origin_key <- rbind(cbind(as.character(z1$id), rep("zone1", nrow(z1))),
                    cbind(as.character(z2$id), rep("zone2", nrow(z2))),
                    cbind(as.character(z3$id), rep("zone3", nrow(z3))),
                    cbind(as.character(z4$id), rep("zone4", nrow(z4))))
colnames(origin_key) <- c("id", "origin")
origin_key <- as.data.frame(origin_key)
new_data <- dplyr::left_join(df,
                             origin_key,
                             by = "id")

Zone of Demise

The code below isolates the last week a ring was observed and uses the associated latitude and longitude during that week as the reference point to calculate its zone of demise.

z4 <-demise %>%
  dplyr::filter(lon <= -60.000, lon > -55.000)
z3 <-demise %>%
  dplyr::filter(lon <= -65.000, lon > -60.000)
z2 <-demise %>%
  dplyr::filter(lon <= -70.000, lon > -65.000)
z1 <-demise %>%
  dplyr::filter(lon <= -75.000, lon > -70.000)
# Number of rings from each zone should sum to 281 (total rings)
test = sum(nrow(z4)+ nrow(z3)+ nrow(z2) + nrow(z1))
test

## [1] 281

# Identify where rings that hit the shelf during the fishing season originated,
# keeping an eye out for those that were born in zones 2 or 1 in particular. 
demise_key <- rbind(cbind(as.character(z1$id), rep("zone1", nrow(z1))),
                    cbind(as.character(z2$id), rep("zone2", nrow(z2))),
                    cbind(as.character(z3$id), rep("zone3", nrow(z3))),
                    cbind(as.character(z4$id), rep("zone4", nrow(z4))))
colnames(demise_key) <- c("id", "demise")
demise_key <- as.data.frame(demise_key)
new_data <- dplyr::left_join(new_data,
                             demise_key,
                             by = "id")

Re-create subsets

Now that zones of origin and demise have been calculated and added to the original dataframe (‘df’) to create a new dataframe (‘new_data’), it’s time to re-create subsets for just rings near shelf (shelf_rings) and rings both near shelf and during the fishing season (shelf_rings_fs).

Specify rings that occur near the shelf and during the fishing season

shelf_rings <-  subset(new_data, lon <= -65.00)
shelf_rings_fs <- subset(shelf_rings, week %in% c(18:44))

Visualize trends in zones of origination and demise

For the next series of figures, the color of each bar represents the geographic (5◦) zone of origin for a given ring. Zones lay between 75 and 55◦ as defined by Silva and colleagues (2020). Zone 1 is represented by orange bars, Zone 2 blue bars, Zone 3 green bars and Zone 4 yellow bars.

A quick comparison of zone of origination for all rings vs only shelf rings across all years and weeks

# Color blind friendly palette :
cbp1 <- c("#999999", "#E69F00", "#56B4E9", "#009E73",
          "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

# Set up color scheme for following geographical zone plots 
zones <- sort(unique(shelf_rings$origin))
zone_cols <- cbp1[2:5]
names(zone_cols)  <- zones

# Generate plots 
p11 <- ggplot(new_data, aes(x = week, fill = origin), 
              binwidth = 40) +
  geom_histogram(alpha = 0.7) +
  scale_fill_manual('Legend', values = zone_cols) +
  xlab('Week') +
  ylab('Frequency')+
  labs(title = 'Origin location: All rings') +
  theme_minimal() +
  theme(legend.position = "none")
p1 <- ggplot(new_data, aes(x = week, fill = origin), 
             binwidth = 40) +
  geom_histogram(alpha = 0.7) +
  scale_fill_manual('Legend', values = zone_cols) +
  xlab('Week') +
  ylab('Frequency')+
  #labs(title = 'Origin location: All rings') +
  theme(legend.position = "none") +
  theme_minimal()+
  facet_wrap(~origin)
p2<- ggplot(shelf_rings, aes(x = week, fill = origin), binwidth = 40) +
  geom_histogram(alpha = 0.7) +
  scale_fill_manual('Legend', values = zone_cols) +
  xlab('Week') +
  ylab('Frequency')+
  labs(title = 'Origin location: Shelf rings') +
  theme_minimal()
p2.2<- ggplot(shelf_rings, aes(x = week, fill = origin), binwidth = 40) +
  geom_histogram(alpha = 0.7) +
  scale_fill_manual('Legend', values = zone_cols) +
  xlab('Week') +
  ylab('Frequency')+
  labs(title = 'Origin location: Shelf rings') +
  theme_minimal() +
  theme(legend.position = "none")
p3 <-ggplot(shelf_rings, aes(x = week, fill = origin), binwidth = 40) +
  geom_histogram(alpha = 0.7) +
  scale_fill_manual('Legend', values = zone_cols) +
  xlab('Week') +
  ylab('Frequency')+
  #labs(title = 'Origin location: Shelf rings') +
  theme_minimal()+
  facet_wrap(~origin)

# Plot the figures 
p11|p1 # All rings

p2.2 |p3 # Shelf rings

p11 | p2 # All rings vs Shelf Rings

Quick take-away from these figures is that a majority of ‘Shelf rings’ (rings near the shelf, or West of -65), originated in Zone 2.

Visualize: Zone of Origination

Next question - does this higher frequency of rings originating in Zone 2 map to a particular system state?

Code below calculates the Zone of origination for rings near the shelf which occur only during the fishing season across good, average and poor years.

dt <- shelf_rings_fs %>%
  dplyr::group_by(score, origin)%>%
  dplyr::tally()%>%
  dplyr::mutate(percent=n/sum(n))
#dt[c(12,13,14),1] <- c('G','G','G')
#dt <- dt[-c(12:14),]
or_plot <- ggplot(data = dt,aes(x= origin, y = n, fill = score))+ geom_bar(stat="identity")+
  geom_text(aes(label=paste0(sprintf("%1.1f", percent*100),"%")),
            position=position_stack(vjust=0.5), 
            colour="white", size = 2) +
  xlab('Zone of origination') +
  ylab('Number of Rings') + 
  #labs(title = 'Origin across') +
  theme_minimal()

a <-ggplot(shelf_rings_fs %>% dplyr::filter(score == 'G'), 
            aes(x = week, fill = origin)) +
  geom_histogram(alpha = 0.7) +
  scale_fill_manual('Legend', values = zone_cols) +
  ylim(0,8) +
  xlab('Week') +
  ylab('Frequency') + 
  labs(title = 'Origin: Good Years') +
  theme_minimal()+
  theme(legend.position = "none") 

aa <- ggplot(shelf_rings_fs %>% dplyr::filter(score == 'G'), 
            aes(x = week, fill = origin)) +
  geom_histogram(alpha = 0.7) +
  scale_fill_manual('Legend', values = zone_cols) +
  ylim(0,8) +
  xlab('Week') +
  ylab('Frequency') + 
  #labs(title = 'Origin: Good Years') +
  theme_minimal()+
  #theme(legend.position = "none") +
  facet_wrap(~origin)


b <- ggplot(shelf_rings_fs %>% dplyr::filter(score == 'A'), 
            aes(x = week, fill = origin)) +
  geom_histogram(alpha = 0.7) +
  scale_fill_manual('Legend', values = zone_cols) +
  ylim(0,8) +
  xlab('Week') +
  ylab('Frequency') + 
  labs(title = 'Origin: Average Years') +
  theme_minimal()+
  theme(legend.position = "none") 
bb <- ggplot(shelf_rings_fs %>% dplyr::filter(score == 'A'), 
             aes(x = week, fill = origin)) +
  geom_histogram(alpha = 0.7) +
  scale_fill_manual('Legend', values = zone_cols) +
  ylim(0,8) +
  xlab('Week') +
  ylab('Frequency') + 
  #labs(title = 'Origin: Average Years') +
  theme_minimal()+
  #theme(legend.position = "none") +
  facet_wrap(~ origin)

c <- ggplot(shelf_rings_fs %>% dplyr::filter(score == 'P'), 
            aes(x = week, fill = origin)) +
  geom_bar(alpha = 0.7) +
  scale_fill_manual('Legend', values = zone_cols) +
  ylim(0,8) +
  xlab('Week') +
  ylab('Frequency') + 
  labs(title = 'Origin: Poor Years') +
  theme_minimal()  +
  theme(legend.position = "none")

  cc <- ggplot(shelf_rings_fs %>% dplyr::filter(score == 'P'), 
               aes(x = week, fill = origin)) +
  geom_histogram(alpha = 0.7) +
  scale_fill_manual('Legend', values = zone_cols) +
  ylim(0,8) +
  xlab('Week') +
  ylab('Frequency') + 
  #labs(title = 'Origin: Poor Years') +
  theme_minimal() +
  facet_wrap(~ origin)

# o_num <-  shelf_rings_fs %>%
#   group_by(year, score, origin) %>%
#   summarise(number = length(unique(id))) 
# o_num$prop <- o_num$number/ring_num$number
# 
# g3<-ggplot(o_num %>% dplyr::filter(score == 'G'),
#        aes(x=year, y=number, fill = origin)) +
#   geom_bar(stat="identity") +
#   scale_fill_manual('Legend', values = zone_cols) +
#   ylim(0,15) +
#   xlab('Year') +
#   ylab('Number of rings') + 
#   labs(title = 'Origin: Good Years') +
#   theme_minimal() +
#    theme(legend.position = "none") +
#   facet_wrap(~ origin)
# 
# g4<-ggplot(o_num %>% dplyr::filter(score == 'P'),
#        aes(x=year, y=number, fill = origin)) +
#   geom_bar(stat="identity", position = 'fill') +
#   scale_fill_manual('Legend', values = zone_cols) +
#   ylim(0,15) +
#   xlab('Year') +
#   ylab('Number of rings') + 
#   labs(title = 'Origin: Poor Years') +
#   theme_minimal() +
#   facet_wrap(~ origin)
# 
# g3|g4

# Plot figures   



or_plot

a|aa

b|bb

c|cc

a | b | c

The above figures suggest that in Poor years - there is a lower frequency of rings born in Zone 2 and and higher frequency of rings born in Zone 1 than in Good years. Rings born in these zones were found to have a higher probability of survival than those born in zones 3 or 4 (Silva_ea_20). This is consistent with previous findings about age.

Quantify differences

shelf_rings_fs$state <- NA
# Add system state
shelf_rings_fs$state <- NA
locs1 <- which(shelf_rings_fs$year == 2011)
shelf_rings_fs[locs1,15] <- '1'
locs2 <- which(shelf_rings_fs$year == 2012)
shelf_rings_fs[locs2,15] <- '1'
locs3 <- which(shelf_rings_fs$year == 2013)
shelf_rings_fs[locs3,15] <- '0'
locs4 <- which(shelf_rings_fs$year == 2014)
shelf_rings_fs[locs4,15] <- '1'
locs5 <- which(shelf_rings_fs$year == 2015)
shelf_rings_fs[locs5,15] <- '0'
locs6 <- which(shelf_rings_fs$year == 2016)
shelf_rings_fs[locs6,15] <- '0'
locs7 <- which(shelf_rings_fs$year == 2017)
shelf_rings_fs[locs7,15] <- '2'
locs8 <- which(shelf_rings_fs$year == 2018)
shelf_rings_fs[locs8,15] <- '2'
locs9 <- which(shelf_rings_fs$year == 2019)
shelf_rings_fs[locs9,15] <- '2'
locs10 <- which(shelf_rings_fs$year == 2020)
shelf_rings_fs[locs10,15] <- '2'

aov3.1 <- aov(year ~ origin, data = shelf_rings_fs)
Anova(aov3.1, type = 'III')

## Anova Table (Type III tests)
## 
## Response: year
##                Sum Sq  Df    F value    Pr(>F)    
## (Intercept) 580104506   1 8.1413e+07 < 2.2e-16 ***
## origin            172   3 8.0316e+00 2.893e-05 ***
## Residuals        4838 679                         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(aov3.1)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = year ~ origin, data = shelf_rings_fs)
## 
## $origin
##                     diff        lwr        upr     p adj
## zone2-zone1  0.933106649  0.2438838  1.6223295 0.0029198
## zone3-zone1  0.927510634  0.1698082  1.6852131 0.0091432
## zone4-zone1 -1.276775856 -2.9554458  0.4018941 0.2048002
## zone3-zone2 -0.005596015 -0.6286085  0.6174165 0.9999955
## zone4-zone2 -2.209882504 -3.8322155 -0.5875495 0.0027037
## zone4-zone3 -2.204286489 -3.8568751 -0.5516979 0.0035076

aov3.2 <- aov(week ~ origin, data = shelf_rings_fs)
Anova(aov3.2, type = 'III')

## Anova Table (Type III tests)
## 
## Response: week
##             Sum Sq  Df   F value Pr(>F)    
## (Intercept) 130084   1 2252.5618 <2e-16 ***
## origin         189   3    1.0913 0.3521    
## Residuals    39212 679                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

aov3.3 <- aov(state ~ origin, data = shelf_rings_fs)
Anova(aov3.3, type = 'III')

## Anova Table (Type III tests)
## 
## Response: state
##             Sum Sq  Df  F value    Pr(>F)    
## (Intercept)  75.64   1 119.2946 < 2.2e-16 ***
## origin       10.91   3   5.7362 0.0007045 ***
## Residuals   430.51 679                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

table(shelf_rings_fs$origin, shelf_rings_fs$score)

##        
##           A   G   P
##   zone1  42  31  70
##   zone2 111 112 104
##   zone3  74  58  62
##   zone4  12   0   7

#chisq.test(table(shelf_rings_fs$origin, shelf_rings_fs$score))$expected

chisq.test(shelf_rings_fs$origin, shelf_rings_fs$score, correct=FALSE)

## 
##  Pearson's Chi-squared test
## 
## data:  shelf_rings_fs$origin and shelf_rings_fs$score
## X-squared = 26.58, df = 6, p-value = 0.0001735

Next, run a Post-hoc Tukey Honest Significant Differences test to determine between which zone of origination those differences lie. The test compares all possible pairs of means.

TukeyHSD(aov3.3)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = state ~ origin, data = shelf_rings_fs)
## 
## $origin
##                    diff         lwr        upr     p adj
## zone2-zone1  0.29719210  0.09159953 0.50278468 0.0012189
## zone3-zone1  0.25210872  0.02608893 0.47812850 0.0217822
## zone4-zone1 -0.09569378 -0.59643469 0.40504713 0.9608346
## zone3-zone2 -0.04508339 -0.23092566 0.14075888 0.9241229
## zone4-zone2 -0.39288588 -0.87682169 0.09104992 0.1570856
## zone4-zone3 -0.34780250 -0.84076345 0.14515846 0.2661213

The results show there is a significant difference in the zone of origination of rings near the shelf between both years and system state. Specifically, there is a difference in rings born in Zone 2 vs Zone 1 that depends on whether the year is classified as a Good, Average or Poor year. The same applies for Rings born in Zones 3 vs Zones 1.

Visualize: Zone of Demise

Zone of demise for shelf rings occuring during the fishing season in good, average and poor years.

The following figures have the same color scheme as before with Zone 1 represented by orange bars, Zone 2 blue, Zone 3 green and Zone 4 yellow bars.

dt <- shelf_rings_fs %>%
  dplyr::group_by(score, demise)%>%
  dplyr::tally()%>%
  dplyr::mutate(percent=n/sum(n))
# dt[c(12,13,14),1] <- c('G','G','G')
# dt <- dt[-c(8:9),]
dm_plot <- ggplot(data = dt, aes(x= demise, y = n, fill = score))+ geom_bar(stat="identity")+
  geom_text(aes(label=paste0(sprintf("%1.1f", percent*100),"%")),
            position=position_stack(vjust=0.5), 
            colour="white", size = 2) +
  xlab('Zone of demise') +
  ylab('Number of Rings') + 
  #labs(title = 'Origin across') +
  theme_minimal()

a <- ggplot(shelf_rings_fs %>% dplyr::filter(score == 'G'), 
       aes(x = week, fill = demise)) +
  geom_histogram(alpha = 0.7) +
  scale_fill_manual('Legend', values = zone_cols) +
  ylim(0,15) +
  xlab('Week') +
  ylab('Frequency') + 
  labs(title = 'Demise: Good Years') +
  theme_minimal() +
  theme(legend.position = "none") 
  
  

b <- ggplot(shelf_rings_fs %>% dplyr::filter(score == 'A'), 
       aes(x = week, fill = demise)) +
  geom_histogram(alpha = 0.7) +
  scale_fill_manual('Legend', values = zone_cols) +
  ylim(0,15) +
  xlab('Week') +
  ylab('Frequency') + 
  labs(title = 'Average Years') +
  theme_minimal() +
  theme(legend.position = "none") 



c <- ggplot(shelf_rings_fs %>% dplyr::filter(score == 'P'), 
       aes(x = week, fill = demise)) +
  geom_histogram(alpha = 0.7) +
  scale_fill_manual('Legend', values = zone_cols) +
  ylim(0,15) +
  xlab('Week') +
  ylab('Frequency') + 
  labs(title = 'Poor Years') +
  theme_minimal()

# Plot figures
dm_plot

# a | b | c 
a|c

Qualitatively it appears that there are more rings near the shelf that demise in Zones 2 in Poor years as compared to Good years.

Quantify differences

aov4.1 <- aov(year ~ demise, data = shelf_rings_fs)
Anova(aov4.1, type = 'III')

## Anova Table (Type III tests)
## 
## Response: year
##                 Sum Sq  Df    F value  Pr(>F)    
## (Intercept) 2042204516   1 2.8012e+08 < 2e-16 ***
## demise              52   2 3.5888e+00 0.02816 *  
## Residuals         4958 680                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(aov4.1)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = year ~ demise, data = shelf_rings_fs)
## 
## $demise
##                   diff       lwr         upr     p adj
## zone2-zone1 -0.6098937 -1.161866 -0.05792161 0.0260938
## zone3-zone1 -1.9562624 -8.304698  4.39217328 0.7494694
## zone3-zone2 -1.3463687 -7.706194  5.01345660 0.8726492

aov4.2 <- aov(days_old ~ demise, data = shelf_rings_fs)
Anova(aov4.2, type = 'III')

## Anova Table (Type III tests)
## 
## Response: days_old
##              Sum Sq  Df F value    Pr(>F)    
## (Intercept) 3401387   1 857.018 < 2.2e-16 ***
## demise        92309   2  11.629 1.081e-05 ***
## Residuals   2698828 680                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(aov4.2)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = days_old ~ demise, data = shelf_rings_fs)
## 
## $demise
##                  diff        lwr       upr     p adj
## zone2-zone1 -25.76333  -38.64199 -12.88467 0.0000094
## zone3-zone1 -75.23260 -223.35483  72.88963 0.4577352
## zone3-zone2 -49.46927 -197.85725  98.91870 0.7135815

aov4.3 <- aov(state ~ demise, data = shelf_rings_fs)
Anova(aov4.3, type = 'III')

## Anova Table (Type III tests)
## 
## Response: state
##             Sum Sq  Df  F value Pr(>F)    
## (Intercept) 473.45   1 733.7322 <2e-16 ***
## demise        2.64   2   2.0466   0.13    
## Residuals   438.78 680                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

table(shelf_rings_fs$demise, shelf_rings_fs$score)

##        
##           A   G   P
##   zone1 172 158 173
##   zone2  67  43  69
##   zone3   0   0   1

ctest <- chisq.test(table(shelf_rings_fs$demise, shelf_rings_fs$score))$expected
fisher.test(table(shelf_rings_fs$demise, shelf_rings_fs$score))

## 
##  Fisher's Exact Test for Count Data
## 
## data:  table(shelf_rings_fs$demise, shelf_rings_fs$score)
## p-value = 0.1828
## alternative hypothesis: two.sided

The results show that there is a significant difference in the Zone of demise for rings near the shelf between years, but not system state. Specifically, there is a difference between rings which demised in Zone 2 versus Zone 1 (p = 0.0260938). Additionally, the ring age differed significantly across weeks for rings which demised in Zone 2 versus those that demised in Zone 1 (p = 0.0000094).

Explore Age of Demise across years and system states

# Add zone of demise to demise_fs data frame 
# (which has age of demise for each ring during fishing season)
demise_fs <- dplyr::left_join(demise_fs,
                             demise_key,
                             by = "id")

write.csv(demise_fs, 'demise_fs.csv', row.names = FALSE)


dz21 <- subset(demise_fs, demise == c('zone1', 'zone2'))


ann_demise_age <- demise_fs %>%
  group_by(year, score) %>%
  summarise(n = sum(days_old),
            mean_aod = mean(days_old),
            sd_aod = sd(days_old),
            max_aod = max(days_old),
            min_aod = min(days_old))
wk_demise_age <- demise_fs %>%
  group_by(year, week, score) %>%
  summarise(n = sum(days_old),
            mean_aod = mean(days_old),
            sd_aod = sd(days_old),
            max_aod = max(days_old),
            min_aod = min(days_old))

g1 = ggplot(data = ann_demise_age,
       mapping = aes(x = year,
                     y = mean_aod)) +
  geom_line(color = 'black') +
  geom_point(aes(fill = factor(score)), shape = 21,  size = 4) +
  scale_x_continuous(breaks = seq(2011,2020, by = 1)) +
  xlab('Year') +
  ylab('Age of Demise (days)') + 
  labs(title = 'Mean Age of Demise', fill = 'System state') +
  geom_errorbar(aes(ymin = mean_aod - sd_aod, ymax = mean_aod + sd_aod),
                width = .2, position = position_dodge(.9)) + 
  theme_minimal()

g2 <- ggplot(data = wk_demise_age,
       mapping = aes(x = week,
                     y = mean_aod)) +
  geom_line(color = 'black') +
  geom_point(aes(fill = factor(score)), shape = 21,  size = 4) +
  scale_x_continuous(breaks = seq(18,44, by = 2)) +
  xlab('Week') +
  ylab('Age of demise (days)') + 
  labs(title = 'Mean Age of Demise', fill = 'System state') +
  geom_errorbar(aes(ymin = mean_aod - sd_aod, ymax = mean_aod + sd_aod),
                width = .2, position = position_dodge(.9)) + 
  theme_minimal()+
  facet_wrap(~score)



g1

g2

aov5.1 <- aov(days_old ~ demise, data = dz21)
Anova(aov5.1, type = 'III')

## Anova Table (Type III tests)
## 
## Response: days_old
##             Sum Sq Df F value    Pr(>F)    
## (Intercept) 169344  1 39.8306 2.383e-07 ***
## demise        4632  1  1.0894    0.3034    
## Residuals   157310 37                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

aov5.2 <- aov(days_old ~ demise, data = demise_fs)
Anova(aov5.2, type = 'III')

## Anova Table (Type III tests)
## 
## Response: days_old
##             Sum Sq  Df  F value    Pr(>F)    
## (Intercept) 315467   1 119.6513 < 2.2e-16 ***
## demise       69555   3   8.7937  1.94e-05 ***
## Residuals   429758 163                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(aov5.2)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = days_old ~ demise, data = demise_fs)
## 
## $demise
##                   diff       lwr        upr     p adj
## zone2-zone1 -24.795699 -56.35546   6.764063 0.1779239
## zone3-zone1 -43.775758 -71.08802 -16.463497 0.0002972
## zone4-zone1 -53.102564 -82.74155 -23.463578 0.0000401
## zone3-zone2 -18.980059 -48.91406  10.953945 0.3559825
## zone4-zone2 -28.306865 -60.37795   3.764219 0.1043274
## zone4-zone3  -9.326807 -37.22834  18.574726 0.8215119

aov5.3 <- aov(days_old ~ score, data = dz21)
Anova(aov5.3, type = 'III')

## Anova Table (Type III tests)
## 
## Response: days_old
##             Sum Sq Df F value   Pr(>F)   
## (Intercept)  54444  1 12.1615 0.001304 **
## score          777  2  0.0868 0.917043   
## Residuals   161164 36                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The age of demise for rings near the shelf that demise in Zone 1 is significantly different than those that demise in Zones 3 and 4.

Determine the number of decaying rings in the region in good vs poor years

Silva_ea_20 found that WCRs are most productive during their decaying stage. If WCRs have the potential to be a food aggregation resource or serve as a hotspot for productivity, I hypothesized that Good years may be associated with a greater number of ‘decaying’ rings.

# Determine the number of decaying rings in the region in good v poor years
num_decay <- demise_fs %>%
  group_by(year,week, score) %>%
  summarise(n = length(unique(id)))

wk_decay <- num_decay %>%
  group_by(year,week, score) %>%
  summarise(dmean = mean(n), 
            dsd = sd(n))

decay_tally <- demise_fs %>%
  dplyr::group_by(week, score)%>%
  dplyr::tally()%>%
  dplyr::mutate(percent=n/sum(n))
# Generate a series of figures
 d1 <- ggplot(data = decay_tally, aes(x = week, y = n, fill = score)) +
  geom_bar(stat="identity") +
  geom_text(aes(label=paste0(sprintf("%1.1f", percent*100),"%")),
            position=position_stack(vjust=0.5), 
            colour="white", size = 2) +
  xlab('Week') +
  ylab('Number of Rings') + 
  labs(title = 'Decaying Rings') +
  theme_minimal()
 
 d2 <- ggplot(data = wk_decay,
       mapping = aes(x = week,
                     y = dmean)) +
  geom_line(color = 'black') +
  geom_point(aes(fill = factor(score)), shape = 21,  size = 4) +
  scale_x_continuous(breaks = seq(18,44, by = 2)) +
  xlab('Week') +
  ylab('Number of rings') + 
  labs(title = 'Decaying Rings by Week', fill = 'System state') +
  #geom_errorbar(aes(ymin = mean_aod - sd_aod, ymax = mean_aod + sd_aod),
  #              width = .2, position = position_dodge(.9)) + 
  theme_minimal()+
  facet_wrap(~year)

 d3 <- ggplot(data = wk_decay,
       mapping = aes(x = week,
                     y = dmean)) +
  geom_line(color = 'black') +
  geom_point(aes(fill = factor(score)), shape = 21,  size = 4) +
  scale_x_continuous(breaks = seq(18,44, by = 2)) +
  xlab('Week') +
  ylab('Number of rings') + 
  labs(title = 'Decaying Rings by Week', fill = 'System state') +
  #geom_errorbar(aes(ymin = mean_aod - sd_aod, ymax = mean_aod + sd_aod),
  #              width = .2, position = position_dodge(.9)) + 
  theme_minimal()+
  facet_wrap(~score)
# Plot figures
d1

#d2
d3

# Compute stats
aov6.1 <- aov(n ~ score * week, data = num_decay)
Anova(aov6.1, type = 'III')

## Anova Table (Type III tests)
## 
## Response: n
##             Sum Sq  Df F value  Pr(>F)  
## (Intercept)  3.005   1  5.6209 0.01943 *
## score        0.528   2  0.4933 0.61189  
## week         0.001   1  0.0010 0.97501  
## score:week   0.625   2  0.5845 0.55906  
## Residuals   60.955 114                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

aov6.2 <- aov(n ~ year * week, data = num_decay)
Anova(aov6.2, type = 'III')

## Anova Table (Type III tests)
## 
## Response: n
##             Sum Sq  Df F value Pr(>F)
## (Intercept)  0.172   1  0.3246 0.5700
## year         0.177   1  0.3337 0.5646
## week         0.213   1  0.4018 0.5274
## year:week    0.214   1  0.4040 0.5263
## Residuals   61.385 116

There is no significant difference in the number of decaying rings in Good vs Bad years. **NOTE I calculated number of ‘decaying’ rings as the sum of the rings that were in their last week of existence (last week of their tracking timeline). I am not sure if this is the best way to do this - is there a known trajectory/decay rate taht is a better way to quantify this?

Visualize: Rings born in Zone 2, Demised in Zone 1

Silva and colleagues (2020) identified rings that originated in zone 2 and demised in zone 1 as having the highest survival probability. The following explores the number and proportion of the rings that are near the shelf during the fishing season are of this category of high survival probability (Z2Z1).

# Calculate the number birthed in z2 and demised in z1 per year
z2z1 <- subset(shelf_rings_fs, 
               origin == 'zone2' & demise == 'zone1')
z2z1_num <-  z2z1 %>%
  group_by(year, score) %>%
  summarise(number = length(unique(id))) 
z2z1_num$prop <- z2z1_num$number/ring_num$number

p1 = ggplot(z2z1_num, aes(x=year, y=number, fill = score)) +
  geom_bar(stat="identity" ) +
  xlab('Year') +
  ylab('Number of Rings') +
  labs(title = 'Rings born in Zone2 Demised in Zone1') +
  geom_text(aes(label=year), vjust=1.6, color="white", size=3.5)+
  theme_minimal() +
  theme(axis.text.x=element_blank()) 


p2 = ggplot(z2z1_num, aes(x=year, y=prop, fill = score)) +
  geom_bar(stat="identity" ) +
  xlab('Year') +
  ylab('Proportion of Rings') +
  #labs(title = 'Proportion of rings born in Zone2 Demised in Zone1') +
  geom_text(aes(label=year), vjust=1.6, color="white", size=3.5)+
  theme_minimal()+
  theme(axis.text.x=element_blank()) 

# Plot
p1/p2

Quantify

table(z2z1_num$number, z2z1_num$score)

##    
##     A G P
##   1 0 0 2
##   2 2 1 1
##   3 1 3 0

aov5.1 <- aov(number ~ score, data = z2z1_num)
Anova(aov5.1, type = 'III')

## Anova Table (Type III tests)
## 
## Response: number
##              Sum Sq Df F value    Pr(>F)    
## (Intercept) 16.3333  1  54.880 0.0001484 ***
## score        3.5167  2   5.908 0.0314049 *  
## Residuals    2.0833  7                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(aov5.1)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = number ~ score, data = z2z1_num)
## 
## $score
##           diff       lwr        upr     p adj
## G-A  0.4166667 -0.810441  1.6437743 0.5999631
## P-A -1.0000000 -2.311833  0.3118333 0.1307893
## P-G -1.4166667 -2.643774 -0.1895590 0.0271497

aov5.2 <- aov(prop ~ score, data = z2z1_num)
Anova(aov5.2, type = 'III')

## Anova Table (Type III tests)
## 
## Response: prop
##                Sum Sq Df F value    Pr(>F)    
## (Intercept) 0.0193094  1 34.4474 0.0006185 ***
## score       0.0020749  2  1.8508 0.2263485    
## Residuals   0.0039238  7                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

In the regions west of -65, there are a greater number of rings that were born in Zone2 and demised in Zone1 in Good years, than in Poor years. This suggests a that near the shelf there is a higher proportion of rings with high survival probability in Good years.

Determine distance between fishing points and ring corridor

dlrlog2 <- dlrlog1 %>% dplyr::filter(is.na(LONGDD) == FALSE, is.na(LATDD) == FALSE)
# make a spatial points object with the coordinates from the dataset
latlon <- SpatialPointsDataFrame(cbind(dlrlog2$LONGDD,dlrlog2$LATDD), 
                        proj4string = CRS("+proj=longlat +datum=NAD83"),
                        data = dlrlog2,
                        match.ID = "VTRIPID") 

UTMcoord <- spTransform(latlon,  CRS("+proj=utm +zone=19 +datum=NAD83 +units=m")) %>% as.data.frame() # transform the gps coordinates to UTM coordinates
# dlrlog1 <- cbind(dlrlog1,as.data.frame(UTMcoord)) %>% # append the UTM coordinates to the data
#   rename(UTMx = coords.x1, UTMy = coords.x2) # rename them to something more intuitive

dlrlog1 <- left_join(dlrlog1,UTMcoord) %>% # append the UTM coordinates to the data
  rename(UTMx = coords.x1, UTMy = coords.x2) # rename them to something more intuitive



dlrlog_shelfdist <- dist2isobath(pt.baty415, x = dlrlog2[,"LONGDD"], y = dlrlog2[,"LATDD"], isobath = -200)
dlrlog2$shelfdist <- dlrlog_shelfdist$distance
dlrlog1 <- left_join(dlrlog1, dlrlog2)
rm(latlon,UTMcoord)

Extracting data from indivdual Warm Core Rings

Create polygon representing WCR in space
Extract chl from ring (and buffering location)
Determine angle from fishing point to ring

library(tidyverse)
library(viridis)
library(mapdata)
library(sf)
library(marmap) # get bathy data
library(ggspatial)
## get coastline from mapdata 
reg = map_data("world2Hires")
reg = subset(reg, region %in% c('Canada', 'USA', 'Mexico'))
reg$long = (360 - reg$long)*-1 
## bring in bathymetry data from marmap package
pt.baty415 <- getNOAA.bathy(lon1 = -66, lon2 = -76, lat1 = 35,
                            lat2 = 41.5, resolution = 1)

### Create polygon 
E2015.27 <- subset(df, id == 'E2015' & week == 27)
E2015.28 <- subset(df, id == 'E2015' & week == 28)
E2015.30 <- subset(df, id == 'E2015' & week == 30)
U2019.31 <- subset(df, id == 'U2019' & week == 31)
U2019.35 <- subset(df, id == 'U2019' & week == 35)
U2020.43 <- subset(df, id == 'U2020' & week == 43)
c14 = circle.polygon(E2015.27$lon,E2015.27$lat, E2015.27$radius,
                     poly.type = "gc.earth", units = 'km')

c15 = circle.polygon(E2015.28$lon,E2015.28$lat, E2015.28$radius,
                     poly.type = "gc.earth", units = 'km')
c16 = circle.polygon(E2015.30[1,5],E2015.30[1,6], E2015.30[1,8],
                     poly.type = "gc.earth", units = 'km')
c3 = circle.polygon(U2019.31$lon,U2019.31$lat, U2019.31$radius,
                    poly.type = "gc.earth", units = 'km')
c7 = circle.polygon(U2019.35$lon,U2019.35$lat, U2019.35$radius,
                    poly.type = "gc.earth", units = 'km')
c4 = circle.polygon(U2020.43$lon,U2020.43$lat, U2020.43$radius,
                    poly.type = "gc.earth", units = 'km')

circ1 = as.data.frame(c4)
circ2 = as.data.frame(c15)
circ3 = as.data.frame(c16)
cir.90 = circ1[nrow(circ1)/4,c(1,2)] # 90 degrees

## make a plot to test angles for reference point 
ggplot() +
  geom_polygon(data = reg, aes(x=long, y = lat, group = group), 
               color = "gray20", fill = "wheat3") +
  coord_sf(xlim = c(-80,-66), ylim = c(33,48)) + 
  geom_contour(data = pt.baty415, aes(x = x, y = y, z = z), 
               breaks=c(-200), colour="gray20", size=0.2) +
  geom_polygon(data = circ1, aes(x=lon, y = lat),
               color = 'goldenrod', fill = NA, alpha = 1) +
  geom_point(data = circ1,
             mapping = aes(x = circ1[1,1], y = circ1[1,2], color = 'red'),
             alpha = 0.5) +
  geom_point(data = circ1,
             mapping = aes(x = circ1[nrow(circ1)/4,1],
                           y = circ1[nrow(circ1)/4,2], color = 'blue'),
             alpha = 0.5) +
  geom_polygon(data = circ2, aes(x=lon, y = lat),
               color = 'goldenrod', fill = NA, alpha = 1) +
  geom_point(data = circ2,
             mapping = aes(x = circ2[1,1], 
                           y = circ2[1,2], color = 'red'),
             alpha = 0.5) +
  geom_point(data = circ2,
             mapping = aes(x = circ2[nrow(circ2)/4,1], 
                           y = circ2[nrow(circ2)/4,2], color = 'blue'),
             alpha = 0.5) +
  geom_polygon(data = circ3, aes(x=lon, y = lat),
               color = 'goldenrod', fill = NA, alpha = 1) +
  geom_point(data = circ3,
             mapping = aes(x = circ3[1,1], y = circ3[1,2], color = 'red'),
             alpha = 0.5) +
  geom_point(data = circ3,
             mapping = aes(x = circ3[nrow(circ3)/4,1],
                           y = circ3[nrow(circ3)/4,2], color = 'blue'),
             alpha = 0.5) +
  # annotation_north_arrow(location = 'tl', which_north = "true",
  #                        pad_x = unit(0.75, "in"), pad_y = unit(0.5, "in"),
  #                        style = north_arrow_minimal) +
  theme_minimal()+
  theme(legend.position = "none")

# *** do I want to use years or ring id? ***
for(i in 1:length(unq_shelf_rings_fs)){
  tmp <- subset(shelf_rings_fs, id == unq_shelf_rings_fs[i])
  weeks <- sort(unique(tmp$week))
  for(k in 1:length(bricks)){
    btmp <- brick(bricks[k])
    btmp@data@names <- as.numeric(c(1:52))
    if(substr(bricks[k], start =  82, stop =  83) == substr(years[i], start =  3, stop =  4)){
      for(j in 1:length(weeks)){
        tmp2 <- subset(tmp, week == weeks[j])
        #coords <- data.frame(x = tmp2$lon, y = tmp2$lat)
        ring = circle.polygon(tmp2$lon, tmp2$lat, tmp2$radius,
                              poly.type = "gc.earth", units = 'km')
        coords <- data.frame(x = ring$lon, y = ring$lat)
        sp_ring <- SpatialPoints(coords)
        rtmp = raster::subset(btmp, which(btmp@data@names == weeks[j]))
        image(rtmp, main = paste0('Week',weeks[j])) # quick visual check
        ring_chl <- extract(rtmp, sp_ring, method='bilinear')
        locs <- which(cbwk$YEAR == years[i] & cbwk$WEEK == weeks[j]) # Change to correct DF
        cbwk$ring_chl[locs] <- ring_chl
        
      }
    }   
  }
  
  
  
 # create a circle
U2020.40 <- subset(df, id == 'U2020' & week == 40)
c1 = circle.polygon(U2020.40$lon,U2020.40$lat, U2020.40$radius,
                   poly.type = "gc.earth", units = 'km')

# Plot on map 
op <- par(mar = c(3, 5, 5, 5) + 0.1, oma = c(1, 1, 1, 1))
maps::map("mapdata::worldHires", fill = TRUE, col = "wheat3",
          xlim = lon.range, ylim = lat.range)
points(E2015.13$lon, E2015.13$lat, pch = 19, col = "goldenrod")
polygon(c1, border = "goldenrod", lwd = 3)
points(fish[,1], fish[,2], pch = 19, col = 'black', cex = 2)
lat.lon.axes(n = 3)
box(lwd = 2)
 

# plot the track of ring
plot(makeLine(tmp), type='l',  main = paste0('WCR_ID :',unq_shelf_rings_fs[i]))
points(tmp)
}

# Plot polygon of ring
for(j in 1:length(years)){
  tmp <- subset(shelf_rings_fs, year == years[j])
  weeks = sort(unique(tmp$week))
  for(i in 1:unq_shelf_rings_fs){
    tmp2 <- subset(shelf_rings_fs, id == unq_shelf_rings_fs[i])
    for(k in 1:length(weeks)){
      tmp3 <- subset(tmp2, week == weeks[k])
      c = circle.polygon(tmp3$lon, tmp3$lat, tmp3$radius,
                         poly.type = "gc.earth", units = 'km')
      maps::map("mapdata::worldHires", fill = TRUE, col = "wheat3",
                xlim = lon.range, ylim = lat.range)
      points(tmp3$lon, tmp3$lat, pch = 19, col = "red")
      polygon(c1, border = "red", lwd = 3)
      lat.lon.axes(n = 3)
      box(lwd = 2)
      
    }
  }
}

##### PLOT circle #####
lon.range = c(-77.99, -62.00)
lat.range = c(34.01, 46.00)


U2020.42 <- subset(df, id == 'U2020' & week == 42)

c3 = circle.polygon(U2020.42$lon,U2020.42$lat, U2020.42$radius,
                   poly.type = "gc.earth", units = 'km')

t3 = c3[nrow(c3)/4,c(1,2)] # 90 degrees

fish = data.frame(lat = c(42.85712167,  42.80302, 42.85641333, 42.854875,   42.78498833, 
                          42.80439, 42.85678333, 42.804105, 42.86433333, 42.83158,                                    42.776335, 42.83301167, 42.852615), 
                  lon = c(-70.652395, -70.64875333, -70.644765, -70.66094667,
                          -70.64431667,-70.65403333,-70.65253,-70.65158667,-70.63724667,
                          -70.66055333,-70.63993,-70.65673333,-70.65805667))

fishy <- fish[1,c(2,1)]
# Plot polygon and fishing points
maps::map("mapdata::worldHires", fill = TRUE, col = "wheat3",
          xlim = lon.range, ylim = lat.range)
par(new=TRUE)
polygon(c3, border = "goldenrod", lwd = 3)
points(fishy$lon, fishy$lat, pch = 17, col = "black", cex = 1.3)
points(c3[1,1], c3[1,2], pch = 19, col = "red") # 270 degrees
points(t3[1], t3[2], pch = 19, col = "blue")
lat.lon.axes(n = 3)
box(lwd = 2)
par(new=TRUE)
contour(grid$lon, grid$lat, isobath,
        levels = c(-200),
        lwd = c(2),
        lty = c(3),
        drawlabels = F, add = TRUE, col = 'darkgray')

contour(bathyLon,bathyLat,bathyZ,
        levels = c(-50, -100, -150, -200, -250),
        lwd = c(1, 1, 2, 2, 3),
        lty = c(3, 1, 3, 1, 3),
        drawlabels = F, add = TRUE, col = 'darkgray')


geom_contour(data = pt.baty415, aes(x = x, y = y, z = z), breaks=c(-200), colour="gray20", size=0.2) 

levelplot(isobath ~ lon * lat, data = grid,  
          pretty = T, col.regions = 'black') +
  xyplot(lat ~ lon, state,  type = 'l', 
         col = '#575757') +
polygon(c3, border = "goldenrod", lwd = 3)+
#points(fishy$lon, fishy$lat, pch = 17, col = "black", cex = 1.3)
points(c3[1,1], c3[1,2], pch = 19, col = "red") + 
points(t3[1], t3[2], pch = 19, col = "blue")



# Generate 
# (i) distance from fishing point to ring 
library(useful)


# Calculate angle between fishing point and ring to determine orientation

To = c3[nrow(c3)/4,c(1,2)]
From = as.matrix(fishy)
df2 = To - From
results = cart2pol(df2[,1], df2[,2], degrees = T)

Calculating Warm Core Ring Metrics

Sarah Salois

10/14/2021

Setting up your environment

Modify data

Calculate first metric : Ages of rings

Focus on rings West of -65 the fishing season

Visualize: age data

Quantify

Ring Occupancy

The following code uses these assumptions about the timing of maturation to explore WCR occupancy during paralarvae and juvenile phases - where rings may serve as transport or food resource aggregation mechanisms.

Looking at Ring occupancy during fishing season

Looking at Ring occupancy during prelarvae and juvenile time periods

Digging deeper

Zone of Origin

Zone of Demise

Re-create subsets

Visualize trends in zones of origination and demise

Visualize: Zone of Origination

Quantify differences

Visualize: Zone of Demise

Quantify differences

Explore Age of Demise across years and system states

Determine the number of decaying rings in the region in good vs poor years

Visualize: Rings born in Zone 2, Demised in Zone 1

Quantify

Determine distance between fishing points and ring corridor

Extracting data from indivdual Warm Core Rings