We previously examined the sensitivity of limiting intensity (LI) to sea surface temerature (SST) for observed and model-generated tropical cyclones (TCs). For the FSU-COAPS and the GFDL-HiRAM, the sensitivity of LI to SST is much lower than it is for observed TCs. The fact that these models do not produce TCs with wind speeds exceeding 50 m s\( ^{-1} \) limits the sensitivity, which is expressed here as the slope from a regression of LI onto SST. Can a higher resolution global climate model (GCM) – the 0.25\( ^{\circ} \) CAM – that generates TCs with wind speeds exceeding 70 m/ s\( ^{-1} \) reproduce the observed relationship between LI and SST?
date()
## [1] "Tue Mar 04 14:47:50 2014"
Below we list the R packages and other R scripts used in this comparison. Please contact Sarah Strazzo at ses09e@fsu.edu if you'd like to view any of the R scripts sourced below.
setwd("C:/Users/Sarah/Dropbox/SpatialExtremes")
require(ggplot2)
require(rgdal)
require(ismev)
require(ncdf)
library(mapproj)
library(mapdata)
library(maptools)
library(maps)
library(grid)
require(plyr)
library(colorRamps)
library(RColorBrewer)
require(spdep)
require(psych)
require(boot)
source("C:/Users/Sarah/Dropbox/SpatialExtremes/selectData.R")
source("C:/Users/Sarah/Dropbox/SpatialExtremes/ismevplus.R")
source("C:/Users/Sarah/Dropbox/SpatialExtremes/Kerrycanes/sensPoly.R")
source("ohcFunctions.R")
source("C:/Users/Sarah/Dropbox/SpatialExtremes/sensitivityUncertainty.R")
We first load the hourly interpolated best-track data, convert the wind speeds and forward motion from knots to m s\( ^{-1} \), and subtract 60% of the forward speed to get an estimate of rotational velocity. Note that we don't subset on hurricanes here as was done previously.
Similarly, we use hourly interpolated CAM (0.25\( ^{\circ} \) resolution) TCs and remove 60% of the forward motion as an estimate of rotational velocity. Note that there are 304 observed and 387 simulated TCs over the 1979–2005 period. Also note that the 60th-80th quantiles of observed and modeled TC wind speeds, which cover the range from approximately 25 to 35 m s\( ^{-1} \), match remarkably well. In Fig. 1, we plot both the wind speed densities and cumulative distributions for best-track and CAM TC wind speeds.
load("C:/Users/Sarah/Dropbox/SpatialExtremes/best.use.RData")
begin = 1979
end = 2005
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))
## 60% 61% 62% 63% 64% 65% 66% 67% 68% 69% 70% 71%
## 28.60 29.11 29.66 30.18 30.65 31.00 31.48 32.08 32.68 33.17 33.42 33.55
## 72% 73% 74% 75% 76% 77% 78% 79% 80%
## 33.86 34.42 35.03 35.68 36.10 36.70 37.41 38.13 38.64
Wind.df$maguv = Wind.df$maguv * 0.5144
Wind.df$WmaxS = Wind.df$WmaxS - Wind.df$maguv * 0.6
load("C:/Users/Sarah/Dropbox/GCMTrackData/CAMTracks_NAtl.RData")
natlwind.df = na.df1[na.df1$Yr >= begin & na.df1$Yr <= end, ]
quantile(natlwind.df$WmaxS, prob = seq(0.6, 0.8, 0.01))
## 60% 61% 62% 63% 64% 65% 66% 67% 68% 69% 70% 71%
## 31.05 31.46 31.88 32.31 32.74 33.20 33.64 34.09 34.52 35.02 35.52 36.04
## 72% 73% 74% 75% 76% 77% 78% 79% 80%
## 36.61 37.25 37.84 38.42 38.96 39.54 40.19 40.79 41.36
natlwind.df$WmaxS = natlwind.df$WmaxS - natlwind.df$maguv * 0.6
length(unique(Wind.df$Sid))
## [1] 304
length(unique(natlwind.df$Sid))
## [1] 387
Figure 1: Comparison of wind speed densities (top) and cumulative distributions (bottom) for best-track and CAM TCs. Top: Frequency of hourly TC records. Rotational intensity shown in m/s. Bottom: Cumulative Distribution of hourly hurricane records
Examining descriptive statistics for observed and modeled rotational wind speed, we find that the mean CAM wind speed of 24.23 m s\( ^{-1} \) compares very well with the mean best-track wind speed of 24.47 m s\( ^{-1} \). Both distributions are positively skewed, though the best-track data are slightly more skewed than the CAM data (1.02 and 0.70 for best-track and CAM respectively).
mC = mean(natlwind.df$WmaxS) # 24.24
mO = mean(Wind.df$WmaxS) # 24.47
sC = skew(natlwind.df$WmaxS) # 0.70
sO = skew(Wind.df$WmaxS) # 1.03
We next add SST data from NOAA's extended reconstructed (ERSST) data set. The data were originally obtained from www.esrl.noaa.gov/psd/data/gridded/data.noaa.ersst.html. The SST data are assigned a lambert conformal conic projection centered at 60\( ^{\circ} \) longitude.
nc = open.ncdf("C:/Users/Sarah/Dropbox/SpatialExtremes/sst.mnmean.nc")
avgsst = reformat.SST(nc, begin, end, 8, 9, 10)
rm(nc)
ll = "+proj=longlat +datum=WGS84"
lcc = "+proj=lcc +lat_1=60 +lat_2=30 +lon_0=-60"
coordinates(avgsst) = c("lon", "lat")
proj4string(avgsst) = CRS(ll)
SST.sdf = spTransform(avgsst, CRS(lcc))
As with the SST data, both TC data sets are assigned a lambert conformal conic projection.
Wind.sdf = Wind.df
coordinates(Wind.sdf) = c("lon", "lat")
proj4string(Wind.sdf) = CRS(ll)
Wind.sdf = spTransform(Wind.sdf, CRS(lcc))
natlwind.sdf = natlwind.df
coordinates(natlwind.sdf) = c("lon", "lat")
proj4string(natlwind.sdf) = CRS(ll)
natlwind.sdf = spTransform(natlwind.sdf, CRS(lcc))
Finally, we generate a tessellation of equal-area hexagons that cover the North Atlantic Basin. Each hexagon has an area of 753,556 km\( ^{2} \). Hexagons are used in place of a rectangular grid because they are slightly more efficient at covering curved TC tracks. The method is particularly useful for binning TC data so that various descriptive statistics can be summarized and easily compared between multiple data sets. As with any method, the hexagon binning comes with some limitations. The specific spatial configuration of the grid (i.e. how large are the hexagons? Where exactly is the grid centered?) does have some influence on per grid statistics. However, this uncertainty can easily be expressed, as is done toward the end of this document.
# bbox(Wind.sdf) bbox(natlwind.sdf) # larger area covered by CAM tracks
set.seed(1025)
bb = bbox(natlwind.sdf)
bb[2, 1] = bb[2, 1] * 0.8
bb[1, 1] = bb[1, 1] * 1.35
bb[1, 2] = bb[1, 2] * 1.08
bb[2, 2] = bb[2, 2] * 1.08
hpt = spsample(natlwind.sdf, type = "hexagonal", n = 150, bb = bb)
hpg = HexPoints2SpatialPolygons(hpt)
# plot(hpg) plot(Wind.sdf, add=TRUE, pch='.', col='red')
# plot(natlwind.sdf, add=TRUE, pch='.', col='red')
np = length(hpg@polygons) # 96
area = hpg@polygons[[1]]@area * 10^-6
# 753556.2 km^2
Once the hexagon lattice has been generated, we overlay the observed and simulated track points onto the grid to obtain per hexagon TC counts and maximum wind speeds. We use the same set of hexagons for both the best-track and CAM data.
hexidOBS = over(x = Wind.sdf, y = hpg)
Wind.sdf$hexid = hexidOBS
hexidCAM = over(x = natlwind.sdf, y = hpg)
natlwind.sdf$hexid = hexidCAM
OBSstats = ddply(Wind.sdf@data, .(hexid, Sid), summarize, nhr = length(Sid),
WmaxS = max(WmaxS))
OBSstats2 = ddply(OBSstats, .(hexid), summarize, countOBS = length(Sid), WmaxOBS = max(WmaxS))
CAMstats = ddply(natlwind.sdf@data, .(hexid, Sid), summarize, nhr = length(Sid),
WmaxS = max(WmaxS))
CAMstats2 = ddply(CAMstats, .(hexid), summarize, countCAM = length(Sid), WmaxCAM = max(WmaxS))
Below, we remove grids with fewer than 15 observed/simulated TCs to ensure we have sufficient data to estimate the limiting intensity model. We also remove grids that contain more than 50% land. We also remove several grids (mostly located in the far northern portion of the basin) for which there are too few strong TCs to obtain an physically meaningful estimate of limiting intensity using the statistical model described below.
hexData = merge(OBSstats2, CAMstats2)
hexData = subset(hexData, countOBS >= 15 & countCAM >= 15)
hexData = subset(hexData, hexid != 19) # No CAM convergence
hexData = subset(hexData, hexid != 39) # Africa
hexData = subset(hexData, hexid != 74) # Canada
hexData = subset(hexData, hexid != 50)
hexData = subset(hexData, hexid != 58) # Over land
hexData = subset(hexData, hexid != 30) # Over land (C. America)
hexData = subset(hexData, hexid != 75) # Newfoundland
hexData = subset(hexData, hexid != 76) # Way up north
hexData = subset(hexData, hexid != 77) # Way up north
hexData = subset(hexData, hexid != 24) # too few strong TCs for LI mod
hexData = subset(hexData, hexid != 23) # too few strong TCs for LI mod
hexData = subset(hexData, hexid != 18) # too few strong TCs for LI mod
row.names(hexData) = paste("ID", hexData$hexid, sep = "")
hspdf = SpatialPolygonsDataFrame(hpg[row.names(hexData)], hexData)
avgsst = over(x = hpg[row.names(hexData)], y = SST.sdf, fn = mean)
hspdf$avgsst = avgsst$avgsst
If we calculate the correlation between the observed and modeled TC counts, we find no significant correlation. However, if we consider only those TCs with wind speeds exceeding 33 m s\( ^{-1} \), then the observed and modeled counts are weakly correlated. In the “extras” section below, we also examine the effects of the specific hexagon configuration on these statistics.
cor.test(hexData$countOBS, hexData$countCAM) # no corr
##
## Pearson's product-moment correlation
##
## data: hexData$countOBS and hexData$countCAM
## t = -1.665, df = 20, p-value = 0.1115
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.67172 0.08522
## sample estimates:
## cor
## -0.3489
cor.test(hexData$WmaxOBS, hexData$WmaxCAM) # 0.758
##
## Pearson's product-moment correlation
##
## data: hexData$WmaxOBS and hexData$WmaxCAM
## t = 4.055, df = 20, p-value = 0.0006189
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.3489 0.8520
## sample estimates:
## cor
## 0.6717
In Fig. 2 we map the per hexagon frequency and intensity for both best-track and CAM TCs. The lack of correlation between observed and modeled TC frequency is clear in this figure. The CAM is very aggressive in developing TCs near the west coast of Africa compared to observations. As with the GFDL-HiRAM and FSU-COAPS, the frequency of CAM-generated TCs over the Gulf of Mexico is much lower than what we observe. Unlike the lower resolution models, however, the CAM does manage to generate TCs with wind speeds in the 70–75 m s\( ^{-1} \) range. As with observed TCs, these stronger model-generated TCs tend to occur over the Gulf of Mexico.
Figure 2: Frequency and intensity maps for grids with at least 15 observed AND modeled TCs.
For each grid with at least 15 observed and CAM-generated TCs, we use a statistical model that combines a Poisson distribution and a generalized Pareto distribution (GPD) to estimate the per grid LI. The Poisson distribution describes the number of TCs with maximum wind speeds that exceed some user-defined threshold, while the GPD describes the wind speed distribution for all TCs with maximum winds exceeding this threshold. We use the set of per storm maximum winds within each grid for this analysis.
The combined Poisson-GPD model for estimating LI is formed for each hexagon and for both observed and model-generated TCs. The model is defined by two parameters: scale (\( \sigma \)) and shape (\( \xi \)). We also must select the maximum wind speed threshold (a free parameter). Selecting an appropriate threshold is difficult and almost always subjective, especially when working with smaller data sets such as ours. The goal is to find a threshold that is small enough that there are sufficient data with which to estimate the model, and large enough that the data fall into the tail of the distribution (i.e. follow a GPD). After some experimenting (some of which can be seen in the “extras” section), we settled on the 69th percentile, which corresponds to 33 m s\( ^{-1} \) from the overall observed TC data set. This threshold essentially represents the highest possible value we could choose without be limited by too few data.
“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.”
oDat = OBSstats
oDat = merge(oDat, hexData, by = "hexid")
oDat = split(oDat, oDat$hexid)
cDat = CAMstats
cDat = merge(cDat, hexData, by = "hexid")
cDat = split(cDat, cDat$hexid)
m = length(begin:end)
lambda = numeric()
sigma = numeric()
xi = numeric()
thresh = numeric()
rate = numeric()
theo = numeric()
for (i in 1:length(oDat)) {
W = oDat[[i]]$WmaxS
rate[i] = length(W)/m
thresh[i] = quantile(W, prob = 0.69)
model = exceed.pp.plot(W, npy = rate[i], threshold = thresh[i], yunits = "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]
}
# length(xi) # 22 sum(xi < 0) # 22 sum(xi < -1) # 14
hspdf$threshOBS = thresh
hspdf$rateOBS = rate
hspdf$lambdaOBS = lambda
hspdf$sigmaOBS = sigma
hspdf$xiOBS = xi
hspdf$theoOBS = theo
range(hspdf@data$theoOBS)
## [1] 47.84 78.70
m = length(begin:end)
lambda = numeric()
sigma = numeric()
xi = numeric()
thresh = numeric()
rate = numeric()
theo = numeric()
for (i in 1:length(cDat)) {
W = cDat[[i]]$WmaxS
rate[i] = length(W)/m
thresh[i] = quantile(W, prob = 0.69)
model = exceed.pp.plot(W, npy = rate[i], threshold = thresh[i], yunits = "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]
}
# length(xi) # 22 sum(xi < 0) # 22 sum(xi < -1) # 13
hspdf$threshCAM = thresh
hspdf$rateCAM = rate
hspdf$lambdaCAM = lambda
hspdf$sigmaCAM = sigma
hspdf$xiCAM = xi
hspdf$theoCAM = theo
range(hspdf@data$theoCAM)
## [1] 47.62 75.78
Now that we have both an average ASO SST and a limiting intensity for each hexagon (for both observed and modeled data), we can finally calculate the “sensitivity” of limiting intensity to SST. We should clarify that in this context, “sensitivity” simply refers to the slope of the regression of LI onto SST. We are not running a climate model forced by a specific set of SSTs, re-running the model with a different set of SSTs, and then analyzing the difference in model output from the two runs as is often done in model sensitivity studies. We are in a sense summarizing the spatial relationship between SST and hurricane intensity. We would expect a positive relationship with intensity increasing as you move farther west in the basin over warmer waters near the gulf stream, Gulf of Mexico, and Caribbean. Indeed, the sensitivity of LI to SST for observed TCs reflects this thinking very well with a highly significant value of 7.46 +/- 1.65 m s\( ^{-1} \) (standard error, s.e.). Interestingly, although the CAM manages to generate fairly strong TCs, the relationship is weaker, with a sensitivity of 3.6 +/- 1.6 (s.e.). The CAM appears to capture this relationship better than lower resolution models, but still falls short of the rather high observed sensitivity. Of course, it is fair to question whether these values change if the specific configuration of hexagons changes. We'll examine this later in the “extras” section.
dat = as.data.frame(hspdf)
dat = subset(dat, avgsst >= 25)
lmOBS.theo = lm(theoOBS ~ avgsst, weights = countOBS, data = dat)
summary(lmOBS.theo)
##
## Call:
## lm(formula = theoOBS ~ avgsst, data = dat, weights = countOBS)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -59.53 -31.16 -3.59 26.65 94.61
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -141.36 41.68 -3.39 0.00373 **
## avgsst 7.46 1.49 5.00 0.00013 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 40.4 on 16 degrees of freedom
## Multiple R-squared: 0.61, Adjusted R-squared: 0.585
## F-statistic: 25 on 1 and 16 DF, p-value: 0.000131
## 7.46 +/- 1.65 (.189)
lmCAM.theo = lm(theoCAM ~ avgsst, weights = countCAM, data = dat)
summary(lmCAM.theo)
##
## Call:
## lm(formula = theoCAM ~ avgsst, data = dat, weights = countCAM)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -87.73 -36.63 4.44 37.60 83.45
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -38.92 48.89 -0.80 0.438
## avgsst 3.61 1.76 2.05 0.057 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 48.5 on 16 degrees of freedom
## Multiple R-squared: 0.208, Adjusted R-squared: 0.159
## F-statistic: 4.21 on 1 and 16 DF, p-value: 0.057
# 3.60 +/- 1.76
lmOBS.max = lm(WmaxOBS ~ avgsst, weights = countOBS, data = dat)
## 8.34 +/- 1.59
lmCAM.max = lm(WmaxCAM ~ avgsst, weights = countCAM, data = dat)
## 3.00 +/- 1.49
Given that the wind speed distributions for observed and modeled TCs do not match perfectly (see Fig. 1), perhaps it is not reasonable to set the same threshold for both sets of data. To account for this, we next try using the same threshold for the observed TCs, but calculating different thresholds for the CAM grids. First, we take the difference between the per grid maximum observed wind speed and the threshold (the 69th percentile from the wind speed distribution for each grid). We then find the threshold for the CAM grids by subtracting that difference from the per grid maximum CAM wind speed. This ensures that the difference between the threshold value and the per grid maximum wind speed is the same for observed and modeled TCs.
thresholds = matrix(0, ncol = 3, nrow = length(oDat))
thresholds = data.frame(thresholds)
for (i in 1:length(oDat)) {
a = oDat[[i]]
b = length(unique(a$Sid))
c = quantile(a$WmaxS, prob = 0.69)
d = max(a$WmaxS) - c
thresholds[i, 1] = b
thresholds[i, 2] = c
thresholds[i, 3] = d
}
colnames(thresholds) = c("countOBS", "threshOBS", "diffOBS")
rownames(thresholds) = rownames(hexData)
thresholds$threshCAM = hexData$WmaxCAM - thresholds$diffOBS
thresholds$DIFF = thresholds$threshOBS - thresholds$threshCAM
m = length(begin:end)
lambda = numeric()
sigma = numeric()
xi = numeric()
thresh = thresholds$threshCAM
rate = numeric()
theo = numeric()
for (i in 1:length(cDat)) {
W = cDat[[i]]$WmaxS
rate[i] = length(W)/m
thr = thresh[i]
model = exceed.pp.plot(W, npy = rate[i], threshold = thr, yunits = "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] = thr - sigma[i]/xi[i]
}
# length(xi) # 22 sum(xi < 0) # 22 sum(xi < -1) # 10
hspdf$threshCAM2 = thresh
hspdf$rateCAM2 = rate
hspdf$lambdaCAM2 = lambda
hspdf$sigmaCAM2 = sigma
hspdf$xiCAM2 = xi
hspdf$theoCAM2 = theo
range(hspdf@data$theoCAM)
## [1] 47.62 75.78
range(hspdf@data$theoCAM2)
## [1] 47.62 79.86
range(hspdf$theoOBS)
## [1] 47.84 78.70
We may now recalculate the sensitivity of LI to SST using the adjusted LI values for t0he CAM. As expected, by allowing lower thresholds for the CAM, we appear to have increased the range of LIs. The sensitivity of LI to SST for CAM-generated TCs is now within the standard error bounds of the observed sensitivity (5.95 +/- 1.50 m s${-1}).
dat = as.data.frame(hspdf)
dat = subset(dat, avgsst >= 25)
lmCAM.theo2 = lm(theoCAM2 ~ avgsst, weights = countCAM, data = dat)
summary(lmCAM.theo2)
##
## Call:
## lm(formula = theoCAM2 ~ avgsst, data = dat, weights = countCAM)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -79.22 -17.94 -6.97 37.43 68.34
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -102.65 41.72 -2.46 0.0256 *
## avgsst 5.95 1.50 3.96 0.0011 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 41.3 on 16 degrees of freedom
## Multiple R-squared: 0.495, Adjusted R-squared: 0.464
## F-statistic: 15.7 on 1 and 16 DF, p-value: 0.00112
# 5.95 +/- 1.50
We have thus estimated the sensitivity of LI to SST using three different methods: 1) a regression of the per hexagon maximum TC wind speed onto SST, 2) a regression of the per hexagon LI defined using the 69th percentile of wind speeds onto SST, and 3) a regression of the per hexagon LI defined using the 69th percentile for observed TCs and an adjusted threshold for modeled TCs onto SST. These regressions are plotted below for both observed and modeled TCs.
Figure 3: Comparing the sensitivity of LI to SST. Top panels use LI obtained by setting a single 30 m s\( ^{-1} \) threshold for both best-track and CAM TCs. The middle panels use the same 30 m s\( ^{-1} \) LI for best-track TCs, but an adjusted LI for the CAM. The bottom panels show the sensitivity of per hexagon maximum intensity to SST
How does the specific configuration of the hexagons affect the per hexagon counts and wind speeds? We generate 100 grids, each with a different spatial offset. For each grid, we calculate the minimum per hexagon TC count (for best-track and the CAM), the maximum per hexagon TC count, the lowest per hexagon maximum intensity, and the highest per hexagon maximum intensity. The ranges of each of these statistics are given in the table below.
allGrids = matrix(0, ncol = 8, nrow = 100)
allGrids = data.frame(allGrids)
set.seed(1025)
bb = bbox(natlwind.sdf)
bb[2, 1] = bb[2, 1] * 0.8
bb[1, 1] = bb[1, 1] * 1.35
bb[1, 2] = bb[1, 2] * 1.08
bb[2, 2] = bb[2, 2] * 1.08
for (i in 1:100) {
print(i)
hpts = spsample(natlwind.sdf, type = "hexagonal", n = 150, bb = bb)
hpg = HexPoints2SpatialPolygons(hpts)
allGrids[i, ] = get.countsWind(Wind.sdf, natlwind.sdf, hpg)
}
colnames(allGrids) = c("oCountMin", "oCountMax", "mCountMin", "mCountMax", "oWmin",
"oWmax", "mWmin", "mWmax")
range(allGrids$oCountMin) # 15 23
range(allGrids$oCountMax) # 63 74
range(allGrids$mCountMin) # 15 24
range(allGrids$mCountMax) # 64 89
range(allGrids$oWmin) # 18.2 34.2
range(allGrids$oWmax) # 78.7 78.7
range(allGrids$mWmin) # 22.5 41.0
range(allGrids$mWmax) # 69.1 73.8
| Type | Min Count | Max Count | Min Wmax | Max Wmax |
|---|---|---|---|---|
| Obs | 15, 23 | 63, 74 | 18.2, 34.1 | 78.7, 78.7 |
| CAM | 15, 24 | 64, 89 | 22.5, 41.0 | 69.1, 73.8 |
How does the specific configuration of the hexagons affect the sensitivity? To test this we generate 100 grids, each one with a slightly different spatial offset. To avoid issues that arise when defining the extreme value model, we estimate LI using the per hexagon maximum intensity.
allSensCAM = matrix(0, ncol = 2, nrow = 100)
allSensCAM = data.frame(allSensCAM)
set.seed(4025)
bb = bbox(natlwind.sdf)
bb[2, 1] = bb[2, 1] * 0.8
bb[1, 1] = bb[1, 1] * 1.35
bb[1, 2] = bb[1, 2] * 1.08
bb[2, 2] = bb[2, 2] * 1.08
for (i in 1:100) {
# print(i)
hpts = spsample(natlwind.sdf, type = "hexagonal", n = 150, bb = bb)
hpg = HexPoints2SpatialPolygons(hpts)
allSensCAM[i, ] = get.sensitivity3(Wind.sdf, natlwind.sdf, hpg, SST.sdf,
1979, 2005)
}
colnames(allSensCAM) = c("OBS", "CAM")
median(allSensCAM$OBS)
## [1] 7.924
median(allSensCAM$CAM)
## [1] 1.717
