R/rvn_rvp_fill_template.R

In RavenR: Raven Hydrological Modelling Framework R Support and Analysis

#' @title Rewrite template rvp file with values
#'
#' @description
#' Rewrites a Raven template rvp file with default parameter values.
#'
#' @details
#' The Raven rvp template file is generated by Raven when the :CreateRVPTemplate command is included in the
#' rvi file (the default when \code{\link{rvn_rvi_write_template}} is used to produce an rvi file). This
#' template file displays the layout of the rvp file with required parameters based on the hydrologic processes
#' in the rvi file, but is not immediately usable. This function uses the soil model information in the rvi file
#' and the HRU class information in the rvh file to rewrite the template file with default parameter values so that
#' it can be used in a model run.
#'
#' If \code{rvp_out} is not provided, Raven will attempt to write to the file prefix of the provided template file with a .rvp extension.
#' If there is a conflict with an existing file and \code{overwrite==FALSE}, the function will automatically overwrite a file
#' with the suffix "_ravenr_generated.rvp".
#'
#' The default parameter values come from the RavenParameters.dat file included with RavenR in the extdata folder. The
#' user may provide their own file with updated values if preferred. Note that the database files held in the RavenR
#' package are unofficial copies of those in the official Raven SVN, and any discrepancies should defer to the Raven SVN versions.
#'
#' Any parameters not found in this file will be written
#' with the value provided by \code{default_param_value}. The default soil thickness for the :SoilProfiles block is provided
#' by the \code{default_soil_thickness} function parameter, which is applied for all soil classes.
#'
#' As an alternative to specifying the three input files (rvi, rvh, rvp_temp.rvp), the \code{fileprefix} of the model (e.g.
#' 'Nith' for Nith.rvi) may be provided instead. If provided, the function will attempt to find all required input files based on
#' the provided fileprefix in the current working directory.
#'
#' Additional commands can be added to the end of the rvp file with \code{extra_commands}, which are not quality controlled but
#' simply appended to the rvp file. This can be useful for non-standard commands such as :RedirectToFile for channel properties rvp files
#' that are not added automatically to base templates rvp files.
#'
#' If you find parameters not found by this function, please open an ticket on Github (\url{https://github.com/rchlumsk/RavenR/issues}).
#'
#' @param rvi_file path to the model *.rvi file
#' @param rvh_file path to the model *.rvh file
#' @param rvp_template_file path to the model rvp template file (*.rvp_temp.rvp)
#' @param fileprefix model name prefix for the main model files; if provided, the function will attempt to find all missing file paths in the current working directory
#' @param rvp_out file path to rewritten rvp file
#' @param overwrite logical for whether to rewrite the *.rvp file if it already exists (default FALSE)
#' @param default_param_value default parameter value to write for any parameters not found in RavenParameters.dat
#' @param default_soil_thickness default soil layer thickness (m) to provide in :SoilProfile block
#' @param use_default_flag writes all soil/land use/vegetation classes with [DEFAULT] flag
#' @param avg_annual_runoff adds a line for :AvgAnnualRunoff if value provided with this parameter and not already in file
#' @param extra_commands additional commands to add to end of rvp file as character vector
#' @param RavenParamsFile path to RavenParameters.dat file (default path points to file included with RavenR installation)
#' @return \code{TRUE} if the function executed successfully
#'
#' @seealso \code{\link{rvn_rvi_getparams}} to get parameter ranges from rvi.
#'
#' @examples
#'
#' ### this section is not run, but illustrates how an rvp template file would be created
#' # -----
#' \dontrun{
#' ## create an rvi file and template file with Raven
#' rvn_rvi_write_template(modelname="HBV-EC",
#'    filename="nithmodel.rvi")
#'
#' ## download Raven.exe if not already downloaded
#' if (!rvn_download(check=TRUE)) {
#'   rvn_download()
#' }
#'
#' ## run Raven to create template file
#' rvn_run(fileprefix="nithmodel")
#' }
#' # -----
#'
#' # load pre-generated template file and other model files
#' nithmodel_template_file <- system.file("extdata","nithmodel.rvp_temp.rvp", package="RavenR")
#' nith_rvi_file <- system.file("extdata","Nith.rvi", package="RavenR")
#' nith_rvh_file <- system.file("extdata","Nith.rvh", package="RavenR")
#' rvp_out_file <- file.path(tempdir(), 'nithmodel.rvp')
#'
#' # rewrite template with parameter values
#' rvn_rvp_fill_template(rvi_file=nith_rvi_file,
#'                       rvh_file=nith_rvh_file,
#'                       rvp_template_file=nithmodel_template_file,
#'                       rvp_out=rvp_out_file)
#'
#' @importFrom dplyr left_join
#' @export rvn_rvp_fill_template
rvn_rvp_fill_template <- function(rvi_file=NULL, rvh_file=NULL, rvp_template_file=NULL, fileprefix=NULL,
                                  rvp_out=NULL, overwrite=FALSE,
                                  default_param_value=0.12345, default_soil_thickness=0.5,
                                  use_default_flag=FALSE, avg_annual_runoff=NULL,
                                  extra_commands=NULL,
                                  RavenParamsFile=system.file("extdata","RavenParameters.dat", package="RavenR")) {

  ## attempt to find all required files if fileprefix provided
  if (!is.null(fileprefix)) {
    if (is.null(rvi_file)) {rvi_file <- paste0(fileprefix,".rvi")}
    if (is.null(rvh_file)) {rvh_file <- paste0(fileprefix,".rvh")}
    if (is.null(rvp_template_file)) {rvp_template_file <- paste0(fileprefix,".rvp_temp.rvp")}
  }

  ## check provided file inputs
  if (!is.character(rvi_file) & file.exists(rvi_file)) {
    stop("rvi_file must be a valid file path")
  }

  if (rvn_substrRight(rvi_file, 4) != ".rvi") {
    stop("rvi_file must be a valid Raven rvi file (*.rvi)")
  }

  if (!is.character(rvh_file) & file.exists(rvh_file)) {
    stop("rvh_file must be a valid file path")
  }

  if (rvn_substrRight(rvh_file, 4) != ".rvh") {
    stop("rvh_file must be a valid Raven rvh file (*.rvh)")
  }

  if (rvn_substrRight(rvp_template_file, 13) != ".rvp_temp.rvp") {
    stop("rvp_template must be a valid Raven template file (.rvp_temp.rvp)")
  }

  if (!is.character(rvp_template_file) & file.exists(rvp_template_file)) {
    stop("rvp_template must be a valid file path")
  }

  ## read in files
  rvi <- rvn_rvi_read(rvi_file)
  rvh <- rvn_rvh_read(rvh_file)
  tt <- readLines(rvp_template_file)

  ## get number of soil layers
  rvi_lines <- readLines(rvi_file)
  soil_model_string <- rvi_lines[grep(pattern=":SoilModel",rvi_lines)] %>%
    trimws() %>%
    strsplit(split="\\s+") %>%
    unlist()

  if (length(soil_model_string) ==3 ) {
    num_soil_classes <- as.numeric(soil_model_string[3])
  } else if (soil_model_string[2] == "SOIL_ONE_LAYER") {
    num_soil_classes <- 1
  } else if (soil_model_string[2] == "SOIL_TWO_LAYER") {
    num_soil_classes <- 2
  } else {
    stop("Unrecognized :SoilModel command.")
  }

  ## set rvp_out
  if (is.null(rvp_out)) {
    rvp_out <- rvn_substrMRight(rvp_template_file, 9)
  }

  ## check if rvp_out exists and overwrite is not enabled
  if (file.exists(rvp_out) & !overwrite) {
    rvp_out <- sprintf("%s_ravenr_generated.rvp",rvn_substrMRight(rvp_out,4))
    overwrite <- TRUE ## overwrite RavenR file is present
  }

  ## read rvh information (profiles)
  soil_classes <- sprintf("SOIL_%02d",seq(1,num_soil_classes))
  soil_profiles <- unique(rvh$HRUtable$SoilProfile)
  land_classes <- unique(rvh$HRUtable$LandUse)
  veg_classes  <- unique(rvh$HRUtable$Vegetation)
  terrain_classes <- unique(rvh$HRUtable$Terrain)

  terrain_classes <-  gsub("[NONE]", replacement=NA, x=terrain_classes)
  ## process classes to remove NA values?


  ## read in required .dat tables
  if (!file.exists(RavenParamsFile)) {
    stop(sprintf("Provided RavenParamsFile does not exist, please check: %s", RavenParamsFile))
  }
  cnames <- c("param","class_type","units","auto","default","min","max")
  RavenParamsTable<-read.table(RavenParamsFile,
                               sep="",
                               col.names=cnames,
                               header=FALSE,
                               blank.lines.skip=TRUE,
                               strip.white=TRUE,
                               stringsAsFactors=FALSE,
                               flush=TRUE,
                               comment.char = "#")

  #### start writing file ----
  fc <- file(rvp_out, open="w+")

  #### write header ----
  writeLines(text=tt[1:2], fc)
  writeLines(sprintf("# Modified by RavenR::rvn_rvp_fill_template (%s)",Sys.Date()),fc)
  writeLines(text=tt[7:8], fc)

  #### write soil classes ----
  # start <- grep(pattern=":SoilClasses",tt)
  # end <- grep(pattern=":EndSoilClasses",tt)
  writeLines("#-----------------------------------------------------------------",fc)
  writeLines("# Soil Classes",fc)
  writeLines("#-----------------------------------------------------------------",fc)
  writeLines(":SoilClasses",fc)
  writeLines("  :Attributes,",fc)
  writeLines("  :Units,",fc)
  writeLines(sprintf("  %s,",soil_classes),fc)
  writeLines(":EndSoilClasses",fc)
  writeLines("",fc)


  #### write land use classes ----
  start <- grep(pattern=":LandUseClasses",tt)
  end <- grep(pattern=":EndLandUseClasses",tt)

  writeLines("#-----------------------------------------------------------------",fc)
  writeLines("# Land Use Classes",fc)
  writeLines("#-----------------------------------------------------------------",fc)
  writeLines(":LandUseClasses",fc)
  writeLines(tt[c(start+1,start+2)],fc) # write attributes and units

  ## get attributes
  temp <- unlist(strsplit(trimws(tt[(start+1)]),split="\\s+"))
  atts <- temp[2:length(temp)]
  atts <- gsub(",","",atts)
  atts <- data.frame("param"=atts, "value"=rep(NA,length(atts)))

  # bring in default values from RavenParamsTable
  atts <- left_join(atts,
                     RavenParamsTable[which(RavenParamsTable$param %in% atts$param), c("param","default")],
                     by="param")
  atts[!is.na(atts$default),]$value <- atts[!is.na(atts$default),]$default

  if (any(is.na(atts$value))) {
    warning(sprintf("Some land use class parameters from rvp file not found in RavenParameters.dat:\n%s",
            paste0(atts[is.na(atts$default),]$param, collapse="\n")))
    atts[is.na(atts$default),]$value <- default_param_value
  }
  if (any(duplicated(atts[,c('param')]))) {
    warning("Some land use class parameters are duplicated in RavenParameters.dat, only the first will be used:\n",
            paste0(atts$param[duplicated(atts[,c('param')])], collapse="\n"))
  }
  atts <- atts[!duplicated(atts[,c('param')]),]

  ## write attribute values for each land use class
  for (i in 1:length(land_classes)) {
    writeLines(paste0(sprintf("%014s,",c(land_classes[i], atts$value)),collapse=" "),fc)
  }
  writeLines(":EndLandUseClasses",fc)
  writeLines("",fc)

  #### write vegetation classes ----
  start <- grep(pattern=":VegetationClasses",tt)
  end <- grep(pattern=":EndVegetationClasses",tt)

  writeLines("#-----------------------------------------------------------------",fc)
  writeLines("# Vegetation Classes",fc)
  writeLines("#-----------------------------------------------------------------",fc)
  writeLines(":VegetationClasses",fc)
  writeLines(tt[c(start+1,start+2)],fc) # write attributes and units

  ## get attributes
  temp <- unlist(strsplit(trimws(tt[(start+1)]),split="\\s+"))
  atts <- temp[2:length(temp)]
  atts <- gsub(",","",atts)
  atts <- data.frame("param"=atts, "value"=rep(NA,length(atts)))

  # bring in default values from RavenParamsTable
  atts <- left_join(atts,
                     RavenParamsTable[which(RavenParamsTable$param %in% atts$param), c("param","default")],
                    by="param")
  atts[!is.na(atts$default),]$value <- atts[!is.na(atts$default),]$default

  if (any(is.na(atts$value))) {
    warning(sprintf("Some vegetation class parameters from rvp file not found in RavenParameters.dat:\n%s",
                    paste0(atts[is.na(atts$default),]$param, collapse="\n")))
    atts[is.na(atts$default),]$value <- default_param_value
  }
  if (any(duplicated(atts[,c('param')]))) {
    warning("Some vegetation class parameters are duplicated in RavenParameters.dat, only the first will be used:\n",
            paste0(atts$param[duplicated(atts[,c('param')])], collapse="\n"))
  }
  atts <- atts[!duplicated(atts[,c('param')]),]

  ## write attribute values for each vegetation class
  for (i in 1:length(veg_classes)) {
    writeLines(paste0(sprintf("  %014s,",c(veg_classes[i], atts$value)),collapse=" "),fc)
  }
  writeLines(":EndVegetationClasses",fc)
  writeLines("",fc)

  ## add in terrain classes (if terrain classes not in NA > 1)


  #### write soil profiles ----
  start <- grep(pattern=":SoilProfiles",tt)
  end <- grep(pattern=":EndSoilProfiles",tt)
  profile_1 <- grep(pattern="PROFILE_1",tt)

  temp <- unlist(strsplit(trimws(tt[(start+1):(profile_1-1)]),split=","))
  existing_profiles <- temp[seq(1,length(temp),by=2)]

  ## misc check
  if (profile_1 < start | profile_1 > end) {
    warning("Search of file for PROFILE_1 in SoilProfiles declaration found inconsistent match. Check results.")
  }

  writeLines("#-----------------------------------------------------------------",fc)
  writeLines("# Soil Profiles",fc)
  writeLines("#-----------------------------------------------------------------",fc)
  writeLines(tt[(start):(profile_1-1)],fc) # write soil profile lines until first profile

  ## write default soil thickness for each layer
  for (i in 1:length(soil_profiles)) {
    if (soil_profiles[i] %notin% existing_profiles) {
      ss <- sprintf("  %014s,%08s,",soil_profiles[i],num_soil_classes)
      for (k in 1:num_soil_classes) {
        ss <- paste0(ss, sprintf("%014s,%014s,", soil_classes[k], default_soil_thickness ))
      }
      writeLines(ss, fc)
    }
  }
  writeLines(":EndSoilProfiles",fc)
  writeLines("",fc)


  #### write global parameters  ----
  global_param_indices <- grep(pattern=":GlobalParameter",tt)

  writeLines("#-----------------------------------------------------------------",fc)
  writeLines("# Global Parameters",fc)
  writeLines("#-----------------------------------------------------------------",fc)

  if (length(global_param_indices) == 0) {
    warning("No global parameters found in %s",rvp_template_file)
  } else {
    for (i in 1:length(global_param_indices)) {

      temp <- unlist(strsplit(trimws(tt[global_param_indices[i]]),split="\\s+"))

      if (length(temp) != 3) {
        warning(sprintf("Error reading global parameter line below:\n%s",tt[global_param_indices[i]]))
      } else {
        param <- temp[2]
        atts <- data.frame("param"=param, "value"=rep(NA,length(param)))

        # bring in default values from RavenParamsTable
        atts <- left_join(atts,
                           RavenParamsTable[which(RavenParamsTable$param %in% atts$param), c("param","default")],
                          by="param")
        atts[!is.na(atts$default),]$value <- atts[!is.na(atts$default),]$default

        if (any(is.na(atts$value))) {
          warning(sprintf("Some global parameters parameters from rvp file not found in RavenParameters.dat:\n%s",
                          paste0(atts[is.na(atts$default),]$param, collapse="\n")))
          atts[is.na(atts$default),]$value <- default_param_value
        }

        if (nrow(atts)>1) {
          warning(sprintf("For param %s, multiple default values found in Raven .dat file. Only the first will be written to file:\n%s",
                          param,
                          paste0(atts$value, collapse=" ")))
          atts <- atts[1,]
        }

        writeLines(sprintf(":GlobalParameter %014s %014s",param, atts$value),fc)
      }
    }
  }
  writeLines("",fc)


  #### write soil parameters  ----
  start <- grep(pattern=":SoilParameterList",tt)
  end <- grep(pattern=":EndSoilParameterList",tt)

  writeLines("#-----------------------------------------------------------------",fc)
  writeLines("# Soil Parameters",fc)
  writeLines("#-----------------------------------------------------------------",fc)
  writeLines(":SoilParameterList",fc)
  writeLines(tt[c(start+1,start+2)],fc) # write attributes and units

  ## get attributes
  temp <- unlist(strsplit(trimws(tt[(start+1)]),split="\\s+"))
  atts <- temp[2:length(temp)]
  atts <- gsub(",","",atts)
  atts <- data.frame("param"=atts, "value"=rep(NA,length(atts)))

  # bring in default values from RavenParamsTable
  atts <- left_join(atts,
                     RavenParamsTable[which(RavenParamsTable$param %in% atts$param), c("param","default")],
                    by="param")
  atts[!is.na(atts$default),]$value <- atts[!is.na(atts$default),]$default

  if (any(is.na(atts$value))) {
    warning(sprintf("Some soil parameters from rvp file not found in RavenParameters.dat:\n%s",
                    paste0(atts[is.na(atts$default),]$param, collapse="\n")))
    atts[is.na(atts$default),]$value <- default_param_value
  }
  if (any(duplicated(atts[,c('param')]))) {
    warning("Some soil parameters are duplicated in RavenParameters.dat, only the first will be used:\n",
            paste0(atts$param[duplicated(atts[,c('param')])], collapse="\n"))
  }
  atts <- atts[!duplicated(atts[,c('param')]),]

  ## write attribute values for each soil class
  if (use_default_flag) {
    writeLines(paste0(sprintf("  %014s,",c("[DEFAULT]", atts$value)),collapse=" "),fc)
    for (i in 1:length(soil_classes)) {
      writeLines(paste0(sprintf("#  %014s,",c(soil_classes[i], atts$value)),collapse=" "),fc)
    }

  } else {
    for (i in 1:length(soil_classes)) {
      writeLines(paste0(sprintf("  %014s,",c(soil_classes[i], atts$value)),collapse=" "),fc)
    }
  }

  writeLines(":EndSoilParameterList",fc)
  writeLines("",fc)


  #### write land use parameters  ----
  start <- grep(pattern=":LandUseParameterList",tt)
  end <- grep(pattern=":EndLandUseParameterList",tt)

  writeLines("#-----------------------------------------------------------------",fc)
  writeLines("# Land Use Parameters",fc)
  writeLines("#-----------------------------------------------------------------",fc)
  writeLines(":LandUseParameterList",fc)
  writeLines(tt[c(start+1,start+2)],fc) # write attributes and units

  ## get attributes
  temp <- unlist(strsplit(trimws(tt[(start+1)]),split="\\s+"))
  atts <- temp[2:length(temp)]
  atts <- gsub(",","",atts)
  atts <- data.frame("param"=atts, "value"=rep(NA,length(atts)))

  # bring in default values from RavenParamsTable
  atts <- left_join(atts,
             RavenParamsTable[which(RavenParamsTable$param %in% atts$param), c("param","default")],
             by="param")
  atts[!is.na(atts$default),]$value <- atts[!is.na(atts$default),]$default

  if (any(is.na(atts$value))) {
    warning(sprintf("Some land use parameters from rvp file not found in RavenParameters.dat:\n%s",
                    paste0(atts[is.na(atts$default),]$param, collapse="\n")))
    atts[is.na(atts$default),]$value <- default_param_value
  }
  if (any(duplicated(atts[,c('param')]))) {
    warning("Some land use parameters are duplicated in RavenParameters.dat, only the first will be used:\n",
            paste0(atts$param[duplicated(atts[,c('param')])], collapse="\n"))
  }
  atts <- atts[!duplicated(atts[,c('param')]),]

  ## write attribute values for each land use class
  if (use_default_flag) {
    writeLines(paste0(sprintf("  %014s,",c("[DEFAULT]", atts$value)),collapse=" "),fc)
    for (i in 1:length(land_classes)) {
      writeLines(paste0(sprintf("#  %014s,",c(land_classes[i], atts$value)),collapse=" "),fc)
    }
  } else {
    for (i in 1:length(land_classes)) {
      writeLines(paste0(sprintf("  %014s,",c(land_classes[i], atts$value)),collapse=" "),fc)
    }
  }

  writeLines(":EndLandUseParameterList",fc)
  writeLines("",fc)


  #### write vegetation parameters  ----
  start <- grep(pattern=":VegetationParameterList",tt)
  end <- grep(pattern=":EndVegetationParameterList",tt)

  writeLines("#-----------------------------------------------------------------",fc)
  writeLines("# Vegetation Parameters",fc)
  writeLines("#-----------------------------------------------------------------",fc)
  writeLines(":VegetationParameterList",fc)
  writeLines(tt[c(start+1,start+2)],fc) # write attributes and units

  ## get attributes
  temp <- unlist(strsplit(trimws(tt[(start+1)]),split="\\s+"))
  atts <- temp[2:length(temp)]
  atts <- gsub(",","",atts)
  atts <- data.frame("param"=atts, "value"=rep(NA,length(atts)))

  # bring in default values from RavenParamsTable
  atts <- left_join(atts,
                    RavenParamsTable[which(RavenParamsTable$param %in% atts$param), c("param","default")],
                    by="param")
  atts[!is.na(atts$default),]$value <- atts[!is.na(atts$default),]$default

  if (any(is.na(atts$value))) {
    warning(sprintf("Some vegetation parameters from rvp file not found in RavenParameters.dat:\n%s",
                    paste0(atts[is.na(atts$default),]$param, collapse="\n")))
    atts[is.na(atts$default),]$value <- default_param_value
  }
  if (any(duplicated(atts[,c('param')]))) {
    warning("Some vegetation parameters are duplicated in RavenParameters.dat, only the first will be used:\n",
            paste0(atts$param[duplicated(atts[,c('param')])], collapse="\n"))
  }
  atts <- atts[!duplicated(atts[,c('param')]),]

  if (any(is.na(atts$value))) {
    warning(sprintf("Some vegetation parameters from rvp file not found in RavenParameters.dat:\n%s",
                    paste0(atts[is.na(atts$default),]$param, collapse="\n")))
    atts[is.na(atts$default),]$value <- default_param_value
  }
  atts <- atts[!duplicated(atts[,c('param')]),]

  ## write attribute values for each vegetation class
  if (use_default_flag) {
    writeLines(paste0(sprintf("  %014s,",c("[DEFAULT]", atts$value)),collapse=" "),fc)
    for (i in 1:length(veg_classes)) {
      writeLines(paste0(sprintf("#  %014s,",c(veg_classes[i], atts$value)),collapse=" "),fc)
    }
  } else {
    for (i in 1:length(veg_classes)) {
      writeLines(paste0(sprintf("  %014s,",c(veg_classes[i], atts$value)),collapse=" "),fc)
    }
  }

  writeLines(":EndVegetationParameterList",fc)
  writeLines("",fc)


  ### add in terrain parameters


  #### write seasonal parameters  ----
  start <- grep(pattern=":SeasonalRelativeLAI",tt)
  end <- grep(pattern=":EndSeasonalRelativeLAI",tt)

  if (length(start) != 0) {
    writeLines(":SeasonalRelativeLAI",fc)
    for (i in 1:length(veg_classes)) {
      writeLines(sprintf("  %s,  %s",
                         veg_classes[i],
                         paste0(sprintf(" %.1f,",(rep(1.0,12))),
                                collapse=" ")),fc)
    }
    writeLines(":EndSeasonalRelativeLAI",fc)
    writeLines("",fc)
  }

  start <- grep(pattern=":SeasonalRelativeHeight",tt)
  end <- grep(pattern=":EndSeasonalRelativeHeight",tt)

  if (length(start) != 0) {
    writeLines(":SeasonalRelativeHeight",fc)
    for (i in 1:length(veg_classes)) {
      writeLines(sprintf("  %s,  %s",
                         veg_classes[i],
                         paste0(sprintf(" %.1f,",(rep(1.0,12))),
                                collapse=" ")),fc)
    }
    writeLines(":EndSeasonalRelativeHeight",fc)
    writeLines("",fc)
  }

  ## write :AvgAnnualRunoff if provided
  if (!is.null(avg_annual_runoff)) {
    start <- grep(pattern=":AvgAnnualRunoff",tt)
    if (length(start) != 0) {
      warning("Trying to add :AvgAnnualRunoff but already found in file, will not be written to file with new value")
    } else if (avg_annual_runoff < 0 | !inherits(avg_annual_runoff,"numeric")) {
      stop("avg_annual_runoff must be a positive double of type numeric.")
    } else {
      writeLines(sprintf(":AvgAnnualRunoff  %.2f",avg_annual_runoff),fc)
      writeLines("",fc)
    }
  }

  ## write extra commands to file
  if (!is.null(extra_commands)) {
    for (i in 1:length(extra_commands)) {
      writeLines(sprintf("%s",extra_commands[i]),fc)
      writeLines("",fc)
    }
  }

  ##### close out file ----
  close(fc)

  return(TRUE)
}