R version 2.15.1 (2012-06-22) – “Roasted Marshmallows”
Copyright © 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)
date()
## [1] "Sat Jul 28 20:04:03 2012"
Set your working directory.
setwd("~/Dropbox/Sensitivity")
Get the required packages.
require(ggplot2)
require(rgdal)
require(ismev)
require(ncdf)
library(mapproj)
library(mapdata)
library(maptools)
library(maps)
library(grid)
library(colorRamps)
library(RColorBrewer)
require(spdep)
require(psych)
require(boot)
Get additional functions for extreme value modeling.
source("selectData.R")
source("ismevplus.R")
Load the hourly interpolated best-track data. Subtract 60% of the forward speed from the best-track speed to get an estimate of rotational velocity. See ftp://texmex.mit.edu/ftp/pub/emanuel/PAPERS/hurr_risk.pdf.
load("best.use.RData")
begin = 1981
end = 2010
Wind.df = best.use[best.use$Yr >= begin & best.use$Yr <= end, ]
Wind.df = subset(Wind.df, Type == "*")
Wind.df$WmaxS = Wind.df$WmaxS * 0.5144
quantile(Wind.df$WmaxS, prob = seq(0.6, 0.8, 0.01)) #33 m/s is the 69th percentile
## 60% 61% 62% 63% 64% 65% 66% 67% 68% 69% 70% 71%
## 28.44 28.87 29.43 29.95 30.48 30.87 31.24 31.80 32.41 32.99 33.34 33.49
## 72% 73% 74% 75% 76% 77% 78% 79% 80%
## 33.76 34.29 34.86 35.53 36.01 36.48 37.21 37.90 38.50
Wind.df$maguv = Wind.df$maguv * 0.5144
Wind.df = subset(Wind.df, WmaxS >= 33)
Wind.df$WmaxS = Wind.df$WmaxS - Wind.df$maguv * 0.6
length(unique(Wind.df$Sid)) # 188 tropical cyclones (1981-2010), 268 (1967-2010)
## [1] 188
Plot a histogram of the wind speeds.
p = ggplot(Wind.df, aes(WmaxS)) + geom_histogram(fill = "gray", binwidth = 2) +
theme_bw()
p + xlab(expression(paste("Hurricane Intensity (rotational) [m ",
s^-1, "]"))) + ylab("Hours")
Next create a SST data frame from the netCDF using data from www.esrl.noaa.gov/psd/data/gridded/data.noaa.ersst.html
nc = open.ncdf("sst.mnmean.nc")
sstvar = nc$var$sst
ncdata = get.var.ncdf(nc, sstvar)
dim(ncdata)
## [1] 180 89 1901
object.size(ncdata)
## 243632368 bytes
Rework the three-dimensional array into a data frame and add longitudes and latitudes.
vals = lapply(nc$var$sst$dim, function(x) as.vector(x$vals))
vals[[1]] = (vals[[1]] + 180)%%360 - 180
vals[[2]] = rev(vals[[2]])
timedate = as.POSIXlt(86400 * vals[[3]], origin = ISOdatetime(1800,
1, 1, 0, 0, 0, tz = "GMT"), tz = "GMT")
timecolumn = paste("Y", 1900 + timedate$year, "M", formatC(as.integer(timedate$mo +
1), 1, flag = "0"), sep = "")
names(vals) = sapply(nc$var$sst$dim, "[", "name")
vals = vals[1:2]
ncdata1 = ncdata[, (dim(ncdata)[2]:1), ]
dims = dim(ncdata1)
dim(ncdata1) = c(dims[1] * dims[2], dims[3])
colnames(ncdata1) = timecolumn
ncdata1 = as.data.frame(ncdata1)
ncdataframe = cbind(expand.grid(vals), ncdata1)
misbyrow = apply(ncdataframe, 1, function(x) sum(is.na(x)))
ncdataframe = ncdataframe[misbyrow == 0, ]
Subset on longitude and latitude to create a domain bounded by the equator and 70N latitude and 100W and 10E longitude. Then compute Aug-Oct SST values on the two-degree grid. Also choose the set of years.
ncdataframe = subset(ncdataframe, lon <= 10 & lon >= -100 & lat >=
0 & lat <= 70)
ll = "+proj=longlat +datum=WGS84"
lcc = "+proj=lcc +lat_1=60 +lat_2=30 +lon_0=-60"
years = begin:end
sehur = paste("Y", years, sep = "")
sstAug = ncdataframe[paste(sehur, "M", formatC(8, 1, flag = "0"),
sep = "")]
sstSep = ncdataframe[paste(sehur, "M", formatC(9, 1, flag = "0"),
sep = "")]
sstOct = ncdataframe[paste(sehur, "M", formatC(10), sep = "")]
sstAugmean = rowMeans(sstAug)
sstSepmean = rowMeans(sstSep)
sstOctmean = rowMeans(sstOct)
ASOsst = cbind(sstAugmean, sstSepmean, sstOctmean)
avgsst = rowMeans(ASOsst)
avgsst = as.data.frame(avgsst)
avgsst$lon = ncdataframe$lon
avgsst$lat = ncdataframe$lat
coordinates(avgsst) = c("lon", "lat")
proj4string(avgsst) = CRS(ll)
SST.sdf = spTransform(avgsst, CRS(lcc))
Now return to the hurricane data and create a projected spatial points data frame from the tracks.
Wind.sdf = Wind.df
coordinates(Wind.sdf) = c("lon", "lat")
proj4string(Wind.sdf) = CRS(ll)
Wind.sdf = spTransform(Wind.sdf, CRS(lcc))
Create a hexagon tessellation of the basin.
hpt = spsample(Wind.sdf, type = "hexagonal", n = 120, bb = bbox(Wind.sdf) *
1.08, offset = c(1, 1))
hpg = HexPoints2SpatialPolygons(hpt)
np = length(hpg@polygons)
area = hpg@polygons[[1]]@area * 10^-6 # area has units of square km
area # 428,496 km^2
## [1] 428496
Function for returning the max cyclone intensity per grid.
get.max.id <- function(x, in.order = FALSE, idfield = "Sid", maxfield = "WmaxS") {
use.max <- x[sapply(split((1:nrow(x)), x[, idfield]), function(i, x) {
xi <- x[i]
ii <- i[max(xi) == xi]
ii[1]
}, x = x[, maxfield]), ]
if (in.order)
use.max <- use.max[rev(order(use.max[, maxfield])), ]
use.max
}
Overlay the hexagon tessellation of the basin on the track points and get the per hurricane maximum.
Wind.hexid = over(x = Wind.sdf, y = hpg)
Wind.sdf$whexid = Wind.hexid
Wind.split = split(Wind.sdf@data, Wind.hexid)
maximum = NULL
for (i in 1:length(Wind.split)) {
maximum[[i]] = get.max.id(Wind.split[[i]])
}
less = which(sapply(maximum, function(x) dim(x)[1] < 15)) # free parameter
maximum = maximum[-less]
id = as.integer(lapply(maximum, function(x) unique(x$whexid)))
hexid = paste("ID", id, sep = "")
dat = over(x = hpg[hexid], y = SST.sdf, fn = mean)
count = unlist(lapply(maximum, function(x) dim(x)[1]))
dat$count = count
hspdf = SpatialPolygonsDataFrame(hpg[hexid], dat, match.ID = TRUE)
Determine the correlation between average SST and hurricane count across the domain. Determine the distribution of hurricanes per grid.
cor(dat$avgsst, dat$count) # -0.194
## [1] -0.194
table(count)
## count
## 15 17 18 20 21 22 24 25 26 27 28 32
## 2 1 3 4 2 2 4 1 1 1 1 2
Determine the per hexagon maximum wind speed and add this value to the spatial polygon data frame.
obsmax = numeric()
for (i in 1:length(maximum)) {
obsmax[i] = max(maximum[[i]]$WmaxS)
}
hspdf$obsmax = obsmax
Obtain map borders.
cl = map("worldHires", xlim = c(-120, 20), ylim = c(-10, 70), plot = TRUE,
resolution = 0.7)
clp = map2SpatialLines(cl, proj4string = CRS(ll))
clp = spTransform(clp, CRS(lcc))
l1 = list("sp.lines", clp, col = "gray")
Plot the hexagons with color scale corresponding to the number of hurricanes and the numbers corresponding to the grid IDs.
centers = coordinates(hspdf)
l2 = list("sp.text", centers, id, col = "red", cex = 1.1)
cr = brewer.pal(4, "Blues")
cr = brewer.pal(7, "Blues")
rng = seq(15, 35, 5)
rng = seq(15, 50, 5)
p0 = spplot(hspdf, "count", col = "white", col.regions = cr, sp.layout = list(l1,
l2), at = rng, colorkey = list(space = "bottom", labels = paste(rng)), sub = expression("Number of hurricanes (1981-2010)"))
p0
Create plots showing the frequency and maximum intensity in each grid.
l2 = list("sp.text", centers, id, col = "black", cex = 1.1)
cr = brewer.pal(4, "Blues")
rng = seq(15, 35, 5)
p1 = spplot(hspdf, "count", col = "white", col.regions = cr, sp.layout = list(l1,
l2), at = rng, colorkey = list(space = "bottom", labels = paste(rng)), sub = expression("Number of hurricanes (1981-2010)"))
cr = brewer.pal(4, "Reds")
rng = seq(40, 80, 10)
p2 = spplot(hspdf, "obsmax", col = "white", col.regions = cr, sp.layout = list(l1,
l2), at = rng, colorkey = list(space = "bottom", labels = paste(rng)), sub = expression(paste("Observed highest intensity [m ",
s^-1, "]")))
p1 = update(p1, main = textGrob("a", x = unit(0.05, "npc")), par.settings = list(fontsize = list(text = 15)))
p2 = update(p2, main = textGrob("b", x = unit(0.05, "npc")), par.settings = list(fontsize = list(text = 15)))
p1
p2
Choose sample grids. Choose one with fewer but stronger and another with more but weaker.
iC = which(id == 18)
iD = which(id == 40)
iE = which(id == 37)
WC = maximum[[iC]]$WmaxS
WD = maximum[[iD]]$WmaxS
WE = maximum[[iE]]$WmaxS
Plot a histogram of the rotational wind speeds from hurricanes in grid E.
WE.df = data.frame(WE = WE)
p21 = ggplot(WE.df, aes(WE)) + geom_histogram(binwidth = 10) + geom_rug()
p21 = p21 + xlab("Wind speed (rotational) [m/s]") + ylab("Frequency") +
theme_bw()
p21
Plot histograms of the rotational wind speeds from hurricanes in grids C and D. Here you use trellis graphics rather than ggplot.
p3 = histogram(WC, breaks = seq(20, 80, 5), ylim = c(-1.5, 35), col = "black",
border = "white", aspect = 0.8, panel = function(x, ...) {
panel.histogram(x, ...)
panel.rug(x, col = "gray")
}, xlab = expression(paste("Wind speed [m ", s^-1, "]")), ylab = "Percent of total")
p4 = histogram(WD, breaks = seq(20, 80, 5), ylim = c(-1.5, 35), col = "black",
border = "white", aspect = 0.8, panel = function(x, ...) {
panel.histogram(x, ...)
panel.rug(x, col = "gray")
}, xlab = expression(paste("Wind speed [m ", s^-1, "]")), ylab = "Percent of total")
p3 = update(p3, main = textGrob("c", x = unit(0.05, "npc")), par.settings = list(fontsize = list(text = 15)))
p4 = update(p4, main = textGrob("d", x = unit(0.05, "npc")), par.settings = list(fontsize = list(text = 15)))
p3
p4
Estimate the empirical return periods for hurricane wind speeds in grids C, D, and E.
m = length(begin:end)
WCs = sort(WC, decreasing = TRUE)
mC = rev(rank(WCs))
rpC = round((m + 1)/mC, 0)
WDs = sort(WD, decreasing = TRUE)
mD = rev(rank(WDs))
rpD = round((m + 1)/mD, 0)
WEs = sort(WE, decreasing = TRUE)
mE = rev(rank(WEs))
rpE = round((m + 1)/mE, 0)
Compute values for the GPD parameter using the method of maximum likelihood and then use a simple functional relationship to link the parameters to the limiting intensity. Estimate the model parameters for each hexagon including the limiting intensity.
threshC = quantile(WC, prob = 0.25)
rateC = length(WC)/m
modelC = exceed.pp.plot(WC, npy = rateC, threshold = threshC, yunit = "m/s",
addtext = "", plotrl = TRUE)
threshD = quantile(WD, prob = 0.25)
rateD = length(WD)/m
modelD = exceed.pp.plot(WD, npy = rateD, threshold = threshD, yunit = "m/s",
addtext = "", plotrl = TRUE)
threshE = quantile(WE, prob = 0.25)
rateE = length(WE)/m
modelE = exceed.pp.plot(WE, npy = rateE, threshold = threshE, yunit = "m/s",
addtext = "", plotrl = TRUE)
sC = modelC$gpd.coef[[2]]
xC = modelC$gpd.coef[[3]]
theoC = threshC - sC/xC
sD = modelD$gpd.coef[[2]]
xD = modelD$gpd.coef[[3]]
theoD = threshD - sD/xD
sE = modelE$gpd.coef[[2]]
xE = modelE$gpd.coef[[3]]
theoE = threshE - sE/xE
Plot the model and the model parameters from grid E.
rpK = c(2, 5, 10, 11, 14, 19)
data.df = modelE$return.level.curves[rpK, ]
data2.df = subset(modelE$empirical.plot, rp > 2)
data2.df$lower = rep(1, length(data2.df$rl))
data2.df$upper = rep(1, length(data2.df$rl))
p22 = ggplot(data = data.df, aes(x = rp, y = q, ymin = lower, ymax = upper)) +
geom_pointrange() + scale_x_log10(breaks = c(2, 5, 10, 20, 50, 100), labels = c(2,
5, 10, 20, 50, 100)) + scale_y_continuous(limits = c(24, 80)) + xlab("Return period [yr]") +
ylab("Wind speed [m/s]")
p22 = p22 + geom_point(aes(x = rp, y = rl), data = data2.df, color = "red") +
theme_bw() + geom_hline(aes(yintercept = theoE), color = "red") + geom_hline(aes(yintercept = threshE),
color = "green")
p22
Create a two panel plot using plot objects p21-p22.
grid.newpage()
pushViewport(viewport(layout = grid.layout(1, 2)))
vplayout <- function(x, y) viewport(layout.pos.row = x, layout.pos.col = y)
print(p21 + theme_bw() + opts(axis.title.x = theme_text(vjust = 0),
axis.title.y = theme_text(size = 10, angle = 90), title = "a", plot.title = theme_text(hjust = 0,
size = 11), aspect.ratio = 1), vp = vplayout(1, 1))
print(p22 + theme_bw() + opts(axis.title.x = theme_text(vjust = 0),
axis.title.y = theme_text(size = 10, angle = 90), title = "b", plot.title = theme_text(hjust = 0,
size = 11), aspect.ratio = 1), vp = vplayout(1, 2))
Create a multiple panel plot using objects p1-p6.
p5 = xyplot(q ~ rp, data = modelC$return.level.curves, aspect = 0.8,
type = "l", lty = c(1, 3, 3), lwd = c(2, 1, 1), col = 1, scales = list(x = list(log = 10,
at = c(1, 10, 100, 1000))), ylim = c(25, 85), panel = function(...) {
panel.abline(h = threshC, lty = 2)
panel.abline(h = theoC, col = 2)
panel.xyplot(x = log10(modelC$empirical.plot$rp), y = modelC$empirical.plot$rl,
pch = 16, col = "grey")
panel.xyplot(x = log10(modelC$return.level.curves$rp), y = modelC$return.level.curves$lower,
type = "l", col = "lightgray")
panel.xyplot(x = log10(modelC$return.level.curves$rp), y = modelC$return.level.curves$upper,
type = "l", col = "lightgray")
panel.xyplot(...)
panel.text(x = log10(1000), y = threshC + 3, labels = expression(paste("u",
{
}[c])), cex = 1.1)
panel.text(x = log10(1), y = theoC + 4, labels = expression(paste("LI",
{
}[c])), cex = 1)
}, xlab = "Return period [yr]", ylab = expression(paste("Return level [m ",
s^-1, "]")))
p6 = xyplot(q ~ rp, data = modelD$return.level.curves, aspect = 0.8,
type = "l", lty = c(1, 3, 3), lwd = c(2, 1, 1), col = 1, scales = list(x = list(log = 10,
at = c(1, 10, 100, 1000))), ylim = c(25, 85), panel = function(...) {
panel.abline(h = threshD, lty = 2)
panel.abline(h = theoD, col = 2)
panel.xyplot(x = log10(modelD$empirical.plot$rp), y = modelD$empirical.plot$rl,
pch = 16, col = "grey")
panel.xyplot(x = log10(modelD$return.level.curves$rp), y = modelD$return.level.curves$lower,
type = "l", col = "lightgray")
panel.xyplot(x = log10(modelD$return.level.curves$rp), y = modelD$return.level.curves$upper,
type = "l", col = "lightgray")
panel.xyplot(...)
panel.text(x = log10(1000), y = threshD + 3, labels = expression(paste("u",
{
}[d])), cex = 1.1)
panel.text(x = log10(1), y = theoD + 4, labels = expression(paste("LI",
{
}[d])), cex = 1)
}, xlab = "Return period [yr]", ylab = expression(paste("Return level [m ",
s^-1, "]")))
p5 = update(p5, main = textGrob("e", x = unit(0.05, "npc")), par.settings = list(fontsize = list(text = 15)))
p6 = update(p6, main = textGrob("f", x = unit(0.05, "npc")), par.settings = list(fontsize = list(text = 15)))
plot(p1, split = c(1, 1, 2, 3), more = TRUE)
plot(p2, split = c(2, 1, 2, 3), more = TRUE)
plot(p3, split = c(1, 2, 2, 3), more = TRUE)
plot(p4, split = c(2, 2, 2, 3), more = TRUE)
plot(p5, split = c(1, 3, 2, 3), more = TRUE)
plot(p6, split = c(2, 3, 2, 3), more = FALSE)
Estimate the statistical model parameters using per hexagon per hurricane maximum wind speed in each grid.
m = length(begin:end)
lambda = numeric()
sigma = numeric()
xi = numeric()
thresh = numeric()
rate = numeric()
theo = numeric()
for (i in 1:length(maximum)) {
W = maximum[[i]]$WmaxS
thresh[i] = quantile(W, prob = 0.25)
rate[i] = length(W)/m
model = exceed.pp.plot(W, npy = rate[i], threshold = thresh[i], yunit = "m/s",
addtext = "", plotrl = FALSE)
lambda[i] = model$gpd.coef[1]
sigma[i] = model$gpd.coef[2]
xi[i] = model$gpd.coef[3]
theo[i] = thresh[i] - sigma[i]/xi[i]
}
hspdf$thresh = thresh
hspdf$rate = rate
hspdf$lambda = lambda
hspdf$sigma = sigma
hspdf$xi = xi
hspdf$theo = theo
Note that with \( \xi \) less than 0 the first derivative of the return level curve decreases with longer return periods. With \( \xi \) greater than 0 the first derivative of the return level increases with longer return periods. With \( \xi \) less than -1 the standard error goes to zero and the predicted maximum is the highest data value in the set. This will tend to under estimate the limiting intensity. It is biased, but probably not by much.
which(xi < -1)
## [1] 1 2 3 4 5 6 7 8 10 13 14 20
Map the model parameters.
clsThr = brewer.pal(4, "Oranges")
clsRate = brewer.pal(6, "Blues")
clsSigma = brewer.pal(5, "Greens")
clsXi = brewer.pal(8, "Purples")
clsLi = brewer.pal(8, "Reds")
rng = seq(25, 45, 5)
p7 = spplot(hspdf["thresh"], sp.layout = list(l1), col.regions = clsThr,
sub = expression(paste("u [m ", s^-1, "]")), colorkey = list(space = "bottom"),
col = "white", at = rng, labels = paste(rng))
rng = seq(0.3, 0.9, 0.1)
p8 = spplot(hspdf["lambda"], sp.layout = list(l1), col.regions = clsRate,
sub = expression(paste(lambda[u], " [yr", {
}^-1, "]")), colorkey = list(space = "bottom"), col = "white", at = rng,
labels = paste(rng))
rng = seq(0, 50, 10)
p8 = spplot(hspdf["sigma"], sp.layout = list(l1), col.regions = clsSigma,
sub = expression(paste(sigma, " [m ", s^-1, "]")), colorkey = list(space = "bottom"),
col = "white", at = rng, labels = paste(rng))
rng = seq(-1.6, 0, 0.2)
p9 = spplot(hspdf["xi"], sp.layout = list(l1), col.regions = clsXi,
sub = expression(paste(xi)), colorkey = list(space = "bottom"), col = "white",
at = rng, labels = paste(rng))
rng = seq(40, 80, 5)
p10 = spplot(hspdf["theo"], sp.layout = list(l1), col.regions = clsLi,
sub = expression(paste("Limiting intensity (LI) [m ", s^-1, "]")), colorkey = list(space = "bottom"),
col = "white", at = rng, labels = paste(rng))
p7 = update(p7, main = textGrob("b", x = unit(0.05, "npc")), par.settings = list(fontsize = list(text = 15)))
p8 = update(p8, main = textGrob("c", x = unit(0.05, "npc")), par.settings = list(fontsize = list(text = 15)))
p9 = update(p9, main = textGrob("d", x = unit(0.05, "npc")), par.settings = list(fontsize = list(text = 15)))
p10 = update(p10, main = textGrob("a", x = unit(0.05, "npc")), par.settings = list(fontsize = list(text = 15)))
plot(p10, split = c(1, 1, 2, 2), more = TRUE)
plot(p7, split = c(2, 1, 2, 2), more = TRUE)
plot(p8, split = c(1, 2, 2, 2), more = TRUE)
plot(p9, split = c(2, 2, 2, 2), more = FALSE)
Map the SST using the same grids.
clsSST = brewer.pal(5, "YlOrRd")
rng = seq(21, 31, 2)
pSST = spplot(hspdf, "avgsst", col = "white", col.regions = clsSST,
sp.layout = list(l1), at = rng, colorkey = list(space = "bottom", labels = paste(rng)),
sub = expression(paste("Sea surface temperature [", degree, "C]")))
pSST
fs1 = numeric()
for (i in 21:24) fs1 = c(fs1, maximum[[i]]$maguv)
fs2 = numeric()
for (i in 1:20) {
fs2 = c(fs2, maximum[[i]]$maguv)
fs3 = maximum[[i]]$maguv
print(t.test(fs1, fs3)$p.value)
}
## [1] 1.619e-10
## [1] 9.547e-11
## [1] 3.405e-08
## [1] 2.605e-09
## [1] 1.413e-09
## [1] 7.875e-10
## [1] 1.291e-11
## [1] 2.656e-11
## [1] 7.036e-06
## [1] 3.63e-11
## [1] 1.224e-05
## [1] 6.69e-09
## [1] 7.231e-05
## [1] 2.076e-08
## [1] 1.07e-11
## [1] 0.00159
## [1] 0.001843
## [1] 0.01091
## [1] 0.0009052
## [1] 0.01135
mean(fs1)
## [1] 11.2
mean(fs2)
## [1] 6.366
The translation speed of hurricanes moving through the four northern grids where SST is cooler than 25C is significantly faster (p-value < .05) than the speed of hurricanes elsewhere. Thus these grids are removed from further analysis.
Extract a data frame of LI and SST and compute the sensitivity. Use weights=count for weighted regression.
dat = hspdf@data
dat$northing = coordinates(hspdf)[, 2]
summary(lm(theo ~ avgsst, data = dat[dat$avgsst >= 25, ]))
##
## Call:
## lm(formula = theo ~ avgsst, data = dat[dat$avgsst >= 25, ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.588 -2.818 -0.707 3.428 7.177
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -165.50 33.11 -5.00 9.3e-05 ***
## avgsst 8.24 1.19 6.94 1.8e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.82 on 18 degrees of freedom
## Multiple R-squared: 0.728, Adjusted R-squared: 0.713
## F-statistic: 48.1 on 1 and 18 DF, p-value: 1.76e-06
##
summary(lm(theo ~ avgsst, weights = count, data = dat[dat$avgsst >=
25, ]))
##
## Call:
## lm(formula = theo ~ avgsst, data = dat[dat$avgsst >= 25, ], weights = count)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -47.5 -14.9 -2.5 15.8 33.7
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -155.54 33.03 -4.71 0.00018 ***
## avgsst 7.89 1.19 6.64 3.1e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 22.8 on 18 degrees of freedom
## Multiple R-squared: 0.71, Adjusted R-squared: 0.694
## F-statistic: 44.1 on 1 and 18 DF, p-value: 3.11e-06
##
summary(lm(theo ~ avgsst, weights = count, data = dat[dat$avgsst >=
23, ]))
##
## Call:
## lm(formula = theo ~ avgsst, data = dat[dat$avgsst >= 23, ], weights = count)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -49.16 -31.74 -1.41 30.42 86.54
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -14.45 28.75 -0.50 0.620
## avgsst 2.85 1.05 2.71 0.013 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 37.9 on 21 degrees of freedom
## Multiple R-squared: 0.259, Adjusted R-squared: 0.223
## F-statistic: 7.33 on 1 and 21 DF, p-value: 0.0132
##
Plot the parameters as a function of SST.
xlabel = expression(paste("Sea surface temperature [", degree, "C]"))
p11 = ggplot(dat[dat$avgsst >= 25, ], aes(avgsst, thresh)) + geom_point() +
geom_smooth(method = "lm", aes(weight = count), color = "red") + theme_bw() +
xlab(xlabel) + ylab(expression(paste("u [m ", s^-1, "]")))
p12 = ggplot(dat[dat$avgsst >= 25, ], aes(avgsst, sigma)) + geom_point() +
geom_smooth(method = "lm", aes(weight = count), color = "red") + theme_bw() +
xlab(xlabel) + ylab(expression(paste(sigma, " [m ", s^-1, "]")))
p13 = ggplot(dat[dat$avgsst >= 25, ], aes(avgsst, xi)) + geom_point() +
geom_smooth(method = "lm", aes(weight = count), color = "red") + theme_bw() +
xlab(xlabel) + ylab(expression(paste(xi)))
p14 = ggplot(dat[dat$avgsst >= 25, ], aes(avgsst, theo)) + geom_point() +
geom_smooth(method = "lm", aes(weight = count), color = "red") + theme_bw() +
xlab(xlabel) + ylab(expression(paste("Limiting intensity [m ", s^-1, "]")))
p11 = p11 + opts(title = "b", plot.title = theme_text(hjust = 0,
size = 15))
p12 = p12 + opts(title = "c", plot.title = theme_text(hjust = 0,
size = 15))
p13 = p13 + opts(title = "d", plot.title = theme_text(hjust = 0,
size = 15))
p14 = p14 + opts(title = "a", plot.title = theme_text(hjust = 0,
size = 15))
grid.newpage()
pushViewport(viewport(layout = grid.layout(2, 2)))
vplayout = function(x, y) viewport(layout.pos.row = x, layout.pos.col = y)
print(p14, vp = vplayout(1, 1))
print(p11, vp = vplayout(1, 2))
print(p12, vp = vplayout(2, 1))
print(p13, vp = vplayout(2, 2))
Check for spatial autocorrelation in the model residuals using Moran's I after removing grids with SST < 25C.
hspdf2 = hspdf[hspdf$avgsst >= 25, ]
hexnb2 = poly2nb(hspdf2)
wts2 = nb2listw(hexnb2, style = "W")
stvty.lm = lm(theo ~ avgsst, weights = count, data = hspdf2)
lm.morantest(stvty.lm, wts2)
##
## Global Moran's I for regression residuals
##
## data:
## model: lm(formula = theo ~ avgsst, data = hspdf2, weights = count)
## weights: wts2
##
## Moran I statistic standard deviate = 0.2005, p-value = 0.4206
## alternative hypothesis: greater
## sample estimates:
## Observed Moran's I Expectation Variance
## -0.05884 -0.09324 0.02946
##
The p-value indicates no significant spatial autocorrelation so it is reasonable to use a non-spatial regression.
Create a bootstrap function to assess uncertainty on the estimates of limiting intensity.
boot.theo = function(data, i, nyear, prob = 0.25) {
x = data[i]
thresh = quantile(x, prob = prob)
rate = length(x)/nyear
model = exceed.pp.plot(x, npy = rate, threshold = thresh, yunit = "m/s",
addtext = "", plotrl = FALSE)
xi = model$gpd.coef[[3]]
sigma = model$gpd.coef[[2]]
if (xi >= 0)
return(Inf) else return(thresh - sigma/xi)
}
Use the function above to compute uncertainty on the limiting intensity. By adjusting minHur, maxGrid, minSST you can test the sensitivity of the results to changes in these settings.
m = length(begin:end)
minHur = 15
maxGrid = 120
minSST = 25
hpt = spsample(Wind.sdf, type = "hexagonal", n = maxGrid, bb = bbox(Wind.sdf) *
1.08, offset = c(1, 1))
hpg = HexPoints2SpatialPolygons(hpt)
np = length(hpg@polygons)
area = hpg@polygons[[1]]@area * 10^-6
Wind.hexid = over(x = Wind.sdf, y = hpg)
Wind.sdf$whexid = Wind.hexid
Wind.split = split(Wind.sdf@data, Wind.hexid)
maximum = NULL
for (i in 1:length(Wind.split)) {
maximum[[i]] = get.max.id(Wind.split[[i]])
}
less = which(sapply(maximum, function(x) dim(x)[1] < minHur))
maximum = maximum[-less]
id = as.integer(lapply(maximum, function(x) unique(x$whexid)))
hexid = paste("ID", id, sep = "")
dat = over(x = hpg[hexid], y = SST.sdf, fn = mean)
count = unlist(lapply(maximum, function(x) dim(x)[1]))
dat$count = count
hspdf = SpatialPolygonsDataFrame(hpg[hexid], dat, match.ID = TRUE)
obsmax = numeric()
for (i in 1:length(maximum)) {
obsmax[i] = max(maximum[[i]]$WmaxS)
}
hspdf$obsmax = obsmax
theo = numeric()
theoL = numeric()
theoU = numeric()
for (i in 1:length(maximum)) {
W = maximum[[i]]$WmaxS
boot.out = boot(W, boot.theo, R = 100, nyear = m)
boot.data = boot.ci(boot.out, conf = 0.8, type = "basic") #.80 gives .10 on upper end
theo[i] = boot.data$t0
theoL[i] = max(max(W), boot.data$basic[4])
theoU[i] = boot.data$basic[5]
}
hspdf$theo = theo
hspdf$theoL = theoL
hspdf$theoU = theoU
hspdf2 = hspdf[hspdf$avgsst >= minSST, ]
hexnb2 = poly2nb(hspdf2)
wts2 = nb2listw(hexnb2, style = "W", zero.policy = TRUE)
stvty.lm = lm(theo ~ avgsst, weights = count, data = hspdf2)
area
## [1] 428496
dim(hspdf)[1]
## [1] 24
dim(hspdf2)[1]
## [1] 20
summary(stvty.lm)$coefficients[2, ]
## Estimate Std. Error t value Pr(>|t|)
## 7.891e+00 1.188e+00 6.642e+00 3.108e-06
Plot the trend line (slope is the sensitivity) with uncertainty lines on the LI estimates.
xlabel = expression(paste("Sea surface temperature [", degree, "C]"))
p = ggplot(hspdf2@data, aes(avgsst, theo)) + geom_point() + theme_bw()
limits = aes(ymax = hspdf2$theoU, ymin = hspdf2$theoL)
p = p + geom_linerange(limits)
p = p + geom_smooth(method = "lm", aes(weight = count))
p + xlab(xlabel) + ylab(expression(paste("Limiting intensity [m ",
s^-1, "]")))
Save the results.
OBS = hspdf2
HiRAM begin and end years.
begin1 = 1981
end1 = 2009
Load the data from gfdl_tracks.Rnw in GCMTrackData folder. Three runs are available r1-r3. Subset by begin and end dates and remove wind speeds less than 10 m/s. Plot a histogram of the wind speeds.
load("gfdlTracks_r1.df.RData")
natlwind.df = natltracks.df1[natltracks.df1$Yr >= begin1 & natltracks.df1$Yr <=
end1, ]
natlwind.df = subset(natlwind.df, WmaxS >= 10)
quantile(natlwind.df$WmaxS, prob = 0.69)
## 69%
## 27.34
p = ggplot(natlwind.df, aes(WmaxS)) + geom_histogram(fill = "gray",
binwidth = 2) + theme_bw()
p + xlab(expression(paste("Hurricane Intensity [m ", s^-1, "]"))) +
ylab("Hours")
Subset on the threshold wind speed to match observed hurricane intensity. Hurricane intensity corresponds to the 69th percentile of all observed TC wind speeds.
qthr = quantile(natlwind.df$WmaxS, prob = 0.69)
natlwind.df = subset(natlwind.df, WmaxS >= qthr)
natlwind.df$WmaxS = natlwind.df$WmaxS - natlwind.df$maguv * 0.6
length(unique(natlwind.df$Sid)) #number of TCs
## [1] 234
Compare rotational wind speeds.
mean(natlwind.df$WmaxS) # 25.37 (r1)
## [1] 25.37
mean(Wind.df$WmaxS) # 41.05
## [1] 41.05
skew(natlwind.df$WmaxS) # -0.240 (r1)
## [1] -0.2398
skew(Wind.df$WmaxS) # 0.812
## [1] 0.8124
plotdat = data.frame(WmaxS = c(Wind.df$WmaxS, natlwind.df$WmaxS),
Type = c(rep("Best-track", length(Wind.df$WmaxS)), rep("HiRAM (r1)", length(natlwind.df$WmaxS))))
p = ggplot(plotdat, aes(WmaxS)) + geom_histogram(fill = "gray", binwidth = 2) +
facet_grid(. ~ Type)
p + theme_bw() + xlab(expression(paste("Rotational speed [m ", s^-1,
"]"))) + ylab("Hours")
Compare translational speeds.
mean(natlwind.df$maguv) # 11.76 (r1)
## [1] 11.76
mean(Wind.df$maguv) # 6.168306
## [1] 6.168
plotdat = data.frame(maguv = c(Wind.df$maguv, natlwind.df$maguv),
Type = c(rep("Best-track", length(Wind.df$maguv)), rep("HiRAM (r1)", length(natlwind.df$maguv))))
p = ggplot(plotdat, aes(maguv)) + geom_histogram(fill = "gray", binwidth = 2) +
facet_grid(. ~ Type)
p + theme_bw() + xlab(expression(paste("Translational speed [m ",
s^-1, "]"))) + ylab("Hours")
Create a projected spatial points data frame from the tracks.
natlwind.sdf = natlwind.df
coordinates(natlwind.sdf) = c("lon", "lat")
proj4string(natlwind.sdf) = CRS(ll)
natlwind.sdf = spTransform(natlwind.sdf, CRS(lcc))
Create a hexagon tessellation of the basin.
GCMhpt = spsample(natlwind.sdf, type = "hexagonal", n = 120, bb = bbox(natlwind.sdf) *
1.08, offset = c(1, 1))
GCMhpg = HexPoints2SpatialPolygons(GCMhpt)
GCMnp = length(GCMhpg@polygons)
GCMarea = GCMhpg@polygons[[1]]@area * 10^-6 # area has units of square km, compare with area of obs
GCMarea # 477929.5 km^2 (r1) 510967.2 km^2 (r2) 516044 km^2 (r3)
## [1] 477929
Overlay the hexagon tessellation of the basin on the track points and get the per hurricane maximum.
GCMWind.hexid = over(x = natlwind.sdf, y = GCMhpg)
natlwind.sdf$whexid = GCMWind.hexid
GCMWind.split = split(natlwind.sdf@data, GCMWind.hexid)
GCMmaximum = NULL
for (i in 1:length(GCMWind.split)) {
GCMmaximum[[i]] = get.max.id(GCMWind.split[[i]])
}
GCMless = which(sapply(GCMmaximum, function(x) dim(x)[1] < 15))
GCMmaximum = GCMmaximum[-GCMless]
GCMid = as.integer(lapply(GCMmaximum, function(x) unique(x$whexid)))
GCMhexid = paste("ID", GCMid, sep = "")
GCMdat = over(x = GCMhpg[GCMhexid], y = SST.sdf, fn = mean)
GCMcount = unlist(lapply(GCMmaximum, function(x) dim(x)[1]))
GCMdat$count = GCMcount
GCMhspdf = SpatialPolygonsDataFrame(GCMhpg[GCMhexid], GCMdat, match.ID = TRUE)
GCMobsmax = numeric()
for (i in 1:length(GCMmaximum)) {
GCMobsmax[i] = max(GCMmaximum[[i]]$WmaxS)
}
GCMhspdf$obsmax = GCMobsmax
Determine the correlation between average SST and hurricane counts across the domain. Determine the distribution of hurricanes per grid.
cor(GCMdat$avgsst, GCMdat$count) # -0.157 (r1)
## [1] -0.1574
table(GCMcount)
## GCMcount
## 16 19 20 21 24 25 28 29 30 32 33 34 36 41 42 46
## 5 1 1 2 1 2 3 1 1 1 1 1 2 2 2 1
Estimate the statistical model parameters using per hexagon per hurricane maximum wind speed in each grid.
GCMm = length(begin1:end1)
GCMlambda = numeric()
GCMsigma = numeric()
GCMxi = numeric()
GCMthresh = numeric()
GCMrate = numeric()
for (i in 1:length(GCMmaximum)) {
GCMW = GCMmaximum[[i]]$WmaxS
GCMthresh[i] = quantile(GCMW, prob = 0.25)
GCMrate[i] = length(GCMW)/GCMm
GCMmodel = exceed.pp.plot(GCMW, npy = GCMrate[i], threshold = GCMthresh[i],
yunit = "m/s", addtext = "", plotrl = FALSE)
GCMlambda[i] = GCMmodel$gpd.coef[1]
GCMsigma[i] = GCMmodel$gpd.coef[2]
GCMxi[i] = GCMmodel$gpd.coef[3]
}
GCMhspdf$thresh = GCMthresh
GCMhspdf$rate = GCMrate
GCMhspdf$lambda = GCMlambda
GCMhspdf$sigma = GCMsigma
GCMhspdf$xi = GCMxi
GCMhspdf$theo = GCMthresh - GCMsigma/GCMxi
Map the model parameters.
clsThr = brewer.pal(7, "Oranges")
clsRate = brewer.pal(5, "Blues")
clsSigma = brewer.pal(4, "Greens")
clsXi = brewer.pal(7, "Purples")
clsLi = brewer.pal(8, "Reds")
l1 = list("sp.lines", clp, col = "gray")
rng = seq(14, 28, 2)
p17 = spplot(GCMhspdf["thresh"], sp.layout = list(l1), col.regions = clsThr,
sub = expression(paste("u [m ", s^-1, "]")), colorkey = list(space = "bottom"),
col = "white", at = rng, labels = paste(rng))
rng = seq(0.3, 1.3, 0.2)
p18 = spplot(GCMhspdf["lambda"], sp.layout = list(l1), col.regions = clsRate,
sub = expression(paste(lambda[u], " [yr", {
}^-1, "]")), colorkey = list(space = "bottom"), col = "white", at = rng,
labels = paste(rng))
rng = seq(5, 25, 5)
p18 = spplot(GCMhspdf["sigma"], sp.layout = list(l1), col.regions = clsSigma,
sub = expression(paste(sigma, " [m ", s^-1, "]")), colorkey = list(space = "bottom"),
col = "white", at = rng, labels = paste(rng))
rng = seq(-1.4, -0.2, 0.2)
p19 = spplot(GCMhspdf["xi"], sp.layout = list(l1), col.regions = clsXi,
sub = expression(paste(xi)), colorkey = list(space = "bottom"), col = "white",
at = rng, labels = paste(rng))
rng = seq(27, 47, 4)
p20 = spplot(GCMhspdf["theo"], sp.layout = list(l1), col.regions = clsLi,
sub = expression(paste("Limiting intensity [m ", s^-1, "]")), colorkey = list(space = "bottom"),
col = "white", at = rng, labels = paste(rng))
p17 = update(p17, main = textGrob("b", x = unit(0.05, "npc")), par.settings = list(fontsize = list(text = 15)))
p18 = update(p18, main = textGrob("c", x = unit(0.05, "npc")), par.settings = list(fontsize = list(text = 15)))
p19 = update(p19, main = textGrob("d", x = unit(0.05, "npc")), par.settings = list(fontsize = list(text = 15)))
p20 = update(p20, main = textGrob("a", x = unit(0.05, "npc")), par.settings = list(fontsize = list(text = 15)))
plot(p20, split = c(1, 1, 2, 2), more = TRUE)
plot(p17, split = c(2, 1, 2, 2), more = TRUE)
plot(p18, split = c(1, 2, 2, 2), more = TRUE)
plot(p19, split = c(2, 2, 2, 2), more = FALSE)
Extract a data frame of LI and SST and compute the sensitivity. Use weights=count for weighted regression.
GCMdat = GCMhspdf@data
GCMdat$northing = coordinates(GCMhspdf)[, 2]
summary(lm(theo ~ avgsst, data = GCMdat[GCMdat$avgsst >= 25, ]))
##
## Call:
## lm(formula = theo ~ avgsst, data = GCMdat[GCMdat$avgsst >= 25,
## ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.783 -2.012 -0.403 1.379 6.987
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.933 11.935 -0.66 0.5138
## avgsst 1.662 0.434 3.82 0.0011 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.64 on 20 degrees of freedom
## Multiple R-squared: 0.422, Adjusted R-squared: 0.394
## F-statistic: 14.6 on 1 and 20 DF, p-value: 0.00106
##
summary(lm(theo ~ avgsst, weights = count, data = GCMdat[GCMdat$avgsst >=
25, ]))
##
## Call:
## lm(formula = theo ~ avgsst, data = GCMdat[GCMdat$avgsst >= 25,
## ], weights = count)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -22.00 -9.03 -2.75 4.81 38.69
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -10.865 11.584 -0.94 0.35947
## avgsst 1.777 0.425 4.18 0.00046 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 14.4 on 20 degrees of freedom
## Multiple R-squared: 0.467, Adjusted R-squared: 0.44
## F-statistic: 17.5 on 1 and 20 DF, p-value: 0.000458
##
Plot the parameters as a function of SST.
xlabel = expression(paste("Sea surface temperature [", degree, "C]"))
p21 = ggplot(GCMdat, aes(avgsst, thresh)) + geom_point() + geom_smooth(method = "lm",
color = "red") + theme_bw() + xlab(xlabel) + ylab(expression(paste("u [m ",
s^-1, "]")))
p22 = ggplot(GCMdat, aes(avgsst, sigma)) + geom_point() + geom_smooth(method = "lm",
color = "red") + theme_bw() + xlab(xlabel) + ylab(expression(paste(sigma,
" [m ", s^-1, "]")))
p23 = ggplot(GCMdat, aes(avgsst, xi)) + geom_point() + geom_smooth(method = "lm",
color = "red") + theme_bw() + xlab(xlabel) + ylab(expression(paste(xi)))
p24 = ggplot(GCMdat, aes(avgsst, theo)) + geom_point() + geom_smooth(method = "lm",
color = "red") + theme_bw() + xlab(xlabel) + ylab(expression(paste("Limiting intensity [m ",
s^-1, "]")))
p21 = p21 + opts(title = "d", plot.title = theme_text(hjust = 0,
size = 15))
p22 = p22 + opts(title = "b", plot.title = theme_text(hjust = 0,
size = 15))
p23 = p23 + opts(title = "c", plot.title = theme_text(hjust = 0,
size = 15))
p24 = p24 + opts(title = "a", plot.title = theme_text(hjust = 0,
size = 15))
grid.newpage()
pushViewport(viewport(layout = grid.layout(2, 2)))
vplayout = function(x, y) viewport(layout.pos.row = x, layout.pos.col = y)
print(p24, vp = vplayout(1, 1))
print(p21, vp = vplayout(1, 2))
print(p22, vp = vplayout(2, 1))
print(p23, vp = vplayout(2, 2))
GCMm = length(begin1:end1)
minHur = 15
maxGrid = 120
minSST = 25
GCMhpt = spsample(natlwind.sdf, type = "hexagonal", n = maxGrid,
bb = bbox(natlwind.sdf) * 1.08, offset = c(1, 1))
GCMhpg = HexPoints2SpatialPolygons(GCMhpt)
GCMnp = length(GCMhpg@polygons)
GCMarea = GCMhpg@polygons[[1]]@area * 10^-6
GCMWind.hexid = over(x = natlwind.sdf, y = GCMhpg)
natlwind.sdf$whexid = GCMWind.hexid
GCMWind.split = split(natlwind.sdf@data, GCMWind.hexid)
GCMmaximum = NULL
for (i in 1:length(GCMWind.split)) {
GCMmaximum[[i]] = get.max.id(GCMWind.split[[i]])
}
GCMless = which(sapply(GCMmaximum, function(x) dim(x)[1] < minHur))
GCMmaximum = GCMmaximum[-GCMless]
GCMid = as.integer(lapply(GCMmaximum, function(x) unique(x$whexid)))
GCMhexid = paste("ID", GCMid, sep = "")
GCMdat = over(x = GCMhpg[GCMhexid], y = SST.sdf, fn = mean)
GCMcount = unlist(lapply(GCMmaximum, function(x) dim(x)[1]))
GCMdat$count = GCMcount
GCMhspdf = SpatialPolygonsDataFrame(GCMhpg[GCMhexid], GCMdat, match.ID = TRUE)
GCMobsmax = numeric()
for (i in 1:length(GCMmaximum)) {
GCMobsmax[i] = max(GCMmaximum[[i]]$WmaxS)
}
GCMhspdf$obsmax = GCMobsmax
GCMtheo = numeric()
GCMtheoL = numeric()
GCMtheoU = numeric()
for (i in 1:length(GCMmaximum)) {
GCMW = GCMmaximum[[i]]$WmaxS
boot.out = boot(GCMW, boot.theo, R = 100, nyear = GCMm)
boot.data = boot.ci(boot.out, conf = 0.8, type = "basic")
GCMtheo[i] = boot.data$t0
GCMtheoL[i] = max(max(GCMW), boot.data$basic[4])
GCMtheoU[i] = boot.data$basic[5]
}
GCMhspdf$theo = GCMtheo
GCMhspdf$theoL = GCMtheoL
GCMhspdf$theoU = GCMtheoU
GCMhspdf2 = GCMhspdf[GCMhspdf$avgsst >= minSST, ]
GCMhexnb2 = poly2nb(GCMhspdf2)
GCMwts2 = nb2listw(GCMhexnb2, style = "W", zero.policy = TRUE)
GCMstvty.lm = lm(theo ~ avgsst, weights = count, data = GCMhspdf2)
GCMarea
## [1] 477929
dim(GCMhspdf)[1]
## [1] 27
dim(GCMhspdf2)[1]
## [1] 22
summary(GCMstvty.lm)$coefficients[2, ]
## Estimate Std. Error t value Pr(>|t|)
## 1.7766591 0.4246383 4.1839356 0.0004576
Plot the trend line (slope is the sensitivity) with uncertainty lines on the LI estimates.
xlabel = expression(paste("Sea surface temperature [", degree, "C]"))
p = ggplot(GCMhspdf2@data, aes(avgsst, theo)) + geom_point() + theme_bw()
limits = aes(ymax = GCMhspdf2$theoU, ymin = GCMhspdf2$theoL)
p = p + geom_linerange(limits)
p = p + geom_smooth(method = "lm", aes(weight = count))
p + xlab(xlabel) + ylab(expression(paste("Limiting intensity [m ",
s^-1, "]")))
Save the results.
HiRAMr1 = GCMhspdf2
load("gfdlTracks_r2.df.RData")
natlwind.df = natltracks.df1[natltracks.df1$Yr >= begin1 & natltracks.df1$Yr <=
end1, ]
natlwind.df = subset(natlwind.df, WmaxS >= 10)
qthr = quantile(natlwind.df$WmaxS, prob = 0.69)
natlwind.df = subset(natlwind.df, WmaxS >= qthr)
natlwind.df$WmaxS = natlwind.df$WmaxS - natlwind.df$maguv * 0.6
length(unique(natlwind.df$Sid))
## [1] 223
mean(natlwind.df$WmaxS)
## [1] 24.93
skew(natlwind.df$WmaxS)
## [1] -0.264
mean(natlwind.df$maguv)
## [1] 11.68
natlwind.sdf = natlwind.df
coordinates(natlwind.sdf) = c("lon", "lat")
proj4string(natlwind.sdf) = CRS(ll)
natlwind.sdf = spTransform(natlwind.sdf, CRS(lcc))
GCMm = length(begin1:end1)
minHur = 15
maxGrid = 120
minSST = 25
GCMhpt = spsample(natlwind.sdf, type = "hexagonal", n = maxGrid,
bb = bbox(natlwind.sdf) * 1.08, offset = c(1, 1))
GCMhpg = HexPoints2SpatialPolygons(GCMhpt)
GCMnp = length(GCMhpg@polygons)
GCMarea = GCMhpg@polygons[[1]]@area * 10^-6
GCMWind.hexid = over(x = natlwind.sdf, y = GCMhpg)
natlwind.sdf$whexid = GCMWind.hexid
GCMWind.split = split(natlwind.sdf@data, GCMWind.hexid)
GCMmaximum = NULL
for (i in 1:length(GCMWind.split)) {
GCMmaximum[[i]] = get.max.id(GCMWind.split[[i]])
}
GCMless = which(sapply(GCMmaximum, function(x) dim(x)[1] < minHur))
GCMmaximum = GCMmaximum[-GCMless]
GCMid = as.integer(lapply(GCMmaximum, function(x) unique(x$whexid)))
GCMhexid = paste("ID", GCMid, sep = "")
GCMdat = over(x = GCMhpg[GCMhexid], y = SST.sdf, fn = mean)
GCMcount = unlist(lapply(GCMmaximum, function(x) dim(x)[1]))
GCMdat$count = GCMcount
GCMhspdf = SpatialPolygonsDataFrame(GCMhpg[GCMhexid], GCMdat, match.ID = TRUE)
GCMobsmax = numeric()
for (i in 1:length(GCMmaximum)) {
GCMobsmax[i] = max(GCMmaximum[[i]]$WmaxS)
}
GCMhspdf$obsmax = GCMobsmax
GCMtheo = numeric()
GCMtheoL = numeric()
GCMtheoU = numeric()
for (i in 1:length(GCMmaximum)) {
GCMW = GCMmaximum[[i]]$WmaxS
boot.out = boot(GCMW, boot.theo, R = 100, nyear = GCMm)
boot.data = boot.ci(boot.out, conf = 0.8, type = "basic")
GCMtheo[i] = boot.data$t0
GCMtheoL[i] = max(max(GCMW), boot.data$basic[4])
GCMtheoU[i] = boot.data$basic[5]
}
GCMhspdf$theo = GCMtheo
GCMhspdf$theoL = GCMtheoL
GCMhspdf$theoU = GCMtheoU
GCMhspdf2 = GCMhspdf[GCMhspdf$avgsst >= minSST, ]
GCMhexnb2 = poly2nb(GCMhspdf2)
GCMwts2 = nb2listw(GCMhexnb2, style = "W", zero.policy = TRUE)
GCMstvty.lm = lm(theo ~ avgsst, weights = count, data = GCMhspdf2)
GCMarea
## [1] 510967
dim(GCMhspdf)[1]
## [1] 22
dim(GCMhspdf2)[1]
## [1] 15
summary(GCMstvty.lm)$coefficients[2, ]
## Estimate Std. Error t value Pr(>|t|)
## 1.68318 0.57667 2.91877 0.01197
HiRAMr2 = GCMhspdf2
load("gfdlTracks_r3.df.RData")
natlwind.df = natltracks.df1[natltracks.df1$Yr >= begin1 & natltracks.df1$Yr <=
end1, ]
natlwind.df = subset(natlwind.df, WmaxS >= 10)
qthr = quantile(natlwind.df$WmaxS, prob = 0.69)
natlwind.df = subset(natlwind.df, WmaxS >= qthr)
natlwind.df$WmaxS = natlwind.df$WmaxS - natlwind.df$maguv * 0.6
length(unique(natlwind.df$Sid))
## [1] 211
mean(natlwind.df$WmaxS)
## [1] 25.91
skew(natlwind.df$WmaxS)
## [1] -0.2723
mean(natlwind.df$maguv)
## [1] 11.96
natlwind.sdf = natlwind.df
coordinates(natlwind.sdf) = c("lon", "lat")
proj4string(natlwind.sdf) = CRS(ll)
natlwind.sdf = spTransform(natlwind.sdf, CRS(lcc))
GCMm = length(begin1:end1)
minHur = 15
maxGrid = 120
minSST = 25
GCMhpt = spsample(natlwind.sdf, type = "hexagonal", n = maxGrid,
bb = bbox(natlwind.sdf) * 1.08, offset = c(1, 1))
GCMhpg = HexPoints2SpatialPolygons(GCMhpt)
GCMnp = length(GCMhpg@polygons)
GCMarea = GCMhpg@polygons[[1]]@area * 10^-6
GCMWind.hexid = over(x = natlwind.sdf, y = GCMhpg)
natlwind.sdf$whexid = GCMWind.hexid
GCMWind.split = split(natlwind.sdf@data, GCMWind.hexid)
GCMmaximum = NULL
for (i in 1:length(GCMWind.split)) {
GCMmaximum[[i]] = get.max.id(GCMWind.split[[i]])
}
GCMless = which(sapply(GCMmaximum, function(x) dim(x)[1] < minHur))
GCMmaximum = GCMmaximum[-GCMless]
GCMid = as.integer(lapply(GCMmaximum, function(x) unique(x$whexid)))
GCMhexid = paste("ID", GCMid, sep = "")
GCMdat = over(x = GCMhpg[GCMhexid], y = SST.sdf, fn = mean)
GCMcount = unlist(lapply(GCMmaximum, function(x) dim(x)[1]))
GCMdat$count = GCMcount
GCMhspdf = SpatialPolygonsDataFrame(GCMhpg[GCMhexid], GCMdat, match.ID = TRUE)
GCMobsmax = numeric()
for (i in 1:length(GCMmaximum)) {
GCMobsmax[i] = max(GCMmaximum[[i]]$WmaxS)
}
GCMhspdf$obsmax = GCMobsmax
GCMtheo = numeric()
GCMtheoL = numeric()
GCMtheoU = numeric()
for (i in 1:length(GCMmaximum)) {
GCMW = GCMmaximum[[i]]$WmaxS
boot.out = boot(GCMW, boot.theo, R = 100, nyear = GCMm)
boot.data = boot.ci(boot.out, conf = 0.8, type = "basic")
GCMtheo[i] = boot.data$t0
GCMtheoL[i] = max(max(GCMW), boot.data$basic[4])
GCMtheoU[i] = boot.data$basic[5]
}
GCMhspdf$theo = GCMtheo
GCMhspdf$theoL = GCMtheoL
GCMhspdf$theoU = GCMtheoU
GCMhspdf2 = GCMhspdf[GCMhspdf$avgsst >= minSST, ]
GCMhexnb2 = poly2nb(GCMhspdf2)
GCMwts2 = nb2listw(GCMhexnb2, style = "W", zero.policy = TRUE)
GCMstvty.lm = lm(theo ~ avgsst, weights = count, data = GCMhspdf2)
GCMarea
## [1] 516044
dim(GCMhspdf)[1]
## [1] 24
dim(GCMhspdf2)[1]
## [1] 17
summary(GCMstvty.lm)$coefficients[2, ]
## Estimate Std. Error t value Pr(>|t|)
## 1.30439 0.46384 2.81212 0.01313
HiRAMr3 = GCMhspdf2
Plot the sensitivity slopes from the observed data and from the three HiRAM runs using a facet grid.
A = OBS
B = HiRAMr1
C = HiRAMr2
D = HiRAMr3
pdat = rbind(A@data, B@data, C@data, D@data)
pdat$type = factor(c(rep("Data (Best-track)", dim(A@data)[1]), rep("GCM (GFDL-HiRAM run 1)",
dim(B@data)[1]), rep("GCM (GFDL-HiRAM run 2)", dim(C@data)[1]), rep("GCM (GFDL-HiRAM run 3)",
dim(D@data)[1])))
xlabel = expression(paste("Sea surface temperature [", degree, "C]"))
ylabel = expression(paste("Limiting intensity [m ", s^-1, "]"))
p25 = ggplot(pdat, aes(x = avgsst, y = theo)) + geom_point() + facet_grid(. ~
type)
p25 = p25 + theme_bw() + xlab(xlabel) + ylab(ylabel)
p25 = p25 + geom_smooth(method = "lm", aes(weight = count))
limits = aes(ymax = c(A$theoU, B$theoU, C$theoU, D$theoU), ymin = c(A$theoL,
B$theoL, C$theoL, D$theoL))
p25 + geom_errorbar(limits) + opts(aspect.ratio = 1) + facet_wrap(~type)
The model sensitivity is lower than the observed sensitivity.
sens = c(1.78, 1.68, 1.3)
sensse = c(0.425, 0.577, 0.464)
sd(sens)
## [1] 0.2532
sd(sens)/sqrt(3)
## [1] 0.1462
mean(sens)
## [1] 1.587
mean(sensse)
## [1] 0.4887
z = mean(sens)/(mean(sensse) + sd(sens)/sqrt(3))
1 - pnorm(z)
## [1] 0.006224
Rerun the sensitivity calculations using wind speeds from the HiRAM that are adjusted according to the scheme proposed in Zhao and Held (2010).
Start with run 1.
load("gfdlTracks_r1.df.RData")
natlwind.df = natltracks.df1[natltracks.df1$Yr >= begin1 & natltracks.df1$Yr <=
end1, ]
natlwind.df = subset(natlwind.df, WmaxS >= 10)
less = natlwind.df$WmaxS <= 34
natlwind.df$WmaxF = 34 + 3.15 * (natlwind.df$WmaxS - 34)
natlwind.df$WmaxF[less] = 34 + 1.22 * (natlwind.df$WmaxS[less] -
34)
natlwind.df$WmaxS = natlwind.df$WmaxF
p = ggplot(natlwind.df, aes(WmaxF)) + geom_histogram(fill = "gray",
binwidth = 2) + theme_bw()
p + xlab(expression(paste("Hurricane Intensity [m ", s^-1, "]"))) +
ylab("Hours")
qthr = quantile(natlwind.df$WmaxS, prob = 0.69)
natlwind.df = subset(natlwind.df, WmaxS >= qthr)
natlwind.df$WmaxS = natlwind.df$WmaxS - natlwind.df$maguv * 0.6
mean(natlwind.df$WmaxS)
## [1] 26.83
skew(natlwind.df$WmaxS)
## [1] 0.9967
natlwind.sdf = natlwind.df
coordinates(natlwind.sdf) = c("lon", "lat")
proj4string(natlwind.sdf) = CRS(ll)
natlwind.sdf = spTransform(natlwind.sdf, CRS(lcc))
GCMm = length(begin1:end1)
minHur = 15
maxGrid = 120
minSST = 25
GCMhpt = spsample(natlwind.sdf, type = "hexagonal", n = maxGrid,
bb = bbox(natlwind.sdf) * 1.08, offset = c(1, 1))
GCMhpg = HexPoints2SpatialPolygons(GCMhpt)
GCMnp = length(GCMhpg@polygons)
GCMarea = GCMhpg@polygons[[1]]@area * 10^-6
GCMWind.hexid = over(x = natlwind.sdf, y = GCMhpg)
natlwind.sdf$whexid = GCMWind.hexid
GCMWind.split = split(natlwind.sdf@data, GCMWind.hexid)
GCMmaximum = NULL
for (i in 1:length(GCMWind.split)) {
GCMmaximum[[i]] = get.max.id(GCMWind.split[[i]])
}
GCMless = which(sapply(GCMmaximum, function(x) dim(x)[1] < minHur))
GCMmaximum = GCMmaximum[-GCMless]
GCMid = as.integer(lapply(GCMmaximum, function(x) unique(x$whexid)))
GCMhexid = paste("ID", GCMid, sep = "")
GCMdat = over(x = GCMhpg[GCMhexid], y = SST.sdf, fn = mean)
GCMcount = unlist(lapply(GCMmaximum, function(x) dim(x)[1]))
GCMdat$count = GCMcount
GCMhspdf = SpatialPolygonsDataFrame(GCMhpg[GCMhexid], GCMdat, match.ID = TRUE)
GCMobsmax = numeric()
for (i in 1:length(GCMmaximum)) {
GCMobsmax[i] = max(GCMmaximum[[i]]$WmaxS)
}
GCMhspdf$obsmax = GCMobsmax
GCMtheo = numeric()
GCMtheoL = numeric()
GCMtheoU = numeric()
for (i in 1:length(GCMmaximum)) {
GCMW = GCMmaximum[[i]]$WmaxS
boot.out = boot(GCMW, boot.theo, R = 100, nyear = GCMm)
boot.data = boot.ci(boot.out, conf = 0.8, type = "basic")
GCMtheo[i] = boot.data$t0
GCMtheoL[i] = max(max(GCMW), boot.data$basic[4])
GCMtheoU[i] = boot.data$basic[5]
}
GCMhspdf$theo = GCMtheo
GCMhspdf$theoL = GCMtheoL
GCMhspdf$theoU = GCMtheoU
GCMhspdf2 = GCMhspdf[GCMhspdf$avgsst >= minSST, ]
GCMhexnb2 = poly2nb(GCMhspdf2)
GCMwts2 = nb2listw(GCMhexnb2, style = "W", zero.policy = TRUE)
GCMstvty.lm = lm(theo ~ avgsst, weights = count, data = GCMhspdf2[GCMhspdf2$theo <
120, ])
GCMarea
## [1] 477929
dim(GCMhspdf)[1]
## [1] 27
dim(GCMhspdf2)[1]
## [1] 22
summary(GCMstvty.lm)$coefficients[2, ]
## Estimate Std. Error t value Pr(>|t|)
## 4.5051 2.0629 2.1839 0.0442
HiRAMr1r = GCMhspdf2
Continue with run 2.
load("gfdlTracks_r2.df.RData")
natlwind.df = natltracks.df1[natltracks.df1$Yr >= begin1 & natltracks.df1$Yr <=
end1, ]
natlwind.df = subset(natlwind.df, WmaxS >= 10)
less = natlwind.df$WmaxS <= 34
natlwind.df$WmaxF = 34 + 3.15 * (natlwind.df$WmaxS - 34)
natlwind.df$WmaxF[less] = 34 + 1.22 * (natlwind.df$WmaxS[less] -
34)
natlwind.df$WmaxS = natlwind.df$WmaxF
p = ggplot(natlwind.df, aes(WmaxF)) + geom_histogram(fill = "gray",
binwidth = 2) + theme_bw()
p + xlab(expression(paste("Hurricane Intensity [m ", s^-1, "]"))) +
ylab("Hours")
qthr = quantile(natlwind.df$WmaxS, prob = 0.69)
natlwind.df = subset(natlwind.df, WmaxS >= qthr) # justification for subset criteria given above
natlwind.df$WmaxS = natlwind.df$WmaxS - natlwind.df$maguv * 0.6
mean(natlwind.df$WmaxS)
## [1] 25.98
skew(natlwind.df$WmaxS)
## [1] 0.7107
natlwind.sdf = natlwind.df
coordinates(natlwind.sdf) = c("lon", "lat")
proj4string(natlwind.sdf) = CRS(ll)
natlwind.sdf = spTransform(natlwind.sdf, CRS(lcc))
GCMm = length(begin1:end1)
minHur = 15
maxGrid = 120
minSST = 25
GCMhpt = spsample(natlwind.sdf, type = "hexagonal", n = maxGrid,
bb = bbox(natlwind.sdf) * 1.08, offset = c(1, 1))
GCMhpg = HexPoints2SpatialPolygons(GCMhpt)
GCMnp = length(GCMhpg@polygons)
GCMarea = GCMhpg@polygons[[1]]@area * 10^-6
GCMWind.hexid = over(x = natlwind.sdf, y = GCMhpg)
natlwind.sdf$whexid = GCMWind.hexid
GCMWind.split = split(natlwind.sdf@data, GCMWind.hexid)
GCMmaximum = NULL
for (i in 1:length(GCMWind.split)) {
GCMmaximum[[i]] = get.max.id(GCMWind.split[[i]])
}
GCMless = which(sapply(GCMmaximum, function(x) dim(x)[1] < minHur))
GCMmaximum = GCMmaximum[-GCMless]
GCMid = as.integer(lapply(GCMmaximum, function(x) unique(x$whexid)))
GCMhexid = paste("ID", GCMid, sep = "")
GCMdat = over(x = GCMhpg[GCMhexid], y = SST.sdf, fn = mean)
GCMcount = unlist(lapply(GCMmaximum, function(x) dim(x)[1]))
GCMdat$count = GCMcount
GCMhspdf = SpatialPolygonsDataFrame(GCMhpg[GCMhexid], GCMdat, match.ID = TRUE)
GCMobsmax = numeric()
for (i in 1:length(GCMmaximum)) {
GCMobsmax[i] = max(GCMmaximum[[i]]$WmaxS)
}
GCMhspdf$obsmax = GCMobsmax
GCMtheo = numeric()
GCMtheoL = numeric()
GCMtheoU = numeric()
for (i in 1:length(GCMmaximum)) {
GCMW = GCMmaximum[[i]]$WmaxS
boot.out = boot(GCMW, boot.theo, R = 100, nyear = GCMm)
boot.data = boot.ci(boot.out, conf = 0.8, type = "basic")
GCMtheo[i] = boot.data$t0
GCMtheoL[i] = max(max(GCMW), boot.data$basic[4])
GCMtheoU[i] = boot.data$basic[5]
}
GCMhspdf$theo = GCMtheo
GCMhspdf$theoL = GCMtheoL
GCMhspdf$theoU = GCMtheoU
GCMhspdf2 = GCMhspdf[GCMhspdf$avgsst >= minSST, ]
GCMhexnb2 = poly2nb(GCMhspdf2)
GCMwts2 = nb2listw(GCMhexnb2, style = "W", zero.policy = TRUE)
GCMstvty.lm = lm(theo ~ avgsst, weights = count, data = GCMhspdf2[GCMhspdf2$theo <
120, ])
GCMarea
## [1] 510967
dim(GCMhspdf)[1]
## [1] 22
dim(GCMhspdf2)[1]
## [1] 15
summary(GCMstvty.lm)$coefficients[2, ]
## Estimate Std. Error t value Pr(>|t|)
## 4.74788 2.45026 1.93770 0.07469
HiRAMr2r = GCMhspdf2
Finally with run 3.
load("gfdlTracks_r3.df.RData")
natlwind.df = natltracks.df1[natltracks.df1$Yr >= begin1 & natltracks.df1$Yr <=
end1, ]
natlwind.df = subset(natlwind.df, WmaxS >= 10)
less = natlwind.df$WmaxS <= 34
natlwind.df$WmaxF = 34 + 3.15 * (natlwind.df$WmaxS - 34)
natlwind.df$WmaxF[less] = 34 + 1.22 * (natlwind.df$WmaxS[less] -
34)
natlwind.df$WmaxS = natlwind.df$WmaxF
p = ggplot(natlwind.df, aes(WmaxF)) + geom_histogram(fill = "gray",
binwidth = 2) + theme_bw()
p + xlab(expression(paste("Hurricane Intensity [m ", s^-1, "]"))) +
ylab("Hours")
qthr = quantile(natlwind.df$WmaxS, prob = 0.69)
natlwind.df = subset(natlwind.df, WmaxS >= qthr) # justification for subset criteria given above
natlwind.df$WmaxS = natlwind.df$WmaxS - natlwind.df$maguv * 0.6
mean(natlwind.df$WmaxS)
## [1] 28.02
skew(natlwind.df$WmaxS)
## [1] 0.8577
natlwind.sdf = natlwind.df
coordinates(natlwind.sdf) = c("lon", "lat")
proj4string(natlwind.sdf) = CRS(ll)
natlwind.sdf = spTransform(natlwind.sdf, CRS(lcc))
GCMm = length(begin1:end1)
minHur = 15
maxGrid = 120
minSST = 25
GCMhpt = spsample(natlwind.sdf, type = "hexagonal", n = maxGrid,
bb = bbox(natlwind.sdf) * 1.08, offset = c(1, 1))
GCMhpg = HexPoints2SpatialPolygons(GCMhpt)
GCMnp = length(GCMhpg@polygons)
GCMarea = GCMhpg@polygons[[1]]@area * 10^-6
GCMWind.hexid = over(x = natlwind.sdf, y = GCMhpg)
natlwind.sdf$whexid = GCMWind.hexid
GCMWind.split = split(natlwind.sdf@data, GCMWind.hexid)
GCMmaximum = NULL
for (i in 1:length(GCMWind.split)) {
GCMmaximum[[i]] = get.max.id(GCMWind.split[[i]])
}
GCMless = which(sapply(GCMmaximum, function(x) dim(x)[1] < minHur))
GCMmaximum = GCMmaximum[-GCMless]
GCMid = as.integer(lapply(GCMmaximum, function(x) unique(x$whexid)))
GCMhexid = paste("ID", GCMid, sep = "")
GCMdat = over(x = GCMhpg[GCMhexid], y = SST.sdf, fn = mean)
GCMcount = unlist(lapply(GCMmaximum, function(x) dim(x)[1]))
GCMdat$count = GCMcount
GCMhspdf = SpatialPolygonsDataFrame(GCMhpg[GCMhexid], GCMdat, match.ID = TRUE)
GCMobsmax = numeric()
for (i in 1:length(GCMmaximum)) {
GCMobsmax[i] = max(GCMmaximum[[i]]$WmaxS)
}
GCMhspdf$obsmax = GCMobsmax
GCMtheo = numeric()
GCMtheoL = numeric()
GCMtheoU = numeric()
for (i in 1:length(GCMmaximum)) {
GCMW = GCMmaximum[[i]]$WmaxS
boot.out = boot(GCMW, boot.theo, R = 100, nyear = GCMm)
boot.data = boot.ci(boot.out, conf = 0.8, type = "basic")
GCMtheo[i] = boot.data$t0
GCMtheoL[i] = max(max(GCMW), boot.data$basic[4])
GCMtheoU[i] = boot.data$basic[5]
}
GCMhspdf$theo = GCMtheo
GCMhspdf$theoL = GCMtheoL
GCMhspdf$theoU = GCMtheoU
GCMhspdf2 = GCMhspdf[GCMhspdf$avgsst >= minSST, ]
GCMhexnb2 = poly2nb(GCMhspdf2)
GCMwts2 = nb2listw(GCMhexnb2, style = "W", zero.policy = TRUE)
GCMstvty.lm = lm(theo ~ avgsst, weights = count, data = GCMhspdf2[GCMhspdf2$theo <
120, ])
GCMarea
## [1] 516044
dim(GCMhspdf)[1]
## [1] 24
dim(GCMhspdf2)[1]
## [1] 17
summary(GCMstvty.lm)$coefficients[2, ]
## Estimate Std. Error t value Pr(>|t|)
## 6.91770 3.33042 2.07712 0.05538
HiRAMr3r = GCMhspdf2
Plot the sensitivity slopes from the observed data and from the three HiRAM runs with the adjusted wind speeds.
A = OBS
B = HiRAMr1r[HiRAMr1r$theoL < 100, ]
C = HiRAMr2r[HiRAMr2r$theoL < 100, ]
D = HiRAMr3r[HiRAMr3r$theoL < 100, ]
pdat = rbind(A@data, B@data, C@data, D@data)
pdat$type = factor(c(rep("Data (Best-track)", dim(A@data)[1]), rep("GCM (GFDL-HiRAM run 1 adjusted)",
dim(B@data)[1]), rep("GCM (GFDL-HiRAM run 2 adjusted)", dim(C@data)[1]),
rep("GCM (GFDL-HiRAM run 3 adjusted)", dim(D@data)[1])))
xlabel = expression(paste("Sea surface temperature [", degree, "C]"))
ylabel = expression(paste("Limiting intensity [m ", s^-1, "]"))
p26 = ggplot(pdat, aes(x = avgsst, y = theo)) + geom_point() + facet_grid(. ~
type)
p26 = p26 + theme_bw() + xlab(xlabel) + ylab(ylabel)
p26 = p26 + geom_smooth(method = "lm", aes(weight = count))
limits = aes(ymax = c(A$theoU, B$theoU, C$theoU, D$theoU), ymin = c(A$theoL,
B$theoL, C$theoL, D$theoL))
p26 + geom_errorbar(limits) + opts(aspect.ratio = 1) + facet_wrap(~type)
Combined sensitivity
sens = c(4.51, 4.75, 6.92)
sensse = c(2.063, 2.45, 3.33)
sd(sens)
## [1] 1.328
sd(sens)/sqrt(3)
## [1] 0.7665
mean(sens)
## [1] 5.393
mean(sensse)
## [1] 2.614
z = mean(sens)/(mean(sensse) + sd(sens)/sqrt(3))
1 - pnorm(z)
## [1] 0.05532
Sensitivity of LI to SST is higher using the adjusted wind speeds, but there is an increase in the standard error making the trends less significant against a null hypothesis of zero sensitivity.
Now repeat the sensitivity calculations using data from the FSU-COAPS model.
begin2 = 1982
end2 = 2009
Run 1
load("fsu_r1i1p1_ATL.RData")
fsuTracks.df = fsuTracks.df[fsuTracks.df$Yr >= begin2 & fsuTracks.df$Yr <=
end2, ]
fsuTracks.df = subset(fsuTracks.df, WmaxS >= 10)
natlwind.df = fsuTracks.df
p = ggplot(natlwind.df, aes(WmaxS)) + geom_histogram(fill = "gray",
binwidth = 2) + theme_bw()
p + xlab(expression(paste("Hurricane Intensity [m ", s^-1, "]"))) +
ylab("Hours")
qthr = quantile(natlwind.df$WmaxS, prob = 0.69)
natlwind.df = subset(natlwind.df, WmaxS >= qthr)
natlwind.df$WmaxS = natlwind.df$WmaxS - natlwind.df$maguv * 0.6
length(unique(natlwind.df$Sid))
## [1] 308
mean(natlwind.df$WmaxS)
## [1] 27.74
skew(natlwind.df$WmaxS)
## [1] -0.5928
mean(natlwind.df$maguv)
## [1] 13.25
natlwind.sdf = natlwind.df
coordinates(natlwind.sdf) = c("lon", "lat")
proj4string(natlwind.sdf) = CRS(ll)
natlwind.sdf = spTransform(natlwind.sdf, CRS(lcc))
GCMm = length(begin1:end1)
minHur = 15
maxGrid = 120
minSST = 25
GCMhpt = spsample(natlwind.sdf, type = "hexagonal", n = maxGrid,
bb = bbox(natlwind.sdf) * 1.08, offset = c(1, 1))
GCMhpg = HexPoints2SpatialPolygons(GCMhpt)
GCMnp = length(GCMhpg@polygons)
GCMarea = GCMhpg@polygons[[1]]@area * 10^-6
GCMWind.hexid = over(x = natlwind.sdf, y = GCMhpg)
natlwind.sdf$whexid = GCMWind.hexid
GCMWind.split = split(natlwind.sdf@data, GCMWind.hexid)
GCMmaximum = NULL
for (i in 1:length(GCMWind.split)) {
GCMmaximum[[i]] = get.max.id(GCMWind.split[[i]])
}
GCMless = which(sapply(GCMmaximum, function(x) dim(x)[1] < minHur))
GCMmaximum = GCMmaximum[-GCMless]
GCMid = as.integer(lapply(GCMmaximum, function(x) unique(x$whexid)))
GCMhexid = paste("ID", GCMid, sep = "")
GCMdat = over(x = GCMhpg[GCMhexid], y = SST.sdf, fn = mean)
GCMcount = unlist(lapply(GCMmaximum, function(x) dim(x)[1]))
GCMdat$count = GCMcount
GCMhspdf = SpatialPolygonsDataFrame(GCMhpg[GCMhexid], GCMdat, match.ID = TRUE)
GCMobsmax = numeric()
for (i in 1:length(GCMmaximum)) {
GCMobsmax[i] = max(GCMmaximum[[i]]$WmaxS)
}
GCMhspdf$obsmax = GCMobsmax
GCMtheo = numeric()
GCMtheoL = numeric()
GCMtheoU = numeric()
for (i in 1:length(GCMmaximum)) {
GCMW = GCMmaximum[[i]]$WmaxS
boot.out = boot(GCMW, boot.theo, R = 100, nyear = GCMm)
boot.data = boot.ci(boot.out, conf = 0.8, type = "basic")
GCMtheo[i] = boot.data$t0
GCMtheoL[i] = max(max(GCMW), boot.data$basic[4])
GCMtheoU[i] = boot.data$basic[5]
}
GCMhspdf$theo = GCMtheo
GCMhspdf$theoL = GCMtheoL
GCMhspdf$theoU = GCMtheoU
GCMhspdf2 = GCMhspdf[GCMhspdf$avgsst >= minSST, ]
GCMhexnb2 = poly2nb(GCMhspdf2)
GCMwts2 = nb2listw(GCMhexnb2, style = "W", zero.policy = TRUE)
GCMstvty.lm = lm(theo ~ avgsst, weights = count, data = GCMhspdf2)
GCMarea
## [1] 444124
dim(GCMhspdf)[1]
## [1] 31
dim(GCMhspdf2)[1]
## [1] 24
summary(GCMstvty.lm)$coefficients[2, ]
## Estimate Std. Error t value Pr(>|t|)
## 2.862 2.640 1.084 0.290
FSUr1 = GCMhspdf2
Run 2
load("fsu_r2i1p1_ATL.RData")
fsuTracks.df = fsuTracks.df[fsuTracks.df$Yr >= begin2 & fsuTracks.df$Yr <=
end2, ]
fsuTracks.df = subset(fsuTracks.df, WmaxS >= 10)
natlwind.df = fsuTracks.df
p = ggplot(natlwind.df, aes(WmaxS)) + geom_histogram(fill = "gray",
binwidth = 2) + theme_bw()
p + xlab(expression(paste("Hurricane Intensity [m ", s^-1, "]"))) +
ylab("Hours")
qthr = quantile(natlwind.df$WmaxS, prob = 0.69)
natlwind.df = subset(natlwind.df, WmaxS >= qthr)
natlwind.df$WmaxS = natlwind.df$WmaxS - natlwind.df$maguv * 0.6
length(unique(natlwind.df$Sid))
## [1] 318
mean(natlwind.df$WmaxS)
## [1] 28.15
skew(natlwind.df$WmaxS)
## [1] -0.5113
mean(natlwind.df$maguv)
## [1] 12.86
natlwind.sdf = natlwind.df
coordinates(natlwind.sdf) = c("lon", "lat")
proj4string(natlwind.sdf) = CRS(ll)
natlwind.sdf = spTransform(natlwind.sdf, CRS(lcc))
GCMm = length(begin1:end1)
minHur = 15
maxGrid = 120
minSST = 25
GCMhpt = spsample(natlwind.sdf, type = "hexagonal", n = maxGrid,
bb = bbox(natlwind.sdf) * 1.08, offset = c(1, 1))
GCMhpg = HexPoints2SpatialPolygons(GCMhpt)
GCMnp = length(GCMhpg@polygons)
GCMarea = GCMhpg@polygons[[1]]@area * 10^-6
GCMWind.hexid = over(x = natlwind.sdf, y = GCMhpg)
natlwind.sdf$whexid = GCMWind.hexid
GCMWind.split = split(natlwind.sdf@data, GCMWind.hexid)
GCMmaximum = NULL
for (i in 1:length(GCMWind.split)) {
GCMmaximum[[i]] = get.max.id(GCMWind.split[[i]])
}
GCMless = which(sapply(GCMmaximum, function(x) dim(x)[1] < minHur))
GCMmaximum = GCMmaximum[-GCMless]
GCMid = as.integer(lapply(GCMmaximum, function(x) unique(x$whexid)))
GCMhexid = paste("ID", GCMid, sep = "")
GCMdat = over(x = GCMhpg[GCMhexid], y = SST.sdf, fn = mean)
GCMcount = unlist(lapply(GCMmaximum, function(x) dim(x)[1]))
GCMdat$count = GCMcount
GCMhspdf = SpatialPolygonsDataFrame(GCMhpg[GCMhexid], GCMdat, match.ID = TRUE)
GCMobsmax = numeric()
for (i in 1:length(GCMmaximum)) {
GCMobsmax[i] = max(GCMmaximum[[i]]$WmaxS)
}
GCMhspdf$obsmax = GCMobsmax
GCMtheo = numeric()
GCMtheoL = numeric()
GCMtheoU = numeric()
for (i in 1:length(GCMmaximum)) {
GCMW = GCMmaximum[[i]]$WmaxS
boot.out = boot(GCMW, boot.theo, R = 100, nyear = GCMm)
boot.data = boot.ci(boot.out, conf = 0.8, type = "basic")
GCMtheo[i] = boot.data$t0
GCMtheoL[i] = max(max(GCMW), boot.data$basic[4])
GCMtheoU[i] = boot.data$basic[5]
}
GCMhspdf$theo = GCMtheo
GCMhspdf$theoL = GCMtheoL
GCMhspdf$theoU = GCMtheoU
GCMhspdf2 = GCMhspdf[GCMhspdf$avgsst >= minSST, ]
GCMhexnb2 = poly2nb(GCMhspdf2)
GCMwts2 = nb2listw(GCMhexnb2, style = "W", zero.policy = TRUE)
GCMstvty.lm = lm(theo ~ avgsst, weights = count, data = GCMhspdf2)
GCMarea
## [1] 450631
dim(GCMhspdf)[1]
## [1] 29
dim(GCMhspdf2)[1]
## [1] 23
summary(GCMstvty.lm)$coefficients[2, ]
## Estimate Std. Error t value Pr(>|t|)
## 0.6417 0.5513 1.1639 0.2575
FSUr2 = GCMhspdf2
Run 3
load("fsu_r3i1p1_ATL.RData")
fsuTracks.df = fsuTracks.df[fsuTracks.df$Yr >= begin2 & fsuTracks.df$Yr <=
end2, ]
fsuTracks.df = subset(fsuTracks.df, WmaxS >= 10)
natlwind.df = fsuTracks.df
p = ggplot(natlwind.df, aes(WmaxS)) + geom_histogram(fill = "gray",
binwidth = 2) + theme_bw()
p + xlab(expression(paste("Hurricane Intensity [m ", s^-1, "]"))) +
ylab("Hours")
qthr = quantile(natlwind.df$WmaxS, prob = 0.69)
natlwind.df = subset(natlwind.df, WmaxS >= qthr)
natlwind.df$WmaxS = natlwind.df$WmaxS - natlwind.df$maguv * 0.6
length(unique(natlwind.df$Sid))
## [1] 283
mean(natlwind.df$WmaxS)
## [1] 28.43
skew(natlwind.df$WmaxS)
## [1] -0.6182
mean(natlwind.df$maguv)
## [1] 12.7
natlwind.sdf = natlwind.df
coordinates(natlwind.sdf) = c("lon", "lat")
proj4string(natlwind.sdf) = CRS(ll)
natlwind.sdf = spTransform(natlwind.sdf, CRS(lcc))
GCMm = length(begin1:end1)
minHur = 15
maxGrid = 120
minSST = 25
GCMhpt = spsample(natlwind.sdf, type = "hexagonal", n = maxGrid,
bb = bbox(natlwind.sdf) * 1.08, offset = c(1, 1))
GCMhpg = HexPoints2SpatialPolygons(GCMhpt)
GCMnp = length(GCMhpg@polygons)
GCMarea = GCMhpg@polygons[[1]]@area * 10^-6
GCMWind.hexid = over(x = natlwind.sdf, y = GCMhpg)
natlwind.sdf$whexid = GCMWind.hexid
GCMWind.split = split(natlwind.sdf@data, GCMWind.hexid)
GCMmaximum = NULL
for (i in 1:length(GCMWind.split)) {
GCMmaximum[[i]] = get.max.id(GCMWind.split[[i]])
}
GCMless = which(sapply(GCMmaximum, function(x) dim(x)[1] < minHur))
GCMmaximum = GCMmaximum[-GCMless]
GCMid = as.integer(lapply(GCMmaximum, function(x) unique(x$whexid)))
GCMhexid = paste("ID", GCMid, sep = "")
GCMdat = over(x = GCMhpg[GCMhexid], y = SST.sdf, fn = mean)
GCMcount = unlist(lapply(GCMmaximum, function(x) dim(x)[1]))
GCMdat$count = GCMcount
GCMhspdf = SpatialPolygonsDataFrame(GCMhpg[GCMhexid], GCMdat, match.ID = TRUE)
GCMobsmax = numeric()
for (i in 1:length(GCMmaximum)) {
GCMobsmax[i] = max(GCMmaximum[[i]]$WmaxS)
}
GCMhspdf$obsmax = GCMobsmax
GCMtheo = numeric()
GCMtheoL = numeric()
GCMtheoU = numeric()
for (i in 1:length(GCMmaximum)) {
GCMW = GCMmaximum[[i]]$WmaxS
boot.out = boot(GCMW, boot.theo, R = 100, nyear = GCMm)
boot.data = boot.ci(boot.out, conf = 0.8, type = "basic")
GCMtheo[i] = boot.data$t0
GCMtheoL[i] = max(max(GCMW), boot.data$basic[4])
GCMtheoU[i] = boot.data$basic[5]
}
GCMhspdf$theo = GCMtheo
GCMhspdf$theoL = GCMtheoL
GCMhspdf$theoU = GCMtheoU
GCMhspdf2 = GCMhspdf[GCMhspdf$avgsst >= minSST, ]
GCMhexnb2 = poly2nb(GCMhspdf2)
GCMwts2 = nb2listw(GCMhexnb2, style = "W", zero.policy = TRUE)
GCMstvty.lm = lm(theo ~ avgsst, weights = count, data = GCMhspdf2)
GCMarea
## [1] 490219
dim(GCMhspdf)[1]
## [1] 25
dim(GCMhspdf2)[1]
## [1] 21
summary(GCMstvty.lm)$coefficients[2, ]
## Estimate Std. Error t value Pr(>|t|)
## 0.4925 0.3715 1.3257 0.2007
FSUr3 = GCMhspdf2
Plot the sensitivity slopes from the observed data and from the three FSU-COAPS runs.
A = OBS
B = FSUr1[FSUr1$theoL < 100, ]
C = FSUr2[FSUr2$theoL < 100, ]
D = FSUr3[FSUr3$theoL < 100, ]
pdat = rbind(A@data, B@data, C@data, D@data)
pdat$type = factor(c(rep("Data (Best-track)", dim(A@data)[1]), rep("GCM (FSU-COAPS run 1)",
dim(B@data)[1]), rep("GCM (FSU-COAPS run 2)", dim(C@data)[1]), rep("GCM (FSU-COAPS run 3)",
dim(D@data)[1])))
xlabel = expression(paste("Sea surface temperature [", degree, "C]"))
ylabel = expression(paste("Limiting intensity [m ", s^-1, "]"))
p26 = ggplot(pdat, aes(x = avgsst, y = theo)) + geom_point() + facet_grid(. ~
type)
p26 = p26 + theme_bw() + xlab(xlabel) + ylab(ylabel)
p26 = p26 + geom_smooth(method = "lm", aes(weight = count))
limits = aes(ymax = c(A$theoU, B$theoU, C$theoU, D$theoU), ymin = c(A$theoL,
B$theoL, C$theoL, D$theoL))
p26 + geom_errorbar(limits) + opts(aspect.ratio = 1) + facet_wrap(~type)
