Review alterations from updated CalCurr and SS3 work

R/SS_write_comps.R

@@ -9,20 +9,28 @@
 #'@param bins a list item, each item is a vector specifying either the age or length bins
 #'@param caal_bool a vector of length \code{ss_data_list} with boolean values for whether the data is conditional age-at-length
 #'@return modified ss_data_list that includes added composition data
-SS_write_comps <- function(ss_data_list, comp_matrix,
-                           data_rows, sampling_month,
-                           data_type, fleet_number,
-                           bins, caal_bool){
+
+#' HAP - 2024 - The following script was edited by Holly Ann Perryman to get updated CalCurr and SAR SS files to run
+
+SS_write_comps <- function(ss_data_list,
+                           comp_matrix,
+                           data_rows,
+                           sampling_month,
+                           data_type,
+                           fleet_number,
+                           bins,
+                           caal_bool){
 
   #Check if age and length comp are in the data_type vector and remove all rows if so
   if("agecomp" %in% data_type){
-    ss_data_list$agecomp <- ss_data_list$agecomp[0,]
+    # HAP - changed this to only take the first 9 columns and thus dropping all a# columns to account for age bins in atlantis differing from the dat file
+    ss_data_list$agecomp <- ss_data_list$agecomp[0,c(1:9)]
     age_comp_ind <- which(data_type=="agecomp")
   }
 
   if("lencomp" %in% data_type){
-
-    ss_data_list$lencomp <- ss_data_list$lencomp[0,c("Yr","Seas","FltSvy","Gender","Part","Nsamp")]
+    # HAP - new column names
+    ss_data_list$lencomp <- ss_data_list$lencomp[0,c("year","month","fleet","sex","part","Nsamp")] # c("Yr","Seas","FltSvy","Gender","Part","Nsamp") # month
     len_comp_ind <- which(data_type=="lencomp")
   }
 
@@ -34,52 +42,46 @@ SS_write_comps <- function(ss_data_list, comp_matrix,
     # if(i==3){
     #   browser()
     # }
-    k <- switch(data_type[i],"lencomp"=k_lencomp,
-                "agecomp"=k_agecomp)
+    k <- switch(data_type[i], "lencomp"=k_lencomp, "agecomp"=k_agecomp)
+    # HAP - new column names
     indices <- (k:(k+length(data_rows[[i]])-1))
-    ss_data_list[[data_type[i]]][indices,"Yr"] <- data_rows[[i]]
-    ss_data_list[[data_type[i]]][indices,"Seas"] <- sampling_month[[i]]
-    ss_data_list[[data_type[i]]][indices,"FltSvy"] <- rep(fleet_number[i], length(data_rows[[i]]))
-
-    ss_data_list[[data_type[i]]][indices,"Gender"] <-
-      rep(0, length(data_rows[[i]]))
-
-    ss_data_list[[data_type[i]]][indices,"Part"] <-
-      rep(0,length(data_rows[[i]]))
-
-    ss_data_list[[data_type[i]]][indices,"Nsamp"] <-
-      comp_matrix[[i]][,"nsamp"]
-
-    type_prefix <- switch(data_type[i],
-                          "agecomp"="a",
-                          "lencomp"="l")
-
+    ss_data_list[[data_type[i]]][indices,"year"]  <- data_rows[[i]]
+    ss_data_list[[data_type[i]]][indices,"month"] <- sampling_month[[i]]
+    ss_data_list[[data_type[i]]][indices,"fleet"] <- rep(fleet_number[i], length(data_rows[[i]]))
+    ss_data_list[[data_type[i]]][indices,"sex"]   <- rep(0, length(data_rows[[i]]))
+    ss_data_list[[data_type[i]]][indices,"part"]  <- rep(0,length(data_rows[[i]]))
+    ss_data_list[[data_type[i]]][indices,"Nsamp"] <- comp_matrix[[i]][,"nsamp"]
+    type_prefix <- switch(data_type[i],"agecomp"="a","lencomp"="l")
 
     if(data_type[i]=="agecomp"){
       if(caal_bool[i]){
         ss_data_list[[data_type[i]]][indices,"Lbin_lo"] <- comp_matrix[[i]][,"lower.bins"]
         ss_data_list[[data_type[i]]][indices,"Lbin_hi"] <- comp_matrix[[i]][,"upper.bins"]
       } else{
-      ss_data_list[[data_type[i]]][indices,"Lbin_lo"] <- rep(-1,length(data_rows[[i]]))
-      ss_data_list[[data_type[i]]][indices,"Lbin_hi"] <- rep(-1,length(data_rows[[i]]))
+        ss_data_list[[data_type[i]]][indices,"Lbin_lo"] <- rep(-1,length(data_rows[[i]]))
+        ss_data_list[[data_type[i]]][indices,"Lbin_hi"] <- rep(-1,length(data_rows[[i]]))
       }
-      ss_data_list[[data_type[i]]][indices, "Ageerr"] <- rep(1,length(data_rows[[i]]))
+      ss_data_list[[data_type[i]]][indices, "ageerr"] <- rep(1,length(data_rows[[i]])) # HAP - new column name
     }
 
-  if(data_type[[i]]=="lencomp"){
-    if(indices[1]==1){
-     ss_data_list[[data_type[i]]][, paste(type_prefix,as.character(bins[[i]]), sep="")] <- 0
+    if(data_type[[i]]=="lencomp"){
+      if(indices[1]==1){
+        ss_data_list[[data_type[i]]][, paste(type_prefix,as.character(bins[[i]]), sep="")] <- 0
+      }
     }
-  }
-    ss_data_list[[data_type[i]]][indices,paste(type_prefix,as.character(bins[[i]]), sep="")] <- comp_matrix[[i]][,as.character(bins[[i]])]
+    # HAP - I had to edit the following line as it was making lencomp data char, rather than num or int, which was giving me issues
+    ss_data_list[[data_type[i]]][indices,paste(type_prefix,as.character(bins[[i]]), sep="")] <- comp_matrix[[i]][,as.character(bins[[i]]), with = FALSE]
+
+    # HAP - have to check for NAs, likely due to mismatching bins across fleets, and turn NAs to zeros
+    ss_data_list[[data_type[i]]][is.na(ss_data_list[[data_type[i]]])] <- 0
+
     if(data_type[i]=="agecomp"){
       k_agecomp <- k_agecomp+length(data_rows[[i]])
     } else{
-
       k_lencomp <- k_lencomp+length(data_rows[[i]])
     }
-    }
-    return(ss_data_list)
-
   }
 
+  return(ss_data_list)
+}
+

R/SS_write_ts.R

@@ -9,24 +9,30 @@
 #'@param  units a vector where each entry must be one of "numbers" or "biomass"
 #'@param fleets a list of vectors, in each vector  each entry is the fleet number
 #'@param data_type a vector with length = \code{length(ts_data)}, each entry must be either "CPUE" or "catch"
-SS_write_ts <- function(ss_data_list, ts_data,
-                 CVs, data_years,
-                 sampling_month, units,
-                 fleets,
-                 data_type){
+
+#' HAP - 2024 - The following script was edited by Holly Ann Perryman to get updated CalCurr and SAR SS files to run
+
+SS_write_ts <- function(ss_data_list,
+                        ts_data,
+                        CVs,
+                        data_years,
+                        sampling_month,
+                        units,
+                        fleets,
+                        data_type){
 
 
   names_cpue <- c("year","seas","index", "obs", "se_log")
   names_catch <- c("year","seas","fleet", "catch", "catch_se")
 
   #Clear existing data
   if("CPUE" %in% data_type){
-  ss_data_list$CPUE <- data.frame(matrix(ncol=5, nrow=0))
-  colnames(ss_data_list$CPUE) <- names_cpue
+    ss_data_list$CPUE <- data.frame(matrix(ncol=5, nrow=0))
+    colnames(ss_data_list$CPUE) <- names_cpue
   }
   if("catch" %in% data_type){
-  ss_data_list$catch <- data.frame(matrix(ncol=5, nrow=0))
-  colnames(ss_data_list$catch) <- names_catch
+    ss_data_list$catch <- data.frame(matrix(ncol=5, nrow=0))
+    colnames(ss_data_list$catch) <- names_catch
   }
 
 
@@ -44,28 +50,26 @@ SS_write_ts <- function(ss_data_list, ts_data,
       ts_data[[i]] <- round(ts_data[[i]]/1000,0)
       ss_data_list$units_of_catch[i] <- 2
       ss_data_list$fleetinfo$units[i] <- 2
-      ss_data_list$CPUEinfo$Units[i] <-2
+      ss_data_list$CPUEinfo$units[i] <-2
     } else{
       ts_data[[i]] <- round(ts_data[[i]],0)
       ss_data_list$units_of_catch[i] <- 1
       ss_data_list$fleetinfo$units[i] <- 1
-      ss_data_list$CPUEinfo$Units[i] <-1
-      }
+      ss_data_list$CPUEinfo$units[i] <-1
+    }
 
     indices <- (k:(k+length(data_years[[i]])-1))
-  ss_data_list[[data_type[i]]][indices,"year"] <- data_years[[i]]
-  ss_data_list[[data_type[i]]][indices, "seas"] <- sampling_month[[i]]
-
-  ss_data_list[[data_type[i]]][indices, col_names[4]] <- ts_data[[i]]
-  ss_data_list[[data_type[i]]][indices, col_names[5]] <- rep(CVs[i], length(indices))
-  ss_data_list[[data_type[i]]][indices, col_names[3]] <- rep(fleets[i], length(indices))
+    ss_data_list[[data_type[i]]][indices,"year"] <- data_years[[i]]
+    ss_data_list[[data_type[i]]][indices, col_names[2]] <- sampling_month[[i]]
+    ss_data_list[[data_type[i]]][indices, col_names[4]] <- ts_data[[i]]
+    ss_data_list[[data_type[i]]][indices, col_names[5]] <- rep(CVs[i], length(indices))
+    ss_data_list[[data_type[i]]][indices, col_names[3]] <- rep(fleets[i], length(indices))
 
-  if(data_type[i]=="CPUE"){
-    k_CPUE <- k_CPUE+length(data_years[[i]])
-  } else{
-
-    k_catch <- k_catch+length(data_years[[i]])
-  }
+    if(data_type[i]=="CPUE"){
+      k_CPUE <- k_CPUE+length(data_years[[i]])
+    } else{
+      k_catch <- k_catch+length(data_years[[i]])
+    }
   }
 
   return(ss_data_list)

R/calc_bio_age.R

@@ -0,0 +1,91 @@
+
+
+##' @title Biomass for age groups
+##' @param fgs data from the groups.csv file
+##' @param nc_out_path ncdf atlantis' output file
+##' @return dataframe: time polygon biomass agecl species
+##' @author Demiurgo (ReactiveAtlantis)
+##' @editor Holly Perryman - Oct. 2024
+
+calc_bio_age <- function(fgs,
+                         nc_out_path,
+                         bboxes)
+{
+  # # --- FOR DEBUGGING
+  # fgs = fgs
+  # nc_out_path = file.path(dir, nc_out)
+  # bboxes = boxes
+  # # ---
+
+  # subset group csv data based on verts (groups with _ResN, _StructN, and _Nums data)
+  fgs_age <- fgs[fgs$NumCohorts  > 1 & !fgs$InvertType %in%
+                   c('PWN', 'PRAWNS', 'PRAWN', 'CEP', 'MOB_EP_OTHER', 'SEAGRASS', 'CORAL', 'MANGROVE', 'MANGROVES', 'SPONGE'), ]
+  # extract nc data
+  nc.out <- ncdf4::nc_open(nc_out_path)
+  Time <- nc.out$dim$t$vals / (60*60*24) # get time, convert from seconds to days
+  #* HAP, functions in atlantisom expect Time to be the number of time steps and not the days
+  Time <- Time / (Time[2] - Time[1])
+  # convertion metrics
+  mg2t   <- 0.00000002  # mgC converted to wet weight in tonnes = 20/1000000000
+  x.cn   <- 5.7   	    # Redfield ratio of C:N 5.7
+  # define placeholder for function output data
+  func_out_bio <- NULL
+  # enter loop over fids
+  for(index_fgs in 1 : nrow(fgs_age)){ # looping through fids --- for debugging: index_fgs = 2
+    # message to user
+    print(paste0(round((index_fgs / nrow(fgs_age)) * 100,2),"% complete computing biomass - ",
+                 "processing group ",
+                 fgs_age[index_fgs,5]))
+    # define placeholder for output age-based output data
+    cohort <- NULL
+    # enter loop over cohorts to compute biomass
+    for(coh in 1 : fgs_age[index_fgs, 'NumCohorts']){ # looping through index_fgs groups --- for debugging: coh = 1
+      # get name for nc output
+      name.fg <- paste0(fgs_age$Name[index_fgs], coh)
+      # get biomass based on the nc output and convert units
+      b.coh   <- (ncdf4::ncvar_get(nc.out, paste0(name.fg, '_ResN'))  +
+                    ncdf4::ncvar_get(nc.out, paste0(name.fg, '_StructN')))  *
+        ncdf4::ncvar_get(nc.out, paste0(name.fg, '_Nums')) * mg2t * x.cn
+      #*
+      #* HAP, I am not sure if this needs layer data or not
+      #*
+      #* (1) provide layer data
+      #*
+      # # convert to data frame
+      # b.coh <- cbind(expand.grid(layer = 1:dim(b.coh)[1],
+      #                            polygon = 1:dim(b.coh)[2],
+      #                            time = 1:dim(b.coh)[3]),
+      #                atoutput = as.vector(b.coh)) # --- head(b.coh)
+      # #* Drop boundary boxes - HAP, I am having issues with NAs with sample_survey_biomass() due to NAs in boundary boxes
+      # b.coh <- b.coh[which(!(b.coh$polygon %in% bboxes)),]
+      # # fix time to Time
+      # b.coh$time <- Time
+      # # add column for age cohort
+      # b.coh$agecl <- coh
+      # # store data
+      # cohort  <- rbind(cohort, b.coh);remove(b.coh) # head(cohort)
+      #*
+      #* (2) collapse by layer
+      #*
+      # collapse by depth, leaving time:polygon
+      b.coh <- apply(b.coh, 2, colSums, na.rm = TRUE)
+      # convert to data frame
+      b.coh <- cbind(expand.grid(time = 1:dim(b.coh)[1],
+                                 polygon = 1:dim(b.coh)[2]),
+                     atoutput = as.vector(b.coh)) # HAP, the header needs to be atoutput as other atlantisom functions expect this
+      #* HAP, drop boundary boxes - I am having issues with NAs with sample_survey_biomass() due to NAs in boundary boxes
+      b.coh <- b.coh[which(!(b.coh$polygon %in% bboxes)),]
+      # fix time to Time
+      b.coh$time <- Time
+      # add column for age cohort
+      b.coh$agecl <- coh
+      # store data
+      cohort  <- rbind(cohort, b.coh);remove(b.coh) # head(cohort)
+    }
+    # generate function output
+    cohort$species <- paste0(as.character(fgs_age$Name[index_fgs]))
+    #
+    func_out_bio <- rbind(func_out_bio, cohort)
+  }
+  return(as.data.frame(func_out_bio))
+}