Skip to content

Instantly share code, notes, and snippets.

@scbrown86
Last active February 27, 2020 21:54
Show Gist options
  • Save scbrown86/993f31a336c971a023db160bbed085d2 to your computer and use it in GitHub Desktop.
Save scbrown86/993f31a336c971a023db160bbed085d2 to your computer and use it in GitHub Desktop.
Calculate modified Hargreaves reference ET with gridded data
# 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]]))
@scbrown86
Copy link
Author

Updated method for calculating solar radiation as a function of latitude. Should be faster on fine res data now as the function iterates through rows of the raster.

@scbrown86
Copy link
Author

Updated so that if a z-index is provided to Tmin (see raster::setZ()), the function accounts for the different date order.

@scbrown86
Copy link
Author

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment