functions for ps and tile plotting added

NAMESPACE

@@ -8,6 +8,8 @@ export(plot_hru_pw_day)
 export(plot_hru_var)
 export(plot_hru_var_aa)
 export(plot_monthly_snow)
+export(plot_ps)
+export(plot_qtile)
 export(plot_variable_at_harvkill)
 export(plot_water_partition)
 export(print_avannual_qtile)
@@ -20,6 +22,9 @@ import(patchwork)
 importFrom(data.table,fread)
 importFrom(dplyr,"%>%")
 importFrom(dplyr,arrange)
+importFrom(dplyr,bind_rows)
+importFrom(dplyr,case_when)
+importFrom(dplyr,distinct)
 importFrom(dplyr,filter)
 importFrom(dplyr,group_by)
 importFrom(dplyr,group_map)
@@ -40,12 +45,17 @@ importFrom(dplyr,summarise_all)
 importFrom(dplyr,tibble)
 importFrom(dplyr,ungroup)
 importFrom(ggplot2,aes)
+importFrom(ggplot2,facet_wrap)
 importFrom(ggplot2,geom_boxplot)
+importFrom(ggplot2,geom_density)
+importFrom(ggplot2,geom_histogram)
 importFrom(ggplot2,geom_hline)
+importFrom(ggplot2,geom_line)
 importFrom(ggplot2,geom_text)
 importFrom(ggplot2,ggplot)
 importFrom(ggplot2,labs)
 importFrom(ggplot2,lims)
+importFrom(ggplot2,theme)
 importFrom(ggplot2,theme_bw)
 importFrom(ggplot2,theme_void)
 importFrom(lubridate,floor_date)

R/globals.R

@@ -6,4 +6,4 @@ utils::globalVariables(c("%&&%", "%//%", ".", "crop", "day", "ecanopy", "eplant"
                          "tile", "time_out", "tmn", "tmpav", "tmx", "topo", "val_max", "val_mean",
                          "val_min", "value", "value_sum", "var", "var1", "var2", "var3", "var4", "var5",
                          "wndspd", "yr", "rhum", "rm_skp", "Date", "Values", "surq_gen", "latq", "perc",
-                         "description"))
+                         "description", "flo", "..density.."))

R/plot_hru_pw.R

@@ -79,38 +79,6 @@ plot_hru_pw_day <- function(sim_verify, hru_id, var, title = "", years = 1900:21
           axis.title.y = element_blank())
 }
 
-#' Print the average annual qtile for HRUs
-#'
-#' print_avannual_qtile prints a table with the average annual qtile in mm
-#' for HRUs that used a tile flow parametrization in landuse.lum
-#'
-#' @param sim_verify Simulation output of the function \code{run_swat_verification()}.
-#'   To plot the heat units at least the output option \code{outputs = 'wb'} must
-#'   be set in  \code{run_swat_verification()}
-#' @param exclude_lum Character vector to define land uses which are excluded
-#'   in the printed table.
-#'
-#' @importFrom dplyr arrange filter left_join rename select %>%
-#'
-#' @return Returns a table with hru ids average annual qtile and attributes.
-#'
-#' @export
-#'
-print_avannual_qtile <- function(sim_verify,
-                                 exclude_lum = c(
-                                   "urhd_lum", "urmd_lum", "urml_lum",
-                                   "urld_lum", "ucom_lum", "uidu_lum",
-                                   "utrn_lum", "uins_lum", "urbn_lum"
-                                 )) {
-
-  sim_verify$hru_wb_aa %>%
-    rename(id = unit) %>%
-    left_join(., sim_verify$lum_mgt, by = "id") %>%
-    filter(tile != 'null') %>%
-    filter(!lu_mgt %in% exclude_lum) %>%
-    select(id, qtile, lu_mgt, mgt, soil) %>%
-    arrange(qtile, id)
-}
 
 #' Aggregate and plot simulated variables saved in hru_pw_day
 #'
