Merge pull request #123 from NOAA-EDAB/Spatial_Biomass_Diagnostic

ATLNTS-373 Spatial biomass diagnostic

R/get_param_seasonal_migration.R

@@ -0,0 +1,63 @@
+#script to get seasonal migration parameters FXXX_SY
+
+get.move = function(bio.file,fgs,nbox = 30){
+  library(dplyr)
+
+  fgs = read.csv(fgs)
+
+  #Get groups and age class
+  groups = fgs$Code
+  is.age = fgs$NumCohorts > 1 
+
+  seasons = 1:4
+
+  #read in lines
+  bio.lines = readLines(bio.file)
+
+  find.line = function(x){
+    line.match = grep(paste0('^',x,'\\b'),bio.lines,ignore.case = T)
+    if(length(line.match) == 0){
+      return(NULL)
+    }else{
+      line.vals = as.numeric(strsplit(bio.lines[line.match + 1],'\t| ')[[1]])[1:nbox]  
+    }
+    return(line.vals)
+  }
+
+  #loop over species
+  spp.ls = list()
+  for(i in 1:length(groups)){
+
+    #construct parameters
+    if(is.age[i] == 1){
+      param.names = c(c(sapply(seasons, function(x) paste0('F',groups[i],'_S',x))),
+                      c(sapply(seasons, function(x) paste0('F',groups[i],'_S',x,'juv'))))
+      season.id = rep(1:4,2)
+    }else{
+      param.names = c(sapply(seasons, function(x) paste0('F',groups[i],'_S',x)))
+      season.id = rep(1:4)
+    }
+
+    spp.vals =  sapply(param.names,find.line)
+
+    if(all(sapply(spp.vals,is.null))){
+      spp.ls[[i]] = NULL
+    }else{
+      spp.ls[[i]] = spp.vals %>%
+        as.data.frame() %>%
+        mutate(box = 0:(nbox-1))%>%
+        reshape2::melt(id.vars = 'box')%>%
+        mutate(group = groups[i],
+               phase = ifelse(grepl('juv$',variable)&is.age[i]==1,'juv','adult'),
+               season = rep(season.id,each = nbox))
+    }
+    # print(groups[i])
+  }
+  return(dplyr::bind_rows(spp.ls))
+}
+
+move.params = get.move(fgs = here::here('currentVersion','neus_groups.csv'),
+                       bio.file = here::here('currentVersion','at_biology.prm')
+                       )
+write.csv(move.params, here::here('data-raw','seasonal_movements.csv'),row.names = F)
+

R/make_atlantis_diagnostic_figures.R

@@ -43,67 +43,6 @@
 #' Author: Ryan Morse, modified by J. Caracappa
 #' 
 
-# atl.dir = 'C:/Users/joseph.caracappa/Documents/Atlantis/ROMS_COBALT/Atlantis_Output_NutrientForcing/'
-# out.dir = 'C:/Users/joseph.caracappa/Documents/Atlantis/ROMS_COBALT/Atlantis_Output_NutrientForcing/Post_Processed/'
-# param.dir = here::here('CurrentVersion')
-# run.prefix = 'neus_output_test'
-# run.name = 'NutrientForcing'
-# 
-# 
-# source(here::here('R','get_atl_paramfiles.R'))
-# source(here::here('R','process_atl_output.R'))
-# param.ls= get_atl_paramfiles(param.dir,atl.dir,include_catch=T)
-# 
-# # process_atl_output( param.dir = here::here('CurrentVersion'),
-# #   atl.dir = 'C:/Users/joseph.caracappa/Documents/Atlantis/ROMS_COBALT/Atlantis_Output_NutrientForcing/',
-# #   out.dir = 'C:/Users/joseph.caracappa/Documents/Atlantis/ROMS_COBALT/Atlantis_Output_NutrientForcing/Post_Processed/',
-# #   run.prefix = 'neus_output_test',
-# #   include_catch = T,
-# #   save.out = T,
-# #   bgm.file = param.ls$bgm,
-# #   groups.file = param.ls$func.groups,
-# #   init.file = param.ls$init.nofill,
-# #   param.ls$biol.prm = param.ls$param.ls$biol.prm,
-# #   run.prm = param.ls$run.prm,
-# #   main.nc = param.ls$main.nc,
-# #   prod.nc = param.ls$prod.nc,
-# #   dietcheck = param.ls$dietcheck,
-# #   ssb = param.ls$ssb,
-# #   yoy = param.ls$yoy,
-# #   catch.file = param.ls$catch,
-# #   totcatch.file = param.ls$catchtot,
-# #   spatial.overlap = F
-# # )
-# 
-# load(paste0(out.dir,'neus_output_test_postprocessed.rdata'))
-# 
-# fig.width = 10
-# fig.height = 8
-# benthic.box = 10
-# benthic.level = 4
-# 
-# bgm.file = param.ls$bgm
-# param.ls$func.groups = param.ls$func.groups
-# run.prm = param.ls$run.prm
-# param.ls$biol.prm = param.ls$biol.prm
-# phytopl.history = here::here('R','phytoplankton_timeseries_biomass_tonnes_1998_2016.csv')
-# zoopl.history = here::here('R','Zooplankton_total_biomass_tonnes_N_20yrs.csv')
-# 
-# plot.benthic = T
-# plot.overall.biomass = T
-# plot.biomass.timeseries = T
-# plot.length.age=T
-# plot.biomass.box=T
-# plot.c.mum=T
-# plot.sn.rn=T
-# plot.recruits=T
-# plot.numbers.timeseries=T
-# plot.physics=T
-# plot.growth.cons=T
-# plot.cohort=T
-# plot.diet=T
-# plot.spatial.biomass=T
-# plot.LTL=T
 
 make_atlantis_diagnostic_figures = function(
   out.dir,
@@ -138,6 +77,7 @@ make_atlantis_diagnostic_figures = function(
   plot.diet,
   plot.consumption,
   plot.spatial.biomass,
+  plot.spatial.biomass.seasonal,
   plot.LTL,
   plot.catch, 
   plot.max.weight = T,
@@ -1129,4 +1069,15 @@ make_atlantis_diagnostic_figures = function(
     dev.off()
 
   }
+
+  if(plot.spatial.biomass.seasonal){
+    source(here::here('R','plot_biomass_box_summary.R'))
+
+    plot_biomass_box_season(bio.box = readRDS(paste0(out.dir,'biomass_box.rds')),
+                           bio.box.invert = readRDS(paste0(out.dir,'biomass_box_invert.rds')),
+                           fig.dir = fig.dir,
+                           species.list = NULL,
+                           plot.presence = T,
+                           save.fig = T)
+  }
 }

R/plot_biomass_box_summary.R

@@ -0,0 +1,231 @@
+#Functions to make box/epu species biomass presence/absence
+
+# Box Gradient ------------------------------------------------------------
+
+plot_biomass_box_season = function(bio.box,
+                                   bio.box.invert,
+                                   plot.presence=F,
+                                   species.list = NULL,
+                                   save.fig = T,
+                                   fig.dir,
+                                   tolerance = 0.1){
+  library(dplyr)
+  library(ggplot2)
+
+  box2epu = read.csv(here::here('Geometry','box2epu.csv'))
+  month2season = read.csv(here::here('data-raw','month2season.csv'))
+  season2name = data.frame(season = 1:4, season.name = unique(month2season$season.name))
+  fgs = read.csv(here::here('currentVersion','neus_groups.csv'))
+
+  move.param = read.csv(here::here('data-raw','seasonal_movements.csv'))%>%
+    group_by(group,box,season)%>%
+    summarise(orig = mean(value,na.rm=T))%>%
+    left_join(season2name)%>%
+    left_join(fgs,by = c('group' = 'Code'))%>%
+    rename(polygon = 'box',species = "LongName")%>%
+    ungroup()%>%
+    left_join(box2epu,by = c('polygon' = 'box'))%>%
+    select(species,polygon,epu,season.name,orig)%>%
+    mutate(orig.presence = ifelse(orig == 0,0,1))%>%
+    filter(!is.na(epu))
+
+
+  bio.all = rbind(bio.box,bio.box.invert)%>%
+    mutate(day = time * 365)
+
+  bio.all = bio.all %>%
+    mutate(date = as.POSIXct(day * 86400, origin = '1964-01-01 00:00:00',tz = 'UTC'),
+           month = as.numeric(format(date,format = '%m')))%>%
+    left_join(month2season)%>%
+    group_by(species,polygon,season.name)%>%
+    summarise(atoutput = mean(atoutput,na.rm=T))
+
+  if(!is.null(species.list)){
+    bio.all =  bio.all %>%
+      filter(species %in% species.list)
+  }
+
+  bio.all = bio.all %>%
+    group_by(species,polygon,season.name)%>%
+    summarise(atoutput = mean(atoutput,na.rm=T))
+
+  bio.max = bio.all %>%
+    group_by(species,season.name)%>%
+    summarise(atoutput.max = sum(atoutput,na.rm=T))
+
+  box.combs = expand.grid(species = unique(bio.all$species),
+                          polygon = sort(unique(bio.all$polygon)),
+                          season.name = unique(month2season$season.name))
+
+  bio.presence = box.combs %>%
+    left_join(bio.all)%>%
+    left_join(box2epu,by = c('polygon' = 'box'))%>%
+    mutate(presence = ifelse(atoutput == 0 | is.na(atoutput),0,1))%>%
+    left_join(bio.max)%>%
+    mutate(atoutput.scaled = atoutput/atoutput.max,
+           atoutput.scaled = ifelse(is.na(atoutput.scaled),0,atoutput.scaled))%>%
+    left_join(move.param)%>%
+    tidyr::replace_na(list(atoutput = 0))%>%
+    mutate(
+      presence.match = ifelse(orig.presence == presence,1,0),
+      prop.match = ifelse(orig >= (atoutput.scaled-tolerance) & orig <= (atoutput.scaled+tolerance), 1,0))
+
+  bio.presence$polygon = factor(bio.presence$polygon,
+                                levels = c(1:7,9,8,12:15,10:11,16:22))
+
+  if(plot.presence){
+
+    # plot.ls = list()
+    plot.name = paste0(fig.dir,'Box_EPU_season_presence.png')
+    p =   ggplot(bio.presence,aes(x=polygon,y=season.name,alpha = presence.match,fill = epu))+
+      geom_tile(color = 'black')+
+      facet_wrap(~species,ncol = 8)+
+      xlab('Box')+
+      guides(alpha = F)+
+      ylab('')+
+      theme_minimal()+
+      theme(panel.grid = element_blank(),
+            legend.position = 'bottom',
+            plot.title = element_text(hjust = 0.5))
+
+  }else{
+     plot.name = paste0(fig.dir,'Box_EPU_Season_scaled.png')
+     p =   ggplot(bio.presence,aes(x=polygon,y=season.name,alpha = atoutput.scaled,fill = epu))+
+       geom_tile(color = 'black')+
+       facet_wrap(~species,ncol = 8)+
+       xlab('Box')+
+       guides(alpha = F)+
+       ylab('')+
+       theme_minimal()+
+       theme(panel.grid = element_blank(),
+             legend.position = 'bottom',
+             plot.title = element_text(hjust = 0.5))
+
+  }
+
+  if(save.fig){
+    p + theme(axis.text.x = element_text(size = 6))
+    ggsave(plot.name,plot = p,width = 24, height = 10, units = 'in',dpi = 350)
+  }else{
+    return(p)
+  }
+
+}
+
+# Box Time Range ------------------------------------------------------------
+
+plot_biomass_box_range = function(bio.box,
+                                     bio.box.invert,
+                                     day.min=NA,
+                                     day.max=NA,
+                                     plot.presence = T,
+                                     species.list = NULL,
+                                     save.fig = T,
+                                     fig.dir){
+
+  library(dplyr)
+  library(ggplot2)
+
+  box2epu = read.csv(here::here('Geometry','box2epu.csv'))
+
+  bio.all = rbind(bio.box,bio.box.invert)%>%
+    mutate(day = time * 365)
+
+  bio.all = bio.all %>%
+    filter(day >= day.min & day <= day.max)
+
+  if(!is.null(species.list)){
+    bio.all =  bio.all %>%
+      filter(species %in% species.list)
+  }
+
+  bio.all = bio.all %>%
+    group_by(species,polygon)%>%
+    summarise(atoutput = mean(atoutput,na.rm=T))
+
+  bio.max = bio.all %>%
+    group_by(species)%>%
+    summarise(atoutput.max = max(atoutput,na.rm=T))
+
+  box.combs = expand.grid(species = unique(bio.all$species),
+                          polygon = sort(unique(bio.all$polygon)))
+
+  bio.presence = box.combs %>%
+    left_join(bio.all)%>%
+    left_join(box2epu,by = c('polygon' = 'box'))%>%
+    mutate(presence = ifelse(atoutput == 0 | is.na(atoutput),0,1))%>%
+    left_join(bio.max)%>%
+    mutate(atoutput.scaled = atoutput/atoutput.max,
+           atoutput.scaled = ifelse(is.na(atoutput.scaled),0,atoutput.scaled))
+
+  bio.presence$polygon = factor(bio.presence$polygon,
+                                levels = c(1:7,9,8,12:15,10:11,16:22))
+
+  spp.order = bio.presence %>%
+    group_by(species)%>%
+    summarise(presence.freq = mean(atoutput.scaled))%>%
+    arrange(desc(presence.freq))
+
+  bio.presence$species = factor(bio.presence$species,
+                                levels = rev(spp.order$species))
+
+  if(plot.presence){
+    p = ggplot(bio.presence,aes(x=polygon,y=species,alpha = presence,fill= epu))+
+      geom_tile(color = 'black')+
+      guides(alpha = F)+
+      xlab('Box')+
+      ylab('')+
+      scale_fill_manual(name = 'EPU', values = RColorBrewer::brewer.pal(3,'Set2'))+
+      theme_minimal()+
+      theme(panel.grid = element_blank(),
+            legend.position = 'bottom')
+
+    plot.name = paste0(fig.dir,'Box_EPU_summary.png')
+  }else{
+
+    p = ggplot(bio.presence,aes(x=polygon,y=species,alpha = atoutput.scaled,fill= epu))+
+      geom_tile(color = 'black')+
+      guides(alpha = F)+
+      xlab('Box')+
+      ylab('')+
+      scale_fill_manual(name = 'EPU', values = RColorBrewer::brewer.pal(3,'Set2'))+
+      theme_minimal()+
+      theme(panel.grid = element_blank(),
+            legend.position = 'bottom')
+
+    plot.name = paste0(fig.dir,'Box_EPU_summary.png')
+  }
+
+  if(save.fig){
+    p + ggsave(plot.name,width = 12, height = 12, units = 'in')
+  }else{
+    return(p)
+  }
+}
+
+
+# Example -----------------------------------------------------------------
+
+# plot_biomass_box_range(bio.box = readRDS('C:/Users/joseph.caracappa/Documents/Atlantis/Obs_Hindcast/Atlantis_Runs/ZM_Spatial_Final/Post_Processed/Data/biomass_box.rds'),
+#                        bio.box.invert = readRDS('C:/Users/joseph.caracappa/Documents/Atlantis/Obs_Hindcast/Atlantis_Runs/ZM_Spatial_Final/Post_Processed/Data/biomass_box_invert.rds'),
+#                        fig.dir = 'C:/Users/joseph.caracappa/Documents/Atlantis/Obs_Hindcast/Atlantis_Runs/ZM_Spatial_Final/Post_Processed/' ,
+#                        day.min = 366,
+#                        day.max = 730,
+#                        # species.list = c('Carnivorous zooplankton','Mesozooplankton','Microzooplankton','Gelatinous zooplankton'),
+#                        species.list = NULL,
+#                        plot.presence = F,
+#                        save.fig = T
+#                        
+# )
+# 
+# plot_biomass_box_season(bio.box = readRDS('C:/Users/joseph.caracappa/Documents/Atlantis/Obs_Hindcast/Atlantis_Runs/ZM_Spatial_Final/Post_Processed/Data/biomass_box.rds'),
+#                        bio.box.invert = readRDS('C:/Users/joseph.caracappa/Documents/Atlantis/Obs_Hindcast/Atlantis_Runs/ZM_Spatial_Final/Post_Processed/Data/biomass_box_invert.rds'),
+#                        fig.dir = 'C:/Users/joseph.caracappa/Documents/Atlantis/Obs_Hindcast/Atlantis_Runs/ZM_Spatial_Final/Post_Processed/' ,
+#                        # species.list = c('Carnivorous zooplankton','Mesozooplankton','Microzooplankton','Gelatinous zooplankton'),
+#                        species.list = NULL,
+#                        plot.presence = T,
+#                        save.fig = T
+# 
+# )
+
+