R/calculate_stats.R
In mccreigh/rwrfhydro: R tools for the WRF Hydro Model
Defines functions
CalcNoahmpWatBudg
CalcNoahmpFluxes
Documented in
CalcNoahmpFluxes
CalcNoahmpWatBudg
#' Calculate water fluxes from NoahMP output
#'
#' \code{CalcNoahmpFluxes} calculates water balance fluxes from accumulated
#' water terms.
#'
#' Read a dataframe derived from NoahMP LDASOUT output (i.e., using
#' \code{\link{GetMultiNcdf}}) and calculate water budget component fluxes from
#' accumulated water variables. NOTE: Currently only works for runs using NoahMP
#' as the LSM.
#'
#' @param ldasoutDf The LDASOUT dataframe
#' @param idCol The optional ID column to use to parse the output;
#' for example, a basin ID column or other index ID. (DEFAULT=NULL)
#' @return The input dataframe with new water flux columns added.
#'
#' @examples
#' ## Take a NoahMP LDASOUT dataframe for a model run of Fourmile Creek
#' ## and generate a dataframe with water fluxes added.
#' \dontrun{
#' modLDASOUT1d.nort.fc <- CalcNoahmpFluxes(modLDASOUT1d.nort.fc)
#' }
#' @keywords manip
#' @concept dataMgmt
#' @family modelEvaluation
#' @export
CalcNoahmpFluxes <- function(ldasoutDf, idCol=NULL) {
ldasoutDf <- ldasoutDf[order(ldasoutDf[,"POSIXct"]),]
if (is.null(idCol)) {
if ("ACCPRCP" %in% colnames(ldasoutDf)) {
ldasoutDf$DEL_ACCPRCP[2:nrow(ldasoutDf)] <- diff(ldasoutDf$ACCPRCP) }
if ("ACCECAN" %in% colnames(ldasoutDf)) {
ldasoutDf$DEL_ACCECAN[2:nrow(ldasoutDf)] <- diff(ldasoutDf$ACCECAN) }
if ("ACCETRAN" %in% colnames(ldasoutDf)) {
ldasoutDf$DEL_ACCETRAN[2:nrow(ldasoutDf)] <- diff(ldasoutDf$ACCETRAN) }
if ("ACCEDIR" %in% colnames(ldasoutDf)) {
ldasoutDf$DEL_ACCEDIR[2:nrow(ldasoutDf)] <- diff(ldasoutDf$ACCEDIR) }
if ("ACCET" %in% colnames(ldasoutDf)) {
ldasoutDf$DEL_ACCET[2:nrow(ldasoutDf)] <- diff(ldasoutDf$ACCET) }
if ("UGDRNOFF" %in% colnames(ldasoutDf)) {
ldasoutDf$DEL_UGDRNOFF[2:nrow(ldasoutDf)] <- diff(ldasoutDf$UGDRNOFF) }
if ("SFCRNOFF" %in% colnames(ldasoutDf)) {
ldasoutDf$DEL_SFCRNOFF[2:nrow(ldasoutDf)] <- diff(ldasoutDf$SFCRNOFF) }
if ("ACSNOM" %in% colnames(ldasoutDf)) {
ldasoutDf$DEL_ACSNOM[2:nrow(ldasoutDf)] <- diff(ldasoutDf$ACSNOM) }
if ("ACSNOW" %in% colnames(ldasoutDf)) {
ldasoutDf$DEL_ACSNOW[2:nrow(ldasoutDf)] <- diff(ldasoutDf$ACSNOW) }
if ("SNEQV" %in% colnames(ldasoutDf)) {
ldasoutDf$DEL_SNEQV[2:nrow(ldasoutDf)] <- diff(ldasoutDf$SNEQV) }
if ("SNOWH" %in% colnames(ldasoutDf)) {
ldasoutDf$DEL_SNOWH[2:nrow(ldasoutDf)] <- diff(ldasoutDf$SNOWH) }
} else {
idList <- unique(ldasoutDf[,idCol])
for (i in 1:length(idList)) {
if ("ACCPRCP" %in% colnames(ldasoutDf)) {
tmp <- subset(ldasoutDf$ACCPRCP, ldasoutDf[,idCol]==idList[i])
tmp <- c(NA, diff(tmp))
ldasoutDf$DEL_ACCPRCP[ldasoutDf[,idCol]==idList[i]] <- tmp }
if ("ACCECAN" %in% colnames(ldasoutDf)) {
tmp <- subset(ldasoutDf$ACCECAN, ldasoutDf[,idCol]==idList[i])
tmp <- c(NA, diff(tmp))
ldasoutDf$DEL_ACCECAN[ldasoutDf[,idCol]==idList[i]] <- tmp }
if ("ACCEDIR" %in% colnames(ldasoutDf)) {
tmp <- subset(ldasoutDf$ACCEDIR, ldasoutDf[,idCol]==idList[i])
tmp <- c(NA, diff(tmp))
ldasoutDf$DEL_ACCEDIR[ldasoutDf[,idCol]==idList[i]] <- tmp }
if ("ACCETRAN" %in% colnames(ldasoutDf)) {
tmp <- subset(ldasoutDf$ACCETRAN, ldasoutDf[,idCol]==idList[i])
tmp <- c(NA, diff(tmp))
ldasoutDf$DEL_ACCETRAN[ldasoutDf[,idCol]==idList[i]] <- tmp }
if ("ACCET" %in% colnames(ldasoutDf)) {
tmp <- subset(ldasoutDf$ACCET, ldasoutDf[,idCol]==idList[i])
tmp <- c(NA, diff(tmp))
ldasoutDf$DEL_ACCET[ldasoutDf[,idCol]==idList[i]] <- tmp }
if ("SFCRNOFF" %in% colnames(ldasoutDf)) {
tmp <- subset(ldasoutDf$SFCRNOFF, ldasoutDf[,idCol]==idList[i])
tmp <- c(NA, diff(tmp))
ldasoutDf$DEL_SFCRNOFF[ldasoutDf[,idCol]==idList[i]] <- tmp }
if ("UGDRNOFF" %in% colnames(ldasoutDf)) {
tmp <- subset(ldasoutDf$UGDRNOFF, ldasoutDf[,idCol]==idList[i])
tmp <- c(NA, diff(tmp))
ldasoutDf$DEL_UGDRNOFF[ldasoutDf[,idCol]==idList[i]] <- tmp }
if ("ACSNOM" %in% colnames(ldasoutDf)) {
tmp <- subset(ldasoutDf$ACSNOM, ldasoutDf[,idCol]==idList[i])
tmp <- c(NA, diff(tmp))
ldasoutDf$DEL_ACSNOM[ldasoutDf[,idCol]==idList[i]] <- tmp }
if ("ACSNOW" %in% colnames(ldasoutDf)) {
tmp <- subset(ldasoutDf$ACSNOW, ldasoutDf[,idCol]==idList[i])
tmp <- c(NA, diff(tmp))
ldasoutDf$DEL_ACSNOW[ldasoutDf[,idCol]==idList[i]] <- tmp }
if ("SNEQV" %in% colnames(ldasoutDf)) {
tmp <- subset(ldasoutDf$SNEQV, ldasoutDf[,idCol]==idList[i])
tmp <- c(NA, diff(tmp))
ldasoutDf$DEL_SNEQV[ldasoutDf[,idCol]==idList[i]] <- tmp }
if ("SNOWH" %in% colnames(ldasoutDf)) {
tmp <- subset(ldasoutDf$SNOWH, ldasoutDf[,idCol]==idList[i])
tmp <- c(NA, diff(tmp))
ldasoutDf$DEL_SNOWH[ldasoutDf[,idCol]==idList[i]] <- tmp }
}
ldasoutDf <- ldasoutDf[order(ldasoutDf[,idCol], ldasoutDf[,"POSIXct"]),]
}
ldasoutDf
}
#' Calculate water balance from WRF-Hydro (w/NoahMP) output
#'
#' \code{CalcNoahmpWatBudg} calculates water budget components from WRF-Hydro
#' (w/NoahMP) model output.
#'
#' \code{CalcNoahmpWatBudg} reads dataframes derived from WRF-Hydro output
#' (i.e., using \code{\link{GetMultiNcdf}}) and calculates water budget
#' partitioning (e.g., surface runoff, evaporation, groundwater). Assumes
#' WRF-Hydro output dataframes have already been masked to the desired basin.
#' See \code{\link{GetMultiNcdf}} documentation for examples of how to do this.
#' NOTE: Currently only works for model runs using NoahMP as the LSM.
#'
#' REQUIRED variables (these terms must be in your input dataframes): \itemize{
#' \item LDASOUT: ACCPRCP, ACCECAN, ACCETRAN, ACCEDIR, SFCRNOFF, UGDRNOFF,
#' SOIL_M (all layers), SNEQV, CANICE, CANLIQ \item RTOUT (optional, use if
#' overland or subsurface routing were activated): QSTRMVOLRT, SFCHEADSUBRT,
#' QBDRY \item GWOUT (optional, use if groundwater bucket model was activated):
#' q_cms, POSIXct }
#'
#' OUTPUT water budget terms (may vary depending on model configuration):
#' \itemize{ \item LSM_PRCP: Total precipitation (mm) \item LSM_ECAN: Total
#' canopy evaporation (mm) \item LSM_ETRAN: Total transpiration (mm) \item
#' LSM_EDIR: Total surface evaporation (mm) \item LSM_DELSWE: Change in snowpack
#' snow water equivalent (mm) \item LSM_DELCANWAT: Change in canopy water
#' storage (liquid + ice) (mm) \item LSM_SFCRNOFF: Surface runoff from LSM
#' \emph{(for an LSM-only run)} (mm) \item LSM_UGDRNOFF: Subsurface runoff from
#' LSM \emph{(for an LSM-only run)} (mm) \item LSM_DELSOILM: Change in total
#' soil moisture storage (mm) \item HYD_QSTRMVOL: Total runoff into channel from
#' land \emph{(routing model only)} (mm) \item HYD_DELSFCHEAD: Change in
#' surface storage \emph{(routing model only)} (mm) \item HYD_QBDRY: Total flow
#' outside of domain \emph{(routing model only)} (mm) \item HYD_GWOUT: Total
#' groundwater outflow \emph{(routing model only)} (mm) \item
#' WB_SFCRNOFF: Total surface runoff used in the water budget calculation
#' (\emph{either} LSM_SFCRNOFF or HYD_QSTRMVOL) (mm) \item WB_GWOUT: Total
#' groundwater outflow used in the water budget calculation (\emph{either}
#' LSM_UGDRNOFF or HYD_GWOUT) (mm) \item WB_DELGWSTOR: Change in groundwater
#' storage (adjusted by surface runoff when no surface routing active) (mm)
#' \item ERROR: Remainder in water budget (mm)
#' \item RUN_FRAC: Runoff fraction, runoff/precipitation \item EVAP_FRAC:
#' Evaporative fraction, evapotranspiration/precipitation \item STOR_FRAC:
#' Change in storage fraction, storagechange/precipitation }
#'
#' @param ldasoutDf The LDASOUT dataframe (required)
#' @param rtoutDf The RTOUT dataframe, if overland or subsurface routing was
#' turned on (default=NULL)
#' @param gwoutDf The GW_OUT dataframe, if groundwater model was turned on
#' (default=NULL)
#' @param sfcrt A flag whether surface overland flow routing was active. All
#' other routing options are determined based on the input dataframes, as
#' needed (e.g., if gwoutDf is provided, it is assumed that the groundwater
#' model was active). (default=FALSE)
#' @param soildeps A list of soil layer depths in mm (top to bottom,
#' default=c(100, 300, 600, 1000))
#' @param basarea The basin area in square km (necessary only if gwoutDf is
#' provided)
#' @return A new dataframe containing the water budget components in mm.
#'
#' @examples
#' ## Take a NoahMP LDASOUT dataframe for a model run of Fourmile Creek with no routing
#' ## options turned on and return a water budget summary.
#'
#' \dontrun{
#' wb.nort.fc <- CalcNoahmpWatBudg(modLDASOUT1d.nort.fc)
#' wb.nort.fc
#' ##
#' ## Take NoahMP LDASOUT, HYDRO model RTOUT, and GW_outflow dataframes for a model
#' ## run of Fourmile Creek with subsurface, overland, and groundwater routing options
#' ## turned on and return a water budget summary. The default soil depths were used
#' ## and the basin is 63.1 km2. NOTE: We MUST specify with the sfcrt flag that overland
#' ## flow routing was turned on. Otherwise the LSM surface runoff term is used.
#'
#' wb.allrt.fc <- CalcNoahmpWatBudg(modLDASOUT1d.allrt.fc,
#' modRTOUT1h.allrt.fc,
#' modGWOUT1h.allrt.fc,
#' sfcrt=TRUE, basarea=63.1)
#' wb.allrt.fc
#' }
#' @keywords manip
#' @concept modelEval
#' @family modelEvaluation
#' @export
CalcNoahmpWatBudg <- function(ldasoutDf, rtoutDf=NULL, gwoutDf=NULL, sfcrt=FALSE,
soildeps=c(100,300,600,1000), basarea=NULL) {
wbDf <- as.data.frame(t(as.matrix(rep(0, 6))))
colnames(wbDf) <- c("LSM_PRCP","LSM_ECAN","LSM_ETRAN","LSM_EDIR",
"LSM_DELSWE","LSM_DELCANWAT")
# LSM water fluxes for all model cases
wbDf[1,1] <- ldasoutDf$ACCPRCP[nrow(ldasoutDf)]-ldasoutDf$ACCPRCP[1]
wbDf[1,2] <- ldasoutDf$ACCECAN[nrow(ldasoutDf)]-ldasoutDf$ACCECAN[1]
wbDf[1,3] <- ldasoutDf$ACCETRAN[nrow(ldasoutDf)]-ldasoutDf$ACCETRAN[1]
wbDf[1,4] <- ldasoutDf$ACCEDIR[nrow(ldasoutDf)]-ldasoutDf$ACCEDIR[1]
wbDf[1,5] <- ldasoutDf$SNEQV[nrow(ldasoutDf)]-ldasoutDf$SNEQV[1]
wbDf[1,6] <- (ldasoutDf$CANICE[nrow(ldasoutDf)] + ldasoutDf$CANLIQ[nrow(ldasoutDf)]) -
(ldasoutDf$CANICE[1] + ldasoutDf$CANLIQ[1])
# LSM overlap water fluxes
wbDf[1,"LSM_SFCRNOFF"] <- ldasoutDf$SFCRNOFF[nrow(ldasoutDf)]-ldasoutDf$SFCRNOFF[1]
wbDf[1,"LSM_UGDRNOFF"] <- ldasoutDf$UGDRNOFF[nrow(ldasoutDf)]-ldasoutDf$UGDRNOFF[1]
numsoil <- length(soildeps)
soilm <- 0.0
for (i in 1:numsoil) {
soilm <- soilm + (ldasoutDf[nrow(ldasoutDf),paste0("SOIL_M",i)]-
ldasoutDf[1,paste0("SOIL_M",i)])*soildeps[i]
}
wbDf[1,"LSM_DELSOILM"] <- soilm
# HYDRO surface/subsurface water fluxes
if (!is.null(rtoutDf)) {
wbDf[1,"HYD_QSTRMVOL"] <- rtoutDf$QSTRMVOLRT[nrow(rtoutDf)] -
rtoutDf$QSTRMVOLRT[1]
wbDf[1,"HYD_DELSFCHEAD"] <- rtoutDf$SFCHEADSUBRT[nrow(rtoutDf)] -
rtoutDf$SFCHEADSUBRT[1]
wbDf[1,"HYD_QBDRY"] <- -(rtoutDf$QBDRYRT[nrow(rtoutDf)] - rtoutDf$QBDRYRT[1])
}
else {
message('Message: No routing output dataframe (rtoutDf) was provided. Proceeding using LSM surface runoff.')
wbDf[1,"HYD_QSTRMVOL"] <- NA
wbDf[1,"HYD_DELSFCHEAD"] <- NA
wbDf[1,"HYD_QBDRY"] <- NA
}
# HYDRO groundwater fluxes
if (!is.null(gwoutDf)) {
if (!is.null(basarea)) {
dt <- as.integer(difftime(gwoutDf$POSIXct[2],gwoutDf$POSIXct[1],units="secs"))
wbDf[1,"HYD_GWOUT"] <- sum(gwoutDf$q_cms/(basarea*1000*1000)*1000*dt, na.rm=T)
# wbDf[1,"HYD_DELGWSTOR"] <- (ldasoutDf$UGDRNOFF[nrow(ldasoutDf)] -
# ldasoutDf$UGDRNOFF[1]) -
# wbDf$HYD_GWOUT
}
else { stop('Error: Groundwater outflow dataframe (gwoutDf) provided but no basin area (basarea). Please provide the basin area.') }
}
else {
message('Message: No groundwater outflow dataframe (gwoutDf) was provided. Proceeding using LSM underground runoff.')
wbDf[1,"HYD_GWOUT"] <- NA
# wbDf[1,"HYD_DELGWSTOR"] <- NA
}
# Overland routing check
if ( sfcrt & is.null(rtoutDf)) {
stop('Error: Surface overland routing (sfcrt) is active (TRUE) but no routing output dataframe (rtoutDf) was specified.')
}
if ( !sfcrt & !is.null(rtoutDf) ) {
message('Message: A routing output file (rtoutDf) was provided but surface overland routing (sfcrt) is inactive (FALSE). Proceeding using LSM surface runoff.')
}
wbDf[1,"WB_SFCRNOFF"] <- ifelse(sfcrt, wbDf[1, "HYD_QSTRMVOL"], wbDf[1, "LSM_SFCRNOFF"])
# Groundwater routing check
wbDf[1,"WB_GWOUT"] <- if (!is.null(gwoutDf)) {
if (!sfcrt) { wbDf[1, "HYD_GWOUT"] - wbDf[1, "LSM_SFCRNOFF"] }
else { wbDf[1, "HYD_GWOUT"] }
}
else { wbDf[1, "LSM_UGDRNOFF"] }
wbDf[1,"WB_DELGWSTOR"] <- (ldasoutDf$UGDRNOFF[nrow(ldasoutDf)]-ldasoutDf$UGDRNOFF[1]) -
wbDf$WB_GWOUT
# Water budget error
wbDf[1,"ERROR"] <- with( wbDf, LSM_PRCP - LSM_ECAN - LSM_ETRAN - LSM_EDIR -
WB_SFCRNOFF -
ifelse(is.na(HYD_QBDRY), 0.0, HYD_QBDRY) -
WB_GWOUT - LSM_DELSOILM -
LSM_DELSWE - LSM_DELCANWAT -
ifelse(is.na(HYD_DELSFCHEAD), 0.0, HYD_DELSFCHEAD) -
ifelse(is.na(WB_DELGWSTOR), 0.0, WB_DELGWSTOR) )
# Calculate various fractions
wbDf[1,"RUN_FRAC"] <- with( wbDf, (WB_SFCRNOFF +
ifelse(is.na(HYD_QBDRY), 0.0, HYD_QBDRY) +
WB_GWOUT) / LSM_PRCP )
wbDf[1,"EVAP_FRAC"] <- with( wbDf, (LSM_ECAN + LSM_ETRAN + LSM_EDIR) / LSM_PRCP )
wbDf[1,"STOR_FRAC"] <- with( wbDf, (LSM_DELSOILM + LSM_DELSWE + LSM_DELCANWAT +
ifelse(is.na(HYD_DELSFCHEAD), 0.0, HYD_DELSFCHEAD) +
ifelse(is.na(WB_DELGWSTOR), 0.0, WB_DELGWSTOR) ) /
LSM_PRCP )
wbDf
}
#' Computes model performance statistics for WRF-Hydro flux output
#'
#' \code{CalcModPerf} calculates model performance statistics for flux output.
#'
#' \code{CalcModPerf} reads a model flux time series (i.e., created using
#' \code{\link{ReadFrxstPts}}) and an observation time series (i.e., created
#' using \code{\link{ReadUsgsGage}}) and calculates model performance statistics
#' (Nash-Sutcliffe Efficiency, Rmse, etc.) at various time scales and for low
#' and high fluxes. The tool will subset data to matching time periods (e.g., if
#' the observed data is at 5-min increments and modelled data is at 1-hr
#' increments, the tool will subset the observed data to select only
#' observations on the matching hour break).
#'
#' Performance Statistics: \cr (mod = model output, obs = observations, n =
#' sample size) \itemize{ \item n: sample size \item nse: Nash-Sutcliffe
#' Efficiency \deqn{nse = 1 - ( sum((obs - mod)^2) / sum((obs - mean(obs))^2) )
#' } \item nselog: log-transformed Nash-Sutcliffe Efficiency \deqn{nselog = 1 -
#' ( sum((log(obs) - log(mod))^2) / sum((log(obs) - mean(log(obs)))^2) ) } \item
#' cor: correlation coefficient \deqn{cor = cor(mod, obs) } \item rmse: root
#' mean squared error \deqn{rmse = sqrt( sum((mod - obs)^2) / n ) } \item
#' rmsenorm: normalized root mean squared error \deqn{rmsenorm = rmse /
#' (max(obs) - min(obs)) } \item bias: percent bias \deqn{bias = sum(mod - obs)
#' / sum(obs) * 100 } \item msd: mean signed deviation \deqn{msd = mean(mod -
#' obs) } \item mae: mean absolute error \deqn{mae = mean(abs(mod -
#' obs)) } \item errcom: error in the center-of-mass of the flux, where
#' center-of-mass is the hour/day when 50\% of daily/monthly/water-year flux has
#' occurred. Reported as number of hours for daily time scale and number of days
#' for monthly and yearly time scales. \item errmaxt: Error in the time of
#' maximum flux. Reported as number of hours for daily time scale and number of
#' days for monthly and yearly time scales). \item errfdc: Error in the
#' integrated flow duration curve between 0.05 and 0.95 exceedance thresholds
#' (in native flow units). }
#'
#' Time scales/Flux types: \itemize{ \item ts = native model/observation time
#' step (e.g., hourly) \item daily = daily time step \item monthly = monthly
#' time step \item yearly = water-year time step \item max10 = high flows;
#' restricted to the portion of the time series where the observed flux is in
#' the highest 10\% (native time step) \item min10 = low flows; restricted to
#' the portion of the time series where the observed flux is in the lowest 10\%
#' (native time step) }
#'
#' @param flxDf.mod The flux output dataframe (required). Assumes only one
#' forecast point per file, so if you have multiple forecast points in your
#' output dataframe, use subset to isolate a single forecast point's data.
#' Also assumes model output and observation both contain POSIXct fields
#' (called "POSIXct").
#' @param flxDf.obs The observed flux dataframe. Assumes only one observation
#' point per file, so if you have multiple observation points in your
#' dataframe, use subset to isolate a single point's data. Also assumes model
#' output and observation both contain POSIXct fields (called "POSIXct").
#' @param flxCol.mod The column name for the flux time series for the MODEL data
#' (default="q_cms")
#' @param flxCol.obs The column name for the flux time series for the OBSERVED
#' data (default="q_cms")
#' @param stdate Start date for statistics (DEFAULT=NULL, all records will be
#' used). Date MUST be specified in POSIXct format with appropriate timezone
#' (e.g., as.POSIXct("2013-05-01 00:00:00", format="\%Y-\%m-\%d \%H:\%M:\%S",
#' tz="UTC"))
#' @param enddate End date for statistics (DEFAULT=NULL, all records will be
#' used). Date MUST be specified in POSIXct format with appropriate timezone
#' (e.g., as.POSIXct("2013-05-01 00:00:00", format="\%Y-\%m-\%d \%H:\%M:\%S",
#' tz="UTC"))
#' @param subdivisions Number of subdivisions used in flow duration curve
#' integration (DEFAULT=1000); increase value if integrate function throws an
#' error.
#' @return A new dataframe containing the model performance statistics.
#'
#' @examples
#' ## Take forecast point model output for Fourmile Creek (modStrh.mod1.fc) and a
#' ## corresponding USGS gage observation file (obsStrh.fc), both at an hourly time step,
#' ## and calculate model performance statistics. The model forecast point data was
#' ## imported using ReadFrxstPts and the gage observation data was imported using
#' ## ReadUsgsGage.
#'
#'
#' \dontrun{
#' CalcModPerf(modStrd.chrt.fc, obsStr5min.fc)
#'
#' > Output:
#' n nse nselog cor rmse rmsenorm bias msd mae errcom errmaxt errfdc
#' units count unitless unitless unitless flux units unitless percent flux units flux units hours|days hours|days flux units
#' t 116 0.81 0.67 0.91 0.211 8.48 7.1 0.0313 0.1385 <NA> <NA> 0.04
#' daily 116 0.81 0.67 0.91 0.211 8.48 7.1 0.0313 0.1385 0 0 0.04
#' monthly 5 0.92 0.82 0.96 0.1073 10.35 9.1 0.0324 0.0833 -1.2 0.2 <NA>
#' yearly 1 <NA> <NA> <NA> <NA> <NA> 7.1 0.0313 0.0313 0 4 <NA>
#' max10 12 0.46 0.58 0.71 0.348 23.63 -2.8 -0.044 0.248 <NA> <NA> <NA>
#' min10 12 -385.99 -14.59 0.35 0.1678 538.78 112.4 0.056 0.0592 <NA> <NA> <NA>
#' }
#' @keywords univar ts
#' @concept modelEval
#' @family modelEvaluation
#' @export
CalcModPerf <- function (flxDf.mod, flxDf.obs, flxCol.mod="q_cms", flxCol.obs="q_cms",
stdate=NULL, enddate=NULL, subdivisions=1000) {
# Internal functions
which.max.dt <- function(dd, qcol, dtcol) { dd[which.max(dd[,qcol]), dtcol] }
CalcCOM.dt <- function(dd, qcol, dtcol) { dd[CalcCOM(dd[,qcol]), dtcol] }
# Prepare data
if (!is.null(stdate) & is.null(enddate)) {
flxDf.obs <- subset(flxDf.obs, POSIXct>=stdate)
flxDf.mod <- subset(flxDf.mod, POSIXct>=stdate)
}
if (is.null(stdate) & !is.null(enddate)) {
flxDf.obs <- subset(flxDf.obs, POSIXct<=enddate)
flxDf.mod <- subset(flxDf.mod, POSIXct<=enddate)
}
if (!is.null(stdate) & !is.null(enddate)) {
flxDf.obs <- subset(flxDf.obs, POSIXct>=stdate & POSIXct<=enddate)
flxDf.mod <- subset(flxDf.mod, POSIXct>=stdate & POSIXct<=enddate)
}
flxDf.mod$qcomp <- flxDf.mod[,flxCol.mod]
flxDf.obs$qcomp <- flxDf.obs[,flxCol.obs]
# model timestep in secs
modT <- as.integer(flxDf.mod$POSIXct[2])-as.integer(flxDf.mod$POSIXct[1])
#if (as.integer(flxDf.obs$POSIXct[2])-as.integer(flxDf.obs$POSIXct[1]) >= 86400) {
# flxDf.obs$POSIXct=as.POSIXct(round(flxDf.obs$POSIXct,"days"), tz="UTC")}
flxDf.mod <- merge(flxDf.mod[c("POSIXct","qcomp")], flxDf.obs[c("POSIXct","qcomp")],
by<-c("POSIXct"), suffixes=c(".mod",".obs"))
flxDf.mod <- subset(flxDf.mod, !is.na(flxDf.mod$qcomp.mod) & !is.na(flxDf.mod$qcomp.obs))
flxDf.mod <- CalcDates(flxDf.mod)
flxDf.mod$date <- as.POSIXct(trunc(flxDf.mod$POSIXct, "days"))
results <- as.data.frame(matrix(nrow = 7, ncol = 12))
colnames(results) = c("n", "nse", "nselog", "cor", "rmse", "rmsenorm", "bias", "msd",
"mae", "errcom", "errmaxt", "errfdc")
rownames(results) = c("units", "t", "daily", "monthly", "yearly", "max10", "min10")
exclvars <- names(flxDf.mod) %in% c("POSIXct", "secs", "timest", "date", "stat")
# Base aggregations
flxDf.mod.d <- aggregate(flxDf.mod[!exclvars], by = list(flxDf.mod$date),
CalcMeanNarm)
flxDf.mod.d$mwy <- paste0(flxDf.mod.d$month, "-", flxDf.mod.d$wy)
flxDf.mod.mwy <- aggregate(flxDf.mod[c("qcomp.mod","qcomp.obs")],
by = list(flxDf.mod$month, flxDf.mod$wy), CalcMeanNarm)
flxDf.mod.wy <- aggregate(flxDf.mod[c("qcomp.mod","qcomp.obs")],
by = list(flxDf.mod$wy), CalcMeanNarm)
# Time of center of mass aggregations
tmp.mod <- plyr::ddply(flxDf.mod, plyr::.(date), CalcCOM.dt, qcol="qcomp.mod",
dtcol="POSIXct")
tmp.obs <- plyr::ddply(flxDf.mod, plyr::.(date), CalcCOM.dt, qcol="qcomp.obs",
dtcol="POSIXct")
flxDf.mod.dcom <- merge(tmp.mod, tmp.obs, by=c("date"))
colnames(flxDf.mod.dcom) <- c("date", "qcomp.mod", "qcomp.obs")
tmp.mod <- plyr::ddply(flxDf.mod.d, plyr::.(month, wy),
CalcCOM.dt, qcol="qcomp.mod", dtcol="Group.1")
tmp.obs <- plyr::ddply(flxDf.mod.d, plyr::.(month, wy),
CalcCOM.dt, qcol="qcomp.obs", dtcol="Group.1")
flxDf.mod.mwycom <- merge(tmp.mod, tmp.obs, by=c("month", "wy"))
colnames(flxDf.mod.mwycom) <- c("month", "wy", "qcomp.mod", "qcomp.obs")
tmp.mod <- plyr::ddply(flxDf.mod.d, plyr::.(wy),
CalcCOM.dt, qcol="qcomp.mod", dtcol="Group.1")
tmp.obs <- plyr::ddply(flxDf.mod.d, plyr::.(wy),
CalcCOM.dt, qcol="qcomp.obs", dtcol="Group.1")
flxDf.mod.wycom <- merge(tmp.mod, tmp.obs, by=c("wy"))
colnames(flxDf.mod.wycom) <- c("wy", "qcomp.mod", "qcomp.obs")
# Time of max aggregations
tmp.mod <- plyr::ddply(flxDf.mod, plyr::.(date),
which.max.dt, qcol="qcomp.mod", dtcol="POSIXct")
tmp.obs <- plyr::ddply(flxDf.mod, plyr::.(date),
which.max.dt, qcol="qcomp.obs", dtcol="POSIXct")
flxDf.mod.dmax <- merge(tmp.mod, tmp.obs, by=c("date"))
colnames(flxDf.mod.dmax) <- c("date", "qcomp.mod", "qcomp.obs")
tmp.mod <- plyr::ddply(flxDf.mod.d, plyr::.(month, wy),
which.max.dt, qcol="qcomp.mod", dtcol="Group.1")
tmp.obs <- plyr::ddply(flxDf.mod.d, plyr::.(month, wy),
which.max.dt, qcol="qcomp.obs", dtcol="Group.1")
flxDf.mod.mwymax <- merge(tmp.mod, tmp.obs, by=c("month", "wy"))
colnames(flxDf.mod.mwymax) <- c("month", "wy", "qcomp.mod", "qcomp.obs")
tmp.mod <- plyr::ddply(flxDf.mod.d, plyr::.(wy),
which.max.dt, qcol="qcomp.mod", dtcol="Group.1")
tmp.obs <- plyr::ddply(flxDf.mod.d, plyr::.(wy),
which.max.dt, qcol="qcomp.obs", dtcol="Group.1")
flxDf.mod.wymax <- merge(tmp.mod, tmp.obs, by=c("wy"))
colnames(flxDf.mod.wymax) <- c("wy", "qcomp.mod", "qcomp.obs")
# NAs cause which.max to return a list, so force back to ints
# (out for now since we pre-screen for NAs)
#flxDf.mod.dmax$qcomp.mod <- as.integer(flxDf.mod.dmax$qcomp.mod)
#flxDf.mod.dmax$qcomp.obs <- as.integer(flxDf.mod.dmax$qcomp.obs)
#flxDf.mod.mwymax$qcomp.mod <- as.integer(flxDf.mod.mwymax$qcomp.mod)
#flxDf.mod.mwymax$qcomp.obs <- as.integer(flxDf.mod.mwymax$qcomp.obs)
#flxDf.mod.wymax$qcomp.mod <- as.integer(flxDf.mod.wymax$qcomp.mod)
#flxDf.mod.wymax$qcomp.obs <- as.integer(flxDf.mod.wymax$qcomp.obs)
# Mins and Maxes
flxDf.mod.max10 <- subset(flxDf.mod, flxDf.mod$qcomp.obs >=
quantile(flxDf.mod$qcomp.obs, 0.90, na.rm=TRUE))
flxDf.mod.min10 <- subset(flxDf.mod, flxDf.mod$qcomp.obs <=
quantile(flxDf.mod$qcomp.obs, 0.10, na.rm=TRUE))
# FDCs
flxDf.mod <- CalcFdc(flxDf.mod, "qcomp.mod")
flxDf.mod <- CalcFdc(flxDf.mod, "qcomp.obs")
flxDf.mod.d <- CalcFdc(flxDf.mod.d, "qcomp.mod")
flxDf.mod.d <- CalcFdc(flxDf.mod.d, "qcomp.obs")
# Compile summary statistics
results["t", "n"] <- length(flxDf.mod$qcomp.mod)
results["t", "nse"] <- round(Nse(flxDf.mod$qcomp.mod, flxDf.mod$qcomp.obs), 2)
results["t", "nselog"] <- round(NseLog(flxDf.mod$qcomp.mod, flxDf.mod$qcomp.obs), 2)
results["t", "cor"] <- round(cor(flxDf.mod$qcomp.mod, flxDf.mod$qcomp.obs,
use="na.or.complete"), 2)
results["t", "rmse"] <- round(Rmse(flxDf.mod$qcomp.mod, flxDf.mod$qcomp.obs), 4)
results["t", "rmsenorm"] <- round(RmseNorm(flxDf.mod$qcomp.mod, flxDf.mod$qcomp.obs), 2)
results["t", "bias"] <- round(sum(flxDf.mod$qcomp.mod-flxDf.mod$qcomp.obs, na.rm=T)/
sum(flxDf.mod$qcomp.obs, na.rm=TRUE) * 100, 1)
results["t", "msd"] <- round(mean(flxDf.mod$qcomp.mod-flxDf.mod$qcomp.obs, na.rm=T), 4)
results["t", "mae"] <- round(mean(abs(flxDf.mod$qcomp.mod-flxDf.mod$qcomp.obs),
na.rm=T), 4)
results["t", "errcom"] <- NA
results["t", "errmaxt"] <- NA
results["t", "errfdc"] <- round(integrate(splinefun(flxDf.mod[,"qcomp.mod.fdc"],
flxDf.mod[,"qcomp.mod"],
method='natural'),
0.05, 0.95, subdivisions=subdivisions)$value -
integrate(splinefun(flxDf.mod[,"qcomp.obs.fdc"],
flxDf.mod[,"qcomp.obs"], method='natural'),
0.05, 0.95, subdivisions=subdivisions)$value, 2 )
results["daily", "n"] <- length(flxDf.mod.d$qcomp.mod)
results["daily", "nse"] <- round(Nse(flxDf.mod.d$qcomp.mod, flxDf.mod.d$qcomp.obs), 2)
results["daily", "nselog"] <- round(NseLog(flxDf.mod.d$qcomp.mod, flxDf.mod.d$qcomp.obs), 2)
results["daily", "cor"] <- round(cor(flxDf.mod.d$qcomp.mod, flxDf.mod.d$qcomp.obs,
use="na.or.complete"), 2)
results["daily", "rmse"] <- round(Rmse(flxDf.mod.d$qcomp.mod, flxDf.mod.d$qcomp.obs), 4)
results["daily", "rmsenorm"] <- round(RmseNorm(flxDf.mod.d$qcomp.mod,
flxDf.mod.d$qcomp.obs), 2)
results["daily", "bias"] <- round(sum(flxDf.mod.d$qcomp.mod-flxDf.mod.d$qcomp.obs, na.rm=T)/
sum(flxDf.mod.d$qcomp.obs, na.rm=TRUE) * 100, 1)
results["daily", "msd"] <- round(mean(flxDf.mod.d$qcomp.mod-flxDf.mod.d$qcomp.obs, na.rm=T), 4)
results["daily", "mae"] <- round(mean(abs(flxDf.mod.d$qcomp.mod-flxDf.mod.d$qcomp.obs),
na.rm=T), 4)
results["daily", "errcom"] <- round(mean(as.numeric(difftime(flxDf.mod.dcom$qcomp.mod,
flxDf.mod.dcom$qcomp.obs,
units="hours")), na.rm=T), 2)
results["daily", "errmaxt"] <- round(mean(as.numeric(difftime(flxDf.mod.dmax$qcomp.mod,
flxDf.mod.dmax$qcomp.obs,
units="hours")), na.rm=T), 2)
results["daily", "errfdc"] <- round(integrate(splinefun(flxDf.mod.d[,"qcomp.mod.fdc"],
flxDf.mod.d[,"qcomp.mod"],
method='natural'), 0.05, 0.95,
subdivisions=subdivisions)$value -
integrate(splinefun(flxDf.mod.d[,"qcomp.obs.fdc"],
flxDf.mod.d[,"qcomp.obs"],
method='natural'), 0.05, 0.95,
subdivisions=subdivisions)$value, 2 )
results["monthly", "n"] <- length(flxDf.mod.mwy$qcomp.mod)
results["monthly", "nse"] <- round(Nse(flxDf.mod.mwy$qcomp.mod, flxDf.mod.mwy$qcomp.obs), 2)
results["monthly", "nselog"] <- round(NseLog(flxDf.mod.mwy$qcomp.mod,
flxDf.mod.mwy$qcomp.obs), 2)
results["monthly", "cor"] <- round(cor(flxDf.mod.mwy$qcomp.mod, flxDf.mod.mwy$qcomp.obs,
use="na.or.complete"), 2)
results["monthly", "rmse"] <- round(Rmse(flxDf.mod.mwy$qcomp.mod, flxDf.mod.mwy$qcomp.obs), 4)
results["monthly", "rmsenorm"] <- round(RmseNorm(flxDf.mod.mwy$qcomp.mod,
flxDf.mod.mwy$qcomp.obs), 2)
results["monthly", "bias"] <- round(sum(flxDf.mod.mwy$qcomp.mod-flxDf.mod.mwy$qcomp.obs,
na.rm=T)/
sum(flxDf.mod.mwy$qcomp.obs, na.rm=TRUE) * 100, 1)
results["monthly", "msd"] <- round(mean(flxDf.mod.mwy$qcomp.mod-flxDf.mod.mwy$qcomp.obs,
na.rm=T), 4)
results["monthly", "mae"] <- round(mean(abs(flxDf.mod.mwy$qcomp.mod-flxDf.mod.mwy$qcomp.obs),
na.rm=T), 4)
results["monthly", "errcom"] <- round(mean(as.numeric(difftime(flxDf.mod.mwycom$qcomp.mod,
flxDf.mod.mwycom$qcomp.obs,
units="days")), na.rm=T), 2)
results["monthly", "errmaxt"] <- round(mean(as.numeric(difftime(flxDf.mod.mwymax$qcomp.mod,
flxDf.mod.mwymax$qcomp.obs,
units="days")), na.rm=T), 2)
results["monthly", "errfdc"] <- NA
results["yearly", "n"] <- length(flxDf.mod.wy$qcomp.mod)
results["yearly", "nse"] <- NA
#round(Nse(flxDf.mod.wy$qcomp.mod, flxDf.mod.wy$qcomp.obs), 2)
results["yearly", "nselog"] <- NA
#round(NseLog(flxDf.mod.wy$qcomp.mod, flxDf.mod.wy$qcomp.obs), 2)
results["yearly", "cor"] <- NA
#round(cor(flxDf.mod.wy$qcomp.mod, flxDf.mod.wy$qcomp.obs, use="na.or.complete"), 2)
results["yearly", "rmse"] <- NA
#round(Rmse(flxDf.mod.wy$qcomp.mod, flxDf.mod.wy$qcomp.obs), 2)
results["yearly", "rmsenorm"] <- NA
#round(RmseNorm(flxDf.mod.wy$qcomp.mod, flxDf.mod.wy$qcomp.obs), 2)
results["yearly", "bias"] <- round(sum(flxDf.mod.wy$qcomp.mod-flxDf.mod.wy$qcomp.obs, na.rm=T)/
sum(flxDf.mod.wy$qcomp.obs, na.rm=TRUE) * 100, 1)
results["yearly", "msd"] <- round(mean(flxDf.mod.wy$qcomp.mod-flxDf.mod.wy$qcomp.obs,
na.rm=T), 4)
results["yearly", "mae"] <- round(mean(abs(flxDf.mod.wy$qcomp.mod-flxDf.mod.wy$qcomp.obs),
na.rm=T), 4)
results["yearly", "errcom"] <- round(mean(as.numeric(difftime(flxDf.mod.wycom$qcomp.mod,
flxDf.mod.wycom$qcomp.obs,
units="days")), na.rm=T), 2)
results["yearly", "errmaxt"] <- round(mean(as.numeric(difftime(flxDf.mod.wymax$qcomp.mod,
flxDf.mod.wymax$qcomp.obs,
units="days")), na.rm=T), 2)
results["yearly", "errfdc"] <- NA
results["max10", "n"] <- length(flxDf.mod.max10$qcomp.mod)
results["max10", "nse"] <- round(Nse(flxDf.mod.max10$qcomp.mod, flxDf.mod.max10$qcomp.obs), 2)
results["max10", "nselog"] <- round(NseLog(flxDf.mod.max10$qcomp.mod,
flxDf.mod.max10$qcomp.obs), 2)
results["max10", "cor"] <- round(cor(flxDf.mod.max10$qcomp.mod, flxDf.mod.max10$qcomp.obs,
use="na.or.complete"), 2)
results["max10", "rmse"] <- round(Rmse(flxDf.mod.max10$qcomp.mod, flxDf.mod.max10$qcomp.obs), 4)
results["max10", "rmsenorm"] <- round(RmseNorm(flxDf.mod.max10$qcomp.mod,
flxDf.mod.max10$qcomp.obs), 2)
results["max10", "bias"] <- round(sum(flxDf.mod.max10$qcomp.mod-flxDf.mod.max10$qcomp.obs,
na.rm=T)/
sum(flxDf.mod.max10$qcomp.obs, na.rm=TRUE) * 100, 1)
results["max10", "msd"] <- round(mean(flxDf.mod.max10$qcomp.mod-flxDf.mod.max10$qcomp.obs,
na.rm=T), 4)
results["max10", "mae"] <- round(mean(abs(flxDf.mod.max10$qcomp.mod-flxDf.mod.max10$qcomp.obs),
na.rm=T), 4)
results["max10", "errcom"] <- NA
results["max10", "errmaxt"] <- NA
results["max10", "errfdc"] <- NA
results["min10", "n"] <- length(flxDf.mod.min10$qcomp.mod)
results["min10", "nse"] <- round(Nse(flxDf.mod.min10$qcomp.mod, flxDf.mod.min10$qcomp.obs), 2)
results["min10", "nselog"] <- round(NseLog(flxDf.mod.min10$qcomp.mod,
flxDf.mod.min10$qcomp.obs), 2)
results["min10", "cor"] <- round(cor(flxDf.mod.min10$qcomp.mod, flxDf.mod.min10$qcomp.obs,
use="na.or.complete"), 2)
results["min10", "rmse"] <- round(Rmse(flxDf.mod.min10$qcomp.mod, flxDf.mod.min10$qcomp.obs), 4)
results["min10", "rmsenorm"] <- round(RmseNorm(flxDf.mod.min10$qcomp.mod,
flxDf.mod.min10$qcomp.obs), 2)
results["min10", "bias"] <- round(sum(flxDf.mod.min10$qcomp.mod-flxDf.mod.min10$qcomp.obs,
na.rm=T)/
sum(flxDf.mod.min10$qcomp.obs, na.rm=TRUE) * 100, 1)
results["min10", "msd"] <- round(mean(flxDf.mod.min10$qcomp.mod-flxDf.mod.min10$qcomp.obs,
na.rm=T), 4)
results["min10", "mae"] <- round(mean(abs(flxDf.mod.min10$qcomp.mod-flxDf.mod.min10$qcomp.obs),
na.rm=T), 4)
results["min10", "errcom"] <- NA
results["min10", "errmaxt"] <- NA
results["min10", "errfdc"] <- NA
# Units
results["units",] <- c("count", "unitless", "unitless", "unitless",
"flux units", "unitless", "percent", "flux units",
"flux units", "hours|days", "hours|days",
"flux units")
results
}
#' Computes model performance statistics for WRF-Hydro flux output
#'
#' \code{CalcModPerfMulti} calculates model performance statistics for flux
#' output.
#'
#' \code{CalcModPerfMulti} reads a model flux time series (i.e., created using
#' \code{\link{ReadFrxstPts}}) and an observation time series (i.e., created
#' using \code{\link{ReadUsgsGage}}) and calculates model performance statistics
#' (Nash-Sutcliffe Efficiency, Rmse, etc.) at various time scales and for low
#' and high fluxes. The tool will subset data to matching time periods (e.g., if
#' the observed data is at 5-min increments and modelled data is at 1-hr
#' increments, the tool will subset the observed data to select only
#' observations on the matching hour break).
#'
#' \code{CalcModPerfMulti} calculates the same statistics as
#' \code{\link{CalcModPerf}}) but returns results as a single-row vs. multi-row
#' dataframe. This is intended to streamline compiling statistics for multiple
#' data runs or multiple sites.
#'
#' Performance Statistics: \cr (mod = model output, obs = observations, n =
#' sample size) \itemize{ \item n: sample size \item nse: Nash-Sutcliffe
#' Efficiency \deqn{nse = 1 - ( sum((obs - mod)^2) / sum((obs - mean(obs))^2) )
#' } \item nselog: log-transformed Nash-Sutcliffe Efficiency \deqn{nselog = 1 -
#' ( sum((log(obs) - log(mod))^2) / sum((log(obs) - mean(log(obs)))^2) ) } \item
#' cor: correlation coefficient \deqn{cor = cor(mod, obs) } \item rmse: root
#' mean squared error \deqn{rmse = sqrt( sum((mod - obs)^2) / n ) } \item
#' rmsenorm: normalized root mean squared error \deqn{rmsenorm = rmse /
#' (max(obs) - min(obs)) } \item bias: percent bias \deqn{bias = sum(mod - obs)
#' / sum(obs) * 100 } \item msd: mean signed deviation \deqn{msd = mean(mod -
#' obs) } \item mae: mean absolute error \deqn{mae = mean(abs(mod -
#' obs)) } \item errcom: error in the center-of-mass of the flux, where
#' center-of-mass is the hour/day when 50\% of daily/monthly/water-year flux has
#' occurred. Reported as number of hours for daily time scale and number of days
#' for monthly and yearly time scales. \item errmaxt: Error in the time of
#' maximum flux. Reported as number of hours for daily time scale and number of
#' days for monthly and yearly time scales). \item errfdc: Error in the
#' integrated flow duration curve between 0.05 and 0.95 exceedance thresholds
#' (in native flow units). \item errflash: Error in the Richards-Baker
#' Flashiness Index \deqn{R-B Index = sum(abs(q_t - q_t-1)) / sum(q)}. }
#'
#' Time scales/Flux types: \itemize{ \item t = native model/observation time
#' step (e.g., hourly) \item dy = daily time step \item mo = monthly time step
#' \item yr = water-year time step \item max10 = high flows; restricted to the
#' portion of the time series where the observed flux is in the highest 10\%
#' (native time step) \item min10 = low flows; restricted to the portion of the
#' time series where the observed flux is in the lowest 10\% (native time step)
#' }
#'
#' @param flxDf.mod The flux output dataframe (required). Assumes only one
#' forecast point per file, so if you have multiple forecast points in your
#' output dataframe, use subset to isolate a single forecast point's data.
#' Also assumes model output and observation both contain POSIXct fields
#' (called "POSIXct").
#' @param flxDf.obs The observed flux dataframe. Assumes only one observation
#' point per file, so if you have multiple observation points in your
#' dataframe, use subset to isolate a single point's data. Also assumes model
#' output and observation both contain POSIXct fields (called "POSIXct").
#' @param flxCol.mod The column name for the flux time series for the MODEL data
#' (default="q_cms")
#' @param flxCol.obs The column name for the flux time series for the OBSERVED
#' data (default="q_cms")
#' @param stdate Start date for statistics (DEFAULT=NULL, all records will be
#' used). Date MUST be specified in POSIXct format with appropriate timezone
#' (e.g., as.POSIXct("2013-05-01 00:00:00", format="\%Y-\%m-\%d \%H:\%M:\%S",
#' tz="UTC"))
#' @param enddate End date for statistics (DEFAULT=NULL, all records will be
#' used). Date MUST be specified in POSIXct format with appropriate timezone
#' (e.g., as.POSIXct("2013-05-01 00:00:00", format="\%Y-\%m-\%d \%H:\%M:\%S",
#' tz="UTC"))
#' @param writeOutPairDf character path/name for file to output the data frame constructed
#' before calculating the statistics and various aggregations. NULL by default gives
#' no output.
#' @param fileTag A to be inserted in the file name, just before the file extension.
#' @return A new dataframe containing the model performance statistics.
#'
#' @examples
#' ## Take forecast point model output for Fourmile Creek (modStrh.mod1.fc) and a
#' ## corresponding USGS gage observation file (obsStrh.fc), both at an hourly time
#' ## step, and calculate model performance statistics. The model forecast point
#' ## data was imported using ReadFrxstPts and the gage observation data was
#' ## imported using ReadUsgsGage.
#'
#'
#' \dontrun{
#' CalcModPerfMulti(modStrd.chrt.fc, obsStr5min.fc)
#'
#' > Output:
#' t_n t_nse t_nselog t_cor t_rmse t_rmsenorm t_bias t_msd t_mae t_errfdc t_errflash
#' 116 0.81 0.67 0.91 0.211 8.48 7.1 0.0313 0.1385 0.04 0.0132
#' dy_n dy_nse dy_nselog dy_cor dy_rmse dy_rmsenorm dy_bias dy_msd dy_mae dy_errcom dy_errmaxt dy_errfdc dy_errflash
#' 116 0.81 0.67 0.91 0.211 8.48 7.1 0.0313 0.1385 0 0 0.04 0.0092
#' mo_n mo_nse mo_nselog mo_cor mo_rmse mo_rmsenorm mo_bias mo_msd mo_mae mo_errcom mo_errmaxt yr_n
#' 5 0.92 0.82 0.96 0.1073 10.35 9.1 0.0324 0.0833 -1.2 0.2 1
#' yr_bias yr_msd yr_mae yr_errcom yr_errmaxt max10_n max10_nse max10_nselog max10_cor max10_rmse max10_rmsenorm max10_bias max10_msd max10_mae
#' 7.1 0.0313 0.0313 0 4 12 0.46 0.58 0.71 0.348 23.63 -2.8 -0.044 0.248
#' min10_n min10_nse min10_nselog min10_cor min10_rmse min10_rmsenorm min10_bias min10_msd min10_mae
#' 12 -385.99 -14.59 0.35 0.1678 538.78 112.4 0.056 0.0592
#' }
#' @keywords univar ts
#' @concept modelEval
#' @family modelEvaluation
#' @export
CalcModPerfMulti <- function (flxDf.mod, flxDf.obs, flxCol.mod="q_cms", flxCol.obs="q_cms",
stdate=NULL, enddate=NULL, writeOutPairDf=NULL, fileTag=NULL) {
# Internal functions
which.max.dt <- function(dd, qcol, dtcol) { dd[which.max(dd[,qcol]), dtcol] }
CalcCOM.dt <- function(dd, qcol, dtcol) { dd[CalcCOM(dd[,qcol]), dtcol] }
# Prepare data
if (!is.null(stdate) & is.null(enddate)) {
flxDf.obs <- subset(flxDf.obs, POSIXct>=stdate)
flxDf.mod <- subset(flxDf.mod, POSIXct>=stdate)
}
if (is.null(stdate) & !is.null(enddate)) {
flxDf.obs <- subset(flxDf.obs, POSIXct<=enddate)
flxDf.mod <- subset(flxDf.mod, POSIXct<=enddate)
}
if (!is.null(stdate) & !is.null(enddate)) {
flxDf.obs <- subset(flxDf.obs, POSIXct>=stdate & POSIXct<=enddate)
flxDf.mod <- subset(flxDf.mod, POSIXct>=stdate & POSIXct<=enddate)
}
flxDf.mod$qcomp <- flxDf.mod[,flxCol.mod]
flxDf.obs$qcomp <- flxDf.obs[,flxCol.obs]
# model timestep in secs
modT <- as.integer(flxDf.mod$POSIXct[2])-as.integer(flxDf.mod$POSIXct[1])
#if (as.integer(flxDf.obs$POSIXct[2])-as.integer(flxDf.obs$POSIXct[1]) >= 86400) {
# flxDf.obs$POSIXct=as.POSIXct(round(flxDf.obs$POSIXct,"days"), tz="UTC")}
flxDf.mod <- merge(flxDf.mod[c("POSIXct","qcomp")], flxDf.obs[c("POSIXct","qcomp")],
by<-c("POSIXct"), suffixes=c(".mod",".obs"))
flxDf.mod <- subset(flxDf.mod, !is.na(flxDf.mod$qcomp.mod) & !is.na(flxDf.mod$qcomp.obs))
if (nrow(flxDf.mod) == 0) {
warning("No data matches.")
return(NA)
}
flxDf.mod <- CalcDates(flxDf.mod)
flxDf.mod$date <- as.POSIXct(trunc(flxDf.mod$POSIXct, "days"))
if(!is.null(writeOutPairDf) && is.character(writeOutPairDf) && writeOutPairDf!='') {
modPerfPairDf <- flxDf.mod
wopdFileName <- if(!is.null(fileTag) && fileTag!='') {
wopdBase <- tools::file_path_sans_ext(writeOutPairDf)
wopdExt <- tools::file_ext(writeOutPairDf)
paste0(wopdBase,'.',fileTag,'.',wopdExt)
} else writeOutPairDf
print(wopdFileName)
save(modPerfPairDf, file=wopdFileName)
rm('modPerfPairDf')
}
results <- as.data.frame(matrix(nrow = 1, ncol = 59))
colnames(results) = c("t_n", "t_nse", "t_nselog", "t_cor", "t_rmse", "t_rmsenorm",
"t_bias", "t_msd", "t_mae", "t_errfdc", "t_errflash",
"dy_n", "dy_nse", "dy_nselog", "dy_cor", "dy_rmse",
"dy_rmsenorm", "dy_bias", "dy_msd", "dy_mae",
"dy_errcom", "dy_errmaxt", "dy_errfdc", "dy_errflash",
"mo_n", "mo_nse", "mo_nselog", "mo_cor", "mo_rmse",
"mo_rmsenorm", "mo_bias", "mo_msd", "mo_mae",
"mo_errcom", "mo_errmaxt",
"yr_n", "yr_bias", "yr_msd", "yr_mae", "yr_errcom", "yr_errmaxt",
"max10_n", "max10_nse", "max10_nselog", "max10_cor",
"max10_rmse", "max10_rmsenorm", "max10_bias",
"max10_msd", "max10_mae",
"min10_n", "min10_nse", "min10_nselog", "min10_cor",
"min10_rmse", "min10_rmsenorm", "min10_bias",
"min10_msd", "min10_mae")
exclvars <- names(flxDf.mod) %in% c("POSIXct", "secs", "timest", "date", "stat")
# Base aggregations
flxDf.mod.d <- aggregate(flxDf.mod[!exclvars], by = list(flxDf.mod$date), CalcMeanNarm)
flxDf.mod.d$mwy <- paste0(flxDf.mod.d$month, "-", flxDf.mod.d$wy)
flxDf.mod.mwy <- aggregate(flxDf.mod[c("qcomp.mod","qcomp.obs")],
by = list(flxDf.mod$month, flxDf.mod$wy), CalcMeanNarm)
flxDf.mod.wy <- aggregate(flxDf.mod[c("qcomp.mod","qcomp.obs")],
by = list(flxDf.mod$wy), CalcMeanNarm)
# Time of center of mass aggregations
tmp.mod <- plyr::ddply(flxDf.mod, plyr::.(date), CalcCOM.dt, qcol="qcomp.mod",
dtcol="POSIXct")
tmp.obs <- plyr::ddply(flxDf.mod, plyr::.(date), CalcCOM.dt, qcol="qcomp.obs",
dtcol="POSIXct")
flxDf.mod.dcom <- merge(tmp.mod, tmp.obs, by=c("date"))
colnames(flxDf.mod.dcom) <- c("date", "qcomp.mod", "qcomp.obs")
tmp.mod <- plyr::ddply(flxDf.mod.d, plyr::.(month, wy), CalcCOM.dt, qcol="qcomp.mod",
dtcol="Group.1")
tmp.obs <- plyr::ddply(flxDf.mod.d, plyr::.(month, wy), CalcCOM.dt, qcol="qcomp.obs",
dtcol="Group.1")
flxDf.mod.mwycom <- merge(tmp.mod, tmp.obs, by=c("month", "wy"))
colnames(flxDf.mod.mwycom) <- c("month", "wy", "qcomp.mod", "qcomp.obs")
tmp.mod <- plyr::ddply(flxDf.mod.d, plyr::.(wy), CalcCOM.dt, qcol="qcomp.mod",
dtcol="Group.1")
tmp.obs <- plyr::ddply(flxDf.mod.d, plyr::.(wy), CalcCOM.dt, qcol="qcomp.obs",
dtcol="Group.1")
flxDf.mod.wycom <- merge(tmp.mod, tmp.obs, by=c("wy"))
colnames(flxDf.mod.wycom) <- c("wy", "qcomp.mod", "qcomp.obs")
# Time of max aggregations
tmp.mod <- plyr::ddply(flxDf.mod, plyr::.(date), which.max.dt, qcol="qcomp.mod",
dtcol="POSIXct")
tmp.obs <- plyr::ddply(flxDf.mod, plyr::.(date), which.max.dt, qcol="qcomp.obs",
dtcol="POSIXct")
flxDf.mod.dmax <- merge(tmp.mod, tmp.obs, by=c("date"))
colnames(flxDf.mod.dmax) <- c("date", "qcomp.mod", "qcomp.obs")
tmp.mod <- plyr::ddply(flxDf.mod.d, plyr::.(month, wy), which.max.dt, qcol="qcomp.mod",
dtcol="Group.1")
tmp.obs <- plyr::ddply(flxDf.mod.d, plyr::.(month, wy), which.max.dt, qcol="qcomp.obs",
dtcol="Group.1")
flxDf.mod.mwymax <- merge(tmp.mod, tmp.obs, by=c("month", "wy"))
colnames(flxDf.mod.mwymax) <- c("month", "wy", "qcomp.mod", "qcomp.obs")
tmp.mod <- plyr::ddply(flxDf.mod.d, plyr::.(wy), which.max.dt, qcol="qcomp.mod",
dtcol="Group.1")
tmp.obs <- plyr::ddply(flxDf.mod.d, plyr::.(wy), which.max.dt, qcol="qcomp.obs",
dtcol="Group.1")
flxDf.mod.wymax <- merge(tmp.mod, tmp.obs, by=c("wy"))
colnames(flxDf.mod.wymax) <- c("wy", "qcomp.mod", "qcomp.obs")
# NAs cause which.max to return a list, so force back to ints
# (out for now since we pre-screen for NAs)
#flxDf.mod.dmax$qcomp.mod <- as.integer(flxDf.mod.dmax$qcomp.mod)
#flxDf.mod.dmax$qcomp.obs <- as.integer(flxDf.mod.dmax$qcomp.obs)
#flxDf.mod.mwymax$qcomp.mod <- as.integer(flxDf.mod.mwymax$qcomp.mod)
#flxDf.mod.mwymax$qcomp.obs <- as.integer(flxDf.mod.mwymax$qcomp.obs)
#flxDf.mod.wymax$qcomp.mod <- as.integer(flxDf.mod.wymax$qcomp.mod)
#flxDf.mod.wymax$qcomp.obs <- as.integer(flxDf.mod.wymax$qcomp.obs)
# Mins and Maxes
flxDf.mod.max10 <- subset(flxDf.mod, flxDf.mod$qcomp.obs>=quantile(flxDf.mod$qcomp.obs,
0.90, na.rm=TRUE))
flxDf.mod.min10 <- subset(flxDf.mod, flxDf.mod$qcomp.obs<=quantile(flxDf.mod$qcomp.obs,
0.10, na.rm=TRUE))
# FDCs
flxDf.mod <- CalcFdc(flxDf.mod, "qcomp.mod")
flxDf.mod <- CalcFdc(flxDf.mod, "qcomp.obs")
flxDf.mod.d <- CalcFdc(flxDf.mod.d, "qcomp.mod")
flxDf.mod.d <- CalcFdc(flxDf.mod.d, "qcomp.obs")
# Compile summary statistics
results["t_n"] <- length(flxDf.mod$qcomp.mod)
results["t_nse"] <- round(Nse(flxDf.mod$qcomp.mod, flxDf.mod$qcomp.obs), 2)
results["t_nselog"] <- round(NseLog(flxDf.mod$qcomp.mod, flxDf.mod$qcomp.obs), 2)
results["t_cor"] <- round(cor(flxDf.mod$qcomp.mod, flxDf.mod$qcomp.obs,
use="na.or.complete"), 2)
results["t_rmse"] <- round(Rmse(flxDf.mod$qcomp.mod, flxDf.mod$qcomp.obs), 4)
results["t_rmsenorm"] <- round(RmseNorm(flxDf.mod$qcomp.mod, flxDf.mod$qcomp.obs), 2)
results["t_bias"] <- round(sum(flxDf.mod$qcomp.mod-flxDf.mod$qcomp.obs, na.rm=T)/
sum(flxDf.mod$qcomp.obs, na.rm=TRUE) * 100, 1)
results["t_msd"] <- round(mean(flxDf.mod$qcomp.mod-flxDf.mod$qcomp.obs, na.rm=T), 4)
results["t_mae"] <- round(mean(abs(flxDf.mod$qcomp.mod-flxDf.mod$qcomp.obs), na.rm=T), 4)
results["t_errfdc"] <- round(integrate(splinefun(flxDf.mod[,"qcomp.mod.fdc"],
flxDf.mod[,"qcomp.mod"], method='natural'),
0.05, 0.95, subdivisions=2000)$value -
integrate(splinefun(flxDf.mod[,"qcomp.obs.fdc"],
flxDf.mod[,"qcomp.obs"], method='natural'),
0.05, 0.95, subdivisions=2000)$value, 2 )
results["t_errflash"] <- round(RBFlash(flxDf.mod$qcomp.obs) - RBFlash(flxDf.mod$qcomp.mod), 2)
results["dy_n"] <- length(flxDf.mod.d$qcomp.mod)
results["dy_nse"] <- round(Nse(flxDf.mod.d$qcomp.mod, flxDf.mod.d$qcomp.obs), 2)
results["dy_nselog"] <- round(NseLog(flxDf.mod.d$qcomp.mod, flxDf.mod.d$qcomp.obs), 2)
results["dy_cor"] <- round(cor(flxDf.mod.d$qcomp.mod, flxDf.mod.d$qcomp.obs,
use="na.or.complete"), 2)
results["dy_rmse"] <- round(Rmse(flxDf.mod.d$qcomp.mod, flxDf.mod.d$qcomp.obs), 4)
results["dy_rmsenorm"] <- round(RmseNorm(flxDf.mod.d$qcomp.mod, flxDf.mod.d$qcomp.obs), 2)
results["dy_bias"] <- round(sum(flxDf.mod.d$qcomp.mod-flxDf.mod.d$qcomp.obs, na.rm=T)/
sum(flxDf.mod.d$qcomp.obs, na.rm=TRUE) * 100, 1)
results["dy_msd"] <- round(mean(flxDf.mod.d$qcomp.mod-flxDf.mod.d$qcomp.obs, na.rm=T), 4)
results["dy_mae"] <- round(mean(abs(flxDf.mod.d$qcomp.mod-flxDf.mod.d$qcomp.obs), na.rm=T), 4)
results["dy_errcom"] <- round(mean(as.numeric(difftime(flxDf.mod.dcom$qcomp.mod,
flxDf.mod.dcom$qcomp.obs,
units="hours")), na.rm=T), 2)
results["dy_errmaxt"] <- round(mean(as.numeric(difftime(flxDf.mod.dmax$qcomp.mod,
flxDf.mod.dmax$qcomp.obs,
units="hours")), na.rm=T), 2)
results["dy_errfdc"] <- round(integrate(splinefun(flxDf.mod.d[,"qcomp.mod.fdc"],
flxDf.mod.d[,"qcomp.mod"],
method='natural'), 0.05, 0.95,
subdivisions=2000)$value -
integrate(splinefun(flxDf.mod.d[,"qcomp.obs.fdc"],
flxDf.mod.d[,"qcomp.obs"],
method='natural'), 0.05, 0.95,
subdivisions=2000)$value, 2 )
results["dy_errflash"] <- round(RBFlash(flxDf.mod.d$qcomp.obs) - RBFlash(flxDf.mod.d$qcomp.mod), 2)
results["mo_n"] <- length(flxDf.mod.mwy$qcomp.mod)
results["mo_nse"] <- round(Nse(flxDf.mod.mwy$qcomp.mod, flxDf.mod.mwy$qcomp.obs), 2)
results["mo_nselog"] <- round(NseLog(flxDf.mod.mwy$qcomp.mod, flxDf.mod.mwy$qcomp.obs), 2)
results["mo_cor"] <- round(cor(flxDf.mod.mwy$qcomp.mod, flxDf.mod.mwy$qcomp.obs,
use="na.or.complete"), 2)
results["mo_rmse"] <- round(Rmse(flxDf.mod.mwy$qcomp.mod, flxDf.mod.mwy$qcomp.obs), 4)
results["mo_rmsenorm"] <- round(RmseNorm(flxDf.mod.mwy$qcomp.mod, flxDf.mod.mwy$qcomp.obs), 2)
results["mo_bias"] <- round(sum(flxDf.mod.mwy$qcomp.mod-flxDf.mod.mwy$qcomp.obs, na.rm=T)/
sum(flxDf.mod.mwy$qcomp.obs, na.rm=TRUE) * 100, 1)
results["mo_msd"] <- round(mean(flxDf.mod.mwy$qcomp.mod-flxDf.mod.mwy$qcomp.obs, na.rm=T), 4)
results["mo_mae"] <- round(mean(abs(flxDf.mod.mwy$qcomp.mod-flxDf.mod.mwy$qcomp.obs), na.rm=T), 4)
results["mo_errcom"] <- round(mean(as.numeric(difftime(flxDf.mod.mwycom$qcomp.mod,
flxDf.mod.mwycom$qcomp.obs,
units="days")), na.rm=T), 2)
results["mo_errmaxt"] <- round(mean(as.numeric(difftime(flxDf.mod.mwymax$qcomp.mod,
flxDf.mod.mwymax$qcomp.obs,
units="days")), na.rm=T), 2)
results["yr_n"] <- length(flxDf.mod.wy$qcomp.mod)
results["yr_bias"] <- round(sum(flxDf.mod.wy$qcomp.mod-flxDf.mod.wy$qcomp.obs,
na.rm=T)/sum(flxDf.mod.wy$qcomp.obs, na.rm=TRUE) * 100, 1)
results["yr_msd"] <- round(mean(flxDf.mod.wy$qcomp.mod-flxDf.mod.wy$qcomp.obs, na.rm=T), 4)
results["yr_mae"] <- round(mean(abs(flxDf.mod.wy$qcomp.mod-flxDf.mod.wy$qcomp.obs), na.rm=T), 4)
results["yr_errcom"] <- round(mean(as.numeric(difftime(flxDf.mod.wycom$qcomp.mod,
flxDf.mod.wycom$qcomp.obs,
units="days")), na.rm=T), 2)
results["yr_errmaxt"] <- round(mean(as.numeric(difftime(flxDf.mod.wymax$qcomp.mod,
flxDf.mod.wymax$qcomp.obs,
units="days")), na.rm=T), 2)
results["max10_n"] <- length(flxDf.mod.max10$qcomp.mod)
results["max10_nse"] <- round(Nse(flxDf.mod.max10$qcomp.mod, flxDf.mod.max10$qcomp.obs), 2)
results["max10_nselog"] <- round(NseLog(flxDf.mod.max10$qcomp.mod, flxDf.mod.max10$qcomp.obs), 2)
results["max10_cor"] <- round(cor(flxDf.mod.max10$qcomp.mod, flxDf.mod.max10$qcomp.obs,
use="na.or.complete"), 2)
results["max10_rmse"] <- round(Rmse(flxDf.mod.max10$qcomp.mod, flxDf.mod.max10$qcomp.obs), 4)
results["max10_rmsenorm"] <- round(RmseNorm(flxDf.mod.max10$qcomp.mod,
flxDf.mod.max10$qcomp.obs), 2)
results["max10_bias"] <- round(sum(flxDf.mod.max10$qcomp.mod-flxDf.mod.max10$qcomp.obs,
na.rm=T)/sum(flxDf.mod.max10$qcomp.obs, na.rm=TRUE) * 100, 1)
results["max10_msd"] <- round(mean(flxDf.mod.max10$qcomp.mod-flxDf.mod.max10$qcomp.obs,
na.rm=T), 4)
results["max10_mae"] <- round(mean(abs(flxDf.mod.max10$qcomp.mod-flxDf.mod.max10$qcomp.obs),
na.rm=T), 4)
results["min10_n"] <- length(flxDf.mod.min10$qcomp.mod)
results["min10_nse"] <- round(Nse(flxDf.mod.min10$qcomp.mod, flxDf.mod.min10$qcomp.obs), 2)
results["min10_nselog"] <- round(NseLog(flxDf.mod.min10$qcomp.mod, flxDf.mod.min10$qcomp.obs), 2)
results["min10_cor"] <- round(cor(flxDf.mod.min10$qcomp.mod, flxDf.mod.min10$qcomp.obs,
use="na.or.complete"), 2)
results["min10_rmse"] <- round(Rmse(flxDf.mod.min10$qcomp.mod, flxDf.mod.min10$qcomp.obs), 4)
results["min10_rmsenorm"] <- round(RmseNorm(flxDf.mod.min10$qcomp.mod,
flxDf.mod.min10$qcomp.obs), 2)
results["min10_bias"] <- round(sum(flxDf.mod.min10$qcomp.mod-flxDf.mod.min10$qcomp.obs,
na.rm=T)/sum(flxDf.mod.min10$qcomp.obs, na.rm=TRUE) * 100, 1)
results["min10_msd"] <- round(mean(flxDf.mod.min10$qcomp.mod-flxDf.mod.min10$qcomp.obs,
na.rm=T), 4)
results["min10_mae"] <- round(mean(abs(flxDf.mod.min10$qcomp.mod-flxDf.mod.min10$qcomp.obs),
na.rm=T), 4)
# Attributes
attr(results, "descrip") <- c(t_n="sample size",
t_nse="nash-sutcliffe efficiency",
t_nselog="log nash-sutcliffe efficiency",
t_cor="correlation",
t_rmse="root mean squared error",
t_rmsenorm="normalized root mean squared error",
t_bias="bias",
t_msd="mean signed deviation",
t_mae="mean absolute error",
t_errfdc="integrated flow duration curve error",
t_errflash="flashiness index error",
dy_n="sample size",
dy_nse="nash-sutcliffe efficiency",
dy_nselog="log nash-sutcliffe efficiency",
dy_cor="correlation",
dy_rmse="root mean squared error",
dy_rmsenorm="normalized root mean squared error",
dy_bias="bias",
dy_msd="mean signed deviation",
dy_mae="mean absolute error",
dy_errcom="error in center-of-mass timing",
dy_errmaxt="error in max timing",
dy_errfdc="integrated flow duration curve error",
dy_errflash="flashiness index error",
mo_n="sample size",
mo_nse="nash-sutcliffe efficiency",
mo_nselog="log nash-sutcliffe efficiency",
mo_cor="correlation",
mo_rmse="root mean squared error",
mo_rmsenorm="normalized root mean squared error",
mo_bias="bias",
mo_msd="mean signed deviation",
mo_mae="mean absolute error",
mo_errcom="error in center-of-mass timing",
mo_errmaxt="error in max timing",
yr_n="sample size",
yr_bias="bias",
yr_msd="mean signed deviation",
yr_mae="mean absolute error",
yr_errcom="error in center-of-mass timing",
yr_errmaxt="error in max timing",
max10_n="sample size",
max10_nse="nash-sutcliffe efficiency",
max10_nselog="log nash-sutcliffe efficiency",
max10_cor="correlation",
max10_rmse="root mean squared error",
max10_rmsenorm="normalized root mean squared error",
max10_bias="bias",
max10_msd="mean signed deviation",
max10_mae="mean absolute error",
min10_n="sample size",
min10_nse="nash-sutcliffe efficiency",
min10_nselog="log nash-sutcliffe efficiency",
min10_cor="correlation",
min10_rmse="root mean squared error",
min10_rmsenorm="normalized root mean squared error",
min10_bias="bias",
min10_msd="mean signed deviation",
min10_mae="mean absolute error")
attr(results, "units") <- c(t_n="count",
t_nse="unitless (-Inf->1)",
t_nselog="unitless (-Inf->1)",
t_cor="unitless (-1->1)",
t_rmse="native flux units",
t_rmsenorm="unitless (0->Inf)",
t_bias="percent (%)",
t_msd="native flux units",
t_mae="native flux units",
t_errfdc="native flux units",
t_errflash="unitless",
dy_n="count",
dy_nse="unitless (-Inf->1)",
dy_nselog="unitless (-Inf->1)",
dy_cor="unitless (-1->1)",
dy_rmse="native flux units",
dy_rmsenorm="native flux units",
dy_bias="percent (%)",
dy_msd="native flux units",
dy_mae="native flux units",
dy_errcom="hours",
dy_errmaxt="hours",
dy_errfdc="native flux units",
dy_errflash="unitless",
mo_n="count",
mo_nse="unitless (-Inf->1)",
mo_nselog="unitless (-Inf->1)",
mo_cor="unitless (-1->1)",
mo_rmse="native flux units",
mo_rmsenorm="native flux units",
mo_bias="percent (%)",
mo_msd="native flux units",
mo_mae="native flux units",
mo_errcom="days",
mo_errmaxt="days",
yr_n="count",
yr_bias="percent (%)",
yr_msd="native flux units",
yr_mae="native flux units",
yr_errcom="days",
yr_errmaxt="days",
max10_n="sample size",
max10_nse="unitless (-Inf->1)",
max10_nselog="unitless (-Inf->1)",
max10_cor="unitless (-1->1)",
max10_rmse="native flux units",
max10_rmsenorm="native flux units",
max10_bias="percent (%)",
max10_msd="native flux units",
max10_mae="native flux units",
min10_n="sample size",
min10_nse="unitless (-Inf->1)",
min10_nselog="unitless (-Inf->1)",
min10_cor="unitless (-1->1)",
min10_rmse="native flux units",
min10_rmsenorm="native flux units",
min10_bias="percent (%)",
min10_msd="native flux units",
min10_mae="native flux units")
results
}
#' Computes flow duration curve statistics for WRF-Hydro streamflow output
#'
#' \code{CalcFdcPerf} calculates flow duration curve statistics for streamflow
#' output.
#'
#' \code{CalcFdcPerf} reads a model forecast point streamflow timeseries (i.e.,
#' created using \code{\link{ReadFrxstPts}}) and a streamflow observation
#' timeseries (i.e., created using \code{\link{ReadUsgsGage}}) and calculates
#' flow duration curve statistics at various exceedance thresholds (e.g., 10\%,
#' 20\%, etc.). The tool will subset data to matching time periods (e.g., if the
#' observed data is at 5-min increments and modelled data is at 1-hr increments,
#' the tool will subset the observed data to select only observations on the
#' matching hour break).
#'
#' Flow Duration Curve Statistics: \cr (mod = model output, obs = observations)
#' \itemize{ \item p.exceed: exceedance threshold (e.g., 0.2 means a flow value
#' that is exceeded 20\% of the time) \item q.mod: MODEL flow value at specified
#' exceedance threshold (in native flow units) \item q.obs: OBSERVED flow value
#' at specified exceedance threshold (in native flow units) \item q.err:
#' difference between model and observed flow values [mod-obs] (in native flow
#' units) \item q.perr: percent error in model flow [(mod-obs)/obs] }
#'
#' @param strDf.mod The forecast point output dataframe (required). Assumes only
#' one forecast point per file, so if you have multiple forecast points in
#' your output dataframe, use subset to isolate a single forecast point's
#' data. Also assumes model output and observation both contain POSIXct fields
#' (called "POSIXct").
#' @param strDf.obs The observed streamflow dataframe. Assumes only one gage per
#' file, so if you have multiple gages in your dataframe, use subset to
#' isolate a single gage's data. Also assumes model output and observation
#' both contain POSIXct fields (called "POSIXct").
#' @param strCol.mod The column name for the streamflow time series for the
#' MODEL data (default="q_cms")
#' @param strCol.obs The column name for the streamflow time series for the
#' OBSERVED data (default="q_cms")
#' @param stdate Start date for statistics (DEFAULT=NULL, all records will be
#' used). Date MUST be specified in POSIXct format with appropriate timezone
#' (e.g., as.POSIXct("2013-05-01 00:00:00", format="\%Y-\%m-\%d \%H:\%M:\%S",
#' tz="UTC"))
#' @param enddate End date for statistics (DEFAULT=NULL, all records will be
#' used). Date MUST be specified in POSIXct format with appropriate timezone
#' (e.g., as.POSIXct("2013-05-01 00:00:00", format="\%Y-\%m-\%d \%H:\%M:\%S",
#' tz="UTC"))
#' @return A new dataframe containing the flow duration curve statistics.
#'
#' @examples
#' ## Take forecast point model output for Fourmile Creek (modStrh.mod1.fc)
#' ## and a corresponding USGS gage observation file (obsStrh.fc), both at an
#' ## hourly time step, and calculate flow duration curve statistics. The
#' ## model forecast point data was imported using ReadFrxstPts and the gage
#' ## observation data was imported using ReadUsgsGage.
#'
#' \dontrun{
#' CalcFdcPerf(modStr1h.allrt.fc, obsStr5min.fc)
#'
#' Output:
#' p.exceed q.mod q.obs
#' 0.1 3.07 5.25
#' 0.2 1.35 2.31
#' 0.3 0.82 1.06
#' 0.4 0.48 0.65
#' 0.5 0.29 0.45
#' 0.6 0.18 0.34
#' 0.7 0.14 0.25
#' 0.8 0.11 0.19
#' 0.9 0.08 0.16
#' }
#' @keywords univar ts
#' @concept modelEval
#' @family modelEvaluation flowDurationCurves
#' @export
CalcFdcPerf <- function (strDf.mod, strDf.obs, strCol.mod="q_cms", strCol.obs="q_cms",
stdate=NULL, enddate=NULL) {
# Prepare data
if (!is.null(stdate) && !is.null(enddate)) {
strDf.obs <- subset(strDf.obs, POSIXct>=stdate & POSIXct<=enddate)
strDf.mod <- subset(strDf.mod, POSIXct>=stdate & POSIXct<=enddate)
}
strDf.mod <- CalcDates(strDf.mod)
strDf.mod$qcomp <- strDf.mod[,strCol.mod]
strDf.obs$qcomp <- strDf.obs[,strCol.obs]
if (as.integer(strDf.obs$POSIXct[2])-as.integer(strDf.obs$POSIXct[1]) >= 86400) {
strDf.obs$POSIXct=as.POSIXct(round(strDf.obs$POSIXct,"days"), tz="UTC")}
strDf.mod <- merge(strDf.mod, strDf.obs[c("POSIXct","qcomp")], by<-c("POSIXct"),
suffixes=c(".mod",".obs"))
#strDf.mod <- subset(strDf.mod, !is.na(strDf$qcomp.mod) & !is.na(strDf$qcomp.obs))
strDf.mod <- CalcFdc(strDf.mod, "qcomp.mod")
strDf.mod <- CalcFdc(strDf.mod, "qcomp.obs")
fdc.mod <- splinefun(strDf.mod[,"qcomp.mod.fdc"], strDf.mod[,"qcomp.mod"],
method='natural')
fdc.obs <- splinefun(strDf.mod[,"qcomp.obs.fdc"], strDf.mod[,"qcomp.obs"],
method='natural')
# Compile summary statistics
results <- as.data.frame(matrix(nrow = 9, ncol = 5))
colnames(results) <- c("p.exceed","q.mod","q.obs","q.err","q.perr")
results[, 1] <- seq(0.1, 0.9, 0.1)
for (i in 1:9) {
results[i, "q.mod"] <- round(fdc.mod(results[i,1]), 2)
results[i, "q.obs"] <- round(fdc.obs(results[i,1]), 2)
results[i, "q.err"] <- results[i, "q.mod"] - results[i, "q.obs"]
results[i, "q.perr"] <- round(results[i, "q.err"] / results[i, "q.obs"] * 100, 1)
}
results
}
mccreigh/rwrfhydro documentation built on Feb. 28, 2021, 1:53 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions CalcNoahmpWatBudg CalcNoahmpFluxes
Documented in CalcNoahmpFluxes CalcNoahmpWatBudg
#' Calculate water fluxes from NoahMP output
#'
#' \code{CalcNoahmpFluxes} calculates water balance fluxes from accumulated
#' water terms.
#'
#' Read a dataframe derived from NoahMP LDASOUT output (i.e., using
#' \code{\link{GetMultiNcdf}}) and calculate water budget component fluxes from
#' accumulated water variables. NOTE: Currently only works for runs using NoahMP
#' as the LSM.
#'
#' @param ldasoutDf The LDASOUT dataframe
#' @param idCol The optional ID column to use to parse the output;
#' for example, a basin ID column or other index ID. (DEFAULT=NULL)
#' @return The input dataframe with new water flux columns added.
#'
#' @examples
#' ## Take a NoahMP LDASOUT dataframe for a model run of Fourmile Creek
#' ## and generate a dataframe with water fluxes added.
#' \dontrun{
#' modLDASOUT1d.nort.fc <- CalcNoahmpFluxes(modLDASOUT1d.nort.fc)
#' }
#' @keywords manip
#' @concept dataMgmt
#' @family modelEvaluation
#' @export
CalcNoahmpFluxes <- function(ldasoutDf, idCol=NULL) {
ldasoutDf <- ldasoutDf[order(ldasoutDf[,"POSIXct"]),]
if (is.null(idCol)) {
if ("ACCPRCP" %in% colnames(ldasoutDf)) {
ldasoutDf$DEL_ACCPRCP[2:nrow(ldasoutDf)] <- diff(ldasoutDf$ACCPRCP) }
if ("ACCECAN" %in% colnames(ldasoutDf)) {
ldasoutDf$DEL_ACCECAN[2:nrow(ldasoutDf)] <- diff(ldasoutDf$ACCECAN) }
if ("ACCETRAN" %in% colnames(ldasoutDf)) {
ldasoutDf$DEL_ACCETRAN[2:nrow(ldasoutDf)] <- diff(ldasoutDf$ACCETRAN) }
if ("ACCEDIR" %in% colnames(ldasoutDf)) {
ldasoutDf$DEL_ACCEDIR[2:nrow(ldasoutDf)] <- diff(ldasoutDf$ACCEDIR) }
if ("ACCET" %in% colnames(ldasoutDf)) {
ldasoutDf$DEL_ACCET[2:nrow(ldasoutDf)] <- diff(ldasoutDf$ACCET) }
if ("UGDRNOFF" %in% colnames(ldasoutDf)) {
ldasoutDf$DEL_UGDRNOFF[2:nrow(ldasoutDf)] <- diff(ldasoutDf$UGDRNOFF) }
if ("SFCRNOFF" %in% colnames(ldasoutDf)) {
ldasoutDf$DEL_SFCRNOFF[2:nrow(ldasoutDf)] <- diff(ldasoutDf$SFCRNOFF) }
if ("ACSNOM" %in% colnames(ldasoutDf)) {
ldasoutDf$DEL_ACSNOM[2:nrow(ldasoutDf)] <- diff(ldasoutDf$ACSNOM) }
if ("ACSNOW" %in% colnames(ldasoutDf)) {
ldasoutDf$DEL_ACSNOW[2:nrow(ldasoutDf)] <- diff(ldasoutDf$ACSNOW) }
if ("SNEQV" %in% colnames(ldasoutDf)) {
ldasoutDf$DEL_SNEQV[2:nrow(ldasoutDf)] <- diff(ldasoutDf$SNEQV) }
if ("SNOWH" %in% colnames(ldasoutDf)) {
ldasoutDf$DEL_SNOWH[2:nrow(ldasoutDf)] <- diff(ldasoutDf$SNOWH) }
} else {
idList <- unique(ldasoutDf[,idCol])
for (i in 1:length(idList)) {
if ("ACCPRCP" %in% colnames(ldasoutDf)) {
tmp <- subset(ldasoutDf$ACCPRCP, ldasoutDf[,idCol]==idList[i])
tmp <- c(NA, diff(tmp))
ldasoutDf$DEL_ACCPRCP[ldasoutDf[,idCol]==idList[i]] <- tmp }
if ("ACCECAN" %in% colnames(ldasoutDf)) {
tmp <- subset(ldasoutDf$ACCECAN, ldasoutDf[,idCol]==idList[i])
tmp <- c(NA, diff(tmp))
ldasoutDf$DEL_ACCECAN[ldasoutDf[,idCol]==idList[i]] <- tmp }
if ("ACCEDIR" %in% colnames(ldasoutDf)) {
tmp <- subset(ldasoutDf$ACCEDIR, ldasoutDf[,idCol]==idList[i])
tmp <- c(NA, diff(tmp))
ldasoutDf$DEL_ACCEDIR[ldasoutDf[,idCol]==idList[i]] <- tmp }
if ("ACCETRAN" %in% colnames(ldasoutDf)) {
tmp <- subset(ldasoutDf$ACCETRAN, ldasoutDf[,idCol]==idList[i])
tmp <- c(NA, diff(tmp))
ldasoutDf$DEL_ACCETRAN[ldasoutDf[,idCol]==idList[i]] <- tmp }
if ("ACCET" %in% colnames(ldasoutDf)) {
tmp <- subset(ldasoutDf$ACCET, ldasoutDf[,idCol]==idList[i])
tmp <- c(NA, diff(tmp))
ldasoutDf$DEL_ACCET[ldasoutDf[,idCol]==idList[i]] <- tmp }
if ("SFCRNOFF" %in% colnames(ldasoutDf)) {
tmp <- subset(ldasoutDf$SFCRNOFF, ldasoutDf[,idCol]==idList[i])
tmp <- c(NA, diff(tmp))
ldasoutDf$DEL_SFCRNOFF[ldasoutDf[,idCol]==idList[i]] <- tmp }
if ("UGDRNOFF" %in% colnames(ldasoutDf)) {
tmp <- subset(ldasoutDf$UGDRNOFF, ldasoutDf[,idCol]==idList[i])
tmp <- c(NA, diff(tmp))
ldasoutDf$DEL_UGDRNOFF[ldasoutDf[,idCol]==idList[i]] <- tmp }
if ("ACSNOM" %in% colnames(ldasoutDf)) {
tmp <- subset(ldasoutDf$ACSNOM, ldasoutDf[,idCol]==idList[i])
tmp <- c(NA, diff(tmp))
ldasoutDf$DEL_ACSNOM[ldasoutDf[,idCol]==idList[i]] <- tmp }
if ("ACSNOW" %in% colnames(ldasoutDf)) {
tmp <- subset(ldasoutDf$ACSNOW, ldasoutDf[,idCol]==idList[i])
tmp <- c(NA, diff(tmp))
ldasoutDf$DEL_ACSNOW[ldasoutDf[,idCol]==idList[i]] <- tmp }
if ("SNEQV" %in% colnames(ldasoutDf)) {
tmp <- subset(ldasoutDf$SNEQV, ldasoutDf[,idCol]==idList[i])
tmp <- c(NA, diff(tmp))
ldasoutDf$DEL_SNEQV[ldasoutDf[,idCol]==idList[i]] <- tmp }
if ("SNOWH" %in% colnames(ldasoutDf)) {
tmp <- subset(ldasoutDf$SNOWH, ldasoutDf[,idCol]==idList[i])
tmp <- c(NA, diff(tmp))
ldasoutDf$DEL_SNOWH[ldasoutDf[,idCol]==idList[i]] <- tmp }
}
ldasoutDf <- ldasoutDf[order(ldasoutDf[,idCol], ldasoutDf[,"POSIXct"]),]
}
ldasoutDf
}
#' Calculate water balance from WRF-Hydro (w/NoahMP) output
#'
#' \code{CalcNoahmpWatBudg} calculates water budget components from WRF-Hydro
#' (w/NoahMP) model output.
#'
#' \code{CalcNoahmpWatBudg} reads dataframes derived from WRF-Hydro output
#' (i.e., using \code{\link{GetMultiNcdf}}) and calculates water budget
#' partitioning (e.g., surface runoff, evaporation, groundwater). Assumes
#' WRF-Hydro output dataframes have already been masked to the desired basin.
#' See \code{\link{GetMultiNcdf}} documentation for examples of how to do this.
#' NOTE: Currently only works for model runs using NoahMP as the LSM.
#'
#' REQUIRED variables (these terms must be in your input dataframes): \itemize{
#' \item LDASOUT: ACCPRCP, ACCECAN, ACCETRAN, ACCEDIR, SFCRNOFF, UGDRNOFF,
#' SOIL_M (all layers), SNEQV, CANICE, CANLIQ \item RTOUT (optional, use if
#' overland or subsurface routing were activated): QSTRMVOLRT, SFCHEADSUBRT,
#' QBDRY \item GWOUT (optional, use if groundwater bucket model was activated):
#' q_cms, POSIXct }
#'
#' OUTPUT water budget terms (may vary depending on model configuration):
#' \itemize{ \item LSM_PRCP: Total precipitation (mm) \item LSM_ECAN: Total
#' canopy evaporation (mm) \item LSM_ETRAN: Total transpiration (mm) \item
#' LSM_EDIR: Total surface evaporation (mm) \item LSM_DELSWE: Change in snowpack
#' snow water equivalent (mm) \item LSM_DELCANWAT: Change in canopy water
#' storage (liquid + ice) (mm) \item LSM_SFCRNOFF: Surface runoff from LSM
#' \emph{(for an LSM-only run)} (mm) \item LSM_UGDRNOFF: Subsurface runoff from
#' LSM \emph{(for an LSM-only run)} (mm) \item LSM_DELSOILM: Change in total
#' soil moisture storage (mm) \item HYD_QSTRMVOL: Total runoff into channel from
#' land \emph{(routing model only)} (mm) \item HYD_DELSFCHEAD: Change in
#' surface storage \emph{(routing model only)} (mm) \item HYD_QBDRY: Total flow
#' outside of domain \emph{(routing model only)} (mm) \item HYD_GWOUT: Total
#' groundwater outflow \emph{(routing model only)} (mm) \item
#' WB_SFCRNOFF: Total surface runoff used in the water budget calculation
#' (\emph{either} LSM_SFCRNOFF or HYD_QSTRMVOL) (mm) \item WB_GWOUT: Total
#' groundwater outflow used in the water budget calculation (\emph{either}
#' LSM_UGDRNOFF or HYD_GWOUT) (mm) \item WB_DELGWSTOR: Change in groundwater
#' storage (adjusted by surface runoff when no surface routing active) (mm)
#' \item ERROR: Remainder in water budget (mm)
#' \item RUN_FRAC: Runoff fraction, runoff/precipitation \item EVAP_FRAC:
#' Evaporative fraction, evapotranspiration/precipitation \item STOR_FRAC:
#' Change in storage fraction, storagechange/precipitation }
#'
#' @param ldasoutDf The LDASOUT dataframe (required)
#' @param rtoutDf The RTOUT dataframe, if overland or subsurface routing was
#' turned on (default=NULL)
#' @param gwoutDf The GW_OUT dataframe, if groundwater model was turned on
#' (default=NULL)
#' @param sfcrt A flag whether surface overland flow routing was active. All
#' other routing options are determined based on the input dataframes, as
#' needed (e.g., if gwoutDf is provided, it is assumed that the groundwater
#' model was active). (default=FALSE)
#' @param soildeps A list of soil layer depths in mm (top to bottom,
#' default=c(100, 300, 600, 1000))
#' @param basarea The basin area in square km (necessary only if gwoutDf is
#' provided)
#' @return A new dataframe containing the water budget components in mm.
#'
#' @examples
#' ## Take a NoahMP LDASOUT dataframe for a model run of Fourmile Creek with no routing
#' ## options turned on and return a water budget summary.
#'
#' \dontrun{
#' wb.nort.fc <- CalcNoahmpWatBudg(modLDASOUT1d.nort.fc)
#' wb.nort.fc
#' ##
#' ## Take NoahMP LDASOUT, HYDRO model RTOUT, and GW_outflow dataframes for a model
#' ## run of Fourmile Creek with subsurface, overland, and groundwater routing options
#' ## turned on and return a water budget summary. The default soil depths were used
#' ## and the basin is 63.1 km2. NOTE: We MUST specify with the sfcrt flag that overland
#' ## flow routing was turned on. Otherwise the LSM surface runoff term is used.
#'
#' wb.allrt.fc <- CalcNoahmpWatBudg(modLDASOUT1d.allrt.fc,
#' modRTOUT1h.allrt.fc,
#' modGWOUT1h.allrt.fc,
#' sfcrt=TRUE, basarea=63.1)
#' wb.allrt.fc
#' }
#' @keywords manip
#' @concept modelEval
#' @family modelEvaluation
#' @export
CalcNoahmpWatBudg <- function(ldasoutDf, rtoutDf=NULL, gwoutDf=NULL, sfcrt=FALSE,
soildeps=c(100,300,600,1000), basarea=NULL) {
wbDf <- as.data.frame(t(as.matrix(rep(0, 6))))
colnames(wbDf) <- c("LSM_PRCP","LSM_ECAN","LSM_ETRAN","LSM_EDIR",
"LSM_DELSWE","LSM_DELCANWAT")
# LSM water fluxes for all model cases
wbDf[1,1] <- ldasoutDf$ACCPRCP[nrow(ldasoutDf)]-ldasoutDf$ACCPRCP[1]
wbDf[1,2] <- ldasoutDf$ACCECAN[nrow(ldasoutDf)]-ldasoutDf$ACCECAN[1]
wbDf[1,3] <- ldasoutDf$ACCETRAN[nrow(ldasoutDf)]-ldasoutDf$ACCETRAN[1]
wbDf[1,4] <- ldasoutDf$ACCEDIR[nrow(ldasoutDf)]-ldasoutDf$ACCEDIR[1]
wbDf[1,5] <- ldasoutDf$SNEQV[nrow(ldasoutDf)]-ldasoutDf$SNEQV[1]
wbDf[1,6] <- (ldasoutDf$CANICE[nrow(ldasoutDf)] + ldasoutDf$CANLIQ[nrow(ldasoutDf)]) -
(ldasoutDf$CANICE[1] + ldasoutDf$CANLIQ[1])
# LSM overlap water fluxes
wbDf[1,"LSM_SFCRNOFF"] <- ldasoutDf$SFCRNOFF[nrow(ldasoutDf)]-ldasoutDf$SFCRNOFF[1]
wbDf[1,"LSM_UGDRNOFF"] <- ldasoutDf$UGDRNOFF[nrow(ldasoutDf)]-ldasoutDf$UGDRNOFF[1]
numsoil <- length(soildeps)
soilm <- 0.0
for (i in 1:numsoil) {
soilm <- soilm + (ldasoutDf[nrow(ldasoutDf),paste0("SOIL_M",i)]-
ldasoutDf[1,paste0("SOIL_M",i)])*soildeps[i]
}
wbDf[1,"LSM_DELSOILM"] <- soilm
# HYDRO surface/subsurface water fluxes
if (!is.null(rtoutDf)) {
wbDf[1,"HYD_QSTRMVOL"] <- rtoutDf$QSTRMVOLRT[nrow(rtoutDf)] -
rtoutDf$QSTRMVOLRT[1]
wbDf[1,"HYD_DELSFCHEAD"] <- rtoutDf$SFCHEADSUBRT[nrow(rtoutDf)] -
rtoutDf$SFCHEADSUBRT[1]
wbDf[1,"HYD_QBDRY"] <- -(rtoutDf$QBDRYRT[nrow(rtoutDf)] - rtoutDf$QBDRYRT[1])
}
else {
message('Message: No routing output dataframe (rtoutDf) was provided. Proceeding using LSM surface runoff.')
wbDf[1,"HYD_QSTRMVOL"] <- NA
wbDf[1,"HYD_DELSFCHEAD"] <- NA
wbDf[1,"HYD_QBDRY"] <- NA
}
# HYDRO groundwater fluxes
if (!is.null(gwoutDf)) {
if (!is.null(basarea)) {
dt <- as.integer(difftime(gwoutDf$POSIXct[2],gwoutDf$POSIXct[1],units="secs"))
wbDf[1,"HYD_GWOUT"] <- sum(gwoutDf$q_cms/(basarea*1000*1000)*1000*dt, na.rm=T)
# wbDf[1,"HYD_DELGWSTOR"] <- (ldasoutDf$UGDRNOFF[nrow(ldasoutDf)] -
# ldasoutDf$UGDRNOFF[1]) -
# wbDf$HYD_GWOUT
}
else { stop('Error: Groundwater outflow dataframe (gwoutDf) provided but no basin area (basarea). Please provide the basin area.') }
}
else {
message('Message: No groundwater outflow dataframe (gwoutDf) was provided. Proceeding using LSM underground runoff.')
wbDf[1,"HYD_GWOUT"] <- NA
# wbDf[1,"HYD_DELGWSTOR"] <- NA
}
# Overland routing check
if ( sfcrt & is.null(rtoutDf)) {
stop('Error: Surface overland routing (sfcrt) is active (TRUE) but no routing output dataframe (rtoutDf) was specified.')
}
if ( !sfcrt & !is.null(rtoutDf) ) {
message('Message: A routing output file (rtoutDf) was provided but surface overland routing (sfcrt) is inactive (FALSE). Proceeding using LSM surface runoff.')
}
wbDf[1,"WB_SFCRNOFF"] <- ifelse(sfcrt, wbDf[1, "HYD_QSTRMVOL"], wbDf[1, "LSM_SFCRNOFF"])
# Groundwater routing check
wbDf[1,"WB_GWOUT"] <- if (!is.null(gwoutDf)) {
if (!sfcrt) { wbDf[1, "HYD_GWOUT"] - wbDf[1, "LSM_SFCRNOFF"] }
else { wbDf[1, "HYD_GWOUT"] }
}
else { wbDf[1, "LSM_UGDRNOFF"] }
wbDf[1,"WB_DELGWSTOR"] <- (ldasoutDf$UGDRNOFF[nrow(ldasoutDf)]-ldasoutDf$UGDRNOFF[1]) -
wbDf$WB_GWOUT
# Water budget error
wbDf[1,"ERROR"] <- with( wbDf, LSM_PRCP - LSM_ECAN - LSM_ETRAN - LSM_EDIR -
WB_SFCRNOFF -
ifelse(is.na(HYD_QBDRY), 0.0, HYD_QBDRY) -
WB_GWOUT - LSM_DELSOILM -
LSM_DELSWE - LSM_DELCANWAT -
ifelse(is.na(HYD_DELSFCHEAD), 0.0, HYD_DELSFCHEAD) -
ifelse(is.na(WB_DELGWSTOR), 0.0, WB_DELGWSTOR) )
# Calculate various fractions
wbDf[1,"RUN_FRAC"] <- with( wbDf, (WB_SFCRNOFF +
ifelse(is.na(HYD_QBDRY), 0.0, HYD_QBDRY) +
WB_GWOUT) / LSM_PRCP )
wbDf[1,"EVAP_FRAC"] <- with( wbDf, (LSM_ECAN + LSM_ETRAN + LSM_EDIR) / LSM_PRCP )
wbDf[1,"STOR_FRAC"] <- with( wbDf, (LSM_DELSOILM + LSM_DELSWE + LSM_DELCANWAT +
ifelse(is.na(HYD_DELSFCHEAD), 0.0, HYD_DELSFCHEAD) +
ifelse(is.na(WB_DELGWSTOR), 0.0, WB_DELGWSTOR) ) /
LSM_PRCP )
wbDf
}
#' Computes model performance statistics for WRF-Hydro flux output
#'
#' \code{CalcModPerf} calculates model performance statistics for flux output.
#'
#' \code{CalcModPerf} reads a model flux time series (i.e., created using
#' \code{\link{ReadFrxstPts}}) and an observation time series (i.e., created
#' using \code{\link{ReadUsgsGage}}) and calculates model performance statistics
#' (Nash-Sutcliffe Efficiency, Rmse, etc.) at various time scales and for low
#' and high fluxes. The tool will subset data to matching time periods (e.g., if
#' the observed data is at 5-min increments and modelled data is at 1-hr
#' increments, the tool will subset the observed data to select only
#' observations on the matching hour break).
#'
#' Performance Statistics: \cr (mod = model output, obs = observations, n =
#' sample size) \itemize{ \item n: sample size \item nse: Nash-Sutcliffe
#' Efficiency \deqn{nse = 1 - ( sum((obs - mod)^2) / sum((obs - mean(obs))^2) )
#' } \item nselog: log-transformed Nash-Sutcliffe Efficiency \deqn{nselog = 1 -
#' ( sum((log(obs) - log(mod))^2) / sum((log(obs) - mean(log(obs)))^2) ) } \item
#' cor: correlation coefficient \deqn{cor = cor(mod, obs) } \item rmse: root
#' mean squared error \deqn{rmse = sqrt( sum((mod - obs)^2) / n ) } \item
#' rmsenorm: normalized root mean squared error \deqn{rmsenorm = rmse /
#' (max(obs) - min(obs)) } \item bias: percent bias \deqn{bias = sum(mod - obs)
#' / sum(obs) * 100 } \item msd: mean signed deviation \deqn{msd = mean(mod -
#' obs) } \item mae: mean absolute error \deqn{mae = mean(abs(mod -
#' obs)) } \item errcom: error in the center-of-mass of the flux, where
#' center-of-mass is the hour/day when 50\% of daily/monthly/water-year flux has
#' occurred. Reported as number of hours for daily time scale and number of days
#' for monthly and yearly time scales. \item errmaxt: Error in the time of
#' maximum flux. Reported as number of hours for daily time scale and number of
#' days for monthly and yearly time scales). \item errfdc: Error in the
#' integrated flow duration curve between 0.05 and 0.95 exceedance thresholds
#' (in native flow units). }
#'
#' Time scales/Flux types: \itemize{ \item ts = native model/observation time
#' step (e.g., hourly) \item daily = daily time step \item monthly = monthly
#' time step \item yearly = water-year time step \item max10 = high flows;
#' restricted to the portion of the time series where the observed flux is in
#' the highest 10\% (native time step) \item min10 = low flows; restricted to
#' the portion of the time series where the observed flux is in the lowest 10\%
#' (native time step) }
#'
#' @param flxDf.mod The flux output dataframe (required). Assumes only one
#' forecast point per file, so if you have multiple forecast points in your
#' output dataframe, use subset to isolate a single forecast point's data.
#' Also assumes model output and observation both contain POSIXct fields
#' (called "POSIXct").
#' @param flxDf.obs The observed flux dataframe. Assumes only one observation
#' point per file, so if you have multiple observation points in your
#' dataframe, use subset to isolate a single point's data. Also assumes model
#' output and observation both contain POSIXct fields (called "POSIXct").
#' @param flxCol.mod The column name for the flux time series for the MODEL data
#' (default="q_cms")
#' @param flxCol.obs The column name for the flux time series for the OBSERVED
#' data (default="q_cms")
#' @param stdate Start date for statistics (DEFAULT=NULL, all records will be
#' used). Date MUST be specified in POSIXct format with appropriate timezone
#' (e.g., as.POSIXct("2013-05-01 00:00:00", format="\%Y-\%m-\%d \%H:\%M:\%S",
#' tz="UTC"))
#' @param enddate End date for statistics (DEFAULT=NULL, all records will be
#' used). Date MUST be specified in POSIXct format with appropriate timezone
#' (e.g., as.POSIXct("2013-05-01 00:00:00", format="\%Y-\%m-\%d \%H:\%M:\%S",
#' tz="UTC"))
#' @param subdivisions Number of subdivisions used in flow duration curve
#' integration (DEFAULT=1000); increase value if integrate function throws an
#' error.
#' @return A new dataframe containing the model performance statistics.
#'
#' @examples
#' ## Take forecast point model output for Fourmile Creek (modStrh.mod1.fc) and a
#' ## corresponding USGS gage observation file (obsStrh.fc), both at an hourly time step,
#' ## and calculate model performance statistics. The model forecast point data was
#' ## imported using ReadFrxstPts and the gage observation data was imported using
#' ## ReadUsgsGage.
#'
#'
#' \dontrun{
#' CalcModPerf(modStrd.chrt.fc, obsStr5min.fc)
#'
#' > Output:
#' n nse nselog cor rmse rmsenorm bias msd mae errcom errmaxt errfdc
#' units count unitless unitless unitless flux units unitless percent flux units flux units hours|days hours|days flux units
#' t 116 0.81 0.67 0.91 0.211 8.48 7.1 0.0313 0.1385 <NA> <NA> 0.04
#' daily 116 0.81 0.67 0.91 0.211 8.48 7.1 0.0313 0.1385 0 0 0.04
#' monthly 5 0.92 0.82 0.96 0.1073 10.35 9.1 0.0324 0.0833 -1.2 0.2 <NA>
#' yearly 1 <NA> <NA> <NA> <NA> <NA> 7.1 0.0313 0.0313 0 4 <NA>
#' max10 12 0.46 0.58 0.71 0.348 23.63 -2.8 -0.044 0.248 <NA> <NA> <NA>
#' min10 12 -385.99 -14.59 0.35 0.1678 538.78 112.4 0.056 0.0592 <NA> <NA> <NA>
#' }
#' @keywords univar ts
#' @concept modelEval
#' @family modelEvaluation
#' @export
CalcModPerf <- function (flxDf.mod, flxDf.obs, flxCol.mod="q_cms", flxCol.obs="q_cms",
stdate=NULL, enddate=NULL, subdivisions=1000) {
# Internal functions
which.max.dt <- function(dd, qcol, dtcol) { dd[which.max(dd[,qcol]), dtcol] }
CalcCOM.dt <- function(dd, qcol, dtcol) { dd[CalcCOM(dd[,qcol]), dtcol] }
# Prepare data
if (!is.null(stdate) & is.null(enddate)) {
flxDf.obs <- subset(flxDf.obs, POSIXct>=stdate)
flxDf.mod <- subset(flxDf.mod, POSIXct>=stdate)
}
if (is.null(stdate) & !is.null(enddate)) {
flxDf.obs <- subset(flxDf.obs, POSIXct<=enddate)
flxDf.mod <- subset(flxDf.mod, POSIXct<=enddate)
}
if (!is.null(stdate) & !is.null(enddate)) {
flxDf.obs <- subset(flxDf.obs, POSIXct>=stdate & POSIXct<=enddate)
flxDf.mod <- subset(flxDf.mod, POSIXct>=stdate & POSIXct<=enddate)
}
flxDf.mod$qcomp <- flxDf.mod[,flxCol.mod]
flxDf.obs$qcomp <- flxDf.obs[,flxCol.obs]
# model timestep in secs
modT <- as.integer(flxDf.mod$POSIXct[2])-as.integer(flxDf.mod$POSIXct[1])
#if (as.integer(flxDf.obs$POSIXct[2])-as.integer(flxDf.obs$POSIXct[1]) >= 86400) {
# flxDf.obs$POSIXct=as.POSIXct(round(flxDf.obs$POSIXct,"days"), tz="UTC")}
flxDf.mod <- merge(flxDf.mod[c("POSIXct","qcomp")], flxDf.obs[c("POSIXct","qcomp")],
by<-c("POSIXct"), suffixes=c(".mod",".obs"))
flxDf.mod <- subset(flxDf.mod, !is.na(flxDf.mod$qcomp.mod) & !is.na(flxDf.mod$qcomp.obs))
flxDf.mod <- CalcDates(flxDf.mod)
flxDf.mod$date <- as.POSIXct(trunc(flxDf.mod$POSIXct, "days"))
results <- as.data.frame(matrix(nrow = 7, ncol = 12))
colnames(results) = c("n", "nse", "nselog", "cor", "rmse", "rmsenorm", "bias", "msd",
"mae", "errcom", "errmaxt", "errfdc")
rownames(results) = c("units", "t", "daily", "monthly", "yearly", "max10", "min10")
exclvars <- names(flxDf.mod) %in% c("POSIXct", "secs", "timest", "date", "stat")
# Base aggregations
flxDf.mod.d <- aggregate(flxDf.mod[!exclvars], by = list(flxDf.mod$date),
CalcMeanNarm)
flxDf.mod.d$mwy <- paste0(flxDf.mod.d$month, "-", flxDf.mod.d$wy)
flxDf.mod.mwy <- aggregate(flxDf.mod[c("qcomp.mod","qcomp.obs")],
by = list(flxDf.mod$month, flxDf.mod$wy), CalcMeanNarm)
flxDf.mod.wy <- aggregate(flxDf.mod[c("qcomp.mod","qcomp.obs")],
by = list(flxDf.mod$wy), CalcMeanNarm)
# Time of center of mass aggregations
tmp.mod <- plyr::ddply(flxDf.mod, plyr::.(date), CalcCOM.dt, qcol="qcomp.mod",
dtcol="POSIXct")
tmp.obs <- plyr::ddply(flxDf.mod, plyr::.(date), CalcCOM.dt, qcol="qcomp.obs",
dtcol="POSIXct")
flxDf.mod.dcom <- merge(tmp.mod, tmp.obs, by=c("date"))
colnames(flxDf.mod.dcom) <- c("date", "qcomp.mod", "qcomp.obs")
tmp.mod <- plyr::ddply(flxDf.mod.d, plyr::.(month, wy),
CalcCOM.dt, qcol="qcomp.mod", dtcol="Group.1")
tmp.obs <- plyr::ddply(flxDf.mod.d, plyr::.(month, wy),
CalcCOM.dt, qcol="qcomp.obs", dtcol="Group.1")
flxDf.mod.mwycom <- merge(tmp.mod, tmp.obs, by=c("month", "wy"))
colnames(flxDf.mod.mwycom) <- c("month", "wy", "qcomp.mod", "qcomp.obs")
tmp.mod <- plyr::ddply(flxDf.mod.d, plyr::.(wy),
CalcCOM.dt, qcol="qcomp.mod", dtcol="Group.1")
tmp.obs <- plyr::ddply(flxDf.mod.d, plyr::.(wy),
CalcCOM.dt, qcol="qcomp.obs", dtcol="Group.1")
flxDf.mod.wycom <- merge(tmp.mod, tmp.obs, by=c("wy"))
colnames(flxDf.mod.wycom) <- c("wy", "qcomp.mod", "qcomp.obs")
# Time of max aggregations
tmp.mod <- plyr::ddply(flxDf.mod, plyr::.(date),
which.max.dt, qcol="qcomp.mod", dtcol="POSIXct")
tmp.obs <- plyr::ddply(flxDf.mod, plyr::.(date),
which.max.dt, qcol="qcomp.obs", dtcol="POSIXct")
flxDf.mod.dmax <- merge(tmp.mod, tmp.obs, by=c("date"))
colnames(flxDf.mod.dmax) <- c("date", "qcomp.mod", "qcomp.obs")
tmp.mod <- plyr::ddply(flxDf.mod.d, plyr::.(month, wy),
which.max.dt, qcol="qcomp.mod", dtcol="Group.1")
tmp.obs <- plyr::ddply(flxDf.mod.d, plyr::.(month, wy),
which.max.dt, qcol="qcomp.obs", dtcol="Group.1")
flxDf.mod.mwymax <- merge(tmp.mod, tmp.obs, by=c("month", "wy"))
colnames(flxDf.mod.mwymax) <- c("month", "wy", "qcomp.mod", "qcomp.obs")
tmp.mod <- plyr::ddply(flxDf.mod.d, plyr::.(wy),
which.max.dt, qcol="qcomp.mod", dtcol="Group.1")
tmp.obs <- plyr::ddply(flxDf.mod.d, plyr::.(wy),
which.max.dt, qcol="qcomp.obs", dtcol="Group.1")
flxDf.mod.wymax <- merge(tmp.mod, tmp.obs, by=c("wy"))
colnames(flxDf.mod.wymax) <- c("wy", "qcomp.mod", "qcomp.obs")
# NAs cause which.max to return a list, so force back to ints
# (out for now since we pre-screen for NAs)
#flxDf.mod.dmax$qcomp.mod <- as.integer(flxDf.mod.dmax$qcomp.mod)
#flxDf.mod.dmax$qcomp.obs <- as.integer(flxDf.mod.dmax$qcomp.obs)
#flxDf.mod.mwymax$qcomp.mod <- as.integer(flxDf.mod.mwymax$qcomp.mod)
#flxDf.mod.mwymax$qcomp.obs <- as.integer(flxDf.mod.mwymax$qcomp.obs)
#flxDf.mod.wymax$qcomp.mod <- as.integer(flxDf.mod.wymax$qcomp.mod)
#flxDf.mod.wymax$qcomp.obs <- as.integer(flxDf.mod.wymax$qcomp.obs)
# Mins and Maxes
flxDf.mod.max10 <- subset(flxDf.mod, flxDf.mod$qcomp.obs >=
quantile(flxDf.mod$qcomp.obs, 0.90, na.rm=TRUE))
flxDf.mod.min10 <- subset(flxDf.mod, flxDf.mod$qcomp.obs <=
quantile(flxDf.mod$qcomp.obs, 0.10, na.rm=TRUE))
# FDCs
flxDf.mod <- CalcFdc(flxDf.mod, "qcomp.mod")
flxDf.mod <- CalcFdc(flxDf.mod, "qcomp.obs")
flxDf.mod.d <- CalcFdc(flxDf.mod.d, "qcomp.mod")
flxDf.mod.d <- CalcFdc(flxDf.mod.d, "qcomp.obs")
# Compile summary statistics
results["t", "n"] <- length(flxDf.mod$qcomp.mod)
results["t", "nse"] <- round(Nse(flxDf.mod$qcomp.mod, flxDf.mod$qcomp.obs), 2)
results["t", "nselog"] <- round(NseLog(flxDf.mod$qcomp.mod, flxDf.mod$qcomp.obs), 2)
results["t", "cor"] <- round(cor(flxDf.mod$qcomp.mod, flxDf.mod$qcomp.obs,
use="na.or.complete"), 2)
results["t", "rmse"] <- round(Rmse(flxDf.mod$qcomp.mod, flxDf.mod$qcomp.obs), 4)
results["t", "rmsenorm"] <- round(RmseNorm(flxDf.mod$qcomp.mod, flxDf.mod$qcomp.obs), 2)
results["t", "bias"] <- round(sum(flxDf.mod$qcomp.mod-flxDf.mod$qcomp.obs, na.rm=T)/
sum(flxDf.mod$qcomp.obs, na.rm=TRUE) * 100, 1)
results["t", "msd"] <- round(mean(flxDf.mod$qcomp.mod-flxDf.mod$qcomp.obs, na.rm=T), 4)
results["t", "mae"] <- round(mean(abs(flxDf.mod$qcomp.mod-flxDf.mod$qcomp.obs),
na.rm=T), 4)
results["t", "errcom"] <- NA
results["t", "errmaxt"] <- NA
results["t", "errfdc"] <- round(integrate(splinefun(flxDf.mod[,"qcomp.mod.fdc"],
flxDf.mod[,"qcomp.mod"],
method='natural'),
0.05, 0.95, subdivisions=subdivisions)$value -
integrate(splinefun(flxDf.mod[,"qcomp.obs.fdc"],
flxDf.mod[,"qcomp.obs"], method='natural'),
0.05, 0.95, subdivisions=subdivisions)$value, 2 )
results["daily", "n"] <- length(flxDf.mod.d$qcomp.mod)
results["daily", "nse"] <- round(Nse(flxDf.mod.d$qcomp.mod, flxDf.mod.d$qcomp.obs), 2)
results["daily", "nselog"] <- round(NseLog(flxDf.mod.d$qcomp.mod, flxDf.mod.d$qcomp.obs), 2)
results["daily", "cor"] <- round(cor(flxDf.mod.d$qcomp.mod, flxDf.mod.d$qcomp.obs,
use="na.or.complete"), 2)
results["daily", "rmse"] <- round(Rmse(flxDf.mod.d$qcomp.mod, flxDf.mod.d$qcomp.obs), 4)
results["daily", "rmsenorm"] <- round(RmseNorm(flxDf.mod.d$qcomp.mod,
flxDf.mod.d$qcomp.obs), 2)
results["daily", "bias"] <- round(sum(flxDf.mod.d$qcomp.mod-flxDf.mod.d$qcomp.obs, na.rm=T)/
sum(flxDf.mod.d$qcomp.obs, na.rm=TRUE) * 100, 1)
results["daily", "msd"] <- round(mean(flxDf.mod.d$qcomp.mod-flxDf.mod.d$qcomp.obs, na.rm=T), 4)
results["daily", "mae"] <- round(mean(abs(flxDf.mod.d$qcomp.mod-flxDf.mod.d$qcomp.obs),
na.rm=T), 4)
results["daily", "errcom"] <- round(mean(as.numeric(difftime(flxDf.mod.dcom$qcomp.mod,
flxDf.mod.dcom$qcomp.obs,
units="hours")), na.rm=T), 2)
results["daily", "errmaxt"] <- round(mean(as.numeric(difftime(flxDf.mod.dmax$qcomp.mod,
flxDf.mod.dmax$qcomp.obs,
units="hours")), na.rm=T), 2)
results["daily", "errfdc"] <- round(integrate(splinefun(flxDf.mod.d[,"qcomp.mod.fdc"],
flxDf.mod.d[,"qcomp.mod"],
method='natural'), 0.05, 0.95,
subdivisions=subdivisions)$value -
integrate(splinefun(flxDf.mod.d[,"qcomp.obs.fdc"],
flxDf.mod.d[,"qcomp.obs"],
method='natural'), 0.05, 0.95,
subdivisions=subdivisions)$value, 2 )
results["monthly", "n"] <- length(flxDf.mod.mwy$qcomp.mod)
results["monthly", "nse"] <- round(Nse(flxDf.mod.mwy$qcomp.mod, flxDf.mod.mwy$qcomp.obs), 2)
results["monthly", "nselog"] <- round(NseLog(flxDf.mod.mwy$qcomp.mod,
flxDf.mod.mwy$qcomp.obs), 2)
results["monthly", "cor"] <- round(cor(flxDf.mod.mwy$qcomp.mod, flxDf.mod.mwy$qcomp.obs,
use="na.or.complete"), 2)
results["monthly", "rmse"] <- round(Rmse(flxDf.mod.mwy$qcomp.mod, flxDf.mod.mwy$qcomp.obs), 4)
results["monthly", "rmsenorm"] <- round(RmseNorm(flxDf.mod.mwy$qcomp.mod,
flxDf.mod.mwy$qcomp.obs), 2)
results["monthly", "bias"] <- round(sum(flxDf.mod.mwy$qcomp.mod-flxDf.mod.mwy$qcomp.obs,
na.rm=T)/
sum(flxDf.mod.mwy$qcomp.obs, na.rm=TRUE) * 100, 1)
results["monthly", "msd"] <- round(mean(flxDf.mod.mwy$qcomp.mod-flxDf.mod.mwy$qcomp.obs,
na.rm=T), 4)
results["monthly", "mae"] <- round(mean(abs(flxDf.mod.mwy$qcomp.mod-flxDf.mod.mwy$qcomp.obs),
na.rm=T), 4)
results["monthly", "errcom"] <- round(mean(as.numeric(difftime(flxDf.mod.mwycom$qcomp.mod,
flxDf.mod.mwycom$qcomp.obs,
units="days")), na.rm=T), 2)
results["monthly", "errmaxt"] <- round(mean(as.numeric(difftime(flxDf.mod.mwymax$qcomp.mod,
flxDf.mod.mwymax$qcomp.obs,
units="days")), na.rm=T), 2)
results["monthly", "errfdc"] <- NA
results["yearly", "n"] <- length(flxDf.mod.wy$qcomp.mod)
results["yearly", "nse"] <- NA
#round(Nse(flxDf.mod.wy$qcomp.mod, flxDf.mod.wy$qcomp.obs), 2)
results["yearly", "nselog"] <- NA
#round(NseLog(flxDf.mod.wy$qcomp.mod, flxDf.mod.wy$qcomp.obs), 2)
results["yearly", "cor"] <- NA
#round(cor(flxDf.mod.wy$qcomp.mod, flxDf.mod.wy$qcomp.obs, use="na.or.complete"), 2)
results["yearly", "rmse"] <- NA
#round(Rmse(flxDf.mod.wy$qcomp.mod, flxDf.mod.wy$qcomp.obs), 2)
results["yearly", "rmsenorm"] <- NA
#round(RmseNorm(flxDf.mod.wy$qcomp.mod, flxDf.mod.wy$qcomp.obs), 2)
results["yearly", "bias"] <- round(sum(flxDf.mod.wy$qcomp.mod-flxDf.mod.wy$qcomp.obs, na.rm=T)/
sum(flxDf.mod.wy$qcomp.obs, na.rm=TRUE) * 100, 1)
results["yearly", "msd"] <- round(mean(flxDf.mod.wy$qcomp.mod-flxDf.mod.wy$qcomp.obs,
na.rm=T), 4)
results["yearly", "mae"] <- round(mean(abs(flxDf.mod.wy$qcomp.mod-flxDf.mod.wy$qcomp.obs),
na.rm=T), 4)
results["yearly", "errcom"] <- round(mean(as.numeric(difftime(flxDf.mod.wycom$qcomp.mod,
flxDf.mod.wycom$qcomp.obs,
units="days")), na.rm=T), 2)
results["yearly", "errmaxt"] <- round(mean(as.numeric(difftime(flxDf.mod.wymax$qcomp.mod,
flxDf.mod.wymax$qcomp.obs,
units="days")), na.rm=T), 2)
results["yearly", "errfdc"] <- NA
results["max10", "n"] <- length(flxDf.mod.max10$qcomp.mod)
results["max10", "nse"] <- round(Nse(flxDf.mod.max10$qcomp.mod, flxDf.mod.max10$qcomp.obs), 2)
results["max10", "nselog"] <- round(NseLog(flxDf.mod.max10$qcomp.mod,
flxDf.mod.max10$qcomp.obs), 2)
results["max10", "cor"] <- round(cor(flxDf.mod.max10$qcomp.mod, flxDf.mod.max10$qcomp.obs,
use="na.or.complete"), 2)
results["max10", "rmse"] <- round(Rmse(flxDf.mod.max10$qcomp.mod, flxDf.mod.max10$qcomp.obs), 4)
results["max10", "rmsenorm"] <- round(RmseNorm(flxDf.mod.max10$qcomp.mod,
flxDf.mod.max10$qcomp.obs), 2)
results["max10", "bias"] <- round(sum(flxDf.mod.max10$qcomp.mod-flxDf.mod.max10$qcomp.obs,
na.rm=T)/
sum(flxDf.mod.max10$qcomp.obs, na.rm=TRUE) * 100, 1)
results["max10", "msd"] <- round(mean(flxDf.mod.max10$qcomp.mod-flxDf.mod.max10$qcomp.obs,
na.rm=T), 4)
results["max10", "mae"] <- round(mean(abs(flxDf.mod.max10$qcomp.mod-flxDf.mod.max10$qcomp.obs),
na.rm=T), 4)
results["max10", "errcom"] <- NA
results["max10", "errmaxt"] <- NA
results["max10", "errfdc"] <- NA
results["min10", "n"] <- length(flxDf.mod.min10$qcomp.mod)
results["min10", "nse"] <- round(Nse(flxDf.mod.min10$qcomp.mod, flxDf.mod.min10$qcomp.obs), 2)
results["min10", "nselog"] <- round(NseLog(flxDf.mod.min10$qcomp.mod,
flxDf.mod.min10$qcomp.obs), 2)
results["min10", "cor"] <- round(cor(flxDf.mod.min10$qcomp.mod, flxDf.mod.min10$qcomp.obs,
use="na.or.complete"), 2)
results["min10", "rmse"] <- round(Rmse(flxDf.mod.min10$qcomp.mod, flxDf.mod.min10$qcomp.obs), 4)
results["min10", "rmsenorm"] <- round(RmseNorm(flxDf.mod.min10$qcomp.mod,
flxDf.mod.min10$qcomp.obs), 2)
results["min10", "bias"] <- round(sum(flxDf.mod.min10$qcomp.mod-flxDf.mod.min10$qcomp.obs,
na.rm=T)/
sum(flxDf.mod.min10$qcomp.obs, na.rm=TRUE) * 100, 1)
results["min10", "msd"] <- round(mean(flxDf.mod.min10$qcomp.mod-flxDf.mod.min10$qcomp.obs,
na.rm=T), 4)
results["min10", "mae"] <- round(mean(abs(flxDf.mod.min10$qcomp.mod-flxDf.mod.min10$qcomp.obs),
na.rm=T), 4)
results["min10", "errcom"] <- NA
results["min10", "errmaxt"] <- NA
results["min10", "errfdc"] <- NA
# Units
results["units",] <- c("count", "unitless", "unitless", "unitless",
"flux units", "unitless", "percent", "flux units",
"flux units", "hours|days", "hours|days",
"flux units")
results
}
#' Computes model performance statistics for WRF-Hydro flux output
#'
#' \code{CalcModPerfMulti} calculates model performance statistics for flux
#' output.
#'
#' \code{CalcModPerfMulti} reads a model flux time series (i.e., created using
#' \code{\link{ReadFrxstPts}}) and an observation time series (i.e., created
#' using \code{\link{ReadUsgsGage}}) and calculates model performance statistics
#' (Nash-Sutcliffe Efficiency, Rmse, etc.) at various time scales and for low
#' and high fluxes. The tool will subset data to matching time periods (e.g., if
#' the observed data is at 5-min increments and modelled data is at 1-hr
#' increments, the tool will subset the observed data to select only
#' observations on the matching hour break).
#'
#' \code{CalcModPerfMulti} calculates the same statistics as
#' \code{\link{CalcModPerf}}) but returns results as a single-row vs. multi-row
#' dataframe. This is intended to streamline compiling statistics for multiple
#' data runs or multiple sites.
#'
#' Performance Statistics: \cr (mod = model output, obs = observations, n =
#' sample size) \itemize{ \item n: sample size \item nse: Nash-Sutcliffe
#' Efficiency \deqn{nse = 1 - ( sum((obs - mod)^2) / sum((obs - mean(obs))^2) )
#' } \item nselog: log-transformed Nash-Sutcliffe Efficiency \deqn{nselog = 1 -
#' ( sum((log(obs) - log(mod))^2) / sum((log(obs) - mean(log(obs)))^2) ) } \item
#' cor: correlation coefficient \deqn{cor = cor(mod, obs) } \item rmse: root
#' mean squared error \deqn{rmse = sqrt( sum((mod - obs)^2) / n ) } \item
#' rmsenorm: normalized root mean squared error \deqn{rmsenorm = rmse /
#' (max(obs) - min(obs)) } \item bias: percent bias \deqn{bias = sum(mod - obs)
#' / sum(obs) * 100 } \item msd: mean signed deviation \deqn{msd = mean(mod -
#' obs) } \item mae: mean absolute error \deqn{mae = mean(abs(mod -
#' obs)) } \item errcom: error in the center-of-mass of the flux, where
#' center-of-mass is the hour/day when 50\% of daily/monthly/water-year flux has
#' occurred. Reported as number of hours for daily time scale and number of days
#' for monthly and yearly time scales. \item errmaxt: Error in the time of
#' maximum flux. Reported as number of hours for daily time scale and number of
#' days for monthly and yearly time scales). \item errfdc: Error in the
#' integrated flow duration curve between 0.05 and 0.95 exceedance thresholds
#' (in native flow units). \item errflash: Error in the Richards-Baker
#' Flashiness Index \deqn{R-B Index = sum(abs(q_t - q_t-1)) / sum(q)}. }
#'
#' Time scales/Flux types: \itemize{ \item t = native model/observation time
#' step (e.g., hourly) \item dy = daily time step \item mo = monthly time step
#' \item yr = water-year time step \item max10 = high flows; restricted to the
#' portion of the time series where the observed flux is in the highest 10\%
#' (native time step) \item min10 = low flows; restricted to the portion of the
#' time series where the observed flux is in the lowest 10\% (native time step)
#' }
#'
#' @param flxDf.mod The flux output dataframe (required). Assumes only one
#' forecast point per file, so if you have multiple forecast points in your
#' output dataframe, use subset to isolate a single forecast point's data.
#' Also assumes model output and observation both contain POSIXct fields
#' (called "POSIXct").
#' @param flxDf.obs The observed flux dataframe. Assumes only one observation
#' point per file, so if you have multiple observation points in your
#' dataframe, use subset to isolate a single point's data. Also assumes model
#' output and observation both contain POSIXct fields (called "POSIXct").
#' @param flxCol.mod The column name for the flux time series for the MODEL data
#' (default="q_cms")
#' @param flxCol.obs The column name for the flux time series for the OBSERVED
#' data (default="q_cms")
#' @param stdate Start date for statistics (DEFAULT=NULL, all records will be
#' used). Date MUST be specified in POSIXct format with appropriate timezone
#' (e.g., as.POSIXct("2013-05-01 00:00:00", format="\%Y-\%m-\%d \%H:\%M:\%S",
#' tz="UTC"))
#' @param enddate End date for statistics (DEFAULT=NULL, all records will be
#' used). Date MUST be specified in POSIXct format with appropriate timezone
#' (e.g., as.POSIXct("2013-05-01 00:00:00", format="\%Y-\%m-\%d \%H:\%M:\%S",
#' tz="UTC"))
#' @param writeOutPairDf character path/name for file to output the data frame constructed
#' before calculating the statistics and various aggregations. NULL by default gives
#' no output.
#' @param fileTag A to be inserted in the file name, just before the file extension.
#' @return A new dataframe containing the model performance statistics.
#'
#' @examples
#' ## Take forecast point model output for Fourmile Creek (modStrh.mod1.fc) and a
#' ## corresponding USGS gage observation file (obsStrh.fc), both at an hourly time
#' ## step, and calculate model performance statistics. The model forecast point
#' ## data was imported using ReadFrxstPts and the gage observation data was
#' ## imported using ReadUsgsGage.
#'
#'
#' \dontrun{
#' CalcModPerfMulti(modStrd.chrt.fc, obsStr5min.fc)
#'
#' > Output:
#' t_n t_nse t_nselog t_cor t_rmse t_rmsenorm t_bias t_msd t_mae t_errfdc t_errflash
#' 116 0.81 0.67 0.91 0.211 8.48 7.1 0.0313 0.1385 0.04 0.0132
#' dy_n dy_nse dy_nselog dy_cor dy_rmse dy_rmsenorm dy_bias dy_msd dy_mae dy_errcom dy_errmaxt dy_errfdc dy_errflash
#' 116 0.81 0.67 0.91 0.211 8.48 7.1 0.0313 0.1385 0 0 0.04 0.0092
#' mo_n mo_nse mo_nselog mo_cor mo_rmse mo_rmsenorm mo_bias mo_msd mo_mae mo_errcom mo_errmaxt yr_n
#' 5 0.92 0.82 0.96 0.1073 10.35 9.1 0.0324 0.0833 -1.2 0.2 1
#' yr_bias yr_msd yr_mae yr_errcom yr_errmaxt max10_n max10_nse max10_nselog max10_cor max10_rmse max10_rmsenorm max10_bias max10_msd max10_mae
#' 7.1 0.0313 0.0313 0 4 12 0.46 0.58 0.71 0.348 23.63 -2.8 -0.044 0.248
#' min10_n min10_nse min10_nselog min10_cor min10_rmse min10_rmsenorm min10_bias min10_msd min10_mae
#' 12 -385.99 -14.59 0.35 0.1678 538.78 112.4 0.056 0.0592
#' }
#' @keywords univar ts
#' @concept modelEval
#' @family modelEvaluation
#' @export
CalcModPerfMulti <- function (flxDf.mod, flxDf.obs, flxCol.mod="q_cms", flxCol.obs="q_cms",
stdate=NULL, enddate=NULL, writeOutPairDf=NULL, fileTag=NULL) {
# Internal functions
which.max.dt <- function(dd, qcol, dtcol) { dd[which.max(dd[,qcol]), dtcol] }
CalcCOM.dt <- function(dd, qcol, dtcol) { dd[CalcCOM(dd[,qcol]), dtcol] }
# Prepare data
if (!is.null(stdate) & is.null(enddate)) {
flxDf.obs <- subset(flxDf.obs, POSIXct>=stdate)
flxDf.mod <- subset(flxDf.mod, POSIXct>=stdate)
}
if (is.null(stdate) & !is.null(enddate)) {
flxDf.obs <- subset(flxDf.obs, POSIXct<=enddate)
flxDf.mod <- subset(flxDf.mod, POSIXct<=enddate)
}
if (!is.null(stdate) & !is.null(enddate)) {
flxDf.obs <- subset(flxDf.obs, POSIXct>=stdate & POSIXct<=enddate)
flxDf.mod <- subset(flxDf.mod, POSIXct>=stdate & POSIXct<=enddate)
}
flxDf.mod$qcomp <- flxDf.mod[,flxCol.mod]
flxDf.obs$qcomp <- flxDf.obs[,flxCol.obs]
# model timestep in secs
modT <- as.integer(flxDf.mod$POSIXct[2])-as.integer(flxDf.mod$POSIXct[1])
#if (as.integer(flxDf.obs$POSIXct[2])-as.integer(flxDf.obs$POSIXct[1]) >= 86400) {
# flxDf.obs$POSIXct=as.POSIXct(round(flxDf.obs$POSIXct,"days"), tz="UTC")}
flxDf.mod <- merge(flxDf.mod[c("POSIXct","qcomp")], flxDf.obs[c("POSIXct","qcomp")],
by<-c("POSIXct"), suffixes=c(".mod",".obs"))
flxDf.mod <- subset(flxDf.mod, !is.na(flxDf.mod$qcomp.mod) & !is.na(flxDf.mod$qcomp.obs))
if (nrow(flxDf.mod) == 0) {
warning("No data matches.")
return(NA)
}
flxDf.mod <- CalcDates(flxDf.mod)
flxDf.mod$date <- as.POSIXct(trunc(flxDf.mod$POSIXct, "days"))
if(!is.null(writeOutPairDf) && is.character(writeOutPairDf) && writeOutPairDf!='') {
modPerfPairDf <- flxDf.mod
wopdFileName <- if(!is.null(fileTag) && fileTag!='') {
wopdBase <- tools::file_path_sans_ext(writeOutPairDf)
wopdExt <- tools::file_ext(writeOutPairDf)
paste0(wopdBase,'.',fileTag,'.',wopdExt)
} else writeOutPairDf
print(wopdFileName)
save(modPerfPairDf, file=wopdFileName)
rm('modPerfPairDf')
}
results <- as.data.frame(matrix(nrow = 1, ncol = 59))
colnames(results) = c("t_n", "t_nse", "t_nselog", "t_cor", "t_rmse", "t_rmsenorm",
"t_bias", "t_msd", "t_mae", "t_errfdc", "t_errflash",
"dy_n", "dy_nse", "dy_nselog", "dy_cor", "dy_rmse",
"dy_rmsenorm", "dy_bias", "dy_msd", "dy_mae",
"dy_errcom", "dy_errmaxt", "dy_errfdc", "dy_errflash",
"mo_n", "mo_nse", "mo_nselog", "mo_cor", "mo_rmse",
"mo_rmsenorm", "mo_bias", "mo_msd", "mo_mae",
"mo_errcom", "mo_errmaxt",
"yr_n", "yr_bias", "yr_msd", "yr_mae", "yr_errcom", "yr_errmaxt",
"max10_n", "max10_nse", "max10_nselog", "max10_cor",
"max10_rmse", "max10_rmsenorm", "max10_bias",
"max10_msd", "max10_mae",
"min10_n", "min10_nse", "min10_nselog", "min10_cor",
"min10_rmse", "min10_rmsenorm", "min10_bias",
"min10_msd", "min10_mae")
exclvars <- names(flxDf.mod) %in% c("POSIXct", "secs", "timest", "date", "stat")
# Base aggregations
flxDf.mod.d <- aggregate(flxDf.mod[!exclvars], by = list(flxDf.mod$date), CalcMeanNarm)
flxDf.mod.d$mwy <- paste0(flxDf.mod.d$month, "-", flxDf.mod.d$wy)
flxDf.mod.mwy <- aggregate(flxDf.mod[c("qcomp.mod","qcomp.obs")],
by = list(flxDf.mod$month, flxDf.mod$wy), CalcMeanNarm)
flxDf.mod.wy <- aggregate(flxDf.mod[c("qcomp.mod","qcomp.obs")],
by = list(flxDf.mod$wy), CalcMeanNarm)
# Time of center of mass aggregations
tmp.mod <- plyr::ddply(flxDf.mod, plyr::.(date), CalcCOM.dt, qcol="qcomp.mod",
dtcol="POSIXct")
tmp.obs <- plyr::ddply(flxDf.mod, plyr::.(date), CalcCOM.dt, qcol="qcomp.obs",
dtcol="POSIXct")
flxDf.mod.dcom <- merge(tmp.mod, tmp.obs, by=c("date"))
colnames(flxDf.mod.dcom) <- c("date", "qcomp.mod", "qcomp.obs")
tmp.mod <- plyr::ddply(flxDf.mod.d, plyr::.(month, wy), CalcCOM.dt, qcol="qcomp.mod",
dtcol="Group.1")
tmp.obs <- plyr::ddply(flxDf.mod.d, plyr::.(month, wy), CalcCOM.dt, qcol="qcomp.obs",
dtcol="Group.1")
flxDf.mod.mwycom <- merge(tmp.mod, tmp.obs, by=c("month", "wy"))
colnames(flxDf.mod.mwycom) <- c("month", "wy", "qcomp.mod", "qcomp.obs")
tmp.mod <- plyr::ddply(flxDf.mod.d, plyr::.(wy), CalcCOM.dt, qcol="qcomp.mod",
dtcol="Group.1")
tmp.obs <- plyr::ddply(flxDf.mod.d, plyr::.(wy), CalcCOM.dt, qcol="qcomp.obs",
dtcol="Group.1")
flxDf.mod.wycom <- merge(tmp.mod, tmp.obs, by=c("wy"))
colnames(flxDf.mod.wycom) <- c("wy", "qcomp.mod", "qcomp.obs")
# Time of max aggregations
tmp.mod <- plyr::ddply(flxDf.mod, plyr::.(date), which.max.dt, qcol="qcomp.mod",
dtcol="POSIXct")
tmp.obs <- plyr::ddply(flxDf.mod, plyr::.(date), which.max.dt, qcol="qcomp.obs",
dtcol="POSIXct")
flxDf.mod.dmax <- merge(tmp.mod, tmp.obs, by=c("date"))
colnames(flxDf.mod.dmax) <- c("date", "qcomp.mod", "qcomp.obs")
tmp.mod <- plyr::ddply(flxDf.mod.d, plyr::.(month, wy), which.max.dt, qcol="qcomp.mod",
dtcol="Group.1")
tmp.obs <- plyr::ddply(flxDf.mod.d, plyr::.(month, wy), which.max.dt, qcol="qcomp.obs",
dtcol="Group.1")
flxDf.mod.mwymax <- merge(tmp.mod, tmp.obs, by=c("month", "wy"))
colnames(flxDf.mod.mwymax) <- c("month", "wy", "qcomp.mod", "qcomp.obs")
tmp.mod <- plyr::ddply(flxDf.mod.d, plyr::.(wy), which.max.dt, qcol="qcomp.mod",
dtcol="Group.1")
tmp.obs <- plyr::ddply(flxDf.mod.d, plyr::.(wy), which.max.dt, qcol="qcomp.obs",
dtcol="Group.1")
flxDf.mod.wymax <- merge(tmp.mod, tmp.obs, by=c("wy"))
colnames(flxDf.mod.wymax) <- c("wy", "qcomp.mod", "qcomp.obs")
# NAs cause which.max to return a list, so force back to ints
# (out for now since we pre-screen for NAs)
#flxDf.mod.dmax$qcomp.mod <- as.integer(flxDf.mod.dmax$qcomp.mod)
#flxDf.mod.dmax$qcomp.obs <- as.integer(flxDf.mod.dmax$qcomp.obs)
#flxDf.mod.mwymax$qcomp.mod <- as.integer(flxDf.mod.mwymax$qcomp.mod)
#flxDf.mod.mwymax$qcomp.obs <- as.integer(flxDf.mod.mwymax$qcomp.obs)
#flxDf.mod.wymax$qcomp.mod <- as.integer(flxDf.mod.wymax$qcomp.mod)
#flxDf.mod.wymax$qcomp.obs <- as.integer(flxDf.mod.wymax$qcomp.obs)
# Mins and Maxes
flxDf.mod.max10 <- subset(flxDf.mod, flxDf.mod$qcomp.obs>=quantile(flxDf.mod$qcomp.obs,
0.90, na.rm=TRUE))
flxDf.mod.min10 <- subset(flxDf.mod, flxDf.mod$qcomp.obs<=quantile(flxDf.mod$qcomp.obs,
0.10, na.rm=TRUE))
# FDCs
flxDf.mod <- CalcFdc(flxDf.mod, "qcomp.mod")
flxDf.mod <- CalcFdc(flxDf.mod, "qcomp.obs")
flxDf.mod.d <- CalcFdc(flxDf.mod.d, "qcomp.mod")
flxDf.mod.d <- CalcFdc(flxDf.mod.d, "qcomp.obs")
# Compile summary statistics
results["t_n"] <- length(flxDf.mod$qcomp.mod)
results["t_nse"] <- round(Nse(flxDf.mod$qcomp.mod, flxDf.mod$qcomp.obs), 2)
results["t_nselog"] <- round(NseLog(flxDf.mod$qcomp.mod, flxDf.mod$qcomp.obs), 2)
results["t_cor"] <- round(cor(flxDf.mod$qcomp.mod, flxDf.mod$qcomp.obs,
use="na.or.complete"), 2)
results["t_rmse"] <- round(Rmse(flxDf.mod$qcomp.mod, flxDf.mod$qcomp.obs), 4)
results["t_rmsenorm"] <- round(RmseNorm(flxDf.mod$qcomp.mod, flxDf.mod$qcomp.obs), 2)
results["t_bias"] <- round(sum(flxDf.mod$qcomp.mod-flxDf.mod$qcomp.obs, na.rm=T)/
sum(flxDf.mod$qcomp.obs, na.rm=TRUE) * 100, 1)
results["t_msd"] <- round(mean(flxDf.mod$qcomp.mod-flxDf.mod$qcomp.obs, na.rm=T), 4)
results["t_mae"] <- round(mean(abs(flxDf.mod$qcomp.mod-flxDf.mod$qcomp.obs), na.rm=T), 4)
results["t_errfdc"] <- round(integrate(splinefun(flxDf.mod[,"qcomp.mod.fdc"],
flxDf.mod[,"qcomp.mod"], method='natural'),
0.05, 0.95, subdivisions=2000)$value -
integrate(splinefun(flxDf.mod[,"qcomp.obs.fdc"],
flxDf.mod[,"qcomp.obs"], method='natural'),
0.05, 0.95, subdivisions=2000)$value, 2 )
results["t_errflash"] <- round(RBFlash(flxDf.mod$qcomp.obs) - RBFlash(flxDf.mod$qcomp.mod), 2)
results["dy_n"] <- length(flxDf.mod.d$qcomp.mod)
results["dy_nse"] <- round(Nse(flxDf.mod.d$qcomp.mod, flxDf.mod.d$qcomp.obs), 2)
results["dy_nselog"] <- round(NseLog(flxDf.mod.d$qcomp.mod, flxDf.mod.d$qcomp.obs), 2)
results["dy_cor"] <- round(cor(flxDf.mod.d$qcomp.mod, flxDf.mod.d$qcomp.obs,
use="na.or.complete"), 2)
results["dy_rmse"] <- round(Rmse(flxDf.mod.d$qcomp.mod, flxDf.mod.d$qcomp.obs), 4)
results["dy_rmsenorm"] <- round(RmseNorm(flxDf.mod.d$qcomp.mod, flxDf.mod.d$qcomp.obs), 2)
results["dy_bias"] <- round(sum(flxDf.mod.d$qcomp.mod-flxDf.mod.d$qcomp.obs, na.rm=T)/
sum(flxDf.mod.d$qcomp.obs, na.rm=TRUE) * 100, 1)
results["dy_msd"] <- round(mean(flxDf.mod.d$qcomp.mod-flxDf.mod.d$qcomp.obs, na.rm=T), 4)
results["dy_mae"] <- round(mean(abs(flxDf.mod.d$qcomp.mod-flxDf.mod.d$qcomp.obs), na.rm=T), 4)
results["dy_errcom"] <- round(mean(as.numeric(difftime(flxDf.mod.dcom$qcomp.mod,
flxDf.mod.dcom$qcomp.obs,
units="hours")), na.rm=T), 2)
results["dy_errmaxt"] <- round(mean(as.numeric(difftime(flxDf.mod.dmax$qcomp.mod,
flxDf.mod.dmax$qcomp.obs,
units="hours")), na.rm=T), 2)
results["dy_errfdc"] <- round(integrate(splinefun(flxDf.mod.d[,"qcomp.mod.fdc"],
flxDf.mod.d[,"qcomp.mod"],
method='natural'), 0.05, 0.95,
subdivisions=2000)$value -
integrate(splinefun(flxDf.mod.d[,"qcomp.obs.fdc"],
flxDf.mod.d[,"qcomp.obs"],
method='natural'), 0.05, 0.95,
subdivisions=2000)$value, 2 )
results["dy_errflash"] <- round(RBFlash(flxDf.mod.d$qcomp.obs) - RBFlash(flxDf.mod.d$qcomp.mod), 2)
results["mo_n"] <- length(flxDf.mod.mwy$qcomp.mod)
results["mo_nse"] <- round(Nse(flxDf.mod.mwy$qcomp.mod, flxDf.mod.mwy$qcomp.obs), 2)
results["mo_nselog"] <- round(NseLog(flxDf.mod.mwy$qcomp.mod, flxDf.mod.mwy$qcomp.obs), 2)
results["mo_cor"] <- round(cor(flxDf.mod.mwy$qcomp.mod, flxDf.mod.mwy$qcomp.obs,
use="na.or.complete"), 2)
results["mo_rmse"] <- round(Rmse(flxDf.mod.mwy$qcomp.mod, flxDf.mod.mwy$qcomp.obs), 4)
results["mo_rmsenorm"] <- round(RmseNorm(flxDf.mod.mwy$qcomp.mod, flxDf.mod.mwy$qcomp.obs), 2)
results["mo_bias"] <- round(sum(flxDf.mod.mwy$qcomp.mod-flxDf.mod.mwy$qcomp.obs, na.rm=T)/
sum(flxDf.mod.mwy$qcomp.obs, na.rm=TRUE) * 100, 1)
results["mo_msd"] <- round(mean(flxDf.mod.mwy$qcomp.mod-flxDf.mod.mwy$qcomp.obs, na.rm=T), 4)
results["mo_mae"] <- round(mean(abs(flxDf.mod.mwy$qcomp.mod-flxDf.mod.mwy$qcomp.obs), na.rm=T), 4)
results["mo_errcom"] <- round(mean(as.numeric(difftime(flxDf.mod.mwycom$qcomp.mod,
flxDf.mod.mwycom$qcomp.obs,
units="days")), na.rm=T), 2)
results["mo_errmaxt"] <- round(mean(as.numeric(difftime(flxDf.mod.mwymax$qcomp.mod,
flxDf.mod.mwymax$qcomp.obs,
units="days")), na.rm=T), 2)
results["yr_n"] <- length(flxDf.mod.wy$qcomp.mod)
results["yr_bias"] <- round(sum(flxDf.mod.wy$qcomp.mod-flxDf.mod.wy$qcomp.obs,
na.rm=T)/sum(flxDf.mod.wy$qcomp.obs, na.rm=TRUE) * 100, 1)
results["yr_msd"] <- round(mean(flxDf.mod.wy$qcomp.mod-flxDf.mod.wy$qcomp.obs, na.rm=T), 4)
results["yr_mae"] <- round(mean(abs(flxDf.mod.wy$qcomp.mod-flxDf.mod.wy$qcomp.obs), na.rm=T), 4)
results["yr_errcom"] <- round(mean(as.numeric(difftime(flxDf.mod.wycom$qcomp.mod,
flxDf.mod.wycom$qcomp.obs,
units="days")), na.rm=T), 2)
results["yr_errmaxt"] <- round(mean(as.numeric(difftime(flxDf.mod.wymax$qcomp.mod,
flxDf.mod.wymax$qcomp.obs,
units="days")), na.rm=T), 2)
results["max10_n"] <- length(flxDf.mod.max10$qcomp.mod)
results["max10_nse"] <- round(Nse(flxDf.mod.max10$qcomp.mod, flxDf.mod.max10$qcomp.obs), 2)
results["max10_nselog"] <- round(NseLog(flxDf.mod.max10$qcomp.mod, flxDf.mod.max10$qcomp.obs), 2)
results["max10_cor"] <- round(cor(flxDf.mod.max10$qcomp.mod, flxDf.mod.max10$qcomp.obs,
use="na.or.complete"), 2)
results["max10_rmse"] <- round(Rmse(flxDf.mod.max10$qcomp.mod, flxDf.mod.max10$qcomp.obs), 4)
results["max10_rmsenorm"] <- round(RmseNorm(flxDf.mod.max10$qcomp.mod,
flxDf.mod.max10$qcomp.obs), 2)
results["max10_bias"] <- round(sum(flxDf.mod.max10$qcomp.mod-flxDf.mod.max10$qcomp.obs,
na.rm=T)/sum(flxDf.mod.max10$qcomp.obs, na.rm=TRUE) * 100, 1)
results["max10_msd"] <- round(mean(flxDf.mod.max10$qcomp.mod-flxDf.mod.max10$qcomp.obs,
na.rm=T), 4)
results["max10_mae"] <- round(mean(abs(flxDf.mod.max10$qcomp.mod-flxDf.mod.max10$qcomp.obs),
na.rm=T), 4)
results["min10_n"] <- length(flxDf.mod.min10$qcomp.mod)
results["min10_nse"] <- round(Nse(flxDf.mod.min10$qcomp.mod, flxDf.mod.min10$qcomp.obs), 2)
results["min10_nselog"] <- round(NseLog(flxDf.mod.min10$qcomp.mod, flxDf.mod.min10$qcomp.obs), 2)
results["min10_cor"] <- round(cor(flxDf.mod.min10$qcomp.mod, flxDf.mod.min10$qcomp.obs,
use="na.or.complete"), 2)
results["min10_rmse"] <- round(Rmse(flxDf.mod.min10$qcomp.mod, flxDf.mod.min10$qcomp.obs), 4)
results["min10_rmsenorm"] <- round(RmseNorm(flxDf.mod.min10$qcomp.mod,
flxDf.mod.min10$qcomp.obs), 2)
results["min10_bias"] <- round(sum(flxDf.mod.min10$qcomp.mod-flxDf.mod.min10$qcomp.obs,
na.rm=T)/sum(flxDf.mod.min10$qcomp.obs, na.rm=TRUE) * 100, 1)
results["min10_msd"] <- round(mean(flxDf.mod.min10$qcomp.mod-flxDf.mod.min10$qcomp.obs,
na.rm=T), 4)
results["min10_mae"] <- round(mean(abs(flxDf.mod.min10$qcomp.mod-flxDf.mod.min10$qcomp.obs),
na.rm=T), 4)
# Attributes
attr(results, "descrip") <- c(t_n="sample size",
t_nse="nash-sutcliffe efficiency",
t_nselog="log nash-sutcliffe efficiency",
t_cor="correlation",
t_rmse="root mean squared error",
t_rmsenorm="normalized root mean squared error",
t_bias="bias",
t_msd="mean signed deviation",
t_mae="mean absolute error",
t_errfdc="integrated flow duration curve error",
t_errflash="flashiness index error",
dy_n="sample size",
dy_nse="nash-sutcliffe efficiency",
dy_nselog="log nash-sutcliffe efficiency",
dy_cor="correlation",
dy_rmse="root mean squared error",
dy_rmsenorm="normalized root mean squared error",
dy_bias="bias",
dy_msd="mean signed deviation",
dy_mae="mean absolute error",
dy_errcom="error in center-of-mass timing",
dy_errmaxt="error in max timing",
dy_errfdc="integrated flow duration curve error",
dy_errflash="flashiness index error",
mo_n="sample size",
mo_nse="nash-sutcliffe efficiency",
mo_nselog="log nash-sutcliffe efficiency",
mo_cor="correlation",
mo_rmse="root mean squared error",
mo_rmsenorm="normalized root mean squared error",
mo_bias="bias",
mo_msd="mean signed deviation",
mo_mae="mean absolute error",
mo_errcom="error in center-of-mass timing",
mo_errmaxt="error in max timing",
yr_n="sample size",
yr_bias="bias",
yr_msd="mean signed deviation",
yr_mae="mean absolute error",
yr_errcom="error in center-of-mass timing",
yr_errmaxt="error in max timing",
max10_n="sample size",
max10_nse="nash-sutcliffe efficiency",
max10_nselog="log nash-sutcliffe efficiency",
max10_cor="correlation",
max10_rmse="root mean squared error",
max10_rmsenorm="normalized root mean squared error",
max10_bias="bias",
max10_msd="mean signed deviation",
max10_mae="mean absolute error",
min10_n="sample size",
min10_nse="nash-sutcliffe efficiency",
min10_nselog="log nash-sutcliffe efficiency",
min10_cor="correlation",
min10_rmse="root mean squared error",
min10_rmsenorm="normalized root mean squared error",
min10_bias="bias",
min10_msd="mean signed deviation",
min10_mae="mean absolute error")
attr(results, "units") <- c(t_n="count",
t_nse="unitless (-Inf->1)",
t_nselog="unitless (-Inf->1)",
t_cor="unitless (-1->1)",
t_rmse="native flux units",
t_rmsenorm="unitless (0->Inf)",
t_bias="percent (%)",
t_msd="native flux units",
t_mae="native flux units",
t_errfdc="native flux units",
t_errflash="unitless",
dy_n="count",
dy_nse="unitless (-Inf->1)",
dy_nselog="unitless (-Inf->1)",
dy_cor="unitless (-1->1)",
dy_rmse="native flux units",
dy_rmsenorm="native flux units",
dy_bias="percent (%)",
dy_msd="native flux units",
dy_mae="native flux units",
dy_errcom="hours",
dy_errmaxt="hours",
dy_errfdc="native flux units",
dy_errflash="unitless",
mo_n="count",
mo_nse="unitless (-Inf->1)",
mo_nselog="unitless (-Inf->1)",
mo_cor="unitless (-1->1)",
mo_rmse="native flux units",
mo_rmsenorm="native flux units",
mo_bias="percent (%)",
mo_msd="native flux units",
mo_mae="native flux units",
mo_errcom="days",
mo_errmaxt="days",
yr_n="count",
yr_bias="percent (%)",
yr_msd="native flux units",
yr_mae="native flux units",
yr_errcom="days",
yr_errmaxt="days",
max10_n="sample size",
max10_nse="unitless (-Inf->1)",
max10_nselog="unitless (-Inf->1)",
max10_cor="unitless (-1->1)",
max10_rmse="native flux units",
max10_rmsenorm="native flux units",
max10_bias="percent (%)",
max10_msd="native flux units",
max10_mae="native flux units",
min10_n="sample size",
min10_nse="unitless (-Inf->1)",
min10_nselog="unitless (-Inf->1)",
min10_cor="unitless (-1->1)",
min10_rmse="native flux units",
min10_rmsenorm="native flux units",
min10_bias="percent (%)",
min10_msd="native flux units",
min10_mae="native flux units")
results
}
#' Computes flow duration curve statistics for WRF-Hydro streamflow output
#'
#' \code{CalcFdcPerf} calculates flow duration curve statistics for streamflow
#' output.
#'
#' \code{CalcFdcPerf} reads a model forecast point streamflow timeseries (i.e.,
#' created using \code{\link{ReadFrxstPts}}) and a streamflow observation
#' timeseries (i.e., created using \code{\link{ReadUsgsGage}}) and calculates
#' flow duration curve statistics at various exceedance thresholds (e.g., 10\%,
#' 20\%, etc.). The tool will subset data to matching time periods (e.g., if the
#' observed data is at 5-min increments and modelled data is at 1-hr increments,
#' the tool will subset the observed data to select only observations on the
#' matching hour break).
#'
#' Flow Duration Curve Statistics: \cr (mod = model output, obs = observations)
#' \itemize{ \item p.exceed: exceedance threshold (e.g., 0.2 means a flow value
#' that is exceeded 20\% of the time) \item q.mod: MODEL flow value at specified
#' exceedance threshold (in native flow units) \item q.obs: OBSERVED flow value
#' at specified exceedance threshold (in native flow units) \item q.err:
#' difference between model and observed flow values [mod-obs] (in native flow
#' units) \item q.perr: percent error in model flow [(mod-obs)/obs] }
#'
#' @param strDf.mod The forecast point output dataframe (required). Assumes only
#' one forecast point per file, so if you have multiple forecast points in
#' your output dataframe, use subset to isolate a single forecast point's
#' data. Also assumes model output and observation both contain POSIXct fields
#' (called "POSIXct").
#' @param strDf.obs The observed streamflow dataframe. Assumes only one gage per
#' file, so if you have multiple gages in your dataframe, use subset to
#' isolate a single gage's data. Also assumes model output and observation
#' both contain POSIXct fields (called "POSIXct").
#' @param strCol.mod The column name for the streamflow time series for the
#' MODEL data (default="q_cms")
#' @param strCol.obs The column name for the streamflow time series for the
#' OBSERVED data (default="q_cms")
#' @param stdate Start date for statistics (DEFAULT=NULL, all records will be
#' used). Date MUST be specified in POSIXct format with appropriate timezone
#' (e.g., as.POSIXct("2013-05-01 00:00:00", format="\%Y-\%m-\%d \%H:\%M:\%S",
#' tz="UTC"))
#' @param enddate End date for statistics (DEFAULT=NULL, all records will be
#' used). Date MUST be specified in POSIXct format with appropriate timezone
#' (e.g., as.POSIXct("2013-05-01 00:00:00", format="\%Y-\%m-\%d \%H:\%M:\%S",
#' tz="UTC"))
#' @return A new dataframe containing the flow duration curve statistics.
#'
#' @examples
#' ## Take forecast point model output for Fourmile Creek (modStrh.mod1.fc)
#' ## and a corresponding USGS gage observation file (obsStrh.fc), both at an
#' ## hourly time step, and calculate flow duration curve statistics. The
#' ## model forecast point data was imported using ReadFrxstPts and the gage
#' ## observation data was imported using ReadUsgsGage.
#'
#' \dontrun{
#' CalcFdcPerf(modStr1h.allrt.fc, obsStr5min.fc)
#'
#' Output:
#' p.exceed q.mod q.obs
#' 0.1 3.07 5.25
#' 0.2 1.35 2.31
#' 0.3 0.82 1.06
#' 0.4 0.48 0.65
#' 0.5 0.29 0.45
#' 0.6 0.18 0.34
#' 0.7 0.14 0.25
#' 0.8 0.11 0.19
#' 0.9 0.08 0.16
#' }
#' @keywords univar ts
#' @concept modelEval
#' @family modelEvaluation flowDurationCurves
#' @export
CalcFdcPerf <- function (strDf.mod, strDf.obs, strCol.mod="q_cms", strCol.obs="q_cms",
stdate=NULL, enddate=NULL) {
# Prepare data
if (!is.null(stdate) && !is.null(enddate)) {
strDf.obs <- subset(strDf.obs, POSIXct>=stdate & POSIXct<=enddate)
strDf.mod <- subset(strDf.mod, POSIXct>=stdate & POSIXct<=enddate)
}
strDf.mod <- CalcDates(strDf.mod)
strDf.mod$qcomp <- strDf.mod[,strCol.mod]
strDf.obs$qcomp <- strDf.obs[,strCol.obs]
if (as.integer(strDf.obs$POSIXct[2])-as.integer(strDf.obs$POSIXct[1]) >= 86400) {
strDf.obs$POSIXct=as.POSIXct(round(strDf.obs$POSIXct,"days"), tz="UTC")}
strDf.mod <- merge(strDf.mod, strDf.obs[c("POSIXct","qcomp")], by<-c("POSIXct"),
suffixes=c(".mod",".obs"))
#strDf.mod <- subset(strDf.mod, !is.na(strDf$qcomp.mod) & !is.na(strDf$qcomp.obs))
strDf.mod <- CalcFdc(strDf.mod, "qcomp.mod")
strDf.mod <- CalcFdc(strDf.mod, "qcomp.obs")
fdc.mod <- splinefun(strDf.mod[,"qcomp.mod.fdc"], strDf.mod[,"qcomp.mod"],
method='natural')
fdc.obs <- splinefun(strDf.mod[,"qcomp.obs.fdc"], strDf.mod[,"qcomp.obs"],
method='natural')
# Compile summary statistics
results <- as.data.frame(matrix(nrow = 9, ncol = 5))
colnames(results) <- c("p.exceed","q.mod","q.obs","q.err","q.perr")
results[, 1] <- seq(0.1, 0.9, 0.1)
for (i in 1:9) {
results[i, "q.mod"] <- round(fdc.mod(results[i,1]), 2)
results[i, "q.obs"] <- round(fdc.obs(results[i,1]), 2)
results[i, "q.err"] <- results[i, "q.mod"] - results[i, "q.obs"]
results[i, "q.perr"] <- round(results[i, "q.err"] / results[i, "q.obs"] * 100, 1)
}
results
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.