Skip to content

Commit

Permalink
Adjust leap year length
Browse files Browse the repository at this point in the history
  • Loading branch information
@jbedia
jbedia committed Apr 24, 2019
1 parent a366372 commit 59516c3
Showing 1 changed file with 102 additions and 73 deletions.
175 changes: 102 additions & 73 deletions R/petGrid.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,41 +39,46 @@
#' @importFrom magrittr %>%
#' @author J Bedia
#' @export
#' @examples
#' @examples \donttest{
#' # Thorthwaite requires monthly mean temperature data as input:
#' data("tas.cru.iberia")
#' thpet <- petGrid(tas = tas.cru.iberia, method = "thornthwaite")
#' require(transformeR)
#' require(magrittr)
#' require(visualizeR)
#' # This is the climatology (by seasons)
#' seasons <- list("DJF" = c(12,1,2), "MAM" = 3:5, "JJA" = 6:8, "SON" = 9:11)
#' season.list <- lapply(1:length(seasons), function(x) {
#' subsetGrid(thpet, season = seasons[[x]]) %>%
#' aggregateGrid(aggr.y = list(FUN = "sum")) %>% climatology()
#' })
#' mg <- makeMultiGrid(season.list, skip.temporal.check = TRUE)
#' plotClimatology(mg, names.attr = names(seasons), at = seq(0, 525, 10),
#' backdrop.theme = "coastline",
#' main = "Mean Potential Evapotranspiration Thornthwaite (1981-2010, mm/season)")
#' spatialPlot(mg, names.attr = names(seasons), at = seq(0, 525, 10),
#' backdrop.theme = "coastline",
#' main = "Mean Potential Evapotranspiration Thornthwaite (1981-2010, mm/season)",
#' rev.colors = TRUE)
#' # A typical operation is computing trends
#' # Here we are interested in the JJA PET trend in Iberia:
#' # We consider the annually aggregated PET series (n = 30 years)
#' jja.pet <- subsetGrid(thpet, season = 6:8) %>% aggregateGrid(aggr.y = list(FUN = "sum"))
#' trend.thpet <- climatology(jja.pet, clim.fun = list(FUN = "trend.1D",
#' dates = getRefDates(jja.pet),
#' method = "kendall"))
#' plotClimatology(trend.thpet, backdrop.theme = "countries",
#' main = "Mann-Kendall Tau PET trend (JJA, 1981-2010)")
#' spatialPlot(trend.thpet, backdrop.theme = "countries",
#' main = "Mann-Kendall Tau PET trend (JJA, 1981-2010)",
#' rev.colors = TRUE)
#' pval.estimate <- climatology(jja.pet, clim.fun = list(FUN = "trend.1D",
#' dates = getRefDates(jja.pet),
#' method = "kendall",
#' return.pvalue = TRUE))
#' sig.points <- map.stippling(clim = pval.estimate, threshold = 0.05, condition = "LT",
#' pch = 19, cex = .5, col = "purple")
#' plotClimatology(trend.thpet, backdrop.theme = "countries",
#' main = "Mann-Kendall Tau PET trend (JJA, 1981-2010)",
#' sp.layout = list(sig.points))
#' spatialPlot(trend.thpet, backdrop.theme = "countries",
#' rev.colors = TRUE,
#' main = "Mann-Kendall Tau PET trend (JJA, 1981-2010)",
#' sp.layout = list(sig.points))
#' # See further examples in 'help(trend.1D)'
#' }


