Nothing
# File hydromod.R
# Part of the hydroPSO R package, https://github.com/hzambran/hydroPSO
# http://cran.r-project.org/web/packages/hydroPSO
# http://www.rforge.net/hydroPSO/
# Copyright 2010-2020 Mauricio Zambrano-Bigiarini & Rodrigo Rojas
# Distributed under GPL 2 or later
################################################################################
# hydromod #
################################################################################
# Purpose : To run a hydrological/environmental model, and get a #
# goodness-of-fit value by comparing the simulated values against #
# observations #
################################################################################
# Output : A list of two elements: #
# 1) sim: simulated values obtained by running the hydrological #
# model #
# 2) GoF: goodness-of fit of the simualted values against observed#
# ones, by using THE USER-DEFINED 'gof' measure #
################################################################################
# Author : Mauricio Zambrano-Bigiarini #
# Started: 14-Dec-2010 at JRC Ispra #
# Updates: 20-Dec-2010 #
# 19-Jan-2011 ; 22-Jan-2011 ; 02-Feb-2011 ; 11-May-2011 #
# 13-Jan-2012 ; 16-Jan-2012 ; 23-Jan-2012 ; 02-May-2012 ; 12-Oct-2012 #
# 15-Oct-2012 ; 22-Nov-2012 #
# 29-May-2013 #
# 10-Jun-2018 #
# 09-Mar-2020 #
################################################################################
hydromod <- function(
param.values, # Numeric vector with the paramter values that will be used in the input files of the hydrological model
param.files="ParamFiles.txt", # Character, with the name of the file (with full path) that stores the name of the files that have to be modified for each parameter
model.drty=getwd(), # Character, with the path of the directory that stores the exe file of the hydrological model and ALL the input files required for the simulation
exe.fname,
stdout=FALSE, # a logical (not NA) indicating whether messages written to 'stdout' should be sent or not. See '?system2'
stderr="", # a logical (not NA) indicating whether messages written to 'stderr' should be sent or not. See '?system2'
verbose= FALSE, # logical, indicating if progress messages have to be printed during the simulations. If \code{verbose=TRUE}, the following messages will appear: i)parameter values for each particle; (ii) model execution; iii) extraction of simulated values; and iv) computation of the goodness-of-fit measure.
###############################
###############################
# function for reading the ouput
# of the hydrological model #
###############################
out.FUN, # Character, with the name of a valid R function to be used for reading the outputs of te model in 'out.fname' and transforming them into a zoo object
out.FUN.args, # list, with the arguments to be passed to 'out.FUN'
###############################
###############################
# function for comparing simulated
# values with observations #
###############################
gof.FUN, # name of the GoF function used to compare simulations and observations
gof.FUN.args=list(), # list, with additional arguments to be passed to 'out.FUN' (additional to 'sim' and 'obs')
gof.Ini, # OPTIONAL. Starting date for the the goodness-of-fit function to be used for comparing observed against simulated values.
# It is used to subset 'obs' (if necessary), AND to define the time period to compare simulated ith observed values
gof.Fin, # OPTIONAL. Ending date for the the goodness-of-fit function to be used for comparing observed against simulated values.
# It is used to subset 'obs' (if necessary), AND to define the time period to compare simulated ith observed values
date.fmt="%Y-%m-%d", # Format used to define 'gof.Ini", "gof.Fin"
obs, # zoo object with the observed values
##############################
# OPTIONAL. for plotting #
# sim vs obs #
do.png=FALSE, # logical, indicating if a png image with the results of the 'ggof' function has to be produced
png.fname, # OPTIONAL, only used when \code{do.png=TRUE}. Name of the PNG file to be produced. By default its name is 'Obs_vs_Sim.png', within the 'model.drty'
width = 1200, # OPTIONAL, only used when \code{do.png=TRUE}. Widht of the output PNG image
height = 600, # OPTIONAL, only used when \code{do.png=TRUE}. Height of the output PNG image
res=90, # OPTIONAL, only used when \code{do.png=TRUE}. Resolution of the output PNG image
main, # OPTIONAL, only used when \code{do.png=TRUE}. character, representing the main title of the plot comparing observed and simulated values
leg.cex=1.2,
tick.tstep= "auto", # character, indicating the time step that have to be used for putting the ticks on the time axis. See 'ggof'
lab.tstep= "auto", #character indicating the time step that have to be used for putting the labels on the time axis. See 'ggof'
lab.fmt=NULL # Character, with the format to be used for drawing the time axis.
###############################
#... # Addtional arguments to be passed to
) {
##############################################################################
# 0) Checkings #
##############################################################################
# Verifying 'param.values'
#if (missing(param.values))
# stop( "Missing argument: 'param.values' has to be given !")
if (!file.exists(param.files))
stop( "Invalid argument: the file '", param.files, "' doesn't exist!" )
if ( missing(exe.fname) )
stop( "Missing argument: 'exe.fname'" )
if ( missing(out.FUN) ) {
stop( "Missing argument: 'out.FUN'" )
} else out.FUN <- match.fun(out.FUN)
if ( missing(gof.FUN) ) {
stop( "Missing argument: 'gof.FUN'" )
} else {gof.name <- gof.FUN
gof.FUN <- match.fun(gof.FUN)
} # ELSE end
setwd(model.drty)
if (!file.exists(exe.fname))
stop( "Invalid argument: the file '", exe.fname, "' does not exist!" )
if ( sessionInfo()[[1]]$os != "linux-gnu") {
dot.pos <- which(strsplit(exe.fname, split=character(0))[[1]] == ".")
ext <- substr(exe.fname,start=dot.pos+1, stop=nchar(exe.fname))
if (ext=="exe") exe.fname <- substr(exe.fname,start=1, stop=dot.pos-1)
} # IF end
if (missing(obs)) stop("Missing argument: 'obs' must be provided")
##############################################################################
# 1) New parameter values -> input files of the hydrological model #
##############################################################################
if (verbose) message(" ")
if (verbose) message("===========================================")
if (verbose) message("[ 1) Writing new parameter values ... ]")
if (verbose) message("===========================================")
ParameterValues2InputFiles(NewValues=param.values, ParamFiles.fname=param.files,
verbose=verbose)
##############################################################################
# 2) Running the hydrological model #
##############################################################################
if (verbose) message("===========================================")
if (verbose) message("[ 2) Running the model ... ]")
if (verbose) message("===========================================")
system2(exe.fname, stdout=stdout, stderr=stderr)
##############################################################################
# 3) Extracting simulated values #
##############################################################################
if (verbose) message("===========================================")
if (verbose) message("[ 3) Extracting simulated values ... ]")
if (verbose) message("===========================================")
out.FUN.argsDefaults <- formals(out.FUN)
out.FUN.args <- modifyList(out.FUN.argsDefaults, out.FUN.args)
sim <- do.call(out.FUN, as.list(out.FUN.args))
##############################################################################
# 4) Goodness of fit #
##############################################################################
if (verbose) message("===========================================")
if (verbose) message("[ 4) Computing the goodness of fit ... ]")
if (verbose) message("===========================================")
# In case different dates are given for 'sim', 'obs', 'gof.Ini', gof.Fin',
# 'sim' and 'obs' are subset to the time period selected for the GoF function
if ( !missing(gof.Ini) | !missing(gof.Fin) )
ifelse ( grepl("%H", date.fmt, fixed=TRUE) | grepl("%M", date.fmt, fixed=TRUE) |
grepl("%S", date.fmt, fixed=TRUE) | grepl("%I", date.fmt, fixed=TRUE) |
grepl("%p", date.fmt, fixed=TRUE) | grepl("%X", date.fmt, fixed=TRUE),
subdaily <- TRUE, subdaily <- FALSE )
if (!missing(gof.Ini) ) {
ifelse(subdaily, gof.Ini <- as.POSIXct(gof.Ini, format=date.fmt),
gof.Ini <- as.Date(gof.Ini, format=date.fmt) )
if (!is.zoo(obs)) { # zoo::is.zoo
stop( "Invalid argument: 'obs' must be a zoo or xts object to use 'gof.Ini' !")
} else obs <- window(obs, start=gof.Ini)
if (!is.zoo(sim)) { # zoo::is.zoo
stop( "Invalid argument: 'sim' must be a zoo or xts object to use 'gof.Ini' !")
} else sim <- window(sim, start=gof.Ini)
} # IF end
if (!missing(gof.Fin) ) {
ifelse(subdaily, gof.Fin <- as.POSIXct(gof.Fin, format=date.fmt),
gof.Fin <- as.Date(gof.Fin, format=date.fmt) )
if (gof.Fin < gof.Ini) {
stop( "Invalid argument: 'gof.Fin < gof.Ini' (", gof.Fin, " < ", gof.Ini, ")" )
} else {
if (!is.zoo(obs)) { # zoo::is.zoo
stop( "Invalid argument: 'obs' must be a zoo or xts object to use 'gof.Fin' !" )
} else obs <- window(obs, end=gof.Fin)
if (!is.zoo(sim)) { # zoo::is.zoo
stop( "Invalid argument: 'sim' must be a zoo or xts object to use 'gof.Fin' !" )
} else sim <- window(sim, end=gof.Fin)
} # IF end
} # IF end
nobs <- length(obs)
nsim <- length(sim)
if (nobs != nsim)
stop( "Invalid argument: number of observations != number of simulations (", nobs, "!=", nsim, ")'" )
gof.FUN.argsDefaults <- formals(gof.FUN)
gof.FUN.args <- modifyList(gof.FUN.argsDefaults, gof.FUN.args)
gof.FUN.args <- modifyList(gof.FUN.args, list(sim=sim, obs=obs))
gof.value <- do.call(gof.FUN, as.list(gof.FUN.args))
if (verbose) message("[", gof.name, "= ", round(gof.value,3), "]")
##############################################################################
# 5) Creating the output object #
##############################################################################
out <- list(2)
out[[1]] <- sim
out[[2]] <- gof.value
names(out) <- c("sim", "GoF")
##############################################################################
# 6) OPTIONAL: GoF plot #
##############################################################################
if (do.png) {
if (verbose) message("===========================================")
if (verbose) message("[ 5) Creating a PNG with 'sim' vs 'obs'...]")
if (verbose) message("===========================================")
if (missing(main)) main <- "Simulations vs Observations"
if (missing(png.fname))
png.fname <- paste(file.path(model.drty), "/Obs_vs_Sim.png", sep="")
png(filename=png.fname, width= width, height= height, res=res)
if ( length(find.package("hydroGOF", quiet=TRUE)) != 0 ) {
hydroGOF::ggof(sim, obs, main=main, cex.main=1.5, leg.cex=leg.cex, tick.tstep=tick.tstep, lab.tstep=lab.tstep, lab.fmt=lab.fmt)
} else plot_out(sim=sim, obs=obs, ptype="corr")
dev.off()
if (verbose) message("===========================================")
if (verbose) message("[ See the file '", basename(png.fname), "' ")
if (verbose) message("===========================================")
} # IF end
return(out)
} # 'hydromod' END
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.