R/water_ET.R

Defines functions dailyET ET24h hourlyET

Documented in dailyET ET24h hourlyET

#' Calculates ET using Penman Monteith hourly formula
#' @param WeatherStation     a data frame with all the needed fields (see 
#' example)
#' @param hours              time of the day in hours in 24hs format
#' @param DOY                day of year
#' @param long.z             longitude for local time
#' @param ET.instantaneous   Logical. True if you want to calculate 
#' instantaneous ET instead of hourly ET. See Details.
#' @param ET       "ETo" for short crops, similar to clipped, cool-season
#' grass; or "ETr" for tall crops, similar to 0.5 m tall full-cover alfalfa.
#' @param height             weather station sensors height in meters
#' @param lat                latitude of weather station in decimal degrees. 
#' Negative values for south latitude
#' @param long               longitude of weather station in decimal degrees. 
#' Negative values for west longitude
#' @param elev               elevation of weather station in meters
#' @details 
#' The only difference on instantaneous ET is how the hour is interpreted. On 
#' FALSE, and for example at 11:00, ET is calculated between 10:00 and 11:00, 
#' on TRUE Et is calculated at 11:00 hs.
#' @return ET hourly in mm.h-1
#' @author Guillermo Federico Olmedo
#' @examples 
#' WeatherStation  <- data.frame(wind=4.72,
#'                               RH=59, 
#'                               temp=24.3,
#'                               radiation=675, 
#'                               height=2.2, 
#'                               lat=-35.37, 
#'                               long=71.5946, 
#'                               elev=124)
#'   hourlyET(WeatherStation, hours=10.5, DOY=363, long.z=71.635)
#' @export
#' @references 
#' Allen 2005 ASCE
hourlyET <- function(WeatherStation, hours, DOY, long.z=WeatherStation$long, 
                     ET.instantaneous=FALSE, ET="ETr", height=2, lat, long, elev){
  if(class(WeatherStation)== "waterWeatherStation"){
    if(!is.null(WeatherStation$at.sat)){
      WeatherStation <- getDataWS(WeatherStation)
      ET.instantaneous = TRUE
      DOY <- WeatherStation$DOY
      hours <-  WeatherStation$hours}
  } else {
    if("datetime" %in% names(WeatherStation)){
      date <- as.POSIXlt(WeatherStation$datetime, format="%Y-%m-%d %H:%M:%S")
      DOY=date$yday+1
      hours=date$hour + date$min/60 + date$sec/3600}
    if(!missing(lat)){WeatherStation$lat <- lat}
    if(!missing(long)){WeatherStation$long <- long}
    if(!missing(elev)){WeatherStation$elev <- elev}
    if(!missing(height)){WeatherStation$height <- height}
    if(length(WeatherStation$DOY) != 0 & length(WeatherStation$hours) != 0 &
       missing(DOY) & missing(hours)){
      DOY <- WeatherStation$DOY
      hours <- WeatherStation$hours
    } 
  }
  tempK <- WeatherStation$temp + 273.16
  Rs <- WeatherStation$radiation * 3600 / 1e6
  P <- 101.3*((293-0.0065*WeatherStation$elev)/293)^5.26
  psi <- 0.000665*P
  Delta <- 2503 * exp((17.27*WeatherStation$temp)/
                  (WeatherStation$temp+237.3))/((WeatherStation$temp+237.3)^2)
  ea.sat <- 0.61078*exp((17.269*WeatherStation$temp)/
                          (WeatherStation$temp+237.3))
  ea <- (WeatherStation$RH/100)*ea.sat
  DPV <- ea.sat - ea
  dr <- 1 + 0.033*(cos(2*pi*DOY/365))
  delta <- 0.409*sin((2*pi*DOY/365)-1.39)
  phi <- WeatherStation$lat*(pi/180)
  b <- 2*pi*(DOY-81)/364
  Sc <- 0.1645*sin(2*b)-0.1255*cos(b)-0.025*sin(b)
  if(ET.instantaneous==TRUE){hours <- hours +0.5}
  hour.angle <- (pi/12)*((hours-0.5+0.06667*(WeatherStation$long*
                                               pi/180-long.z*pi/180)+Sc)-12)
  w1 <- hour.angle-((pi)/24)
  w2 <- hour.angle+((pi)/24)
  hour.angle.s <- acos(-tan(phi)*tan(delta))
  w1c <- ifelse(w1< -hour.angle.s, -hour.angle.s, 
                ifelse(w1>hour.angle.s, hour.angle.s, ifelse(w1>w2, w2, w1)))
  w2c <- ifelse(w2< -hour.angle.s, -hour.angle.s, 
                ifelse(w2>hour.angle.s, hour.angle.s, w2))
  Beta <- asin((sin(phi)*sin(delta)+cos(phi)*cos(delta)*cos(hour.angle)))
  Ra <- ifelse(Beta <= 0, 1e-45, ((12/pi)*4.92*dr)*
                 (((w2c-w1c)*sin(phi)*sin(delta))+(cos(phi)*cos(delta)*
                                                     (sin(w2)-sin(w1)))))
  Rso <- (0.75+2e-5*WeatherStation$elev)*Ra
  Rs.Rso <- ifelse(Rs/Rso<=0.3, 0, ifelse(Rs/Rso>=1, 1, Rs/Rso))
  fcd <- ifelse(1.35*Rs.Rso-0.35<=0.05, 0.05, 
                ifelse(1.35*Rs.Rso-0.35<1, 1.35*Rs.Rso-0.35,1))
  Rn.a <- ((1-0.23)*Rs) - (2.042e-10*fcd*(0.34-0.14*(ea^0.5))*tempK^4)
  Ghr <- c(0.1, 0.5)
  if(ET=="ETr"){Ghr <- c(0.04, 0.2)}
  G.day <- Rn.a * ifelse(Rn.a > 0, Ghr[1], Ghr[2])
  wind.2 <- WeatherStation$wind *(4.87/(log(67.8*WeatherStation$height-5.42)))
  windfactors <- c(0.24, 0.96)
  if(ET=="ETr"){windfactors <-  c(0.25, 1.7)}
  windf <- ifelse(WeatherStation$radiation > 0, windfactors[1], windfactors[2])
  ET.hourly <- ((0.408*Delta*(Rn.a-G.day))+(psi*(37/tempK)*wind.2*(DPV)))/
    (Delta+(psi*(1+(windf*wind.2))))
  if(ET=="ETr"){
    ET.hourly <- ((0.408*Delta*(Rn.a-G.day))+(psi*(66/tempK)*wind.2*(DPV)))/
      (Delta+(psi*(1+(windf*wind.2))))
  }
  return(ET.hourly)
}

