date()
## [1] "Tue Aug 7 16:02:40 2012"
require(sp)
require(rgdal)
require(maps)
require(maptools)
require(colorRamps)
require(grid)
require(RColorBrewer)
require(pROC)
require(ggplot2)
require(plyr)
Data are from the International Best Track Archive for Climate Stewardship (IBTrACS). The 6-hourly data has been interpolated to hourly data using the method outlined in Elsner and Jagger (2012).
First subset on the years 1981–2008 and convert the wind speed to m/s.
# 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
ch.i = nrow(ib.df)
Next create a spatial data frame from this data and assign a mollweide projection with a reference longitude of 150 degrees and spatial units of kilometers.
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 +units=km"
ib.sdf = spTransform(ib.sdf, CRS(moll))
Extract three cyclones.
ib.samples = subset(ib.sdf@data, Yr == 1985 & name == "KATE" | Yr ==
1992 & name == "ANDREW" | Yr == 2005 & name == "KATRINA")
The R functions spsample and HexPoints2SpatialPolygons are used to first create a grid of points and then connect those points to form a grid of hexagons. Track data from both observations and GCMs will be overlaid onto this hexagon grid.
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] 745377
The area is in square kilometers.
Cyclone counts represent the number of storms that “passed through” each hexagon from the grid created above for the period 1981–2008. For example, since TCs don't form over the equator, a hexagon centered over the equator may contain no data. Conversely, a hexagon centered in the West Pacific may contain hundreds of storms. We sum the number of cyclones contained in each hexagon and store this information in a polygon data frame hspdf.
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 = 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)
Using the GCM-generated track data, we follow the same steps from above to calculate per hexagon cyclone counts and cyclone hours.
setwd("~/Dropbox/GCMTrackData")
load("track.df1.RData")
gt.df = subset(track.df1, Yr > 1981 & Yr < 2009)
gt.df$Wmax = gt.df$Wmax * 0.5144
ch.g = nrow(gt.df)
gt.sdf = gt.df[, 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.cn = ddply(gt.sdf@data, .(hexid), summarize, cn = length(unique(Sid)))
ID = paste("ID", gt.ch$hexid, sep = "")
gt.chn = data.frame(nhours = gt.ch$nrow, ncyclones = gt.cn$cn)
rownames(gt.chn) = ID
hspdfHiRAMr1 = SpatialPolygonsDataFrame(hpg.i[ID], gt.chn, match.ID = TRUE)
load("track.df2.RData")
gt.df = subset(track.df2, Yr > 1981 & Yr < 2009)
gt.df$Wmax = gt.df$Wmax * 0.5144
ch.g = nrow(gt.df)
gt.sdf = gt.df[, 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.cn = ddply(gt.sdf@data, .(hexid), summarize, cn = length(unique(Sid)))
ID = paste("ID", gt.ch$hexid, sep = "")
gt.chn = data.frame(nhours = gt.ch$nrow, ncyclones = gt.cn$cn)
rownames(gt.chn) = ID
hspdfHiRAMr2 = SpatialPolygonsDataFrame(hpg.i[ID], gt.chn, match.ID = TRUE)
load("track.df3.RData")
gt.df = subset(track.df3, Yr > 1981 & Yr < 2009)
gt.df$Wmax = gt.df$Wmax * 0.5144
ch.g = nrow(gt.df)
gt.sdf = gt.df[, 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.cn = ddply(gt.sdf@data, .(hexid), summarize, cn = length(unique(Sid)))
ID = paste("ID", gt.ch$hexid, sep = "")
gt.chn = data.frame(nhours = gt.ch$nrow, ncyclones = gt.cn$cn)
rownames(gt.chn) = ID
hspdfHiRAMr3 = SpatialPolygonsDataFrame(hpg.i[ID], gt.chn, match.ID = TRUE)
With the hexagons set up for both observed and modeled tracks, we can visually compare them. First map cyclone hours hexagons.
l1.i = list("sp.points", ib.sdf, pch = 20, col = "black", cex = 0.1)
cl.i = map("world", plot = FALSE, wrap = TRUE)
clp.i = map2SpatialLines(cl.i, proj4string = CRS(ll))
clp.i = spTransform(clp.i, CRS(moll))
l2.i = list("sp.lines", clp.i, lwd = 0.2)
cr = brewer.pal(7, "Reds")
p1 = spplot(hspdf, "nhours", col = "white", col.regions = cr, sp.layout = list(l2.i),
at = seq(0, 14000, 2000), xlim = bbox(hspdf)[1, ], ylim = bbox(hspdf)[2,
], colorkey = list(space = "bottom", labels = paste(seq(0, 14000, 2000))),
sub = expression("Cyclone Hours"), main = "IBTrACS Cyclone Hours: 1981-2008")
p2 = spplot(hspdfHiRAMr1, "nhours", col = "white", col.regions = cr,
sp.layout = list(l2.i), at = seq(0, 14000, 2000), xlim = bbox(hspdfHiRAMr1)[1,
], ylim = bbox(hspdfHiRAMr1)[2, ], colorkey = list(space = "bottom",
labels = paste(seq(0, 14000, 2000))), sub = expression("Cyclone Hours"),
main = "GFDL-HiRAM (R1) Cyclone Hours: 1981-2008")
print(p1, split = c(1, 1, 1, 2), more = TRUE)
print(p2, split = c(1, 2, 1, 2), more = FALSE)
Next map cyclone count hexagons.
cr2 = brewer.pal(7, "Reds")
p3 = spplot(hspdf, "ncyclones", col = "white", col.regions = cr2,
sp.layout = list(l2.i), at = seq(0, 350, 50), xlim = bbox(hspdf)[1, ], ylim = bbox(hspdf)[2,
], colorkey = list(space = "bottom", labels = paste(seq(0, 350, 50))),
sub = expression("Number of Cyclones"), main = "IBTrACS Cyclone Counts: 1981-2008")
p4 = spplot(hspdfHiRAMr1, "ncyclones", col = "white", col.regions = cr2,
sp.layout = list(l2.i), at = seq(0, 350, 50), xlim = bbox(hspdfHiRAMr1)[1,
], ylim = bbox(hspdfHiRAMr1)[2, ], colorkey = list(space = "bottom",
labels = paste(seq(0, 350, 50))), sub = expression("Number of Cyclones"),
main = "GFDL-HiRAM (R1) Cyclone Counts: 1981-2008")
print(p3, split = c(1, 1, 1, 2), more = TRUE)
print(p4, split = c(1, 2, 1, 2), more = FALSE)
