SIRATES_2022

Inital data visualizations

Pulling in shape files and bathymetry for maps

# Bringing in shape files and bathymetry for maps

wd = ("C:/Users/sarah.salois/Documents/github/ssalois1/calc_wcr_metrics")
US.areas <- rgdal::readOGR(wd, 'USA')

## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\sarah.salois\Documents\github\ssalois1\calc_wcr_metrics", layer: "USA"
## with 24 features
## It has 4 fields

proj4string(US.areas) <- CRS("+init=epsg:4326")
## bring in bathymetry
atl <- marmap::getNOAA.bathy(-73,-69, 39, 42)
atl = fortify.bathy(atl)
blues = colorRampPalette(brewer.pal(9,'Blues'))(25)
depths <- c(0,50,100,200,300,500,Inf)
blues2 <- blues[c(5,7,9,11,13,24)]
coast = map_data("world2Hires")
coast = subset(coast, region %in% c('USA'))
coast$lon = (360 - coast$long)*-1 

#plot(US.areas)

Data Wrangling

Manually wrangled lengths and net metrics to be in the long format. Note: I will need to go back and write code to wrangle the excel data from R for efficiency and reproducibility. After loading the pre-wrangled data, I convert bottom temperatures from F to C and depth in fathoms to meters. A quick QC-check on names of unique species.

# Still need to figure out pulling in a workbook 
#df <- read_excel('Illex_biological_data fm_2022_Salt_Water_Intrusion.xlsx')
lengths <- read.csv('sirates_lengths.csv')
net_metrics <- read.csv('sirates_net_metrics.csv')
wtml <- read.csv('sirates_length_weights_illex.csv')
unique(lengths$species)

##  [1] "loligo"              "butterfish"          "illex"              
##  [4] "rough_scad"          "whiting"             "spiny_dogfish"      
##  [7] "spotted_hake"        "unknown"             "monkfish"           
## [10] "moonfish"            "fourspot_flounder"   "gulfstream_flounder"
## [13] "golden_tilefish"     "tilefish"

# error in station G, guessing bt was actually 50, may want to make NA
net_metrics[26,10] <- 50.0
# convert to metrics and celcius 
net_metrics$bottom_temp_c <- (net_metrics$bottom_temp - 32) * 0.5556
net_metrics$headrope_depth_m <- (net_metrics$headrope_depth * 1.8288) 
net_metrics$footrope_depth_m <- (net_metrics$footrope_depth * 1.8288) 
# Subset based on two net positions 
net_set <- net_metrics %>% filter(net_position == 'net_set')
haul_back <- net_metrics %>% filter(net_position == 'haul_back')

Map of Stations

The following charts out a map of the location of each of the 8 tows.

# Set up the dataframe as sf object 
#net_set_sp = coordinates(net_set) <- c('lon', 'lat')
#net_set_sp = proj4string(net_set) <- CRS("+init=epsg:4326")
net_set_sf <- st_as_sf(net_set, coords = c(4,3))
# Create a color palette 
wespal <- wes_palette("Zissou1")

p = ggplot() + 
  geom_contour_filled(data = atl,
                      aes(x=x,y=y,z=-1*z),
                      breaks=c(0,50,100,250,500,Inf),
                      size=c(0.3)) + 
  scale_fill_manual(values = blues2, # 5:20
                    name = paste("Depth (m)"),
                    labels = depths) +
  geom_contour(data = atl,
               aes(x=x,y=y,z=-1*z),
               breaks=c(0,50,100,250,500,Inf),
               size=c(0.3), 
               col = 'lightgrey') +
  new_scale_fill() +
  geom_sf(data = net_set_sf, 
                aes(color = bottom_temp_c), size = 5, 
                show.legend = 'point') +
  scale_color_gradientn(colours = wespal) +
  ggrepel::geom_label_repel(data = net_set_sf,
    aes(label = station, geometry = geometry),
    stat = "sf_coordinates",
    min.segment.length = 0,
    label.size = NA, 
    #color = text_color
  ) + 
  geom_polygon(data = coast, aes(x=lon, y = lat, group = group), 
                             color = 'gray20', fill = "wheat3") +
  coord_sf(xlim = c(-72,-69.5), ylim = c(39.5,41.8), 
           datum = sf::st_crs(4326)) +
  xlab('Longitude') + 
  ylab('Latitude')+
 # guides(fill = 'none') + #color = FALSE, this was attempt to get rid of dots
  theme_bw()