#' Calculates daily ET from a surface energy balance and Weather Station 
#' @param Rn             Net radiation. See netRadiation()
#' @param G              Soil Heat Flux. See soilHeatFlux()
#' @param H              Sensible Heat Flux. See calcH()
#' @param Ts             Land surface temperature. See surfaceTemperature()
#' @param WeatherStation WeatherStation data at the flyby from the satellite. 
#' Can be a waterWeatherStation object calculate using read.WSdata and MTL file
#' @param ETr.daily      hourly ETr for every hour of the day. See dailyET()
#' @param C.rad          correction term used in sloping terrain to correct for 
#' variation in 24 h versus instantaneous energy availability. See Allen (2007)
#' @author Guillermo Federico Olmedo
#' @references 
#' R. G. Allen, M. Tasumi, and R. Trezza, "Satellite-based energy balance for mapping evapotranspiration with internalized calibration (METRIC) - Model" Journal of Irrigation and Drainage Engineering, vol. 133, p. 380, 2007
#' @importFrom grDevices colorRampPalette
#' @export
ET24h <- function(Rn, G, H, Ts, WeatherStation, ETr.daily, C.rad=1){
  if(class(WeatherStation)== "waterWeatherStation"){
    WeatherStation <- getDataWS(WeatherStation)
  }
  LE = Rn - G - H
  ET.inst <- 3600*LE/((2.501 - 0.00236 * (Ts - 273.15)) * (1e6))
  ETo.hourly <- hourlyET(WeatherStation, WeatherStation$hours, WeatherStation$DOY)
  ETr.Fr <- ET.inst/ETo.hourly
  ET.24 <- ETr.Fr * ETr.daily * C.rad
  ET.24[ET.24 < 0]  <- 0
  #ET.24[ET.24 > quantile(ET.24, 0.9)] <- quantile(ET.24, 0.9)
  rgb.palette <- grDevices::colorRampPalette(c("red3","snow2","blue"),  space = "rgb")
  print(spplot(ET.24, scales=list(draw=TRUE), xlab="Easting", ylab="Northing", col.regions=rgb.palette, main= "24-Hour Evapotranspiration (mm/day)",
               colorkey=list(height=1), at=seq(0,ceiling(ETr.daily*1.5),length.out=50), maxpixels=ncell(ET.24) * 0.3))
  saveLoadClean(imagestack = ET.24, 
                file = "ET24", overwrite=TRUE)
}


