Newer
Older
StatSeasAtlHurr <- 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
}