Load packages

require(dplyr)
## Warning: package 'dplyr' was built under R version 4.1.1
require(zooper) # devtools::install_github("InteragencyEcologicalProgram/zooper", ref="v2.4.1")
require(lubridate)
## Warning: package 'lubridate' was built under R version 4.1.1
require(readr)
## Warning: package 'readr' was built under R version 4.1.1
require(tidyr)
## Warning: package 'tidyr' was built under R version 4.1.1
require(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.1.1
require(readxl)
## Warning: package 'readxl' was built under R version 4.1.1
require(stringr)
## Warning: package 'stringr' was built under R version 4.1.1
require(sf)
## Warning: package 'sf' was built under R version 4.1.3
source("functions/region_assigner.R")

Load core stations

download.file("https://portal.edirepository.org/nis/dataviewer?packageid=edi.522.7&entityid=71dd301f30a2bc2e40f5da573dde9f97", destfile = file.path(tempdir(), "zoop_station_lookup.csv"))

zoop_stations<-read_csv(file.path(tempdir(), "zoop_station_lookup.csv"))%>%
  filter(Core%in%c(1, 2))

Load zoop data

if(utils::packageVersion('zooper') != '2.4.1') {
    stop('zooper version 2.4.1 is required, please update zooper to the correct version: 
         devtools::install_github("InteragencyEcologicalProgram/zooper", ref="v2.4.1")')
}

zoop_data<-Zoopsynther(Data_type="Community", Sources="EMP", Time_consistency = TRUE)%>%
  mutate(Month=month(Date))
## [1] "No disclaimers here! Enjoy the clean data!"

Read in zoop mass conversions and groupings Also Read in the energy density data and join to the mass grouping

zoop_energy_data<-read.csv("data/zoop_energy_data/taxa_energy_density_dry.csv")

zoop_mass_group<-read_excel("data/data_in/Zoop Categories for SEM.xlsx", na = "NA")%>%
  mutate(Taxlifestage=paste(Taxname, Lifestage))%>%
  rename(Group=`SEM Category`)%>%
  mutate(Group=str_replace_all(Group, " ", "_"))%>%
  filter(!is.na(Carbon_mass_micrograms)) %>% 
  left_join(zoop_energy_data, by = "Taxname")%>%
  mutate(Taxlifestage=if_else(Group=="Mysids", "Mysida", Taxlifestage))%>%
  select(Group, Taxlifestage, Carbon_mass_micrograms, energy_density_j_per_g_dry_mass)

# List all the zoop groups
unique(zoop_mass_group$Group)
## [1] "Cladoceran"           "Herbivorous_Copepods" "Predatory_Copepods"  
## [4] "Rotifers"

Load Mysid biomass data

mysid_energy = zoop_energy_data %>% 
  filter(Taxonomic.Category == "Mysids") %>% 
  summarise(energy = mean(energy_density_j_per_g_dry_mass)) %>% 
  unlist()
  

zoop_mysid<-read_excel("data/data_in/1972-2020MysidBPUEMatrix.xlsx",
                       sheet="Mysid_BPUE_matrix_1972-2020", na = "NA",
                       col_types = c(rep("numeric", 4), "date", "text", "text", rep("text", 7), rep("numeric", 8)))%>%
  select(Date=SampleDate, Station=StationNZ, `Acanthomysis aspera`:Unidentified)%>%
  mutate(BPUE_mysid=rowSums(select(., -Date, -Station), na.rm=T))%>%
  mutate(BPUE_mysid=BPUE_mysid*1000, # Convert to ug
         Group="Mysids")%>%
  mutate(SampleID=paste("EMP", Station, Date))%>%
  select(SampleID, Group, BPUE_mysid)

Start processing the zoop data

zoop_data_mass<-zoop_data%>%
  mutate(Taxlifestage=str_remove(Taxlifestage, fixed("_UnID")))%>%
  filter(
    !(SizeClass=="Meso" & #eliminating species which are counted in meso and micro and retained better in the micro net from the meso calcs
        
        Taxlifestage%in%c("Asplanchna Adult", "Copepoda Larva","Cyclopoida Juvenile", "Eurytemora affinis Larva", "Harpacticoida Undifferentiated",
                          "Keratella Adult", "Limnoithona Adult", "Limnoithona Juvenile", "Limnoithona sinenesis Adult", "Limnoithona tetraspina
                                    Adult", "Oithona Adult", "Oithona Juvenile", "Oithona davisae Adult", "Polyarthra Adult","Pseudodiaptomus Larva", 
                          "Rotifera Adult", "Sinocalanus doerrii Larva", "Synchaeta Adult", "Synchaeta bicornis Adult", "Trichocerca Adult")) &
      
      !(SizeClass=="Micro" & #removing categories better retained in meso net from micro net matrix
          
          Taxlifestage%in%c("Cirripedia Larva", "Cyclopoida Adult", "Oithona similis")))%>%
  mutate(Taxlifestage=if_else(Order=="Mysida", "Mysida", Taxlifestage), # Need to summarise Mysid CPUE first before joining to Mysid BPUE
         Taxname=if_else(Order=="Mysida", "Mysida", Taxname))%>%
  group_by(across(-c(Species, Genus, Family, CPUE)))%>%
  summarise(CPUE=sum(CPUE), .groups="drop")%>%
  left_join(zoop_mass_group, by=c("Taxlifestage"))%>%
  mutate(BPUE=CPUE*Carbon_mass_micrograms,
         Group=if_else(Order=="Mysida", "Mysids", Group))%>%
  left_join(zoop_mysid, 
            by=c("Group", "SampleID"))%>%
  filter(!is.na(Group))%>% # This removes anyone without an assigned group
  mutate(BPUE=if_else(is.na(BPUE), BPUE_mysid, BPUE),
         energy_density_j_per_g_dry_mass = if_else(Group == "Mysids",
                                                   mysid_energy,
                                                   energy_density_j_per_g_dry_mass),
         JPUE=(BPUE/1000)*energy_density_j_per_g_dry_mass)

Summarize zoop effort by month and year

zoop_data_sum<-zoop_data_mass%>%
  select(Month, Year, SampleID)%>%
  distinct()%>%
  group_by(Month, Year)%>%
  summarise(N=n(), .groups="drop")

Plot heat map

ggplot(zoop_data_sum, aes(y=Month, x=Year, fill=N))+
  geom_tile()+
  scale_fill_viridis_c()+
  scale_y_continuous(breaks=1:12)+
  coord_cartesian(expand=FALSE)+
  theme_bw()

Summarize zoop effort by date and station

zoop_data_station<-zoop_data_mass%>%
  select(Month, Year, Station)%>%
  distinct()%>%
  mutate(type="present")%>%
  complete(Month, Year, Station, fill=list(type="missing"))%>%
  mutate(Date=dmy(paste("1", Month, Year)))

Plot zoop effort by date and station for annual metrics

ggplot(zoop_data_station%>%filter(Station%in%zoop_stations$StationNZ), aes(x=Date, y=Station, color=type))+
  geom_point()+
  scale_color_manual(values=c("firebrick3", "black"))+
  theme_bw()

Find zoop stations that can be used for monthly metrics

zoop_stations_month<-zoop_data_mass%>%
  filter(Year>=1995)%>%
  group_by(Station, Month)%>%
  summarise(N_years=n_distinct(Year), .groups="drop")%>%
  filter(N_years>=22)%>%
  group_by(Station)%>%
  summarise(N_months=n_distinct(Month))%>%
  filter(N_months==12)

Plot zoop effort by date and station for monthly metrics

ggplot(zoop_data_station%>%filter(Station%in%zoop_stations_month$Station), aes(x=Date, y=Station, color=type))+
  geom_point()+
  scale_color_manual(values=c("firebrick3", "black"))+
  theme_bw()

Create annual data set with final set of annual stations

zoop_data_annual<-zoop_data_mass%>%
  filter(Station%in%zoop_stations$StationNZ)%>%
  filter(!Station%in%c("NZEZ6", "NZEZ2", "NZD16", "NZD06", "NZ080", "NZ042"))%>% # Renove stations that aren't continuous in time
  filter(Month%in%3:11)

Create monthly data set with final set of monthly stations

zoop_data_month<-zoop_data_mass%>%
  filter(Station%in%zoop_stations_month$Station)

Export zoop data for analysis

Annual data

zoop_data_annual_final<-zoop_data_annual%>%
  region_assigner(analysis="annual")

zoop_data_annual_final_noregions<-zoop_data_annual_final%>%
  group_by(Year, Group)%>%
  summarise(across(c(BPUE, CPUE, JPUE), ~sum(.x, na.rm=T)), .groups="drop")%>%
  pivot_wider(names_from = Group, values_from = c(BPUE, CPUE, JPUE),
              names_glue = "{Group}_{.value}")

str(zoop_data_annual_final_noregions)
## tibble [49 x 16] (S3: tbl_df/tbl/data.frame)
##  $ Year                     : num [1:49] 1972 1973 1974 1975 1976 ...
##  $ Cladoceran_BPUE          : num [1:49] 358380 155219 545954 225884 206525 ...
##  $ Herbivorous_Copepods_BPUE: num [1:49] 2203974 1191548 1521217 705526 1260426 ...
##  $ Mysids_BPUE              : num [1:49] 3247402 3945872 5040697 2011625 3374640 ...
##  $ Predatory_Copepods_BPUE  : num [1:49] 221896 66391 59958 37254 34763 ...
##  $ Rotifers_BPUE            : num [1:49] 2426765 3264591 2589798 1829724 658250 ...
##  $ Cladoceran_CPUE          : num [1:49] 213259 168612 375922 200851 212726 ...
##  $ Herbivorous_Copepods_CPUE: num [1:49] 633291 377878 507249 315244 449380 ...
##  $ Mysids_CPUE              : num [1:49] 8151 14237 15419 11508 11202 ...
##  $ Predatory_Copepods_CPUE  : num [1:49] 66040 19759 17845 11088 10346 ...
##  $ Rotifers_CPUE            : num [1:49] 19070589 23030366 20429422 13870210 5455748 ...
##  $ Cladoceran_JPUE          : num [1:49] 6916387 3338574 11205737 4727240 4679148 ...
##  $ Herbivorous_Copepods_JPUE: num [1:49] 58115448 31377746 39902377 18090502 32585420 ...
##  $ Mysids_JPUE              : num [1:49] 6.70e+07 8.14e+07 1.04e+08 4.15e+07 6.96e+07 ...
##  $ Predatory_Copepods_JPUE  : num [1:49] 5080981 1520228 1372919 853044 796001 ...
##  $ Rotifers_JPUE            : num [1:49] 39865768 35116977 43948077 31147852 4574664 ...
write_csv(zoop_data_annual_final_noregions, "data/annual_averages/zoop_annual_noregions.csv")

zoop_data_annual_final_regions<-zoop_data_annual_final%>%
  group_by(Year, Group, Region)%>%
  summarise(across(c(BPUE, CPUE, JPUE), ~sum(.x, na.rm=T)), .groups="drop")%>%
  pivot_wider(names_from = Group, values_from = c(BPUE, CPUE, JPUE),
              names_glue = "{Group}_{.value}")

str(zoop_data_annual_final_regions)
## tibble [147 x 17] (S3: tbl_df/tbl/data.frame)
##  $ Year                     : num [1:147] 1972 1972 1972 1973 1973 ...
##  $ Region                   : chr [1:147] "North" "South" "West" "North" ...
##  $ Cladoceran_BPUE          : num [1:147] 6071 351399 910 8005 139483 ...
##  $ Herbivorous_Copepods_BPUE: num [1:147] 648704 305727 1249542 283215 160355 ...
##  $ Mysids_BPUE              : num [1:147] 1555159 591244 1100998 884064 508617 ...
##  $ Predatory_Copepods_BPUE  : num [1:147] 52847 150460 18589 15604 35538 ...
##  $ Rotifers_BPUE            : num [1:147] 391741 1497800 537224 372596 1343326 ...
##  $ Cladoceran_CPUE          : num [1:147] 5951 206360 948 8141 152309 ...
##  $ Herbivorous_Copepods_CPUE: num [1:147] 182914 87687 362690 89924 53217 ...
##  $ Mysids_CPUE              : num [1:147] 3486 1435 3229 2536 646 ...
##  $ Predatory_Copepods_CPUE  : num [1:147] 15728 44780 5532 4644 10577 ...
##  $ Rotifers_CPUE            : num [1:147] 3356471 11233529 4480588 2661045 8298719 ...
##  $ Cladoceran_JPUE          : num [1:147] 117790 6780188 18409 161694 3069959 ...
##  $ Herbivorous_Copepods_JPUE: num [1:147] 17204469 8025960 32885019 7466324 4196202 ...
##  $ Mysids_JPUE              : num [1:147] 32086385 12198675 22716038 18240197 10493898 ...
##  $ Predatory_Copepods_JPUE  : num [1:147] 1210089 3445240 425652 357302 813741 ...
##  $ Rotifers_JPUE            : num [1:147] 6994533 28596900 4274336 5167977 24776936 ...
write_csv(zoop_data_annual_final_regions, "data/annual_averages/zoop_annual_regions.csv")

Export annual station list to data/stations folder

## Use zoop_data_annual_final so the stations saved here match those used 
## in the annual and annual regional datasets saved above.
stations_final_annual <- zoop_data_annual_final %>% 
  select(Station, Latitude, Longitude) %>%
  distinct() %>%
  mutate(Survey="EMP zooplantkon")

stations_final_annual %>%
  write.csv(file=file.path("data/stations/stations_zoop_annual.csv"), row.names=FALSE)

Plot the final set of annual stations

ggplot()+
  geom_sf(data=deltamapr::WW_Delta%>%st_transform(crs=4326))+
  geom_point(data=stations_final_annual, aes(x=Longitude, y=Latitude), color="red")

Monthly data

zoop_data_month_final<-zoop_data_month%>%
  region_assigner(analysis="monthly")

zoop_data_month_final_regions <- zoop_data_month_final %>%
  group_by(Year, Month, Group, Region)%>%
  summarise(across(c(BPUE, CPUE, JPUE), ~sum(.x, na.rm=T)), .groups="drop")%>%
  pivot_wider(names_from = Group, values_from = c(BPUE, CPUE, JPUE),
              names_glue = "{Group}_{.value}")

str(zoop_data_month_final_regions)
## tibble [1,861 x 18] (S3: tbl_df/tbl/data.frame)
##  $ Year                     : num [1:1861] 1972 1972 1972 1972 1972 ...
##  $ Month                    : num [1:1861] 1 1 1 2 2 2 3 3 3 4 ...
##  $ Region                   : chr [1:1861] "North" "South" "West" "North" ...
##  $ Cladoceran_BPUE          : num [1:1861] 265.1 411.7 12.2 300.4 43.5 ...
##  $ Herbivorous_Copepods_BPUE: num [1:1861] 220 965 5103 281 207 ...
##  $ Mysids_BPUE              : num [1:1861] 16089 15963 8214 4509 15744 ...
##  $ Predatory_Copepods_BPUE  : num [1:1861] 4576 7612 3345 4668 11090 ...
##  $ Rotifers_BPUE            : num [1:1861] 2659 7671 1871 2988 10424 ...
##  $ Cladoceran_CPUE          : num [1:1861] 358.52 635.81 6.21 236.46 66.67 ...
##  $ Herbivorous_Copepods_CPUE: num [1:1861] 61.2 272.9 1861.3 75.2 57.3 ...
##  $ Mysids_CPUE              : num [1:1861] 44.62 10.41 5.82 6.02 18.26 ...
##  $ Predatory_Copepods_CPUE  : num [1:1861] 1362 2266 996 1389 3301 ...
##  $ Rotifers_CPUE            : num [1:1861] 13529 41765 8824 24118 46471 ...
##  $ Cladoceran_JPUE          : num [1:1861] 5268 8862 249 6397 916 ...
##  $ Herbivorous_Copepods_JPUE: num [1:1861] 5776 25309 133537 7132 5416 ...
##  $ Mysids_JPUE              : num [1:1861] 331944 329356 169473 93035 324829 ...
##  $ Predatory_Copepods_JPUE  : num [1:1861] 104789 174302 76598 106894 253949 ...
##  $ Rotifers_JPUE            : num [1:1861] 51774 149366 36425 58189 202973 ...
write_csv(zoop_data_month_final_regions, "data/monthly_averages/zoop_month.csv")

Export monthly station list to data/stations folder

## Use zoop_data_month_final so the stations saved here match those used 
## in the monthly dataset saved above.
stations_final_monthly <- zoop_data_month_final %>% 
  select(Station, Latitude, Longitude) %>%
  distinct() %>%
  mutate(Survey="EMP zooplantkon")

stations_final_monthly %>%
  write.csv(file=file.path("data/stations/stations_zoop_month.csv"), row.names=FALSE)

Plot the final set of monthly stations

ggplot()+
  geom_sf(data=deltamapr::WW_Delta%>%st_transform(crs=4326))+
  geom_point(data=stations_final_monthly, aes(x=Longitude, y=Latitude), color="red")