#' Calculates daily ET using Penman Monteith hourly formula for every hour
#' @param WeatherStation a data frame with all the needed fields (see example)
#' @param DOY      day of year
#' @param height   weather station sensors height in meters
#' @param lat      latitude in decimal degrees of the weather station
#' @param long     longitude in decimal degrees of the weather station
#' @param elev     elevation in meters of the weather station
#' @param ET       "ETo" for short crops, similar to clipped, cool-season
#' grass; or "ETr" for tall crops, similar to 0.5 m tall full-cover alfalfa.
#' @param long.z   longitude for local time
#' @param MTL      Metadata file. If not provided will look for one on
#' working directory. If provided or present will calculate weather conditions
#' on satellite overpass.
#' @param date     if date == "auto" will use a MTL file provided or present in 
#' the working folder to select the date.
#' @return ET      daily in mm.h-1
#' @author Guillermo Federico Olmedo
#' @export
#' @references 
#' Allen 2005 ASCE
#' @examples 
#' csvfile <- system.file("extdata", "apples.csv", package="water")
#' 
#' WeatherStation <- read.WSdata(WSdata = csvfile, date.format = "%d/%m/%Y", 
#' lat=-35.42222, long= -71.38639, elev=201, height= 2.2, cf=c(1,0.2777778,1,1))
#' 
#' dailyET(WeatherStation = WeatherStation, lat=-35.422, long=-71.386, elev=124, 
#' ET="ETo")
#' 
dailyET <- function(WeatherStation, DOY, height, lat, long, elev, ET="ETr", 
                    long.z=WeatherStation$long, date = "auto", MTL){
  if(class(WeatherStation)== "waterWeatherStation"){
    if(missing(height)){height <- WeatherStation$location$height}
    if(missing(lat)){lat <- WeatherStation$location$lat}
    if(missing(long)){long <- WeatherStation$location$long}
    if(missing(elev)){elev <- WeatherStation$location$elev}
    
    ## Join with satellite data
    if(missing(MTL)){MTL <- list.files(pattern = "MTL.txt", full.names = T)}
    if(length(MTL)!=0){  ## TODO: Check if MTL date is inside WS data
      MTL <- readLines(MTL, warn=FALSE)
      date.line <- grep("DATE_ACQUIRED",MTL,value=TRUE)
      sat.date <-regmatches(date.line,regexec(text=date.line,
                                              pattern="([0-9]{4})(-)([0-9]{2})(-)([0-9]{2})"))[[1]][1]
      sat.date <- strptime(sat.date,  
                           format = "%Y-%m-%d", tz="GMT")
      if(date == "auto"){
        data <- WeatherStation$hourly[as.character(trunc(WeatherStation$hourly$datetime, "days")) == 
                                        as.character(as.Date(sat.date)),]
      }  else {
        data <- WeatherStation$hourly[as.Date(WeatherStation$hourly$datetime) == as.Date(date),]
      }
      ET.daily <- vector()
      for(i in 1:24){
          ET.daily <- c(ET.daily, hourlyET(data[i,], lat=lat, 
                                           long = long, elev=elev, ET=ET, 
                                           height = height))
        }
    }
    if(length(MTL)==0){
      ET.daily <- vector()
      for(i in 1:24){
        date <- as.POSIXlt(WeatherStation$hourly[i,1], format="%Y-%m-%d %H:%M:%S")
        ET.daily <- c(ET.daily, hourlyET(WeatherStation$hourly[i,], lat=lat, 
                                         long = long, elev=elev, ET=ET, 
                                         height = height))
      }
  }} else {
    print("not implemented yet")}
  return(sum(ET.daily, na.rm=T))
}
midraed/water documentation built on Feb. 20, 2022, 7:30 a.m.