date()
## [1] "Fri Nov 30 15:07:04 2012"
setwd("~/Dropbox/GCMTrackData/Sarah")
require(sp)
require(rgdal)
require(maps)
require(maptools)
require(colorRamps)
require(grid)
require(RColorBrewer)
require(pROC)
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.
Observational data used in this research come from the International Best Track Archive for Climate Stewardship (IBTrACS, available online at http://www.ncdc.noaa.gov/oa/ibtracs/) (Knapp et al. 2010). For computational purposes, the 6-hourly data have been interpolated to hourly data using the method outlined in Elsner and Jagger (2012). Although IBTrACS includes more than a century's worth of track data, we subset the data from the period 1982–2008. This time period matches the time period over which the models were run. Additionally, because these data were obtained during the satellite era, we have more confidence in their accuracy.
Model-derived track data are obtained from a private website courtesy of the U.S. Climate Variability and Predictability Research Program (CLIVAR) Hurricane Working Group. We use data from two different atmospheric (non-couupled) GCMs. As with the observational data, the modeled track data are provided at 6-hourly intervals and have been interpolated to hourly intervals using the same method.
We first use model-derived track data from the GFDL HiRAM, version 2.2 (Zhao et al. 2012). The model data come from a control simulation forced with prescribed SSTs from the Hadley Centre Global Sea Ice and Sea Surface Temperature (Rayner et al. 2003). The model features a 50 km horizontal resolution and 32 vertical levels. As described by Zhao et al. (2009), TC-like vortices are detected and tracked using an algorithm similar to that used by Vitart et al. (2003). The algorithm searches for a coinciding (within 2 degrees latitude and longitude) relative vorticity maximum at 850 hPa, a sea-level pressure minimum, and a maximum in the 300-500 hPa averaged temperature field. The vortex trajectories were considered TC tracks when the modeled maximum surface winds exceeded 15.2 m/s during at least 2 (not necessarily consecutive) days (Zhao et al. 2009). We use track data from three separate runs of the HiRAM, referred to here as r1, r2, and r3.
In addition to the HiRAM track data, we also consider tracks derived from the FSU-COAPS global spectral model (Cocke and LaRow 2000; LaRow et al. 2008). As with the GFDL HiRAM, the FSU-COAPS model is not coupled and is forced with prescribed SSTs from the HadISST dataset. The spectral model has 27 vertical levels and a T126 horizontal resolution, which corresponds to roughly 0.94 degrees. The simulated TC tracks were also obtained using an algorithm based on Vitart et al. (2003). Because global track data are not yet available for this model, FSU-COAPS tracks are compared with observations and HiRAM tracks for the North Atlantic Basin only. Again, we use tracks from three model runs.
For the analysis, we convert both observational and model maximum wind speed measurements to meters per second. Additionally, we consider only those observations and modeled points that exceed 17 m/s (tropical storm strength winds).
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('/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)
We create a spatial data frame (ib.sdf) from the observational data (ib.df). Next, a planar coordinate system is assigned wtih a mollweide projection.
ib.sdf = ib.df[, c("name", "lon", "lat", "Sid", "Sn", "Yr", "Mo", "Da", "Wmax",
"DWmaxDt")]
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))
Additional details on the spatial framework used here can be found in Elsner et al. (2012). We first define a common grid of equal area hexagons that cover the global tropical and subtropical region to approximately 70 degrees N/S latitude. Each hexagon has an area of 727.9 thousand square kilometers, an area slightly larger than the state of Texas. The selected area is sufficiently small to capture regional variability. Once the hexagonal grid is defined, we populate the cells with track attributes from the observational and model datasets.
hpt.i = spsample(ib.sdf, type = "hexagonal", n = 1000, bb = bbox(ib.sdf) * 1.08,
offset = c(0, 0))
hpg.i = HexPoints2SpatialPolygons(hpt.i)
Extract the number of hexagons and the area of each hexagon.
np.i = length(hpg.i@polygons)
area.i = hpg.i@polygons[[1]]@area
area.i
## [1] 727859
The area is in square kilometers.
As we are interested in comparing observed and modeled spatial distribustion of storms, we calculate cyclone counts within the hexagonal grid framework. Cyclone counts represent the number of storms that “passed through” each hexagon in the spatial grid. For example, a hexagon centered in the western North Pacific may contain hundreds of storms for the period 1982–2008. Figure 1 illustrates the hexagonal framework for the year 2005 in the Gulf of Mexico. We sum the number of cyclones contained in each hexagon and store this information in a polygon data frame for further calculations. This procedure is carried out for both observed and modeled data.
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 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")
Defining hte hexagonal spatial data frame from observations, calculating per hexagon storm counts
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)
Calculating per hexagon storm counts for the HiRAM, r1
load("~/Dropbox/GCMTrackData/track.df1.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.g = nrow(gt.df1)
gt.sdf = gt.df1[, c("lon", "lat", "Sid", "Sn", "Yr", "Mo", "Da", "Wmax", "DWmaxDt")]
coordinates(gt.sdf) = c("lon", "lat")
proj4string(gt.sdf) = CRS(ll)
gt.sdf = spTransform(gt.sdf, CRS(moll))
l1.g1 = list("sp.points", gt.sdf, pch = 20, col = "black", cex = 0.2)
gt.hexid = over(gt.sdf, hpg.i)
gt.sdf$hexid = gt.hexid
gt.ch = ddply(gt.sdf@data, .(hexid), "nrow")
gt.cn1 = ddply(gt.sdf@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.ch$hexid, sep = "")
gt.chn1 = data.frame(nhours = gt.ch$nrow, ncyclones = gt.cn1$cn)
rownames(gt.chn1) = ID
hspdfHiRAMr1 = SpatialPolygonsDataFrame(hpg.i[ID], gt.chn1, match.ID = TRUE)
rm(track.df1)
Calculating per hexagon storm counts for the HiRAM, r2
load("~/Dropbox/GCMTrackData/track.df2.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.g = nrow(gt.df2)
gt.sdf = gt.df2[, c("lon", "lat", "Sid", "Sn", "Yr", "Mo", "Da", "Wmax", "DWmaxDt")]
coordinates(gt.sdf) = c("lon", "lat")
proj4string(gt.sdf) = CRS(ll)
gt.sdf = spTransform(gt.sdf, CRS(moll))
l1.g2 = list("sp.points", gt.sdf, pch = 20, col = "black", cex = 0.2)
gt.hexid = over(gt.sdf, hpg.i)
gt.sdf$hexid = gt.hexid
gt.ch = ddply(gt.sdf@data, .(hexid), "nrow")
gt.cn2 = ddply(gt.sdf@data, .(hexid), summarize, cn = length(unique(Sid)))
id.mod2 = gt.cn2$hexid # for difference mapping later
ID = paste("ID", gt.ch$hexid, sep = "")
gt.chn2 = data.frame(nhours = gt.ch$nrow, ncyclones = gt.cn2$cn)
rownames(gt.chn2) = ID
hspdfHiRAMr2 = SpatialPolygonsDataFrame(hpg.i[ID], gt.chn2, match.ID = TRUE)
rm(track.df2)
Calculating per hexagon storm counts for the HiRAM, r3
load("~/Dropbox/GCMTrackData/track.df3.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.g = nrow(gt.df3)
gt.sdf = gt.df3[, c("lon", "lat", "Sid", "Sn", "Yr", "Mo", "Da", "Wmax", "DWmaxDt")]
coordinates(gt.sdf) = c("lon", "lat")
proj4string(gt.sdf) = CRS(ll)
gt.sdf = spTransform(gt.sdf, CRS(moll))
l1.g3 = list("sp.points", gt.sdf, pch = 20, col = "black", cex = 0.2)
gt.hexid = over(gt.sdf, hpg.i)
gt.sdf$hexid = gt.hexid
gt.ch = ddply(gt.sdf@data, .(hexid), "nrow")
gt.cn3 = ddply(gt.sdf@data, .(hexid), summarize, cn = length(unique(Sid)))
id.mod3 = gt.cn3$hexid # for difference mapping later
ID = paste("ID", gt.ch$hexid, sep = "")
gt.chn3 = data.frame(nhours = gt.ch$nrow, ncyclones = gt.cn3$cn)
rownames(gt.chn3) = ID
hspdfHiRAMr3 = SpatialPolygonsDataFrame(hpg.i[ID], gt.chn3, match.ID = TRUE)
rm(track.df3)
mean(ib.df$Wmax)
## [1] 31.18
min(ib.df$Wmax)
## [1] 17
max(ib.df$Wmax)
## [1] 82.76
mean(gt.df1$Wmax)
## [1] 25.86
min(gt.df1$Wmax)
## [1] 17
max(gt.df1$Wmax)
## [1] 50.82
mean(gt.df2$Wmax)
## [1] 25.44
min(gt.df2$Wmax)
## [1] 17
max(gt.df2$Wmax)
## [1] 52.03
mean(gt.df3$Wmax)
## [1] 26.01
min(gt.df3$Wmax)
## [1] 17
max(gt.df3$Wmax)
## [1] 52.69
With the hexagons set up for both observed and modeled tracks, we can visually compare them. Figure 2 contains a map of cyclone counts for (a) observations and (b) the GFDL HiRAM. Darker reds indicate areas with higher storm frequencies over the 1982–2008 time period. In both maps, local maxima in TC counts are present in the western and eastern Pacific basins.
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)
To facilitate comparison of the observed and modeled TCs, we create difference maps. For cyclone count hexagons, this is accomplished by simply subtracting the values contained in the modeled storm count hexagons from those in the observed storm count hexagons.
z1 = merge(ib.cn, gt.cn1, by = "hexid", all = TRUE)
z1[is.na(z1)] = 0
z1$Diff = z1$cn.x - z1$cn.y
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)
Plotting difference maps
blue2white2red = colorRampPalette(c("blue", "white", "red"), space = "Lab")
p3 = spplot(hspdfDiffz1, "Diff", col = "white", col.regions = blue2white2red(8),
sp.layout = list(l2.i), at = seq(-100, 100, 25), xlim = bbox(hspdf)[1, ],
ylim = bbox(hspdf)[2, ], colorkey = list(space = "bottom", labels = paste(seq(-100,
100, 25))), sub = expression("Cyclone Count Difference"))
p4 = spplot(hspdfDiffz2, "Diff", col = "white", col.regions = blue2white2red(8),
sp.layout = list(l2.i), at = seq(-100, 100, 25), xlim = bbox(hspdf)[1, ],
ylim = bbox(hspdf)[2, ], colorkey = list(space = "bottom", labels = paste(seq(-100,
100, 25))), sub = expression("Cyclone Count Difference"))
p5 = spplot(hspdfDiffz3, "Diff", col = "white", col.regions = blue2white2red(8),
sp.layout = list(l2.i), at = seq(-100, 100, 25), xlim = bbox(hspdf)[1, ],
ylim = bbox(hspdf)[2, ], colorkey = list(space = "bottom", labels = paste(seq(-100,
100, 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)
From the difference maps, we can qualitatively asses model performance. We see that in some regions the model tends to generate too many storms, while in others, it generates fewer storms than are observed. To summarize model performance in a more quantitative manner, we use receiver operator character (ROC) analysis. ROC analysis offers a means of viewing the tradeoff between the false alarm rate and the hit rate of the model, with the fale alarm being the case that the model predicts a storm, but none actually occurs, and a hit being the case that the model correctly predicts a storm. Here, we consider the set of hexagons containing observational data to be the “truth” dataset. Because the observed and modeled hexagons are based on the same grid, we can easily implement a ROC analysis. For a hexagon populated with at least one observed storm, we test whether or not that hexagon is also populated by at least one modeled storm. A ROC curve is then created by varyng the threshold to find the optimal combination of hit rate and false alarm rate. For example, rather than using at least one observed storm per hexagon, we may be most likely to have a hit (a modeled storm and observed storm in the same hexagon) when the hexagon is populated with at least 3 or 4 observed storms. From the ROC analysis we obtain a ROC curve, which displays the tradeoff between false alarm rate (x-axis) and hit rate (y-axis). Ideally, we want the optimal false-alarm rate/hit rate combination to be in the upper left corner of the ROC space. In this ideal scenario, the area under the ROC curve is exactly one. If the area under the curve (AUC) is less than 50%, the model has no skill. However, if the AUC value is above 50%, the model has some skill.
Using the observational and GFDL-HiRAM data with the hexagon grid, we implement a ROC analysis for each year during the 1982–2008 time period. The ROC curve for 1982 is shown in Fig. X as an example. For a 7% false alarm rate, the model correctly identifies hexagons with observed activity as having activity 67% of the time. We are interested in quantitatively showing whether or not hexagons that contain observed storms also contain modeled storms, or stated differently, how well the observed and modeled storms match spatially. If the model does fairly well, the AUC values should be above 50%. Fig. X shows the AUC values plotted as a time series for the period 1982–2008. The AUC values are all higher than 50%, which indicates that at the resolution of this hexagon grid, the model matches the spatial distribution of observed storms fairly well.
(NOTE: The execution of the ROC analysis rely on several additional scripts: track_hexagons.R, which contains several useful functions for calculating per hexagon storm counts by year and for storing the yearly per hexagon storm counts in a matrix; and finally, get.AUC, which calculates the area under the ROC curve (with 95% confidence intervals) for the yearly per hexagon storm counts.)
# source('track_hexagons.R')# obs_counts = get.count.totals(ib.df, hpg.i) GFDL_counts1 =# get.count.totals(gt.df1, hpg.i) GFDL_counts2 = get.count.totals(gt.df2,# hpg.i) GFDL_counts3 = get.count.totals(gt.df3, hpg.i)# all_auc_r1 = get.auc(obs_counts, GFDL_counts1) all_auc_r1 =# cbind(1982:2008, all_auc_r1) colnames(all_auc_r1) = c('year', 'ci_lo',# 'auc', 'ci_hi') write.table(all_auc_r1, file='Gr1AUC_1000.txt')# all_auc_r2 = get.auc(obs_counts, GFDL_counts2) all_auc_r2 =# cbind(1982:2008, all_auc_r2) colnames(all_auc_r2) = c('year', 'ci_lo',# 'auc', 'ci_hi') write.table(all_auc_r2, file='Gr2AUC_1000.txt')# all_auc_r3 = get.auc(obs_counts, GFDL_counts3) all_auc_r3 =# cbind(1982:2008, all_auc_r3) colnames(all_auc_r3) = c('year', 'ci_lo',# 'auc', 'ci_hi') write.table(all_auc_r3, file='Gr3AUC_1000.txt')
obCounts = get.hex(ib.df, 1982, hpg.i)
obCounts = obCounts$count
moCounts = get.hex(gt.df1, 1982, hpg.i)
moCounts = moCounts$count
obCounts = obCounts > 0
rocEx = roc(obCounts, moCounts, percent = TRUE, ci = TRUE)
plot.roc(rocEx, legacy.axes = TRUE, print.auc = TRUE, xlab = "False Alarm Rate (%)",
ylab = "Hit Rate (%)", main = "ROC Analysis: HiRAM (r1) 1982")
##
## Call:
## roc.default(response = obCounts, predictor = moCounts, percent = TRUE, ci = TRUE)
##
## Data: moCounts in 708 controls (obCounts FALSE) < 174 cases (obCounts TRUE).
## Area under the curve: 84.8%
## 95% CI: 81.4%-88.2% (DeLong)
Plotting ROC AUC values with error bars using ggplot2
multiplot(p.a, p.b, p.c, cols = 1)
Because the results of the ROC analysis are partially dependent on the area of each cell within the hexagonal lattice, we test these results by varying the “resolution” (bad word?) of our hexagonal grid. The initial grid is obtained using a spatial sampling algorithm using n=1000 sampling points. This resulted in 882 hexagons each with an area of 727.9 thousand kilometers. We now compare the results for sampling the data using n=500, n=750, n=1500, and n=2000 points, which correspond to areas of 1,455.7 thousand, 970,5, 485.2, and 363.9 thousand square kilometers, repectively. We observe that while the differences in hexagon area do indeed affect the AUC values, these numbers nonetheless remain above the 50% threshold. As expected, smaller areas are associated with slightly lower AUC values, while larger areas area associated with higher AUC values. It should be noted that at this global scale, this type of ROC analysis is not as useful for diagnosing regional and subregional observation-model discrepancies. Rather, the ROC analysis is simply a means of quantitatively demonstrating how well the model generates a realistic global spatial distribution of TCs. These results indicate that, in general, the model is able to capture the known spatial distribution of storms.
# hex500.pts = spsample(ib.sdf, type='hexagonal', n=500,# bb=bbox(ib.sdf)*1.08, offset=c(0, 0)) hex500 =# HexPoints2SpatialPolygons(hex500.pts)# obs_counts500 = get.count.totals(ib.df, hex500) GFDL_counts500 =# get.count.totals(gt.df1, hex500) all_auc500 = get.auc(obs_counts500,# GFDL_counts500) all_auc500 = cbind(1982:2008, all_auc500)# colnames(all_auc500) = c('year', 'ci_lo', 'auc', 'ci_hi')# write.table(all_auc500, file='Gr1AUC_500.txt')
# hex750.pts = spsample(ib.sdf, type='hexagonal', n=750,# bb=bbox(ib.sdf)*1.08, offset=c(0, 0)) hex750 =# HexPoints2SpatialPolygons(hex750.pts)# obs_counts750 = get.count.totals(ib.df, hex750) GFDL_counts750 =# get.count.totals(gt.df1, hex750) all_auc750 = get.auc(obs_counts750,# GFDL_counts750) all_auc750 = cbind(1982:2008, all_auc750)# colnames(all_auc750) = c('year', 'ci_lo', 'auc', 'ci_hi')# write.table(all_auc750, file='Gr1AUC_750.txt')
# hex1500.pts = spsample(ib.sdf, type='hexagonal', n=1500,# bb=bbox(ib.sdf)*1.08, offset=c(0, 0)) hex1500 =# HexPoints2SpatialPolygons(hex1500.pts)# obs_counts1500 = get.count.totals(ib.df, hex1500) GFDL_counts1500 =# get.count.totals(gt.df1, hex1500) all_auc1500 = get.auc(obs_counts1500,# GFDL_counts1500) all_auc1500 = cbind(1982:2008, all_auc1500)# colnames(all_auc1500) = c('year', 'ci_lo', 'auc', 'ci_hi')# write.table(all_auc1500, file='Gr1AUC_1500.txt')
# hex2000.pts = spsample(ib.sdf, type='hexagonal', n=2000,# bb=bbox(ib.sdf)*1.08, offset=c(0, 0)) hex2000 =# HexPoints2SpatialPolygons(hex2000.pts)# obs_counts2000 = get.count.totals(ib.df, hex2000) GFDL_counts2000 =# get.count.totals(gt.df1, hex2000) all_auc2000 = get.auc(obs_counts2000,# GFDL_counts2000) all_auc2000 = cbind(1982:2008, all_auc2000)# colnames(all_auc2000) = c('year', 'ci_lo', 'auc', 'ci_hi')# write.table(all_auc2000, file='Gr1AUC_2000.txt')
Plotting ROC AUC time series with varying resolution hexagonal lattices
multiplot(p.500, p.750, p.1000, p.1500, p.2000, cols = 2)
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.
# 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
percentile.atl = length(ib.atl$Wmax[ib.atl$Wmax <= 17])/length(ib.atl$Wmax)
ib.atl = subset(ib.atl, Wmax >= 17)
ib.sATL = ib.atl[, c("name", "lon", "lat", "Sid", "Sn", "Yr", "Mo", "Da", "Wmax",
"DWmaxDt")]
coordinates(ib.sATL) = c("lon", "lat")
ll = "+proj=longlat +datum=WGS84"
proj4string(ib.sATL) = CRS(ll)
# moll.atl = '+proj=moll +lon_0=360 +units=km'
lcc = "+proj=lcc +lat_1=60 +lat_2=30 +lon_0=-60 +units=km"
# ib.sATL = spTransform(ib.sATL, CRS(moll.atl))
ib.sATL = spTransform(ib.sATL, CRS(lcc))
hpt.atl = spsample(ib.sATL, type = "hexagonal", n = 500, bb = bbox(ib.sATL) *
1.08, offset = c(0, 0), iter = 100)
hpg.atl = HexPoints2SpatialPolygons(hpt.atl)
ib.idATL = over(ib.sATL, hpg.atl)
ib.sATL$hexid = ib.idATL
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)
## repeat above for FSU tracks (r1, r2, & r3)
load("~/Dropbox/GCMTrackData/fsu_r1i1p1_ATL.RData")
fsu.r1 = fsuTracks.df
rm(fsuTracks.df)
fsu.r1 = subset(fsu.r1, Yr > 1981 & Yr < 2009)
fsu.r1 = subset(fsu.r1, Wmax >= 17)
fsu.sdf1 = fsu.r1[, c("name", "lon", "lat", "Sid", "Sn", "Yr", "Mo", "Da", "Wmax",
"DWmaxDt")]
coordinates(fsu.sdf1) = c("lon", "lat")
proj4string(fsu.sdf1) = CRS(ll)
# fsu.sdf1 = spTransform(fsu.sdf1, CRS(moll.atl))
fsu.sdf1 = spTransform(fsu.sdf1, CRS(lcc))
fsu.hexid1 = over(fsu.sdf1, hpg.atl)
fsu.sdf1$hexid = fsu.hexid1
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)
## repeat above for FSU tracks (r1, r2, & r3)
load("~/Dropbox/GCMTrackData/fsu_r2i1p1_ATL.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.sdf2, CRS(moll.atl))
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)
## repeat above for FSU tracks (r1, r2, & r3)
load("~/Dropbox/GCMTrackData/fsu_r3i1p1_ATL.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.sdf3, CRS(moll.atl))
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)
## repeat above for GFDL tracks (r1, r2, & r3)
load("~/Dropbox/GCMTrackData/gfdlTracks_r1.df.RData")
gfdl.na1 = natltracks.df1
rm(natltracks.df1)
gfdl.na1 = subset(gfdl.na1, Yr > 1981 & Yr < 2009)
gfdl.na1 = subset(gfdl.na1, Wmax >= 17)
gfdl.sna1 = gfdl.na1[, c("name", "lon", "lat", "Sid", "Sn", "Yr", "Mo", "Da",
"Wmax", "DWmaxDt")]
coordinates(gfdl.sna1) = c("lon", "lat")
proj4string(gfdl.sna1) = CRS(ll)
# gfdl.sna1 = spTransform(gfdl.sna1, CRS(moll.atl))
gfdl.sna1 = spTransform(gfdl.sna1, CRS(lcc))
gfdl.hexid1 = over(gfdl.sna1, hpg.atl)
gfdl.sna1$hexid = gfdl.hexid1
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)
## repeat above for GFDL tracks (r1, r2, & r3)
load("~/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)
## repeat above for GFDL tracks (r1, r2, & r3)
load("~/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)
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)
From this figure it is apparent that there are some discrepancies between the FSU-COAPs model and observations. The FSU-COAPS model appears to generate more storms in the southwestern portion of the north Atlantic. Figure Xb also suggests that most of these model-generated storms recurve to the north before impacting the United States. This is consistent with LaRow et al. (2008) who attribute this to the model's large-scale steering flow during the first half of the hurricane season. The framework also allows us to directly compare the FSU-COAPS model with the HiRAM. In contrast to the FSU-COAPS model, which tends to overgenerate storms in the north Atlantic, the GFDL HiRAM does not generate as many storms as are observed. The spatial pattern of the HiRAM storms is more consistent with observations, although there are fewer modeled storms in the Gulf of Mexico and Caribbean than are observed. This latter point is also true of the FSU-COAPS model.
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)
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)
If we consider only those hexagons which contain storm data for both the FSU-COAPS model and observations, we can compute some useful descriptive statistics. For example, we compute the correlation between between the modeled number of cyclones and observed number of cyclones for the hexagonal grid. For the north Atlantic basin, the correlation is about 0.47.
Using this common grid, we can also examine the factor by which the modeled storm frequency exceeds that of the observed frequency. We first calculate the relative ratios for both obserations and model by dividing the number of cyclones in each hexagon by the total number of cyclones for the entire grid. We then divide the the model relative ratio by the observed relative ratio to obtain the factor by which modeled storm frequency exceeds observed storm frequency. These logarithm of these factors are shown in Figure ?. The colorbar shows the logarithm of these factors. Positive numbers indicate overprediction by the model, while negative numbers indicate underprediction. The largest exceedance factor is 6.3 and the smallest is 0.03. Of the 198 hexagons in the grid, 129 contained more observed than modeled storms, while 69 contained more modeled than observed storms. The overprediction region is clearly visible in the center of FigX. It is also apparent that few modeled storms are present in the Gulf of Mexico, at higher latitudes, and near the Cape Verde Islands.
Qualitatively comparing 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))