# Plot map with labels
p + geom_text()

This is a zoomed in version. Note: Something changes when I add the second layer of color - depth legend adds points, have tried to fix by specifiying color = z but did not work + Is potentially addressed here may look into for future plotting

p2 = ggplot() + 
  geom_contour_filled(data = atl,
                      aes(x=x,y=y,z=-1*z, color = 1*z),
                      breaks=c(0,50,100,250,500,Inf),
                      size=c(0.3)) + 
  scale_fill_manual(values = blues2, # 5:20
                    name = paste("Depth (m)"),
                    labels = depths) +
  geom_contour(data = atl,
               aes(x=x,y=y,z=-1*z),
               breaks=c(0,50,100,250,500,Inf),
               size=c(0.3), 
               col = 'lightgrey') +
  new_scale_fill() +
  geom_sf(data = net_set_sf, 
          aes(color = bottom_temp_c), size = 5, 
          show.legend = 'point') +
  scale_color_gradientn(colours = wespal, 
                        name = paste("Bottom Temp (°C)")) +
  ggrepel::geom_label_repel(data = net_set_sf,
                            aes(label = station, geometry = geometry),
                            stat = "sf_coordinates",
                            min.segment.length = 0,
                            label.size = NA) + 
  geom_polygon(data = coast, aes(x=lon, y = lat, group = group), 
               color = 'gray20', fill = "wheat3") +
  coord_sf(xlim = c(-71.4,-70.5), ylim = c(40.0,41.0), 
           datum = sf::st_crs(4326)) +
  xlab('Longitude') + 
  ylab('Latitude')+
  guides(fill = 'none') + #color = FALSE, 
  theme_bw()

p2 + geom_text()

Interactive maps with Leaflet

Net position : Net set

## -- Leaflet maps -- ##

netset_xyz <- net_set %>% 
  st_as_sf(coords = c('lon', 'lat'), crs = 4326)
pal <- colorNumeric(palette = "magma", netset_xyz$bottom_temp_c)

# Net set map
netset_xyz %>% 
  st_transform(4326) %>% 
  as("Spatial") %>% 
  leaflet() %>% 
  #addTiles() %>% 
  addProviderTiles("Esri.WorldStreetMap") %>%
  addCircles(color = ~pal(bottom_temp_c),
             popup = paste0(
               "<strong>Station: </strong>", 
               netset_xyz$station, "<br>",
               "<strong>Net Position: </strong>", 
               netset_xyz$net_position, "<br>",
               "<strong>Headrope Depth (m): </strong>",
               netset_xyz$headrope_depth_m, "<br>",
               "<strong>Footrope Depth (m): </strong>", 
               netset_xyz$footrope_depth_m, "<br>",
               "<strong>Temperature (°C): </strong>", 
               netset_xyz$bottom_temp_c, "<br>"
             )) %>% 
  addLegend(pal = pal, values = ~bottom_temp_c, 
            title = "Bottom temperature (°C)", opacity = 1)

Net position : Haul back

# Haul back map
haulback_xyz <- haul_back %>% 
  st_as_sf(coords = c('lon', 'lat'), crs = 4326)
pal <- colorNumeric(palette = "magma", haulback_xyz$bottom_temp_c)

