Newer
Older
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
81
82
83
84
85
86
87
88
StatSeasAtlTC <- function(atlano=NULL,tropano=NULL,hrvar='HR') {
#
# Verify that variables are either TC, HR or PDI.
# -----------------------------------------------
#
if (hrvar != "HR" && hrvar != "TC" && hrvar != "PDI") {
stop("Hurricane variable not recognized.")
}
#
# Verify that both Atl and Trop SSTA are present.
# -----------------------------------------------
#
if (is.null(atlano) ) {
stop("Atlantic SST missing")
}
if (is.null(tropano)) {
stop("Tropical SST missing")
}
#
# Verify that Atl and Trop SSTA are of the same dimensions.
# ---------------------------------------------------------
#
if(length(dim(atlano)) != length(dim(tropano)) ){
stop("Input arrays are of different dimensions.")
} else {
for (i in 1:length(dim(atlano))){
if (dim(atlano)[i] != dim(tropano)[i]) {
stop("Input arrays are of different sizes.")
}
}
}
#
# Get the values of the betas according to the hurricane activity measure we specified.
# -------------------------------------------------------------------------------------
#
if(hrvar=='HR'){
# beta's are derived from Villarini et al. (2012), Mon Wea Rev, 140, 44-65.
# beta's are for corrected hurricane data + ERSST with SBC criteria (table 2)
beta0=1.85
betaAtl=1.05
betaTrop=-1.17
} else if (hrvar=='TC'){
# beta's are from Villarini et al. (2010), Mon Wea Rev, 138, 2681-2705.
# beta's are for corrected TC data (lifetime >=48h) + ERSST (table 5)
beta0= 2.10
betaAtl= 1.02
betaTrop= -1.05
} else if (hrvar=='PDI'){
# beta's are from Villarini et al. (2012), J Clim, 25, 625-637.
# beta's are from ERSST, with SBC penalty criterion (table 1)
beta0=0.76
betaAtl=1.94
betaTrop=-1.78
}
#
# Create matrix of similar dimension as atlano for beta0.
# -------------------------------------------------------
#
intercept <- array(beta0, dim(atlano) )
#
# Compute statistical relationship b/w SSTAs and mean hurricane activity.
# -----------------------------------------------------------------------
#
atl <- betaAtl * atlano
trop <- betaTrop * tropano
#
temp <- intercept + atl + trop
#
statval <- list(mean=array(NaN,dim(atl)),var=array(NaN,dim(atl)) )
statval$mean[] <- vapply(X=temp, FUN=exp, numeric(1))
#
# Compute the variance of the distribution.
# TC and HR follow a Poisson distribution, so the variance is equal to the mean.
# PDI follows a gamma distribution, with sigma=-0.57. (variance= sigma^2 * mean^2).
# ----------------------------------------------------------------------------------
#
if(hrvar == "HR" && hrvar == "TC"){
statval$var <- statval$mean
} else {
sigma=-0.57
statval$var[] <- sigma^2 * vapply(X=statval$mean, FUN=function(x) x^2, numeric(1))
}
#
# Output
# ~~~~~~~~
#
statval
}