diff --git a/R/petGrid.R b/R/petGrid.R index 7f93e96..ccd8b0b 100644 --- a/R/petGrid.R +++ b/R/petGrid.R @@ -39,12 +39,13 @@ #' @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) { @@ -52,9 +53,10 @@ #' 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) @@ -62,18 +64,21 @@ #' 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, @@ -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. @@ -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") { @@ -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") } - -