@@ -173,7 +141,11 @@ plot_hru_var <- function(sim_verify, hru_id, var, period = "day", fn_summarize =
 #' }
 
 plot_hru_var_aa <- function(sim_verify, lum = NULL, mgt = NULL, soil = NULL){
-  p <- paste(lum, mgt, soil)
+  p <- if(is.null(lum) & is.null(mgt) & is.null(soil)){
+    p <- "all"
+  } else {
+    p <- paste0(lum,"|",mgt,"|",soil)
+  }
   id <- get_hru_id_by_attribute(sim_verify, lum, mgt, soil)
   fig <- sim_verify$hru_wb_aa[sim_verify$hru_wb_aa$unit %in% id$id, -c(1:7)] %>%
     .[, colSums(.!= 0) > 0] %>%
@@ -182,7 +154,7 @@ plot_hru_var_aa <- function(sim_verify, lum = NULL, mgt = NULL, soil = NULL){
     group_by(p, var) %>%
     group_map(~plot_ly(., y=~Values, color = ~var, colors = "cyan4", type = 'box'), keep = TRUE) %>%
     subplot(nrows = 5) %>%
-    layout(title = paste("HRUs selected by", p))
+    layout(title = paste0("HRUs selected: ", p))
   return(fig)
 }
 
@@ -197,6 +169,8 @@ plot_hru_var_aa <- function(sim_verify, lum = NULL, mgt = NULL, soil = NULL){
 #' @param soil Optional character vector with soil type labels as defined in the soil data.
 #' @param exclude_lum Character vector to define land uses which are excluded
 #'   in the printed table.
+#' @param boxpoints Optional Boolean TRUE for displaying outliers, FALSE for hiding them.
+#'  \code{Default = TRUE}
 #' @return plotly figure object
 #' @importFrom dplyr %>% mutate group_by rename left_join summarise_all filter select
 #' @importFrom tidyr pivot_longer
@@ -211,7 +185,7 @@ plot_hru_var_aa <- function(sim_verify, lum = NULL, mgt = NULL, soil = NULL){
 plot_water_partition <- function(sim_verify, tile = NULL, lum = NULL, mgt = NULL, soil = NULL, exclude_lum = c(
   "urhd_lum", "urmd_lum", "urml_lum",
   "urld_lum", "ucom_lum", "uidu_lum",
-  "utrn_lum", "uins_lum", "urbn_lum")){
+  "utrn_lum", "uins_lum", "urbn_lum"), boxpoints = TRUE){
   df <- sim_verify$hru_wb_aa %>%
     rename(id = unit) %>%
     left_join(., sim_verify$lum_mgt, by = "id") %>%
@@ -261,7 +235,7 @@ plot_water_partition <- function(sim_verify, tile = NULL, lum = NULL, mgt = NULL
     add_pie(hole = 0.3)
   ##Preparing box plot
   box_pl <- plot_ly(df[c("var", "Values")], x=~Values,  color = ~var,  type = "box",  colors = pal,
-                    showlegend = F) %>%
+                    showlegend = F,  boxpoints = boxpoints) %>%
     layout(yaxis = list(autorange = "reversed"))
   ##Putting into one figure and annotations
   fig <- subplot(box_pl, pie_pl, nrows = 1, margin = 0.05) %>%

R/plot_ps_tile.R

+#' Plot yearly simulated point source values (as yearly loads or average concentrations)
+#'
+#' @param sim_verify Simulation output of the function \code{run_swat_verification()}.
+#'   To plot the heat units at least the output option \code{outputs = 'wb'} must
+#'   be set in  \code{run_swat_verification()}
+#' @param conc Boolean, TRUE to provide figure for point source average yearly pollutant concentrations,
+#' FALSE - for average yearly pollutant loads.
+#' @return ggplot object for point source yearly simulated data
+#' @importFrom dplyr %>% mutate group_by group_map select bind_rows case_when
+#' @importFrom ggplot2 ggplot geom_line facet_wrap labs theme_bw theme aes
+#' @importFrom tidyr pivot_longer
+#' @export
+#'
+#' @examples
+#' \dontrun{
+#' plot_ps(sim_nostress, TRUE)
+#' }
+
+plot_ps <- function(sim_verify, conc = FALSE){
+  df <- sim_verify$recall_yr[, -c(1:3,6)] %>%
+    .[, colSums(.!= 0) > 0] %>%
+    mutate(yr = as.factor(yr),
+           name = gsub("hru00", "ps", name))
+  if(conc){
+    df <- df %>%
+      pivot_longer(-c(yr, name, flo), names_to = 'var', values_to = 'Values') %>%
+      mutate(Values = ifelse(var == "sed", (Values/flo)*(1000000/(24*69*60*365.25)),
+                             (Values/flo)*(1000/(24*69*60*365.25)))) %>%
+      select(yr, name, var, Values) %>%
+      bind_rows(df[c("yr", "name", "flo")] %>% mutate(var = "flo") %>% rename(Values = flo)) %>%
+      mutate(var = case_when(var == 'flo' ~ "flo m3/s",
+                             var %in% c("orgn", "no3", "nh3", "no2") ~ paste(var, "N mg/y"),
+                             var %in% c("sedp", "solp") ~ paste(var, "P mg/y"),
+                             TRUE ~  paste(var, "mg/l")))
+  } else {
+    df <- df %>%
+      pivot_longer(-c(yr, name), names_to = 'var', values_to = 'Values') %>%
+      mutate(var = case_when(var == 'flo' ~ "flo m3/s",
+                             var == 'sed' ~ "sed t/y",
+                             var %in% c("orgn", "no3", "nh3", "no2") ~ paste(var, "N kg/y"),
+                             var %in% c("sedp", "solp") ~ paste(var, "P kg/y"),
+                             TRUE ~  paste(var, "kg/y")))
+  }
+  fig <- ggplot(df, aes(x=yr, y=Values,  group=name, colour=name))+
+    geom_line(size=1.5)+
+    facet_wrap(~var, scales = "free_y")+
+    labs(color='Point sources', x = 'Year') +
+    theme_bw()+
+    theme(strip.background = element_rect(fill = "deepskyblue3", colour = "azure3"),
+          strip.text = element_text(color = "white", face="bold"),
+          panel.border = element_rect(colour = "azure3"),
+          axis.text.x = element_text(angle = 25, hjust=1))
+  return(fig)
+}
+
+#' Print the average annual qtile for HRUs
+#'
+#' print_avannual_qtile prints a table with the average annual qtile in mm
+#' for HRUs that used a tile flow parametrization in landuse.lum
+#'
+#' @param sim_verify Simulation output of the function \code{run_swat_verification()}.
+#'   To plot the heat units at least the output option \code{outputs = 'wb'} must
+#'   be set in  \code{run_swat_verification()}
+#' @param exclude_lum Character vector to define land uses which are excluded
+#'   in the printed table.
+#'
+#' @importFrom dplyr arrange filter left_join rename select %>%
+#'
+#' @return Returns a table with hru ids average annual qtile and attributes.
+#'
+#' @export
+#'
+print_avannual_qtile <- function(sim_verify,
+                                 exclude_lum = c(
+                                   "urhd_lum", "urmd_lum", "urml_lum",
+                                   "urld_lum", "ucom_lum", "uidu_lum",
+                                   "utrn_lum", "uins_lum", "urbn_lum"
+                                 )) {
+
+  sim_verify$hru_wb_aa %>%
+    rename(id = unit) %>%
+    left_join(., sim_verify$lum_mgt, by = "id") %>%
+    filter(tile != 'null') %>%
+    filter(!lu_mgt %in% exclude_lum) %>%
+    select(id, qtile, lu_mgt, mgt, soil) %>%
+    arrange(qtile, id)
+}
+
+#' Plot tile drain flow histogram and distribution curve
+#'
+#' @param sim_verify Simulation output of the function \code{run_swat_verification()}.
+#'   To plot the heat units at least the output option \code{outputs = 'wb'} must
+#'   be set in  \code{run_swat_verification()}
+#' @param exclude_lum Character vector to define land uses which are excluded
+#'   in the printed table.
+#'
+#' @importFrom ggplot2 ggplot geom_histogram labs theme_bw theme geom_density aes
+#' @return Returns ggplot object with histogram and density curve for tile drain from distribution.
+#' @export
+#'
+#' @examples
+#' \dontrun{
+#' plot_qtile(sim_nostress)
+#' }
+
+plot_qtile <- function(sim_verify, exclude_lum = c("urhd_lum", "urmd_lum", "urml_lum",
+                                                   "urld_lum", "ucom_lum", "uidu_lum",
+                                                   "utrn_lum", "uins_lum", "urbn_lum")){
+  df <- print_avannual_qtile(sim_verify, exclude_lum)
+
+  fig <- ggplot(df, aes(x=qtile)) +
+    geom_histogram(aes(y=..density..), color="black", fill="blue", breaks = seq(min(df$qtile), max(df$qtile), 10))+
+    geom_density(alpha=.3, fill="white", linewidth = 1, color = "grey25", linetype = "twodash")+
+    labs(title = "Tile drain flow density mm/year")+
+    theme_bw()+
+    theme(panel.border = element_blank(),
+          axis.line = element_line(color='black'),
+          axis.title.x=element_blank())
+  return(fig)
+}

R/run_swat_verify.R

@@ -31,7 +31,7 @@
 #' @return Returns the simulation results for the defined output variables as a
 #'   list of tibbles.
 #'
-#' @importFrom dplyr %>%
+#' @importFrom dplyr %>% distinct
 #' @importFrom processx run
 #' @importFrom stringr str_split
 #' @importFrom tibble tibble
@@ -77,6 +77,7 @@ run_swat_verification <- function(project_path, outputs = c('wb', 'mgt', 'plt'),
       model_output$basin_wb_day <- read_tbl('basin_wb_day.txt', run_path, 3)
       model_output$basin_pw_day <- read_tbl('basin_pw_day.txt', run_path, 3)
       model_output$hru_wb_aa <- read_tbl('hru_wb_aa.txt', run_path, 3)
+      model_output$recall_yr <- read_tbl('recall_yr.txt', run_path, 3)
     }
     if ('mgt' %in% outputs) {
       model_output$mgt_out <- read_mgt(run_path)
@@ -86,6 +87,7 @@ run_swat_verification <- function(project_path, outputs = c('wb', 'mgt', 'plt'),
       landuse_lum <- read_tbl('landuse.lum', run_path, 2)
       model_output$lum_mgt <- left_join(hru_data, landuse_lum,
                                         by = c("lu_mgt" = 'name')) %>%
+        distinct() %>%
         select(id, topo, hydro, soil, lu_mgt, plnt_com, mgt, tile)
     }
   }

man/plot_ps.Rd

+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot_ps_tile.R
+\name{plot_ps}
+\alias{plot_ps}
+\title{Plot yearly simulated point source values (as yearly loads or average concentrations)}
+\usage{
+plot_ps(sim_verify, conc = FALSE)
+}
+\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 = 'wb'} must
+be set in  \code{run_swat_verification()}}
+
+\item{conc}{Boolean, TRUE to provide figure for point source average yearly pollutant concentrations,
+FALSE - for average yearly pollutant loads.}
+}
+\value{
+ggplot object for point source yearly simulated data
+}
+\description{
+Plot yearly simulated point source values (as yearly loads or average concentrations)
+}
+\examples{
+\dontrun{
+plot_ps(sim_nostress, TRUE)
+}
+}

man/plot_qtile.Rd

+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot_ps_tile.R
+\name{plot_qtile}
+\alias{plot_qtile}
+\title{Plot tile drain flow histogram and distribution curve}
+\usage{
+plot_qtile(
+  sim_verify,
+  exclude_lum = c("urhd_lum", "urmd_lum", "urml_lum", "urld_lum", "ucom_lum", "uidu_lum",
+    "utrn_lum", "uins_lum", "urbn_lum")
+)
+}
+\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 = 'wb'} must
+be set in  \code{run_swat_verification()}}
+
+\item{exclude_lum}{Character vector to define land uses which are excluded
+in the printed table.}
+}
+\value{
+Returns ggplot object with histogram and density curve for tile drain from distribution.
+}
+\description{
+Plot tile drain flow histogram and distribution curve
+}
+\examples{
+\dontrun{
+plot_qtile(sim_nostress)
+}
+}

man/plot_water_partition.Rd

@@ -6,12 +6,13 @@
 \usage{
 plot_water_partition(
   sim_verify,
-  tile = TRUE,
+  tile = NULL,
   lum = NULL,
   mgt = NULL,
   soil = NULL,
   exclude_lum = c("urhd_lum", "urmd_lum", "urml_lum", "urld_lum", "ucom_lum", "uidu_lum",
-    "utrn_lum", "uins_lum", "urbn_lum")
+    "utrn_lum", "uins_lum", "urbn_lum"),
+  boxpoints = TRUE
 )
 }
 \arguments{
@@ -29,6 +30,9 @@ be set in  \code{run_swat_verification()}}
 
 \item{exclude_lum}{Character vector to define land uses which are excluded
 in the printed table.}
+
+\item{boxpoints}{Optional Boolean TRUE for displaying outliers, FALSE for hiding them.
+\code{Default = TRUE}}
 }
 \value{
 plotly figure object

man/print_avannual_qtile.Rd

 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/plot_hru_pw.R
+% Please edit documentation in R/plot_ps_tile.R
 \name{print_avannual_qtile}
 \alias{print_avannual_qtile}
 \title{Print the average annual qtile for HRUs}