haulback_xyz %>% 
  st_transform(4326) %>% 
  as("Spatial") %>% 
  leaflet() %>% 
  #addTiles() %>% 
  addProviderTiles("Esri.WorldStreetMap") %>%
  addCircles(color = ~pal(bottom_temp_c),
             popup = paste0(
               "<strong>Station: </strong>", 
               haulback_xyz$station, "<br>",
               "<strong>Net Position: </strong>", 
               haulback_xyz$net_position, "<br>",
               "<strong>Headrope Depth (m): </strong>", 
               haulback_xyz$headrope_depth_m, "<br>",
               "<strong>Footrope Depth (m): </strong>", 
               haulback_xyz$footrope_depth_m, "<br>",
               "<strong>Temperature (°C): </strong>", 
               haulback_xyz$bottom_temp_c, "<br>"
             )) %>% 
  addLegend(pal = pal, values = ~bottom_temp_c, 
            title = "Bottom temperature (°C)", opacity = 1)

Boxplots

All species:

## -- Boxplots -- ##
# All species
p = ggplot(lengths, aes(x = factor(station), y = length, col = species))+ 
  labs(x = 'Station', y = 'Lengths (mm)')+
  scale_color_brewer(palette = "Dark2") +
  theme_classic() + 
  #facet_wrap(~factor(fleet))
  geom_boxplot()
p

Illex vs Loligo

# Illex vs Loligo
p = ggplot(lengths %>% filter(species %in% c('loligo', 'illex')), 
           aes(x = factor(station), y = length, col = species))+ 
  labs(x = 'Station', y = 'Lengths (mm)') +
  scale_color_brewer(palette = "Dark2") +
  theme_classic() + 
  #facet_wrap(~factor(fleet))
  geom_boxplot()


p = p +  theme(text = element_text(size = 15), 
               axis.text=element_text(size = 16)) 
p

Just Illex

# Just Illex 
p = ggplot(lengths %>% filter(species == 'illex'), 
           aes(x = factor(station), y = length, col = species))+ 
  labs(x = 'Station', y = 'Lengths (mm)') +
  scale_color_brewer(palette = "Dark2") +
  theme_classic() + 
  #facet_wrap(~factor(fleet))
  geom_boxplot()

p = p +  theme(text = element_text(size = 15), 
               axis.text=element_text(size = 16)) 
p

# Just Illex 
p = ggplot(wtml , 
           aes(x = factor(station), y = weight*1000, col = factor(station)))+ 
  labs(x = 'Station', y = 'Weights (g)') +
  scale_color_brewer(palette = "Dark2") +
  theme_classic() + 
  #facet_wrap(~factor(fleet))
  geom_boxplot()

p = p +  theme(text = element_text(size = 15), 
               axis.text=element_text(size = 16)) 
p

Illex catch length/weight distributions

Lengths

Ridge plot for quick glance of lengths across stations

ggplot(lengths %>% filter(species == 'illex'), 
       aes(x = length, y = factor(station), fill =  factor(station))) +
  geom_density_ridges(alpha = .8, color = 'white',
                      scale = 2.5, rel_min_height = .01) +
  labs(title = '2022 Illex Lengths', x = 'Length (mm)', y = 'Station', fill = 'Station') +
  guides(color = guide_legend(title = 'Station')) +
  theme_ridges() +
  theme_classic()

ggplot(lengths %>% filter(species == 'illex'),
       aes(x = length, y = factor(station), 
                          fill = stat(quantile))) +
  stat_density_ridges(quantile_lines = FALSE,
                      calc_ecdf = TRUE, rel_min_height = 0.0001,
                      geom = "density_ridges_gradient") +
  labs(title = 'Illex lengths',
       subtitle = 'SIRATES',
       x = 'Length (mm)', y = 'Station') +
  scale_fill_brewer(name = '', palette = 'Greys') +
  ecodata::theme_ts()

## Picking joint bandwidth of 3.15

Weights

Ridge plot for quick glance of weights across stations

