NMFS Trawl EDA

NMFS Trawl Data

Select columns of interest
Select stratum for full offshore survey
Eliminate strata not sampled for full time span
Select fall and spring surveys only
Select years 1970-2017
Remove biomass NAs

trawl.data<-survdat %>%
  select(ID, EST_YEAR,SEASON, STRATUM, DECDEG_BEGLAT,DECDEG_BEGLON,
         SVSPP, COMNAME, CATCHSEX,BIOMASS, AVGDEPTH, ABUNDANCE) %>%
  filter(STRATUM >= 01010 & STRATUM <= 01760) %>%
  filter(STRATUM!=1310 & STRATUM!=1320 & STRATUM!=1330 & STRATUM!=1350 &
           STRATUM!=1410 & STRATUM!=1420 & STRATUM!=1490) %>%
  filter(SEASON == "SPRING" | SEASON == "FALL") %>%
  filter(EST_YEAR >= 1970 & EST_YEAR < 2018) %>%
  filter(!is.na(BIOMASS))

Generate template of all sampling occasions by selecting distinct cases by ID, EST_YEAR, SEASON, STRATUM, and AVGDEPTH.

template<-trawl.data %>%
  distinct(ID, EST_YEAR, SEASON, STRATUM, AVGDEPTH,.keep_all = TRUE) %>%
  select(ID, EST_YEAR, SEASON, STRATUM,DECDEG_BEGLAT,DECDEG_BEGLON,AVGDEPTH)

Summarize biomass for each spp in each sample

samples<-trawl.data %>%
  distinct(ID, SVSPP, CATCHSEX,.keep_all = TRUE)

samples.biomass<-samples %>%
  group_by(ID,SVSPP) %>%
  summarise(BIOMASS = sum(BIOMASS))

Calculate center lat/long.

trawl.spp<-sort(unique(samples$SVSPP))
year <- seq(1970, 2017)
center7017 <- data.frame()

for (s in 1:length(trawl.spp)){ # loop for species
  
  spp.sub <- samples[samples$SVSPP==trawl.spp[s],]   # subset data for species
  
 
  for (y in 1:length(year)){  # loop for years
    
    spp.yr<- spp.sub[spp.sub$EST_YEAR == year[y], ]    # subset data for time
    
    biomass.total.year <- sum(spp.yr$BIOMASS)
    
    center.lat <- sum(spp.yr$DECDEG_BEGLAT*(spp.yr$BIOMASS/biomass.total.year))
    center.lon <- sum(spp.yr$DECDEG_BEGLON*(spp.yr$BIOMASS/biomass.total.year))
    center7017[(s-1)*(length(year))+y, 1]<- trawl.spp[s]
    center7017[(s-1)*(length(year))+y, 2]<- year[y]
    center7017[(s-1)*(length(year))+y, 3] <- center.lat
    center7017[(s-1)*(length(year))+y, 4] <- center.lon
    
    colnames(center7017)<- c("SVSPP","YEAR","CENTER_LAT", "CENTER_LON")
  }
  
}

Screen species.

screen7017 <- data.frame()

for (s in 1:length(trawl.spp)){ # loop for species
  spp.sub <- center7017[center7017$SVSPP==trawl.spp[s],] 
  count<-nrow(spp.sub[spp.sub$CENTER_LAT<=1,])  #THIS INDICATES YEARS IN WHICH THEY WEREN'T OBSERVED
  screen7017[s, 1]<- trawl.spp[s]
  screen7017[s, 2]<- count
  
  colnames(screen7017)<- c("SVSPP","MISSING_YRS")
}

nearfull7017<-screen7017[screen7017$MISSING_YRS==0,]

ts7017<-merge(center7017,nearfull7017,by="SVSPP",all.x=FALSE)
names<-unique(trawl.data[,c('SVSPP', 'COMNAME')])
ts7017<-merge(ts7017,names, by="SVSPP",all.y=FALSE)
ts7017<-ts7017[with(ts7017,order(SVSPP, YEAR)),]

Plot ts of faster moving species center.

##  1970s  1980s  1990s  2000s  2010s   NA's 
## 321957 294232 323834 565113  54704 300924

Distribution by species/decade

Create figures (biomass by latitude over decades, biomass weighted density by decade, and slope of center, leading, trailing edge).

Top plot: This barplot (not histogram) illustrates the relationship between latitude and biomass over time. This does not only include center of biomass, but latitude for all catch by species over each decade. I used a color scale to define the center 50% of the data, as well as the 80, 85, 90, and 95 percentiles.

Bottom left: This ggridge shows the biomass weighted density by decade. This is more like a histogram (counts), but I weighted it by biomass so it is more representative of actual catch.

Bottom right: In this figure, the points represent the annual center of biomass, the black line is the lm of center of biomass latitude ~ year (with the 99% CI of the lm). The dark grey dashed line represents the IQR of the species’ latitude distribution. The light grey dashed line represents mean plus and minus two standard deviations

Quartiles and standard deviations were calculated by using the summarize function from the FSA package in R on the lm of (lat ~ year) for each species (all data, not just center).

The dashed lines are a lm for q1.lat~year, q3.lat~year, and so on.

Potential “edges” by year

## 1970s 1980s 1990s 2000s 2010s  NA's 
## 60114 57446 65302 93485  8509 56611

Annual differences between CI values.

diffdf<-do.call(rbind,difflist) %>%
  as_tibble() %>%
  select(EST_YEAR, SPECIES, starts_with("mean")) %>%
  gather(key = "range","difference",-EST_YEAR, -SPECIES) %>%
  group_by(EST_YEAR, range) 

meanlabels<-as_labeller(c('mean1' = "Difference between 99%-95%", 'mean5' = 
                            "Difference between 95%-90%", 'mean10' = "Difference between 99%-90%"))

ggplot(diffdf, aes(EST_YEAR,difference, group = EST_YEAR)) + stat_boxplot(geom ='errorbar', width = 0.5) +
  geom_boxplot(outlier.shape = NA) + facet_wrap(~fct_relevel(range,"mean1","mean5","mean10"), ncol = 1, labeller = meanlabels) +   labs(x = "Year", y = "Degrees latitude", title = "All species by year") +
  theme(panel.grid = element_blank())

uniqueonly<-trawl.data %>% 
  distinct(ID, .keep_all = TRUE) %>% 
  filter(EST_YEAR == 2017)

bbox<-c(left = -76, bottom = 35, right = -65, top = 45)

library(ggmap)

## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.

## Please cite ggmap if you use it! See citation("ggmap") for details.

#create same density map, but animated by monthly density
ggmap(get_stamenmap(bbox, zoom = 5,maptype = "toner-lite")) + geom_point(aes(x = DECDEG_BEGLON, y = DECDEG_BEGLAT, colour = as.factor(STRATUM)), data = uniqueonly, size = 1.25, alpha = 0.9) + theme(legend.position = "none")

## Source : http://tile.stamen.com/toner-lite/5/9/11.png

## Source : http://tile.stamen.com/toner-lite/5/10/11.png

## Source : http://tile.stamen.com/toner-lite/5/9/12.png

## Source : http://tile.stamen.com/toner-lite/5/10/12.png

ggsave("trawl2017.pdf",plot = last_plot())

## Saving 7 x 5 in image

NMFS Trawl EDA

LGCarlson

5/20/2019

NMFS Trawl Data

Distribution by species/decade

Potential “edges” by year