Code/text in support of our paper to be submitted to Journal of Advances in Modeling Earth Systems
http://james.agu.org/index.php/JAMES/article/view/31
Obtain the required packages and define a function for obtaining the maximum cyclone intensity per grid:
date()
## [1] "Thu Jul 04 10:11:09 2013"
# setwd('~/Dropbox/Sensitivity/Kerrycanes')
setwd("C:/Users/Sarah/Dropbox/SpatialExtremes/Kerrycanes")
require(R.matlab)
require(R.utils)
require(ggplot2)
require(rgdal)
require(ismev)
require(ncdf)
require(plyr)
library(mapproj)
library(mapdata)
library(maptools)
library(maps)
library(grid)
library(colorRamps)
library(RColorBrewer)
require(spdep)
require(psych)
require(boot)
require(spgwr)
source("selectData.R")
source("ismevplus.R")
source("sensPoly.R")
Best track data for the North Atlantic, Gulf of Mexico, and Caribbean sea are available from the NOAA NHC HURDAT database (http://www.nhc.noaa.gov/pastall.shtml). The data are provided in 6-hourly intervals. Wind speed are given in knots. More information can be found in Jarvinen et al. (1984).
Interpolate raw best-track to 1-hr intervals (Elsner and Jagger 2013), subract 60% of forward speed to obtain maximum rotational velocity):
begin2 = 1967
end2 = 2010
WindOBS2.df = best.use[best.use$Yr >= begin & best.use$Yr <= end, ]
WindOBS2.df = subset(WindOBS2.df, Type == "*")
WindOBS2.df$WmaxS = WindOBS2.df$WmaxS * 0.5144
WindOBS2.df$maguv = WindOBS2.df$maguv * 0.5144
WindOBS2.df$WmaxS = WindOBS2.df$WmaxS - WindOBS2.df$maguv * 0.6
evOnly = WindOBS2.df$hr/2 == trunc(WindOBS2.df$hr/2)
WindOBS2.df = WindOBS2.df[evOnly, ]
WindOBS2.df = subset(WindOBS2.df, M == 0)
The simulated track data are available in Matlab format from ftp://texmex.mit.edu/pub/emanuel/ForJim/ncep850b.zip.
###Note: The file was replaced on 8/21/12 at 2:06:00 PM with ncep850b_rev.zip. Email from Kerry: ''I placed a revised zip file in the same directory as before; this is identical to the old one but contains an additional array, vnetstore, that contains the peak surface winds in knots that includes the translation effect plus our baroclinic correction. The latter is unpublished, but basically adds a fraction of the background wind shear to the surface wind speed, based on isallobaric arguments. I am not sure how this might affect sensitivity.''
# tmp = download.file(url =
# 'ftp://texmex.mit.edu/pub/emanuel/ForJim/ncep850b_rev.zip', destfile =
# 'ncep.zip', mode='wb') ncep = unzip('ncep.zip')
ncep = readMat("ncep850b.mat")
# names(ncep) str(ncep)
Track information is stored in arrays with the first index representing the cyclone number and the second indicating the center fix along the track. Additional vectors record the cyclone year and the annual frequency. The maximum rotational velocity is given in units of knots.
Reformat the data to match best-track:
## Sid SYear Mo Da hr lon lat WmaxS P YrRate
## 1 1 1980 9 24 20 -68.50 17.76 8.315 1004 9.4
## 2 1 1980 9 24 22 -68.66 17.81 8.229 1004 9.4
## 3 1 1980 9 25 0 -68.82 17.85 8.466 1004 9.4
## 4 1 1980 9 25 2 -68.98 17.90 8.855 1004 9.4
## 5 1 1980 9 25 4 -69.13 17.95 9.360 1004 9.4
## 6 1 1980 9 25 6 -69.29 18.01 9.930 1003 9.4
## [1] 6200
Figure 1: Comparison of wind speed densities (top) and cumulative distributions (bottom). Top: Frequency of 2-hour hurricane records. Rotational intensity shown in m/s. Bottom: Cumulative Distribution of 2-hour hurricane records
Wind = c(WindOBS.df$WmaxS, WindSIM.df$WmaxS)
Label = c(rep("Best-track", length(WindOBS.df$WmaxS)), rep("Simulated", length(WindSIM.df$WmaxS)))
WindALL.df = data.frame(Label = Label, Wind = Wind)
p = ggplot(WindALL.df, aes(Wind)) + facet_grid(~Label) + geom_density(fill = "grey",
alpha = 0.2) + theme_bw()
p = p + xlab(expression(paste("Rotational Wind Speed [m s", {
}^-1, "]"))) + ylab("Relative Frequency\n Two Hour Records (1980-2010)")
p = p + ggtitle("a") + theme(plot.title = element_text(hjust = 0))
cumdist = ddply(WindALL.df, .(Label), summarize, Wind = unique(Wind), ecdf = ecdf(Wind)(unique(Wind)))
p1 = ggplot(cumdist, aes(Wind, ecdf, color = Label)) + geom_step() + theme_bw()
p1 = p1 + xlab(expression(paste("Rotational Wind Speed [m s", {
}^-1, "]"))) + ylab("Cumulative distribution of\n 2-hour records (1980-2010)")
p1 = p1 + ggtitle("b") + theme(plot.title = element_text(hjust = 0))
multiplot(p, p1, cols = 1)
Calculate skewness and kurtosis for best-track and simulated data:
sOBS = skew(WindALL.df$Wind[WindALL.df$Label == "Best-track"])
sOBS
## [1] 0.9836
sSIM = skew(WindALL.df$Wind[WindALL.df$Label == "Simulated"])
sSIM
## [1] 1.114
kOBS = kurtosi(WindALL.df$Wind[WindALL.df$Label == "Best-track"])
kOBS
## [1] 0.5181
kSIM = kurtosi(WindALL.df$Wind[WindALL.df$Label == "Simulated"])
kSIM
## [1] 1.604
SST data available from www.esrl.noaa.gov/psd/data/gridded/data.noaa.ersst.html in netCDF format.
Extract SST data and assign an lambert conformal conic projection (lcc).
Specify track points using a lcc projection:
Create a hexagonal tessellation of the basin bounding the set of best-track hurricanes. A rectangular domain encompassing the set of observed hurricane locations is gridded into equal-area hexagons (Elsner et al. 2012).
## [1] 833379
Overlay the hexagon tesselation on the observed and simulated track points to obtain per grid TC count and per grid max wind speed.
hexidOBS = over(x = WindOBS.sdf, y = hpg)
WindOBS.sdf$hexid = hexidOBS
hexidSIM = over(x = WindSIM.sdf, y = hpg)
WindSIM.sdf$hexid = hexidSIM
hexidOBS2 = over(x = WindOBS2.sdf, y = hpg)
WindOBS2.sdf$hexid = hexidOBS2
# hexidOBS3 = over(x=WindOBS3.sdf, y=hpg) WindOBS3.sdf$hexid = hexidOBS3
How many simulated storm hours are outside the tessellation boundaries?
sum(is.na(hexidOBS))
## [1] 0
sum(is.na(hexidSIM))
## [1] 11437
sum(is.na(hexidSIM))/dim(WindSIM.df)[1] * 100 #1.56% of hexagons outside (SIM)
## [1] 1.577
Compute per hexagon statistics.
OBSstats = ddply(WindOBS.sdf@data, .(hexid, Sid), summarize, nhr = length(Sid),
WmaxS = max(WmaxS))
OBSstats2 = ddply(OBSstats, .(hexid), summarize, countOBS = length(Sid), WmaxOBS = max(WmaxS))
SIMstats = ddply(WindSIM.sdf@data, .(hexid, Sid), summarize, nhr = length(Sid),
WmaxS = max(WmaxS))
SIMstats2 = ddply(SIMstats, .(hexid), summarize, countSIM = length(Sid), WmaxSIM = max(WmaxS))
OBS2stats = ddply(WindOBS2.sdf@data, .(hexid, Sid), summarize, nhr = length(Sid),
WmaxS = max(WmaxS))
OBS2stats2 = ddply(OBS2stats, .(hexid), summarize, countOBS2 = length(Sid),
WmaxOBS2 = max(WmaxS))
How many simulated storms are outside the grid?
sum(is.na(SIMstats))
## [1] 615
sum(is.na(SIMstats))/length(unique(WindSIM.df$Sid)) * 100
## [1] 9.919
Merge the per hexagon observed and simulated data. Remove grids with countOBS less than 15. Remove grids over land.
hexData = merge(OBSstats2, SIMstats2)
hexData = merge(hexData, OBS2stats2)
# hexData = merge(hexData, OBS3stats2)
hexData = subset(hexData, countOBS >= 20)
# hexData = subset(hexData, hexid !=49 & hexid !=50) #grids over land
row.names(hexData) = paste("ID", hexData$hexid, sep = "")
Correlations.
cor(hexData$countOBS, hexData$countSIM)
## [1] 0.6888
cor(hexData$WmaxOBS, hexData$WmaxSIM)
## [1] 0.7739
Create a spatial polygons data frame.
hspdf = SpatialPolygonsDataFrame(hpg[row.names(hexData)], hexData)
Compute grid average SST.
avgsst = over(x = hpg[row.names(hexData)], y = SST.sdf, fn = mean)
avgsst2 = over(x = hpg[row.names(hexData)], y = SST.sdf2, fn = mean)
hspdf$avgsst = avgsst$avgsst
hspdf$avgsst2 = avgsst2$avgsst2
hspdf = hspdf[!hspdf$hexid == 9, ]
Obtain map boundaries for plotting.
Figure 2: Frequency and intensity maps. The shading gives per hexagon TC counts and the red text denotes the per hexagon maximum wind speed
centers = coordinates(hspdf)
l2 = list("sp.text", centers, round(hspdf$WmaxOBS), col = "red", cex = 1.1)
cr = brewer.pal(6, "Blues")
rng = seq(20, 80, 10)
p0 = spplot(hspdf, "countOBS", col = "white", col.regions = cr, sp.layout = list(l1,
l2), at = rng, colorkey = list(space = "bottom", labels = paste(rng)), sub = expression("Observed number of hurricanes (1980-2010)"))
l2 = list("sp.text", centers, round(hspdf$WmaxSIM), col = "red", cex = 1.1)
cr = brewer.pal(6, "Blues")
rng = seq(400, 2300, 300)
p1 = spplot(hspdf, "countSIM", col = "white", col.regions = cr, sp.layout = list(l1,
l2), at = rng, colorkey = list(space = "bottom", labels = paste(rng)), sub = expression("Simulated number of hurricanes (1980-2010)"))
p0 = update(p0, main = textGrob("a", x = unit(0.05, "npc")), par.settings = list(fontsize = list(text = 15)))
p1 = update(p1, main = textGrob("b", x = unit(0.05, "npc")), par.settings = list(fontsize = list(text = 15)))
plot(p0, split = c(1, 1, 1, 2), more = TRUE)
plot(p1, split = c(1, 2, 1, 2), more = FALSE)
NOTE: maximum intensities (red text) for the simulated TCs are significantly higher in the eastern portion of the basin. Do these storms simply intensify faster? See gen and max plots later…
Figure 3: Relative risk of a simulated hurricane. Factor by which the simulated relative TC frequency exceeds that of the observed relative TC frequency.
hspdf$IntDiff = hspdf$WmaxSIM - hspdf$WmaxOBS
hspdf$RelativeRisk = (hspdf$countSIM/sum(hspdf$countSIM))/(hspdf$countOBS/sum(hspdf$countOBS))
cr = brewer.pal(7, "RdBu")
cr2 = c(rev(cr[1:5]), rev(cr[6:11]))
rng = seq(0.4, 1.8, 0.2)
p3 = spplot(hspdf, "RelativeRisk", col = "white", col.regions = cr, sp.layout = list(l1),
at = rng, colorkey = list(space = "bottom", labels = paste(rng)), sub = expression("Relative Risk (Simulation/Best-track)"))
p3
Map SST throughout the basin: Figure 4: Map of SST averaged over the grids.
clsSST = brewer.pal(7, "Reds")
l3 = list("sp.text", centers, id, col = "blue", cex = 1.1)
rng = seq(9, 30, 3)
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
Compute correlation between average SST and hurricane counts:
Bottom line: very little correspondence between storm counts and SST for Observed and for Simulated storms.
Probably not used (give cor.test results instead…maybe table it?)
Calculate correlation between SST and per grid maximum wind speed:
Plot WmaxOBS vs. WmaxSIM:
First, obtain LI using Poisson-EV model:
hspdf$thresh = thresh
hspdf$rate = rate
hspdf$lambda = lambda
hspdf$sigma = sigma
hspdf$xi = xi
hspdf$theo = theo
range(hspdf@data$theo)
## [1] 32.34 103.01
head(hspdf@data)
## hexid countOBS WmaxOBS countSIM WmaxSIM countOBS2 WmaxOBS2 avgsst
## ID5 5 20 67.43 606 84.02 20 67.43 28.67
## ID6 6 33 59.00 1112 79.04 36 59.00 28.81
## ID7 7 37 64.84 1282 80.65 41 64.84 28.64
## ID8 8 32 57.15 1021 87.09 37 57.15 28.43
## ID15 15 46 78.09 1100 86.55 47 78.09 28.95
## ID16 16 44 73.19 1617 92.70 47 73.19 28.91
## avgsst2 IntDiff RelativeRisk thresh rate lambda sigma xi theo
## ID5 28.56 16.592 1.1878 30 0.6452 0.1290 72.16 -1.928 67.43
## ID6 28.73 20.037 1.3210 30 1.0645 0.1935 39.58 -1.365 59.00
## ID7 28.55 15.805 1.3583 30 1.1935 0.1935 42.47 -1.219 64.84
## ID8 28.34 29.938 1.2508 30 1.0323 0.2258 36.08 -1.329 57.15
## ID15 28.87 8.461 0.9374 30 1.4839 0.7419 49.97 -1.039 78.09
## ID16 28.83 19.508 1.4406 30 1.4194 0.5484 51.10 -1.183 73.19
sum(hspdf$xi < -1) # 15/31 == Wmax
## [1] 15
sum(hspdf$xi < 0) #31 out of 31 asymptote to LI
## [1] 31
# test = hspdf@data test[test$xi > 0,] ## Hexid for xi > 0 = 9 (farthest
# SE)
Calculate sensitivity using LI from model:
dat = as.data.frame(hspdf)
dat$northing = coordinates(hspdf)[, 2]
lmOBS.theo = summary(lm(theo ~ avgsst, weights = countOBS, data = dat[dat$avgsst >=
25, ]))
sensOBS.theo = coef(lm(theo ~ avgsst, weights = countOBS, data = dat[dat$avgsst >=
25, ]))[2]
Extract a data frame of LI and SST and compute the sensitivity, Use weights=count for weighted regression, test if northing is an important predictor:
Plot SST vs. LI (same as previous plot…): Figure 6: Sensitivity of limiting intensity to SST. The slope of the red line is the sensitivity of limiting intensity to spatial variation in SST.
dat = subset(dat, avgsst >= 25)
avgsst = rep(dat$avgsst, 2)
type = c(rep("Best-track", nrow(dat)), rep("Simulated", nrow(dat)))
LI = c(dat$WmaxOBS, dat$WmaxSIM)
count = c(dat$countOBS, dat$countSIM)
dat.df = data.frame(avgsst = avgsst, type = type, LI = LI, count = count)
xlabel = expression(paste("Sea surface temperature [", degree, "C]"))
ylabel = expression(paste("Limiting intensity [m ", s^-1, "]"))
ggplot(dat.df, aes(x = avgsst, y = LI)) + geom_point() + geom_smooth(method = "lm",
aes(weight = count), color = "red") + facet_wrap(~type) + xlab(xlabel) +
ylab(ylabel) + theme_bw()
d = hspdf@data[hspdf$avgsst >= 25, ]
d2 = subset(d, WmaxOBS < 40)
d2 #hexid = 11 & 23
## hexid countOBS WmaxOBS countSIM WmaxSIM countOBS2 WmaxOBS2 avgsst
## ID22 22 34 34.12 454 59.91 35 34.12 27.87
## avgsst2 IntDiff RelativeRisk thresh rate lambda sigma xi theo
## ID22 27.75 25.79 0.5235 30 1.097 0.06452 5.18 -1.259 34.12
Estimate Moran's I on the model residuals.
hspdf2 = hspdf[hspdf$avgsst >= 25, ]
hexnb2 = poly2nb(hspdf2)
wts2 = nb2listw(hexnb2, style = "W")
stvty.lmOBS = lm(WmaxOBS ~ avgsst, weights = countOBS, data = hspdf2)
stvty.lmSIM = lm(WmaxSIM ~ avgsst, weights = countSIM, data = hspdf2)
lm.morantest(stvty.lmOBS, wts2)
##
## Global Moran's I for regression residuals
##
## data:
## model: lm(formula = WmaxOBS ~ avgsst, data = hspdf2, weights =
## countOBS)
## weights: wts2
##
## Moran I statistic standard deviate = 2.417, p-value = 0.007819
## alternative hypothesis: greater
## sample estimates:
## Observed Moran's I Expectation Variance
## 0.25042 -0.07095 0.01767
lm.morantest(stvty.lmSIM, wts2)
##
## Global Moran's I for regression residuals
##
## data:
## model: lm(formula = WmaxSIM ~ avgsst, data = hspdf2, weights =
## countSIM)
## weights: wts2
##
## Moran I statistic standard deviate = 4.427, p-value = 4.77e-06
## alternative hypothesis: greater
## sample estimates:
## Observed Moran's I Expectation Variance
## 0.51948 -0.07020 0.01774
Perform a geographically weighted regression:
bwOBS = gwr.sel(WmaxOBS ~ avgsst, weights = countOBS, data = hspdf2)
## Bandwidth: 3156792 CV score: 2363777
## Bandwidth: 5102697 CV score: 2578424
## Bandwidth: 1954156 CV score: 2137657
## Bandwidth: 1210887 CV score: 2110550
## Bandwidth: 1348390 CV score: 2086513
## Bandwidth: 1530241 CV score: 2082566
## Bandwidth: 1461950 CV score: 2081044
## Bandwidth: 1467327 CV score: 2081046
## Bandwidth: 1463882 CV score: 2081042
## Bandwidth: 1463928 CV score: 2081042
## Bandwidth: 1463925 CV score: 2081042
## Bandwidth: 1463925 CV score: 2081042
## Bandwidth: 1463925 CV score: 2081042
## Bandwidth: 1463926 CV score: 2081042
## Bandwidth: 1463925 CV score: 2081042
## Bandwidth: 1463925 CV score: 2081042
## Bandwidth: 1463925 CV score: 2081042
## Bandwidth: 1463925 CV score: 2081042
## Bandwidth: 1463925 CV score: 2081042
bwSIM = gwr.sel(WmaxSIM ~ avgsst, weights = countSIM, data = hspdf2)
## Bandwidth: 3156792 CV score: 1.293e+09
## Bandwidth: 5102697 CV score: 1.441e+09
## Bandwidth: 1954156 CV score: 1.064e+09
## Bandwidth: 1210887 CV score: 809190798
## Bandwidth: 751521 CV score: 637824374
## Bandwidth: 467617 CV score: 623658127
## Bandwidth: 552185 CV score: 6.08e+08
## Bandwidth: 588362 CV score: 607598337
## Bandwidth: 574608 CV score: 607310868
## Bandwidth: 574441 CV score: 607310714
## Bandwidth: 574202 CV score: 607310636
## Bandwidth: 574210 CV score: 607310635
## Bandwidth: 574210 CV score: 607310635
## Bandwidth: 574210 CV score: 607310635
## Bandwidth: 574210 CV score: 607310635
## Bandwidth: 574210 CV score: 607310635
## Bandwidth: 574210 CV score: 607310635
gwrOBS = gwr(WmaxOBS ~ avgsst, weights = countOBS, data = hspdf2, bandwidth = bwOBS)
gwrSIM = gwr(WmaxSIM ~ avgsst, weights = countSIM, data = hspdf2, bandwidth = bwSIM)
rSensOBSSIM = cor(gwrOBS$SDF$avgsst, gwrSIM$SDF$avgsst)
range(gwrOBS$SDF$avgsst)
## [1] -0.02293 9.24954
range(gwrSIM$SDF$avgsst)
## [1] -5.712 20.742
Figure t: Geographically weighted regression of LI on SST. Colors represent the sensitivity in m s1 K1. Gray represents sensitivity values less than 1
gwrOBS$SDF$avgsst[gwrOBS$SDF$avgsst < 0] = -1
gwrSIM$SDF$avgsst[gwrSIM$SDF$avgsst < 0] = -1
cls = brewer.pal(6, "Oranges")
cls = c("gray", cls)
rng = seq(-4, 24, 4)
p0 = spplot(gwrOBS$SDF, "avgsst", col = "white", col.regions = cls, sp.layout = list(l1),
at = rng, colorkey = list(space = "bottom", labels = c("", rng[2:length(rng)])),
sub = expression(paste("Sensitivity (observed) [m ", s^-1, K^-1, "]")))
p1 = spplot(gwrSIM$SDF, "avgsst", col = "white", col.regions = cls, sp.layout = list(l1),
at = rng, colorkey = list(space = "bottom", labels = c("", rng[2:length(rng)])),
sub = expression(paste("Sensitivity (simulated) [m ", s^-1, K^-1, "]")))
p0 = update(p0, main = textGrob("a", x = unit(0.05, "npc")), par.settings = list(fontsize = list(text = 15)))
p1 = update(p1, main = textGrob("b", x = unit(0.05, "npc")), par.settings = list(fontsize = list(text = 15)))
plot(p0, split = c(1, 1, 1, 2), more = TRUE)
plot(p1, split = c(1, 2, 1, 2), more = FALSE)
rng = seq(-2, 10, 2)
cls = brewer.pal(5, "Oranges")
cls = c("gray", cls)
spplot(gwrOBS$SDF, "avgsst", col = "white", col.regions = cls, sp.layout = list(l1),
at = rng, colorkey = list(space = "bottom", labels = c("", rng[2:length(rng)])),
sub = expression(paste("Sensitivity (observed) [m ", s^-1, K^-1, "]")))
Compare observed and simulated sensitivities:
We test these results by generating several different grids:
set.seed(2013)
hpt.i1 = spsample(WindOBS.sdf, type = "hexagonal", n = 120, bb = bb)
hpg.i1 = HexPoints2SpatialPolygons(hpt.i1)
hpg.i1@polygons[[1]]@area * 10^-6 # square
## [1] 833379
hpt.i2 = spsample(WindOBS.sdf, type = "hexagonal", n = 120, bb = bb)
hpg.i2 = HexPoints2SpatialPolygons(hpt.i2)
hpg.i2@polygons[[1]]@area * 10^-6 # square km
## [1] 833379
hpt.i3 = spsample(WindOBS.sdf, type = "hexagonal", n = 120, bb = bb)
hpg.i3 = HexPoints2SpatialPolygons(hpt.i3)
hpg.i3@polygons[[1]]@area * 10^-6 # square km
## [1] 833379
hpt.i4 = spsample(WindOBS.sdf, type = "hexagonal", n = 120, bb = bb)
hpg.i4 = HexPoints2SpatialPolygons(hpt.i4)
hpg.i4@polygons[[1]]@area * 10^-6 # square km
## [1] 833379
hpt.i5 = spsample(WindOBS.sdf, type = "hexagonal", n = 120, bb = bb)
hpg.i5 = HexPoints2SpatialPolygons(hpt.i5)
hpg.i5@polygons[[1]]@area * 10^-6 # square km
## [1] 833379
hspdf.i1 = get.hspdf(WindOBS.sdf, WindSIM.sdf, hpg.i1, SST.sdf)
hspdf.i2 = get.hspdf(WindOBS.sdf, WindSIM.sdf, hpg.i2, SST.sdf)
hspdf.i3 = get.hspdf(WindOBS.sdf, WindSIM.sdf, hpg.i3, SST.sdf)
hspdf.i4 = get.hspdf(WindOBS.sdf, WindSIM.sdf, hpg.i4, SST.sdf)
hspdf.i5 = get.hspdf(WindOBS.sdf, WindSIM.sdf, hpg.i5, SST.sdf)
hspdf.i1 = hspdf.i1[!is.na(hspdf.i1$avgsst), ]
hspdf.i1 = hspdf.i1[hspdf.i1$avgsst >= 25, ]
hspdf.i2 = hspdf.i2[!is.na(hspdf.i2$avgsst), ]
hspdf.i2 = hspdf.i2[hspdf.i2$avgsst >= 25, ]
hspdf.i3 = hspdf.i3[!is.na(hspdf.i3$avgsst), ]
hspdf.i3 = hspdf.i3[hspdf.i3$avgsst >= 25, ]
hspdf.i4 = hspdf.i4[!is.na(hspdf.i4$avgsst), ]
hspdf.i4 = hspdf.i4[hspdf.i4$avgsst >= 25, ]
hspdf.i5 = hspdf.i5[!is.na(hspdf.i5$avgsst), ]
hspdf.i5 = hspdf.i5[hspdf.i5$avgsst >= 25, ]
GWR.i1 = get.gwr(hspdf.i1)
## Bandwidth: 3335771 CV score: 3265831
## Bandwidth: 5392003 CV score: 3668844
## Bandwidth: 2064950 CV score: 2771571
## Bandwidth: 1279540 CV score: 2364031
## Bandwidth: 794130 CV score: 2189149
## Bandwidth: 494129 CV score: 2723981
## Bandwidth: 970815 CV score: 2212746
## Bandwidth: 791539 CV score: 2189500
## Bandwidth: 837977 CV score: 2186987
## Bandwidth: 888717 CV score: 2191916
## Bandwidth: 829335 CV score: 2186891
## Bandwidth: 830435 CV score: 2186890
## Bandwidth: 830195 CV score: 2186890
## Bandwidth: 830190 CV score: 2186890
## Bandwidth: 830190 CV score: 2186890
## Bandwidth: 830190 CV score: 2186890
## Bandwidth: 830190 CV score: 2186890
## Bandwidth: 830190 CV score: 2186890
## Bandwidth: 830190 CV score: 2186890
## Bandwidth: 3335771 CV score: 1.572e+09
## Bandwidth: 5392003 CV score: 1.715e+09
## Bandwidth: 2064950 CV score: 1.334e+09
## Bandwidth: 1279540 CV score: 1.069e+09
## Bandwidth: 794130 CV score: 851838016
## Bandwidth: 494129 CV score: 974196625
## Bandwidth: 831460 CV score: 867238712
## Bandwidth: 727980 CV score: 830712589
## Bandwidth: 638657 CV score: 833138234
## Bandwidth: 689410 CV score: 825425732
## Bandwidth: 687514 CV score: 825370864
## Bandwidth: 683048 CV score: 825331659
## Bandwidth: 683896 CV score: 825329180
## Bandwidth: 683924 CV score: 825329179
## Bandwidth: 683918 CV score: 825329179
## Bandwidth: 683918 CV score: 825329179
## Bandwidth: 683918 CV score: 825329179
## Bandwidth: 683918 CV score: 825329179
GWR.i2 = get.gwr(hspdf.i2)
## Bandwidth: 3335771 CV score: 3632295
## Bandwidth: 5392003 CV score: 4134528
## Bandwidth: 2064950 CV score: 2992435
## Bandwidth: 1279540 CV score: 2410678
## Bandwidth: 794130 CV score: 2105537
## Bandwidth: 494129 CV score: 2178051
## Bandwidth: 753193 CV score: 2102837
## Bandwidth: 745898 CV score: 2102908
## Bandwidth: 752621 CV score: 2102837
## Bandwidth: 752618 CV score: 2102837
## Bandwidth: 752619 CV score: 2102837
## Bandwidth: 752619 CV score: 2102837
## Bandwidth: 752619 CV score: 2102837
## Bandwidth: 752619 CV score: 2102837
## Bandwidth: 752619 CV score: 2102837
## Bandwidth: 3335771 CV score: 1.633e+09
## Bandwidth: 5392003 CV score: 1.817e+09
## Bandwidth: 2064950 CV score: 1.379e+09
## Bandwidth: 1279540 CV score: 1.148e+09
## Bandwidth: 794130 CV score: 1.038e+09
## Bandwidth: 494129 CV score: 1.06e+09
## Bandwidth: 738317 CV score: 1.038e+09
## Bandwidth: 763695 CV score: 1.038e+09
## Bandwidth: 764757 CV score: 1.038e+09
## Bandwidth: 765109 CV score: 1.038e+09
## Bandwidth: 765088 CV score: 1.038e+09
## Bandwidth: 765088 CV score: 1.038e+09
## Bandwidth: 765088 CV score: 1.038e+09
## Bandwidth: 765088 CV score: 1.038e+09
## Bandwidth: 765088 CV score: 1.038e+09
## Bandwidth: 765088 CV score: 1.038e+09
## Bandwidth: 765088 CV score: 1.038e+09
## Bandwidth: 765088 CV score: 1.038e+09
## Bandwidth: 765088 CV score: 1.038e+09
GWR.i3 = get.gwr(hspdf.i3)
## Bandwidth: 3156792 CV score: 2318072
## Bandwidth: 5102697 CV score: 2661541
## Bandwidth: 1954156 CV score: 1949020
## Bandwidth: 1210887 CV score: 1835163
## Bandwidth: 612753 CV score: 1957485
## Bandwidth: 1295288 CV score: 1833016
## Bandwidth: 1300007 CV score: 1833033
## Bandwidth: 1292096 CV score: 1833013
## Bandwidth: 1292112 CV score: 1833013
## Bandwidth: 1292113 CV score: 1833013
## Bandwidth: 1292113 CV score: 1833013
## Bandwidth: 1292113 CV score: 1833013
## Bandwidth: 1292113 CV score: 1833013
## Bandwidth: 1292113 CV score: 1833013
## Bandwidth: 1292113 CV score: 1833013
## Bandwidth: 1292113 CV score: 1833013
## Bandwidth: 3156792 CV score: 2.018e+09
## Bandwidth: 5102697 CV score: 2.232e+09
## Bandwidth: 1954156 CV score: 1.706e+09
## Bandwidth: 1210887 CV score: 1.43e+09
## Bandwidth: 751521 CV score: 1.245e+09
## Bandwidth: 467617 CV score: 1.17e+09
## Bandwidth: 292155 CV score: 1.17e+09
## Bandwidth: 379664 CV score: 1.17e+09
## Bandwidth: 372426 CV score: 1.17e+09
## Bandwidth: 413259 CV score: 1.17e+09
## Bandwidth: 434022 CV score: 1.17e+09
## Bandwidth: 444916 CV score: 1.17e+09
## Bandwidth: 429468 CV score: 1.17e+09
## Bandwidth: 431372 CV score: 1.17e+09
## Bandwidth: 431294 CV score: 1.17e+09
## Bandwidth: 431330 CV score: 1.17e+09
## Bandwidth: 431330 CV score: 1.17e+09
## Bandwidth: 431331 CV score: 1.17e+09
## Bandwidth: 431331 CV score: 1.17e+09
## Bandwidth: 431330 CV score: 1.17e+09
## Bandwidth: 431330 CV score: 1.17e+09
## Bandwidth: 431330 CV score: 1.17e+09
## Bandwidth: 431330 CV score: 1.17e+09
## Bandwidth: 431330 CV score: 1.17e+09
GWR.i4 = get.gwr(hspdf.i4)
## Bandwidth: 3335771 CV score: 3138665
## Bandwidth: 5392003 CV score: 3627265
## Bandwidth: 2064950 CV score: 2524709
## Bandwidth: 1279540 CV score: 2039808
## Bandwidth: 794130 CV score: 1762502
## Bandwidth: 494129 CV score: 1713943
## Bandwidth: 488874 CV score: 1714395
## Bandwidth: 544499 CV score: 1711081
## Bandwidth: 639850 CV score: 1717506
## Bandwidth: 573287 CV score: 1711193
## Bandwidth: 556351 CV score: 1710948
## Bandwidth: 556704 CV score: 1710948
## Bandwidth: 556861 CV score: 1710948
## Bandwidth: 556854 CV score: 1710948
## Bandwidth: 556854 CV score: 1710948
## Bandwidth: 556854 CV score: 1710948
## Bandwidth: 556857 CV score: 1710948
## Bandwidth: 556855 CV score: 1710948
## Bandwidth: 556855 CV score: 1710948
## Bandwidth: 556854 CV score: 1710948
## Bandwidth: 556854 CV score: 1710948
## Bandwidth: 556854 CV score: 1710948
## Bandwidth: 556854 CV score: 1710948
## Bandwidth: 556854 CV score: 1710948
## Bandwidth: 556854 CV score: 1710948
## Bandwidth: 3335771 CV score: 1.971e+09
## Bandwidth: 5392003 CV score: 2.234e+09
## Bandwidth: 2064950 CV score: 1.608e+09
## Bandwidth: 1279540 CV score: 1.288e+09
## Bandwidth: 794130 CV score: 1.102e+09
## Bandwidth: 494129 CV score: 1.039e+09
## Bandwidth: 178357 CV score: NA
## Warning: NA/Inf replaced by maximum positive value
## Bandwidth: 373515 CV score: 1.028e+09
## Bandwidth: 259922 CV score: 1.027e+09
## Bandwidth: 311360 CV score: 1.027e+09
## Bandwidth: 288028 CV score: 1.027e+09
## Bandwidth: 322490 CV score: 1.027e+09
## Bandwidth: 302937 CV score: 1.027e+09
## Bandwidth: 305159 CV score: 1.027e+09
## Bandwidth: 305823 CV score: 1.027e+09
## Bandwidth: 306085 CV score: 1.027e+09
## Bandwidth: 306041 CV score: 1.027e+09
## Bandwidth: 306040 CV score: 1.027e+09
## Bandwidth: 306040 CV score: 1.027e+09
## Bandwidth: 306039 CV score: 1.027e+09
## Bandwidth: 305957 CV score: 1.027e+09
## Bandwidth: 306008 CV score: 1.027e+09
## Bandwidth: 306027 CV score: 1.027e+09
## Bandwidth: 306035 CV score: 1.027e+09
## Bandwidth: 306038 CV score: 1.027e+09
## Bandwidth: 306039 CV score: 1.027e+09
## Bandwidth: 306039 CV score: 1.027e+09
## Bandwidth: 306039 CV score: 1.027e+09
## Bandwidth: 306039 CV score: 1.027e+09
## Bandwidth: 306039 CV score: 1.027e+09
## Bandwidth: 306039 CV score: 1.027e+09
## Bandwidth: 306039 CV score: 1.027e+09
## Bandwidth: 306039 CV score: 1.027e+09
## Bandwidth: 306039 CV score: 1.027e+09
## Bandwidth: 306039 CV score: 1.027e+09
## Bandwidth: 306039 CV score: 1.027e+09
GWR.i5 = get.gwr(hspdf.i5)
## Bandwidth: 3335771 CV score: 3193933
## Bandwidth: 5392003 CV score: 3692575
## Bandwidth: 2064950 CV score: 2565233
## Bandwidth: 1279540 CV score: 2081537
## Bandwidth: 794130 CV score: 1811527
## Bandwidth: 494129 CV score: 1757531
## Bandwidth: 456276 CV score: 1760397
## Bandwidth: 525221 CV score: 1755714
## Bandwidth: 627935 CV score: 1760366
## Bandwidth: 547367 CV score: 1755082
## Bandwidth: 561691 CV score: 1755078
## Bandwidth: 554695 CV score: 1755036
## Bandwidth: 554679 CV score: 1755036
## Bandwidth: 554754 CV score: 1755036
## Bandwidth: 554755 CV score: 1755036
## Bandwidth: 554755 CV score: 1755036
## Bandwidth: 554755 CV score: 1755036
## Bandwidth: 554755 CV score: 1755036
## Bandwidth: 554755 CV score: 1755036
## Bandwidth: 554755 CV score: 1755036
## Bandwidth: 554755 CV score: 1755036
## Bandwidth: 3335771 CV score: 1.991e+09
## Bandwidth: 5392003 CV score: 2.255e+09
## Bandwidth: 2064950 CV score: 1.628e+09
## Bandwidth: 1279540 CV score: 1.313e+09
## Bandwidth: 794130 CV score: 1.14e+09
## Bandwidth: 494129 CV score: 1.083e+09
## Bandwidth: 195875 CV score: 1.075e+09
## Bandwidth: 300840 CV score: 1.075e+09
## Bandwidth: 250094 CV score: 1.075e+09
## Bandwidth: 374670 CV score: 1.076e+09
## Bandwidth: 329041 CV score: 1.075e+09
## Bandwidth: 346470 CV score: 1.075e+09
## Bandwidth: 329481 CV score: 1.075e+09
## Bandwidth: 332564 CV score: 1.075e+09
## Bandwidth: 337876 CV score: 1.075e+09
## Bandwidth: 333496 CV score: 1.075e+09
## Bandwidth: 333305 CV score: 1.075e+09
## Bandwidth: 333288 CV score: 1.075e+09
## Bandwidth: 333289 CV score: 1.075e+09
## Bandwidth: 333289 CV score: 1.075e+09
## Bandwidth: 333295 CV score: 1.075e+09
## Bandwidth: 333292 CV score: 1.075e+09
## Bandwidth: 333290 CV score: 1.075e+09
## Bandwidth: 333290 CV score: 1.075e+09
## Bandwidth: 333289 CV score: 1.075e+09
## Bandwidth: 333289 CV score: 1.075e+09
## Bandwidth: 333289 CV score: 1.075e+09
## Bandwidth: 333289 CV score: 1.075e+09
## Bandwidth: 333289 CV score: 1.075e+09
## Bandwidth: 333289 CV score: 1.075e+09
## Bandwidth: 333289 CV score: 1.075e+09
gwrOBS$SDF$avgsst[gwrOBS$SDF$avgsst < 0] = 0
gwrSIM$SDF$avgsst[gwrSIM$SDF$avgsst < 0] = 0
gwrOBS$SDF$diff = gwrOBS$SDF$avgsst - gwrSIM$SDF$avgsst
Conclude: No need to use spatial regression. This is confirmed with the small Moran's I on the model residuals.
Modeled sensitivity.
Choose a spatial error model.
Note the increase is standard error for the sensitivity estimate. 6 +/- 2.12 * 1.96. The estimate is 6.0 (1.8, 10.2) m/s/K (95% CI). This compares with a sensitivity of 7.8 (4.2, 11.3).
Create grids of varying hexagon size:
bb = bbox(WindOBS.sdf) * 1.2
bb[2, 1] = bbox(WindOBS.sdf)[2, 1] * 0.7
set.seed(3042)
hpt3 = spsample(WindOBS.sdf, type = "hexagonal", n = 60, bb = bb, offset = c(0.5,
0.5))
hpg3 = HexPoints2SpatialPolygons(hpt3)
area3 = hpg3@polygons[[1]]@area * 10^-6 # square km
area3 #1.67E6
## [1] 1666757
set.seed(3042)
hpt4 = spsample(WindOBS.sdf, type = "hexagonal", n = 90, bb = bb, offset = c(0.5,
0.5))
hpg4 = HexPoints2SpatialPolygons(hpt4)
area4 = hpg4@polygons[[1]]@area * 10^-6 # square km
area4 #1.11E6
## [1] 1111171
set.seed(3042)
hpt5 = spsample(WindOBS.sdf, type = "hexagonal", n = 200, bb = bb, offset = c(0.5,
0.5))
hpg5 = HexPoints2SpatialPolygons(hpt5)
area5 = hpg5@polygons[[1]]@area * 10^-6 # square km
area5 #5.00E5
## [1] 5e+05
set.seed(3042)
hpt6 = spsample(WindOBS.sdf, type = "hexagonal", n = 240, bb = bb, offset = c(0.5,
0.5))
hpg6 = HexPoints2SpatialPolygons(hpt6)
area6 = hpg6@polygons[[1]]@area * 10^-6 # square km
area6 #4.17E5
## [1] 416689
The function “get.hspdf” does everything the “overlayWindpointsonHexagons” chunk does for each of the grids created above:
hspdf3 = get.hspdf(WindOBS.sdf, WindSIM.sdf, hpg3, SST.sdf)
hspdf4 = get.hspdf(WindOBS.sdf, WindSIM.sdf, hpg4, SST.sdf)
hspdf5 = get.hspdf(WindOBS.sdf, WindSIM.sdf, hpg5, SST.sdf)
hspdf6 = get.hspdf(WindOBS.sdf, WindSIM.sdf, hpg6, SST.sdf)
Compute correlations:
rsstOBS3 = cor.test(hspdf3$avgsst[hspdf3$avgsst >= 25], hspdf3$WmaxOBS[hspdf3$avgsst >=
25], use = "complete") #0.517
rsstSIM3 = cor.test(hspdf3$avgsst[hspdf3$avgsst >= 25], hspdf3$WmaxSIM[hspdf3$avgsst >=
25], use = "complete") #0.548
rsstOBS4 = cor.test(hspdf4$avgsst[hspdf4$avgsst >= 25], hspdf4$WmaxOBS[hspdf4$avgsst >=
25], use = "complete") #0.699
rsstSIM4 = cor.test(hspdf4$avgsst[hspdf4$avgsst >= 25], hspdf4$WmaxSIM[hspdf4$avgsst >=
25], use = "complete") #0.715
rsstOBS5 = cor.test(hspdf5$avgsst[hspdf5$avgsst >= 25], hspdf5$WmaxOBS[hspdf5$avgsst >=
25], use = "complete") #0.562
rsstSIM5 = cor.test(hspdf5$avgsst[hspdf5$avgsst >= 25], hspdf5$WmaxSIM[hspdf5$avgsst >=
25], use = "complete") #0.614
rsstOBS6 = cor.test(hspdf6$avgsst[hspdf6$avgsst >= 25], hspdf6$WmaxOBS[hspdf6$avgsst >=
25], use = "complete") #0.652
rsstSIM6 = cor.test(hspdf6$avgsst[hspdf6$avgsst >= 25], hspdf6$WmaxSIM[hspdf6$avgsst >=
25], use = "complete") #0.560
Compute Sensitivity for different grid sizes:
dat3 = hspdf3@data
dat4 = hspdf4@data
dat5 = hspdf5@data
dat6 = hspdf6@data
dat3$northing = coordinates(hspdf3)[, 2]
sensO3 = summary(lm(WmaxOBS ~ avgsst, weights = countOBS, data = dat3[dat3$avgsst >=
25, ])) ## 7.76 +/- 2.45
sensS3 = summary(lm(WmaxSIM ~ avgsst, weights = countOBS, data = dat3[dat3$avgsst >=
25, ])) ## 6.93.6 +/- 1.90
dat4$northing = coordinates(hspdf4)[, 2]
sensO4 = summary(lm(WmaxOBS ~ avgsst, weights = countOBS, data = dat4[dat4$avgsst >=
25, ])) ## 8.86 +/- 2.82
sensS4 = summary(lm(WmaxSIM ~ avgsst, weights = countOBS, data = dat4[dat4$avgsst >=
25, ])) ## 10.6 +/- 1.92
dat5$northing = coordinates(hspdf5)[, 2]
sensO5 = summary(lm(WmaxOBS ~ avgsst, weights = countOBS, data = dat5[dat5$avgsst >=
25, ])) ## 7.76 +/- 1.54
sensS5 = summary(lm(WmaxSIM ~ avgsst, weights = countOBS, data = dat5[dat5$avgsst >=
25, ])) ## 7.02 +/- 1.40
dat6$northing = coordinates(hspdf6)[, 2]
sensO6 = summary(lm(WmaxOBS ~ avgsst, weights = countOBS, data = dat6[dat6$avgsst >=
25, ])) ## 7.88 +/- 1.46
sensS6 = summary(lm(WmaxSIM ~ avgsst, weights = countOBS, data = dat6[dat6$avgsst >=
25, ])) ## 6.33 +/- 1.34