ggplot(wtml,
       aes(x = weight, y = factor(station), 
                          fill = stat(quantile))) +
  stat_density_ridges(quantile_lines = FALSE,
                      calc_ecdf = TRUE, rel_min_height = 0.0001,
                      geom = "density_ridges_gradient") +
  labs(title = 'Illex weights',
        x = 'Weights (kg)', y = 'Station') +
  scale_fill_brewer(name = '', palette = 'Greys') +
  ecodata::theme_ts()

## Picking joint bandwidth of 0.00301

Histograms of lengths across locations

# Just Illex 
ggplot(lengths %>% filter(species == 'illex'), 
       aes(x = length)) +
  geom_histogram(aes(y = ..density..), 
                 colour = "black",
                 fill = "white") +
  geom_density(alpha = 0.2, fill = "#FF6666") +
  #labs(title = '2021') +
  xlab('Length (mm)') +
  ylab('Density') + 
  theme_classic() + 
  facet_wrap(~factor(station))

Histogram of weights across locations

# Just Illex 
ggplot(wtml,
       aes(x = weight)) +
  geom_histogram(aes(y = ..density..), 
                 colour = "black",
                 fill = "white") +
  geom_density(alpha = 0.2, fill = "#FF6666") +
  xlab('Weight (kg)') +
  ylab('Density') + 
  theme_classic() + 
  facet_wrap(~factor(station))

Inspecting trends

Across all locations

#Plot series with trend 

g <- ggplot(wtml, aes(length, weight)) +           
  geom_point() + 
  theme_classic() 

g + stat_smooth(method = "lm",
              formula = y ~ x,
              geom = "smooth")

g + stat_smooth(method = "loess",
              formula = y ~ x,
              geom = "smooth") 

g <- ggplot(wtml, aes(length, weight)) +           
  geom_point() + 
  theme_classic() + 
  facet_wrap(~factor(station))

g+ stat_smooth(method = "loess",
              formula = y ~ x,
              geom = "smooth") 

Relative condition (Kn) = W/Wwhere W is the predicted length-specific mean weight for the fish population in a given region

# Condition
wtml <- wtml %>% mutate(weight = weight * 1000)
wtml$W <- 0.04810*(wtml$length*.10)^2.71990
wtml$Kn <- wtml$weight/wtml$W
wtml$W21 <- 0.082888*(wtml$length*.10)^2.530143
wtml$Kn21 <- wtml$weight/wtml$W21

ggplot(wtml, aes(x = factor(station), y = Kn, fill = factor(station))) + 
  labs(x = 'Station', y = 'Relative Condition (Kn)',
       title = 'Relative Condition:',
       subtitle = 'Based on 1970s coefficients') +
  theme_classic() +
  guides(fill=guide_legend(title=NULL))+
  geom_hline(yintercept = 1, lty = 2) +
  ylim(0,2) +
  geom_boxplot()

## Warning: Removed 4 rows containing non-finite values (stat_boxplot).

ggplot(wtml, 
            aes(x = factor(station), y = Kn21, fill = factor(station))) + 
  labs(x = 'Station', y = 'Relative Condition (Kn)',
       title = 'Relative Condition:',
       subtitle = 'Based on 2021/2022 coefficients') +
  theme_classic() +
  guides(fill=guide_legend(title=NULL))+
  geom_hline(yintercept = 1, lty = 2) +
  ylim(0,2) +
  geom_boxplot()

## Warning: Removed 1 rows containing non-finite values (stat_boxplot).

# Ordered by median 
ggplot(wtml,
            aes(x = reorder(station, Kn, median), y = Kn, fill = factor(station))) +
  labs(x = 'Station', y = 'Relative Condition (Kn)',
       title = 'Relative Condition: ordered by median',
       subtitle = 'Based on 1970s coefficients') +
  theme_classic() +
  guides(fill=guide_legend(title=NULL))+
  geom_hline(yintercept = 1, lty = 2) +
  ylim(0,2) +
  geom_boxplot()

## Warning: Removed 4 rows containing non-finite values (stat_boxplot).

