Added sel_sc argument to plot_traj function

R/plot_traj.R

@@ -10,6 +10,7 @@
 #' if age_class_use = "pairs", to determine which age classes are mature (thus contribute to the number of pairs).
 #' @param onset_year a vector containing the years of each wind farm start being active
 #' (thus, the year at whihc each fatality value starts kicking in)
+#' @param sel_sc scenario to display on the plot. Either "all" or the ID number of a given scenario.
 #' @param xlab a character string. Label for the x axis.
 #' @param ylab a character string. Label for the y axis.
 #' @param Legend a vector of character strings. The legend to show on the side of the plot.
@@ -25,9 +26,13 @@
 #'
 #'
 #'
-plot_traj <- function(N, age_class_use = "NotJuv0", fecundities = NULL, onset_year = NULL,
+plot_traj <- function(N, age_class_use = "NotJuv0", fecundities = NULL, onset_year = NULL, sel_sc = "all",
                       xlab = "Year", ylab = "Population size", Legend = NULL, ylim = NULL, ...){
 
+  # select subset of legends, if needed
+  if(sel_sc != "all") sel_sc = as.numeric(sel_sc)
+  if(sel_sc != "all") Legend = Legend[c(1, sel_sc+1)]
+
   # Get metrics and dimensions
   TH <- dim(N)[2]
   nsc <- dim(N)[3]
@@ -45,7 +50,7 @@ plot_traj <- function(N, age_class_use = "NotJuv0", fecundities = NULL, onset_ye
       N <- N[mature,,,]/2
     }
   }
-  dim(N)
+
 
   if(is.null(onset_year)) onset_year <- 1
   years <- min(onset_year) + (1:TH) - 1
@@ -57,14 +62,20 @@ plot_traj <- function(N, age_class_use = "NotJuv0", fecundities = NULL, onset_ye
   }else{
     out <- colSums(N)
   }
-  dim(out)
+
   N_avg <- apply(out, c(1,2), median)
   N_lci <- apply(out, c(1,2), quantile, prob = pnorm(0.5))
   N_uci <- apply(out, c(1,2), quantile, prob = pnorm(-0.5))
 
   # Build dataframe
-  df <- as.data.frame(cbind(year = years, N_avg = N_avg[,1], N_lci = N_lci[,1], N_uci = N_uci[,1], scenario = 1))
-  for(j in 2:nsc) df <- rbind(df, cbind(year = years, N_avg = N_avg[,j], N_lci = N_lci[,j], N_uci = N_uci[,j], scenario = j))
+  if(sel_sc == "all"){
+    df <- as.data.frame(cbind(year = years, N_avg = N_avg[,1], N_lci = N_lci[,1], N_uci = N_uci[,1], scenario = 1))
+    for(j in 2:nsc) df <- rbind(df, cbind(year = years, N_avg = N_avg[,j], N_lci = N_lci[,j], N_uci = N_uci[,j], scenario = j))
+  }else{
+    df <- as.data.frame(cbind(year = years, N_avg = N_avg[,1], N_lci = N_lci[,1], N_uci = N_uci[,1], scenario = 1))
+    df <- rbind(df, cbind(year = years, N_avg = N_avg[,sel_sc+1], N_lci = N_lci[,sel_sc+1], N_uci = N_uci[,sel_sc+1], scenario = sel_sc+1))
+  }
+
 
   ## Define Graphic Parameters
   size = 1.5
@@ -79,9 +90,10 @@ plot_traj <- function(N, age_class_use = "NotJuv0", fecundities = NULL, onset_ye
   #  geom_ribbon(data = dplyr::filter(df, .data$scenario > 1), aes(ymin = .data$N_lci, ymax = .data$N_uci, fill = factor(.data$scenario)), linetype = 0, alpha = 0.100)
 
   # change color palette (we want sc0 in black)
+  if(sel_sc == "all") ColoR <- custom_palette_c25()[1:nsc] else ColoR <- custom_palette_c25()[c(1, sel_sc + 1)]
+
   p <- p +
-    scale_color_manual(breaks = 1:nsc,
-                       values = custom_palette_c25()[1:nsc],
+    scale_color_manual(values = ColoR,
                        labels = Legend, aesthetics = c("colour", "fill"))
 
 

man/plot_traj.Rd

@@ -9,6 +9,7 @@ plot_traj(
   age_class_use = "NotJuv0",
   fecundities = NULL,
   onset_year = NULL,
+  sel_sc = "all",
   xlab = "Year",
   ylab = "Population size",
   Legend = NULL,
@@ -28,6 +29,8 @@ if age_class_use = "pairs", to determine which age classes are mature (thus cont
 \item{onset_year}{a vector containing the years of each wind farm start being active
 (thus, the year at whihc each fatality value starts kicking in)}
 
+\item{sel_sc}{scenario to display on the plot. Either "all" or the ID number of a given scenario.}
+
 \item{xlab}{a character string. Label for the x axis.}
 
 \item{ylab}{a character string. Label for the y axis.}

run_analysis.R

@@ -157,6 +157,13 @@ get_metrics(N = out$run$N)$scenario$impact[time_horizon, ,-1] %>% round(.,2)
 res = get_metrics(N = out$run$N, cumulated_impacts = cumulated_impacts)
 names(res)
 
+N_0 = N
+
+N=N_0
+
+
+
+
 
 plot_impact(N, sel_sc = "3", show_CI = 0.999, Legend = paste("sc", (1:length(fatalities_mean))-1))