Nothing
#' Frequency of frosts
#'
#' @description This function calculates the frequency of frost occurrence, by month.
#' @param mn vector containing the daily minimum temperature series.
#' @param dates sequence of daily dates coinciding with temperature data series.
#' @param thres temperature threshold considered to trigger frost occurrence (0 by default).
#' @param out system path to the generated output file ('.pdf').
#' @param iniday first day of the year ("dd-mm") to consider for first and last occurrence of frost.
#' @param endday last day of the year ("dd-mm") to consider for first and last occurrence of frost.
#' @details Despite the logical threshold of temperature is 0 ÂșC to determine frost occurrence, the argument "thres" is open to change in case of different units of temperature.
#' @return If the output path is defined, a pdf file is created.
#' @examples
#'\donttest{
#' frostFreqs(mn = daily_tmin, dates = seq.Date(as.Date('1981-01-01'),
#' as.Date('2010-12-31'), by ='day'), thres = 0)
#'}
#' @import zoo
#' @import stats
#' @importFrom reshape melt
#' @importFrom grDevices dev.off pdf
#' @export
frostFreqs <- function(mn, dates, thres = 0, out = NULL, iniday = '07-01', endday = '06-30'){
#global variables
value <- month <- days <- NULL
x <- zoo(mn, dates)
nfrosts <- function(x, thres){
freqs <- matrix(0, ncol = 31, nrow = 12)
colnames(freqs) <- 1:31
for(i in 1:12){
d <- as.numeric(x[which(as.numeric(substr(time(x), 6, 7)) == i)])
d[d <= thres] <- 'frost'
d[d != 'frost'] <- 'no frost'
r <- rle(d)
r <- r$lengths[which(r$values == 'frost')]
r <- table(r)
freqs[i,match(as.numeric(names(r)), 1:31)] <- as.numeric(r)
}
freqs
}
years <- unique(substr(dates, 1, 4))
lf <- list()
for(y in 1:length(years)){
lf[[y]] <- nfrosts(x = x[which(years[y] == substr(time(x), 1, 4))],
thres = thres)
}
allf <- abind::abind(lf, along = 3)
allf <- apply(allf, c(1,2), sum, na.rm = T)
# what percentage represents each day for its corresponding month
allf <- t(apply(allf, 1, function(x) {x * 100 / sum(x, na.rm = T)}))
allf <- melt(allf)
names(allf) <- c('month', 'days', 'value')
# frost probability, by month
pf <- numeric()
for(i in 1:12){
w <- which(as.numeric(format(time(x), '%m')) == i)
ec <- ecdf(as.numeric(x[w]))
pf[i] <- ec(thres)
}
for(i in 1:12){
w <- which(allf$month == i)
allf$value[w] <- allf$value[w] * pf[i]
}
cols <- rev(c('#a50026','#d73027','#f46d43','#fdae61','#fee090','#ffffbf','#e0f3f8','#abd9e9','#74add1','#4575b4','#313695'))
colfunc <- colorRampPalette(cols)
allf$days <- factor(allf$days)
for(i in 2:31) allf$days <- relevel(allf$days, i) #reorder levels
aord <- allf[order(allf$days, decreasing = T), ]
aord$month <- factor(aord$month, levels = c(7:12, 1:6))
pg <- ggplot(aord, aes(y = value, x = month, fill = days)) +
# geom_bar(stat = "identity", position = position_nudge(x = 0.5)) +
geom_bar(stat = "identity") +
scale_fill_manual(name = "Frost spell length in days (frequency of occurrence)", values=rev(colfunc(31))) +
scale_x_discrete(name = "", labels = month.abb[c(7:12, 1:6)]) +
scale_y_continuous(name = "Probability (%)", breaks = seq(0,100, 10), limits = c(0, 100)) +
geom_vline(xintercept=seq(0.5, 12.5, 1),
lwd=0.2, colour="grey") +
theme_bw() +
theme(legend.position = "bottom",
legend.spacing.x = unit(-0.1,"cm"),
legend.key.size = unit(1,"lines"),
legend.title=element_text(size = 10),
legend.text=element_text(size = 7),
axis.title.x=element_blank(),
panel.grid.major.x = element_blank(),
axis.ticks.x = element_blank(),
legend.margin=margin(0,0,0,0)) +
guides(fill=guide_legend(ncol = 31,
label.position = 'bottom',
label.hjust = 0,
title.position = "top",
reverse = TRUE))
############
# plot of first and last frost occurrence probability
frostyears <- aggregate(x, by = list(substr(dates, 1, 4)), function(x) length(which(x < thres)))
frostprob <- length(which(frostyears > 0)) / length(frostyears)
# first frost probability
ff <- firstFrost(mn, dates, iniday = iniday, endday = endday, type = 'doy', thres = thres)
ffun <- sort(unique(ff))
ec <- ecdf(ff)
# plot(ec)
pff <- sapply(ffun, ec) * frostprob
# last frost probability
lf <- lastFrost(mn, dates, iniday = iniday, endday = endday, type = 'doy', thres = 0)
lfun <- sort(unique(lf))
ec <- ecdf(lf)
# plot(ec)
plf <- (1-sapply(lfun, ec)) * frostprob
pfrost <- data.frame(day = seq.Date(as.Date('2019-07-01'),
as.Date('2020-06-30'), by ='day'),
first = NA, last = NA)
pfrost$first[ffun] <- pff
pfrost$last[lfun] <- plf
pfrost <- melt(pfrost, id = c('day'))
pf <- ggplot(data=pfrost, aes(x=day, y=value, colour=variable)) +
geom_point() +
geom_line(data = pfrost[!is.na(pfrost$value),]) +
theme_void() +
# theme_bw() +
scale_color_manual(values = c('blue2', 'darkviolet'),
labels = c("Before", "After"),
name="Date of frost probability") +
ylim(0,1) +
# adjust x-axis
xlim(as.Date(c(min(pfrost$day)+15,max(pfrost$day)-15), format="%Y/%m/%d")) +
theme(legend.position = c(0.17, 0.9),
legend.title=element_text(size=10),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
axis.text.y=element_blank(),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()) +
guides(colour=guide_legend(ncol=2))
gt <- combine_plots(pg, pf)
if(!is.null(out)){
pdf(out, height = 6, width = 7, useDingbats = F)
plot(gt)
dev.off()
} else{
plot(gt)
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.