petGrid <- function(tasmin = NULL,
Expand Down Expand Up @@ -200,7 +205,8 @@ petGrid.har <- function(tasmin, tasmax, pr, ...) {
#' @importFrom magrittr extract extract2 %>% %<>%
#' @keywords internal
#' @note This function is intended for daily data only.
#' For hourly or shorter periods, see FAO, eq. 28 p47 (not implemented yet)
#' For hourly or shorter periods, see FAO, eq. 28 p47 (not implemented yet).
#' Radiation will be returned in MJ.m-2.day-1.
#' @references Allen, R.G., Pereira, L.S., Raes, D., Smith, M., 2006. FAO Irrigation and Drainage Paper. Crop Evapotranspiration (guidelines for computing crop water requirements) (No. 56). FAO.


Expand All @@ -223,13 +229,17 @@ petGrid.hs <- function(tasmin, tasmax, what) {
lats %<>% expand.grid() %>% extract(,2)
}
lats <- lats * (pi / 180) ## Conversion decimal degrees --> radians
J <- getRefDates(tasmin) %>% as.Date() %>% format("%j") %>% as.integer()
Gsc <- 0.082 ## Solar constant (MJ.m-2.min-1)
ds <- 0.409 * sin(2 * pi * J / 365 - 1.39) ## Solar declination, FAO, eq. 24 (p. 46)
refdates <- getRefDates(tasmin)
J <- refdates %>% as.Date() %>% format("%j") %>% as.integer()
ind <- getYearsAsINDEX(tasmin) %>% which.leap()
year.len <- rep(365, length(refdates))
year.len[ind] <- 366 ## Number of days per year (controlling for leap years)
ds <- 0.409 * sin(2 * pi * J / year.len - 1.39) ## Solar declination, FAO, eq. 24 (p. 46)
n.mem <- getShape(tasmin, "member")
aux.lat <- matrix(lats, nrow = length(J), ncol = length(lats), byrow = TRUE)
ws <- acos(-tan(aux.lat) * tan(ds)) ## Sunset hour angle, FAO, eq. 25 (p. 46)
dr <- 1 + 0.033 * cos(2 * pi * J / 365) ## inverse relative distance Earth-Sun , FAO, eq. 23 (p. 46)
dr <- 1 + 0.033 * cos(2 * pi * J / year.len) ## inverse relative distance Earth-Sun , FAO, eq. 23 (p. 46)
Gsc <- 0.082 ## Solar constant (MJ.m-2.min-1)
Ra <- (24 * 60 * Gsc * dr * (ws * sin(aux.lat) * sin(ds) + cos(aux.lat) * cos(ds) * sin(ws))) / pi ## FAO, eq. 21 (p. 46)
et0.list <- lapply(1:n.mem, function(x) {
if (what == "rad") {
Expand All @@ -251,67 +261,86 @@ petGrid.hs <- function(tasmin, tasmax, what) {



#' @import transformeR
#' @author J. Bedia, M. Iturbide
#' @import transformeR
#' @importFrom magrittr extract extract2 %>% %<>%
#' @keywords internal
#' Allen, R.G., Pereira, L.S., Raes, D., Smith, M., 2006. FAO Irrigation and Drainage Paper. Crop Evapotranspiration (guidelines for computing crop water requirements) (No. 56). FAO.
# #' @import transformeR
# #' @author J. Bedia, M. Iturbide
# #' @import transformeR
# #' @importFrom magrittr extract extract2 %>% %<>%
# #' @keywords internal
# #' Allen, R.G., Pereira, L.S., Raes, D., Smith, M., 2006. FAO Irrigation and Drainage Paper. Crop Evapotranspiration (guidelines for computing crop water requirements) (No. 56). FAO.
#
# # load("~/workspace/gitRepo/fireDanger/ignore/fffdi_vignette/fffi_data.Rdata", verbose = TRUE)
# # tasmin <- tasmin.fin
# # tasmax <- tasmax.fin
#
# petGrid.pm <- function(tasmin, tasmax, elev, u2, HRm, what) {
# tasmin %<>% redim(member = TRUE)
# tasmax %<>% redim(member = TRUE)
# suppressMessages(checkDim(tasmin, tasmax))
# checkTemporalConsistency(tasmin, tasmax)
# if (isFALSE(getTimeResolution(tasmin) == "DD")) {
# stop("Penman-Monteith method is meant for daily data", call. = FALSE)
# }
# if (typeofGrid(tasmin) == "rotated_grid") {
# stop("Rotated grids are unsupported. Please consider regridding", call. = FALSE)
# }
# what = match.arg(what, choices = c("PET", "rad"))
# lats <- getCoordinates(tasmin)
# if (is.data.frame(lats)) {
# lats %<>% extract2("y")
# } else {
# lats %<>% expand.grid() %>% extract(,2)
# }
# lats <- lats * (pi / 180) # Conversion to radians
# refdates <- getRefDates(tasmin)
# J <- refdates %>% as.Date() %>% format("%j") %>% as.integer()
# nmonthdays <- sapply(refdates, FUN = "ndays")
# ind <- getYearsAsINDEX(tasmin) %>% which.leap()
# year.len <- rep(365, length(refdates))
# year.len[ind] <- 366 ## Number of days per year (controlling for leap years)
# Gsc <- 0.082
# ds <- 0.409 * sin(2 * pi * J / year.len - 1.39) ## Solar declination FAO, eq. 24 (p. 46)
# n.mem <- getShape(tasmin, "member")
# aux.lat <- matrix(lats, nrow = length(J), ncol = length(lats), byrow = TRUE)
# ws <- acos(-tan(aux.lat) * tan(ds)) ## Sunset hour angle, FAO eq. 25 (p. 46)
# N <- (24 / pi) * ws ## Daylight hours (maximum possible duration of sunshine hours), FAO, eq. 34 (p. 48)
# dr <- 1 + 0.033 * cos(2 * pi * J / year.len)
# Ra <- (24 * 60 * Gsc * dr * (ws * sin(aux.lat) * sin(ds) + cos(aux.lat) * cos(ds) * sin(ws))) / pi ## FAO, eq. 21 (p. 46)
# et0.list <- lapply(1:n.mem, function(x) {
# if (what == "rad") {
# return(Ra)
# } else {
# P <- 101.3 * ((293 - 0.0065 * elev) / 293) ^ 5.26
# Rs <- (0.25 + 0.5 * (n/N)) * Ra ## Solar radiation using Angstrom formula, FAO, eq. 35 (p. 50)
# Rso <- (0.75 + 2/100000) * Ra
# Rns <- 0.77 * Rs
# o <- 4.903/1000000000
# tn <- subsetGrid(tasmin, members = x, drop = TRUE) %>% extract2("Data") %>% array3Dto2Dmat()
# tx <- subsetGrid(tasmax, members = x, drop = TRUE) %>% extract2("Data") %>% array3Dto2Dmat()
# tm <- (tx + tn) / 2
# lam <- 2.501 - (2.361 / 1000) * tm
# y <- 0.00163 * (P / lam)
# es <- (0.6108 * exp(17.27 * tn / (tn + 237.3)) + 0.6108 * exp(17.27 * tx / (tx + 237.3))) / 2
# ea <- es * HRm / 100
# Rnl <- o * (tm + 273.16) * (0.34 - 0.14 * ea^0.5) * (1.35 * Rs / Rso - 0.35)
# Rn <- Rns - Rnl
# A <- (4098 * es) / (tm + 237.3)^2
# ((0.408 * A * Rn) + (y * 900 * u2 * (es - ea)) / (tm + 273)) / (A + y * (1 + 0.34 * u2)) %>% return()
# }
# })
# return(list("et0.list" = et0.list, "ref.grid" = tasmin))
# }

#' Calculate the number of days of the current month
#' @param d A date (character) in format YYYY-MM-DD...
#' @return The number of days of the current month
#' @references
#' \url{http://stackoverflow.com/questions/6243088/find-out-the-number-of-days-of-a-month-in-r}
#' @keywords internal
#' @note This function is an internal helper of loadeR, from wehere it has been duplicated (bad practice...)
#' @importFrom utils tail

petGrid.pm <- function(tasmin, tasmax, elev, u2, HRm, what) {
tasmin %<>% redim(member = TRUE)
tasmax %<>% redim(member = TRUE)
suppressMessages(checkDim(tasmin, tasmax))
checkTemporalConsistency(tasmin, tasmax)
if (isFALSE(getTimeResolution(tasmin) == "DD")) {
stop("Penman-Monteith method is meant for daily data", call. = FALSE)
}
if (typeofGrid(tasmin) == "rotated_grid") {
stop("Rotated grids are unsupported. Please consider regridding", call. = FALSE)
}
what = match.arg(what, choices = c("PET", "rad"))
lats <- getCoordinates(tasmin)
if (is.data.frame(lats)) {
lats %<>% extract2("y")
} else {
lats %<>% expand.grid() %>% extract(,2)
}
lats <- lats * (pi / 180)
J <- getRefDates(tasmin) %>% as.Date() %>% format("%j") %>% as.integer()
Gsc <- 0.082
ds <- 0.409 * sin(2 * pi * J / 365 - 1.39)
n.mem <- getShape(tasmin, "member")
aux.lat <- matrix(lats, nrow = length(J), ncol = length(lats), byrow = TRUE)
ws <- acos(-tan(aux.lat) * tan(ds))
N <- 24 / pi * ws ## Daylight hours, FAO, eq. 34 (p. 48)
dr <- 1 + 0.033 * cos(2 * pi * J / 365)
Ra <- (24 * 60 * Gsc * dr * (ws * sin(aux.lat) * sin(ds) + cos(aux.lat) * cos(ds) * sin(ws))) / pi
P <- 101.3 * ((293 - 0.0065 * elev) / 293) ^ 5.26
Rs <- (0.25 + 0.5 * (n/N)) * Ra ## Solar radiation using Angstron formula, FAO, eq. 35 (p. 50)
Rso <- (0.75 + 2/100000) * Ra
Rns <- 0.77 * Rs
o <- 4.903/1000000000
et0.list <- lapply(1:n.mem, function(x) {
if (what == "rad") {
return(Ra)
} else {
tn <- subsetGrid(tasmin, members = x, drop = TRUE) %>% extract2("Data") %>% array3Dto2Dmat()
tx <- subsetGrid(tasmax, members = x, drop = TRUE) %>% extract2("Data") %>% array3Dto2Dmat()
tm <- (tx + tn) / 2
lam <- 2.501 - (2.361 / 1000) * tm
y <- 0.00163 * (P / lam)
es <- (0.6108 * exp(17.27 * tn / (tn + 237.3)) + 0.6108 * exp(17.27 * tx / (tx + 237.3))) / 2
ea <- es * HRm / 100
Rnl <- o * (tm + 273.16) * (0.34 - 0.14 * ea^0.5) * (1.35 * Rs / Rso - 0.35)
Rn <- Rns - Rnl
A <- (4098 * es) / (tm + 237.3)^2
((0.408 * A * Rn) + (y * 900 * u2 * (es - ea)) / (tm + 273)) / (A + y * (1 + 0.34 * u2)) %>% return()
}
})
return(list("et0.list" = et0.list, "ref.grid" = tasmin))
ndays <- function(d) {
as.difftime(tail((28:31)[which(!is.na(as.Date(paste0(substr(d, 1, 8), 28:31), '%Y-%m-%d')))], 1), units = "days")
}




0 comments on commit 59516c3

Please sign in to comment.