Nothing
R/Plotting.R
In Evapotranspiration: Modelling Actual, Potential and Reference Crop Evapotranspiration
Defines functions
ETForcings
ETComparison
ETPlot
Documented in
ETComparison
ETForcings
ETPlot
ETPlot <- function(results, type = "Aggregation", OBS, OBSplot = FALSE, Sdate = time(results$ET.Daily)[1], Edate = time(results$ET.Daily)[length(results$ET.Daily)]) { # plot estimations and observations
# Aggregation plot (default)
if (type == "Aggregation") {
if (is.null(Sdate)) {
Sdate = time(results$ET.Daily)[1]
}
if (is.null(Edate)) {
Edate = time(results$ET.Daily)[length(results$ET.Daily)]
}
# Plots - Aggregation
par(ask=FALSE)
plot.new()
x.Date <- as.Date(time(results$ET.Daily[which(time(results$ET.Daily) == as.Date(Sdate)):which(time(results$ET.Daily) == as.Date(Edate))]))
if (!is.null(results$ET.Daily)) {
plot(results$ET.Daily[which(as.Date(time(results$ET.Daily)) == x.Date[1]) : which(as.Date(time(results$ET.Daily)) == x.Date[length(x.Date)])], main = paste("Daily", results$ET_formulation, results$ET_type), xlab = "Year", ylab = list(c(results$ET_type, "mm/day")))
if (OBSplot == TRUE) {
if (!is.null(OBS)) {
if (!is.null(OBS$E_obs.Daily)) {
lines(OBS$E_obs.Daily[which(as.Date(time(results$ET.Daily)) == x.Date[1]) : which(as.Date(time(results$ET.Daily)) == x.Date[length(x.Date)])], main = "Observed evaporation", type = "o", pch = ".", col = "RED", xlab = "Year", ylab = "Observed evaporation mm/day")
legend("topright", inset = .05, c(paste(results$ET_formulation,results$ET_type),"Observed Class-A pan evaporation"), cex = 0.8, col = c("BLACK","RED"), lty = 1)
}
} else {
warning("No observed data available for plotting, only the estimated values are plotted")
}
}
} else {
warning("No daily results obtained from ", results$ET_formulation)
}
par(ask=TRUE)
plot(results$ET.Monthly[which(as.yearmon(time(results$ET.Monthly)) == as.yearmon(x.Date[1])):which(as.yearmon(time(results$ET.Monthly)) == as.yearmon(x.Date[length(x.Date)]))], ylim = c(0,max(results$ET.Monthly)*1.5), main = paste("Monthly", results$ET_formulation, results$ET_type), type = "o", pch = ".", xlab = "Year", ylab = list(c(results$ET_type, "mm/month")))
if (OBSplot == TRUE) {
if (!is.null(OBS)) {
lines(OBS$E_obs.Monthly[which(as.yearmon(time(results$ET.Monthly)) == as.yearmon(x.Date[1])):which(as.yearmon(time(results$ET.Monthly)) == as.yearmon(x.Date[length(x.Date)]))], main = "Observed evaporation", type = "o", pch = ".", ylim = c(0,400), col = "RED", xlab = "Year", ylab = "Observed evaporation mm/month")
legend("topright", inset = .05, c(paste(results$ET_formulation,results$ET_type),"Observed Class-A pan evaporation"), cex = 0.8, col = c("BLACK","RED"), lty = 1)
} else {
warning("No observed data available for plotting, only the estimated values are plotted")
}
}
plot(results$ET.Annual[which(floor(as.numeric(as.yearmon(time(results$ET.Annual)))) == floor(as.numeric(as.yearmon(x.Date[1])))):which(floor(as.numeric(as.yearmon(time(results$ET.Annual)))) == floor(as.numeric(as.yearmon(x.Date[length(x.Date)]))))], main = paste("Annual", results$ET_formulation, results$ET_type), type = "o", xlab = "Year", ylim = c(0,3000), ylab = list(c(results$ET_type, "mm/year")))
if (OBSplot == TRUE) {
if (!is.null(OBS)) {
lines(OBS$E_obs.Annual[which(floor(as.numeric(as.yearmon(time(results$ET.Annual)))) == floor(as.numeric(as.yearmon(x.Date[1])))):which(floor(as.numeric(as.yearmon(time(results$ET.Annual)))) == floor(as.numeric(as.yearmon(x.Date[length(x.Date)]))))], main = "Observed evaporation", type = "o", col = "RED", ylim = c(0,3000), xlab = "Year", ylab = "Observed evaporation mm/year")
legend("topright", inset = .05, c(paste(results$ET_formulation,results$ET_type),"Observed Class-A pan evaporation"), cex = 0.8, col = c("BLACK","RED"), lty = 1)
} else {
warning("No observed data available for plotting, only the estimated values are plotted")
}
}
par(ask=FALSE)
paste("Completed plotting aggregations for results calculated by", results$ET_formulation, "formulation")
}
# Aggregation plot (default)
if (type == "Average") {
Sdate = time(results$ET.Daily)[1]
Edate = time(results$ET.Daily)[length(results$ET.Daily)]
par(ask=FALSE)
plot.new()
# Plots - Average
plot(results$ET.MonthlyAve~unique(1:12), main = paste("Monthly average", results$ET_formulation, results$ET_type), type = "o", ylim = c(0, max(results$ET.MonthlyAve)*1.5), xlab = "Month", ylab = list(c(paste("Monthly average", results$ET_type), "mm/day")))
if (OBSplot == TRUE) {
if (!is.null(OBS)) {
lines(OBS$E_obs.MonthlyAve~unique(1:12), main = "Observed evaporation", type = "o", ylim = c(0,10), col = "RED", xlab = "Month", ylab = "Monthly average observed evaporation mm")
legend("topright", inset = .05, c(paste(results$ET_formulation,results$ET_type),"Observed Class-A pan evaporation"), cex = 0.8, col = c("BLACK","RED"), lty = 1)
} else {
warning("No observed data available for plotting, only the estimated values are plotted")
}
}
par(ask=TRUE)
plot(results$ET.AnnualAve~unique(as.POSIXlt(data$Date.daily)$year + 1900), main = paste("Annual average", results$ET_formulation, results$ET_type), type = "o", ylim = c(0, max(results$ET.AnnualAve)*1.5), xlab = "Year", ylab = list(c(paste("Annual average", results$ET_type), "mm/day")))
if (OBSplot == TRUE) {
if (!is.null(OBS)) {
lines(OBS$E_obs.AnnualAve~unique(as.POSIXlt(OBS$Date.OBS)$year + 1900), main = "Observed evaporation", type = "o", ylim = c(0,10), col = "RED", xlab = "Year", ylab = "Annual average observed evaporation mm")
legend("topright", inset = .05, c(paste(results$ET_formulation,results$ET_type),"Observed Class-A pan evaporation"), cex = 0.8, col = c("BLACK","RED"), lty = 1)
} else {
warning("No observed data available for plotting, only the estimated values are plotted")
}
}
par(ask=FALSE)
paste("Completed plotting averages for results calculated by", results$ET_formulation, "formulation")
}
par(ask=FALSE)
}
#-------------------------------------------------------------------------------------
ETComparison <- function(results1, results2, results3 = NULL, results4 = NULL, results5 = NULL, results6 = NULL, results7 = NULL,
labs, Sdate=NULL, Edate=NULL, type = "Monthly", ylim=rep(NA,2)) { # plot up to 7 estimations
# check number of sets of results to compare and compile list of legend text
if (exists("results1")==F | exists("results2")==F) {
stop("Please provide at least results1 and results2 for producing comparison plot")
}
if (is.null(results1) | is.null(results2)) {
stop("Please provide at least results1 and results2 for producing comparison plot")
}
Ncomp <- 2
legtext <- vector(mode = "character", length = Ncomp)
legtext[1] <- paste(results1$ET_formulation, results1$ET_type, labs[1])
legtext[2] <- paste(results2$ET_formulation, results2$ET_type, labs[2])
if (!is.null(results3)) {
Ncomp <- 3
length(legtext) <- 3
legtext[3] <- paste(results3$ET_formulation, results3$ET_type, labs[3])
}
if (!is.null(results4)) {
Ncomp <- 4
length(legtext) <- 4
legtext[4] <- paste(results4$ET_formulation, results4$ET_type, labs[4])
}
if (!is.null(results5)) {
Ncomp <- 5
length(legtext) <- 5
legtext[5] <- paste(results5$ET_formulation, results5$ET_type, labs[5])
}
if (!is.null(results6)) {
Ncomp <- 6
length(legtext) <- 6
legtext[6] <- paste(results6$ET_formulation, results6$ET_type, labs[6])
}
if (!is.null(results7)) {
Ncomp <- 7
length(legtext) <- 7
legtext[7] <- paste(results7$ET_formulation, results7$ET_type, labs[7])
}
if (length(labs)!= Ncomp) {
stop("Please ensure that the number of labels is the same as the number of sets of results to compare")
}
# Time series
if (type == "Daily") {
# Daily
par(ask=FALSE)
if (is.null(results1$ET.Daily) | is.null(results2$ET.Daily)) {
stop("Unable to compare results because estimations of daily ET are not available in both results1 and results2")
}
if (is.null(Sdate)) {
Sdate = time(results1$ET.Daily)[1]
}
if (is.null(Edate)) {
Edate = time(results1$ET.Daily)[length(results1$ET.Daily)]
}
x.Date <- as.Date(time(results1$ET.Daily[which(time(results1$ET.Daily) == as.Date(Sdate)):which(time(results1$ET.Daily) == as.Date(Edate))]))
comp <- matrix(nrow=length(x.Date),ncol=Ncomp+1)
comp[,1] <- x.Date
comp[,2] <- results1$ET.Daily[which(as.Date(time(results1$ET.Daily)) == x.Date[1]) : which(as.Date(time(results1$ET.Daily)) == x.Date[length(x.Date)])]
comp[,3] <- results2$ET.Daily[which(as.Date(time(results2$ET.Daily)) == x.Date[1]) : which(as.Date(time(results2$ET.Daily)) == x.Date[length(x.Date)])]
if (Ncomp >= 3) {
if (!is.null(results3)) {
if (!is.null(results3$ET.Daily)) {
comp[,4] <- results3$ET.Daily[which(as.Date(time(results3$ET.Daily)) == x.Date[1]) : which(as.Date(time(results3$ET.Daily)) == x.Date[length(x.Date)])]
} else {
warning("No estimated daily ET available for plotting from results4")
comp[,4] <- NULL
}
}
}
if (Ncomp >= 4) {
if (!is.null(results4)) {
if (!is.null(results4$ET.Daily)) {
comp[,5] <- results4$ET.Daily[which(as.Date(time(results4$ET.Daily)) == x.Date[1]) : which(as.Date(time(results4$ET.Daily)) == x.Date[length(x.Date)])]
} else {
warning("No estimated daily ET available for plotting from results5")
comp[,5] <- NULL
}
}
}
if (Ncomp >= 5) {
if (!is.null(results5)) {
if (!is.null(results5$ET.Daily)) {
comp[,6] <- results5$ET.Daily[which(as.Date(time(results5$ET.Daily)) == x.Date[1]) : which(as.Date(time(results5$ET.Daily)) == x.Date[length(x.Date)])]
} else {
warning("No estimated daily ET available for plotting from results6")
comp[,6] <- NULL
}
}
}
if (Ncomp >= 6) {
if (!is.null(results6)) {
if (!is.null(results6$ET.Daily)) {
comp[,7] <- results6$ET.Daily[which(as.Date(time(results6$ET.Daily)) == x.Date[1]) : which(as.Date(time(results6$ET.Daily)) == x.Date[length(x.Date)])]
} else {
warning("No estimated daily ET available for plotting from results7")
comp[,7] <- NULL
}
}
}
if (Ncomp >= 7) {
if (!is.null(results7)) {
if (!is.null(results7$ET.Daily)) {
comp[,8] <- results7$ET.Daily[which(as.Date(time(results7$ET.Daily)) == x.Date[1]) : which(as.Date(time(results7$ET.Daily)) == x.Date[length(x.Date)])]
} else {
warning("No estimated daily ET available for plotting from results8")
comp[,8] <- NULL
}
}
}
if (is.na(ylim[1])) {
ylim[1] <- 0
}
if (is.na(ylim[2])) {
ylim[2] <- max(results1$ET.Daily)*1.5
}
plot(comp[,2]~as.Date(comp[,1]), type = "l", main = "Daily estimations of ET", ylim = ylim, col = 2, xlab = "Year", ylab = "Estimated ET mm/day", cex.lab = 1.3, cex.main = 2.1, par(mar = c(5.1,4.7,4.7,2.1)))
for (i in 1:(Ncomp-1)) {
lines(comp[,i+2]~as.Date(comp[,1]), type = "o", pch = ".", col = i+2)
}
legend("topright", inset = .01, c(paste(legtext[1:Ncomp])), cex = (1), col = c(seq(from=2,to=Ncomp+1)), lty = 1)
# Non-exceedance probability
par(ask=TRUE)
sortcomp <- matrix(nrow=length(x.Date),ncol=Ncomp+2)
sortcomp[,1] <- comp[,1]
for (j in 1:Ncomp) {
sortcomp[,j+1] <- sort(comp[,j+1])
}
sortcomp[,Ncomp+2] <- seq(from = 1/nrow(sortcomp), to = 1, by = 1/nrow(sortcomp))
plot(sortcomp[,2]~sortcomp[,Ncomp+2], type = "l", main = "Non-exceedance probability of the daily ET", xlim = c(0,1), ylim = ylim, col = 2, xlab = "Non-exceedance probability", ylab = "Estimated daily ET mm/day", cex.lab = 1.3, cex.main = 2.1, par(mar = c(5.1,4.7,4.7,2.1)))
for (k in 1:(Ncomp-1)) {
lines(sortcomp[,k+2]~sortcomp[,Ncomp+2], type = "o", pch = ".", col = k+2)
}
legend("topright", inset = .01, c(paste(legtext[1:Ncomp])), cex = (1), col = c(seq(from=2,to=Ncomp+1)), lty = 1)
# Box plot
par(ask=TRUE)
colnames(comp)[2:(Ncomp+1)] <- legtext
boxplot(comp[,2:(Ncomp+1)], col = c(seq(from=2,to=Ncomp+1)), boxwex = 0.5, range = 0, main = "Box plot of the daily estimations of ET", xlab = "", ylab = "Distribution of daily PET in mm", show.names = FALSE, pars = list(ylim = ylim, cex.lab = 1.3, cex.main = 2.1, par(mar = c(5.1,4.7,4.7,2.1))))
legend("topright", inset = .01, c(paste(legtext[1:Ncomp])), cex = (1), fill = c(seq(from=2,to=Ncomp+1)), bty="o")
}
if (type == "Monthly") {
if (is.null(Sdate)) {
Sdate = time(results1$ET.Monthly)[1]
}
if (is.null(Edate)) {
Edate = time(results1$ET.Monthly)[length(results1$ET.Monthly)]
}
# Monthly
par(ask=FALSE)
if (is.null(results1$ET.Monthly) | is.null(results2$ET.Monthly)) {
stop("Unable to compare results because estimations of monthly ET are not available in both results1 and results2")
}
x.Date <- as.yearmon(time(results1$ET.Monthly[which(time(results1$ET.Monthly) == as.yearmon(Sdate)):which(time(results1$ET.Monthly) == as.yearmon(Edate))]))
comp <- matrix(nrow=length(x.Date),ncol=Ncomp+1)
comp[,1] <- x.Date
comp[,2] <- results1$ET.Monthly[which(as.yearmon(time(results1$ET.Monthly)) == x.Date[1]) : which(as.yearmon(time(results1$ET.Monthly)) == x.Date[length(x.Date)])]
comp[,3] <- results2$ET.Monthly[which(as.yearmon(time(results2$ET.Monthly)) == x.Date[1]) : which(as.yearmon(time(results2$ET.Monthly)) == x.Date[length(x.Date)])]
if (Ncomp >= 3) {
if (!is.null(results3)) {
if (!is.null(results3$ET.Monthly)) {
comp[,4] <- results3$ET.Monthly[which(as.yearmon(time(results3$ET.Monthly)) == x.Date[1]) : which(as.yearmon(time(results3$ET.Monthly)) == x.Date[length(x.Date)])]
} else {
warning("No estimated monthly ET available for plotting from results4")
comp[,4] <- NULL
}
}
}
if (Ncomp >= 4) {
if (!is.null(results4)) {
if (!is.null(results4$ET.Monthly)) {
comp[,5] <- results4$ET.Monthly[which(as.yearmon(time(results4$ET.Monthly)) == x.Date[1]) : which(as.yearmon(time(results4$ET.Monthly)) == x.Date[length(x.Date)])]
} else {
warning("No estimated monthly ET available for plotting from results5")
comp[,5] <- NULL
}
}
}
if (Ncomp >= 5) {
if (!is.null(results5)) {
if (!is.null(results5$ET.Monthly)) {
comp[,6] <- results5$ET.Monthly[which(as.yearmon(time(results5$ET.Monthly)) == x.Date[1]) : which(as.yearmon(time(results5$ET.Monthly)) == x.Date[length(x.Date)])]
} else {
warning("No estimated monthly ET available for plotting from results6")
comp[,6] <- NULL
}
}
}
if (Ncomp >= 6) {
if (!is.null(results6)) {
if (!is.null(results6$ET.Monthly)) {
comp[,7] <- results6$ET.Monthly[which(as.yearmon(time(results6$ET.Monthly)) == x.Date[1]) : which(as.yearmon(time(results6$ET.Monthly)) == x.Date[length(x.Date)])]
} else {
warning("No estimated monthly ET available for plotting from results7")
comp[,7] <- NULL
}
}
}
if (Ncomp >= 7) {
if (!is.null(results7)) {
if (!is.null(results7$ET.Monthly)) {
comp[,8] <- results7$ET.Monthly[which(as.yearmon(time(results7$ET.Monthly)) == x.Date[1]) : which(as.yearmon(time(results7$ET.Monthly)) == x.Date[length(x.Date)])]
} else {
warning("No estimated monthly ET available for plotting from results8")
comp[,8] <- NULL
}
}
}
if (is.null(ylim[1])) {
ylim[1] <- 0
}
if (is.null(ylim[2])) {
ylim[2] <- max(results1$ET.Monthly)*1.5
}
plot(comp[,2]~as.yearmon(comp[,1]), type = "l", main = "Monthly estimations of ET", ylim = ylim, col = 2, xlab = "Year", ylab = "Monthly aggregations of estimated ET mm/month", cex.lab = 1.3, cex.main = 2.1, par(mar = c(5.1,4.7,4.7,2.1)))
for (i in 1:(Ncomp-1)) {
lines(comp[,i+2]~as.yearmon(comp[,1]), type = "o", pch = ".", col = i+2)
}
legend("topright", inset = .01, c(paste(legtext[1:Ncomp])), cex = (1), col = c(seq(from=2,to=Ncomp+1)), lty = 1)
# Non-exceedance probability
par(ask=TRUE)
sortcomp <- matrix(nrow=length(x.Date),ncol=Ncomp+2)
sortcomp[,1] <- comp[,1]
for (j in 1:Ncomp) {
sortcomp[,j+1] <- sort(comp[,j+1])
}
sortcomp[,Ncomp+2] <- seq(from = 1/nrow(sortcomp), to = 1, by = 1/nrow(sortcomp))
plot(sortcomp[,2]~sortcomp[,Ncomp+2], type = "l", main = "Non-exceedance probability of the monthly ET", xlim = c(0,1), ylim = ylim, col = 2, xlab = "Non-exceedance probability", ylab = "Monthly aggregations of estimated ET mm/month", cex.lab = 1.3, cex.main = 2.1, par(mar = c(5.1,4.7,4.7,2.1)))
for (k in 1:(Ncomp-1)) {
lines(sortcomp[,k+2]~sortcomp[,Ncomp+2], type = "o", pch = ".", col = k+2)
}
legend("topright", inset = .01, c(paste(legtext[1:Ncomp])), cex = (1), col = c(seq(from=2,to=Ncomp+1)), lty = 1)
# Box plot
par(ask=TRUE)
colnames(comp)[2:(Ncomp+1)] <- legtext
boxplot(comp[,2:(Ncomp+1)], col = c(seq(from=2,to=Ncomp+1)), main = "Box plot of the monthly estimations of ET", boxwex = 0.5, range = 0, xlab = "", ylab = "Distribution of monthly PET in mm", show.names = FALSE, pars=list(ylim = ylim, cex.lab = 1.3, cex.main = 2.1, par(mar = c(5.1,4.7,4.7,2.1))))
legend("topright", inset = .01, c(paste(legtext[1:Ncomp])), cex = (1), fill = c(seq(from=2,to=Ncomp+1)), bty="o")
}
if (type == "Annual") {
if (is.null(Sdate)) {
Sdate = time(results1$ET.Annual)[1]
}
if (is.null(Edate)) {
Edate = time(results1$ET.Annual)[length(results1$ET.Annual)]
}
# Annual
par(ask=FALSE)
if (is.null(results1$ET.Annual) | is.null(results2$ET.Annual)) {
stop("Unable to compare results because estimations of annual ET are not available in both results1 and results2")
}
x.Date <- floor(as.numeric(as.yearmon(time(results1$ET.Annual[which(time(results1$ET.Annual) == floor(as.numeric(as.yearmon(Sdate)))):which(time(results1$ET.Annual) == floor(as.numeric(as.yearmon(Edate))))]))))
comp <- matrix(nrow=length(x.Date),ncol=Ncomp+1)
comp[,1] <- x.Date
comp[,2] <- results1$ET.Annual[which(floor(as.numeric(as.yearmon(time(results1$ET.Annual)))) == x.Date[1]) : which(floor(as.numeric(as.yearmon(time(results1$ET.Annual)))) == x.Date[length(x.Date)])]
comp[,3] <- results2$ET.Annual[which(floor(as.numeric(as.yearmon(time(results2$ET.Annual)))) == x.Date[1]) : which(floor(as.numeric(as.yearmon(time(results2$ET.Annual)))) == x.Date[length(x.Date)])]
if (Ncomp >= 3) {
if (!is.null(results3)) {
if (!is.null(results3$ET.Annual)) {
comp[,4] <- results3$ET.Annual[which(as.Date(time(results3$ET.Annual)) == x.Date[1]) : which(as.Date(time(results3$ET.Annual)) == x.Date[length(x.Date)])]
} else {
warning("No estimated annual ET available for plotting from results4")
comp[,4] <- NULL
}
}
}
if (Ncomp >= 4) {
if (!is.null(results4)) {
if (!is.null(results4$ET.Annual)) {
comp[,5] <- results4$ET.Annual[which(floor(as.numeric(as.yearmon(time(results4$ET.Annual)))) == x.Date[1]) : which(floor(as.numeric(as.yearmon(time(results4$ET.Annual)))) == x.Date[length(x.Date)])]
} else {
warning("No estimated annual ET available for plotting from results5")
comp[,5] <- NULL
}
}
}
if (Ncomp >= 5) {
if (!is.null(results5)) {
if (!is.null(results5$ET.Annual)) {
comp[,6] <- results5$ET.Annual[which(floor(as.numeric(as.yearmon(time(results5$ET.Annual)))) == x.Date[1]) : which(floor(as.numeric(as.yearmon(time(results5$ET.Annual)))) == x.Date[length(x.Date)])]
} else {
warning("No estimated annual ET available for plotting from results6")
comp[,6] <- NULL
}
}
}
if (Ncomp >= 6) {
if (!is.null(results6)) {
if (!is.null(results6$ET.Annual)) {
comp[,7] <- results6$ET.Annual[which(floor(as.numeric(as.yearmon(time(results6$ET.Annual)))) == x.Date[1]) : which(floor(as.numeric(as.yearmon(time(results6$ET.Annual)))) == x.Date[length(x.Date)])]
} else {
warning("No estimated annual ET available for plotting from results7")
comp[,7] <- NULL
}
}
}
if (Ncomp >= 7) {
if (!is.null(results7)) {
if (!is.null(results7$ET.Annual)) {
comp[,8] <- results7$ET.Annual[which(floor(as.numeric(as.yearmon(time(results7$ET.Annual)))) == x.Date[1]) : which(floor(as.numeric(as.yearmon(time(results7$ET.Annual)))) == x.Date[length(x.Date)])]
} else {
warning("No estimated annual ET available for plotting from results8")
comp[,8] <- NULL
}
}
}
if (is.null(ylim[1])) {
ylim[1] <- 0
}
if (is.null(ylim[2])) {
ylim[2] <- max(results1$ET.Annual)*1.5
}
plot(comp[,2]~floor(as.numeric(as.yearmon(comp[,1]))), type = "l", main = "Annual estimations of ET", ylim = ylim, col = 2, xlab = "Year", ylab = "Annual aggregations of estimated ET mm/year", cex.lab = 1.3, cex.main = 2.1, par(mar = c(5.1,4.7,4.7,2.1)))
for (i in 1:(Ncomp-1)) {
lines(comp[,i+2]~floor(as.numeric(as.yearmon(comp[,1]))), type = "l", col = i+2)
}
legend("topright", inset = .01, c(paste(legtext[1:Ncomp])), cex = (1), fill = c(seq(from=2,to=Ncomp+1)), bty="o")
# Non-exceedance probability
par(ask=TRUE)
sortcomp <- matrix(nrow=length(x.Date),ncol=Ncomp+2)
sortcomp[,1] <- comp[,1]
for (j in 1:Ncomp) {
sortcomp[,j+1] <- sort(comp[,j+1])
}
sortcomp[,Ncomp+2] <- seq(from = 1/nrow(sortcomp), to = 1, by = 1/nrow(sortcomp))
plot(sortcomp[,2]~sortcomp[,Ncomp+2], type = "l", main = "Non-exceedance probability of the annual ET", xlim = c(0,1), ylim = ylim, col = 2, xlab = "Non-exceedance probability", ylab = "Annual aggregations of estimated ET mm/year", cex.lab = 1.3, cex.main = 2.1, par(mar = c(5.1,4.7,4.7,2.1)))
for (k in 1:(Ncomp-1)) {
lines(sortcomp[,k+2]~sortcomp[,Ncomp+2], type = "o", pch = ".", col = k+2)
}
legend("topright", inset = .01, c(paste(legtext[1:Ncomp])), cex = (1), col = c(seq(from=2,to=Ncomp+1)), lty = 1)
# Box plot
par(ask=TRUE)
colnames(comp)[2:(Ncomp+1)] <- legtext
boxplot(comp[,2:(Ncomp+1)], col = c(seq(from=2,to=Ncomp+1)), boxwex = 0.5, range = 0, main = "Box plot of the annual estimations of ET", xlab = "", ylab = "Distribution of annual PET in mm", show.names = FALSE, pars = list(ylim = ylim, cex.lab = 1.3, cex.main = 2.1, par(mar = c(5.1,4.7,4.7,2.1))))
legend("topright", inset = .01, c(paste(legtext[1:Ncomp])), cex = (1), fill = c(seq(from=2,to=Ncomp+1)), bty="o")
}
if (type == "MonthlyAve") {
if (is.null(Sdate)) {
Sdate = time(results1$ET.MonthlyAve)[1]
}
if (is.null(Edate)) {
Edate = time(results1$ET.MonthlyAve)[length(results1$ET.MonthlyAve)]
}
# Monthly average
par(ask=FALSE)
if (is.null(results1$ET.MonthlyAve) | is.null(results2$ET.MonthlyAve)) {
stop("Unable to compare results because estimations of monthly averaged ET are not available in both results1 and results2")
}
x.Date <- floor(as.numeric(as.yearmon(time(results1$ET.MonthlyAve[which(time(results1$ET.MonthlyAve) == floor(as.numeric(as.yearmon(Sdate)))):which(time(results1$ET.MonthlyAve) == floor(as.numeric(as.yearmon(Edate))))]))))
comp <- matrix(nrow=length(results1$ET.MonthlyAve),ncol=Ncomp+1)
comp[,1] <- x.Date
comp[,2] <- results1$ET.MonthlyAve
comp[,3] <- results2$ET.MonthlyAve
if (Ncomp >= 3) {
if (!is.null(results3)) {
if (!is.null(results3$ET.MonthlyAve)) {
comp[,4] <- results3$ET.MonthlyAve
} else {
warning("No estimated monthly averaged ET available for plotting from results4")
comp[,4] <- NULL
}
}
}
if (Ncomp >= 4) {
if (!is.null(results4)) {
if (!is.null(results4$ET.MonthlyAve)) {
comp[,5] <- results4$ET.MonthlyAve
} else {
warning("No estimated monthly averaged ET available for plotting from results5")
comp[,5] <- NULL
}
}
}
if (Ncomp >= 5) {
if (!is.null(results5)) {
if (!is.null(results5$ET.MonthlyAve)) {
comp[,6] <- results5$ET.MonthlyAve
} else {
warning("No estimated monthly averaged ET available for plotting from results6")
comp[,6] <- NULL
}
}
}
if (Ncomp >= 6) {
if (!is.null(results6)) {
if (!is.null(results6$ET.MonthlyAve)) {
comp[,7] <- results6$ET.MonthlyAve
} else {
warning("No estimated monthly averaged ET available for plotting from results7")
comp[,7] <- NULL
}
}
}
if (Ncomp >= 7) {
if (!is.null(results7)) {
if (!is.null(results7$ET.MonthlyAve)) {
comp[,8] <- results7$ET.MonthlyAve
} else {
warning("No estimated monthly averaged ET available for plotting from results8")
comp[,8] <- NULL
}
}
}
if (is.null(ylim[1])) {
ylim[1] <- 0
}
if (is.null(ylim[2])) {
ylim[2] <- max(results1$ET.MonthlyAve)*1.5
}
plot(comp[,2]~unique(1:12), type = "l", main = "Monthly average ET", ylim = ylim, col = 2, xlab = "Month", ylab = "Monthly averages of estimated daily ET mm/day", cex.lab = 1.3, cex.main = 2.1, par(mar = c(5.1,4.7,4.7,2.1)))
for (i in 1:(Ncomp-1)) {
lines(comp[,i+2]~unique(1:12), type = "o", pch = ".", col = i+2)
}
legend("topright", inset = .01, c(paste(legtext[1:Ncomp])), cex = (1), col = c(seq(from=2,to=Ncomp+1)), lty = 1)
# Non-exceedance probability
par(ask=TRUE)
sortcomp <- matrix(nrow=length(results1$ET.MonthlyAve),ncol=Ncomp+2)
sortcomp[,1] <- comp[,1]
for (j in 1:Ncomp) {
sortcomp[,j+1] <- sort(comp[,j+1])
}
sortcomp[,Ncomp+2] <- seq(from = 1/nrow(sortcomp), to = 1, by = 1/nrow(sortcomp))
plot(sortcomp[,2]~sortcomp[,Ncomp+2], type = "l", main = "Non-exceedance probability of the monthly average ET", xlim = c(0,1), ylim = ylim, col = 2, xlab = "Non-exceedance probability", ylab = "Monthly averages of estimated daily ET mm/day", cex.lab = 1.3, cex.main = 2.1, par(mar = c(5.1,4.7,4.7,2.1)))
for (k in 1:(Ncomp-1)) {
lines(sortcomp[,k+2]~sortcomp[,Ncomp+2], type = "o", pch = ".", col = k+2)
}
legend("topright", inset = .01, c(paste(legtext[1:Ncomp])), cex = (1), col = c(seq(from=2,to=Ncomp+1)), lty = 1)
# Box plot
par(ask=TRUE)
colnames(comp)[2:(Ncomp+1)] <- legtext
boxplot(comp[,2:(Ncomp+1)], col = c(seq(from=2,to=Ncomp+1)), boxwex = 0.5, range = 0, main = "Box plot of the monthly average ET", xlab = "", ylab = "Distribution of monthly averaged PET in mm/day", show.names = FALSE, pars = list(ylim = ylim, cex.lab = 1.3, cex.main = 2.1, par(mar = c(5.1,4.7,4.7,2.1))))
legend("topright", inset = .01, c(paste(legtext[1:Ncomp])), cex = (1), fill = c(seq(from=2,to=Ncomp+1)), bty="o")
}
if (type == "AnnualAve") {
# Annual average
par(ask=FALSE)
if (is.null(results1$ET.AnnualAve) | is.null(results2$ET.AnnualAve)) {
stop("Unable to compare results because estimations of annually averaged ET are not available in both results1 and results2")
}
x.Date <- floor(as.numeric(as.yearmon(time(results1$ET.AnnualAve[which(time(results1$ET.AnnualAve) == floor(as.numeric(as.yearmon(Sdate)))):which(time(results1$ET.AnnualAve) == floor(as.numeric(as.yearmon(Edate))))]))))
comp <- matrix(nrow=length(results1$ET.AnnualAve),ncol=Ncomp+1)
comp[,1] <- x.Date
comp[,2] <- results1$ET.AnnualAve
comp[,3] <- results2$ET.AnnualAve
if (Ncomp >= 3) {
if (!is.null(results3)) {
if (!is.null(results3$ET.AnnualAve)) {
comp[,4] <- results3$ET.AnnualAve
} else {
warning("No estimated annually averaged ET available for plotting from results4")
comp[,4] <- NULL
}
}
}
if (Ncomp >= 4) {
if (!is.null(results4)) {
if (!is.null(results4$ET.AnnualAve)) {
comp[,5] <- results4$ET.AnnualAve
} else {
warning("No estimated annually averaged ET available for plotting from results5")
comp[,5] <- NULL
}
}
}
if (Ncomp >= 5) {
if (!is.null(results5)) {
if (!is.null(results5$ET.AnnualAve)) {
comp[,6] <- results5$ET.AnnualAve
comp[,6] <- NULL
}
}
}
if (Ncomp >= 6) {
if (!is.null(results6)) {
if (!is.null(results6$ET.AnnualAve)) {
comp[,7] <- results6$ET.AnnualAve
} else {
warning("No estimated annually averaged ET available for plotting from results7")
comp[,7] <- NULL
}
}
}
if (Ncomp >= 7) {
if (!is.null(results7)) {
if (!is.null(results7$ET.AnnualAve)) {
comp[,8] <- results7$ET.AnnualAve
} else {
warning("No estimated annually averaged ET available for plotting from results8")
comp[,8] <- NULL
}
}
}
if (is.null(ylim[1])) {
ylim[1] <- 0
}
if (is.null(ylim[2])) {
ylim[2] <- max(results1$ET.AnnualAve)*1.5
}
plot(comp[,2]~floor(as.numeric(as.yearmon(time(results1$ET.Annual)))), type = "l", main = "Annually average ET", ylim = ylim, col = 2, xlab = "Month", ylab = "Annual averages of estimated daily ET mm/day", cex.lab = 1.3, cex.main = 2.1, par(mar = c(5.1,4.7,4.7,2.1)))
for (i in 1:(Ncomp-1)) {
lines(comp[,i+2]~floor(as.numeric(as.yearmon(time(results1$ET.Annual)))), type = "l", col = i+2)
}
legend("topright", inset = .01, c(paste(legtext[1:Ncomp])), cex = (1), col = c(seq(from=2,to=Ncomp+1)), lty = 1)
# Non-exceedance probability
par(ask=TRUE)
sortcomp <- matrix(nrow=length(results1$ET.AnnualAve),ncol=Ncomp+2)
sortcomp[,1] <- comp[,1]
for (j in 1:Ncomp) {
sortcomp[,j+1] <- sort(comp[,j+1])
}
sortcomp[,Ncomp+2] <- seq(from = 1/nrow(sortcomp), to = 1, by = 1/nrow(sortcomp))
plot(sortcomp[,2]~sortcomp[,Ncomp+2], type = "l", main = "Non-exceedance probability of the annual average ET", xlim = c(0,1), ylim = ylim, col = 2, xlab = "Non-exceedance probability", ylab = "Annual averages of estimated daily ET mm/day", cex.lab = 1.3, cex.main = 2.1, par(mar = c(5.1,4.7,4.7,2.1)))
for (k in 1:(Ncomp-1)) {
lines(sortcomp[,k+2]~sortcomp[,Ncomp+2], type = "o", pch = ".", col = k+2)
}
legend("topright", inset = .01, c(paste(legtext[1:Ncomp])), cex = (1), col = c(seq(from=2,to=Ncomp+1)), lty = 1)
# Box plot
par(ask=TRUE)
colnames(comp)[2:(Ncomp+1)] <- legtext
boxplot(comp[,2:(Ncomp+1)], col = c(seq(from=2,to=Ncomp+1)), boxwex = 0.5, range = 0, main = "Box plot of the annual average ET", xlab = "", ylab = "Distribution ofannually averaged PET in mm/day", show.names = FALSE, pars = list(ylim = ylim, cex.lab = 1.3, cex.main = 2.1, par(mar = c(5.1,4.7,4.7,2.1))))
legend("topright", inset = .01, c(paste(legtext[1:Ncomp])), cex = (1), fill = c(seq(from=2,to=Ncomp+1)), bty="o")
}
par(ask=FALSE)
}
#-------------------------------------------------------------------------------------
ETForcings <- function(data, results, forcing) {
par(ask=FALSE)
plot.new()
# define forcing names
Fnames <- list(v=c("Tmax","Tmin","u2","uz","Rs","n","Precip","Epan","RHmax","RHmin","Tdew"),
n=c("maximum temperature","minimum temperature", "average wind speed at 2m", "average wind speed", "daily solar radiation", "daily sunshine hours",
"precipitation", "Class-A pan evaporation", "maximum relative humidity", "minimum relative humidity", "average dew point temeprature"),
u=c("degree Celcius","degree Celcius", "m/s", "m/s", "MJ.per.sqm.per.day", "hours", "mm/day", "mm/day", "%", "%", "degree Celcius"))
# check forcing value
if (forcing %in% Fnames$v) {
par(mfrow=c(2,2))
# Aggregated PET vs forcing
# Plot daily relationship
if (!is.null(results$ET.Daily)) {
interval <- "daily"
if (!is.null(data[[forcing]])) {
plot(results$ET.Daily~data[[forcing]], main = results$ET_formulation, xlab = paste(interval, Fnames$n[Fnames$v == forcing], Fnames$u[Fnames$v == forcing]), ylab = list(c(results$ET_type, "mm/day")))
} else {
plot(1, type="n", axes=F, xlab="", ylab="")
text(1, paste("no", forcing, "data to plot"))
}
} else {
plot(1, type="n", axes=F, xlab="", ylab="")
text(1, "no daily ET estimation to plot")
}
# Plot monthly relationship
if (!is.null(results$ET.Monthly)) {
interval <- "monthly"
if (!is.null(data[[forcing]])) {
if (forcing != "Precip") {
Fmonthly <- aggregate(data[[forcing]], as.yearmon(data$Date.daily, "%m/%y"), mean)
} else {
Fmonthly <- aggregate(data[[forcing]], as.yearmon(data$Date.daily, "%m/%y"), sum)
}
plot(results$ET.Monthly~Fmonthly, main = results$ET_formulation, xlab = paste(interval, Fnames$n[Fnames$v == forcing], Fnames$u[Fnames$v == forcing]), ylab = list(c(results$ET_type, "mm/month")))
} else {
plot(1, type="n", axes=F, xlab="", ylab="")
text(1, paste("no", forcing, "data to plot"))
}
} else {
plot(1, type="n", axes=F, xlab="", ylab="")
text(1, "no monthly ET estimation to plot")
}
# Plot annual relationship
if (!is.null(results$ET.Annual)) {
interval <- "annual"
if (!is.null(data[[forcing]])) {
if (forcing != "Precip") {
Fannual <- aggregate(Fmonthly, floor(as.numeric(as.yearmon(data$Date.monthly, "%m/%y"))), mean)
} else {
Fannual <- aggregate(Fmonthly, floor(as.numeric(as.yearmon(data$Date.monthly, "%m/%y"))), sum)
}
plot(results$ET.Annual~Fannual, main = results$ET_formulation, xlab = paste(interval, Fnames$n[Fnames$v == forcing], Fnames$u[Fnames$v == forcing]), ylab = list(c(results$ET_type, "mm/year")))
} else {
plot(1, type="n", axes=F, xlab="", ylab="")
text(1, paste("no", forcing, "data to plot"))
}
} else {
plot(1, type="n", axes=F, xlab="", ylab="")
text(1, "no annual ET estimation to plot")
}
} else {
stop("Please select forcing from 'Tmax','Tmin','u2','uz','Rs','n','Precip','Epan','RHmax','RHmin','Tdew'") # forcing is not defined
}
}
#-------------------------------------------------------------------------------------
Try the Evapotranspiration package in your browser
Any scripts or data that you put into this service are public.
Evapotranspiration documentation built on Jan. 10, 2022, 9:06 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions ETForcings ETComparison ETPlot
Documented in ETComparison ETForcings ETPlot
ETPlot <- function(results, type = "Aggregation", OBS, OBSplot = FALSE, Sdate = time(results$ET.Daily)[1], Edate = time(results$ET.Daily)[length(results$ET.Daily)]) { # plot estimations and observations
# Aggregation plot (default)
if (type == "Aggregation") {
if (is.null(Sdate)) {
Sdate = time(results$ET.Daily)[1]
}
if (is.null(Edate)) {
Edate = time(results$ET.Daily)[length(results$ET.Daily)]
}
# Plots - Aggregation
par(ask=FALSE)
plot.new()
x.Date <- as.Date(time(results$ET.Daily[which(time(results$ET.Daily) == as.Date(Sdate)):which(time(results$ET.Daily) == as.Date(Edate))]))
if (!is.null(results$ET.Daily)) {
plot(results$ET.Daily[which(as.Date(time(results$ET.Daily)) == x.Date[1]) : which(as.Date(time(results$ET.Daily)) == x.Date[length(x.Date)])], main = paste("Daily", results$ET_formulation, results$ET_type), xlab = "Year", ylab = list(c(results$ET_type, "mm/day")))
if (OBSplot == TRUE) {
if (!is.null(OBS)) {
if (!is.null(OBS$E_obs.Daily)) {
lines(OBS$E_obs.Daily[which(as.Date(time(results$ET.Daily)) == x.Date[1]) : which(as.Date(time(results$ET.Daily)) == x.Date[length(x.Date)])], main = "Observed evaporation", type = "o", pch = ".", col = "RED", xlab = "Year", ylab = "Observed evaporation mm/day")
legend("topright", inset = .05, c(paste(results$ET_formulation,results$ET_type),"Observed Class-A pan evaporation"), cex = 0.8, col = c("BLACK","RED"), lty = 1)
}
} else {
warning("No observed data available for plotting, only the estimated values are plotted")
}
}
} else {
warning("No daily results obtained from ", results$ET_formulation)
}
par(ask=TRUE)
plot(results$ET.Monthly[which(as.yearmon(time(results$ET.Monthly)) == as.yearmon(x.Date[1])):which(as.yearmon(time(results$ET.Monthly)) == as.yearmon(x.Date[length(x.Date)]))], ylim = c(0,max(results$ET.Monthly)*1.5), main = paste("Monthly", results$ET_formulation, results$ET_type), type = "o", pch = ".", xlab = "Year", ylab = list(c(results$ET_type, "mm/month")))
if (OBSplot == TRUE) {
if (!is.null(OBS)) {
lines(OBS$E_obs.Monthly[which(as.yearmon(time(results$ET.Monthly)) == as.yearmon(x.Date[1])):which(as.yearmon(time(results$ET.Monthly)) == as.yearmon(x.Date[length(x.Date)]))], main = "Observed evaporation", type = "o", pch = ".", ylim = c(0,400), col = "RED", xlab = "Year", ylab = "Observed evaporation mm/month")
legend("topright", inset = .05, c(paste(results$ET_formulation,results$ET_type),"Observed Class-A pan evaporation"), cex = 0.8, col = c("BLACK","RED"), lty = 1)
} else {
warning("No observed data available for plotting, only the estimated values are plotted")
}
}
plot(results$ET.Annual[which(floor(as.numeric(as.yearmon(time(results$ET.Annual)))) == floor(as.numeric(as.yearmon(x.Date[1])))):which(floor(as.numeric(as.yearmon(time(results$ET.Annual)))) == floor(as.numeric(as.yearmon(x.Date[length(x.Date)]))))], main = paste("Annual", results$ET_formulation, results$ET_type), type = "o", xlab = "Year", ylim = c(0,3000), ylab = list(c(results$ET_type, "mm/year")))
if (OBSplot == TRUE) {
if (!is.null(OBS)) {
lines(OBS$E_obs.Annual[which(floor(as.numeric(as.yearmon(time(results$ET.Annual)))) == floor(as.numeric(as.yearmon(x.Date[1])))):which(floor(as.numeric(as.yearmon(time(results$ET.Annual)))) == floor(as.numeric(as.yearmon(x.Date[length(x.Date)]))))], main = "Observed evaporation", type = "o", col = "RED", ylim = c(0,3000), xlab = "Year", ylab = "Observed evaporation mm/year")
legend("topright", inset = .05, c(paste(results$ET_formulation,results$ET_type),"Observed Class-A pan evaporation"), cex = 0.8, col = c("BLACK","RED"), lty = 1)
} else {
warning("No observed data available for plotting, only the estimated values are plotted")
}
}
par(ask=FALSE)
paste("Completed plotting aggregations for results calculated by", results$ET_formulation, "formulation")
}
# Aggregation plot (default)
if (type == "Average") {
Sdate = time(results$ET.Daily)[1]
Edate = time(results$ET.Daily)[length(results$ET.Daily)]
par(ask=FALSE)
plot.new()
# Plots - Average
plot(results$ET.MonthlyAve~unique(1:12), main = paste("Monthly average", results$ET_formulation, results$ET_type), type = "o", ylim = c(0, max(results$ET.MonthlyAve)*1.5), xlab = "Month", ylab = list(c(paste("Monthly average", results$ET_type), "mm/day")))
if (OBSplot == TRUE) {
if (!is.null(OBS)) {
lines(OBS$E_obs.MonthlyAve~unique(1:12), main = "Observed evaporation", type = "o", ylim = c(0,10), col = "RED", xlab = "Month", ylab = "Monthly average observed evaporation mm")
legend("topright", inset = .05, c(paste(results$ET_formulation,results$ET_type),"Observed Class-A pan evaporation"), cex = 0.8, col = c("BLACK","RED"), lty = 1)
} else {
warning("No observed data available for plotting, only the estimated values are plotted")
}
}
par(ask=TRUE)
plot(results$ET.AnnualAve~unique(as.POSIXlt(data$Date.daily)$year + 1900), main = paste("Annual average", results$ET_formulation, results$ET_type), type = "o", ylim = c(0, max(results$ET.AnnualAve)*1.5), xlab = "Year", ylab = list(c(paste("Annual average", results$ET_type), "mm/day")))
if (OBSplot == TRUE) {
if (!is.null(OBS)) {
lines(OBS$E_obs.AnnualAve~unique(as.POSIXlt(OBS$Date.OBS)$year + 1900), main = "Observed evaporation", type = "o", ylim = c(0,10), col = "RED", xlab = "Year", ylab = "Annual average observed evaporation mm")
legend("topright", inset = .05, c(paste(results$ET_formulation,results$ET_type),"Observed Class-A pan evaporation"), cex = 0.8, col = c("BLACK","RED"), lty = 1)
} else {
warning("No observed data available for plotting, only the estimated values are plotted")
}
}
par(ask=FALSE)
paste("Completed plotting averages for results calculated by", results$ET_formulation, "formulation")
}
par(ask=FALSE)
}
#-------------------------------------------------------------------------------------
ETComparison <- function(results1, results2, results3 = NULL, results4 = NULL, results5 = NULL, results6 = NULL, results7 = NULL,
labs, Sdate=NULL, Edate=NULL, type = "Monthly", ylim=rep(NA,2)) { # plot up to 7 estimations
# check number of sets of results to compare and compile list of legend text
if (exists("results1")==F | exists("results2")==F) {
stop("Please provide at least results1 and results2 for producing comparison plot")
}
if (is.null(results1) | is.null(results2)) {
stop("Please provide at least results1 and results2 for producing comparison plot")
}
Ncomp <- 2
legtext <- vector(mode = "character", length = Ncomp)
legtext[1] <- paste(results1$ET_formulation, results1$ET_type, labs[1])
legtext[2] <- paste(results2$ET_formulation, results2$ET_type, labs[2])
if (!is.null(results3)) {
Ncomp <- 3
length(legtext) <- 3
legtext[3] <- paste(results3$ET_formulation, results3$ET_type, labs[3])
}
if (!is.null(results4)) {
Ncomp <- 4
length(legtext) <- 4
legtext[4] <- paste(results4$ET_formulation, results4$ET_type, labs[4])
}
if (!is.null(results5)) {
Ncomp <- 5
length(legtext) <- 5
legtext[5] <- paste(results5$ET_formulation, results5$ET_type, labs[5])
}
if (!is.null(results6)) {
Ncomp <- 6
length(legtext) <- 6
legtext[6] <- paste(results6$ET_formulation, results6$ET_type, labs[6])
}
if (!is.null(results7)) {
Ncomp <- 7
length(legtext) <- 7
legtext[7] <- paste(results7$ET_formulation, results7$ET_type, labs[7])
}
if (length(labs)!= Ncomp) {
stop("Please ensure that the number of labels is the same as the number of sets of results to compare")
}
# Time series
if (type == "Daily") {
# Daily
par(ask=FALSE)
if (is.null(results1$ET.Daily) | is.null(results2$ET.Daily)) {
stop("Unable to compare results because estimations of daily ET are not available in both results1 and results2")
}
if (is.null(Sdate)) {
Sdate = time(results1$ET.Daily)[1]
}
if (is.null(Edate)) {
Edate = time(results1$ET.Daily)[length(results1$ET.Daily)]
}
x.Date <- as.Date(time(results1$ET.Daily[which(time(results1$ET.Daily) == as.Date(Sdate)):which(time(results1$ET.Daily) == as.Date(Edate))]))
comp <- matrix(nrow=length(x.Date),ncol=Ncomp+1)
comp[,1] <- x.Date
comp[,2] <- results1$ET.Daily[which(as.Date(time(results1$ET.Daily)) == x.Date[1]) : which(as.Date(time(results1$ET.Daily)) == x.Date[length(x.Date)])]
comp[,3] <- results2$ET.Daily[which(as.Date(time(results2$ET.Daily)) == x.Date[1]) : which(as.Date(time(results2$ET.Daily)) == x.Date[length(x.Date)])]
if (Ncomp >= 3) {
if (!is.null(results3)) {
if (!is.null(results3$ET.Daily)) {
comp[,4] <- results3$ET.Daily[which(as.Date(time(results3$ET.Daily)) == x.Date[1]) : which(as.Date(time(results3$ET.Daily)) == x.Date[length(x.Date)])]
} else {
warning("No estimated daily ET available for plotting from results4")
comp[,4] <- NULL
}
}
}
if (Ncomp >= 4) {
if (!is.null(results4)) {
if (!is.null(results4$ET.Daily)) {
comp[,5] <- results4$ET.Daily[which(as.Date(time(results4$ET.Daily)) == x.Date[1]) : which(as.Date(time(results4$ET.Daily)) == x.Date[length(x.Date)])]
} else {
warning("No estimated daily ET available for plotting from results5")
comp[,5] <- NULL
}
}
}
if (Ncomp >= 5) {
if (!is.null(results5)) {
if (!is.null(results5$ET.Daily)) {
comp[,6] <- results5$ET.Daily[which(as.Date(time(results5$ET.Daily)) == x.Date[1]) : which(as.Date(time(results5$ET.Daily)) == x.Date[length(x.Date)])]
} else {
warning("No estimated daily ET available for plotting from results6")
comp[,6] <- NULL
}
}
}
if (Ncomp >= 6) {
if (!is.null(results6)) {
if (!is.null(results6$ET.Daily)) {
comp[,7] <- results6$ET.Daily[which(as.Date(time(results6$ET.Daily)) == x.Date[1]) : which(as.Date(time(results6$ET.Daily)) == x.Date[length(x.Date)])]
} else {
warning("No estimated daily ET available for plotting from results7")
comp[,7] <- NULL
}
}
}
if (Ncomp >= 7) {
if (!is.null(results7)) {
if (!is.null(results7$ET.Daily)) {
comp[,8] <- results7$ET.Daily[which(as.Date(time(results7$ET.Daily)) == x.Date[1]) : which(as.Date(time(results7$ET.Daily)) == x.Date[length(x.Date)])]
} else {
warning("No estimated daily ET available for plotting from results8")
comp[,8] <- NULL
}
}
}
if (is.na(ylim[1])) {
ylim[1] <- 0
}
if (is.na(ylim[2])) {
ylim[2] <- max(results1$ET.Daily)*1.5
}
plot(comp[,2]~as.Date(comp[,1]), type = "l", main = "Daily estimations of ET", ylim = ylim, col = 2, xlab = "Year", ylab = "Estimated ET mm/day", cex.lab = 1.3, cex.main = 2.1, par(mar = c(5.1,4.7,4.7,2.1)))
for (i in 1:(Ncomp-1)) {
lines(comp[,i+2]~as.Date(comp[,1]), type = "o", pch = ".", col = i+2)
}
legend("topright", inset = .01, c(paste(legtext[1:Ncomp])), cex = (1), col = c(seq(from=2,to=Ncomp+1)), lty = 1)
# Non-exceedance probability
par(ask=TRUE)
sortcomp <- matrix(nrow=length(x.Date),ncol=Ncomp+2)
sortcomp[,1] <- comp[,1]
for (j in 1:Ncomp) {
sortcomp[,j+1] <- sort(comp[,j+1])
}
sortcomp[,Ncomp+2] <- seq(from = 1/nrow(sortcomp), to = 1, by = 1/nrow(sortcomp))
plot(sortcomp[,2]~sortcomp[,Ncomp+2], type = "l", main = "Non-exceedance probability of the daily ET", xlim = c(0,1), ylim = ylim, col = 2, xlab = "Non-exceedance probability", ylab = "Estimated daily ET mm/day", cex.lab = 1.3, cex.main = 2.1, par(mar = c(5.1,4.7,4.7,2.1)))
for (k in 1:(Ncomp-1)) {
lines(sortcomp[,k+2]~sortcomp[,Ncomp+2], type = "o", pch = ".", col = k+2)
}
legend("topright", inset = .01, c(paste(legtext[1:Ncomp])), cex = (1), col = c(seq(from=2,to=Ncomp+1)), lty = 1)
# Box plot
par(ask=TRUE)
colnames(comp)[2:(Ncomp+1)] <- legtext
boxplot(comp[,2:(Ncomp+1)], col = c(seq(from=2,to=Ncomp+1)), boxwex = 0.5, range = 0, main = "Box plot of the daily estimations of ET", xlab = "", ylab = "Distribution of daily PET in mm", show.names = FALSE, pars = list(ylim = ylim, cex.lab = 1.3, cex.main = 2.1, par(mar = c(5.1,4.7,4.7,2.1))))
legend("topright", inset = .01, c(paste(legtext[1:Ncomp])), cex = (1), fill = c(seq(from=2,to=Ncomp+1)), bty="o")
}
if (type == "Monthly") {
if (is.null(Sdate)) {
Sdate = time(results1$ET.Monthly)[1]
}
if (is.null(Edate)) {
Edate = time(results1$ET.Monthly)[length(results1$ET.Monthly)]
}
# Monthly
par(ask=FALSE)
if (is.null(results1$ET.Monthly) | is.null(results2$ET.Monthly)) {
stop("Unable to compare results because estimations of monthly ET are not available in both results1 and results2")
}
x.Date <- as.yearmon(time(results1$ET.Monthly[which(time(results1$ET.Monthly) == as.yearmon(Sdate)):which(time(results1$ET.Monthly) == as.yearmon(Edate))]))
comp <- matrix(nrow=length(x.Date),ncol=Ncomp+1)
comp[,1] <- x.Date
comp[,2] <- results1$ET.Monthly[which(as.yearmon(time(results1$ET.Monthly)) == x.Date[1]) : which(as.yearmon(time(results1$ET.Monthly)) == x.Date[length(x.Date)])]
comp[,3] <- results2$ET.Monthly[which(as.yearmon(time(results2$ET.Monthly)) == x.Date[1]) : which(as.yearmon(time(results2$ET.Monthly)) == x.Date[length(x.Date)])]
if (Ncomp >= 3) {
if (!is.null(results3)) {
if (!is.null(results3$ET.Monthly)) {
comp[,4] <- results3$ET.Monthly[which(as.yearmon(time(results3$ET.Monthly)) == x.Date[1]) : which(as.yearmon(time(results3$ET.Monthly)) == x.Date[length(x.Date)])]
} else {
warning("No estimated monthly ET available for plotting from results4")
comp[,4] <- NULL
}
}
}
if (Ncomp >= 4) {
if (!is.null(results4)) {
if (!is.null(results4$ET.Monthly)) {
comp[,5] <- results4$ET.Monthly[which(as.yearmon(time(results4$ET.Monthly)) == x.Date[1]) : which(as.yearmon(time(results4$ET.Monthly)) == x.Date[length(x.Date)])]
} else {
warning("No estimated monthly ET available for plotting from results5")
comp[,5] <- NULL
}
}
}
if (Ncomp >= 5) {
if (!is.null(results5)) {
if (!is.null(results5$ET.Monthly)) {
comp[,6] <- results5$ET.Monthly[which(as.yearmon(time(results5$ET.Monthly)) == x.Date[1]) : which(as.yearmon(time(results5$ET.Monthly)) == x.Date[length(x.Date)])]
} else {
warning("No estimated monthly ET available for plotting from results6")
comp[,6] <- NULL
}
}
}
if (Ncomp >= 6) {
if (!is.null(results6)) {
if (!is.null(results6$ET.Monthly)) {
comp[,7] <- results6$ET.Monthly[which(as.yearmon(time(results6$ET.Monthly)) == x.Date[1]) : which(as.yearmon(time(results6$ET.Monthly)) == x.Date[length(x.Date)])]
} else {
warning("No estimated monthly ET available for plotting from results7")
comp[,7] <- NULL
}
}
}
if (Ncomp >= 7) {
if (!is.null(results7)) {
if (!is.null(results7$ET.Monthly)) {
comp[,8] <- results7$ET.Monthly[which(as.yearmon(time(results7$ET.Monthly)) == x.Date[1]) : which(as.yearmon(time(results7$ET.Monthly)) == x.Date[length(x.Date)])]
} else {
warning("No estimated monthly ET available for plotting from results8")
comp[,8] <- NULL
}
}
}
if (is.null(ylim[1])) {
ylim[1] <- 0
}
if (is.null(ylim[2])) {
ylim[2] <- max(results1$ET.Monthly)*1.5
}
plot(comp[,2]~as.yearmon(comp[,1]), type = "l", main = "Monthly estimations of ET", ylim = ylim, col = 2, xlab = "Year", ylab = "Monthly aggregations of estimated ET mm/month", cex.lab = 1.3, cex.main = 2.1, par(mar = c(5.1,4.7,4.7,2.1)))
for (i in 1:(Ncomp-1)) {
lines(comp[,i+2]~as.yearmon(comp[,1]), type = "o", pch = ".", col = i+2)
}
legend("topright", inset = .01, c(paste(legtext[1:Ncomp])), cex = (1), col = c(seq(from=2,to=Ncomp+1)), lty = 1)
# Non-exceedance probability
par(ask=TRUE)
sortcomp <- matrix(nrow=length(x.Date),ncol=Ncomp+2)
sortcomp[,1] <- comp[,1]
for (j in 1:Ncomp) {
sortcomp[,j+1] <- sort(comp[,j+1])
}
sortcomp[,Ncomp+2] <- seq(from = 1/nrow(sortcomp), to = 1, by = 1/nrow(sortcomp))
plot(sortcomp[,2]~sortcomp[,Ncomp+2], type = "l", main = "Non-exceedance probability of the monthly ET", xlim = c(0,1), ylim = ylim, col = 2, xlab = "Non-exceedance probability", ylab = "Monthly aggregations of estimated ET mm/month", cex.lab = 1.3, cex.main = 2.1, par(mar = c(5.1,4.7,4.7,2.1)))
for (k in 1:(Ncomp-1)) {
lines(sortcomp[,k+2]~sortcomp[,Ncomp+2], type = "o", pch = ".", col = k+2)
}
legend("topright", inset = .01, c(paste(legtext[1:Ncomp])), cex = (1), col = c(seq(from=2,to=Ncomp+1)), lty = 1)
# Box plot
par(ask=TRUE)
colnames(comp)[2:(Ncomp+1)] <- legtext
boxplot(comp[,2:(Ncomp+1)], col = c(seq(from=2,to=Ncomp+1)), main = "Box plot of the monthly estimations of ET", boxwex = 0.5, range = 0, xlab = "", ylab = "Distribution of monthly PET in mm", show.names = FALSE, pars=list(ylim = ylim, cex.lab = 1.3, cex.main = 2.1, par(mar = c(5.1,4.7,4.7,2.1))))
legend("topright", inset = .01, c(paste(legtext[1:Ncomp])), cex = (1), fill = c(seq(from=2,to=Ncomp+1)), bty="o")
}
if (type == "Annual") {
if (is.null(Sdate)) {
Sdate = time(results1$ET.Annual)[1]
}
if (is.null(Edate)) {
Edate = time(results1$ET.Annual)[length(results1$ET.Annual)]
}
# Annual
par(ask=FALSE)
if (is.null(results1$ET.Annual) | is.null(results2$ET.Annual)) {
stop("Unable to compare results because estimations of annual ET are not available in both results1 and results2")
}
x.Date <- floor(as.numeric(as.yearmon(time(results1$ET.Annual[which(time(results1$ET.Annual) == floor(as.numeric(as.yearmon(Sdate)))):which(time(results1$ET.Annual) == floor(as.numeric(as.yearmon(Edate))))]))))
comp <- matrix(nrow=length(x.Date),ncol=Ncomp+1)
comp[,1] <- x.Date
comp[,2] <- results1$ET.Annual[which(floor(as.numeric(as.yearmon(time(results1$ET.Annual)))) == x.Date[1]) : which(floor(as.numeric(as.yearmon(time(results1$ET.Annual)))) == x.Date[length(x.Date)])]
comp[,3] <- results2$ET.Annual[which(floor(as.numeric(as.yearmon(time(results2$ET.Annual)))) == x.Date[1]) : which(floor(as.numeric(as.yearmon(time(results2$ET.Annual)))) == x.Date[length(x.Date)])]
if (Ncomp >= 3) {
if (!is.null(results3)) {
if (!is.null(results3$ET.Annual)) {
comp[,4] <- results3$ET.Annual[which(as.Date(time(results3$ET.Annual)) == x.Date[1]) : which(as.Date(time(results3$ET.Annual)) == x.Date[length(x.Date)])]
} else {
warning("No estimated annual ET available for plotting from results4")
comp[,4] <- NULL
}
}
}
if (Ncomp >= 4) {
if (!is.null(results4)) {
if (!is.null(results4$ET.Annual)) {
comp[,5] <- results4$ET.Annual[which(floor(as.numeric(as.yearmon(time(results4$ET.Annual)))) == x.Date[1]) : which(floor(as.numeric(as.yearmon(time(results4$ET.Annual)))) == x.Date[length(x.Date)])]
} else {
warning("No estimated annual ET available for plotting from results5")
comp[,5] <- NULL
}
}
}
if (Ncomp >= 5) {
if (!is.null(results5)) {
if (!is.null(results5$ET.Annual)) {
comp[,6] <- results5$ET.Annual[which(floor(as.numeric(as.yearmon(time(results5$ET.Annual)))) == x.Date[1]) : which(floor(as.numeric(as.yearmon(time(results5$ET.Annual)))) == x.Date[length(x.Date)])]
} else {
warning("No estimated annual ET available for plotting from results6")
comp[,6] <- NULL
}
}
}
if (Ncomp >= 6) {
if (!is.null(results6)) {
if (!is.null(results6$ET.Annual)) {
comp[,7] <- results6$ET.Annual[which(floor(as.numeric(as.yearmon(time(results6$ET.Annual)))) == x.Date[1]) : which(floor(as.numeric(as.yearmon(time(results6$ET.Annual)))) == x.Date[length(x.Date)])]
} else {
warning("No estimated annual ET available for plotting from results7")
comp[,7] <- NULL
}
}
}
if (Ncomp >= 7) {
if (!is.null(results7)) {
if (!is.null(results7$ET.Annual)) {
comp[,8] <- results7$ET.Annual[which(floor(as.numeric(as.yearmon(time(results7$ET.Annual)))) == x.Date[1]) : which(floor(as.numeric(as.yearmon(time(results7$ET.Annual)))) == x.Date[length(x.Date)])]
} else {
warning("No estimated annual ET available for plotting from results8")
comp[,8] <- NULL
}
}
}
if (is.null(ylim[1])) {
ylim[1] <- 0
}
if (is.null(ylim[2])) {
ylim[2] <- max(results1$ET.Annual)*1.5
}
plot(comp[,2]~floor(as.numeric(as.yearmon(comp[,1]))), type = "l", main = "Annual estimations of ET", ylim = ylim, col = 2, xlab = "Year", ylab = "Annual aggregations of estimated ET mm/year", cex.lab = 1.3, cex.main = 2.1, par(mar = c(5.1,4.7,4.7,2.1)))
for (i in 1:(Ncomp-1)) {
lines(comp[,i+2]~floor(as.numeric(as.yearmon(comp[,1]))), type = "l", col = i+2)
}
legend("topright", inset = .01, c(paste(legtext[1:Ncomp])), cex = (1), fill = c(seq(from=2,to=Ncomp+1)), bty="o")
# Non-exceedance probability
par(ask=TRUE)
sortcomp <- matrix(nrow=length(x.Date),ncol=Ncomp+2)
sortcomp[,1] <- comp[,1]
for (j in 1:Ncomp) {
sortcomp[,j+1] <- sort(comp[,j+1])
}
sortcomp[,Ncomp+2] <- seq(from = 1/nrow(sortcomp), to = 1, by = 1/nrow(sortcomp))
plot(sortcomp[,2]~sortcomp[,Ncomp+2], type = "l", main = "Non-exceedance probability of the annual ET", xlim = c(0,1), ylim = ylim, col = 2, xlab = "Non-exceedance probability", ylab = "Annual aggregations of estimated ET mm/year", cex.lab = 1.3, cex.main = 2.1, par(mar = c(5.1,4.7,4.7,2.1)))
for (k in 1:(Ncomp-1)) {
lines(sortcomp[,k+2]~sortcomp[,Ncomp+2], type = "o", pch = ".", col = k+2)
}
legend("topright", inset = .01, c(paste(legtext[1:Ncomp])), cex = (1), col = c(seq(from=2,to=Ncomp+1)), lty = 1)
# Box plot
par(ask=TRUE)
colnames(comp)[2:(Ncomp+1)] <- legtext
boxplot(comp[,2:(Ncomp+1)], col = c(seq(from=2,to=Ncomp+1)), boxwex = 0.5, range = 0, main = "Box plot of the annual estimations of ET", xlab = "", ylab = "Distribution of annual PET in mm", show.names = FALSE, pars = list(ylim = ylim, cex.lab = 1.3, cex.main = 2.1, par(mar = c(5.1,4.7,4.7,2.1))))
legend("topright", inset = .01, c(paste(legtext[1:Ncomp])), cex = (1), fill = c(seq(from=2,to=Ncomp+1)), bty="o")
}
if (type == "MonthlyAve") {
if (is.null(Sdate)) {
Sdate = time(results1$ET.MonthlyAve)[1]
}
if (is.null(Edate)) {
Edate = time(results1$ET.MonthlyAve)[length(results1$ET.MonthlyAve)]
}
# Monthly average
par(ask=FALSE)
if (is.null(results1$ET.MonthlyAve) | is.null(results2$ET.MonthlyAve)) {
stop("Unable to compare results because estimations of monthly averaged ET are not available in both results1 and results2")
}
x.Date <- floor(as.numeric(as.yearmon(time(results1$ET.MonthlyAve[which(time(results1$ET.MonthlyAve) == floor(as.numeric(as.yearmon(Sdate)))):which(time(results1$ET.MonthlyAve) == floor(as.numeric(as.yearmon(Edate))))]))))
comp <- matrix(nrow=length(results1$ET.MonthlyAve),ncol=Ncomp+1)
comp[,1] <- x.Date
comp[,2] <- results1$ET.MonthlyAve
comp[,3] <- results2$ET.MonthlyAve
if (Ncomp >= 3) {
if (!is.null(results3)) {
if (!is.null(results3$ET.MonthlyAve)) {
comp[,4] <- results3$ET.MonthlyAve
} else {
warning("No estimated monthly averaged ET available for plotting from results4")
comp[,4] <- NULL
}
}
}
if (Ncomp >= 4) {
if (!is.null(results4)) {
if (!is.null(results4$ET.MonthlyAve)) {
comp[,5] <- results4$ET.MonthlyAve
} else {
warning("No estimated monthly averaged ET available for plotting from results5")
comp[,5] <- NULL
}
}
}
if (Ncomp >= 5) {
if (!is.null(results5)) {
if (!is.null(results5$ET.MonthlyAve)) {
comp[,6] <- results5$ET.MonthlyAve
} else {
warning("No estimated monthly averaged ET available for plotting from results6")
comp[,6] <- NULL
}
}
}
if (Ncomp >= 6) {
if (!is.null(results6)) {
if (!is.null(results6$ET.MonthlyAve)) {
comp[,7] <- results6$ET.MonthlyAve
} else {
warning("No estimated monthly averaged ET available for plotting from results7")
comp[,7] <- NULL
}
}
}
if (Ncomp >= 7) {
if (!is.null(results7)) {
if (!is.null(results7$ET.MonthlyAve)) {
comp[,8] <- results7$ET.MonthlyAve
} else {
warning("No estimated monthly averaged ET available for plotting from results8")
comp[,8] <- NULL
}
}
}
if (is.null(ylim[1])) {
ylim[1] <- 0
}
if (is.null(ylim[2])) {
ylim[2] <- max(results1$ET.MonthlyAve)*1.5
}
plot(comp[,2]~unique(1:12), type = "l", main = "Monthly average ET", ylim = ylim, col = 2, xlab = "Month", ylab = "Monthly averages of estimated daily ET mm/day", cex.lab = 1.3, cex.main = 2.1, par(mar = c(5.1,4.7,4.7,2.1)))
for (i in 1:(Ncomp-1)) {
lines(comp[,i+2]~unique(1:12), type = "o", pch = ".", col = i+2)
}
legend("topright", inset = .01, c(paste(legtext[1:Ncomp])), cex = (1), col = c(seq(from=2,to=Ncomp+1)), lty = 1)
# Non-exceedance probability
par(ask=TRUE)
sortcomp <- matrix(nrow=length(results1$ET.MonthlyAve),ncol=Ncomp+2)
sortcomp[,1] <- comp[,1]
for (j in 1:Ncomp) {
sortcomp[,j+1] <- sort(comp[,j+1])
}
sortcomp[,Ncomp+2] <- seq(from = 1/nrow(sortcomp), to = 1, by = 1/nrow(sortcomp))
plot(sortcomp[,2]~sortcomp[,Ncomp+2], type = "l", main = "Non-exceedance probability of the monthly average ET", xlim = c(0,1), ylim = ylim, col = 2, xlab = "Non-exceedance probability", ylab = "Monthly averages of estimated daily ET mm/day", cex.lab = 1.3, cex.main = 2.1, par(mar = c(5.1,4.7,4.7,2.1)))
for (k in 1:(Ncomp-1)) {
lines(sortcomp[,k+2]~sortcomp[,Ncomp+2], type = "o", pch = ".", col = k+2)
}
legend("topright", inset = .01, c(paste(legtext[1:Ncomp])), cex = (1), col = c(seq(from=2,to=Ncomp+1)), lty = 1)
# Box plot
par(ask=TRUE)
colnames(comp)[2:(Ncomp+1)] <- legtext
boxplot(comp[,2:(Ncomp+1)], col = c(seq(from=2,to=Ncomp+1)), boxwex = 0.5, range = 0, main = "Box plot of the monthly average ET", xlab = "", ylab = "Distribution of monthly averaged PET in mm/day", show.names = FALSE, pars = list(ylim = ylim, cex.lab = 1.3, cex.main = 2.1, par(mar = c(5.1,4.7,4.7,2.1))))
legend("topright", inset = .01, c(paste(legtext[1:Ncomp])), cex = (1), fill = c(seq(from=2,to=Ncomp+1)), bty="o")
}
if (type == "AnnualAve") {
# Annual average
par(ask=FALSE)
if (is.null(results1$ET.AnnualAve) | is.null(results2$ET.AnnualAve)) {
stop("Unable to compare results because estimations of annually averaged ET are not available in both results1 and results2")
}
x.Date <- floor(as.numeric(as.yearmon(time(results1$ET.AnnualAve[which(time(results1$ET.AnnualAve) == floor(as.numeric(as.yearmon(Sdate)))):which(time(results1$ET.AnnualAve) == floor(as.numeric(as.yearmon(Edate))))]))))
comp <- matrix(nrow=length(results1$ET.AnnualAve),ncol=Ncomp+1)
comp[,1] <- x.Date
comp[,2] <- results1$ET.AnnualAve
comp[,3] <- results2$ET.AnnualAve
if (Ncomp >= 3) {
if (!is.null(results3)) {
if (!is.null(results3$ET.AnnualAve)) {
comp[,4] <- results3$ET.AnnualAve
} else {
warning("No estimated annually averaged ET available for plotting from results4")
comp[,4] <- NULL
}
}
}
if (Ncomp >= 4) {
if (!is.null(results4)) {
if (!is.null(results4$ET.AnnualAve)) {
comp[,5] <- results4$ET.AnnualAve
} else {
warning("No estimated annually averaged ET available for plotting from results5")
comp[,5] <- NULL
}
}
}
if (Ncomp >= 5) {
if (!is.null(results5)) {
if (!is.null(results5$ET.AnnualAve)) {
comp[,6] <- results5$ET.AnnualAve
comp[,6] <- NULL
}
}
}
if (Ncomp >= 6) {
if (!is.null(results6)) {
if (!is.null(results6$ET.AnnualAve)) {
comp[,7] <- results6$ET.AnnualAve
} else {
warning("No estimated annually averaged ET available for plotting from results7")
comp[,7] <- NULL
}
}
}
if (Ncomp >= 7) {
if (!is.null(results7)) {
if (!is.null(results7$ET.AnnualAve)) {
comp[,8] <- results7$ET.AnnualAve
} else {
warning("No estimated annually averaged ET available for plotting from results8")
comp[,8] <- NULL
}
}
}
if (is.null(ylim[1])) {
ylim[1] <- 0
}
if (is.null(ylim[2])) {
ylim[2] <- max(results1$ET.AnnualAve)*1.5
}
plot(comp[,2]~floor(as.numeric(as.yearmon(time(results1$ET.Annual)))), type = "l", main = "Annually average ET", ylim = ylim, col = 2, xlab = "Month", ylab = "Annual averages of estimated daily ET mm/day", cex.lab = 1.3, cex.main = 2.1, par(mar = c(5.1,4.7,4.7,2.1)))
for (i in 1:(Ncomp-1)) {
lines(comp[,i+2]~floor(as.numeric(as.yearmon(time(results1$ET.Annual)))), type = "l", col = i+2)
}
legend("topright", inset = .01, c(paste(legtext[1:Ncomp])), cex = (1), col = c(seq(from=2,to=Ncomp+1)), lty = 1)
# Non-exceedance probability
par(ask=TRUE)
sortcomp <- matrix(nrow=length(results1$ET.AnnualAve),ncol=Ncomp+2)
sortcomp[,1] <- comp[,1]
for (j in 1:Ncomp) {
sortcomp[,j+1] <- sort(comp[,j+1])
}
sortcomp[,Ncomp+2] <- seq(from = 1/nrow(sortcomp), to = 1, by = 1/nrow(sortcomp))
plot(sortcomp[,2]~sortcomp[,Ncomp+2], type = "l", main = "Non-exceedance probability of the annual average ET", xlim = c(0,1), ylim = ylim, col = 2, xlab = "Non-exceedance probability", ylab = "Annual averages of estimated daily ET mm/day", cex.lab = 1.3, cex.main = 2.1, par(mar = c(5.1,4.7,4.7,2.1)))
for (k in 1:(Ncomp-1)) {
lines(sortcomp[,k+2]~sortcomp[,Ncomp+2], type = "o", pch = ".", col = k+2)
}
legend("topright", inset = .01, c(paste(legtext[1:Ncomp])), cex = (1), col = c(seq(from=2,to=Ncomp+1)), lty = 1)
# Box plot
par(ask=TRUE)
colnames(comp)[2:(Ncomp+1)] <- legtext
boxplot(comp[,2:(Ncomp+1)], col = c(seq(from=2,to=Ncomp+1)), boxwex = 0.5, range = 0, main = "Box plot of the annual average ET", xlab = "", ylab = "Distribution ofannually averaged PET in mm/day", show.names = FALSE, pars = list(ylim = ylim, cex.lab = 1.3, cex.main = 2.1, par(mar = c(5.1,4.7,4.7,2.1))))
legend("topright", inset = .01, c(paste(legtext[1:Ncomp])), cex = (1), fill = c(seq(from=2,to=Ncomp+1)), bty="o")
}
par(ask=FALSE)
}
#-------------------------------------------------------------------------------------
ETForcings <- function(data, results, forcing) {
par(ask=FALSE)
plot.new()
# define forcing names
Fnames <- list(v=c("Tmax","Tmin","u2","uz","Rs","n","Precip","Epan","RHmax","RHmin","Tdew"),
n=c("maximum temperature","minimum temperature", "average wind speed at 2m", "average wind speed", "daily solar radiation", "daily sunshine hours",
"precipitation", "Class-A pan evaporation", "maximum relative humidity", "minimum relative humidity", "average dew point temeprature"),
u=c("degree Celcius","degree Celcius", "m/s", "m/s", "MJ.per.sqm.per.day", "hours", "mm/day", "mm/day", "%", "%", "degree Celcius"))
# check forcing value
if (forcing %in% Fnames$v) {
par(mfrow=c(2,2))
# Aggregated PET vs forcing
# Plot daily relationship
if (!is.null(results$ET.Daily)) {
interval <- "daily"
if (!is.null(data[[forcing]])) {
plot(results$ET.Daily~data[[forcing]], main = results$ET_formulation, xlab = paste(interval, Fnames$n[Fnames$v == forcing], Fnames$u[Fnames$v == forcing]), ylab = list(c(results$ET_type, "mm/day")))
} else {
plot(1, type="n", axes=F, xlab="", ylab="")
text(1, paste("no", forcing, "data to plot"))
}
} else {
plot(1, type="n", axes=F, xlab="", ylab="")
text(1, "no daily ET estimation to plot")
}
# Plot monthly relationship
if (!is.null(results$ET.Monthly)) {
interval <- "monthly"
if (!is.null(data[[forcing]])) {
if (forcing != "Precip") {
Fmonthly <- aggregate(data[[forcing]], as.yearmon(data$Date.daily, "%m/%y"), mean)
} else {
Fmonthly <- aggregate(data[[forcing]], as.yearmon(data$Date.daily, "%m/%y"), sum)
}
plot(results$ET.Monthly~Fmonthly, main = results$ET_formulation, xlab = paste(interval, Fnames$n[Fnames$v == forcing], Fnames$u[Fnames$v == forcing]), ylab = list(c(results$ET_type, "mm/month")))
} else {
plot(1, type="n", axes=F, xlab="", ylab="")
text(1, paste("no", forcing, "data to plot"))
}
} else {
plot(1, type="n", axes=F, xlab="", ylab="")
text(1, "no monthly ET estimation to plot")
}
# Plot annual relationship
if (!is.null(results$ET.Annual)) {
interval <- "annual"
if (!is.null(data[[forcing]])) {
if (forcing != "Precip") {
Fannual <- aggregate(Fmonthly, floor(as.numeric(as.yearmon(data$Date.monthly, "%m/%y"))), mean)
} else {
Fannual <- aggregate(Fmonthly, floor(as.numeric(as.yearmon(data$Date.monthly, "%m/%y"))), sum)
}
plot(results$ET.Annual~Fannual, main = results$ET_formulation, xlab = paste(interval, Fnames$n[Fnames$v == forcing], Fnames$u[Fnames$v == forcing]), ylab = list(c(results$ET_type, "mm/year")))
} else {
plot(1, type="n", axes=F, xlab="", ylab="")
text(1, paste("no", forcing, "data to plot"))
}
} else {
plot(1, type="n", axes=F, xlab="", ylab="")
text(1, "no annual ET estimation to plot")
}
} else {
stop("Please select forcing from 'Tmax','Tmin','u2','uz','Rs','n','Precip','Epan','RHmax','RHmin','Tdew'") # forcing is not defined
}
}
#-------------------------------------------------------------------------------------
Try the Evapotranspiration package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.