A little data wrangling to map mean weight/condition by location

# Add location based on net set lat lon
wtml$lat <- NA
wtml$lon <- NA

wtml <- wtml %>% relocate(lat, .before=length)
wtml <- wtml %>% relocate(lon, .before=length)
a <- subset(net_set, station == 'A',  select = c('lat', 'lon', 'station'))
b <- subset(net_set, station == 'B', select = c('lat', 'lon', 'station'))
c <- subset(net_set, station == 'C', select = c('lat', 'lon', 'station'))
d <- subset(net_set, station == 'D', select = c('lat', 'lon', 'station'))
e <- subset(net_set, station == 'E', select = c('lat', 'lon', 'station'))
f <- subset(net_set, station == 'F', select = c('lat', 'lon', 'station'))
g <- subset(net_set, station == 'G', select = c('lat', 'lon', 'station'))
h <- subset(net_set, station == 'H', select = c('lat', 'lon', 'station'))
locs1 <- which(wtml$station %in% a$station)
locs2 <- which(wtml$station %in% b$station)
locs3 <- which(wtml$station %in% c$station)
locs4 <- which(wtml$station %in% d$station)
locs5 <- which(wtml$station %in% e$station)
locs6 <- which(wtml$station %in% f$station)
locs7 <- which(wtml$station %in% g$station)
locs8 <- which(wtml$station %in% h$station)
wtml[locs1,2] <- a[1,1]
wtml[locs1,3] <- a[1,2]
wtml[locs2,2] <- b[1,1]
wtml[locs2,3] <- b[1,2]
wtml[locs3,2] <- c[1,1]
wtml[locs3,3] <- c[1,2]
wtml[locs4,2] <- d[1,1]
wtml[locs4,3] <- d[1,2]
wtml[locs5,2] <- e[1,1]
wtml[locs5,3] <- e[1,2]
wtml[locs6,2] <- f[1,1]
wtml[locs6,3] <- f[1,2]
wtml[locs7,2] <- g[1,1]
wtml[locs7,3] <- g[1,2]
wtml[locs8,2] <- h[1,1]
wtml[locs8,3] <- h[1,2]

wtml_stats2 <- wtml %>% 
  group_by(station) %>% 
  summarise(n = n(),
         wt_mean = mean(weight),
            ml_mean = mean(length),
            wt_median = median(weight), 
            ml_median = median(length),
            wt_sd = sd(weight), 
            ml_sd = sd(length), 
            wt_max = max(weight), 
            ml_max = max(length), 
            wt_min = min(weight), 
            ml_min = min(length), 
            kn_mean = mean(Kn), 
            kn_median = median(Kn), 
            lat = mean(lat), 
            lon = mean(lon)) 

Plotting condtion and mean weight in space

wespal2 <- wes_palette("Zissou1", 3, type = "continuous")


p3 = ggplot() + 
  geom_contour_filled(data = atl,
                      aes(x=x,y=y,z=-1*z, color = 1*z),
                      breaks=c(0,50,100,250,500,Inf),
                      size=c(0.3)) + 
  scale_fill_manual(values = blues2, # 5:20
                    name = paste("Depth (m)"),
                    labels = depths) +
  geom_contour(data = atl,
               aes(x=x,y=y,z=-1*z),
               breaks=c(0,50,100,250,500,Inf),
               size=c(0.3), 
               col = 'lightgrey') +
  new_scale_fill() +
  geom_point(data = wtml_stats2, 
            aes(x = lon, y = lat, color = kn_mean, size = wt_mean),
          show.legend = 'point') +
  scale_color_gradientn(colours = wespal2, 
                        name = ('Mean condition')) +
  ggrepel::geom_label_repel(data = net_set_sf,
                            aes(label = station, geometry = geometry),
                            stat = "sf_coordinates",
                            min.segment.length = 0,
                            label.size = NA) + 
  geom_polygon(data = coast, aes(x=lon, y = lat, group = group), 
               color = 'gray20', fill = "wheat3") +
  coord_sf(xlim = c(-71.4,-70.5), ylim = c(40.0,41.0), 
           datum = sf::st_crs(4326)) +
  xlab('Longitude') + 
  ylab('Latitude')+
  guides(fill = 'none') + #color = FALSE, 
  theme_bw()

