diff --git a/DESCRIPTION b/DESCRIPTION
index cce0b6bc9efd4e8b22ab6c767da13a3c5ca7523a..c61f4a403e3c93798c75730d25a658ee547b7827 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
Package: SWATdoctR
Type: Package
Title: Finding the right diagnoses and treatments for SWAT+ models
-Version: 0.1.14
+Version: 0.1.15
Author: c(person("Christoph", "Schürz",
email = "christoph.schuerz@ufz.de",
role = c("aut", "cre")),
diff --git a/R/plot_mgt_harv.R b/R/plot_mgt_harv.R
index d70fe95ee16c05bb065cf5fc58f8ec5c800ef789..3c8b9eb91e2cf8295915c0ccee1d78b91c327aab 100644
--- a/R/plot_mgt_harv.R
+++ b/R/plot_mgt_harv.R
@@ -34,89 +34,78 @@
#'
plot_variable_at_harvkill <- function(sim_verify, variable, years = 1900:2100) {
##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]])){
+ is_single_sim <- "mgt_out" %in% names(sim_verify)
+ is_list_sim <- "mgt_out" %in% names(sim_verify[[1]])
+ if(is_single_sim){
+ years <- years[years >= min(sim_verify$basin_wb_day$yr)]
+
+ tbl_harv <- prepare_var_at_harvkill(sim_verify$mgt_out, years, variable)
+ } else if (is_list_sim){
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 <- prepare_var_at_harvkill(mgt_out, years, variable)
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, 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, 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, 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, name) %>%
- set_names(c('crop', 'water', 'aeration',
- 'temperature', 'nitrogen', 'phosphorus', 'simulation')) %>%
- pivot_longer(., cols = -c(crop, simulation)) %>%
- mutate(., name = factor(name,
- levels = c('water', 'aeration', 'temperature',
- 'nitrogen', 'phosphorus')))
-
- y_lbl <- 'Plant stress at harvest/kill'
} else {
- stop("Variable must be one of: 'phu', 'yield', 'bioms'")
+ stop('Incorrect input for sim_verify')
}
+ y_lbl <- case_when(variable == 'phu' ~ 'Crop HU fraction before kill',
+ variable == 'yield' ~ 'Crop yield per growing period',
+ variable == 'bioms' ~ 'Biomass at harvest/kill',
+ variable == 'stress' ~ 'Plant stress in growing period')
+
gg <- ggplot(data = tbl_harv) +
labs(x = 'Crops', y = y_lbl) +
theme_bw()
if (variable == 'stress') {
gg <- gg +
- geom_boxplot(aes(x = crop, y = value, fill = name),
+ geom_boxplot(aes(x = crop, y = var, fill = stress),
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')+
- facet_wrap(~simulation)+
- theme(axis.text.x=element_text(angle=45, hjust=1))
+ scale_fill_manual(values = c('deepskyblue4', 'lightcyan3', 'indianred3',
+ 'palegreen4', 'orange3')) +
+ theme(legend.position = 'bottom')
+
+ if (is_list_sim) {
+ gg <- gg + facet_wrap(~simulation)
+ }
} else {
- 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 (is_single_sim) {
+ gg <- gg +
+ geom_boxplot(aes(x = crop, y = var),
+ color = "grey30", outlier.size=1, outlier.colour = "grey50")
+ } else {
+ 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
+#' @param mgt_out Data frame with management output data
+#' @param years Integer vector of years for which values should be calculated
+#' @param variable String to define which variable is selected for plotting
#'
-#' @return dataframe for plotting
+#' @return data frame for plotting
#' @importFrom dplyr %>% filter group_by mutate ungroup
#' @keywords internal
#'
@@ -125,29 +114,229 @@ plot_variable_at_harvkill <- function(sim_verify, variable, years = 1900:2100) {
#' var_at_harv(sim_nostress$mgt_out, 1900:2100)
#' }
-var_at_harv <- function(mgt_out, years){
- tbl_harv_b <- mgt_out %>%
+prepare_var_at_harvkill <- function(mgt_out, years, variable){
+ if(variable == 'phu') {
+ tbl_var <- prepare_phu(mgt_out, years)
+
+ # y_lbl <- 'Crop HU fraction before kill'
+ }else if(variable == 'yield') {
+ tbl_var <- prepare_yield(mgt_out, years)
+
+ # y_lbl <- 'Crop yield per growing period'
+ } else if(variable == 'bioms') {
+ tbl_var <- prepare_biomass(mgt_out, years)
+
+ # y_lbl <- 'Biomass at harvest/kill'
+ } else if(variable == 'stress') {
+ tbl_var <- prepare_stress(mgt_out, years)
+
+ # y_lbl <- 'Plant stress at harvest/kill'
+ } else {
+ stop("Variable must be one of: 'phu', 'yield', 'bioms', 'stress'")
+ }
+
+ return(tbl_var)
+}
+
+#' Prepare mgt_out data for plotting yields
+#'
+#' For static land uses (no plant and kill operations) annual yields are calculated.
+#' For crops with defined growing periods the yield sums of that period is calculated.
+#'
+#' @param mgt_out Data frame with management output data
+#' @param years Integer vector of years for which values should be calculated
+#'
+#' @return data frame with crop names and yields
+#'
+#' @importFrom dplyr bind_rows filter group_by mutate select summarise ungroup %>%
+#' @importFrom purrr set_names
+#'
+#' @keywords internal
+#'
+prepare_yield <- function(mgt_out, years) {
+ yield_tbl <- mgt_out %>%
+ filter(operation %in% c('PLANT','HARVEST', 'KILL')) %>%
+ select(hru, year, operation, op_typ, op_var) %>%
+ group_by(hru, op_typ) %>%
+ mutate(grw_group = ifelse(operation == 'KILL', 1, 0),
+ grw_group = cumsum(grw_group)) %>%
+ ungroup() %>%
+ group_by(hru) %>%
+ mutate(no_plntkill = max(grw_group) == 0)
+
+ yield_mono <- yield_tbl %>%
+ filter(no_plntkill) %>%
+ filter(year %in% years) %>%
+ group_by(hru, year, op_typ) %>%
+ summarise(op_var = sum(op_var), .groups = 'drop') %>%
+ select(-hru, -year)
+
+ yield_rota <- yield_tbl %>%
+ filter(!no_plntkill) %>%
+ filter(operation == 'HARVEST') %>%
+ group_by(hru, op_typ, grw_group) %>%
+ summarise(op_var = sum(op_var), year = max(year), .groups = 'drop') %>%
filter(year %in% years) %>%
- filter(operation %in% c('HARVEST', 'KILL'))
+ select(-hru, -grw_group, -year)
- tbl_harv <- tbl_harv_b %>%
+ yield <- bind_rows(yield_mono, yield_rota) %>%
+ set_names(c('crop', 'var'))
+
+ return(yield)
+}
+
+#' Prepare mgt_out data for plotting biomass
+#'
+#' Biomass is only extracted for harvest/kill operations.
+#'
+#' @param mgt_out Data frame with management output data
+#' @param years Integer vector of years for which values should be calculated
+#'
+#' @return data frame with crop names and bio mass
+#'
+#' @importFrom dplyr filter group_by mutate n select summarise ungroup %>%
+#' @importFrom purrr set_names
+#'
+#' @keywords internal
+#'
+prepare_biomass <- function(mgt_out, years) {
+ bioms <- 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){
- 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')
+ filter(n == 2, operation == 'HARVEST') %>%
+ select(op_typ, plant_bioms) %>%
+ set_names(c('crop', 'var'))
+
+ return(bioms)
+}
+
+#' Prepare mgt_out data for plotting PHU fractions
+#'
+#' PHU fractions are extracted for the last operation before a kill operation is set.
+#'
+#' @param mgt_out Data frame with management output data
+#' @param years Integer vector of years for which values should be calculated
+#'
+#' @return data frame with crop names and phu fractions
+#'
+#' @importFrom dplyr arrange filter group_by lead mutate n select summarise ungroup %>%
+#' @importFrom purrr set_names
+#'
+#' @keywords internal
+#'
+prepare_phu <- function(mgt_out, years) {
+ phu_tbl <- mgt_out %>%
+ mutate(date = 1e4*year + 100*mon + day) %>%
+ select(hru, date, year, operation, op_typ, phuplant) %>%
+ group_by(hru, op_typ, year) %>%
+ mutate(is_kill = ifelse(operation == 'KILL', TRUE, FALSE)) %>%
+ ungroup() %>%
+ arrange(hru, date) %>%
+ mutate(is_kill = is_kill | lead(is_kill)) %>%
+ filter(is_kill) %>%
+ group_by(hru, op_typ, year) %>%
+ mutate(date_diff = diff(date)) %>%
+ ungroup() %>%
+ filter(year %in% years)
+
+ if(any(phu_tbl$date_diff > 1)) {
+ crop_diff <- phu_tbl %>%
+ filter(operation == 'KILL') %>%
+ filter(date_diff > 1) %>%
+ .$op_typ %>%
+ unique(.)
+
+ message('For the following crops the difference between the kill operation\n',
+ 'and the operation before that was more than 1 day:\n',
+ paste(crop_diff, collapse = ', '), '\n',
+ 'The plotted PHU fractions may therefore not be correct.')
}
- return(tbl_harv)
+
+ phu <- phu_tbl %>%
+ filter(operation != 'KILL') %>%
+ select(op_typ, phuplant) %>%
+ set_names(c('crop', 'var'))
+
+ return(phu)
}
+#' Prepare mgt_out data for plotting stress factors
+#'
+#' For static land uses (no plant and kill operations) annual changes in stress
+#' values are calculated (annual increase of a stress factor).
+#' For crops with defined growing periods the stress values for the last harvest
+#' operation is extracted.
+#'
+#' @param mgt_out Data frame with management output data
+#' @param years Integer vector of years for which values should be calculated
+#'
+#' @return data frame with crop names, stress names and values
+#'
+#' @importFrom dplyr filter group_by mutate select summarise ungroup %>%
+#' @importFrom purrr set_names
+#' @importFrom tidyr pivot_longer
+#'
+#' @keywords internal
+#'
+prepare_stress <- function(mgt_out, years) {
+ stress_tbl <- mgt_out %>%
+ filter(operation %in% c('PLANT','HARVEST', 'KILL')) %>%
+ select(hru, year, operation, op_typ, var4, var5, var3, var1, var2) %>%
+ group_by(hru, op_typ) %>%
+ mutate(grw_group = ifelse(operation == 'KILL', 1, 0),
+ grw_group = cumsum(grw_group)) %>%
+ ungroup() %>%
+ group_by(hru) %>%
+ mutate(no_plntkill = max(grw_group) == 0) %>%
+ ungroup()
+
+ stress_mono <- stress_tbl %>%
+ filter(no_plntkill) %>%
+ group_by(hru, year, op_typ) %>%
+ summarise(var4 = max(var4), var5 = max(var5), var3 = max(var3),
+ var2 = max(var2), var1 = max(var1), .groups = 'drop') %>%
+ group_by(hru, op_typ) %>%
+ mutate(var4 = diff1(var4), var5 = diff1(var5), var3 = diff1(var3),
+ var2 = diff1(var2), var1 = diff1(var1)) %>%
+ ungroup() %>%
+ filter(year %in% years) %>%
+ select(-hru, -year)
+
+ stress_rota <- stress_tbl %>%
+ filter(year %in% years) %>%
+ filter(!no_plntkill) %>%
+ filter(operation == 'HARVEST') %>%
+ group_by(hru, op_typ, grw_group) %>%
+ summarise(var4 = max(var4), var5 = max(var5), var3 = max(var3),
+ var2 = max(var2), var1 = max(var1), .groups = 'drop') %>%
+ ungroup(.) %>%
+ select(-hru, -grw_group)
+
+ stress <- bind_rows(stress_mono, stress_rota) %>%
+ set_names(c('crop', 'water', 'aeration',
+ 'temperature', 'nitrogen', 'phosphorus')) %>%
+ pivot_longer(., cols = - crop, names_to = 'stress', values_to = 'var') %>%
+ mutate(., stress = factor(stress,
+ levels = c('water', 'aeration', 'temperature',
+ 'nitrogen', 'phosphorus')))
+
+ return(stress)
+}
+
+#' Concatenate the first value of a vector and the differences of following values.
+#'
+#' @param x numeric vector
+#'
+#' @return numeric vector
+#'
+#' @keywords internal
+#'
+diff1 <- function(x) {
+ c(x[1], diff(x))
+}
#' Filter HRUs by their variable values at harvest/kill
#'
@@ -187,4 +376,3 @@ filter_hru_at_harvkill <- function(sim_verify, ...) {
filter(., ...) %>%
select(-year)
}
-
diff --git a/man/diff1.Rd b/man/diff1.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..6b17d2c5ba086e6561a618860601653f0807153f
--- /dev/null
+++ b/man/diff1.Rd
@@ -0,0 +1,18 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot_mgt_harv.R
+\name{diff1}
+\alias{diff1}
+\title{Concatenate the first value of a vector and the differences of following values.}
+\usage{
+diff1(x)
+}
+\arguments{
+\item{x}{numeric vector}
+}
+\value{
+numeric vector
+}
+\description{
+Concatenate the first value of a vector and the differences of following values.
+}
+\keyword{internal}
diff --git a/man/prepare_biomass.Rd b/man/prepare_biomass.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..e6f451066270a524b181d9b8ae1d3a7079439715
--- /dev/null
+++ b/man/prepare_biomass.Rd
@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot_mgt_harv.R
+\name{prepare_biomass}
+\alias{prepare_biomass}
+\title{Prepare mgt_out data for plotting biomass}
+\usage{
+prepare_biomass(mgt_out, years)
+}
+\arguments{
+\item{mgt_out}{Data frame with management output data}
+
+\item{years}{Integer vector of years for which values should be calculated}
+}
+\value{
+data frame with crop names and bio mass
+}
+\description{
+Biomass is only extracted for harvest/kill operations.
+}
+\keyword{internal}
diff --git a/man/prepare_phu.Rd b/man/prepare_phu.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..b587e7481ae53207832ee8b115b25699f54f5179
--- /dev/null
+++ b/man/prepare_phu.Rd
@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot_mgt_harv.R
+\name{prepare_phu}
+\alias{prepare_phu}
+\title{Prepare mgt_out data for plotting PHU fractions}
+\usage{
+prepare_phu(mgt_out, years)
+}
+\arguments{
+\item{mgt_out}{Data frame with management output data}
+
+\item{years}{Integer vector of years for which values should be calculated}
+}
+\value{
+data frame with crop names and phu fractions
+}
+\description{
+PHU fractions are extracted for the last operation before a kill operation is set.
+}
+\keyword{internal}
diff --git a/man/prepare_stress.Rd b/man/prepare_stress.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..ca4959958822dc1233bcdd1a5df88692cfe27bcd
--- /dev/null
+++ b/man/prepare_stress.Rd
@@ -0,0 +1,23 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot_mgt_harv.R
+\name{prepare_stress}
+\alias{prepare_stress}
+\title{Prepare mgt_out data for plotting stress factors}
+\usage{
+prepare_stress(mgt_out, years)
+}
+\arguments{
+\item{mgt_out}{Data frame with management output data}
+
+\item{years}{Integer vector of years for which values should be calculated}
+}
+\value{
+data frame with crop names, stress names and values
+}
+\description{
+For static land uses (no plant and kill operations) annual changes in stress
+values are calculated (annual increase of a stress factor).
+For crops with defined growing periods the stress values for the last harvest
+operation is extracted.
+}
+\keyword{internal}
diff --git a/man/var_at_harv.Rd b/man/prepare_var_at_harvkill.Rd
similarity index 50%
rename from man/var_at_harv.Rd
rename to man/prepare_var_at_harvkill.Rd
index 303e994683d281b76641af01d8b69fe9c18fc7b6..708e2916b03bee01e02f5b486fdb79535ef2f90b 100644
--- a/man/var_at_harv.Rd
+++ b/man/prepare_var_at_harvkill.Rd
@@ -1,18 +1,20 @@
% 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}
+\name{prepare_var_at_harvkill}
+\alias{prepare_var_at_harvkill}
\title{Prepare dataframe for plot_variable_at_harvkill function}
\usage{
-var_at_harv(mgt_out, years)
+prepare_var_at_harvkill(mgt_out, years, variable)
}
\arguments{
-\item{mgt_out}{dataframe with model output data}
+\item{mgt_out}{Data frame with management output data}
-\item{years}{integer vector to filter dataframe}
+\item{years}{Integer vector of years for which values should be calculated}
+
+\item{variable}{String to define which variable is selected for plotting}
}
\value{
-dataframe for plotting
+data frame for plotting
}
\description{
Prepare dataframe for plot_variable_at_harvkill function
diff --git a/man/prepare_yield.Rd b/man/prepare_yield.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..f54cf8296304b1eb559be17955f5c34f66c26eb1
--- /dev/null
+++ b/man/prepare_yield.Rd
@@ -0,0 +1,21 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot_mgt_harv.R
+\name{prepare_yield}
+\alias{prepare_yield}
+\title{Prepare mgt_out data for plotting yields}
+\usage{
+prepare_yield(mgt_out, years)
+}
+\arguments{
+\item{mgt_out}{Data frame with management output data}
+
+\item{years}{Integer vector of years for which values should be calculated}
+}
+\value{
+data frame with crop names and yields
+}
+\description{
+For static land uses (no plant and kill operations) annual yields are calculated.
+For crops with defined growing periods the yield sums of that period is calculated.
+}
+\keyword{internal}