Calculate modified Hargreaves reference ET with gridded data

modhargreaves.R

# https://doi.org/10.1023/A:1015508322413
# The refet() function assumes that the rasters are ordered Jan-Dec
# if a different order is required, make sure a z-index is provided for
# the Tmin raster. The refet() function will then calculate
# julian dates and month lengths off this z-index.

library(raster)
library(sirad)

Tmin <- resample(crop(x = getData("worldclim", res = 10, var="tmin"),
                      y = extent(c(-20,160,50,80))),
                 y = raster(res = 1, ext = extent(-20, 160, 50, 80)))
Tmin <- Tmin / 10
Tmax <- resample(crop(x = getData("worldclim", res = 10, var="tmax"),
                      y = extent(c(-20,160,50,80))),
                 y = raster(res = 1, ext = extent(-20, 160, 50, 80)))
Tmax <- Tmax / 10
Prec <- resample(crop(x = getData("worldclim", res = 10, var="prec"),
                      y = extent(c(-20,160,50,80))),
                 y = raster(res = 1, ext = extent(-20, 160, 50, 80)))

refet <- function (Tmin, Tmax, Prec, SolarRad = NULL, seasonal = TRUE) {
  stopifnot(exprs = {nlayers(Tmin) == 12
    nlayers(Tmax) == 12
    nlayers(Prec) == 12})
  stopifnot(exprs = {inherits(Tmin, c("RasterBrick", "RasterStack"))
    inherits(Tmax, c("RasterBrick", "RasterStack"))
    inherits(Prec, c("RasterBrick", "RasterStack"))})
  if (is.null(getZ(Tmin))) {
    warning("No z-index for Tmin. Assuming Jan - Dec")
  }
  require(sirad)
  require(zoo)
  ET0 = Tmin * NA
  Tavg = (Tmin + Tmax)/2
  Trange = Tmax - Tmin
  Trange[Trange < 0] = 0
  ## estimate monthly solar radiation from latitude if not provided
  if (is.null(SolarRad)) {
    ## set julian day if no z-index is provided
    jd = if (is.null(getZ(Tmin))) {
      c(15, 46, 75, 106, 136, 167, 197, 228, 258, 289, 320, 350)
    } else {
      ## else calc julian day from dates
      as.POSIXlt(getZ(Tmin))$yday
    }
    ## raster of latitudes
    SolarRad = brick(nl = 12, nrow = nrow(Tmin), ncols = ncol(Tmin),
                     crs = crs(Tmin), 
                     xmn = xmin(Tmin), xmx = xmax(Tmin), 
                     ymn = ymin(Tmin), ymx = ymax(Tmin))
    SolarRad[] <- NA
    ## for each julian day
    for (i in jd) {
      ## for each row in the raster (same lat!)
      for (j in 1:nrow(SolarRad)) {
        ## extract the cell numbers
        cellNums <- cellFromRow(SolarRad, j)
        ## calculate the latitude
        lat <- xyFromCell(SolarRad, cellNums[1])[2]
        ## calculate the solar radiation
        ra <- as.numeric(sirad::extrat(i = i, lat = sirad::radians(lat))[1])
        ## set the value of all cells in that row to the solar radiation value
        SolarRad[[which(i == jd)]][cellNums][] <- ra
      }
    }
    names(SolarRad) = paste0("SolarRad_", month.abb)
  }
  ab <- Trange - 0.0123 * Prec
  ET0 <- 0.0013 * 0.408 * SolarRad * (Tavg + 17) * ab^0.76
  ET0[is.nan(ab^0.76)] <- 0
  ET0[ET0 < 0] <- 0
  ## month lengths
  ndays <- function(d) {
    x <- zoo::as.yearmon(d)
    d <- as.numeric(zoo::as.Date(x + 1/12) - zoo::as.Date(x))
    d[which(d == 28)] <- 28.25
    return(d)
  }
  ## if no z-index assume Jan-Dec year
  mlen = if (is.null(getZ(Tmin))) {
    c(31, 28.25, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31)
  } else {
    ## else get month lengths from zInd
    sapply(getZ(Tmin), ndays)
  }
  ## Multiply by month lengths to get monthly ET0
  ET0 <- ET0 * mlen
  month.names = if (is.null(getZ(Tmin))) {
    month.abb
  } else {
    format(getZ(Tmin),"%b")
  }
  names(ET0) <- paste0("ET0_", month.names)
  ET0[["AnnAvgET0"]] <- calc(ET0, mean, na.rm = TRUE)
  if (seasonal) {
    djf <- which(month.names %in% c("Dec", "Jan", "Feb"))
    mam <- which(month.names %in% c("Mar", "Apr", "May"))
    jja <- which(month.names %in% c("Jun", "Jul", "Aug"))
    son <- which(month.names %in% c("Sep", "Oct", "Nov"))
    ET0[["DJF"]] <- calc(ET0[[djf]], mean, na.rm = TRUE)
    ET0[["MAM"]] <- calc(ET0[[mam]], mean, na.rm = TRUE)
    ET0[["JJA"]] <- calc(ET0[[jja]], mean, na.rm = TRUE)
    ET0[["SON"]] <- calc(ET0[[son]], mean, na.rm = TRUE)
  }
  return(ET0)
}