p3 + geom_text()

## Warning: `show.legend` must be a logical vector.

Estimating age at weight

Code below computes age estimate for each weight Parameters of weight at age estimation per Hendrickson and Hart 2006 for males and females combined. The formula for Weight as a function of Age is: weight = alpha * age^beta, therefore age estimate is the inverse: (weight/alpha)^(1/beta)

Frequency tables

Generate weight/length frequencies

tabw = table(wtml$weight) 
wtml %>% group_by(weight) %>% tally()

## # A tibble: 10 x 2
##    weight     n
##     <dbl> <int>
##  1      5    23
##  2     10   294
##  3     20   290
##  4     30    33
##  5     40     5
##  6     50     2
##  7     70     1
##  8     80     1
##  9     90     1
## 10    120     1

#count(wtml, 'length') 
tabl = table(wtml$length) 
wtml %>% group_by(length) %>% tally()

## # A tibble: 20 x 2
##    length     n
##     <int> <int>
##  1     19     1
##  2     55     2
##  3     60     3
##  4     65     2
##  5     70    23
##  6     75    21
##  7     80    98
##  8     85    60
##  9     90   126
## 10     92     1
## 11     95    78
## 12    100   113
## 13    105    39
## 14    110    54
## 15    115    13
## 16    120    10
## 17    130     3
## 18    140     2
## 19    150     1
## 20    170     1

# CrossTable in R
CrossTable(wtml$weight, wtml$length, prop.t=TRUE, prop.r=TRUE, prop.c=TRUE)

Age of 70 days is approximately 10 weeks old, the xaxis then increases by 2 weeks Birth weeks of 18:22 correspond to May 2:May 30, 2022 23-25 June 6, June 20th

library(patchwork)

## 
## Attaching package: 'patchwork'

## The following object is masked from 'package:interp':
## 
##     area

## The following object is masked from 'package:raster':
## 
##     area

# estimate age in days
alpha<-1.12e-6
beta<-3.60
wtml$age.est<-(wtml$weight/alpha)^(1/beta)
tab.w = table(wtml$weight)
barplot(tab.w, main=paste("Weight frequency for Drysten sample in week 35 "),
        xlab = "Weight", ylab = "Frequency", col='deeppink')

tab.age = table(round(wtml$age.est))
barplot(tab.age, main=paste("Age frequency for Drysten sample in week 35 "),
        xlab = "Age in days", ylab = "Frequency", col='blue')

# Analysis of putative  Birth weeks
# Compute birth week as function of age at capture
wtml$week <- 35
wtml$birth.week<-round(wtml$week-wtml$age.est/7)  
b.week.table<-table(round(wtml$birth.week))
barplot(b.week.table, main=paste("Birthweek frequency for Drysten sample in week 35 "),
             xlab="Birth Week", ylab="Frequency", col="red")

ml_wt_22 <- read.csv('ml_wt_21_22.csv')
wt_22 <- ml_wt_22 %>%
  filter(PARAM_TYPE == 'WT') %>%
  rename(weight = PARAM_VALUE_NUM)



w1 = ggplot(wt_22 %>% filter(year == 2022),
       aes(x = weight, y = factor(week), fill =  factor(week))) +
  geom_density_ridges(alpha = .8, color = 'white',
                      scale = 2.5, rel_min_height = .01) +
  labs(title = '2022 Illex Weights', x = 'Weight (gm)', y = 'Week', fill = 'Week') +
  guides(color = guide_legend(title = 'Week')) +
  theme_ridges() +
  theme_classic()
w1

## Picking joint bandwidth of 4.34