-
Notifications
You must be signed in to change notification settings - Fork 34
/
thornthwaite.R
executable file
·80 lines (73 loc) · 2.18 KB
/
thornthwaite.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
#' @title Computation of potential evapotranspiration.
#'
#'
#' @description See hargreaves
#'
#'
#' @details See hargreaves
#'
#'
#' @return A time series with the values of monthly potential or reference evapotranspiration, in mm.
#' If the input is a matrix or a multivariate time series each column will be treated as independent
#' data (e.g., diferent observatories), and the output will be a multivariate time series.
#'
#'
#' @rdname Potential-evapotranspiration
#'
#'
#' @importFrom stats ts cycle frequency start
#'
#'
#' @export
#'
thornthwaite <-
function(Tave, lat, na.rm=FALSE) {
if (anyNA(Tave) && na.rm==FALSE) {
stop('Data must not contain NAs')
}
if (is.ts(Tave) && frequency(Tave)!=12) {
stop('Data should be a monthly time series (frequency = 12)')
}
if (!is.null(ncol(Tave))) {
if (length(lat)!=ncol(Tave)) {
stop('Longitude of latitudes vector does not coincide with the number of columns in Tave')
}
}
if (!is.ts(Tave)) {
Tave <- ts(as.matrix(Tave),frequency=12)
} else {
Tave <- ts(as.matrix(Tave),frequency=frequency(Tave),start=start(Tave))
}
PE <- Tave*NA
n <- nrow(Tave)
m <- ncol(Tave)
c <- cycle(Tave)
# Monthly correction factor, depending on latitude (K)
# tangens of the average solar declination angle for each month of the year
tanLat <- tan(lat/57.2957795)
tanDelta <- c(-0.37012566,-0.23853358,-0.04679872,0.16321764,
0.32930908,0.40677729,0.3747741,0.239063,
0.04044485,-0.16905776,-0.33306377,-0.40743608)
tanLatM <-matrix(tanLat,nrow=12,ncol=length(tanLat),byrow=TRUE)*tanDelta
tanLatM[tanLatM< (-1)] <- -1
tanLatM[tanLatM>1] <- 1
omega <- acos(-tanLatM)
# montly average number of sun hours
N <- 24*omega/pi
#
days <- c(31,28,31,30,31,30,31,31,30,31,30,31)
K <- N/12 * days/30
# Annual temperature efficiency index (J)
T <- matrix(unlist(by(Tave,c,colMeans,na.rm=na.rm)),ncol=m,byrow=TRUE)
T[T<0] <- 0
J <- colSums((T/5)^1.514,na.rm=na.rm)
J <- matrix(J,n,m,TRUE)
# Empirical exponent (c)
J2 <- J*J; J3 <- J2*J
q <- 0.000000675*J3 - 0.0000771*J2 + 0.01792*J + 0.49239
# Potential evapotranspiration series (PE)
Tave[Tave<0] <- 0
PE <- K[c,] * 16 * (10*Tave/J)^q
colnames(PE) <- rep('PET_tho',m)
return(PE)
}