From 02169aaaf31d83ac284cbb65c205ed1fcf4482ec Mon Sep 17 00:00:00 2001 From: andybeet <22455149+andybeet@users.noreply.github.com> Date: Wed, 8 Jul 2026 07:43:58 -0400 Subject: [PATCH] style(air): format all files in R folder --- R/agg_custom.R | 12 +- R/compare_spatial_vars.R | 122 +++---- R/diag_fleet_catch.R | 28 +- R/diag_maxsize.R | 10 +- R/get_forcing_temperature.R | 2 +- R/get_param_FXXX_SY.R | 2 +- R/get_param_move_temp.R | 6 +- R/get_param_recruit_hdistrib.R | 2 +- R/get_param_recruit_vdistrib.R | 2 +- R/get_param_vert.R | 2 +- R/get_rawdata.R | 10 +- R/helperfiles.r | 8 +- R/make_atlantis_diagnostic_figures.R | 491 ++++++++++++++------------- R/plot_biomass_box_summary.R | 68 ++-- R/plot_overall_predation.R | 104 +++--- R/process_atl_output.R | 392 ++++++++++----------- R/process_catch_fleet.R | 44 +-- 17 files changed, 657 insertions(+), 648 deletions(-) diff --git a/R/agg_custom.R b/R/agg_custom.R index 919a697..0c833a1 100644 --- a/R/agg_custom.R +++ b/R/agg_custom.R @@ -7,7 +7,7 @@ #' @noRd agg_custom <- function(data, groups, fun, agg.scale) { - out = atlantistools::agg_data(data = data, groups = groups, fun = fun) %>% + out <- atlantistools::agg_data(data = data, groups = groups, fun = fun) %>% mutate( date = as.POSIXct( time * 365 * 86400, @@ -17,15 +17,15 @@ agg_custom <- function(data, groups, fun, agg.scale) { ) if (agg.scale == 'month') { - out$time.agg = as.numeric(factor(format(out$date, format = '%m'))) + out$time.agg <- as.numeric(factor(format(out$date, format = '%m'))) } else if (agg.scale == 'year') { - out$time.agg = as.numeric(factor(format(out$date, format = '%Y'))) + out$time.agg <- as.numeric(factor(format(out$date, format = '%Y'))) } else { - out$time.agg = out$time * 365 + out$time.agg <- out$time * 365 } - match.cols = c(groups[-which(groups == 'time')], 'time.agg') - out.agg = out %>% + match.cols <- c(groups[-which(groups == 'time')], 'time.agg') + out.agg <- out %>% group_by_at(match.cols) %>% dplyr::summarise( time = floor(min(time)), diff --git a/R/compare_spatial_vars.R b/R/compare_spatial_vars.R index 6902443..9267fcf 100644 --- a/R/compare_spatial_vars.R +++ b/R/compare_spatial_vars.R @@ -23,7 +23,7 @@ #' #' @export -compare_spatial_vars = function( +compare_spatial_vars <- function( param.dir, run.dirs, run.names, @@ -41,10 +41,10 @@ compare_spatial_vars = function( #Pull in reference data # ref.data = read.csv(ref.file, as.ist =T)%>% - ref.data = ref.data %>% + ref.data <- ref.data %>% dplyr::mutate(polygon = as.factor(polygon)) - init.data = init.data %>% + init.data <- init.data %>% dplyr::mutate(polygon = as.factor(polygon)) ##Check if ref.data$statistic == data.type @@ -58,23 +58,23 @@ compare_spatial_vars = function( #Pull in spatial data - boxes = atlantistools::convert_bgm(param.ls$bgm.file) %>% + boxes <- atlantistools::convert_bgm(param.ls$bgm.file) %>% dplyr::mutate(polygon = as.factor(polygon)) # ggplot2::ggplot(boxes,ggplot2::aes(x = long, y = lat, fill = polygon,group = polygon))+ # ggplot2::geom_polygon() - run.data.all.ls = list() + run.data.all.ls <- list() #Loop through each run.dirs for (i in 1:length(run.dirs)) { #Pull in run data if (tolower(ref.data$var.name[1]) == 'biomass') { - run.data = readRDS(paste0( + run.data <- readRDS(paste0( run.dirs[i], '/Post_Processed/Data/biomass_box.rds' )) - run.data.yr = run.data %>% + run.data.yr <- run.data %>% dplyr::mutate(polygon = as.factor(polygon)) %>% ##Calculate proportion by box over ref.years dplyr::filter(time >= ref.years[1] & time <= ref.years[2]) %>% @@ -85,11 +85,11 @@ compare_spatial_vars = function( dplyr::mutate(model.total = sum(atoutput)) %>% dplyr::ungroup() } else if (tolower(ref.data$var.name[1]) == 'numbers') { - run.data = readRDS(paste0( + run.data <- readRDS(paste0( run.dirs[i], '/Post_Processed/Data/numbers_box.rds' )) - run.data.yr = run.data %>% + run.data.yr <- run.data %>% dplyr::mutate(polygon = as.factor(polygon)) %>% ##Calculate proportion by box over ref.years dplyr::filter(time >= ref.years[1] & time <= ref.years[2]) %>% @@ -100,8 +100,8 @@ compare_spatial_vars = function( dplyr::mutate(model.total = sum(atoutput)) %>% dplyr::ungroup() } else if (tolower(ref.data$var.name[1]) == 'catch.total') { - run.data = readRDS(paste0(run.dirs[i], '/Post_Processed/Data/catch.rds')) - run.data.yr = run.data %>% + run.data <- readRDS(paste0(run.dirs[i], '/Post_Processed/Data/catch.rds')) + run.data.yr <- run.data %>% dplyr::mutate(polygon = as.factor(polygon)) %>% ##Calculate proportion by box over ref.years dplyr::filter(time >= ref.years[1] & time <= ref.years[2]) %>% @@ -114,11 +114,11 @@ compare_spatial_vars = function( dplyr::mutate(model.total = sum(atoutput)) %>% dplyr::ungroup() } else if (tolower(ref.data$var.name[1]) == 'catch_fleet') { - run.data = readRDS(paste0( + run.data <- readRDS(paste0( run.dirs[i], '/Post_Processed/Data/catch_fleet.rds' )) - run.data.yr = run.data %>% + run.data.yr <- run.data %>% dplyr::mutate(polygon = as.factor(polygon)) %>% ##Calculate proportion by box over ref.years dplyr::filter(time >= ref.years[1] & time <= ref.years[2]) %>% @@ -130,14 +130,14 @@ compare_spatial_vars = function( dplyr::ungroup() } - run.data.yr$polygon = as.factor(run.data.yr$polygon) + run.data.yr$polygon <- as.factor(run.data.yr$polygon) #Calculate proportion or value for model.val if (data.type == 'value') { - run.data.yr = run.data.yr %>% + run.data.yr <- run.data.yr %>% dplyr::mutate(model.val = atoutput) } else if (data.type == 'proportion') { - run.data.yr = run.data.yr %>% + run.data.yr <- run.data.yr %>% dplyr::mutate(model.val = atoutput / model.total) } else { warning("data.type must be either 'absolute' or 'proportion'") @@ -145,10 +145,10 @@ compare_spatial_vars = function( #Join with ref.data if (tolower(ref.data$var.name[1]) %in% c('catch', 'catch_fleet')) { - init.data2 = init.data %>% + init.data2 <- init.data %>% dplyr::select(species, polygon, init.value) %>% dplyr::mutate(polygon = as.factor(polygon)) - run.data.yr = run.data.yr %>% + run.data.yr <- run.data.yr %>% dplyr::left_join(ref.data) %>% dplyr::left_join(init.data2) %>% dplyr::mutate( @@ -156,7 +156,7 @@ compare_spatial_vars = function( statistic = ref.data$statistic[1] ) } else { - run.data.yr = run.data.yr %>% + run.data.yr <- run.data.yr %>% dplyr::left_join(ref.data) %>% dplyr::left_join(init.data) } @@ -165,22 +165,22 @@ compare_spatial_vars = function( #Do difference if (comparison.type == 'difference') { - run.data.compare = run.data.yr %>% + run.data.compare <- run.data.yr %>% dplyr::mutate(compare.val = model.val - ref.value) #Do scalars } else if (comparison.type == 'scalar') { - run.data.compare = run.data.yr %>% + run.data.compare <- run.data.yr %>% dplyr::mutate(compare.val = (model.val - ref.value) / ref.value) } else { warning('comparison.type must be either "difference" or "scalar"') } - run.data.all.ls[[i]] = run.data.compare %>% + run.data.all.ls[[i]] <- run.data.compare %>% dplyr::mutate(run.name = run.names[i]) } #Join to final dataframe with (species|polygon|atoutput|model.val|var.name|statistic|ref.value|compare.val|run.name|data.type|comparison.type) - run.data.all = dplyr::bind_rows(run.data.all.ls) %>% + run.data.all <- dplyr::bind_rows(run.data.all.ls) %>% dplyr::mutate(data.type = data.type, comparison.type = comparison.type) # out.name = paste(run.names,collapse = '_') @@ -189,10 +189,10 @@ compare_spatial_vars = function( if (plot == T) { #If plotting create PDF with the following plots - fgs = read.csv(param.ls$groups.file) - spp.names = sort(unique(run.data.all$species)) - spp.codes = fgs$Code[match(spp.names, fgs$LongName)] - box.id = sort(unique(boxes$polygon)) + fgs <- read.csv(param.ls$groups.file) + spp.names <- sort(unique(run.data.all$species)) + spp.codes <- fgs$Code[match(spp.names, fgs$LongName)] + box.id <- sort(unique(boxes$polygon)) # filter species of interest if (!is.null(speciesCodes)) { @@ -214,19 +214,19 @@ compare_spatial_vars = function( #Get species data #species ref box - ref.data.box = boxes %>% + ref.data.box <- boxes %>% dplyr::left_join(dplyr::filter(ref.data, species == spp.names[s])) - init.data.box = boxes %>% + init.data.box <- boxes %>% dplyr::left_join(dplyr::filter(init.data, species == spp.names[s])) - plot.data.ls = list() + plot.data.ls <- list() for (r in 1:length(run.names)) { - plot.data.spp = run.data.all %>% + plot.data.spp <- run.data.all %>% dplyr::filter(species == spp.names[s] & run.name == run.names[r]) #get missing boxes - missing.df = data.frame( + missing.df <- data.frame( species = spp.names[s], polygon = box.id[which(!(box.id %in% plot.data.spp$polygon))], atoutput = NA, @@ -242,17 +242,17 @@ compare_spatial_vars = function( dplyr::left_join(ref.data) %>% dplyr::left_join(init.data) - plot.data.spp = dplyr::bind_rows(plot.data.spp, missing.df) %>% + plot.data.spp <- dplyr::bind_rows(plot.data.spp, missing.df) %>% dplyr::mutate(polygon = as.factor(polygon)) - plot.data.ls[[r]] = boxes %>% + plot.data.ls[[r]] <- boxes %>% dplyr::left_join(plot.data.spp, by = 'polygon') } - plot.data = dplyr::bind_rows(plot.data.ls) + plot.data <- dplyr::bind_rows(plot.data.ls) - plot.spp.ls = list() + plot.spp.ls <- list() #1: Maps of ref values - p1 = ggplot2::ggplot( + p1 <- ggplot2::ggplot( ref.data.box, ggplot2::aes(x = long, y = lat, group = polygon, fill = ref.value) ) + @@ -269,7 +269,7 @@ compare_spatial_vars = function( ) #2: Map of init values - p2 = ggplot2::ggplot( + p2 <- ggplot2::ggplot( init.data.box, ggplot2::aes(x = long, y = lat, group = polygon, fill = init.value) ) + @@ -287,7 +287,7 @@ compare_spatial_vars = function( #3: Maps of run values - p3 = ggplot2::ggplot( + p3 <- ggplot2::ggplot( plot.data, ggplot2::aes(x = long, y = lat, group = polygon, fill = model.val) ) + @@ -305,7 +305,7 @@ compare_spatial_vars = function( ) #4: Maps of comparisons between runs and ref values - p4 = ggplot2::ggplot( + p4 <- ggplot2::ggplot( plot.data, ggplot2::aes(x = long, y = lat, group = polygon, fill = compare.val) ) + @@ -323,7 +323,7 @@ compare_spatial_vars = function( plot.title = ggplot2::element_text(hjust = 0.5) ) - plot.layout = matrix( + plot.layout <- matrix( c( 1, 1, @@ -357,9 +357,9 @@ compare_spatial_vars = function( #If Catch Fleet } else { - fisheries = read.csv(param.ls$fishery.prm) + fisheries <- read.csv(param.ls$fishery.prm) - fleet.combs = run.data.all %>% + fleet.combs <- run.data.all %>% dplyr::distinct(species, fleet) if (!is.null(speciesCodes)) { @@ -379,21 +379,21 @@ compare_spatial_vars = function( #Get species data #species ref box - ref.data.box = boxes %>% + ref.data.box <- boxes %>% dplyr::left_join(dplyr::filter( ref.data, species == fleet.combs$species[sf] & fleet == fleet.combs$fleet[sf] )) - init.data.box = boxes %>% + init.data.box <- boxes %>% dplyr::left_join(dplyr::filter( init.data, species == fleet.combs$species[sf] )) - plot.data.ls = list() + plot.data.ls <- list() for (r in 1:length(run.names)) { - plot.data.spp = run.data.all %>% + plot.data.spp <- run.data.all %>% dplyr::filter( species == fleet.combs$species[sf] & fleet == fleet.combs$fleet[sf] & @@ -401,9 +401,9 @@ compare_spatial_vars = function( ) #get missing boxes - missing.box = box.id[which(!(box.id %in% plot.data.spp$polygon))] + missing.box <- box.id[which(!(box.id %in% plot.data.spp$polygon))] if (length(missing.box) > 0) { - missing.df = data.frame( + missing.df <- data.frame( species = fleet.combs$species[sf], fleet = fleet.combs$fleet[sf], polygon = box.id[which(!(box.id %in% plot.data.spp$polygon))], @@ -420,21 +420,21 @@ compare_spatial_vars = function( dplyr::left_join(ref.data) %>% dplyr::left_join(init.data) - plot.data.spp = dplyr::bind_rows(plot.data.spp, missing.df) %>% + plot.data.spp <- dplyr::bind_rows(plot.data.spp, missing.df) %>% dplyr::mutate(polygon = as.factor(polygon)) } else { - plot.data.spp = plot.data.spp %>% + plot.data.spp <- plot.data.spp %>% dplyr::mutate(polygon = as.factor(polygon)) } - plot.data.ls[[r]] = boxes %>% + plot.data.ls[[r]] <- boxes %>% dplyr::left_join(plot.data.spp, by = 'polygon') } - plot.data = dplyr::bind_rows(plot.data.ls) + plot.data <- dplyr::bind_rows(plot.data.ls) - plot.spp.ls = list() + plot.spp.ls <- list() #1: Maps of ref values - p1 = ggplot2::ggplot( + p1 <- ggplot2::ggplot( ref.data.box, ggplot2::aes(x = long, y = lat, group = polygon, fill = ref.value) ) + @@ -451,7 +451,7 @@ compare_spatial_vars = function( ) #2: Map of init values - p2 = ggplot2::ggplot( + p2 <- ggplot2::ggplot( init.data.box, ggplot2::aes(x = long, y = lat, group = polygon, fill = init.value) ) + @@ -469,7 +469,7 @@ compare_spatial_vars = function( #3: Maps of run values - p3 = ggplot2::ggplot( + p3 <- ggplot2::ggplot( plot.data, ggplot2::aes(x = long, y = lat, group = polygon, fill = model.val) ) + @@ -487,7 +487,7 @@ compare_spatial_vars = function( ) #4: Maps of comparisons between runs and ref values - p4 = ggplot2::ggplot( + p4 <- ggplot2::ggplot( plot.data, ggplot2::aes(x = long, y = lat, group = polygon, fill = compare.val) ) + @@ -505,8 +505,8 @@ compare_spatial_vars = function( plot.title = ggplot2::element_text(hjust = 0.5) ) - spp.match = fgs$Code[which(fgs$LongName == fleet.combs$species[sf])] - plot.name = paste0( + spp.match <- fgs$Code[which(fgs$LongName == fleet.combs$species[sf])] + plot.name <- paste0( fleet.combs$fleet[sf], ':', fleet.combs$species[sf], @@ -517,7 +517,7 @@ compare_spatial_vars = function( ' ', data.type ) - plot.layout = matrix( + plot.layout <- matrix( c( 1, 1, diff --git a/R/diag_fleet_catch.R b/R/diag_fleet_catch.R index 675a36c..f1eab49 100644 --- a/R/diag_fleet_catch.R +++ b/R/diag_fleet_catch.R @@ -44,41 +44,41 @@ diag_fleet_catch <- function( min.dist = 100, relChangeThreshold = 0.01 ) { - boxes = atlantistools::convert_bgm(bgm) %>% + boxes <- atlantistools::convert_bgm(bgm) %>% dplyr::distinct(polygon, inside_lat, inside_long) - catch.fleet = process_catch_fleet( + catch.fleet <- process_catch_fleet( fishery.prm = fishery.prm, catch = catch.file, groups.file = fgs ) if (!is.null(speciesCodes)) { - fgs.df = read.csv(fgs, as.is = T) + fgs.df <- read.csv(fgs, as.is = T) - spp.match = fgs.df$LongName[which(fgs.df$Code %in% speciesCodes)] + spp.match <- fgs.df$LongName[which(fgs.df$Code %in% speciesCodes)] - catch.fleet = catch.fleet %>% + catch.fleet <- catch.fleet %>% dplyr::filter(species %in% spp.match) - catch.ref = catch.ref %>% + catch.ref <- catch.ref %>% dplyr::filter(species %in% spp.match) } - max.yr = max(catch.fleet$time) + max.yr <- max(catch.fleet$time) #mnagnitude - catch.mag.model = catch.fleet %>% + catch.mag.model <- catch.fleet %>% dplyr::filter(time >= (max.yr - nYrs)) %>% dplyr::group_by(species, fleet, time) %>% dplyr::summarise(catch.model = sum(atoutput, na.rm = T)) %>% dplyr::group_by(species, fleet) %>% dplyr::summarise(catch.model = mean(catch.model, na.rm = T)) - catch.mag.ref = catch.ref %>% + catch.mag.ref <- catch.ref %>% dplyr::group_by(species, fleet) %>% dplyr::summarise(catch.ref = mean(ref.value, na.rm = T)) - catch.mag.all = catch.mag.model %>% + catch.mag.all <- catch.mag.model %>% dplyr::left_join(catch.mag.ref) %>% dplyr::mutate( catch.rel = catch.ref / catch.model, @@ -91,7 +91,7 @@ diag_fleet_catch <- function( ) #distance - catch.cog.model = catch.fleet %>% + catch.cog.model <- catch.fleet %>% dplyr::filter(time >= (max.yr - nYrs)) %>% dplyr::left_join(boxes) %>% dplyr::mutate( @@ -111,7 +111,7 @@ diag_fleet_catch <- function( cog.y.model = mean(cog.y, na.rm = T) ) - catch.cog.ref = catch.ref %>% + catch.cog.ref <- catch.ref %>% dplyr::left_join(boxes) %>% dplyr::mutate( catch.wgt.x = inside_long * ref.value, @@ -125,7 +125,7 @@ diag_fleet_catch <- function( ) %>% dplyr::mutate(cog.x.ref = cog.x / catch.tot, cog.y.ref = cog.y / catch.tot) - catch.cog.all = catch.cog.model %>% + catch.cog.all <- catch.cog.model %>% dplyr::left_join(catch.cog.ref) %>% dplyr::mutate( dist = sqrt((cog.x.ref - cog.x.model)^2 + (cog.y.ref - cog.y.model)^2), @@ -133,7 +133,7 @@ diag_fleet_catch <- function( ) #combine for output - catch.diag.all = catch.mag.all %>% + catch.diag.all <- catch.mag.all %>% dplyr::left_join(catch.cog.all) %>% dplyr::select( species, diff --git a/R/diag_maxsize.R b/R/diag_maxsize.R index 01d2cd2..77f2999 100644 --- a/R/diag_maxsize.R +++ b/R/diag_maxsize.R @@ -80,9 +80,9 @@ diag_maxsize <- function( # list of variables to pull from main nc file. # Needed for biomass calculation. Each variable resides in list element - vars = list('Nums', 'StructN', 'ResN', 'N') - group.types = list(groups.age, groups.age, groups.age, groups.bp) - rawdata.main = Map( + vars <- list('Nums', 'StructN', 'ResN', 'N') + group.types <- list(groups.age, groups.age, groups.age, groups.bp) + rawdata.main <- Map( atlantistools::load_nc, select_variable = vars, select_groups = group.types, @@ -96,7 +96,7 @@ diag_maxsize <- function( ) # calculate biomass for species,age, polygon, layer, time - spatial.biomass = atlantistools::calculate_biomass_spatial( + spatial.biomass <- atlantistools::calculate_biomass_spatial( nums = rawdata.main[[1]], sn = rawdata.main[[2]], rn = rawdata.main[[3]], @@ -106,7 +106,7 @@ diag_maxsize <- function( bps = bio.pools ) # grab numbers in time and space - spatialNumbers = rawdata.main[[1]] %>% + spatialNumbers <- rawdata.main[[1]] %>% dplyr::rename(numbers = .data$atoutput) # filter biomass for species with 10 cohorts and convert to kilograms spatialBiomass <- spatial.biomass %>% diff --git a/R/get_forcing_temperature.R b/R/get_forcing_temperature.R index 724ec6f..4e22f13 100644 --- a/R/get_forcing_temperature.R +++ b/R/get_forcing_temperature.R @@ -98,7 +98,7 @@ get_forcing_temperature <- function(param.ls, plotFigs = F) { # round up to nearest decade for plotting maxTime <- max(tempD$time) + (10 - max(tempD$time) %% 10) for (ilayer in sort(unique(tempD$layer))) { - data = tempD |> + data <- tempD |> dplyr::filter(layer == ilayer) |> dplyr::select(-variable, -layer) if (nrow(data) == 0) { diff --git a/R/get_param_FXXX_SY.R b/R/get_param_FXXX_SY.R index 4ca839a..6270835 100644 --- a/R/get_param_FXXX_SY.R +++ b/R/get_param_FXXX_SY.R @@ -21,7 +21,7 @@ get_param_FXXX_SY <- function(bio.prm) { values <- gsub("\t", " ", values) # create a null dataframe - out.df = data.frame( + out.df <- data.frame( group = NULL, season = NULL, cohort = NULL, diff --git a/R/get_param_move_temp.R b/R/get_param_move_temp.R index 5490839..6780728 100644 --- a/R/get_param_move_temp.R +++ b/R/get_param_move_temp.R @@ -12,14 +12,14 @@ get_param_move_temp <- function(bio.prm) { # read n parameter file - bio.lines = readLines(bio.prm) + bio.lines <- readLines(bio.prm) # identify lines that contain the min_move_temp and max_move_temp values - line = grep(paste0("_move_temp"), bio.lines) + line <- grep(paste0("_move_temp"), bio.lines) # extract the lines lines <- bio.lines[line] # create a null dataframe - out.df = data.frame(group = NULL, value = NULL, limit = NULL) + out.df <- data.frame(group = NULL, value = NULL, limit = NULL) # loop through the lines and extract the values, groups and the min/max for (i in 1:length(lines)) { linesi <- gsub("_move_temp", "", lines[i]) diff --git a/R/get_param_recruit_hdistrib.R b/R/get_param_recruit_hdistrib.R index 44e2ec6..aed4941 100644 --- a/R/get_param_recruit_hdistrib.R +++ b/R/get_param_recruit_hdistrib.R @@ -21,7 +21,7 @@ get_param_recruit_hdistrib <- function(bio.prm) { values <- gsub("\t", " ", values) # create a null dataframe - out.df = data.frame(group = NULL, polygon = NULL, value = NULL) + out.df <- data.frame(group = NULL, polygon = NULL, value = NULL) # loop through the lines and extract the values, groups and the min/max for (i in 1:length(lines)) { # remove trailing part of the line where the number of values on next line is listed diff --git a/R/get_param_recruit_vdistrib.R b/R/get_param_recruit_vdistrib.R index e816b02..3dfed9a 100644 --- a/R/get_param_recruit_vdistrib.R +++ b/R/get_param_recruit_vdistrib.R @@ -23,7 +23,7 @@ get_param_recruit_vdistrib <- function(bio.prm) { values <- gsub("\t", " ", values) # create a null dataframe - out.df = data.frame(group = NULL, layer = NULL, value = NULL) + out.df <- data.frame(group = NULL, layer = NULL, value = NULL) # loop through the lines and extract the values, groups and the min/max for (i in 1:length(lines)) { # remove trailing part of the line where the number of values on next line is listed diff --git a/R/get_param_vert.R b/R/get_param_vert.R index 79b6cfb..a7f9ac9 100644 --- a/R/get_param_vert.R +++ b/R/get_param_vert.R @@ -39,7 +39,7 @@ get_param_vert <- function(bio.prm) { valuesD <- gsub("\t", " ", valuesD) # create a null dataframe - outdf = data.frame( + outdf <- data.frame( group = NULL, daynight = NULL, cohort = NULL, diff --git a/R/get_rawdata.R b/R/get_rawdata.R index d1d23c0..5f5ccfe 100644 --- a/R/get_rawdata.R +++ b/R/get_rawdata.R @@ -6,9 +6,9 @@ get_rawdata <- function(group, group.type) { if (group.type == 'age') { - main.vars = age.vars + main.vars <- age.vars } else { - main.vars = bp.vars + main.vars <- bp.vars } if ( @@ -17,7 +17,7 @@ get_rawdata <- function(group, group.type) { )] == 0 ) { - rawdata.spp.f = list(data.frame( + rawdata.spp.f <- list(data.frame( species = group.types$species[i], polygon = NA, agecl = NA, @@ -26,7 +26,7 @@ get_rawdata <- function(group, group.type) { atoutput = NA )) } else { - rawdata.spp.f = Map( + rawdata.spp.f <- Map( load_nc_temp, select_variable = main.vars, select_groups = group, @@ -41,7 +41,7 @@ get_rawdata <- function(group, group.type) { } if (group.types$group[i] == 'bp') { - rawdata.spp.f[[1]]$agecl = 1 + rawdata.spp.f[[1]]$agecl <- 1 } return(rawdata.spp.f) } diff --git a/R/helperfiles.r b/R/helperfiles.r index a957042..5006606 100644 --- a/R/helperfiles.r +++ b/R/helperfiles.r @@ -18,15 +18,15 @@ run.filename <- function(command, code) { #Utility function -bind.save = function(x, name, out.dir) { - x2 = dplyr::bind_rows(x) +bind.save <- function(x, name, out.dir) { + x2 <- dplyr::bind_rows(x) saveRDS(x2, file.path(out.dir, paste0(name, '.rds'))) } #Utility function -add.title = function(p, title) { - p = p + +add.title <- function(p, title) { + p <- p + ggplot2::ggtitle(title) + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) return(p) diff --git a/R/make_atlantis_diagnostic_figures.R b/R/make_atlantis_diagnostic_figures.R index dedc9f2..9424740 100644 --- a/R/make_atlantis_diagnostic_figures.R +++ b/R/make_atlantis_diagnostic_figures.R @@ -43,7 +43,7 @@ #' Author: Ryan Morse, modified by J. Caracappa #' @export -make_atlantis_diagnostic_figures = function( +make_atlantis_diagnostic_figures <- function( out.dir, fig.dir, atl.dir, @@ -81,31 +81,31 @@ make_atlantis_diagnostic_figures = function( fig.dir <- check_string(fig.dir) #Load groups data - group.code = atlantistools::get_age_acronyms(param.ls$groups.file) - group.data = atlantistools::load_fgs(param.ls$groups.file) - group.index = dplyr::select(group.data, c(Code, LongName)) + group.code <- atlantistools::get_age_acronyms(param.ls$groups.file) + group.data <- atlantistools::load_fgs(param.ls$groups.file) + group.index <- dplyr::select(group.data, c(Code, LongName)) #Load BGM data - box.bgm = atlantistools::load_box(param.ls$bgm) + box.bgm <- atlantistools::load_box(param.ls$bgm) #plot parameters - plot.labels = list(x = 'Time (years)', y = 'Biomass (tonnes)') + plot.labels <- list(x = 'Time (years)', y = 'Biomass (tonnes)') # Benthic box timeseries -------------------------------------------------- #Select box for timeseries of all benthic groups if (plot.benthic | plot.all) { print("benthic") - biomass.spatial.stanza = readRDS(file.path( + biomass.spatial.stanza <- readRDS(file.path( out.dir, 'biomass_spatial_stanza.rds' )) - benthic.biomass.spatial = dplyr::filter( + benthic.biomass.spatial <- dplyr::filter( biomass.spatial.stanza, layer == benthic.level & polygon == benthic.box ) - benthic.spp = unique(benthic.biomass.spatial$species) - box.area = box.bgm$boxes[[benthic.box + 1]]$area + benthic.spp <- unique(benthic.biomass.spatial$species) + box.area <- box.bgm$boxes[[benthic.box + 1]]$area pdf( file = file.path( @@ -114,8 +114,8 @@ make_atlantis_diagnostic_figures = function( ) ) for (i in 1:length(benthic.spp)) { - spp.ind = benthic.spp[i] - benthic.biomass.spp = dplyr::filter( + spp.ind <- benthic.spp[i] + benthic.biomass.spp <- dplyr::filter( benthic.biomass.spatial, species == spp.ind ) @@ -136,47 +136,47 @@ make_atlantis_diagnostic_figures = function( if (plot.catch | plot.all) { print("Catch") - catchmt = readRDS(file.path(out.dir, 'catchmt.rds')) + catchmt <- readRDS(file.path(out.dir, 'catchmt.rds')) #Catch by species time series (metric tonnes) - temp.plot.1 = atlantistools::plot_line(catchmt) - temp.plot.1 = ggplot2::update_labels(temp.plot.1, labels = plot.labels) - temp.plot.1 = ggplot2::update_labels( + temp.plot.1 <- atlantistools::plot_line(catchmt) + temp.plot.1 <- ggplot2::update_labels(temp.plot.1, labels = plot.labels) + temp.plot.1 <- ggplot2::update_labels( temp.plot.1, labels = list(x = 'Time (years)', y = 'Metric Tonnes') ) - temp.plot.1 = add.title(temp.plot.1, 'Catch') + temp.plot.1 <- add.title(temp.plot.1, 'Catch') #Catch at age time series (numbers) - totcatch = readRDS(file.path(out.dir, 'totcatch.rds')) + totcatch <- readRDS(file.path(out.dir, 'totcatch.rds')) - temp.plot.2 = atlantistools::plot_line(totcatch, col = 'agecl') - temp.plot.2 = ggplot2::update_labels( + temp.plot.2 <- atlantistools::plot_line(totcatch, col = 'agecl') + temp.plot.2 <- ggplot2::update_labels( p = temp.plot.2, labels = c(plot.labels, list(colour = 'Ageclas')) ) - temp.plot.2 = ggplot2::update_labels( + temp.plot.2 <- ggplot2::update_labels( temp.plot.2, labels = list(x = 'Time (years)', y = 'Numbers') ) - temp.plot.2 = add.title(temp.plot.2, 'Catch at Age') + temp.plot.2 <- add.title(temp.plot.2, 'Catch at Age') #Catch at age - percent - catch.age.pct = atlantistools::agg_perc( + catch.age.pct <- atlantistools::agg_perc( totcatch, groups = c('time', 'species') ) - temp.plot.6 = atlantistools::plot_bar( + temp.plot.6 <- atlantistools::plot_bar( catch.age.pct, fill = 'agecl', wrap = 'species' ) - temp.plot.6 = ggplot2::update_labels( + temp.plot.6 <- ggplot2::update_labels( temp.plot.6, labels = list(x = 'Time (years)', y = 'Numbers (%)') ) - temp.plot.6 = add.title(temp.plot.6, 'Catch at age - Percent') + temp.plot.6 <- add.title(temp.plot.6, 'Catch at age - Percent') pdf( file.path(fig.dir, paste0(run.name, ' Catch Timeseries.pdf')), @@ -197,16 +197,16 @@ make_atlantis_diagnostic_figures = function( if (plot.mortality | plot.all) { print("mortality") # plot mortality from Mort.txt - mort = readRDS(file.path(out.dir, 'mort.rds')) + mort <- readRDS(file.path(out.dir, 'mort.rds')) itype <- 1 plotMort <- list() # Annual Mortality time series M, F by species on same plot - temp.plot.1 = atlantistools::plot_line(mort, col = "source") - temp.plot.1 = ggplot2::update_labels( + temp.plot.1 <- atlantistools::plot_line(mort, col = "source") + temp.plot.1 <- ggplot2::update_labels( temp.plot.1, labels = list(x = 'Time (years)', y = 'Mortality') ) - temp.plot.1 = add.title(temp.plot.1, 'Mortality (F & M2)') + temp.plot.1 <- add.title(temp.plot.1, 'Mortality (F & M2)') plotMort[[itype]] <- temp.plot.1 # plot mortaliy from specificMort.txt @@ -217,15 +217,18 @@ make_atlantis_diagnostic_figures = function( itype <- itype + 1 mort <- specificmort %>% dplyr::filter(mort == atype) - temp.plot = atlantistools::plot_line(mort, col = 'agecl') - temp.plot = ggplot2::update_labels( + temp.plot <- atlantistools::plot_line(mort, col = 'agecl') + temp.plot <- ggplot2::update_labels( p = temp.plot, labels = c( list(x = 'Time (years)', y = 'Mortality'), list(colour = 'Ageclas') ) ) - temp.plot = add.title(temp.plot, paste0('Mortality at Age (', atype, ')')) + temp.plot <- add.title( + temp.plot, + paste0('Mortality at Age (', atype, ')') + ) plotMort[[itype]] <- temp.plot } @@ -247,16 +250,16 @@ make_atlantis_diagnostic_figures = function( dplyr::filter(code %in% allCodes) %>% dplyr::filter(agecl == iage) - pct = atlantistools::agg_perc(mortality, groups = c('time', 'species')) - temp.plot = atlantistools::plot_bar(pct, fill = 'mort', wrap = 'species') - temp.plot = ggplot2::update_labels( + pct <- atlantistools::agg_perc(mortality, groups = c('time', 'species')) + temp.plot <- atlantistools::plot_bar(pct, fill = 'mort', wrap = 'species') + temp.plot <- ggplot2::update_labels( temp.plot, labels = list(x = 'Time (years)', y = 'Rate (proportion)') ) + ggplot2::scale_y_continuous( labels = scales::label_number(accuracy = 0.01) ) - temp.plot = add.title( + temp.plot <- add.title( temp.plot, paste0( "Relative Mortality Rates for species with 10 age classes (Age ", @@ -280,16 +283,16 @@ make_atlantis_diagnostic_figures = function( ) %>% dplyr::filter(agecl == i2age) - pct = atlantistools::agg_perc(mortality, groups = c('time', 'species')) - temp.plot = atlantistools::plot_bar(pct, fill = 'mort', wrap = 'species') - temp.plot = ggplot2::update_labels( + pct <- atlantistools::agg_perc(mortality, groups = c('time', 'species')) + temp.plot <- atlantistools::plot_bar(pct, fill = 'mort', wrap = 'species') + temp.plot <- ggplot2::update_labels( temp.plot, labels = list(x = 'Time (years)', y = 'Rate (proportion)') ) + ggplot2::scale_y_continuous( labels = scales::label_number(accuracy = 0.01) ) - temp.plot = add.title( + temp.plot <- add.title( temp.plot, paste0( "Relative Mortality Rates for species with 2 age classes (Age ", @@ -312,16 +315,16 @@ make_atlantis_diagnostic_figures = function( ) %>% dplyr::filter(agecl == 1) - pct = atlantistools::agg_perc(mortality, groups = c('time', 'species')) - temp.plot = atlantistools::plot_bar(pct, fill = 'mort', wrap = 'species') - temp.plot = ggplot2::update_labels( + pct <- atlantistools::agg_perc(mortality, groups = c('time', 'species')) + temp.plot <- atlantistools::plot_bar(pct, fill = 'mort', wrap = 'species') + temp.plot <- ggplot2::update_labels( temp.plot, labels = list(x = 'Time (years)', y = 'Rate (proportion)') ) + ggplot2::scale_y_continuous( labels = scales::label_number(accuracy = 0.01) ) - temp.plot = add.title( + temp.plot <- add.title( temp.plot, paste0("Relative Mortality Rates for species with 1 age class (Age 1) ") ) @@ -350,27 +353,27 @@ make_atlantis_diagnostic_figures = function( #Make overall biomass plot (stacked barplot of total biomass domain-wide) if (plot.overall.biomass | plot.all) { print("biomass overall") - biomass = readRDS(file.path(out.dir, 'biomass.rds')) + biomass <- readRDS(file.path(out.dir, 'biomass.rds')) #combine threshold = 10 - biomass.df.10 = atlantistools::combine_groups( + biomass.df.10 <- atlantistools::combine_groups( biomass, group_col = 'species', combine_thresh = 10 ) - temp.plot.1 = atlantistools::plot_bar(biomass.df.10) - temp.plot.1 = temp.plot.1 + ggplot2::ggtitle('Top 10 Groups') - temp.plot.1 = ggplot2::update_labels(temp.plot.1, labels = plot.labels) + temp.plot.1 <- atlantistools::plot_bar(biomass.df.10) + temp.plot.1 <- temp.plot.1 + ggplot2::ggtitle('Top 10 Groups') + temp.plot.1 <- ggplot2::update_labels(temp.plot.1, labels = plot.labels) #Combine threshold = 20 - biomass.df.20 = atlantistools::combine_groups( + biomass.df.20 <- atlantistools::combine_groups( biomass, group_col = 'species', combine_thresh = 20 ) - temp.plot.2 = atlantistools::plot_bar(biomass.df.20) - temp.plot.2 = temp.plot.2 + ggplot2::ggtitle('Top 20 Groups') - temp.plot.2 = ggplot2::update_labels(temp.plot.2, labels = plot.labels) + temp.plot.2 <- atlantistools::plot_bar(biomass.df.20) + temp.plot.2 <- temp.plot.2 + ggplot2::ggtitle('Top 20 Groups') + temp.plot.2 <- ggplot2::update_labels(temp.plot.2, labels = plot.labels) pdf( file = file.path(fig.dir, paste0(run.name, ' overall biomass.pdf')), @@ -388,71 +391,71 @@ make_atlantis_diagnostic_figures = function( #Make biomass timeseries plots if (plot.biomass.timeseries | plot.all) { print("biomass ") - biomass = readRDS(file.path(out.dir, 'biomass.rds')) + biomass <- readRDS(file.path(out.dir, 'biomass.rds')) #biomass by species timeseries - temp.plot.1 = atlantistools::plot_line(biomass) - temp.plot.1 = ggplot2::update_labels(temp.plot.1, labels = plot.labels) - temp.plot.1 = add.title(temp.plot.1, 'Biomass') + temp.plot.1 <- atlantistools::plot_line(biomass) + temp.plot.1 <- ggplot2::update_labels(temp.plot.1, labels = plot.labels) + temp.plot.1 <- add.title(temp.plot.1, 'Biomass') #biomass at age timeseries - biomass.age = readRDS(file.path(out.dir, 'biomass_age.rds')) + biomass.age <- readRDS(file.path(out.dir, 'biomass_age.rds')) - temp.plot.2 = atlantistools::plot_line(biomass.age, col = 'agecl') - temp.plot.2 = ggplot2::update_labels( + temp.plot.2 <- atlantistools::plot_line(biomass.age, col = 'agecl') + temp.plot.2 <- ggplot2::update_labels( p = temp.plot.2, labels = c(plot.labels, list(colour = 'Ageclas')) ) - temp.plot.2 = add.title(temp.plot.2, 'Biomass at Age') + temp.plot.2 <- add.title(temp.plot.2, 'Biomass at Age') #biomass at age relative to initial biomass timeseries - rel.biomass.age = atlantistools::convert_relative_initial(biomass.age) - temp.plot.3 = atlantistools::plot_line(rel.biomass.age, col = 'agecl') - temp.plot.3 = ggplot2::update_labels( + rel.biomass.age <- atlantistools::convert_relative_initial(biomass.age) + temp.plot.3 <- atlantistools::plot_line(rel.biomass.age, col = 'agecl') + temp.plot.3 <- ggplot2::update_labels( temp.plot.3, list(x = 'Time (years)', y = expression(biomass / bio[init])) ) - temp.plot.3 = atlantistools::plot_add_box(temp.plot.3) - temp.plot.3 = add.title( + temp.plot.3 <- atlantistools::plot_add_box(temp.plot.3) + temp.plot.3 <- add.title( temp.plot.3, 'Biomass at Age Relatative to Initial Biomass' ) #Biomass vs Bio init - rel.biomass = atlantistools::convert_relative_initial(biomass) - temp.plot.4 = atlantistools::plot_line(rel.biomass) - temp.plot.4 = ggplot2::update_labels( + rel.biomass <- atlantistools::convert_relative_initial(biomass) + temp.plot.4 <- atlantistools::plot_line(rel.biomass) + temp.plot.4 <- ggplot2::update_labels( temp.plot.4, list(x = 'Time (years)', y = expression(Biomass / Biomass[init])) ) - temp.plot.4 = atlantistools::plot_add_box(temp.plot.4) - temp.plot.4 = add.title( + temp.plot.4 <- atlantistools::plot_add_box(temp.plot.4) + temp.plot.4 <- add.title( temp.plot.4, 'Biomass Relatative to Initial Biomass' ) #Invert bio timeseries - biomass.age.invert = readRDS(file.path(out.dir, 'biomass_age_invert.rds')) + biomass.age.invert <- readRDS(file.path(out.dir, 'biomass_age_invert.rds')) - temp.plot.5 = atlantistools::plot_line(biomass.age.invert) - temp.plot.5 = ggplot2::update_labels(temp.plot.5, labels = plot.labels) - temp.plot.5 = add.title(temp.plot.5, 'Invert Biomass') + temp.plot.5 <- atlantistools::plot_line(biomass.age.invert) + temp.plot.5 <- ggplot2::update_labels(temp.plot.5, labels = plot.labels) + temp.plot.5 <- add.title(temp.plot.5, 'Invert Biomass') #Bio at age - bio.age.pct = atlantistools::agg_perc( + bio.age.pct <- atlantistools::agg_perc( biomass.age, groups = c('time', 'species') ) - temp.plot.6 = atlantistools::plot_bar( + temp.plot.6 <- atlantistools::plot_bar( bio.age.pct, fill = 'agecl', wrap = 'species' ) - temp.plot.6 = ggplot2::update_labels( + temp.plot.6 <- ggplot2::update_labels( temp.plot.6, labels = list(x = 'Time (years)', y = 'Numbers (%)') ) - temp.plot.6 = add.title(temp.plot.6, 'Biomass at age - Percent') + temp.plot.6 <- add.title(temp.plot.6, 'Biomass at age - Percent') pdf( file.path(fig.dir, paste0(run.name, ' biomass timeseries.pdf')), @@ -476,19 +479,19 @@ make_atlantis_diagnostic_figures = function( #Make length.age plots if (plot.length.age | plot.all) { print("length age") - length.age = readRDS(file.path(out.dir, 'length_age.rds')) + length.age <- readRDS(file.path(out.dir, 'length_age.rds')) #Length at age ts by spp # init.length.old = read.csv(paste0(param.dir,'/vertebrate_init_length_cm.csv'),header =T, stringsAsFactors = F) - init.length = read.csv( + init.length <- read.csv( file.path(param.dir, 'vertebrate_init_length_cm_Adjusted.csv'), header = T, stringsAsFactors = F ) %>% dplyr::select(Code, species, agecl, new.length.ref) %>% tidyr::spread(agecl, new.length.ref) - init.length = init.length[order(init.length$species), ] - spp.names = unique(length.age$species) + init.length <- init.length[order(init.length$species), ] + spp.names <- unique(length.age$species) pdf( file = file.path( @@ -500,9 +503,9 @@ make_atlantis_diagnostic_figures = function( onefile = T ) for (x in 1:length(spp.names)) { - spp.id = spp.names[x] - length.age.spp = dplyr::filter(length.age, species == spp.id & time > 0) - init.length.age.spp = dplyr::filter( + spp.id <- spp.names[x] + length.age.spp <- dplyr::filter(length.age, species == spp.id & time > 0) + init.length.age.spp <- dplyr::filter( init.length[, 2:12], species == spp.id ) @@ -518,22 +521,22 @@ make_atlantis_diagnostic_figures = function( dev.off() #Length age together - temp.plot.1 = atlantistools::plot_line(length.age, col = 'agecl') - temp.plot.1 = ggplot2::update_labels( + temp.plot.1 <- atlantistools::plot_line(length.age, col = 'agecl') + temp.plot.1 <- ggplot2::update_labels( temp.plot.1, labels = c(x = 'Time (years)', y = 'Length (cm)', colour = 'Ageclass') ) - temp.plot.1 = add.title(temp.plot.1, 'Length-at-age') + temp.plot.1 <- add.title(temp.plot.1, 'Length-at-age') #Length at age vs. length init - rel.length.age = atlantistools::convert_relative_initial(length.age) - temp.plot.2 = atlantistools::plot_line(rel.length.age, col = 'agecl') - temp.plot.2 = ggplot2::update_labels( + rel.length.age <- atlantistools::convert_relative_initial(length.age) + temp.plot.2 <- atlantistools::plot_line(rel.length.age, col = 'agecl') + temp.plot.2 <- ggplot2::update_labels( temp.plot.2, list(x = 'Time (years)', y = expression(length / length[init])) ) atlantistools::plot_add_box(temp.plot.2) - temp.plot.2 = add.title(temp.plot.2, 'Length at Age vs. Initial Length') + temp.plot.2 <- add.title(temp.plot.2, 'Length at Age vs. Initial Length') pdf( file = file.path( @@ -601,27 +604,27 @@ make_atlantis_diagnostic_figures = function( #Make Biomass Box plots if (plot.biomass.box | plot.all) { print("biomass box") - biomass.box = readRDS(file.path(out.dir, 'biomass_box.rds')) + biomass.box <- readRDS(file.path(out.dir, 'biomass_box.rds')) #Bio per box - temp.plot.1 = atlantistools::plot_line(biomass.box) - temp.plot.1 = ggplot2::update_labels(temp.plot.1, plot.labels) - temp.plot.1 = atlantistools::custom_grid( + temp.plot.1 <- atlantistools::plot_line(biomass.box) + temp.plot.1 <- ggplot2::update_labels(temp.plot.1, plot.labels) + temp.plot.1 <- atlantistools::custom_grid( temp.plot.1, grid_x = 'polygon', grid_y = 'species' ) - temp.plot.1 = add.title(temp.plot.1, 'Biomass by Box') + temp.plot.1 <- add.title(temp.plot.1, 'Biomass by Box') #Invert bio by box - biomass.box.invert = readRDS(file.path(out.dir, 'biomass_box_invert.rds')) - temp.plot.2 = atlantistools::plot_line(biomass.box.invert) - temp.plot.2 = ggplot2::update_labels(temp.plot.2, plot.labels) - temp.plot.2 = atlantistools::custom_grid( + biomass.box.invert <- readRDS(file.path(out.dir, 'biomass_box_invert.rds')) + temp.plot.2 <- atlantistools::plot_line(biomass.box.invert) + temp.plot.2 <- ggplot2::update_labels(temp.plot.2, plot.labels) + temp.plot.2 <- atlantistools::custom_grid( temp.plot.2, grid_x = 'polygon', grid_y = 'species' ) - temp.plot.2 = add.title(temp.plot.2, 'Invert Biomass by Box') + temp.plot.2 <- add.title(temp.plot.2, 'Invert Biomass by Box') pdf( file = file.path( @@ -776,51 +779,51 @@ make_atlantis_diagnostic_figures = function( #SN/RN plots if (plot.sn.rn | plot.all) { print("sn,rn") - SN.box = readRDS(file.path(out.dir, 'SN_box.rds')) - RN.box = readRDS(file.path(out.dir, 'RN_box.rds')) - RN.age = readRDS(file.path(out.dir, 'RN_age.rds')) - SN.age = readRDS(file.path(out.dir, 'SN_age.rds')) + SN.box <- readRDS(file.path(out.dir, 'SN_box.rds')) + RN.box <- readRDS(file.path(out.dir, 'RN_box.rds')) + RN.age <- readRDS(file.path(out.dir, 'RN_age.rds')) + SN.age <- readRDS(file.path(out.dir, 'SN_age.rds')) #SN per box - temp.plot.1 = atlantistools::plot_line(SN.box) - temp.plot.1 = atlantistools::custom_grid( + temp.plot.1 <- atlantistools::plot_line(SN.box) + temp.plot.1 <- atlantistools::custom_grid( temp.plot.1, grid_x = 'polygon', grid_y = 'species' ) - temp.plot.1 = add.title(temp.plot.1, 'SN by Box') + temp.plot.1 <- add.title(temp.plot.1, 'SN by Box') #RN per box - temp.plot.2 = atlantistools::plot_line(RN.box) - temp.plot.2 = atlantistools::custom_grid( + temp.plot.2 <- atlantistools::plot_line(RN.box) + temp.plot.2 <- atlantistools::custom_grid( temp.plot.2, grid_x = 'polygon', grid_y = 'species' ) - temp.plot.2 = add.title(temp.plot.2, 'RN by Box') + temp.plot.2 <- add.title(temp.plot.2, 'RN by Box') #SN vs SN init - SN.rel = atlantistools::convert_relative_initial(SN.age) - temp.plot.3 = atlantistools::plot_line(SN.rel, col = 'agecl') - temp.plot.3 = ggplot2::update_labels( + SN.rel <- atlantistools::convert_relative_initial(SN.age) + temp.plot.3 <- atlantistools::plot_line(SN.rel, col = 'agecl') + temp.plot.3 <- ggplot2::update_labels( temp.plot.3, list(x = 'Time (years)', y = expression(SN / SN[init])) ) - temp.plot.3 = atlantistools::plot_add_box(temp.plot.3) - temp.plot.3 = add.title(temp.plot.3, 'SN vs SN Init') + temp.plot.3 <- atlantistools::plot_add_box(temp.plot.3) + temp.plot.3 <- add.title(temp.plot.3, 'SN vs SN Init') #RN vs RN init - RN.rel = atlantistools::convert_relative_initial(RN.age) - temp.plot.4 = atlantistools::plot_line(RN.rel, col = 'agecl') - temp.plot.4 = ggplot2::update_labels( + RN.rel <- atlantistools::convert_relative_initial(RN.age) + temp.plot.4 <- atlantistools::plot_line(RN.rel, col = 'agecl') + temp.plot.4 <- ggplot2::update_labels( temp.plot.4, list(x = 'Time (years)', y = expression(RN / RN[init])) ) - temp.plot.4 = atlantistools::plot_add_box(temp.plot.4) - temp.plot.4 = add.title(temp.plot.4, 'RN vs RN Init') + temp.plot.4 <- atlantistools::plot_add_box(temp.plot.4) + temp.plot.4 <- add.title(temp.plot.4, 'RN vs RN Init') #SN/RN domain-wide - RN.SN = SN.box %>% + RN.SN <- SN.box %>% dplyr::rename('SN' = atoutput) %>% dplyr::left_join(RN.box, by = c("species", "polygon", "time")) %>% dplyr::rename('RN' = atoutput) %>% @@ -828,7 +831,7 @@ make_atlantis_diagnostic_figures = function( dplyr::summarize(SN = sum(SN, na.rm = T), RN = sum(RN, na.rm = T)) %>% dplyr::mutate(RN.SN = RN / SN) - temp.plot.5 = ggplot2::ggplot(RN.SN, ggplot2::aes(x = time, y = RN.SN)) + + temp.plot.5 <- ggplot2::ggplot(RN.SN, ggplot2::aes(x = time, y = RN.SN)) + ggplot2::geom_line() + ggplot2::geom_hline(yintercept = 2.65, lty = 2) + ggplot2::facet_wrap(~species) + @@ -863,35 +866,35 @@ make_atlantis_diagnostic_figures = function( if (plot.recruits | plot.all) { print("recruits") # Recruits TS - ssb.recruits = readRDS(file.path(out.dir, 'ssb_recruits.rds')) + ssb.recruits <- readRDS(file.path(out.dir, 'ssb_recruits.rds')) - temp.plot.1 = atlantistools::plot_line(ssb.recruits, y = 'rec') - temp.plot.1 = ggplot2::update_labels( + temp.plot.1 <- atlantistools::plot_line(ssb.recruits, y = 'rec') + temp.plot.1 <- ggplot2::update_labels( temp.plot.1, labels = list(x = 'Time (days)', y = 'Numbers') ) - temp.plot.1 = add.title(temp.plot.1, 'Recruits') + temp.plot.1 <- add.title(temp.plot.1, 'Recruits') # SSB TS - temp.plot.2 = atlantistools::plot_line(ssb.recruits, y = 'ssb') - temp.plot.2 = ggplot2::update_labels( + temp.plot.2 <- atlantistools::plot_line(ssb.recruits, y = 'ssb') + temp.plot.2 <- ggplot2::update_labels( temp.plot.2, labels = list(x = 'Time (days)', y = 'Numbers') ) - temp.plot.2 = add.title(temp.plot.2, 'SSB') + temp.plot.2 <- add.title(temp.plot.2, 'SSB') # Recruit per SSB - ssb.recruits$rec.per.sbb = ssb.recruits$rec / ssb.recruits$ssb - temp.plot.3 = atlantistools::plot_line( + ssb.recruits$rec.per.sbb <- ssb.recruits$rec / ssb.recruits$ssb + temp.plot.3 <- atlantistools::plot_line( ssb.recruits, y = 'rec.per.sbb', yexpand = T ) - temp.plot.3 = ggplot2::update_labels( + temp.plot.3 <- ggplot2::update_labels( temp.plot.3, labels = list(x = 'Time (days)', y = 'Numbers') ) - temp.plot.3 = add.title(temp.plot.3, 'Recruits per SSB') + temp.plot.3 <- add.title(temp.plot.3, 'Recruits per SSB') pdf( file = file.path( @@ -916,60 +919,60 @@ make_atlantis_diagnostic_figures = function( if (plot.numbers.timeseries | plot.all) { print("numbers") #Numbers TS - numbers = readRDS(file.path(out.dir, 'numbers.rds')) + numbers <- readRDS(file.path(out.dir, 'numbers.rds')) - temp.plot.1 = atlantistools::plot_line(numbers) - temp.plot.1 = ggplot2::update_labels( + temp.plot.1 <- atlantistools::plot_line(numbers) + temp.plot.1 <- ggplot2::update_labels( temp.plot.1, labels = list(x = 'Time (years)', y = 'Numbers') ) - temp.plot.1 = add.title(temp.plot.1, 'Numbers') + temp.plot.1 <- add.title(temp.plot.1, 'Numbers') #Numbers at age - numbers.age = readRDS(file.path(out.dir, 'numbers_age.rds')) + numbers.age <- readRDS(file.path(out.dir, 'numbers_age.rds')) - temp.plot.2 = atlantistools::plot_line(numbers.age, col = 'agecl') - temp.plot.2 = ggplot2::update_labels( + temp.plot.2 <- atlantistools::plot_line(numbers.age, col = 'agecl') + temp.plot.2 <- ggplot2::update_labels( temp.plot.2, labels = list(x = 'Time (years)', y = 'Numbers', colour = 'Ageclass') ) - temp.plot.2 = add.title(temp.plot.2, 'Numbers at age') + temp.plot.2 <- add.title(temp.plot.2, 'Numbers at age') #Num age vs num init - nums.rel = atlantistools::convert_relative_initial(numbers.age) - temp.plot.3 = atlantistools::plot_line(nums.rel, col = 'agecl') - temp.plot.3 = ggplot2::update_labels( + nums.rel <- atlantistools::convert_relative_initial(numbers.age) + temp.plot.3 <- atlantistools::plot_line(nums.rel, col = 'agecl') + temp.plot.3 <- ggplot2::update_labels( temp.plot.3, list(x = 'Time (years)', y = expression(Numbers / Numbers[init])) ) - temp.plot.3 = atlantistools::plot_add_box(temp.plot.3) - temp.plot.3 = add.title(temp.plot.3, 'Numbers vs. Initial Numbers') + temp.plot.3 <- atlantistools::plot_add_box(temp.plot.3) + temp.plot.3 <- add.title(temp.plot.3, 'Numbers vs. Initial Numbers') #Numbers per ageclass, used to scale recruitment values from initial conditions ageclass 1 - nums.rel = dplyr::group_by(nums.rel, species, agecl) - nums.rel = dplyr::summarise(nums.rel, avg = mean(atoutput)) - nums.rel = tidyr::spread(nums.rel, agecl, avg) + nums.rel <- dplyr::group_by(nums.rel, species, agecl) + nums.rel <- dplyr::summarise(nums.rel, avg = mean(atoutput)) + nums.rel <- tidyr::spread(nums.rel, agecl, avg) - nums.scale = 1 / rowMeans(nums.rel[, 2]) - nums.c = data.frame(nums.rel[, 1]) - nums.c$scale = nums.scale + nums.scale <- 1 / rowMeans(nums.rel[, 2]) + nums.c <- data.frame(nums.rel[, 1]) + nums.c$scale <- nums.scale - nums.init = dplyr::filter(numbers.age, time == 0 & agecl == 1) + nums.init <- dplyr::filter(numbers.age, time == 0 & agecl == 1) - RN.age = readRDS(file.path(out.dir, 'RN_age.rds')) - SN.age = readRDS(file.path(out.dir, 'SN_age.rds')) + RN.age <- readRDS(file.path(out.dir, 'RN_age.rds')) + SN.age <- readRDS(file.path(out.dir, 'SN_age.rds')) - RN.age = dplyr::filter(RN.age, time == 0 & agecl == 1) - SN.age = dplyr::filter(SN.age, time == 0 & agecl == 1) + RN.age <- dplyr::filter(RN.age, time == 0 & agecl == 1) + SN.age <- dplyr::filter(SN.age, time == 0 & agecl == 1) - nums.RN.SN = dplyr::left_join( + nums.RN.SN <- dplyr::left_join( nums.init, group.index, by = c('species' = 'LongName') ) - nums.RN.SN$SN_RN = RN.age$atoutput + SN.age$atoutput - nums.RN.SN$totalN = nums.RN.SN$atoutput * nums.RN.SN$SN_RN - numscale.f = dplyr::left_join(nums.RN.SN, nums.c, by = 'species') + nums.RN.SN$SN_RN <- RN.age$atoutput + SN.age$atoutput + nums.RN.SN$totalN <- nums.RN.SN$atoutput * nums.RN.SN$SN_RN + numscale.f <- dplyr::left_join(nums.RN.SN, nums.c, by = 'species') write.csv( numscale.f, @@ -981,20 +984,20 @@ make_atlantis_diagnostic_figures = function( ) #num at age % - num.pct = atlantistools::agg_perc( + num.pct <- atlantistools::agg_perc( numbers.age, groups = c('time', 'species') ) - temp.plot.4 = atlantistools::plot_bar( + temp.plot.4 <- atlantistools::plot_bar( num.pct, fill = 'agecl', wrap = 'species' ) - temp.plot.4 = ggplot2::update_labels( + temp.plot.4 <- ggplot2::update_labels( temp.plot.4, labels = list(x = "Time (years0", y = 'Numbers (%)') ) - temp.plot.4 = add.title(temp.plot.4, 'Numbers at age - Percent') + temp.plot.4 <- add.title(temp.plot.4, 'Numbers at age - Percent') #Biomass Pool Grazers #Biomass pool grazing @@ -1026,88 +1029,88 @@ make_atlantis_diagnostic_figures = function( if (plot.physics | plot.all) { print("physics") #Physics snapshot - physics.statevars = readRDS(file.path(out.dir, 'physics_statevars.rds')) + physics.statevars <- readRDS(file.path(out.dir, 'physics_statevars.rds')) - temp.plot.1 = atlantistools::plot_line(physics.statevars, wrap = NULL) - temp.plot.1 = atlantistools::custom_grid( + temp.plot.1 <- atlantistools::plot_line(physics.statevars, wrap = NULL) + temp.plot.1 <- atlantistools::custom_grid( temp.plot.1, grid_x = 'polygon', grid_y = 'variable' ) - temp.plot.1 = ggplot2::update_labels(temp.plot.1, list(y = '')) - temp.plot.1 = add.title(temp.plot.1, 'Physics Snapshot') + temp.plot.1 <- ggplot2::update_labels(temp.plot.1, list(y = '')) + temp.plot.1 <- add.title(temp.plot.1, 'Physics Snapshot') #Phys plots - physics = atlantistools::flip_layers(physics.statevars) - physics = split(physics, physics$variable) - phys.plots = list() + physics <- atlantistools::flip_layers(physics.statevars) + physics <- split(physics, physics$variable) + phys.plots <- list() for (v in 1:length(physics)) { - phys.plots[[v]] = atlantistools::plot_line(physics[[v]], wrap = NULL) - phys.plots[[v]] = atlantistools::custom_grid( + phys.plots[[v]] <- atlantistools::plot_line(physics[[v]], wrap = NULL) + phys.plots[[v]] <- atlantistools::custom_grid( phys.plots[[v]], grid_x = 'polygon', grid_y = 'layer' ) - phys.plots[[v]] = add.title(phys.plots[[v]], names(physics)[v]) - phys.plots[[v]] = ggplot2::update_labels( + phys.plots[[v]] <- add.title(phys.plots[[v]], names(physics)[v]) + phys.plots[[v]] <- ggplot2::update_labels( phys.plots[[v]], labels = list(x = 'time', y = names(phys.plots)[v]) ) } #fluxes 1 - flux = readRDS(file.path(out.dir, 'flux.rds')) - temp.plot.2 = atlantistools::flip_layers(flux) - temp.plot.2 = atlantistools::plot_line( + flux <- readRDS(file.path(out.dir, 'flux.rds')) + temp.plot.2 <- atlantistools::flip_layers(flux) + temp.plot.2 <- atlantistools::plot_line( temp.plot.2, wrap = NULL, col = 'variable' ) - temp.plot.2 = atlantistools::custom_grid( + temp.plot.2 <- atlantistools::custom_grid( temp.plot.2, grid_x = 'polygon', grid_y = 'layer' ) - temp.plot.2 = ggplot2::update_labels(temp.plot.2, list(y = '')) - temp.plot.2 = add.title(temp.plot.2, 'Fluxes') + temp.plot.2 <- ggplot2::update_labels(temp.plot.2, list(y = '')) + temp.plot.2 <- add.title(temp.plot.2, 'Fluxes') #fluxes 2 - source.sink = readRDS(file.path(out.dir, 'source_sink.rds')) - temp.plot.3 = atlantistools::flip_layers(source.sink) - temp.plot.3 = atlantistools::plot_line( + source.sink <- readRDS(file.path(out.dir, 'source_sink.rds')) + temp.plot.3 <- atlantistools::flip_layers(source.sink) + temp.plot.3 <- atlantistools::plot_line( temp.plot.3, wrap = NULL, col = 'variable' ) - temp.plot.3 = atlantistools::custom_grid( + temp.plot.3 <- atlantistools::custom_grid( temp.plot.3, grid_x = 'polygon', grid_y = 'layer' ) - temp.plot.3 = ggplot2::update_labels(temp.plot.3, list(y = '')) - temp.plot.3 = add.title(temp.plot.3, 'Source-Sink') + temp.plot.3 <- ggplot2::update_labels(temp.plot.3, list(y = '')) + temp.plot.3 <- add.title(temp.plot.3, 'Source-Sink') #Changes in wc w/ rel dz - dz = readRDS(file.path(out.dir, 'dz.rds')) - nominal.dz = readRDS(file.path(out.dir, 'nominal_dz.rds')) + dz <- readRDS(file.path(out.dir, 'dz.rds')) + nominal.dz <- readRDS(file.path(out.dir, 'nominal_dz.rds')) - check.dz = dplyr::left_join(dz, nominal.dz, by = c('polygon', 'layer')) - check.dz = dplyr::mutate(check.dz, check.dz = atoutput.x / atoutput.y) - check.dz = dplyr::filter(check.dz, !is.na(check.dz)) + check.dz <- dplyr::left_join(dz, nominal.dz, by = c('polygon', 'layer')) + check.dz <- dplyr::mutate(check.dz, check.dz = atoutput.x / atoutput.y) + check.dz <- dplyr::filter(check.dz, !is.na(check.dz)) - temp.plot.4 = atlantistools::plot_line( + temp.plot.4 <- atlantistools::plot_line( check.dz, x = 'time', y = 'check.dz', wrap = 'polygon', col = 'layer' ) - temp.plot.4 = ggplot2::update_labels( + temp.plot.4 <- ggplot2::update_labels( temp.plot.4, list(x = 'Time (years)', y = expression(dz / nominal_dz)) ) - temp.plot.4 = add.title(temp.plot.4, 'Change in Water Column Height') - temp.plot.4 = ggplot2::update_labels(temp.plot.4, list(y = '')) + temp.plot.4 <- add.title(temp.plot.4, 'Change in Water Column Height') + temp.plot.4 <- ggplot2::update_labels(temp.plot.4, list(y = '')) # Additional Temp plots (boxplot by layer and polygon) tempD <- physics.statevars |> @@ -1155,50 +1158,50 @@ make_atlantis_diagnostic_figures = function( if (plot.growth.cons | plot.all) { print("growth") #Growth at ageclass v growth init - growth.age = readRDS(file.path(out.dir, 'growth_age.rds')) - growth.rel = atlantistools::convert_relative_initial(growth.age) - temp.plot.1 = atlantistools::plot_line(growth.rel, col = 'agecl') - temp.plot.1 = ggplot2::update_labels( + growth.age <- readRDS(file.path(out.dir, 'growth_age.rds')) + growth.rel <- atlantistools::convert_relative_initial(growth.age) + temp.plot.1 <- atlantistools::plot_line(growth.rel, col = 'agecl') + temp.plot.1 <- ggplot2::update_labels( temp.plot.1, list(x = 'Time (years)', y = expression(Growth / Growth[init])) ) - temp.plot.1 = atlantistools::plot_add_box(temp.plot.1) - temp.plot.1 = add.title(temp.plot.1, 'Growth at Age vs. Iniital Growth') + temp.plot.1 <- atlantistools::plot_add_box(temp.plot.1) + temp.plot.1 <- add.title(temp.plot.1, 'Growth at Age vs. Iniital Growth') #Grwoth vs growth init - growth.rel.init = readRDS(file.path(out.dir, 'growth_rel_init.rds')) - temp.plot.2 = atlantistools::plot_line( + growth.rel.init <- readRDS(file.path(out.dir, 'growth_rel_init.rds')) + temp.plot.2 <- atlantistools::plot_line( growth.rel.init, y = 'gr_rel', col = 'agecl' ) - temp.plot.2 = ggplot2::update_labels( + temp.plot.2 <- ggplot2::update_labels( temp.plot.2, list(y = expression((Growth - Growth[req] / Growth[req]))) ) - temp.plot.2 = add.title(temp.plot.2, 'Growth vs. Initial Growth') + temp.plot.2 <- add.title(temp.plot.2, 'Growth vs. Initial Growth') #Consumptions at age vs. initial - eat_age = readRDS(file.path(out.dir, 'eat_age.rds')) - cons.rel = atlantistools::convert_relative_initial(eat_age) - temp.plot.3 = atlantistools::plot_line(cons.rel, col = 'agecl') - temp.plot.3 = ggplot2::update_labels( + eat_age <- readRDS(file.path(out.dir, 'eat_age.rds')) + cons.rel <- atlantistools::convert_relative_initial(eat_age) + temp.plot.3 <- atlantistools::plot_line(cons.rel, col = 'agecl') + temp.plot.3 <- ggplot2::update_labels( temp.plot.3, list(x = 'Time (years)', y = expression(Consumption / Consumption[init])) ) - temp.plot.3 = atlantistools::plot_add_box(temp.plot.3) - temp.plot.3 = add.title( + temp.plot.3 <- atlantistools::plot_add_box(temp.plot.3) + temp.plot.3 <- add.title( temp.plot.3, 'Consumption at Age vs. Initial Consumption' ) #Consumption at age timeseries - temp.plot.4 = atlantistools::plot_line(eat_age, col = 'agecl') - temp.plot.4 = ggplot2::update_labels( + temp.plot.4 <- atlantistools::plot_line(eat_age, col = 'agecl') + temp.plot.4 <- ggplot2::update_labels( temp.plot.4, list(x = 'Time (years)', y = 'Biomass (tonnes)', color = 'Ageclass') ) - temp.plot.4 = add.title(temp.plot.4, 'Consumption at Age') + temp.plot.4 <- add.title(temp.plot.4, 'Consumption at Age') pdf( file = file.path( @@ -1223,7 +1226,7 @@ make_atlantis_diagnostic_figures = function( #plot cohort timeseries if (plot.cohort | plot.all) { print("cohort") - numbers.age = readRDS(file.path(out.dir, 'numbers_age.rds')) + numbers.age <- readRDS(file.path(out.dir, 'numbers_age.rds')) pdf( file = file.path(fig.dir, paste0(run.name, ' Cohort Timeseries.pdf')), width = 24, @@ -1231,13 +1234,13 @@ make_atlantis_diagnostic_figures = function( onefile = T ) for (i in 1:10) { - age.sub = dplyr::filter(numbers.age, agecl == i) - temp.plot = atlantistools::plot_line(age.sub) - temp.plot = ggplot2::update_labels( + age.sub <- dplyr::filter(numbers.age, agecl == i) + temp.plot <- atlantistools::plot_line(age.sub) + temp.plot <- ggplot2::update_labels( temp.plot, list(x = 'Time (years)', y = 'Numbers') ) - temp.plot = add.title(temp.plot, paste0('Age-', i)) + temp.plot <- add.title(temp.plot, paste0('Age-', i)) gridExtra::grid.arrange(temp.plot) } dev.off() @@ -1250,12 +1253,12 @@ make_atlantis_diagnostic_figures = function( #Diet figures if (plot.diet | plot.all) { print("diet") - bio_consumed = readRDS(file.path(out.dir, 'biomass_consumed.rds')) + bio_consumed <- readRDS(file.path(out.dir, 'biomass_consumed.rds')) if (nrow(bio_consumed) > 0) { # diet.plots = atlantistools::plot_diet(result$biomass.consumed, wrap_col = 'agecl', combine_thresh = 3) - wrap_col = 'agecl' - combine_thresh = 3 - species = NULL + wrap_col <- 'agecl' + combine_thresh <- 3 + species <- NULL print("1") atlantistools::check_df_names( data = bio_consumed, @@ -1349,11 +1352,11 @@ make_atlantis_diagnostic_figures = function( print('Output file size too large to generate consumption plots') } else { #source(here::here('R','plot_overall_predation.R')) - consumption = get_consumption( + consumption <- get_consumption( prod.file = param.ls$prod.nc, fgs.file = param.ls$groups.file ) - data.sub = subset_diet( + data.sub <- subset_diet( diet.file = param.ls$dietcheck, consumption = consumption, spp.names = group.index$Code @@ -1378,9 +1381,9 @@ make_atlantis_diagnostic_figures = function( out.dir, 'biomass_spatial_stanza.rds' )) - volume = readRDS(file.path(out.dir, 'volume.rds')) + volume <- readRDS(file.path(out.dir, 'volume.rds')) - temp.plots = atlantistools::plot_spatial_box( + temp.plots <- atlantistools::plot_spatial_box( bio_spatial = biomass.spatial.stanza, bgm_as_df = atlantistools::convert_bgm(bgm = param.ls$bgm), timesteps = 7 @@ -1403,7 +1406,7 @@ make_atlantis_diagnostic_figures = function( biomass.spatial.stanza, !is.na(layer) ) - temp.plots.2 = atlantistools::plot_spatial_ts( + temp.plots.2 <- atlantistools::plot_spatial_ts( bio_spatial = biomass.spatial.stanza, bgm_as_df = atlantistools::convert_bgm(bgm = param.ls$bgm), vol = volume @@ -1440,31 +1443,31 @@ make_atlantis_diagnostic_figures = function( if (plot.spatial.catch | plot.all) { print("spatial catch") - bgm = atlantistools::convert_bgm(bgm = param.ls$bgm) + bgm <- atlantistools::convert_bgm(bgm = param.ls$bgm) - biomass.box = readRDS(file.path(out.dir, 'biomass_box.rds')) %>% + biomass.box <- readRDS(file.path(out.dir, 'biomass_box.rds')) %>% dplyr::filter(time >= (max(time) - 10)) %>% dplyr::group_by(species, polygon) %>% dplyr::summarise(biomass = mean(atoutput, na.rm = T)) - catch = readRDS(file.path(out.dir, 'catch.rds')) %>% + catch <- readRDS(file.path(out.dir, 'catch.rds')) %>% dplyr::filter(time >= (max(time) - 10)) %>% dplyr::group_by(species, polygon) %>% dplyr::summarise(catch = mean(atoutput, na.rm = T)) - biomass.catch.box = biomass.box %>% + biomass.catch.box <- biomass.box %>% dplyr::left_join(catch) - i = 1 + i <- 1 pdf(paste0(fig.dir, '/spatial_biomass_catch.pdf')) for (i in 1:nrow(group.index)) { - biomass.catch.spp = biomass.catch.box %>% + biomass.catch.spp <- biomass.catch.box %>% dplyr::filter(species == group.index$LongName[i]) - biomass.catch.spp.polygon = bgm %>% + biomass.catch.spp.polygon <- bgm %>% dplyr::left_join(biomass.catch.spp) - p1 = ggplot2::ggplot( + p1 <- ggplot2::ggplot( biomass.catch.spp.polygon, ggplot2::aes(x = long, y = lat, fill = biomass, group = polygon) ) + @@ -1478,7 +1481,7 @@ make_atlantis_diagnostic_figures = function( legend.position = 'bottom', legend.key.width = ggplot2::unit(0.4, 'in') ) - p2 = ggplot2::ggplot( + p2 <- ggplot2::ggplot( biomass.catch.spp.polygon, ggplot2::aes(x = long, y = lat, fill = catch, group = polygon) ) + diff --git a/R/plot_biomass_box_summary.R b/R/plot_biomass_box_summary.R index 0af4963..0877e1e 100644 --- a/R/plot_biomass_box_summary.R +++ b/R/plot_biomass_box_summary.R @@ -6,7 +6,7 @@ # Box Gradient ------------------------------------------------------------ # This needs to be reworked -plot_biomass_box_season = function( +plot_biomass_box_season <- function( bio.box, bio.box.invert, plot.presence = F, @@ -15,15 +15,15 @@ plot_biomass_box_season = function( fig.dir, tolerance = 0.1 ) { - box2epu = read.csv(here::here('Geometry', 'box2epu.csv')) - month2season = read.csv(here::here('data-raw', 'month2season.csv')) - season2name = data.frame( + 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')) + fgs <- read.csv(here::here('currentVersion', 'neus_groups.csv')) - move.param = read.csv(here::here('data-raw', 'seasonal_movements.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) %>% @@ -35,10 +35,10 @@ plot_biomass_box_season = function( mutate(orig.presence = ifelse(orig == 0, 0, 1)) %>% filter(!is.na(epu)) - bio.all = rbind(bio.box, bio.box.invert) %>% + bio.all <- rbind(bio.box, bio.box.invert) %>% mutate(day = time * 365) - bio.all = bio.all %>% + bio.all <- bio.all %>% mutate( date = as.POSIXct( day * 86400, @@ -52,25 +52,25 @@ plot_biomass_box_season = function( summarise(atoutput = mean(atoutput, na.rm = T)) if (!is.null(species.list)) { - bio.all = bio.all %>% + bio.all <- bio.all %>% filter(species %in% species.list) } - bio.all = bio.all %>% + bio.all <- bio.all %>% group_by(species, polygon, season.name) %>% summarise(atoutput = mean(atoutput, na.rm = T)) - bio.max = bio.all %>% + bio.max <- bio.all %>% group_by(species, season.name) %>% summarise(atoutput.max = sum(atoutput, na.rm = T)) - box.combs = expand.grid( + 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 %>% + 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)) %>% @@ -91,15 +91,15 @@ plot_biomass_box_season = function( ) ) - bio.presence$polygon = factor( + 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( + 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) ) + @@ -115,8 +115,8 @@ plot_biomass_box_season = function( plot.title = element_text(hjust = 0.5) ) } else { - plot.name = paste0(fig.dir, 'Box_EPU_Season_scaled.png') - p = ggplot( + 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) ) + @@ -156,7 +156,7 @@ plot_biomass_box_season = function( #' @noRd # Box Time Range ------------------------------------------------------------ -plot_biomass_box_range = function( +plot_biomass_box_range <- function( bio.box, bio.box.invert, day.min = NA, @@ -166,33 +166,33 @@ plot_biomass_box_range = function( save.fig = T, fig.dir ) { - box2epu = read.csv(here::here('Geometry', 'box2epu.csv')) + box2epu <- read.csv(here::here('Geometry', 'box2epu.csv')) - bio.all = rbind(bio.box, bio.box.invert) %>% + bio.all <- rbind(bio.box, bio.box.invert) %>% mutate(day = time * 365) - bio.all = bio.all %>% + bio.all <- bio.all %>% filter(day >= day.min & day <= day.max) if (!is.null(species.list)) { - bio.all = bio.all %>% + bio.all <- bio.all %>% filter(species %in% species.list) } - bio.all = bio.all %>% + bio.all <- bio.all %>% group_by(species, polygon) %>% summarise(atoutput = mean(atoutput, na.rm = T)) - bio.max = bio.all %>% + bio.max <- bio.all %>% group_by(species) %>% summarise(atoutput.max = max(atoutput, na.rm = T)) - box.combs = expand.grid( + box.combs <- expand.grid( species = unique(bio.all$species), polygon = sort(unique(bio.all$polygon)) ) - bio.presence = box.combs %>% + 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)) %>% @@ -202,23 +202,23 @@ plot_biomass_box_range = function( atoutput.scaled = ifelse(is.na(atoutput.scaled), 0, atoutput.scaled) ) - bio.presence$polygon = factor( + bio.presence$polygon <- factor( bio.presence$polygon, levels = c(1:7, 9, 8, 12:15, 10:11, 16:22) ) - spp.order = bio.presence %>% + spp.order <- bio.presence %>% group_by(species) %>% summarise(presence.freq = mean(atoutput.scaled)) %>% arrange(desc(presence.freq)) - bio.presence$species = factor( + bio.presence$species <- factor( bio.presence$species, levels = rev(spp.order$species) ) if (plot.presence) { - p = ggplot( + p <- ggplot( bio.presence, aes(x = polygon, y = species, alpha = presence, fill = epu) ) + @@ -233,9 +233,9 @@ plot_biomass_box_range = function( theme_minimal() + theme(panel.grid = element_blank(), legend.position = 'bottom') - plot.name = paste0(fig.dir, 'Box_EPU_summary.png') + plot.name <- paste0(fig.dir, 'Box_EPU_summary.png') } else { - p = ggplot( + p <- ggplot( bio.presence, aes(x = polygon, y = species, alpha = atoutput.scaled, fill = epu) ) + @@ -250,7 +250,7 @@ plot_biomass_box_range = function( theme_minimal() + theme(panel.grid = element_blank(), legend.position = 'bottom') - plot.name = paste0(fig.dir, 'Box_EPU_summary.png') + plot.name <- paste0(fig.dir, 'Box_EPU_summary.png') } if (save.fig) { diff --git a/R/plot_overall_predation.R b/R/plot_overall_predation.R index c9b1540..4d7da72 100644 --- a/R/plot_overall_predation.R +++ b/R/plot_overall_predation.R @@ -4,24 +4,24 @@ # Gets total consumption for all groups ----------------------------------- get_consumption <- function(prod.file, fgs.file) { - prod.nc = ncdf4::nc_open(prod.file) - prod.vars = names(prod.nc$var) - groups = read.csv(fgs.file, as.is = T) - time.vals = as.Date(as.POSIXct( + prod.nc <- ncdf4::nc_open(prod.file) + prod.vars <- names(prod.nc$var) + groups <- read.csv(fgs.file, as.is = T) + time.vals <- as.Date(as.POSIXct( prod.nc$dim$t$vals, origin = '1964-01-01 00:00:00', tz = 'UTC' )) - time.days = as.numeric(difftime( + time.days <- as.numeric(difftime( time.vals, as.Date('1964-01-01'), units = 'd' )) - consumption.all.ls = list() + consumption.all.ls <- list() for (i in 1:nrow(groups)) { - prod.vars.group = prod.vars[grep(groups$Name[i], prod.vars)] - group.eat = prod.vars.group[c( + prod.vars.group <- prod.vars[grep(groups$Name[i], prod.vars)] + group.eat <- prod.vars.group[c( grep('Eat', prod.vars.group), grep('Grazing', prod.vars.group) )] @@ -29,16 +29,16 @@ get_consumption <- function(prod.file, fgs.file) { if (length(group.eat) == 0) { next() } - group.df.ls = list() + group.df.ls <- list() for (v in 1:length(group.eat)) { - age.var = colSums(ncdf4::ncvar_get(prod.nc, group.eat[v])) + age.var <- colSums(ncdf4::ncvar_get(prod.nc, group.eat[v])) if (grepl('_Eat', group.eat[v])) { - group.cohort = as.numeric(strsplit( + group.cohort <- as.numeric(strsplit( group.eat[v], paste0(groups$Name[i], '|_Eat') )[[1]][2]) - 1 - group.df.ls[[v]] = data.frame( + group.df.ls[[v]] <- data.frame( Predator = groups$Code[i], Cohort = group.cohort, Time = time.days, @@ -47,7 +47,7 @@ get_consumption <- function(prod.file, fgs.file) { stringsAsFactors = F ) } else { - group.df.ls[[v]] = data.frame( + group.df.ls[[v]] <- data.frame( Predator = groups$Code[i], Cohort = 0, Time = time.days, @@ -57,9 +57,9 @@ get_consumption <- function(prod.file, fgs.file) { ) } } - consumption.all.ls[[i]] = do.call('rbind', group.df.ls) + consumption.all.ls[[i]] <- do.call('rbind', group.df.ls) } - consumption.all = do.call('rbind', consumption.all.ls) + consumption.all <- do.call('rbind', consumption.all.ls) return(consumption.all) } @@ -68,10 +68,10 @@ get_consumption <- function(prod.file, fgs.file) { #Reads in diet data, transforms to longform, adds total consumption, and subsets subset_diet <- function(diet.file, consumption, spp.names) { - data = data.table::fread(diet.file) %>% + data <- data.table::fread(diet.file) %>% dplyr::select(-Stock, -Updated) #Convert Data to long format - data.long = reshape2::melt( + data.long <- reshape2::melt( data, id.vars = c('Time', 'Predator', 'Cohort'), variable.name = 'Prey', @@ -102,11 +102,11 @@ plot_overall_predation <- function( file.prefix ) { #Collapse small contributors into "Rest" - data.tot = data %>% + data.tot <- data %>% dplyr::group_by(Time, date, Prey) %>% dplyr::summarize(consumed.tot = sum(consumed.prey, na.rm = T)) - data.new = data %>% + data.new <- data %>% dplyr::left_join(data.tot) %>% # dplyr::group_by(Time,date,Prey,Predator) %>% dplyr::mutate( @@ -114,14 +114,14 @@ plot_overall_predation <- function( less.min = consumed.pct < min.fract ) - data.small.pct = data.new %>% + data.small.pct <- data.new %>% dplyr::filter(less.min == T) %>% dplyr::group_by(Time, date, Prey) %>% dplyr::summarize(consumed.pct = sum(consumed.pct)) %>% dplyr::mutate(Predator = 'Rest') %>% dplyr::arrange(Time, Prey, Predator, consumed.pct) - data.final = data.new %>% + data.final <- data.new %>% dplyr::filter(less.min == F) %>% dplyr::select(-less.min) %>% dplyr::bind_rows(data.small.pct) %>% @@ -134,87 +134,87 @@ plot_overall_predation <- function( ) #Get Biomass Data - biomass.data = read.table(bioindex.file, header = T, stringsAsFactors = F) - biomass.colnames = colnames(biomass.data) + biomass.data <- read.table(bioindex.file, header = T, stringsAsFactors = F) + biomass.colnames <- colnames(biomass.data) #Get catch Data if (file.exists(catch.file)) { - catch.data = read.table(catch.file, header = T, stringsAsFactors = F) - catch.colnames = colnames(catch.data) + catch.data <- read.table(catch.file, header = T, stringsAsFactors = F) + catch.colnames <- colnames(catch.data) } else { - catch.data = NULL - catch.colnames = NULL + catch.data <- NULL + catch.colnames <- NULL } #Loop through species - plot.cols = c( + plot.cols <- c( RColorBrewer::brewer.pal(12, 'Set3'), RColorBrewer::brewer.pal(8, 'Dark2'), RColorBrewer::brewer.pal(8, 'Set2'), RColorBrewer::brewer.pal(9, 'Set1') ) - plot.spp = sort(unique(as.character(data.new$Prey))) + plot.spp <- sort(unique(as.character(data.new$Prey))) - filename = file.path(fig.dir, paste0(file.prefix, '_TotalConsumption.pdf')) + filename <- file.path(fig.dir, paste0(file.prefix, '_TotalConsumption.pdf')) pdf(file = filename, width = 16, height = 8, onefile = T) for (i in 1:length(plot.spp)) { - data.spp = dplyr::filter(data.final, Prey == plot.spp[i]) + data.spp <- dplyr::filter(data.final, Prey == plot.spp[i]) #Identify all groups who have zero consumption values across all times/box/layers - which.zero = data.spp %>% + which.zero <- data.spp %>% dplyr::group_by(Predator) %>% dplyr::summarize(tot = sum(consumed.pct, na.rm = T)) %>% dplyr::mutate(all.zero = ifelse(tot == 0, T, F)) %>% dplyr::filter(all.zero == T) - which.zero = as.character(which.zero$Predator) + which.zero <- as.character(which.zero$Predator) #Remove zero consumption spp - data.spp = data.spp %>% + data.spp <- data.spp %>% dplyr::filter(!(Predator %in% which.zero)) #Total consumption - data.tot.spp = data.tot %>% + data.tot.spp <- data.tot %>% dplyr::filter(Prey == plot.spp[i]) #Biomass Timeseries - biomass.spp = biomass.data[, c( + biomass.spp <- biomass.data[, c( 1, grep(paste0('\\b', plot.spp[i], '\\b'), biomass.colnames) )] - colnames(biomass.spp)[2] = 'value' - biomass.spp$Metric = 'biomass' - biomass.spp$date = as.POSIXct( + colnames(biomass.spp)[2] <- 'value' + biomass.spp$Metric <- 'biomass' + biomass.spp$date <- as.POSIXct( biomass.spp$Time * 86400, origin = '1964-01-01 00:00:00', tz = 'UTC' ) #Catch Timeseries - catch.match = grep(paste0('\\b', plot.spp[i], '\\b'), catch.colnames) + catch.match <- grep(paste0('\\b', plot.spp[i], '\\b'), catch.colnames) if (length(catch.match) == 0) { - dum.dat = biomass.spp - dum.dat$value = NA - dum.dat$Metric = 'catch' - bio.catch.data = rbind(biomass.spp, dum.dat) + dum.dat <- biomass.spp + dum.dat$value <- NA + dum.dat$Metric <- 'catch' + bio.catch.data <- rbind(biomass.spp, dum.dat) } else { - catch.spp = catch.data[, c(1, catch.match)] - colnames(catch.spp)[2] = 'value' - catch.spp$Metric = 'catch' - catch.spp$date = as.POSIXct( + catch.spp <- catch.data[, c(1, catch.match)] + colnames(catch.spp)[2] <- 'value' + catch.spp$Metric <- 'catch' + catch.spp$date <- as.POSIXct( catch.spp$Time * 86400, origin = '1964-01-01 00:00:00', tz = 'UTC' ) #Combine Bio Catch - bio.catch.data = rbind(biomass.spp, catch.spp) + bio.catch.data <- rbind(biomass.spp, catch.spp) } #Plot prey - f1 = ggplot2::ggplot( + f1 <- ggplot2::ggplot( data.spp, ggplot2::aes(x = date, y = consumed.pct, fill = Predator) ) + @@ -228,7 +228,7 @@ plot_overall_predation <- function( ggplot2::guides(fill = ggplot2::guide_legend(nrow = 1)) + ggplot2::theme(legend.position = 'bottom') - f2 = ggplot2::ggplot( + f2 <- ggplot2::ggplot( data.tot.spp, ggplot2::aes(x = date, y = consumed.tot), size = 1.5 @@ -238,7 +238,7 @@ plot_overall_predation <- function( ggplot2::ylab('Total consumption (mg N m-3 d-1)') + ggplot2::theme_classic() - f3 = ggplot2::ggplot( + f3 <- ggplot2::ggplot( bio.catch.data, ggplot2::aes(x = date, y = value, lty = Metric) ) + diff --git a/R/process_atl_output.R b/R/process_atl_output.R index 283f38f..3c634f7 100644 --- a/R/process_atl_output.R +++ b/R/process_atl_output.R @@ -22,7 +22,7 @@ #' #' @export -process_atl_output = function( +process_atl_output <- function( param.dir, atl.dir, out.dir = file.path(atl.dir, 'Post_Processed/Data/'), @@ -67,52 +67,52 @@ process_atl_output = function( } #Read in groups file - fgs = atlantistools::load_fgs(param.ls$groups.file) %>% + fgs <- atlantistools::load_fgs(param.ls$groups.file) %>% dplyr::select(Code, Name, LongName) #Get boundary box - bboxes = atlantistools::get_boundary( + bboxes <- atlantistools::get_boundary( boxinfo = atlantistools::load_box(param.ls$bgm.file) ) #Get epibenthic biopool groups - bio.pools = atlantistools::load_bps(param.ls$groups.file, param.ls$init.file) + bio.pools <- atlantistools::load_bps(param.ls$groups.file, param.ls$init.file) #Get biomass conversion scalar - bio.conv = atlantistools::get_conv_mgnbiot(param.ls$biol.prm) + bio.conv <- atlantistools::get_conv_mgnbiot(param.ls$biol.prm) #All groups extracted (names, age-structured, biopools, and codes) - group.names = atlantistools::get_groups(param.ls$groups.file) - groups.age = atlantistools::get_age_groups(param.ls$groups.file) - groups.bp = group.names[!group.names %in% groups.age] - codes.age = atlantistools::get_age_acronyms(param.ls$groups.file) - groups.data = atlantistools::load_fgs(param.ls$groups.file) + group.names <- atlantistools::get_groups(param.ls$groups.file) + groups.age <- atlantistools::get_age_groups(param.ls$groups.file) + groups.bp <- group.names[!group.names %in% groups.age] + codes.age <- atlantistools::get_age_acronyms(param.ls$groups.file) + groups.data <- atlantistools::load_fgs(param.ls$groups.file) # Read Physics ------------------------------------------------------------ #Always make volume objects - vol.dz = atlantistools::load_nc_physics( + vol.dz <- atlantistools::load_nc_physics( nc = param.ls$main.nc, select_physics = c('volume', 'dz'), prm_run = param.ls$run.prm, bboxes = bboxes ) - dz = dplyr::filter(vol.dz, variable == 'dz') - vol = dplyr::filter(vol.dz, variable == 'volume') + dz <- dplyr::filter(vol.dz, variable == 'dz') + vol <- dplyr::filter(vol.dz, variable == 'volume') #Aggregate volume vertically - vol.ts = atlantistools::agg_data( + vol.ts <- atlantistools::agg_data( vol, groups = c('time', 'polygon'), fun = sum, out = 'volume' ) - nominal.dz = as.data.frame(atlantistools::load_init( + nominal.dz <- as.data.frame(atlantistools::load_init( init = param.ls$init.file, vars = 'nominal_dz' )) - nominal.dz = dplyr::filter(nominal.dz, !is.na(layer)) + nominal.dz <- dplyr::filter(nominal.dz, !is.na(layer)) saveRDS(vol.ts, file = file.path(out.dir, 'volume.rds')) saveRDS(dz, file = file.path(out.dir, 'dz.rds')) @@ -120,19 +120,19 @@ process_atl_output = function( rm(vol.ts) if (plot.physics | plot.all | process.all) { - flux = atlantistools::load_nc_physics( + flux <- atlantistools::load_nc_physics( nc = param.ls$main.nc, select_physics = c('eflux', 'vflux'), prm_run = param.ls$run.prm, bboxes = bboxes ) - source.sink = atlantistools::load_nc_physics( + source.sink <- atlantistools::load_nc_physics( nc = param.ls$main.nc, select_physics = c('hdsource', 'hdsink'), prm_run = param.ls$run.prm, bboxes = bboxes ) - phys.statevars = atlantistools::load_nc_physics( + phys.statevars <- atlantistools::load_nc_physics( nc = param.ls$main.nc, select_physics = c( 'salt', @@ -148,7 +148,7 @@ process_atl_output = function( ) #Exclude sediment from salinity - phys.statevars = dplyr::filter( + phys.statevars <- dplyr::filter( phys.statevars, !(variable == 'salt' & layer == max(layer) & time == min(time)) ) @@ -163,7 +163,7 @@ process_atl_output = function( # Other Parameter Objects ------------------------------------------------- #Read in age matrix - data.age.mat = atlantistools::prm_to_df( + data.age.mat <- atlantistools::prm_to_df( prm_biol = param.ls$biol.prm, fgs = param.ls$groups.file, group = codes.age, @@ -172,7 +172,7 @@ process_atl_output = function( saveRDS(data.age.mat, file = file.path(out.dir, 'data_age_mat.rds')) #Read in diet matrix - data.diet.mat = atlantistools::load_dietmatrix( + data.diet.mat <- atlantistools::load_dietmatrix( prm_biol = param.ls$biol.prm, fgs = param.ls$groups.file, convert_names = T @@ -180,53 +180,56 @@ process_atl_output = function( # saveRDS(data.diet.mat,file = file.path(out.dir,'diet_matrix.rds')) #length.age tempmat - biol.prm.lines = read.table( + biol.prm.lines <- read.table( param.ls$biol.prm, col.name = 1:100, comment.char = '', fill = T, header = F ) - lia.match = biol.prm.lines[grep('li_a_', biol.prm.lines[, 1]), 1:20] - tempmat = matrix(NA, nrow = nrow(lia.match), ncol = 3) + lia.match <- biol.prm.lines[grep('li_a_', biol.prm.lines[, 1]), 1:20] + tempmat <- matrix(NA, nrow = nrow(lia.match), ncol = 3) for (igroup in 1:nrow(tempmat)) { - tempmat[igroup, 1] = strsplit(as.character(lia.match[igroup, 1]), 'li_a_')[[ + tempmat[igroup, 1] <- strsplit( + as.character(lia.match[igroup, 1]), + 'li_a_' + )[[ 1 ]][2] } - tempmat[, 2] = as.numeric(as.character(lia.match[, 2])) - lib.match = grep('li_b_', biol.prm.lines[, 1]) - tempmat[, 3] = as.numeric(as.character(biol.prm.lines[lib.match, 2])) + tempmat[, 2] <- as.numeric(as.character(lia.match[, 2])) + lib.match <- grep('li_b_', biol.prm.lines[, 1]) + tempmat[, 3] <- as.numeric(as.character(biol.prm.lines[lib.match, 2])) - groups.data2 = groups.data[, c('Code', 'LongName')] - tempmat2 = as.data.frame(tempmat[2:dim(tempmat)[1], ]) - colnames(tempmat2) = c('Code', 'li_a', 'li_b') - tempmat3 = dplyr::left_join(tempmat2, groups.data2, by = 'Code') + groups.data2 <- groups.data[, c('Code', 'LongName')] + tempmat2 <- as.data.frame(tempmat[2:dim(tempmat)[1], ]) + colnames(tempmat2) <- c('Code', 'li_a', 'li_b') + tempmat3 <- dplyr::left_join(tempmat2, groups.data2, by = 'Code') ##Growth relative to initial conditions - recruit.weight = atlantistools::prm_to_df( + recruit.weight <- atlantistools::prm_to_df( prm_biol = param.ls$biol.prm, fgs = param.ls$groups.file, group = codes.age, parameter = c('KWRR', 'KWSR', 'AgeClassSize') ) - pd = atlantistools::load_init_weight( + pd <- atlantistools::load_init_weight( init = param.ls$init.nofill, fgs = param.ls$groups.file, bboxes = bboxes ) - pd = dplyr::left_join(pd, recruit.weight, by = "species") - pd = split(pd, pd$species) + pd <- dplyr::left_join(pd, recruit.weight, by = "species") + pd <- split(pd, pd$species) #Calculate weight difference from one ageclass to the next for (i in seq_along(pd)) { - pd[[i]]$wdiff = c( + pd[[i]]$wdiff <- c( (pd[[i]]$rn[1] + pd[[i]]$sn[1]) - (pd[[i]]$kwrr[1] + pd[[i]]$kwsr[1]), diff(pd[[i]]$rn + pd[[i]]$sn) ) } - pd = do.call(rbind, pd) - pd$growth_req = pd$wdiff / (365 * pd$ageclasssize) + pd <- do.call(rbind, pd) + pd$growth_req <- pd$wdiff / (365 * pd$ageclasssize) if (any(pd$growth_req < 0)) { warning( "Required growth negative for some groups. Please check your initial conditions files." @@ -235,24 +238,24 @@ process_atl_output = function( } unique(pd$species[which(pd$growth_req < 0)]) - growth.required = dplyr::select(pd, c(species, agecl, growth_req)) + growth.required <- dplyr::select(pd, c(species, agecl, growth_req)) # Process DietCheck ------------------------------------------------------- if (plot.diet | plot.all | process.all) { if (large.file == F) { - data.dietcheck.orig = atlantistools::load_dietcheck( + data.dietcheck.orig <- atlantistools::load_dietcheck( dietcheck = param.ls$dietcheck, fgs = param.ls$groups.file, prm_run = param.ls$run.prm, convert_names = T ) #Normalize proprotions so they always sum to 1 - dietcheck.tot = data.dietcheck.orig %>% + dietcheck.tot <- data.dietcheck.orig %>% dplyr::group_by(time, pred, agecl) %>% dplyr::summarise(atoutput.tot = sum(atoutput, na.rm = T)) - data.dietcheck = data.dietcheck.orig %>% + data.dietcheck <- data.dietcheck.orig %>% dplyr::left_join(dietcheck.tot, by = c("time", "pred", "agecl")) %>% dplyr::rename(atoutput.old = 'atoutput') %>% dplyr::mutate(atoutput = atoutput.old / atoutput.tot) %>% @@ -261,53 +264,53 @@ process_atl_output = function( saveRDS(data.dietcheck, file = file.path(out.dir, 'data_dietcheck.rds')) } else { ##NEUS only 613 diet obs per timestep - nsteps = 365 / + nsteps <- 365 / extract_prm(prm_biol = param.ls$run.prm, variables = "toutinc") - line.incr = 613 * nsteps + line.incr <- 613 * nsteps if (system == 'windows') { - nline.str = system( + nline.str <- system( paste0('find /c /v "" ', param.ls$dietcheck), intern = T )[2] - nline = as.numeric(strsplit(nline.str, ' ')[[1]][3]) + nline <- as.numeric(strsplit(nline.str, ' ')[[1]][3]) } else { - nline.str = system(paste0('wc -l ', param.ls$dietcheck), intern = T) - nline = as.numeric(strsplit(nline.str, ' ')[[1]][1]) + nline.str <- system(paste0('wc -l ', param.ls$dietcheck), intern = T) + nline <- as.numeric(strsplit(nline.str, ' ')[[1]][1]) } - line.seq = c(seq(0, nline, line.incr), nline) + line.seq <- c(seq(0, nline, line.incr), nline) - diet.agg = list() + diet.agg <- list() - diet.colnames = colnames(data.table::fread(param.ls$dietcheck, nrow = 1)) + diet.colnames <- colnames(data.table::fread(param.ls$dietcheck, nrow = 1)) for (i in 1:(length(line.seq) - 1)) { #read chunk and aggregate by specified interval - lines2read = line.seq[i + 1] - line.seq[i] - diet.slice = data.table::fread( + lines2read <- line.seq[i + 1] - line.seq[i] + diet.slice <- data.table::fread( param.ls$dietcheck, skip = line.seq[i] + 1, nrow = lines2read ) - colnames(diet.slice) = diet.colnames + colnames(diet.slice) <- diet.colnames - diet.slice.dates = as.POSIXct( + diet.slice.dates <- as.POSIXct( diet.slice$Time * 86400, origin = '1964-01-01' ) if (agg.scale == 'month') { - diet.slice$time.agg = as.numeric(factor(format( + diet.slice$time.agg <- as.numeric(factor(format( diet.slice.dates, format = '%m' ))) } else if (agg.scale == 'year') { - diet.slice$time.agg = as.numeric(factor(format( + diet.slice$time.agg <- as.numeric(factor(format( diet.slice.dates, format = '%Y' ))) } else { - diet.slice$time.agg = 1:nrow(diet.slice) + diet.slice$time.agg <- 1:nrow(diet.slice) } - diet.slice = diet.slice %>% + diet.slice <- diet.slice %>% # dplyr::mutate(Time = Time)%>% dplyr::select(-Stock, -Updated) %>% tidyr::gather( @@ -351,17 +354,17 @@ process_atl_output = function( # diet.slice[,6:ncol(diet.slice)] = diet.slice[,6:ncol(diet.slice)]/pred.sum # diet.slice = diet.slice[which(pred.sum !=0),] - diet.agg[[i]] = diet.slice + diet.agg[[i]] <- diet.slice print(i) } - data.dietcheck = dplyr::bind_rows(diet.agg) %>% + data.dietcheck <- dplyr::bind_rows(diet.agg) %>% dplyr::mutate(atoutput = as.numeric(atoutput)) - pred.sum = data.dietcheck %>% + pred.sum <- data.dietcheck %>% dplyr::group_by(time, pred, agecl) %>% dplyr::summarise(atoutput.sum = sum(atoutput, na.rm = T)) - data.dietcheck = data.dietcheck %>% + data.dietcheck <- data.dietcheck %>% dplyr::left_join(pred.sum) %>% dplyr::mutate(atoutput = atoutput / atoutput.sum) %>% dplyr::select(-atoutput.sum) @@ -382,12 +385,12 @@ process_atl_output = function( # Main NetCDF objects ----------------------------------------------------- #Set up biological variable groups - group.types = dplyr::bind_rows(list( + group.types <- dplyr::bind_rows(list( data.frame(species = groups.age, group = 'age'), data.frame(species = groups.bp, group = 'bp') )) - age.vars = c('Nums', 'StructN', 'ResN', 'N') - bp.vars = 'N' + age.vars <- c('Nums', 'StructN', 'ResN', 'N') + bp.vars <- 'N' if ( plot.overall.biomass | @@ -407,31 +410,31 @@ process_atl_output = function( plot.all | process.all ) { - numbers = list() - numbers.age = list() - numbers.box = list() - spatial.numbers = list() - RN.box = list() - SN.box = list() - RN.age = list() - SN.age = list() - RN.age.mean = list() - SN.age.mean = list() - biomass.age = list() - biomass.age.invert = list() - spatial.biomass = list() - spatial.numbers = list() - spatial.biomass.stanza = list() - biomass = list() - biomass.box = list() - sp.overlap = list() - biomass.box.invert = list() - length.age = list() + numbers <- list() + numbers.age <- list() + numbers.box <- list() + spatial.numbers <- list() + RN.box <- list() + SN.box <- list() + RN.age <- list() + SN.age <- list() + RN.age.mean <- list() + SN.age.mean <- list() + biomass.age <- list() + biomass.age.invert <- list() + spatial.biomass <- list() + spatial.numbers <- list() + spatial.biomass.stanza <- list() + biomass <- list() + biomass.box <- list() + sp.overlap <- list() + biomass.box.invert <- list() + length.age <- list() message("Reading in 'Nums', 'StructN', 'ResN', 'N'") if (large.file == F) { - vars = list('Nums', 'StructN', 'ResN', 'N') - group.types = list(groups.age, groups.age, groups.age, groups.bp) - rawdata.main = Map( + vars <- list('Nums', 'StructN', 'ResN', 'N') + group.types <- list(groups.age, groups.age, groups.age, groups.bp) + rawdata.main <- Map( atlantistools::load_nc, select_variable = vars, select_groups = group.types, @@ -456,7 +459,7 @@ process_atl_output = function( process.all | plot.all ) { - spatial.biomass = atlantistools::calculate_biomass_spatial( + spatial.biomass <- atlantistools::calculate_biomass_spatial( nums = rawdata.main[[1]], sn = rawdata.main[[2]], rn = rawdata.main[[3]], @@ -493,21 +496,21 @@ process_atl_output = function( #Biomass timeseries objects if (plot.biomass.timeseries | process.all | plot.all) { - biomass.age = dplyr::filter( + biomass.age <- dplyr::filter( spatial.biomass, species %in% data.age.mat$species ) - biomass.age = atlantistools::agg_data( + biomass.age <- atlantistools::agg_data( biomass.age, groups = c('species', 'agecl', 'time'), fun = sum ) - biomass.age.invert = dplyr::filter( + biomass.age.invert <- dplyr::filter( spatial.biomass, !(species %in% data.age.mat$species) ) - biomass.age.invert = atlantistools::agg_data( + biomass.age.invert <- atlantistools::agg_data( biomass.age.invert, groups = c('species', 'time'), fun = sum @@ -528,7 +531,7 @@ process_atl_output = function( process.all | plot.all ) { - biomass.box = atlantistools::agg_data( + biomass.box <- atlantistools::agg_data( spatial.biomass, groups = c('species', 'polygon', 'time'), fun = sum @@ -537,7 +540,7 @@ process_atl_output = function( #Biomass Seasonal objects if (plot.spatial.biomass.seasonal | process.all | plot.all) { - biomass.box.invert = dplyr::filter( + biomass.box.invert <- dplyr::filter( biomass.box, !(species %in% data.age.mat$species) ) @@ -547,7 +550,7 @@ process_atl_output = function( } if (plot.spatial.biomass | process.all | plot.all) { - spatial.biomass.stanza = atlantistools::combine_ages( + spatial.biomass.stanza <- atlantistools::combine_ages( spatial.biomass, grp_col = 'species', agemat = data.age.mat @@ -558,7 +561,7 @@ process_atl_output = function( } if (plot.weight | plot.all | process.all) { - spatialNumbers = rawdata.main[[1]] %>% + spatialNumbers <- rawdata.main[[1]] %>% dplyr::rename(numbers = atoutput) # filter biomass for species with 10 cohorts and convert to kilograms @@ -596,23 +599,23 @@ process_atl_output = function( process.all | plot.all ) { - RN.age = atlantistools::agg_data( + RN.age <- atlantistools::agg_data( data = rawdata.main[[3]], groups = c('species', 'agecl', 'time'), fun = sum ) - SN.age = atlantistools::agg_data( + SN.age <- atlantistools::agg_data( data = rawdata.main[[2]], groups = c('species', 'agecl', 'time'), fun = sum ) - SN.age.mean = atlantistools::agg_data( + SN.age.mean <- atlantistools::agg_data( data = rawdata.main[[2]], groups = c('species', 'time', 'agecl'), fun = mean ) - RN.age.mean = atlantistools::agg_data( + RN.age.mean <- atlantistools::agg_data( data = rawdata.main[[3]], groups = c('species', 'time', 'agecl'), fun = mean @@ -629,17 +632,17 @@ process_atl_output = function( #Numbers only objects if (plot.numbers.timeseries | plot.cohort | plot.all | process.all) { #numbers - numbers = atlantistools::agg_data( + numbers <- atlantistools::agg_data( data = rawdata.main[[1]], groups = c('species', 'time'), fun = sum ) - numbers.age = atlantistools::agg_data( + numbers.age <- atlantistools::agg_data( data = rawdata.main[[1]], groups = c('species', 'agecl', 'time'), fun = sum ) - numbers.box = atlantistools::agg_data( + numbers.box <- atlantistools::agg_data( data = rawdata.main[[1]], groups = c('species', 'polygon', 'time'), fun = sum @@ -656,34 +659,34 @@ process_atl_output = function( #length-age only objects if (plot.length.age | plot.c.mum | process.all | plot.all) { #Use mean RN+SN per age for each species, convert to weight, get length w/Von Bert. - RN.SN.age = dplyr::left_join( + RN.SN.age <- dplyr::left_join( RN.age.mean, SN.age.mean, by = c('species', 'agecl', 'time') ) - colnames(RN.SN.age) = c('species', 'time', 'agecl', 'RN', 'SN') + colnames(RN.SN.age) <- c('species', 'time', 'agecl', 'RN', 'SN') - biomass.age2 = dplyr::left_join( + biomass.age2 <- dplyr::left_join( RN.SN.age, tempmat3[, 2:4], by = c('species' = 'LongName') ) - biomass.age2$grams_N_Ind = (biomass.age2$RN + biomass.age2$SN) * + biomass.age2$grams_N_Ind <- (biomass.age2$RN + biomass.age2$SN) * 5.7 * 20 / 1000 - biomass.age2$length_age = (as.numeric(as.character( + biomass.age2$length_age <- (as.numeric(as.character( biomass.age2$grams_N_Ind )) / as.numeric(as.character(biomass.age2$li_a)))^(1 / as.numeric(as.character(biomass.age2$li_b))) - length.age = biomass.age2[, c( + length.age <- biomass.age2[, c( 'species', 'agecl', 'time', 'length_age' )] - colnames(length.age)[4] = 'atoutput' + colnames(length.age)[4] <- 'atoutput' bind.save(length.age, 'length_age', out.dir) rm('length.age') @@ -691,12 +694,12 @@ process_atl_output = function( #SN.RN only objects if (plot.sn.rn | process.all | plot.all) { - RN.box = atlantistools::agg_data( + RN.box <- atlantistools::agg_data( data = rawdata.main[[3]], groups = c('species', 'polygon', 'time'), fun = sum ) - SN.box = atlantistools::agg_data( + SN.box <- atlantistools::agg_data( data = rawdata.main[[2]], groups = c('species', 'polygon', 'time'), fun = sum @@ -727,62 +730,62 @@ process_atl_output = function( if (group.types$group[i] == 'age') { #numbers - rawdata.spp = get_rawdata( + rawdata.spp <- get_rawdata( group = group.types$species[i], group.type = group.types$group[i] ) - numbers[[i]] = agg_custom( + numbers[[i]] <- agg_custom( data = rawdata.spp[[1]], groups = c('species', 'time'), fun = sum, agg.scale ) - numbers.age[[i]] = agg_custom( + numbers.age[[i]] <- agg_custom( data = rawdata.spp[[1]], groups = c('species', 'agecl', 'time'), fun = sum, agg.scale ) - numbers.box[[i]] = agg_custom( + numbers.box[[i]] <- agg_custom( data = rawdata.spp[[1]], groups = c('species', 'polygon', 'time'), fun = sum, agg.scale ) - RN.box[[i]] = agg_custom( + RN.box[[i]] <- agg_custom( data = rawdata.spp[[3]], groups = c('species', 'polygon', 'time'), fun = sum, agg.scale ) - SN.box[[i]] = agg_custom( + SN.box[[i]] <- agg_custom( data = rawdata.spp[[2]], groups = c('species', 'polygon', 'time'), fun = sum, agg.scale ) - RN.age[[i]] = agg_custom( + RN.age[[i]] <- agg_custom( data = rawdata.spp[[3]], groups = c('species', 'agecl', 'time'), fun = sum, agg.scale ) - SN.age[[i]] = agg_custom( + SN.age[[i]] <- agg_custom( data = rawdata.spp[[2]], groups = c('species', 'agecl', 'time'), fun = sum, agg.scale ) - SN.age.mean[[i]] = agg_custom( + SN.age.mean[[i]] <- agg_custom( data = rawdata.spp[[2]], groups = c('species', 'time', 'agecl'), fun = mean, agg.scale ) - RN.age.mean[[i]] = agg_custom( + RN.age.mean[[i]] <- agg_custom( data = rawdata.spp[[3]], groups = c('species', 'time', 'agecl'), fun = mean, @@ -791,36 +794,39 @@ process_atl_output = function( #make length.age #Use mean RN+SN per age for each species, convert to weight, get length w/Von Bert. - RN.SN.age = dplyr::left_join( + RN.SN.age <- dplyr::left_join( RN.age[[i]], SN.age[[i]], by = c('species', 'agecl', 'time') ) - colnames(RN.SN.age) = c('species', 'agecl', 'time', 'RN', 'SN') + colnames(RN.SN.age) <- c('species', 'agecl', 'time', 'RN', 'SN') - biomass.age2 = dplyr::left_join( + biomass.age2 <- dplyr::left_join( RN.SN.age, tempmat3[, 2:4], by = c('species' = 'LongName') ) - biomass.age2$grams_N_Ind = (biomass.age2$RN + biomass.age2$SN) * + biomass.age2$grams_N_Ind <- (biomass.age2$RN + biomass.age2$SN) * 5.7 * 20 / 1000 - biomass.age2$length_age = (as.numeric(as.character( + biomass.age2$length_age <- (as.numeric(as.character( biomass.age2$grams_N_Ind )) / as.numeric(as.character(biomass.age2$li_a)))^(1 / as.numeric(as.character(biomass.age2$li_b))) - length.age[[i]] = biomass.age2[, c( + length.age[[i]] <- biomass.age2[, c( 'species', 'agecl', 'time', 'length_age' )] - colnames(length.age[[i]])[4] = 'atoutput' + colnames(length.age[[i]])[4] <- 'atoutput' - spatial.biomass[[i]] = rename(rawdata.spp[[1]], nums = 'atoutput') %>% + spatial.biomass[[i]] <- rename( + rawdata.spp[[1]], + nums = 'atoutput' + ) %>% left_join(rename(rawdata.spp[[2]], sn = 'atoutput')) %>% left_join(rename(rawdata.spp[[3]], rn = 'atoutput')) %>% mutate(atoutput = nums * (sn + rn) * bio.conv) %>% @@ -833,7 +839,7 @@ process_atl_output = function( # vol_dz = NA, # bio_conv = NA, # bps = NA) - spatial.numbers[[i]] = agg_custom( + spatial.numbers[[i]] <- agg_custom( data = rawdata.spp[[1]], groups = c('species', 'agecl', 'polygon', 'layer', 'time'), fun = mean, @@ -841,7 +847,7 @@ process_atl_output = function( ) %>% dplyr::rename(numbers = 'atoutput') - spatial.biomass[[i]] = agg_custom( + spatial.biomass[[i]] <- agg_custom( data = spatial.biomass[[i]], groups = c('species', 'agecl', 'polygon', 'layer', 'time'), fun = mean, @@ -849,13 +855,13 @@ process_atl_output = function( ) #Aggregate spatial biomass based on stanzas - spatial.biomass.stanza[[i]] = atlantistools::combine_ages( + spatial.biomass.stanza[[i]] <- atlantistools::combine_ages( spatial.biomass[[i]], grp_col = 'species', agemat = data.age.mat ) } else { - blank.df = data.frame( + blank.df <- data.frame( species = fgs$LongName[which(fgs$Name == group.types$species[i])], polygon = NA, agecl = NA, @@ -864,12 +870,12 @@ process_atl_output = function( atoutput = NA ) - rawdata.spp = get_rawdata( + rawdata.spp <- get_rawdata( group = group.types$species[i], group.type = group.types$group[i] ) - spatial.biomass[[i]] = atlantistools::calculate_biomass_spatial( + spatial.biomass[[i]] <- atlantistools::calculate_biomass_spatial( nums = blank.df, sn = blank.df, rn = blank.df, @@ -900,34 +906,34 @@ process_atl_output = function( ) #Aggregate biomass by box - biomass.box[[i]] = agg_custom( + biomass.box[[i]] <- agg_custom( spatial.biomass[[i]], groups = c('species', 'polygon', 'time'), fun = sum, agg.scale ) - biomass.box.invert[[i]] = dplyr::filter( + biomass.box.invert[[i]] <- dplyr::filter( biomass.box[[i]], !(species %in% data.age.mat$species) ) #Aggregate by ageclass - biomass.age[[i]] = dplyr::filter( + biomass.age[[i]] <- dplyr::filter( spatial.biomass[[i]], species %in% data.age.mat$species ) - biomass.age[[i]] = agg_custom( + biomass.age[[i]] <- agg_custom( biomass.age[[i]], groups = c('species', 'agecl', 'time'), fun = sum, agg.scale ) - biomass.age.invert[[i]] = dplyr::filter( + biomass.age.invert[[i]] <- dplyr::filter( spatial.biomass[[i]], !(species %in% data.age.mat$species) ) - biomass.age.invert[[i]] = agg_custom( + biomass.age.invert[[i]] <- agg_custom( biomass.age.invert[[i]], groups = c('species', 'time'), fun = sum, @@ -939,7 +945,7 @@ process_atl_output = function( # max age (mean) - biomass / numbers - spatialNumbers = dplyr::bind_rows(spatial.numbers) + spatialNumbers <- dplyr::bind_rows(spatial.numbers) # filter biomass for species with 10 cohorts and convert to kilograms spatialBiomass <- dplyr::bind_rows(spatial.biomass) %>% dplyr::filter(species %in% unique(spatialNumbers$species)) %>% @@ -969,21 +975,21 @@ process_atl_output = function( #Read in raw untransformed data from prod.nc file #Set up biological variable groups - age.vars.prod = c('Eat', 'Growth') - bp.vars.prod = 'Grazing' + age.vars.prod <- c('Eat', 'Growth') + bp.vars.prod <- 'Grazing' #Setup spaces for objects - eat.age = list() - grazing = list() - growth.age = list() - growth.rel.init = list() - bio.consumed = list() + eat.age <- list() + grazing <- list() + growth.age <- list() + growth.rel.init <- list() + bio.consumed <- list() message("Reading in 'Eat', 'Grazing', 'Growth'") if (large.file == F) { - vars = list('Eat', 'Grazing', 'Growth') - group.types = list(groups.age, groups.bp, groups.age) - rawdata.prod = Map( + vars <- list('Eat', 'Grazing', 'Growth') + group.types <- list(groups.age, groups.bp, groups.age) + rawdata.prod <- Map( atlantistools::load_nc, select_variable = vars, select_groups = group.types, @@ -997,22 +1003,22 @@ process_atl_output = function( ) ##Recreate bio.consumed manually without full_join - data_eat = dplyr::bind_rows(rawdata.prod[[1]], rawdata.prod[[2]]) - ts_eat = sort(unique(data_eat$time)) - ts_dm = sort(unique(data.dietcheck$time)) - matching = sum(ts_eat %in% ts_dm) / length(ts_eat) - boxvol = atlantistools::agg_data( + data_eat <- dplyr::bind_rows(rawdata.prod[[1]], rawdata.prod[[2]]) + ts_eat <- sort(unique(data_eat$time)) + ts_dm <- sort(unique(data.dietcheck$time)) + matching <- sum(ts_eat %in% ts_dm) / length(ts_eat) + boxvol <- atlantistools::agg_data( vol, groups = c('polygon', 'time'), out = 'vol', fun = sum ) - pred.names = unique(data.dietcheck$pred) - bio.consumed = list() + pred.names <- unique(data.dietcheck$pred) + bio.consumed <- list() for (i in 1:length(pred.names)) { - consumed_bio = data_eat %>% + consumed_bio <- data_eat %>% dplyr::filter(species == pred.names[i]) %>% dplyr::left_join(boxvol, by = c('polygon', 'time')) %>% dplyr::mutate(atoutput = atoutput * vol) %>% @@ -1023,7 +1029,7 @@ process_atl_output = function( ) %>% dplyr::filter(time %in% ts_eat) %>% dplyr::rename(pred = species) - bio.consumed[[i]] = consumed_bio %>% + bio.consumed[[i]] <- consumed_bio %>% dplyr::filter(!is.na(atoutput.x)) %>% dplyr::filter(!is.na(atoutput.y)) %>% dplyr::mutate(atoutput = atoutput.x * atoutput.y) %>% @@ -1035,11 +1041,11 @@ process_atl_output = function( print(pred.names[i]) gc() } - bio.consumed = dplyr::bind_rows(bio.consumed) + bio.consumed <- dplyr::bind_rows(bio.consumed) saveRDS(bio.consumed, file.path(out.dir, 'biomass_consumed.rds')) rm(bio.consumed) gc() - bio.consumed = atlantistools::calculate_consumed_biomass( + bio.consumed <- atlantistools::calculate_consumed_biomass( eat = rawdata.prod[[1]], grazing = rawdata.prod[[2]], dm = data.dietcheck, @@ -1047,17 +1053,17 @@ process_atl_output = function( bio_conv = bio.conv ) - eat.age = atlantistools::agg_data( + eat.age <- atlantistools::agg_data( data = rawdata.prod[[1]], groups = c('species', 'time', 'agecl'), fun = mean ) - grazing = atlantistools::agg_data( + grazing <- atlantistools::agg_data( data = rawdata.prod[[2]], groups = c('species', 'time'), fun = mean ) - growth.age = atlantistools::agg_data( + growth.age <- atlantistools::agg_data( data = rawdata.prod[[3]], groups = c('species', 'time', 'agecl'), fun = mean @@ -1074,13 +1080,13 @@ process_atl_output = function( } if (group.types$group[i] == 'age') { - prod.vars = age.vars.prod + prod.vars <- age.vars.prod } else { - prod.vars = bp.vars.prod + prod.vars <- bp.vars.prod } ## Process PROD Data by spp - proddata.spp = Map( + proddata.spp <- Map( load_nc_temp, select_variable = prod.vars, select_groups = group.types$species[i], @@ -1097,7 +1103,7 @@ process_atl_output = function( print(i) next() } else if (group.types$group[i] != 'age') { - grazing[[i]] = agg_custom( + grazing[[i]] <- agg_custom( proddata.spp[[1]], groups = c('species', 'agecl', 'polygon', 'time'), fun = mean, @@ -1105,14 +1111,14 @@ process_atl_output = function( ) %>% filter(time %in% unique(data.dietcheck$time)) } else { - eat.age[[i]] = agg_custom( + eat.age[[i]] <- agg_custom( proddata.spp[[1]], groups = c('species', 'agecl', 'polygon', 'time'), fun = mean, agg.scale ) %>% filter(time %in% unique(data.dietcheck$time)) - growth.age[[i]] = agg_custom( + growth.age[[i]] <- agg_custom( proddata.spp[[2]], groups = c('species', 'agecl', 'polygon', 'time'), fun = mean, @@ -1124,12 +1130,12 @@ process_atl_output = function( print(i) } - grazing = dplyr::bind_rows(grazing) - eat.age = dplyr::bind_rows(eat.age) + grazing <- dplyr::bind_rows(grazing) + eat.age <- dplyr::bind_rows(eat.age) - vol.temp = dplyr::filter(vol, time %in% unique(data.dietcheck$time)) + vol.temp <- dplyr::filter(vol, time %in% unique(data.dietcheck$time)) - bio.consumed = calculate_consumed_biomass( + bio.consumed <- calculate_consumed_biomass( eat = eat.age, grazing = grazing, dm = data.dietcheck, @@ -1138,26 +1144,26 @@ process_atl_output = function( ) } - growth.age = dplyr::bind_rows(growth.age) - growth.age = atlantistools::agg_data( + growth.age <- dplyr::bind_rows(growth.age) + growth.age <- atlantistools::agg_data( data = growth.age, groups = c('species', 'time', 'agecl'), fun = mean ) #make growth.rel.init - growth.rel.init = dplyr::left_join( + growth.rel.init <- dplyr::left_join( growth.age, growth.required, by = c("species", "agecl") ) - growth.rel.init = dplyr::mutate( + growth.rel.init <- dplyr::mutate( growth.rel.init, gr_rel = (atoutput - growth_req) / growth_req ) - which.inf = which(growth.rel.init$gr_rel == 'Inf') - growth.rel.init$gr_rel[which.inf] = 1 + which.inf <- which(growth.rel.init$gr_rel == 'Inf') + growth.rel.init$gr_rel[which.inf] <- 1 saveRDS(growth.rel.init, file.path(out.dir, 'growth_rel_init.rds')) saveRDS(growth.age, file.path(out.dir, 'growth_age.rds')) @@ -1165,12 +1171,12 @@ process_atl_output = function( rm(growth.age, growth.rel.init) #aggregate other prod objects - grazing = atlantistools::agg_data( + grazing <- atlantistools::agg_data( data = grazing, groups = c('species', 'time'), fun = mean ) - eat.age = atlantistools::agg_data( + eat.age <- atlantistools::agg_data( data = eat.age, groups = c('species', 'time', 'agecl'), fun = mean @@ -1186,7 +1192,7 @@ process_atl_output = function( # Do Recruitment ---------------------------------------------------------- if (plot.recruits | process.all | plot.all) { - ssb.recruits = atlantistools::load_rec( + ssb.recruits <- atlantistools::load_rec( yoy = param.ls$yoy, ssb = param.ls$ssb, prm_biol = param.ls$biol.prm @@ -1198,7 +1204,7 @@ process_atl_output = function( # Do catch ------------------------------------------------------------------- message("Reading Catch") if (plot.catch | plot.spatial.catch | process.all | plot.all) { - catch = atlantistools::load_nc( + catch <- atlantistools::load_nc( param.ls$catch, fgs = param.ls$groups.file, bps = bio.pools, @@ -1208,7 +1214,7 @@ process_atl_output = function( bboxes = bboxes, check_acronyms = F ) - totcatch = atlantistools::agg_data( + totcatch <- atlantistools::agg_data( catch, groups = c('species', 'time', 'agecl'), fun = sum @@ -1229,7 +1235,7 @@ process_atl_output = function( } if (plot.catch.fleet | process.all | plot.all) { - catch.fleet = process_catch_fleet( + catch.fleet <- process_catch_fleet( fishery.prm = param.ls$fishery.prm, catch = param.ls$catch, groups.file = param.ls$groups.file diff --git a/R/process_catch_fleet.R b/R/process_catch_fleet.R index 196c2c3..f82c887 100644 --- a/R/process_catch_fleet.R +++ b/R/process_catch_fleet.R @@ -20,49 +20,49 @@ #' #' @export -process_catch_fleet = function(fishery.prm, catch, groups.file) { - fisheries = read.csv(fishery.prm) - fgs = read.csv(groups.file) +process_catch_fleet <- function(fishery.prm, catch, groups.file) { + fisheries <- read.csv(fishery.prm) + fgs <- read.csv(groups.file) - catch.nc = ncdf4::nc_open(catch) + catch.nc <- ncdf4::nc_open(catch) - catch.varname = names(catch.nc$var) + catch.varname <- names(catch.nc$var) - times = catch.nc$dim$t$vals / 86400 + times <- catch.nc$dim$t$vals / 86400 - get_fleet_num = function(x) { - a = strsplit(x, '_')[[1]][3] - b = strsplit(a, 'FC')[[1]][2] + get_fleet_num <- function(x) { + a <- strsplit(x, '_')[[1]][3] + b <- strsplit(a, 'FC')[[1]][2] return(as.numeric(b) - 1) } - spp.catch.out.ls = list() - s = 1 + spp.catch.out.ls <- list() + s <- 1 for (s in 1:nrow(fgs)) { - spp.vars = grep(paste0('^', fgs$Code[s], '_'), catch.varname, value = T) + spp.vars <- grep(paste0('^', fgs$Code[s], '_'), catch.varname, value = T) if (length(spp.vars) == 0) { next() } - spp.catch.vars = grep('Catch', spp.vars, value = T) + spp.catch.vars <- grep('Catch', spp.vars, value = T) if (length(spp.catch.vars) == 0) { next() } #get the fleet ID - fleet.nums = sapply(spp.catch.vars, get_fleet_num, USE.NAMES = F) - fleet.names = fisheries$Code[match(fleet.nums, fisheries$Index)] + fleet.nums <- sapply(spp.catch.vars, get_fleet_num, USE.NAMES = F) + fleet.names <- fisheries$Code[match(fleet.nums, fisheries$Index)] - fleet.out.ls = list() - f = 1 + fleet.out.ls <- list() + f <- 1 for (f in 1:length(spp.catch.vars)) { - data.fleet = ncdf4::ncvar_get(catch.nc, spp.catch.vars[f]) %>% + data.fleet <- ncdf4::ncvar_get(catch.nc, spp.catch.vars[f]) %>% as.data.frame() - colnames(data.fleet) = times - fleet.out.ls[[f]] = dplyr::bind_cols( + colnames(data.fleet) <- times + fleet.out.ls[[f]] <- dplyr::bind_cols( data.frame(polygon = (1:nrow(data.fleet)) - 1), data.fleet ) %>% @@ -72,10 +72,10 @@ process_catch_fleet = function(fishery.prm, catch, groups.file) { dplyr::select(species, fleet, polygon, time, atoutput) } - spp.catch.out.ls[[s]] = dplyr::bind_rows(fleet.out.ls) + spp.catch.out.ls[[s]] <- dplyr::bind_rows(fleet.out.ls) } - catch.fleet = dplyr::bind_rows(spp.catch.out.ls) + catch.fleet <- dplyr::bind_rows(spp.catch.out.ls) ncdf4::nc_close(catch.nc)