#' Create sinusoidal time series of potential evapotranspiration input
#'
#' @description Generate series of potential evapotranspiration
#'
#' @param ts as vector of POSIXct data/times
#' @param eMin Minimum daily PE total (m or mm)
#' @param eMax Maximum daily PE total (m or mm)
#'
#' @details Dynamic TOPMODEL requires a time series of potential
#' evapotranspiration in order to calculate and remove actual
#' evapotranspiration from the root zone during a run. Many sophisticated
#' physical models have been developed for estimating potential and actual evapotranspiration, including the
#' Priestly-Taylor (Priestley and Taylor, 1972) and Penman-Monteith (Montieth,
#' 1965) methods. These, however, require detailed meteorological data such as
#' radiation input and relative humidities that are, in general, difficult to
#' obtain. Calder (1983) demonstrated that a simple approximation using a
#' sinusoidal variation in potential evapotranspiration to be a good
#' approximation to more complex schemes.
#'
#' If the insolation is also taken to vary sinusoidally through the daylight
#' hours then, ignoring diurnal meteorological variations, the potential
#' evapotranspiration during daylight hours for each year day number can be
#' calculated (for the catchment's latitude). Integration over the daylight
#' hours allows the daily maximum to be calculated and thus a sub-daily series
#' generated.
#'
#' @return Time series (xts) of potential evapotranspiration totals for the time steps given in same units as eMin and eMax
#'
#' @references Beven, K. J. (2012). Rainfall-runoff modelling : the primer. Chichester, UK, Wiley-Blackwell.
#' @references Calder, I. R. (1986). A stochastic model of rainfall interception. Journal of Hydrology, 89(1), 65-71.
#'
#' @examples
#' ## Generating daily PET data for 1970
#' ## the values of eMin and eMax may not by not be realistic
#' st <- as.POSIXct("1970-01-02 00:00:00",tz='GMT')
#' fn <- as.POSIXct("1971-01-01 00:00:00",tz='GMT')
#' daily_ts <- seq(st,fn,by=24*60*60)
#' dpet <- evap_est(daily_ts,0,1)
#'
#' ## create hourly data for the same period
#' st <- as.POSIXct("1970-01-01 01:00:00",tz='GMT')
#' fn <- as.POSIXct("1971-01-01 00:00:00",tz='GMT')
#' hour_ts <- seq(st,fn,by=1*60*60)
#' hpet <- evap_est(hour_ts,0,1)
#'
#' ## the totals should eb the same...
#' stopifnot(all.equal(sum(hpet), sum(dpet)))
#' @export
evap_est <- function(ts, eMin=0, eMax=0){
## Check min and max
if(!(eMin < eMax)){
stop("eMin should be less then eMax")
}
## Check timestep
dt <- diff(as.numeric(ts))
if(!all(dt[]==dt[1])){
stop("Irregularly spaced time series supplied")
}else{
dt <- dt[1]
}
if(dt > 24*60*60){ stop("Time step can be no mroe then a day") }
## create a series of daily PET values based on eMin and eMax
## day 0 is Jan 1, 31 Dec is day 364 or day 365 depending on if leap year
yday <- 0:365
fact <- 1+sin(2*pi*yday/365-pi/2)
daily_pet <- eMin + 0.5*(eMax-eMin)*fact
dawn <- (10 - 2.5*fact)*60*60 # in seconds from start of day
dayLength <- (6 + 4*fact) * 60*60 # in sec from start fo day
## make some elements for computation
sts <- ts-dt # this is the start of the timestep
dts <- as.POSIXct(3600*24*floor(as.numeric(ts)/(3600*24)),origin="1970-01-01",tz='GMT') # start of the day on which the timestep finishes
dsts <- as.POSIXct(3600*24*floor(as.numeric(sts)/(3600*24)),origin="1970-01-01",tz='GMT') # start of the day on which the timestep starts
##
pet <- rep(0,length(ts))
tmp <- dsts
idx <- tmp<dts
while( any(idx) ){
pet[idx] <- pet[idx] + daily_pet[ as.POSIXlt(tmp[idx])$yday + 1]
tmp <- tmp + 24*60*60
idx <- tmp<dts
}
## take away the fraction from the start
sc <- as.numeric(sts) - as.numeric(dsts)
frc <- (sc - dawn[ as.POSIXlt(sts)$yday +1 ])/ dayLength[ as.POSIXlt(sts)$yday +1]
frc <- pmin(1,pmax(0,frc))
pet <- pet - daily_pet[ as.POSIXlt(sts)$yday + 1]*0.5*(1-cos(frc*pi))
## add fraction from start fo currect day
sc <- as.numeric(ts) - as.numeric(dts)
frc <- (sc - dawn[ as.POSIXlt(ts)$yday +1 ])/ dayLength[ as.POSIXlt(ts)$yday +1]
frc <- pmin(1,pmax(0,frc))
pet <- pet + daily_pet[ as.POSIXlt(ts)$yday + 1]*0.5*(1-cos(frc*pi))
return(xts::xts(pet,order.by=ts))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.