R/f_read_output.R

In biopsichas/swattools: SWAT output assessment and visualising of modelling results

#' Function to read directly SWAT model output files from folder paths
#'
#'This function allows reading all output.rch, output.sub, output.hru files generated by SWAT2012
#'
#' @param folder String giving path to the output.*** folder. (example "data-raw/monthly")
#' @param starting_date Sting with time series starting date (example "1997-01-01")
#' @param ending_date Sting with time series ending date (example "2019-12-31")
#' @param output_type Number the same as in IPRINT (0 - Monthly, 1 - Daily, 2 - Year)
#' @param output_file String showing ending of output file. Currently implemented are "rch", "sub", "hru" files.
#' @return Data.frame imported from SWAT output files
#' @importFrom readr read_table2 col_character col_skip cols
#' @importFrom lubridate month year mday
#' @importFrom dplyr bind_cols arrange mutate filter %>%
#' @importFrom tidyr separate
#' @export

read_swat_output <- function(folder, output_type, output_file, starting_date, ending_date){

  ##Reading text files and writing column names for .rch, .sub and .hru SWAT output files
  if (output_file == "rch"){
    ##Reading in the output.rch data
    output_raw <- read_table2(paste0(folder, "/", "output.", output_file), skip = 9, col_types = cols(X3 = col_skip(), X4 = col_character()), col_names = FALSE) %>%
      mutate(X4 = suppressWarnings(as.numeric(X4)))
    ## column names for output.rch file
    col_names <- c("FILE","RCH","MON","AREAkm2","FLOW_INcms","FLOW_OUTcms","EVAPcms","TLOSScms",
                       "SED_INtons","SED_OUTtons","SEDCONCmg_kg","ORGN_INkg","ORGN_OUTkg","ORGP_INkg","ORGP_OUTkg",
                       "NO3_INkg","NO3_OUTkg","NH4_INkg","NH4_OUTkg","NO2_INkg","NO2_OUTkg","MINP_INkg","MINP_OUTkg",
                       "CHLA_INkg","CHLA_OUTkg","CBOD_INkg","CBOD_OUTkg","DISOX_INkg","DISOX_OUTkg","SOLPST_INmg",
                       "SOLPST_OUTmg","SORPST_INmg","SORPST_OUTmg","REACTPSTmg","VOLPSTmg","SETTLPSTmg","RESUSP_PSTmg",
                       "DIFFUSEPSTmg","REACBEDPSTmg","BURYPSTmg","BED_PSTmg","BACTP_OUTct","BACTLP_OUTct","CMETAL_1kg",
                       "CMETAL_2kg","CMETAL_3kg","TOTNkg","TOTPkg","NO3_mg_l","WTMPdegc", "RCH2")

  } else if(output_file == "sub"){
    ##Reading in the output.sub data
    output_raw <- read_table2(paste0(folder, "/", "output.", output_file), skip = 9, col_types = cols(X3 = col_skip(), X4 = col_character()), col_names = FALSE) %>%
      separate(3, into = c('M', 'A'), sep = "\\.") %>% ##One column have to be separated as is jointed
      mutate(A = as.numeric(paste0(".", A)),
             M = suppressWarnings(as.numeric(M)))
    ##Column names for output.sub file
    col_names <- c("SUB", "GIS", "MON", "AREAkm2", "PRECIPmm", "SNOMELTmm", "PETmm", "ETmm", "SWmm",
                       "PERCmm", "SURQmm", "GW_Qmm", "WYLDmm", "SYLDt/ha", "ORGNkg/ha", "ORGPkg/ha",
                       "NSURQkg/ha", "SOLPkg/ha", "SEDPkg/ha", "LATQmm", "LATNO3kg/ha", "GWNO3kg/ha",
                       "CHOLAmic/L", "CBODUmg/L", "DOXQmg/L", "TNO3kg/ha", "RCH")

  } else if(output_file == "hru"){
    ##Reading in the output.hru data
    output_raw <- read_table2(paste0(folder, "/", "output.", output_file), skip = 9, col_types = cols(X5 = col_skip(), X6 = col_character()), col_names = FALSE) %>%
      separate(5, into = c('M', 'A'), sep = "\\.") %>% ##One column have to be seperated as is jointed
      mutate(A = as.numeric(paste0(".", A)),
             M = suppressWarnings(as.numeric(M)))
    ##Column names for output.hru file
    col_names <- c("LULC", "HRU", "GIS", "SUB", "MON", "AREAkm2", "WYLDmm", "ORGNkg_ha", "ORGPkg_ha",
                       "SEDPkg_ha", "NSURQkg_ha", "NLATQkg_ha", "NO3GWkg_ha", "SOLPkg_ha", "P_GWkg_ha", "YLDt_ha", "TNO3kg_ha")
  } else {
    stop("Output file is not supported!!! Only 'rch', 'sub', 'hru' are available!!!" )
  }
  ##Assigning column names
  colnames(output_raw) <- col_names

  ##Generating daily time for all modelled period
  date_daily <- data.frame(date = rep(seq(from = as.Date(starting_date), to = as.Date(ending_date), by = 1),
                                      max(unique(output_raw[, 2])))) %>%
    arrange(date) %>%
    mutate(Year = year(date),
           Month = month(date),
           Day = mday(date))

  ##Dealing with output of different timescale
  if (output_type == 1){ ##For daily output
    output <- bind_cols(date_daily, output_raw) ##Binding dates to the dataframe

  } else if (output_type == 2){ ##For yearly output
    output <- bind_cols(filter(date_daily, Month == 1 & Day == 1), filter(output_raw, MON > 1900))

  } else if (output_type == 0){ ##For monthly output
    output <- bind_cols(filter(date_daily, Day == 15), filter(output_raw, MON <= 12))

  } else {
    stop("Output type could only be a numbers 0 for monthly, 1 for daily and 2 for yearly output!!!")
  }

  ##Returning result table
  return(output)
}