diff --git a/R/globals.R b/R/globals.R
index 351308e6c2cc1bc5f12617ef84d630ec6f9dbf72..ce423d3caa3fcf0de2ee2b4beb9359b3f887b765 100644
--- a/R/globals.R
+++ b/R/globals.R
@@ -8,4 +8,4 @@ utils::globalVariables(c("%&&%", "%//%", ".", "crop", "day", "ecanopy", "eplant"
"wndspd", "yr", "rhum", "rm_skp", "Date", "Values", "surq_gen", "latq", "perc",
"description", "flo", "..density..", "wateryld", "cn", "sw_final", "esr", "mmdd",
"lu_mgt_ini", "has_suffix", "suffix_upd", "lu_mgt_upd", "schedule_upd",
- "plnt_com_upd", "n_chr"))
+ "plnt_com_upd", "n_chr", "simulation"))
diff --git a/R/plot_mgt_harv.R b/R/plot_mgt_harv.R
index c5d6db7ffca527ba8656f56a6e4ad771397d4e8c..fdc9d7102b5f3e67c9fd75793e383e171ba1c863 100644
--- a/R/plot_mgt_harv.R
+++ b/R/plot_mgt_harv.R
@@ -4,10 +4,11 @@
#' fractions ('phu'), crop yields ('yield'), or plant biomass ('bioms') for crops
#' at harvest-kill of a crop separated for all identified crops.
#'
-#' @param sim_verify Simulation output of the function
-#' \code{run_swat_verification()}.
+#' @param sim_verify Simulation output of the function \code{run_swat_verification()}
+#' as a single object or as a list of objects (representing \code{run_swat_verification()}
+#' function runs with different settings).
#' To plot the heat units at least the output option \code{outputs = 'mgt'}
-#' must be set in \code{run_swat_verification()}
+#' must be set in \code{run_swat_verification()}.
#' @param variable Selected variable to be plotted. Must be one of: 'phu',
#' 'yield', 'bioms'
#' @param years Simulated years which are aggregated in the boxplot
@@ -19,51 +20,62 @@
#' @importFrom purrr set_names
#'
#' @export
+#' @examples
+#' \dontrun{
+#' ##As single object
+#' plot_variable_at_harvkill(sim_nostress, 'yield')
+#' ##As a list
+#' sim_list<- list(
+#' no_stress = sim_nostress,
+#' only_nutrient = sim_nutrientstress,
+#' all_stress = sim_allstress)
+#' plot_variable_at_harvkill(sim_list, 'yield')
+#' }
#'
plot_variable_at_harvkill <- function(sim_verify, variable, years = 1900:2100) {
- tbl_harv <- sim_verify$mgt_out %>%
- filter(year %in% years) %>%
- filter(operation %in% c('HARVEST', 'KILL')) %>%
- group_by(hru,year, mon, day) %>%
- mutate(n = n()) %>%
- ungroup() %>%
- filter(n == 2, operation == 'HARVEST')
-
- ##Just in case dates of harvest and kill differ a bit
- if(dim(tbl_harv)[1] == 0){
- tbl_harv <- sim_verify$mgt_out %>%
- filter(year %in% years) %>%
- filter(operation %in% c('HARVEST', 'KILL')) %>%
- group_by(hru,year, mon) %>%
- mutate(n = n()) %>%
- ungroup() %>%
- filter(n == 2, operation == 'HARVEST')
+ ##Identifying if it is single object of a list and generating df for plotting
+ if("mgt_out" %in% names(sim_verify)){
+ tbl_harv <- var_at_harv(sim_verify$mgt_out, years)
+ tbl_harv$name <- deparse(substitute(sim_verify))
+ } else if ("mgt_out" %in% names(sim_verify[[1]])){
+ df <- NULL
+ nn <- c()
+ for(n in names(sim_verify)){
+ mgt_out <- sim_verify[[n]][["mgt_out"]]
+ tbl_harv <- var_at_harv(mgt_out, years)
+ tbl_harv$name <- n
+ nn <- c(nn, n)
+ if(is.null(df)){df <- tbl_harv} else {df <- bind_rows(df, tbl_harv)}
+ }
+ tbl_harv <- df %>%
+ mutate(name = factor(name, levels = nn))
}
+
if(variable == 'phu') {
tbl_harv <- tbl_harv %>%
- select(op_typ, phuplant) %>%
- set_names(c('crop', 'var'))
+ select(op_typ, phuplant, name) %>%
+ set_names(c('crop', 'var', 'simulation'))
y_lbl <- 'Crop HUs at harvest/kill'
}else if(variable == 'yield') {
tbl_harv <- tbl_harv %>%
- select(op_typ, op_var) %>%
- set_names(c('crop', 'var'))
+ select(op_typ, op_var, name) %>%
+ set_names(c('crop', 'var', 'simulation'))
y_lbl <- 'Crop yield at harvest/kill'
} else if(variable == 'bioms') {
tbl_harv <- tbl_harv %>%
- select(op_typ, plant_bioms) %>%
- set_names(c('crop', 'var'))
+ select(op_typ, plant_bioms, name) %>%
+ set_names(c('crop', 'var', 'simulation'))
y_lbl <- 'Biomass at harvest/kill'
} else if(variable == 'stress') {
tbl_harv <- tbl_harv %>%
- select(op_typ, var4, var5, var3, var1, var2) %>%
+ select(op_typ, var4, var5, var3, var1, var2, name) %>%
set_names(c('crop', 'water', 'aeration',
- 'temperature', 'nitrogen', 'phosphorus')) %>%
- pivot_longer(., cols = -crop) %>%
+ 'temperature', 'nitrogen', 'phosphorus', 'simulation')) %>%
+ pivot_longer(., cols = -c(crop, simulation)) %>%
mutate(., name = factor(name,
levels = c('water', 'aeration', 'temperature',
'nitrogen', 'phosphorus')))
@@ -79,19 +91,59 @@ plot_variable_at_harvkill <- function(sim_verify, variable, years = 1900:2100) {
if (variable == 'stress') {
gg <- gg +
- geom_boxplot(aes(x = crop, y = value, fill = name)) +
+ geom_boxplot(aes(x = crop, y = value, fill = name),
+ color = "grey30", outlier.size=1, outlier.colour = "grey50") +
labs(fill = 'Stress factor') +
scale_fill_manual(values = c('deepskyblue4', 'lightcyan3',
- 'indianred3', 'palegreen4', 'orange3')) +
- theme(legend.position = 'bottom')
+ 'indianred3', 'palegreen4', 'orange3')) +
+ theme(legend.position = 'bottom')+
+ facet_wrap(~simulation)+
+ theme(axis.text.x=element_text(angle=45, hjust=1))
} else {
- gg <- gg + geom_boxplot(aes(x = crop, y = var))
+ gg <- gg + geom_boxplot(aes(x = crop, y = var, fill = simulation),
+ color = "grey30", outlier.size=1, outlier.colour = "grey50")+
+ scale_fill_manual(values = c('white', 'palegreen4', 'orange3', 'lightcyan3', 'indianred3'))
}
-
if(variable == 'phu') {
gg <- gg + geom_hline(yintercept = 1, lty = 2)
}
-
return(gg)
}
+
+#' Prepare dataframe for plot_variable_at_harvkill function
+#'
+#' @param mgt_out dataframe with model output data
+#' @param years integer vector to filter dataframe
+#'
+#' @return dataframe for plotting
+#' @importFrom dplyr %>% filter group_by mutate ungroup
+#' @keywords internal
+#'
+#' @examples
+#' \dontrun{
+#' var_at_harv(sim_nostress$mgt_out, 1900:2100)
+#' }
+
+var_at_harv <- function(mgt_out, years){
+ tbl_harv_b <- mgt_out %>%
+ filter(year %in% years) %>%
+ filter(operation %in% c('HARVEST', 'KILL'))
+
+ tbl_harv <- tbl_harv_b %>%
+ group_by(hru, year, mon, day) %>%
+ mutate(n = n()) %>%
+ ungroup() %>%
+ filter(n == 2, operation == 'HARVEST')
+
+ ##Just in case dates of harvest and kill differ a bit
+ if(dim(tbl_harv)[1] == 0){
+ print("harv and kill are not on the same day. Switching to closest on same month.")
+ tbl_harv <- tbl_harv_b %>%
+ group_by(hru, year, mon) %>%
+ mutate(n = n()) %>%
+ ungroup() %>%
+ filter(n == 2, operation == 'HARVEST')
+ }
+ return(tbl_harv)
+}
diff --git a/man/plot_variable_at_harvkill.Rd b/man/plot_variable_at_harvkill.Rd
index 90e2ca83a66b605edb43814513fb81d8d3361684..0f29b8a9e55d2bbbaa4a4ad684169549096417eb 100644
--- a/man/plot_variable_at_harvkill.Rd
+++ b/man/plot_variable_at_harvkill.Rd
@@ -7,10 +7,11 @@
plot_variable_at_harvkill(sim_verify, variable, years = 1900:2100)
}
\arguments{
-\item{sim_verify}{Simulation output of the function
-\code{run_swat_verification()}.
-To plot the heat units at least the output option \code{outputs = 'mgt'}
-must be set in \code{run_swat_verification()}}
+\item{sim_verify}{Simulation output of the function \code{run_swat_verification()}
+as a single object or as a list of objects (representing \code{run_swat_verification()}
+function runs with different settings).
+ To plot the heat units at least the output option \code{outputs = 'mgt'}
+ must be set in \code{run_swat_verification()}.}
\item{variable}{Selected variable to be plotted. Must be one of: 'phu',
'yield', 'bioms'}
@@ -25,3 +26,16 @@ plot_variable_at_harvkill plots boxplots of one of the variables crop heat unit
fractions ('phu'), crop yields ('yield'), or plant biomass ('bioms') for crops
at harvest-kill of a crop separated for all identified crops.
}
+\examples{
+\dontrun{
+##As single object
+plot_variable_at_harvkill(sim_nostress, 'yield')
+##As a list
+sim_list<- list(
+no_stress = sim_nostress,
+only_nutrient = sim_nutrientstress,
+all_stress = sim_allstress)
+plot_variable_at_harvkill(sim_list, 'yield')
+}
+
+}
diff --git a/man/var_at_harv.Rd b/man/var_at_harv.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..303e994683d281b76641af01d8b69fe9c18fc7b6
--- /dev/null
+++ b/man/var_at_harv.Rd
@@ -0,0 +1,25 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot_mgt_harv.R
+\name{var_at_harv}
+\alias{var_at_harv}
+\title{Prepare dataframe for plot_variable_at_harvkill function}
+\usage{
+var_at_harv(mgt_out, years)
+}
+\arguments{
+\item{mgt_out}{dataframe with model output data}
+
+\item{years}{integer vector to filter dataframe}
+}
+\value{
+dataframe for plotting
+}
+\description{
+Prepare dataframe for plot_variable_at_harvkill function
+}
+\examples{
+\dontrun{
+var_at_harv(sim_nostress$mgt_out, 1900:2100)
+}
+}
+\keyword{internal}