date()
## [1] "Mon Mar 18 10:49:16 2013"
setwd("C:/Users/Sarah/Dropbox/GCMTrackData/Sarah")
# setwd('~/Dropbox/GCMTrackData/Sarah')
require(sp)
require(rgdal)
require(maps)
require(maptools)
require(colorRamps)
require(grid)
require(RColorBrewer)
require(verification)
require(ggplot2)
require(plyr)
source("track_hexagons.R")
Of broad scientific and public interest is the reliability of general circulation models (GCMs) to simulate future regional and local tropical cyclone (TC) occurrences. Atmospheric GCMs are now able to generate vortices resembling actual tropical cyclones (TCs), but questions remain about their fidelity to actual TCs. Here the authors demonstrate a spatial lattice approach for comparing actual with simulated TC occurrences regionally using actual TCs from the IBTrACS data set and GCM-generated TCs from the GFDL HiRAM and FSU-COAPS models over the common period 1982–2008. Results show that the spatial distribution of TCs generated by the GFDL model compare well with observations globally, although there are areas of over and under prediction, particularly in parts of the Pacific. Difference maps using the spatial lattice highlight these discrepancies. Additionally, comparisons focusing on the North Atlantic basin are made using both models. Results confirm a large area of over prediction by the FSU-COAPS model in the south-central portion of the basin. Relevant to projections of future U.S. hurricane activity is the fact that both models under predict TC activity in the Gulf of Mexico.
Load IBTrACS data, subset on data from years 1982–2008, convert wind speeds from knots to ms-1, subset on wind speeds greater than 17 ms-1 (tropical storm strength): (Note: use world.interp.v02r01B.R - includes basins defined in “define_basins.R”)
# con = url('http://myweb.fsu.edu/jelsner/world.interp.v02r01.R')
# load(con) load('G:/Research/world.interp.v02r01.R')
load("G:/Research/world.interp.v02r01B.R")
# load('/Users/jelsner/Desktop/Misc Data/world.interp.v02r01.R')
ib.df = subset(world.interp, Yr > 1981 & Yr < 2009)
ib.df$Wmax = ib.df$Wmax * 0.5144
mean(ib.df$Wmax)
## [1] 22.87
# percentile = length(ib.df$Wmax[ib.df$Wmax <= 17])/length(ib.df$Wmax)
ib.df = subset(ib.df, Wmax >= 17)
ch.i = nrow(ib.df)
rm(world.interp)
Create spatial data frame, assign planar coordinate system with a mollweide projection:
ib.sdf = ib.df[, c("name", "lon", "lat", "Sid", "Sn", "Yr", "Mo", "Da", "Wmax",
"DWmaxDt", "basins")]
coordinates(ib.sdf) = c("lon", "lat")
ll = "+proj=longlat +ellps=WGS84"
proj4string(ib.sdf) = CRS(ll)
moll = "+proj=moll +lon_0=150 +datum=WGS84 +units=km"
ib.sdf = spTransform(ib.sdf, CRS(moll))
Set up hexagonal grid:
set.seed(2032)
hpt.i = spsample(ib.sdf, type = "hexagonal", n = 1000, bb = bbox(ib.sdf) * 1.08)
hpg.i = HexPoints2SpatialPolygons(hpt.i)
Extract the number of hexagons and the area of each hexagon (in square km):
np.i = length(hpg.i@polygons)
area.i = hpg.i@polygons[[1]]@area
area.i
## [1] 727859
To obtain per hexagon cyclone counts, we sum the number of unique storm IDs (Sid) occuring in each hexagon.
Compute per hexagon cyclone counts for the Gulf of Mexico in 2005 as an example to demonstrate the spatial lattice framework:
ib.gulf = subset(ib.df, Yr == 2005)
ib.gulf = subset(ib.gulf, lat > 20 & lat < 45)
ib.gulf = subset(ib.gulf, lon > -100 & lon < -80)
ib.Sgulf = ib.gulf[, c("lon", "lat", "Sid", "Sn", "Yr", "Mo", "Da", "Wmax",
"DWmaxDt")]
coordinates(ib.Sgulf) = c("lon", "lat")
proj4string(ib.Sgulf) = CRS(ll)
moll2 = "+proj=moll +lon_0=340 +datum=WGS84"
# lcc = '+proj=lcc +lat_1=30 +lat_2=30 +lon_0=-85 +datum=WGS84'
ib.Sgulf = spTransform(ib.Sgulf, CRS(moll2))
# ib.Sgulf = spTransform(ib.Sgulf, CRS(lcc))
hpt.gulf = spsample(ib.Sgulf, type = "hexagonal", n = 250, offset = c(1, -1))
hpg.gulf = HexPoints2SpatialPolygons(hpt.gulf)
gulfHSPDF = gulf.counts(ib.Sgulf, hpg.gulf)
l1.gulf = list("sp.points", ib.Sgulf, pch = 20, col = "dark gray", cex = 1)
# require(mapdata)
cl.i = map("world", plot = FALSE, wrap = TRUE)
# cl.i = map('worldHires', xlim=c(-110, -70), ylim=c(15, 50), plot=FALSE,
# resolution=.7)
clp.i = map2SpatialLines(cl.i, proj4string = CRS(ll))
clp.gulf = spTransform(clp.i, CRS(moll2))
# clp.gulf = spTransform(clp.i, CRS(lcc))
l2.gulf = list("sp.lines", clp.gulf, lwd = 1, col = "black")
cr = brewer.pal(6, "Reds")
Plot Gulf of Mexico example:
spplot(gulfHSPDF, "tracks.count", col = "white", col.regions = cr, sp.layout = list(l1.gulf,
l2.gulf), at = seq(0, 6, 1), xlim = bbox(gulfHSPDF)[1, ], ylim = bbox(gulfHSPDF)[2,
], colorkey = list(space = "bottom", labels = paste(seq(0, 6, 1))), sub = expression("Number of Cyclones"),
main = "Observed TCs in the Gulf of Mexico in 2005")
Calculating per hexagon cyclone counts for the global grid of IBTrACS data:
ib.hexid = over(ib.sdf, hpg.i)
ib.sdf$hexid = ib.hexid
ib.ch = ddply(ib.sdf@data, .(hexid), "nrow")
ib.cn = ddply(ib.sdf@data, .(hexid), summarize, cn = length(unique(Sid)))
id.obs = ib.cn$hexid # for difference mapping later
ID = paste("ID", ib.ch$hexid, sep = "")
ib.chn = data.frame(nhours = ib.ch$nrow, ncyclones = ib.cn$cn)
rownames(ib.chn) = ID
hspdf = SpatialPolygonsDataFrame(hpg.i[ID], ib.chn, match.ID = TRUE)
Calculate per hexagon cyclone counts for the HiRAM, r1:
# load('C:/Users/Sarah/Dropbox/GCMTrackData/gfdlPmin_r1.RData')
# load('C:/Users/Sarah/Dropbox/GCMTrackData/track.df1.RData')
load("C:/Users/Sarah/Dropbox/GCMTrackData/track.df1B.RData")
gt.df1 = subset(track.df1, Yr > 1981 & Yr < 2009)
mean(gt.df1$Wmax)
## [1] 24.57
gt.df1 = subset(gt.df1, Wmax >= 17)
ch.g1 = nrow(gt.df1)
gt.sdf1 = gt.df1[, c("lon", "lat", "Sid", "Sn", "Yr", "Mo", "Da", "Wmax", "DWmaxDt",
"basins")]
coordinates(gt.sdf1) = c("lon", "lat")
proj4string(gt.sdf1) = CRS(ll)
gt.sdf1 = spTransform(gt.sdf1, CRS(moll))
l1.g1 = list("sp.points", gt.sdf1, pch = 20, col = "black", cex = 0.2)
gt.hexid1 = over(gt.sdf1, hpg.i)
gt.sdf1$hexid = gt.hexid1
gt.ch1 = ddply(gt.sdf1@data, .(hexid), "nrow")
gt.cn1 = ddply(gt.sdf1@data, .(hexid), summarize, cn = length(unique(Sid)))
id.mod1 = gt.cn1$hexid # for difference mapping later, should save gt.cn for each model instead
ID = paste("ID", gt.ch1$hexid, sep = "")
gt.chn1 = data.frame(nhours = gt.ch1$nrow, ncyclones = gt.cn1$cn)
rownames(gt.chn1) = ID
hspdfHiRAMr1 = SpatialPolygonsDataFrame(hpg.i[ID], gt.chn1, match.ID = TRUE)
rm(track.df1)
Calculate per hexagon storm counts for the HiRAM, r2
# load('C:/Users/Sarah//Dropbox/GCMTrackData/track.df2.RData')
load("C:/Users/Sarah//Dropbox/GCMTrackData/track.df2B.RData")
# load('C:/Users/Sarah/Dropbox/GCMTrackData/gfdlPmin_r2.RData')
gt.df2 = subset(track.df2, Yr > 1981 & Yr < 2009)
mean(gt.df2$Wmax)
## [1] 24.15
gt.df2 = subset(gt.df2, Wmax >= 17)
ch.g2 = nrow(gt.df2)
gt.sdf2 = gt.df2[, c("lon", "lat", "Sid", "Sn", "Yr", "Mo", "Da", "Wmax", "DWmaxDt",
"basins")]
coordinates(gt.sdf2) = c("lon", "lat")
proj4string(gt.sdf2) = CRS(ll)
gt.sdf2 = spTransform(gt.sdf2, CRS(moll))
l1.g2 = list("sp.points", gt.sdf2, pch = 20, col = "black", cex = 0.2)
gt.hexid2 = over(gt.sdf2, hpg.i)
gt.sdf2$hexid = gt.hexid2
gt.ch2 = ddply(gt.sdf2@data, .(hexid), "nrow")
gt.cn2 = ddply(gt.sdf2@data, .(hexid), summarize, cn = length(unique(Sid)))
id.mod2 = gt.cn2$hexid # for difference mapping later
ID = paste("ID", gt.ch2$hexid, sep = "")
gt.chn2 = data.frame(nhours = gt.ch2$nrow, ncyclones = gt.cn2$cn)
rownames(gt.chn2) = ID
hspdfHiRAMr2 = SpatialPolygonsDataFrame(hpg.i[ID], gt.chn2, match.ID = TRUE)
rm(track.df2)
Calculate per hexagon storm counts for the HiRAM, r3:
# load('C:/Users/Sarah/Dropbox/GCMTrackData/track.df3.RData')
load("C:/Users/Sarah/Dropbox/GCMTrackData/track.df3B.RData")
# load('C:/Users/Sarah/Dropbox/GCMTrackData/gfdlPmin_r3.RData')
gt.df3 = subset(track.df3, Yr > 1981 & Yr < 2009)
mean(gt.df3$Wmax)
## [1] 24.72
gt.df3 = subset(gt.df3, Wmax >= 17)
ch.g3 = nrow(gt.df3)
gt.sdf3 = gt.df3[, c("lon", "lat", "Sid", "Sn", "Yr", "Mo", "Da", "Wmax", "DWmaxDt",
"basins")]
coordinates(gt.sdf3) = c("lon", "lat")
proj4string(gt.sdf3) = CRS(ll)
gt.sdf3 = spTransform(gt.sdf3, CRS(moll))
l1.g3 = list("sp.points", gt.sdf3, pch = 20, col = "black", cex = 0.2)
gt.hexid3 = over(gt.sdf3, hpg.i)
gt.sdf3$hexid = gt.hexid3
gt.ch3 = ddply(gt.sdf3@data, .(hexid), "nrow")
gt.cn3 = ddply(gt.sdf3@data, .(hexid), summarize, cn = length(unique(Sid)))
id.mod3 = gt.cn3$hexid # for difference mapping later
ID = paste("ID", gt.ch3$hexid, sep = "")
gt.chn3 = data.frame(nhours = gt.ch3$nrow, ncyclones = gt.cn3$cn)
rownames(gt.chn3) = ID
hspdfHiRAMr3 = SpatialPolygonsDataFrame(hpg.i[ID], gt.chn3, match.ID = TRUE)
rm(track.df3)
Compute mean, min (17 ms-1), and max winds for obs and HiRAM:
mean(ib.df$Wmax)
min(ib.df$Wmax)
max(ib.df$Wmax)
mean(gt.df1$Wmax)
min(gt.df1$Wmax)
max(gt.df1$Wmax)
mean(gt.df2$Wmax)
min(gt.df2$Wmax)
max(gt.df2$Wmax)
mean(gt.df3$Wmax)
min(gt.df3$Wmax)
max(gt.df3$Wmax)
Set up global map:
Map cyclone hours for IBTrACS and GFDL HiRAM, r1:
p1 = spplot(hspdf, "ncyclones", col = "white", col.regions = cr2, sp.layout = list(l2.i),
at = seq(0, 300, 50), xlim = bbox(hspdf)[1, ], ylim = bbox(hspdf)[2, ],
colorkey = list(space = "bottom", labels = paste(seq(0, 300, 50))), sub = expression("Number of Cyclones"))
p2 = spplot(hspdfHiRAMr1, "ncyclones", col = "white", col.regions = cr2, sp.layout = list(l2.i),
at = seq(0, 300, 50), xlim = bbox(hspdf)[1, ], ylim = bbox(hspdf)[2, ],
colorkey = list(space = "bottom", labels = paste(seq(0, 300, 50))), sub = expression("Number of Cyclones"))
p1 = update(p1, main = textGrob("(a)", x = unit(0.05, "npc")), par.settings = list(fontsize = list(text = 14)))
p2 = update(p2, main = textGrob("(b)", x = unit(0.05, "npc")), par.settings = list(fontsize = list(text = 14)))
plot(p1, split = c(1, 1, 1, 2), more = TRUE)
plot(p2, split = c(1, 2, 1, 2), more = FALSE)
Calculate difference between IBTrACS and HiRAM cyclones counts:
z1 = merge(ib.cn, gt.cn1, by = "hexid", all = TRUE)
z1[is.na(z1)] = 0
z1$Diff = z1$cn.x - z1$cn.y # cn.x = obs, cn.y = mod
rownames(z1) = paste("ID", z1$hexid, sep = "")
hspdfDiffz1 = SpatialPolygonsDataFrame(hpg.i[rownames(z1)], z1, match.ID = TRUE)
z2 = merge(ib.cn, gt.cn2, by = "hexid", all = TRUE)
z2[is.na(z2)] = 0
z2$Diff = z2$cn.x - z2$cn.y
rownames(z2) = paste("ID", z2$hexid, sep = "")
hspdfDiffz2 = SpatialPolygonsDataFrame(hpg.i[rownames(z2)], z2, match.ID = TRUE)
z3 = merge(ib.cn, gt.cn3, by = "hexid", all = TRUE)
z3[is.na(z3)] = 0
z3$Diff = z3$cn.x - z3$cn.y
rownames(z3) = paste("ID", z3$hexid, sep = "")
hspdfDiffz3 = SpatialPolygonsDataFrame(hpg.i[rownames(z3)], z3, match.ID = TRUE)
Plot difference maps:
blue2white2red = colorRampPalette(c("blue", "white", "red"), space = "Lab")
p3 = spplot(hspdfDiffz1, "Diff", col = "white", col.regions = blue2white2red(10),
sp.layout = list(l2.i), at = seq(-125, 125, 25), xlim = bbox(hspdf)[1, ],
ylim = bbox(hspdf)[2, ], colorkey = list(space = "bottom", labels = paste(seq(-125,
125, 25))), sub = expression("Cyclone Count Difference"))
p4 = spplot(hspdfDiffz2, "Diff", col = "white", col.regions = blue2white2red(10),
sp.layout = list(l2.i), at = seq(-125, 125, 25), xlim = bbox(hspdf)[1, ],
ylim = bbox(hspdf)[2, ], colorkey = list(space = "bottom", labels = paste(seq(-125,
125, 25))), sub = expression("Cyclone Count Difference"))
p5 = spplot(hspdfDiffz3, "Diff", col = "white", col.regions = blue2white2red(10),
sp.layout = list(l2.i), at = seq(-125, 125, 25), xlim = bbox(hspdf)[1, ],
ylim = bbox(hspdf)[2, ], colorkey = list(space = "bottom", labels = paste(seq(-125,
125, 25))), sub = expression("Cyclone Count Difference"))
p3 = update(p3, main = textGrob("(a)", x = unit(0.05, "npc")), par.settings = list(fontsize = list(text = 14)))
p4 = update(p4, main = textGrob("(b)", x = unit(0.05, "npc")), par.settings = list(fontsize = list(text = 14)))
p5 = update(p5, main = textGrob("(c)", x = unit(0.05, "npc")), par.settings = list(fontsize = list(text = 14)))
plot(p3, split = c(1, 1, 1, 3), more = TRUE)
plot(p4, split = c(1, 2, 1, 3), more = TRUE)
plot(p5, split = c(1, 3, 1, 3), more = FALSE)
Bias = (Hits + False Positives) / (Hits + False Negatives). Overall bias (all years, all basins):
countsR1 = hspdfDiffz1@data
colnames(countsR1) = c("hexid", "OBS", "MOD", "DIFF")
countsR1$obsBin1 = countsR1$OBS > 200
countsR1$modBin1 = countsR1$MOD > 200
countsR1$hits = countsR1$obsBin1 == countsR1$modBin1
countsR1$FP = countsR1$obsBin1 < countsR1$modBin1
countsR1$FN = countsR1$obsBin1 > countsR1$modBin1
biasR1 = (sum(countsR1$hits) + sum(countsR1$FP))/(sum(countsR1$hits) + sum(countsR1$FN)) # 1.10 --> overprediction (in terms of area w/TCs)
CSI = (sum(countsR1$hits))/(sum(countsR1$hits) + sum(countsR1$FN) + sum(countsR1$FP))
# How much overprediction in general? length(unique(ib.df$Sid))
# length(unique(gt.df1$Sid)) length(unique(gt.df2$Sid))
# length(unique(gt.df3$Sid))
Calculate Probability of Detection (POD) 100 times for all basins (using grids the same size as those used in maps):
basin = c("NA", "EP", "WP", "SP", "NI", "SI", "SA")
###############################################################
pod1000B = matrix(0, nrow = 100, ncol = 7)
pod1000B = data.frame(pod1000B)
for (j in 1:length(basin)) {
ba = basin[j]
set.seed(2032)
for (i in 1:100) {
hpts = spsample(ib.sdf, type = "hexagonal", n = 1000, bb = bbox(ib.sdf) *
1.08)
hpg = HexPoints2SpatialPolygons(hpts)
pod1000B[i, j] = get.pod(ib.df, gt.df1, hpg, ba)
}
}
podPerform = pod1000B
save(podPerform, file = "podPerform.RData")
###############################################################
podAllB = rep(0, 100)
set.seed(2032)
for (i in 1:100) {
hpts = spsample(ib.sdf, type = "hexagonal", n = 1000, bb = bbox(ib.sdf) *
1.08)
hpg = HexPoints2SpatialPolygons(hpts)
podAllB[i] = grid.pod.tot(ib.df, gt.df1, hpg)
}
load("podPerform.RData")
podPerform$ALL = podAllB
save(podPerform, file = "podPerform.RData")
Calculate 95% C.I. for probability of detection:
load("podPerform.RData")
ciPerform.pod = matrix(0, nrow = 2, ncol = 8)
mPerform.pod = rep(0, 8)
for (i in 1:length(podPerform[1, ])) {
ciPerform.pod[1, i] = quantile(podPerform[, i], 0.025)
ciPerform.pod[2, i] = quantile(podPerform[, i], 0.975)
mPerform.pod[i] = median(podPerform[, i])
}
Calculate success rate 100 times for all basins (using same size hexagons as those used in maps):
basin = c("NA", "EP", "WP", "SP", "NI", "SI", "SA")
###############################################################
srPerform = matrix(0, nrow = 100, ncol = 7)
srPerform = data.frame(srPerform)
for (j in 1:length(basin)) {
ba = basin[j]
set.seed(2032)
for (i in 1:100) {
hpts = spsample(ib.sdf, type = "hexagonal", n = 1000, bb = bbox(ib.sdf) *
1.08)
hpg = HexPoints2SpatialPolygons(hpts)
srPerform[i, j] = get.sr(ib.df, gt.df1, hpg, ba)
}
}
save(srPerform, file = "srPerform.RData")
###############################################################
srAllB = rep(0, 100)
set.seed(2032)
for (i in 1:100) {
hpts = spsample(ib.sdf, type = "hexagonal", n = 1000, bb = bbox(ib.sdf) *
1.08)
hpg = HexPoints2SpatialPolygons(hpts)
srAllB[i] = grid.sr.tot(ib.df, gt.df1, hpg)
}
load("srPerform.RData")
srPerform$ALL = srAllB
save(srPerform, file = "srPerform.RData")
Calculate 95% C.I. for SR:
load("srPerform.RData")
ciPerform.sr = matrix(0, nrow = 2, ncol = 8)
mPerform.sr = rep(0, 8)
for (i in 1:length(srPerform[1, ])) {
ciPerform.sr[1, i] = quantile(srPerform[, i], 0.025)
ciPerform.sr[2, i] = quantile(srPerform[, i], 0.975)
mPerform.sr[i] = median(srPerform[, i])
}
Make performance diagram6
require(verification)
performance.diagram()
points(mPerform.sr[1], mPerform.pod[1], pch = 19, cex = 1)
text(mPerform.sr[1] + 0.025, mPerform.pod[1] - 0.02, labels = "NA")
segments(mPerform.sr[1], mPerform.pod[1], x1 = ciPerform.sr[1, 1])
segments(mPerform.sr[1], mPerform.pod[1], x1 = ciPerform.sr[2, 1])
segments(mPerform.sr[1], mPerform.pod[1], y1 = ciPerform.pod[1, 1])
segments(mPerform.sr[1], mPerform.pod[1], y1 = ciPerform.pod[2, 1])
points(mPerform.sr[2], mPerform.pod[2], pch = 19, cex = 1)
text(mPerform.sr[2] - 0.025, mPerform.pod[2] + 0.02, labels = "EP")
segments(mPerform.sr[2], mPerform.pod[2], x1 = ciPerform.sr[1, 2])
segments(mPerform.sr[2], mPerform.pod[2], x1 = ciPerform.sr[2, 2])
segments(mPerform.sr[2], mPerform.pod[2], y1 = ciPerform.pod[1, 2])
segments(mPerform.sr[2], mPerform.pod[2], y1 = ciPerform.pod[2, 2])
points(mPerform.sr[3], mPerform.pod[3], pch = 19, cex = 1)
text(mPerform.sr[3] + 0.025, mPerform.pod[3] + 0.02, labels = "WP")
segments(mPerform.sr[3], mPerform.pod[3], x1 = ciPerform.sr[1, 3])
segments(mPerform.sr[3], mPerform.pod[3], x1 = ciPerform.sr[2, 3])
segments(mPerform.sr[3], mPerform.pod[3], y1 = ciPerform.pod[1, 3])
segments(mPerform.sr[3], mPerform.pod[3], y1 = ciPerform.pod[2, 3])
points(mPerform.sr[4], mPerform.pod[4], pch = 19, cex = 1)
text(mPerform.sr[4] + 0.025, mPerform.pod[4] - 0.02, labels = "SP")
segments(mPerform.sr[4], mPerform.pod[4], x1 = ciPerform.sr[1, 4])
segments(mPerform.sr[4], mPerform.pod[4], x1 = ciPerform.sr[2, 4])
segments(mPerform.sr[4], mPerform.pod[4], y1 = ciPerform.pod[1, 4])
segments(mPerform.sr[4], mPerform.pod[4], y1 = ciPerform.pod[2, 4])
points(mPerform.sr[5], mPerform.pod[5], pch = 19, cex = 1)
text(mPerform.sr[5] + 0.025, mPerform.pod[5] + 0.02, labels = "NI")
segments(mPerform.sr[5], mPerform.pod[5], x1 = ciPerform.sr[1, 5])
segments(mPerform.sr[5], mPerform.pod[5], x1 = ciPerform.sr[2, 5])
segments(mPerform.sr[5], mPerform.pod[5], y1 = ciPerform.pod[1, 5])
segments(mPerform.sr[5], mPerform.pod[5], y1 = ciPerform.pod[2, 5])
points(mPerform.sr[6], mPerform.pod[6], pch = 19, cex = 1)
text(mPerform.sr[6] - 0.025, mPerform.pod[6] - 0.02, labels = "SI")
segments(mPerform.sr[6], mPerform.pod[6], x1 = ciPerform.sr[1, 6])
segments(mPerform.sr[6], mPerform.pod[6], x1 = ciPerform.sr[2, 6])
segments(mPerform.sr[6], mPerform.pod[6], y1 = ciPerform.pod[1, 6])
segments(mPerform.sr[6], mPerform.pod[6], y1 = ciPerform.pod[2, 6])
points(mPerform.sr[7], mPerform.pod[7], pch = 19, cex = 1)
text(mPerform.sr[7] + 0.025, mPerform.pod[7] - 0.02, labels = "SA")
segments(mPerform.sr[7], mPerform.pod[7], x1 = ciPerform.sr[1, 7])
segments(mPerform.sr[7], mPerform.pod[7], x1 = ciPerform.sr[2, 7])
segments(mPerform.sr[7], mPerform.pod[7], y1 = ciPerform.pod[1, 7])
segments(mPerform.sr[7], mPerform.pod[7], y1 = ciPerform.pod[2, 7])
points(mPerform.sr[8], mPerform.pod[8], pch = 5, cex = 1.5)
text(mPerform.sr[8] + 0.025, mPerform.pod[8] - 0.02, labels = "ALL")
segments(mPerform.sr[8], mPerform.pod[8], x1 = ciPerform.sr[1, 8], lwd = 4)
segments(mPerform.sr[8], mPerform.pod[8], x1 = ciPerform.sr[2, 8], lwd = 4)
segments(mPerform.sr[8], mPerform.pod[8], y1 = ciPerform.pod[1, 8], lwd = 4)
segments(mPerform.sr[8], mPerform.pod[8], y1 = ciPerform.pod[2, 8], lwd = 4)
Estimate uncertainty for all basins, n=500 hexagons:
basin = c("NA", "EP", "WP", "SP", "NI", "SI", "SA")
b500B = matrix(0, nrow = 100, ncol = 7)
b500B = data.frame(b500B)
for (j in 1:length(basin)) {
ba = basin[j]
set.seed(2032)
for (i in 1:100) {
hpts = spsample(ib.sdf, type = "hexagonal", n = 500, bb = bbox(ib.sdf) *
1.08)
hpg = HexPoints2SpatialPolygons(hpts)
b500B[i, j] = grid.bias(ib.df, gt.df1, ba, hpg)
}
}
save(b500B, file = "b500B.RData")
Estimate uncertainty for all basins, n = 750 hexagons
b750B = matrix(0, nrow = 100, ncol = 7)
b750B = data.frame(b750B)
for (j in 1:length(basin)) {
ba = basin[j]
set.seed(2032)
for (i in 1:100) {
hpts = spsample(ib.sdf, type = "hexagonal", n = 750, bb = bbox(ib.sdf) *
1.08)
hpg = HexPoints2SpatialPolygons(hpts)
b750B[i, j] = grid.bias(ib.df, gt.df1, ba, hpg)
}
}
save(b750B, file = "b750B.RData")
Estimate uncertainty for all basins, n = 1000 hexagons
b1000B = matrix(0, nrow = 100, ncol = 7)
b1000B = data.frame(b1000B)
for (j in 1:length(basin)) {
ba = basin[j]
set.seed(2032)
for (i in 1:100) {
hpts = spsample(ib.sdf, type = "hexagonal", n = 1000, bb = bbox(ib.sdf) *
1.08)
hpg = HexPoints2SpatialPolygons(hpts)
b1000B[i, j] = grid.bias(ib.df, gt.df1, ba, hpg)
}
}
save(b1000B, file = "b1000B.RData")
Estimate uncertainty for all basins, n = 1500 hexagons
b1500B = matrix(0, nrow = 100, ncol = 7)
b1500B = data.frame(b1500B)
for (j in 1:length(basin)) {
ba = basin[j]
set.seed(2032)
for (i in 1:100) {
hpts = spsample(ib.sdf, type = "hexagonal", n = 1500, bb = bbox(ib.sdf) *
1.08)
hpg = HexPoints2SpatialPolygons(hpts)
b1500B[i, j] = grid.bias(ib.df, gt.df1, ba, hpg)
}
}
save(b1500B, file = "b1500B.RData")
Estimate uncertainty for all basins, n = 2000 hexagons
b2000B = matrix(0, nrow = 100, ncol = 7)
b2000B = data.frame(b2000B)
for (j in 1:length(basin)) {
ba = basin[j]
set.seed(2032)
for (i in 1:100) {
hpts = spsample(ib.sdf, type = "hexagonal", n = 2000, bb = bbox(ib.sdf) *
1.08)
hpg = HexPoints2SpatialPolygons(hpts)
b2000B[i, j] = grid.bias(ib.df, gt.df1, ba, hpg)
}
}
save(b2000B, file = "b2000B.RData")
Get uncertainty, global grids
bAllB = matrix(0, nrow = 100, ncol = 5)
bAllB = data.frame(bAllB)
grids = c(500, 750, 1000, 1500, 2000)
for (j in 1:length(grids)) {
num = grids[j]
set.seed(2032)
for (i in 1:100) {
hpts = spsample(ib.sdf, type = "hexagonal", n = num, bb = bbox(ib.sdf) *
1.08)
hpg = HexPoints2SpatialPolygons(hpts)
bAllB[i, j] = grid.bias.tot(ib.df, gt.df1, hpg)
}
}
save(bAllB, file = "bAllB.RData")
Calculate 90% C.I.:
load("b500B.RData")
load("b750B.RData")
load("b1000B.RData")
load("b1500B.RData")
load("b2000B.RData")
load("bAllB.RData")
ci500 = matrix(0, nrow = 3, ncol = 7)
for (i in 1:length(b500B[1, ])) {
ci500[1, i] = quantile(b500B[, i], 0.5)
ci500[2, i] = quantile(b500B[, i], 0.025)
ci500[3, i] = quantile(b500B[, i], 0.975)
}
ci750 = matrix(0, nrow = 3, ncol = 7)
for (i in 1:length(b750B[1, ])) {
ci750[1, i] = quantile(b750B[, i], 0.5)
ci750[2, i] = quantile(b750B[, i], 0.025)
ci750[3, i] = quantile(b750B[, i], 0.975)
}
ci1000 = matrix(0, nrow = 3, ncol = 7)
for (i in 1:length(b1000B[1, ])) {
ci1000[1, i] = quantile(b1000B[, i], 0.5)
ci1000[2, i] = quantile(b1000B[, i], 0.025)
ci1000[3, i] = quantile(b1000B[, i], 0.975)
}
ci1500 = matrix(0, nrow = 3, ncol = 7)
for (i in 1:length(b1500B[1, ])) {
ci1500[1, i] = quantile(b1500B[, i], 0.5)
ci1500[2, i] = quantile(b1500B[, i], 0.025)
ci1500[3, i] = quantile(b1500B[, i], 0.975)
}
ci2000 = matrix(0, nrow = 3, ncol = 7)
for (i in 1:length(b2000B[1, ])) {
ci2000[1, i] = quantile(b2000B[, i], 0.5)
ci2000[2, i] = quantile(b2000B[, i], 0.025)
ci2000[3, i] = quantile(b2000B[, i], 0.975)
}
ciALL = matrix(0, nrow = 5, ncol = 3)
for (i in 1:length(bAllB[1, ])) {
ciALL[i, 1] = quantile(bAllB[, i], 0.5)
ciALL[i, 2] = quantile(bAllB[, i], 0.025)
ciALL[i, 3] = quantile(bAllB[, i], 0.975)
}
Estimate (CSI) uncertainty for all basins, n=500 hexagons:
basin = c("NA", "EP", "WP", "SP", "NI", "SI", "SA")
c500B = matrix(0, nrow = 100, ncol = 7)
c500B = data.frame(c500B)
for (j in 1:length(basin)) {
ba = basin[j]
set.seed(2032)
for (i in 1:100) {
hpts = spsample(ib.sdf, type = "hexagonal", n = 500, bb = bbox(ib.sdf) *
1.08)
hpg = HexPoints2SpatialPolygons(hpts)
c500B[i, j] = get.csi(ib.df, gt.df1, ba, hpg)
}
}
save(c500B, file = "c500B.RData")
Estimate uncertainty for all basins, n = 750 hexagons
c750B = matrix(0, nrow = 100, ncol = 7)
c750B = data.frame(c750B)
for (j in 1:length(basin)) {
ba = basin[j]
set.seed(2032)
for (i in 1:100) {
hpts = spsample(ib.sdf, type = "hexagonal", n = 750, bb = bbox(ib.sdf) *
1.08)
hpg = HexPoints2SpatialPolygons(hpts)
c750B[i, j] = get.csi(ib.df, gt.df1, ba, hpg)
}
}
save(c750B, file = "c750B.RData")
Estimate uncertainty for all basins, n = 1000 hexagons
c1000B = matrix(0, nrow = 100, ncol = 7)
c1000B = data.frame(c1000B)
for (j in 1:length(basin)) {
ba = basin[j]
set.seed(2032)
for (i in 1:100) {
hpts = spsample(ib.sdf, type = "hexagonal", n = 1000, bb = bbox(ib.sdf) *
1.08)
hpg = HexPoints2SpatialPolygons(hpts)
c1000B[i, j] = get.csi(ib.df, gt.df1, ba, hpg)
}
}
save(c1000B, file = "c1000B.RData")
Estimate uncertainty for all basins, n = 1500 hexagons
c1500B = matrix(0, nrow = 100, ncol = 7)
c1500B = data.frame(c1500B)
for (j in 1:length(basin)) {
ba = basin[j]
set.seed(2032)
for (i in 1:100) {
hpts = spsample(ib.sdf, type = "hexagonal", n = 1500, bb = bbox(ib.sdf) *
1.08)
hpg = HexPoints2SpatialPolygons(hpts)
c1500B[i, j] = get.csi(ib.df, gt.df1, ba, hpg)
}
}
save(c1500B, file = "c1500B.RData")
Estimate uncertainty for all basins, n = 2000 hexagons
c2000B = matrix(0, nrow = 100, ncol = 7)
c2000B = data.frame(c2000B)
for (j in 1:length(basin)) {
ba = basin[j]
set.seed(2032)
for (i in 1:100) {
hpts = spsample(ib.sdf, type = "hexagonal", n = 2000, bb = bbox(ib.sdf) *
1.08)
hpg = HexPoints2SpatialPolygons(hpts)
c2000B[i, j] = get.csi(ib.df, gt.df1, ba, hpg)
}
}
save(c2000B, file = "c2000B.RData")
Get uncertainty, global grids
cAllB = matrix(0, nrow = 100, ncol = 5)
cAllB = data.frame(cAllB)
grids = c(500, 750, 1000, 1500, 2000)
for (j in 1:length(grids)) {
num = grids[j]
set.seed(2032)
for (i in 1:100) {
hpts = spsample(ib.sdf, type = "hexagonal", n = num, bb = bbox(ib.sdf) *
1.08)
hpg = HexPoints2SpatialPolygons(hpts)
cAllB[i, j] = grid.csi.tot(ib.df, gt.df1, hpg)
}
}
save(cAllB, file = "cAllB.RData")
Calculate 90% C.I. for CSI:
load("c500B.RData")
load("c750B.RData")
load("c1000B.RData")
load("c1500B.RData")
load("c2000B.RData")
load("cAllB.RData")
ci500CSI = matrix(0, nrow = 3, ncol = 7)
for (i in 1:length(c500B[1, ])) {
ci500CSI[1, i] = quantile(c500B[, i], 0.5)
ci500CSI[2, i] = quantile(c500B[, i], 0.025)
ci500CSI[3, i] = quantile(c500B[, i], 0.975)
}
ci750CSI = matrix(0, nrow = 3, ncol = 7)
for (i in 1:length(c750B[1, ])) {
ci750CSI[1, i] = quantile(c750B[, i], 0.5)
ci750CSI[2, i] = quantile(c750B[, i], 0.025)
ci750CSI[3, i] = quantile(c750B[, i], 0.975)
}
ci1000CSI = matrix(0, nrow = 3, ncol = 7)
for (i in 1:length(c1000B[1, ])) {
ci1000CSI[1, i] = quantile(c1000B[, i], 0.5)
ci1000CSI[2, i] = quantile(c1000B[, i], 0.025)
ci1000CSI[3, i] = quantile(c1000B[, i], 0.975)
}
ci1500CSI = matrix(0, nrow = 3, ncol = 7)
for (i in 1:length(c1500B[1, ])) {
ci1500CSI[1, i] = quantile(c1500B[, i], 0.5)
ci1500CSI[2, i] = quantile(c1500B[, i], 0.025)
ci1500CSI[3, i] = quantile(c1500B[, i], 0.975)
}
ci2000CSI = matrix(0, nrow = 3, ncol = 7)
for (i in 1:length(c2000B[1, ])) {
ci2000CSI[1, i] = quantile(c2000B[, i], 0.5)
ci2000CSI[2, i] = quantile(c2000B[, i], 0.025)
ci2000CSI[3, i] = quantile(c2000B[, i], 0.975)
}
ciALLCSI = matrix(0, nrow = 5, ncol = 3)
for (i in 1:length(cAllB[1, ])) {
ciALLCSI[i, 1] = quantile(cAllB[, i], 0.5)
ciALLCSI[i, 2] = quantile(cAllB[, i], 0.025)
ciALLCSI[i, 3] = quantile(cAllB[, i], 0.975)
}
Look at convergence to 1.0 as set higher thresholds…
We next compare the GFDL HiRAM simulated tracks with those generated from the FSU-COAPS spectral model. Because global tracks are not yet available for the FSU-COAPS model, we focus our comparison on the North Atlantic Basin.
Load observations, assign a lambert conformal conic projection, and calculate per hexagon TC counts:
# con = url('http://myweb.fsu.edu/jelsner/world.interp.v02r01.R')
# load(con) load('C:/Documents and Settings/sstrazzo/My
# Documents/Research/world.interp.v02r01.R')
load("G:Research/world.interp.v02r01.R")
ib.atl = subset(world.interp, Yr > 1981 & Yr < 2009)
rm(world.interp)
ib.atl = subset(ib.atl, Basin == "NA")
ib.atl$Wmax = ib.atl$Wmax * 0.5144
ib.gen = ib.atl
ib.atl = subset(ib.atl, Wmax >= 17)
ib.sATL = ib.atl[, c("name", "lon", "lat", "Sid", "Sn", "Yr", "Mo", "Da", "Wmax",
"DWmaxDt")]
ib.sGen = ib.gen[, c("name", "lon", "lat", "Sid", "Sn", "Yr", "Mo", "Da", "Wmax",
"DWmaxDt")]
coordinates(ib.sATL) = c("lon", "lat")
coordinates(ib.sGen) = c("lon", "lat")
ll = "+proj=longlat +datum=WGS84"
proj4string(ib.sATL) = CRS(ll)
proj4string(ib.sGen) = CRS(ll)
lcc = "+proj=lcc +lat_1=60 +lat_2=30 +lon_0=-60 +units=km"
ib.sATL = spTransform(ib.sATL, CRS(lcc))
ib.sGen = spTransform(ib.sGen, CRS(lcc))
hpt.atl = spsample(ib.sATL, type = "hexagonal", n = 500, bb = bbox(ib.sATL) *
1.08, offset = c(0, 0), iter = 100)
hpt.gen = spsample(ib.sGen, type = "hexagonal", n = 500, bb = bbox(ib.sATL) *
1.08, offset = c(0, 0), iter = 100)
hpg.atl = HexPoints2SpatialPolygons(hpt.atl)
hpg.gen = HexPoints2SpatialPolygons(hpt.gen)
ib.idATL = over(ib.sATL, hpg.atl)
ib.idGen = over(ib.sGen, hpg.gen)
ib.sATL$hexid = ib.idATL
ib.sGen$hexid = ib.idGen
ib.ch = ddply(ib.sATL@data, .(hexid), "nrow")
ib.ch = ib.ch[!is.na(ib.ch$hexid), ]
ib.cn = ddply(ib.sATL@data, .(hexid), summarize, cn = length(unique(Sid)))
ib.cn = ib.cn[!is.na(ib.cn$hexid), ]
id.obsATL = ib.cn$hexid # for difference mapping later
ID = paste("ID", ib.ch$hexid, sep = "")
ib.chn = data.frame(nhours = ib.ch$nrow, ncyclones = ib.cn$cn)
rownames(ib.chn) = ID
hspdf.atl = SpatialPolygonsDataFrame(hpg.atl[ID], ib.chn, match.ID = TRUE)
Load FSU, r1, assign projection, and compute per hexagon cyclone counts:
## repeat above for FSU tracks (r1, r2, & r3)
load("C:/Users/Sarah/Dropbox/GCMTrackData/fsu_r1i1p1_ATL.RData")
# load('C:/Users/Sarah/Dropbox/GCMTrackData/fsuPmin_r1.Rdata')
fsu.r1 = fsuTracks.df
rm(fsuTracks.df)
fsu.r1 = subset(fsu.r1, Yr > 1981 & Yr < 2009)
fsu.gen = fsu.r1
fsu.r1 = subset(fsu.r1, Wmax >= 17)
fsu.sdf1 = fsu.r1[, c("name", "lon", "lat", "Sid", "Sn", "Yr", "Mo", "Da", "Wmax",
"DWmaxDt")]
fsu.sGen = fsu.gen[, c("name", "lon", "lat", "Sid", "Sn", "Yr", "Mo", "Da",
"Wmax", "DWmaxDt")]
coordinates(fsu.sdf1) = c("lon", "lat")
coordinates(fsu.sGen) = c("lon", "lat")
proj4string(fsu.sdf1) = CRS(ll)
proj4string(fsu.sGen) = CRS(ll)
fsu.sdf1 = spTransform(fsu.sdf1, CRS(lcc))
fsu.sGen = spTransform(fsu.sGen, CRS(lcc))
fsu.hexid1 = over(fsu.sdf1, hpg.atl)
fsu.idGen = over(fsu.sGen, hpg.gen)
fsu.sdf1$hexid = fsu.hexid1
fsu.sGen$hexid = fsu.idGen
fsu.ch = ddply(fsu.sdf1@data, .(hexid), "nrow")
fsu.ch = fsu.ch[!is.na(fsu.ch$hexid), ]
fsu.cn = ddply(fsu.sdf1@data, .(hexid), summarize, cn = length(unique(Sid)))
fsu.cn = fsu.cn[!is.na(fsu.cn$hexid), ]
id.fsu1 = fsu.cn$hexid # for difference mapping later
ID = paste("ID", fsu.ch$hexid, sep = "")
fsu.chn = data.frame(nhours = fsu.ch$nrow, ncyclones = fsu.cn$cn)
rownames(fsu.chn) = ID
hspdf.fsu1 = SpatialPolygonsDataFrame(hpg.atl[ID], fsu.chn, match.ID = TRUE)
Load FSU, r2, assign projection, and compute per hexagon cyclone counts:
## repeat above for FSU tracks (r1, r2, & r3)
load("C:/Users/Sarah/Dropbox/GCMTrackData/fsu_r2i1p1_ATL.RData")
# load('C:/Users/Sarah/Dropbox/GCMTrackData/fsuPmin_r2.RData')
fsu.r2 = fsuTracks.df
rm(fsuTracks.df)
fsu.r2 = subset(fsu.r2, Yr > 1981 & Yr < 2009)
fsu.r2 = subset(fsu.r2, Wmax >= 17)
fsu.sdf2 = fsu.r2[, c("name", "lon", "lat", "Sid", "Sn", "Yr", "Mo", "Da", "Wmax",
"DWmaxDt")]
coordinates(fsu.sdf2) = c("lon", "lat")
proj4string(fsu.sdf2) = CRS(ll)
fsu.sdf2 = spTransform(fsu.sdf1, CRS(lcc))
fsu.hexid2 = over(fsu.sdf2, hpg.atl)
fsu.sdf2$hexid = fsu.hexid2
fsu.ch = ddply(fsu.sdf2@data, .(hexid), "nrow")
fsu.ch = fsu.ch[!is.na(fsu.ch$hexid), ]
fsu.cn = ddply(fsu.sdf2@data, .(hexid), summarize, cn = length(unique(Sid)))
fsu.cn = fsu.cn[!is.na(fsu.cn$hexid), ]
id.fsu2 = fsu.cn$hexid # for difference mapping later
ID = paste("ID", fsu.ch$hexid, sep = "")
fsu.chn = data.frame(nhours = fsu.ch$nrow, ncyclones = fsu.cn$cn)
rownames(fsu.chn) = ID
hspdf.fsu2 = SpatialPolygonsDataFrame(hpg.atl[ID], fsu.chn, match.ID = TRUE)
Load FSU, r3, assign projection, and compute per hexagon cyclone counts:
## repeat above for FSU tracks (r1, r2, & r3)
load("C:/Users/Sarah/Dropbox/GCMTrackData/fsu_r3i1p1_ATL.RData")
# load('C:/Users/Sarah/Dropbox/GCMTrackData/fsuPmin_r3.RData')
fsu.r3 = fsuTracks.df
rm(fsuTracks.df)
fsu.r3 = subset(fsu.r3, Yr > 1981 & Yr < 2009)
fsu.r3 = subset(fsu.r3, Wmax >= 17)
fsu.sdf3 = fsu.r3[, c("name", "lon", "lat", "Sid", "Sn", "Yr", "Mo", "Da", "Wmax",
"DWmaxDt")]
coordinates(fsu.sdf3) = c("lon", "lat")
proj4string(fsu.sdf3) = CRS(ll)
fsu.sdf3 = spTransform(fsu.sdf1, CRS(lcc))
fsu.hexid3 = over(fsu.sdf3, hpg.atl)
fsu.sdf3$hexid = fsu.hexid3
fsu.ch = ddply(fsu.sdf3@data, .(hexid), "nrow")
fsu.ch = fsu.ch[!is.na(fsu.ch$hexid), ]
fsu.cn = ddply(fsu.sdf3@data, .(hexid), summarize, cn = length(unique(Sid)))
fsu.cn = fsu.cn[!is.na(fsu.cn$hexid), ]
id.fsu3 = fsu.cn$hexid # for difference mapping later
ID = paste("ID", fsu.ch$hexid, sep = "")
fsu.chn = data.frame(nhours = fsu.ch$nrow, ncyclones = fsu.cn$cn)
rownames(fsu.chn) = ID
hspdf.fsu3 = SpatialPolygonsDataFrame(hpg.atl[ID], fsu.chn, match.ID = TRUE)
Load GFDL-HiRAM, r1, assign projection, and compute per hexagon cyclone counts:
## repeat above for GFDL tracks (r1, r2, & r3)
load("C:/Users/Sarah/Dropbox/GCMTrackData/gfdlTracks_r1.df.RData")
gfdl.na1 = natltracks.df1
rm(natltracks.df1)
gfdl.na1 = subset(gfdl.na1, Yr > 1981 & Yr < 2009)
gfdl.gen = gfdl.na1
gfdl.na1 = subset(gfdl.na1, Wmax >= 17)
gfdl.sna1 = gfdl.na1[, c("name", "lon", "lat", "Sid", "Sn", "Yr", "Mo", "Da",
"Wmax", "DWmaxDt")]
gfdl.sGen = gfdl.gen[, c("name", "lon", "lat", "Sid", "Sn", "Yr", "Mo", "Da",
"Wmax", "DWmaxDt")]
coordinates(gfdl.sna1) = c("lon", "lat")
coordinates(gfdl.sGen) = c("lon", "lat")
proj4string(gfdl.sna1) = CRS(ll)
proj4string(gfdl.sGen) = CRS(ll)
gfdl.sna1 = spTransform(gfdl.sna1, CRS(lcc))
gfdl.sGen = spTransform(gfdl.sGen, CRS(lcc))
gfdl.hexid1 = over(gfdl.sna1, hpg.atl)
gfdl.idGen = over(gfdl.sGen, hpg.gen)
gfdl.sna1$hexid = gfdl.hexid1
gfdl.sGen$hexid = gfdl.idGen
gfdl.ch = ddply(gfdl.sna1@data, .(hexid), "nrow")
gfdl.ch = gfdl.ch[!is.na(gfdl.ch$hexid), ]
gfdl.cn = ddply(gfdl.sna1@data, .(hexid), summarize, cn = length(unique(Sid)))
gfdl.cn = gfdl.cn[!is.na(gfdl.cn$hexid), ]
id.gfdl1 = gfdl.cn$hexid # for difference mapping later
ID = paste("ID", gfdl.ch$hexid, sep = "")
gfdl.chn = data.frame(nhours = gfdl.ch$nrow, ncyclones = gfdl.cn$cn)
rownames(gfdl.chn) = ID
hspdf.na.gfdl1 = SpatialPolygonsDataFrame(hpg.atl[ID], gfdl.chn, match.ID = TRUE)
Load GFDL-HiRAM, r2, assign projection, and compute per hexagon cyclone counts:
## repeat above for GFDL tracks (r1, r2, & r3)
load("C:/Users/Sarah/Dropbox/GCMTrackData/gfdlTracks_r2.df.RData")
gfdl.na2 = natltracks.df1
rm(natltracks.df1)
gfdl.na2 = subset(gfdl.na2, Yr > 1981 & Yr < 2009)
gfdl.na2 = subset(gfdl.na2, Wmax >= 17)
gfdl.sna2 = gfdl.na2[, c("name", "lon", "lat", "Sid", "Sn", "Yr", "Mo", "Da",
"Wmax", "DWmaxDt")]
coordinates(gfdl.sna2) = c("lon", "lat")
proj4string(gfdl.sna2) = CRS(ll)
# gfdl.sna2 = spTransform(gfdl.sna2, CRS(moll.atl))
gfdl.sna2 = spTransform(gfdl.sna2, CRS(lcc))
gfdl.hexid2 = over(gfdl.sna2, hpg.atl)
gfdl.sna2$hexid = gfdl.hexid2
gfdl.ch = ddply(gfdl.sna2@data, .(hexid), "nrow")
gfdl.ch = gfdl.ch[!is.na(gfdl.ch$hexid), ]
gfdl.cn = ddply(gfdl.sna2@data, .(hexid), summarize, cn = length(unique(Sid)))
gfdl.cn = gfdl.cn[!is.na(gfdl.cn$hexid), ]
id.gfdl2 = gfdl.cn$hexid # for difference mapping later
ID = paste("ID", gfdl.ch$hexid, sep = "")
gfdl.chn = data.frame(nhours = gfdl.ch$nrow, ncyclones = gfdl.cn$cn)
rownames(gfdl.chn) = ID
hspdf.na.gfdl2 = SpatialPolygonsDataFrame(hpg.atl[ID], gfdl.chn, match.ID = TRUE)
Load GFDL-HiRAM, r3, assign projection, and compute per hexagon cyclone counts:
## repeat above for GFDL tracks (r1, r2, & r3)
load("C:/Users/Sarah/Dropbox/GCMTrackData/gfdlTracks_r3.df.RData")
gfdl.na3 = natltracks.df1
rm(natltracks.df1)
gfdl.na3 = subset(gfdl.na3, Yr > 1981 & Yr < 2009)
gfdl.na3 = subset(gfdl.na3, Wmax >= 17)
gfdl.sna3 = gfdl.na3[, c("name", "lon", "lat", "Sid", "Sn", "Yr", "Mo", "Da",
"Wmax", "DWmaxDt")]
coordinates(gfdl.sna3) = c("lon", "lat")
proj4string(gfdl.sna3) = CRS(ll)
# gfdl.sna3 = spTransform(gfdl.sna3, CRS(moll.atl))
gfdl.sna3 = spTransform(gfdl.sna3, CRS(lcc))
gfdl.hexid3 = over(gfdl.sna3, hpg.atl)
gfdl.sna3$hexid = gfdl.hexid3
gfdl.ch = ddply(gfdl.sna3@data, .(hexid), "nrow")
gfdl.ch = gfdl.ch[!is.na(gfdl.ch$hexid), ]
gfdl.cn = ddply(gfdl.sna3@data, .(hexid), summarize, cn = length(unique(Sid)))
gfdl.cn = gfdl.cn[!is.na(gfdl.cn$hexid), ]
id.gfdl3 = gfdl.cn$hexid # for difference mapping later
ID = paste("ID", gfdl.ch$hexid, sep = "")
gfdl.chn = data.frame(nhours = gfdl.ch$nrow, ncyclones = gfdl.cn$cn)
rownames(gfdl.chn) = ID
hspdf.na.gfdl3 = SpatialPolygonsDataFrame(hpg.atl[ID], gfdl.chn, match.ID = TRUE)
Plot TC counts for the North Atlantic (obs, FSU-COAPS, GFDL-HiRAM):
cl.atl = map("world", xlim = c(-120, 10), ylim = c(0, 65), plot = FALSE, wrap = TRUE) #xlim,ylim reduces comp
clp.atl = map2SpatialLines(cl.atl, proj4string = CRS(ll))
# clp.atl = spTransform(clp.atl, CRS(moll.atl))
clp.atl = spTransform(clp.atl, CRS(lcc))
l2.atl = list("sp.lines", clp.atl, lwd = 0.2)
# l1.atl = list('sp.points', ib.sATL, pch=20, col='black', cex=.1)
cr = brewer.pal(6, "Reds")
p1 = spplot(hspdf.atl, "ncyclones", col = "white", col.regions = cr, sp.layout = list(l2.atl),
at = seq(0, 150, 25), xlim = bbox(hspdf.fsu1)[1, ], ylim = bbox(hspdf.fsu1)[2,
], colorkey = list(space = "bottom", labels = paste(seq(0, 150, 25))),
sub = expression("Number of Cyclones"))
p2 = spplot(hspdf.fsu1, "ncyclones", col = "white", col.regions = cr, sp.layout = list(l2.atl),
at = seq(0, 150, 25), xlim = bbox(hspdf.fsu1)[1, ], ylim = bbox(hspdf.fsu1)[2,
], colorkey = list(space = "bottom", labels = paste(seq(0, 150, 25))),
sub = expression("Number of Cyclones"))
p3 = spplot(hspdf.na.gfdl1, "ncyclones", col = "white", col.regions = cr, sp.layout = list(l2.atl),
at = seq(0, 150, 25), xlim = bbox(hspdf.fsu1)[1, ], ylim = bbox(hspdf.fsu1)[2,
], colorkey = list(space = "bottom", labels = paste(seq(0, 150, 25))),
sub = expression("Number of Cyclones"))
p1 = update(p1, main = textGrob("(a)", x = unit(0.05, "npc")), par.settings = list(fontsize = list(text = 14)))
p2 = update(p2, main = textGrob("(b)", x = unit(0.05, "npc")), par.settings = list(fontsize = list(text = 14)))
p3 = update(p3, main = textGrob("(c)", x = unit(0.05, "npc")), par.settings = list(fontsize = list(text = 14)))
plot(p1, split = c(1, 1, 1, 3), more = TRUE)
plot(p2, split = c(1, 2, 1, 3), more = TRUE)
plot(p3, split = c(1, 3, 1, 3), more = FALSE)
Calculate the factor by which modeled TC counts exceed observed counts (“relative risk”):
Obs = hspdf.atl@data
Mod = hspdf.fsu1@data
Obs$hexid = rownames(Obs)
Mod$hexid = rownames(Mod)
z = merge(Obs, Mod, by = "hexid", all = TRUE)
z[is.na(z)] = 0
rownames(z) = z$hexid
# Frequency of TCs relative to climatology
z$fObs = z$ncyclones.x/length(unique(ib.atl$Sid))
z$fMod = z$ncyclones.y/length(unique(fsu.r1$Sid))
cor(z$fMod, z$fObs)
## [1] 0.6228
# Factor by which the modeled TC frequency exceeds the observed TC
# frequency relative to the observed climatology
z$rr = z$fMod/z$fObs
z$l2rr = log(z$rr, base = 2)
labs = c(-32, -16, -8, -4, -2, 0, 2, 4, 8, 16, 32)
hspdfRR = SpatialPolygonsDataFrame(hpg.atl[rownames(z)], z, match.ID = TRUE)
Mod2 = hspdf.na.gfdl1@data
Mod2$hexid = rownames(Mod2)
z2 = merge(Obs, Mod2, by = "hexid", all = TRUE)
z2[is.na(z2)] = 0
rownames(z2) = z2$hexid
# Frequency of TCs relative to climatology
z2$fObs = z2$ncyclones.x/length(unique(ib.atl$Sid))
z2$fMod2 = z2$ncyclones.y/length(unique(gfdl.na1$Sid))
cor(z2$fMod2, z2$fObs)
## [1] 0.7956
# Factor by which the modeled TC frequency exceeds the observed TC
# frequency relative to the observed climatology
z2$rr = z2$fMod/z2$fObs
z2$l2rr = log(z2$rr, base = 2)
labs = c(-32, -16, -8, -4, -2, 0, 2, 4, 8, 16, 32)
hspdfRR2 = SpatialPolygonsDataFrame(hpg.atl[rownames(z2)], z2, match.ID = TRUE)
Plot relative risk:
blue2white2red = colorRampPalette(c("blue", "white", "red"), space = "Lab")
p1 = spplot(hspdfRR, "l2rr", col = "white", col.regions = blue2white2red(10),
sp.layout = list(l2.atl), at = seq(-5, 5, 1), colorkey = list(space = "bottom",
labels = paste(labs)), sub = expression("Factor by which model frequency exceeds observed frequency"))
p2 = spplot(hspdfRR2, "l2rr", col = "white", col.regions = blue2white2red(10),
sp.layout = list(l2.atl), at = seq(-5, 5, 1), colorkey = list(space = "bottom",
labels = paste(labs)), sub = expression("Factor by which model frequency exceeds observed frequency"))
p1 = update(p1, main = textGrob("(a)", x = unit(0.05, "npc")), par.settings = list(fontsize = list(text = 14)))
p2 = update(p2, main = textGrob("(b)", x = unit(0.05, "npc")), par.settings = list(fontsize = list(text = 14)))
plot(p1, split = c(1, 1, 1, 2), more = TRUE)
plot(p2, split = c(1, 2, 1, 2), more = FALSE)
Qualitatively compare the GFDL and FSU-COAPS models
m = merge(z, z2, by = "hexid", all = TRUE)
m$groups = ifelse(m$rr.x > 1 & m$rr.y > 1, 4, ifelse(m$rr.x > 1 & m$rr.y < 1,
3, ifelse(m$rr.x < 1 & m$rr.y > 1, 2, ifelse(m$rr.x < 1 & m$rr.y < 1, 1,
0))))
m = subset(m, groups != "NA")
rownames(m) = m$hexid
hspdfCOMP = SpatialPolygonsDataFrame(hpg.atl[rownames(m)], m, match.ID = TRUE)
mycolors = c("blue", "cyan", "magenta", "red")
txt = list(c("FSU under, GFDL Under", "FSU Under, GFDL Over", "FSU Over, GFDL Under",
"FSU Over, GFDL Over"))
spplot(hspdfCOMP, "groups", col = "white", col.regions = mycolors, sp.layout = list(l2.atl),
colorkey = F, key = list(space = "bottom", rect = list(col = mycolors),
text = txt, cex = 0.9))
Mention using pressure (model) to estimate windspeeds using a p-w relationship:
load("G:/Research/world.interp.v02r01.R")
load("C:/Users/Sarah/Dropbox/GCMTrackData/gfdlPmin_r1.RData")
load("C:/Users/Sarah/Dropbox/GCMTrackData/gfdlPmin_r2.RData")
load("C:/Users/Sarah/Dropbox/GCMTrackData/gfdlPmin_r3.RData")
obs = subset(world.interp, Yr > 1981 & Yr < 2009)
r1 = subset(gfdlPmin_r1, Yr > 1981 & Yr < 2009)
r2 = subset(gfdlPmin_r2, Yr > 1981 & Yr < 2009)
r3 = subset(mod, Yr > 1981 & Yr < 2009)
obs$Wmax = 0.5144 * obs$Wmax
length(unique(obs$Sid)) # 2926
length(unique(r1$Sid)) # 2946
length(unique(r2$Sid)) # 2866
length(unique(r3$Sid)) # 2981
obs = subset(obs, Wmax > 17)
r1 = subset(r1, Wmax > 17)
r2 = subset(r2, Wmax > 17)
r3 = subset(r3, Wmax > 17)
length(unique(obs$Sid)) # 2512
length(unique(r1$Sid)) # 2945
length(unique(r2$Sid)) # 2863
length(unique(r3$Sid)) # 2981
(use data without taking 17m/s threshold to get “real” genesis point)
dat1 = ib.sGen@data
# dat1 = ib.sATL@data
dat1 = dat1[!duplicated(dat1$Sid), ]
dat1 = dat1[!is.na(dat1$hexid), ]
dat1 = ddply(dat1, .(hexid), summarize, genPts = length(unique(Sid)))
dat2 = fsu.sGen@data
# dat2 = fsu.sdf1@data
dat2 = dat2[!duplicated(dat2$Sid), ]
dat2 = dat2[!is.na(dat2$hexid), ]
dat2 = ddply(dat2, .(hexid), summarize, genPts = length(unique(Sid)))
dat3 = gfdl.sGen@data
# dat3 = gfdl.sna1@data
dat3 = dat3[!duplicated(dat3$Sid), ]
dat3 = dat3[!is.na(dat3$hexid), ]
dat3 = ddply(dat3, .(hexid), summarize, genPts = length(unique(Sid)))
gen = merge(dat1, dat2, by = "hexid", all = TRUE)
gen = merge(gen, dat3, by = "hexid", all = TRUE)
gen[is.na(gen)] = 0
colnames(gen) = c("hexid", "OBS", "FSU", "GFDL")
gen = gen[!gen$hexid == 106, ]
gen = gen[!gen$hexid == 170, ]
gen = gen[!gen$hexid == 212, ]
gen = gen[!gen$hexid == 44, ]
gen = gen[!gen$hexid == 22, ]
gen = gen[!gen$hexid == 23, ]
rownames(gen) = paste("ID", gen$hexid, sep = "")
hspdfGEN = SpatialPolygonsDataFrame(hpg.atl[rownames(gen)], gen, match.ID = TRUE)
cl.atl = map("world", xlim = c(-120, 10), ylim = c(0, 65), plot = FALSE, wrap = TRUE) #xlim,ylim reduces comp
clp.atl = map2SpatialLines(cl.atl, proj4string = CRS(ll))
clp.atl = spTransform(clp.atl, CRS(lcc))
l2.atl = list("sp.lines", clp.atl, lwd = 0.2)
centers = coordinates(hspdfGEN)
lids = list("sp.text", centers, hspdfGEN$hexid, col = "blue", cex = 1.1)
hspdfGEN$FSU[hspdfGEN$FSU > 14] = 16
cls = brewer.pal(8, "Oranges")
rng = seq(1, 17, 2)
p1 = spplot(hspdfGEN, "OBS", col = "white", col.regions = cls, sp.layout = list(l2.atl),
at = rng, xlim = bbox(hspdfGEN)[1, ], ylim = bbox(hspdfGEN)[2, ], colorkey = list(space = "bottom",
labels = c(rng[1:length(rng) - 1], "> 15")), sub = expression("Number of Genesis Points"))
p2 = spplot(hspdfGEN, "FSU", col = "white", col.regions = cls, sp.layout = list(l2.atl),
at = rng, xlim = bbox(hspdfGEN)[1, ], ylim = bbox(hspdfGEN)[2, ], colorkey = list(space = "bottom",
labels = c(rng[1:length(rng) - 1], "> 15")), sub = expression("Number of Genesis Points"))
p3 = spplot(hspdfGEN, "GFDL", col = "white", col.regions = cls, sp.layout = list(l2.atl),
at = rng, xlim = bbox(hspdfGEN)[1, ], ylim = bbox(hspdfGEN)[2, ], colorkey = list(space = "bottom",
labels = c(rng[1:length(rng) - 1], "> 15")), sub = expression("Number of Genesis Points"))
p1 = update(p1, main = textGrob("(a)", x = unit(0.05, "npc")), par.settings = list(fontsize = list(text = 14)))
p2 = update(p2, main = textGrob("(b)", x = unit(0.05, "npc")), par.settings = list(fontsize = list(text = 14)))
p3 = update(p3, main = textGrob("(c)", x = unit(0.05, "npc")), par.settings = list(fontsize = list(text = 14)))
plot(p1, split = c(1, 1, 1, 3), more = TRUE)
plot(p2, split = c(1, 2, 1, 3), more = TRUE)
plot(p3, split = c(1, 3, 1, 3), more = FALSE)