# # If using refet() function in a loop it is quicker to
# # pre-compute solar radiation and provide it in as a raster brick
# jd <- as.POSIXlt(seq(as.Date("1980-01-16", format = "%Y-%m-%d"),
#                      as.Date("1980-12-16", format = "%Y-%m-%d"),
#                      by = "month"))$yday
# sr <- brick(nl = 12, nrow = nrow(Tmin), ncols = ncol(Tmin),
#             crs = crs(Tmin), 
#             xmn = xmin(Tmin), xmx = xmax(Tmin), 
#             ymn = ymin(Tmin), ymx = ymax(Tmin))
# sr[] <- NA
# for (i in jd) {
#   for (j in 1:nrow(sr)) {
#     cellNums <- cellFromRow(sr, j)
#     lat <- xyFromCell(sr, cellNums[1])[2]
#     ra <- as.numeric(extrat(i = i, lat = radians(lat))[1])
#     sr[[which(i == jd)]][cellNums][] <- ra
#   }
# }
# names(sr) <- paste0("SolarRad_", month.abb)

ET0 <- refet(Tmin = Tmin, Tmax = Tmax, Prec = Prec, SolarRad = NULL,
             seasonal = TRUE)
spplot(ET0)
testCell <- cbind(c(56),c(55.5))
cellFromXY(ET0, testCell)

## have to subset to 1:12, other layers are annual averages and seasonal
ras_ET0 <- as.numeric(ET0[[1:12]][4397])

calc_ET0 <- as.numeric(SPEI::hargreaves(Tmin = as.numeric(Tmin[4397]),
                                        Tmax = as.numeric(Tmax[4397]), 
                                        # Ra = as.numeric(sr[4397]),
                                        Pre = as.numeric(Prec[4397]),
                                        lat = 55.5))

# Cor = 0.99997. Difference is due to slight difference in 
# solar radiation calc. 
# Calculate sr, and then uncomment Ra line above and cor = 1
cor.test(ras_ET0, calc_ET0)

# # Testing a different month order
# order <- c(6,7,8,9,10,11,12,1,2,3,4,5)
# Tmin <- Tmin[[order]]
# Tmax <- Tmax[[order]]
# Prec <- Prec[[order]]
# Tmin <- setZ(Tmin, seq(as.Date("1980-06-16", format = "%Y-%m-%d"),
#                        as.Date("1981-05-16", format = "%Y-%m-%d"),
#                        by = "month"), "Date")
# Tmin
# ET0_zInd <- refet(Tmin, Tmax, Prec, SolarRad = NULL, seasonal = TRUE)
# spplot(ET0_zInd)
# 
# identical(ET0[["ET0_Jul"]][], ET0_zInd[["ET0_Jul"]][]) ## TRUE
# identical(ET0[["JJA"]][], ET0_zInd[["JJA"]][]) ## TRUE
# spplot(stack(ET0[[7]], ET0_zInd[[2]]))

Updated to return annual averages, and seasonal averages if required.