require(sp)
require(rgdal)
require(maps)
require(maptools)
require(colorRamps)
require(grid)
require(RColorBrewer)
require(pROC)
require(ggplot2)
Data are from the International Best Track Archive for Climate Stewardship (IBTrACS). The 6-hourly data has been interpolated to hourly data.
First we subset on the years 1981–2008 and convert the windspeed to m/s
load("C:/Documents and Settings/sstrazzo/My Documents/Research/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)
We now create a spatial data frame from this data and assign a mollweide projection with a reference longitude of 150 degrees
ib.sdf = ib.df[, c("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"
ib.sdf = spTransform(ib.sdf, CRS(moll))
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, offset = c(1, -1))
hpg.i = HexPoints2SpatialPolygons(hpt.i)
We can 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] 6.39e+11
First, we calculate the per hexagon storm counts. Storm 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 storms contained in each hexagon and store this information in a polygon data frame, count.pdf.
ib.hexid = overlay(ib.sdf, hpg.i)
ib.split = split(ib.sdf@data, ib.hexid)
names(ib.split) = sapply(hpg.i@polygons, function(x) x@ID)[as.numeric(names(ib.split))]
ib.count = sapply(ib.split, function(x) length(table(x$Sid)))
ib.count = data.frame(ib.count)
count.pdf = SpatialPolygonsDataFrame(hpg.i[rownames(ib.count)], ib.count)
We next calculate the per hexagon storm hours. Since our original 6-hourly observations have been interpolated to hourly data, per hexagon storm hours are simply the number of observation points contained within each hexagon.
int.i = overlay(ib.sdf, hpg.i, fn = max)
int.i = data.frame(int.i)
rownam.i = row.names(int.i)
hn.i = as.integer(rownam.i)
hex.i = cbind(hn.i, int.i)
row.names(int.i) = paste("ID", hn.i, sep = "")
int.i = as.data.frame(int.i)
hspdf.i = SpatialPolygonsDataFrame(hpg.i[hn.i, ], int.i, match.ID = TRUE)
ib.sdf@data = data.frame(num = rep(1, ch.i))
co.i = overlay(ib.sdf, hspdf.i, fn = sum)
hspdf.i$hours = co.i[, 1]
Using the GCM-generated track data, we follow the same steps from above to calculate per hexagon storm counts and storm hours.
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))
gt.hexid = overlay(gt.sdf, hpg.i)
gt.split = split(gt.sdf@data, gt.hexid)
names(gt.split) = sapply(hpg.i@polygons, function(x) x@ID)[as.numeric(names(gt.split))]
gt.count = sapply(gt.split, function(x) length(table(x$Sid)))
gt.count = data.frame(gt.count)
countGT.pdf = SpatialPolygonsDataFrame(hpg.i[rownames(gt.count)], gt.count)
int.g = overlay(gt.sdf, hpg.i, fn = max)
int.g = data.frame(int.g)
rownam.g = row.names(int.g)
hn.g = as.integer(rownam.g)
hex.g = cbind(hn.g, int.g)
row.names(int.g) = paste("ID", hn.g, sep = "")
int.g = as.data.frame(int.g)
hspdf.g = SpatialPolygonsDataFrame(hpg.i[hn.g, ], int.g, match.ID = TRUE)
gt.sdf@data = data.frame(num = rep(1, ch.g))
co.g = overlay(gt.sdf, hspdf.g, fn = sum)
hspdf.g$hours = co.g[, 1]
l1.g = list("sp.points", gt.sdf, pch = 20, col = "black", cex = 0.2)
Using the GCM-generated track data, we follow the same steps from above to calculate per hexagon storm counts and storm hours.
load("track.df2.RData")
gt.df2 = subset(track.df2, Yr > 1981 & Yr < 2009)
gt.df2$Wmax = gt.df2$Wmax * 0.5144
ch.g2 = nrow(gt.df2)
gt.sdf2 = gt.df2[, c("lon", "lat", "Sid", "Sn", "Yr", "Mo", "Da", "Wmax", "DWmaxDt")]
coordinates(gt.sdf2) = c("lon", "lat")
proj4string(gt.sdf2) = CRS(ll)
gt.sdf2 = spTransform(gt.sdf2, CRS(moll))
gt.hexid2 = overlay(gt.sdf2, hpg.i)
gt.split2 = split(gt.sdf2@data, gt.hexid2)
names(gt.split2) = sapply(hpg.i@polygons, function(x) x@ID)[as.numeric(names(gt.split2))]
gt.count2 = sapply(gt.split2, function(x) length(table(x$Sid)))
gt.count2 = data.frame(gt.count2)
countGT.pdf2 = SpatialPolygonsDataFrame(hpg.i[rownames(gt.count2)], gt.count2)
int.g2 = overlay(gt.sdf2, hpg.i, fn = max)
int.g2 = data.frame(int.g2)
rownam.g2 = row.names(int.g2)
hn.g2 = as.integer(rownam.g2)
hex.g2 = cbind(hn.g2, int.g2)
row.names(int.g2) = paste("ID", hn.g2, sep = "")
int.g2 = as.data.frame(int.g2)
hspdf.g2 = SpatialPolygonsDataFrame(hpg.i[hn.g2, ], int.g2, match.ID = TRUE)
gt.sdf2@data = data.frame(num = rep(1, ch.g2))
co.g2 = overlay(gt.sdf2, hspdf.g2, fn = sum)
hspdf.g2$hours = co.g2[, 1]
l1.g2 = list("sp.points", gt.sdf2, pch = 20, col = "black", cex = 0.2)
Using the GCM-generated track data, we follow the same steps from above to calculate per hexagon storm counts and storm hours.
load("track.df3.RData")
gt.df3 = subset(track.df3, Yr > 1981 & Yr < 2009)
gt.df3$Wmax = gt.df3$Wmax * 0.5144
ch.g3 = nrow(gt.df3)
gt.sdf3 = gt.df3[, c("lon", "lat", "Sid", "Sn", "Yr", "Mo", "Da", "Wmax", "DWmaxDt")]
coordinates(gt.sdf3) = c("lon", "lat")
proj4string(gt.sdf3) = CRS(ll)
gt.sdf3 = spTransform(gt.sdf3, CRS(moll))
gt.hexid3 = overlay(gt.sdf3, hpg.i)
gt.split3 = split(gt.sdf3@data, gt.hexid3)
names(gt.split3) = sapply(hpg.i@polygons, function(x) x@ID)[as.numeric(names(gt.split3))]
gt.count3 = sapply(gt.split3, function(x) length(table(x$Sid)))
gt.count3 = data.frame(gt.count3)
countGT.pdf3 = SpatialPolygonsDataFrame(hpg.i[rownames(gt.count3)], gt.count3)
int.g3 = overlay(gt.sdf3, hpg.i, fn = max)
## Warning: Reached total allocation of 1535Mb: see help(memory.size)
## Warning: Reached total allocation of 1535Mb: see help(memory.size)
int.g3 = data.frame(int.g3)
rownam.g3 = row.names(int.g3)
hn.g3 = as.integer(rownam.g3)
hex.g3 = cbind(hn.g3, int.g3)
row.names(int.g3) = paste("ID", hn.g3, sep = "")
int.g3 = as.data.frame(int.g3)
hspdf.g3 = SpatialPolygonsDataFrame(hpg.i[hn.g3, ], int.g3, match.ID = TRUE)
gt.sdf3@data = data.frame(num = rep(1, ch.g3))
co.g3 = overlay(gt.sdf3, hspdf.g3, fn = sum)
## Warning: Reached total allocation of 1535Mb: see help(memory.size)
## Warning: Reached total allocation of 1535Mb: see help(memory.size)
## Warning: Reached total allocation of 1535Mb: see help(memory.size)
## Warning: Reached total allocation of 1535Mb: see help(memory.size)
hspdf.g3$hours = co.g3[, 1]
l1.g3 = list("sp.points", gt.sdf3, pch = 20, col = "black", cex = 0.2)
With the hexagons set up for both observed and modeled tracks, we can visually compare them.
We 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(6, "Reds")
p1 = spplot(hspdf.i, "hours", col = "white", col.regions = cr, sp.layout = list(l2.i),
at = seq(0, 12000, 2000), xlim = bbox(hspdf.i)[1, ], ylim = bbox(hspdf.i)[2,
], colorkey = list(space = "bottom", labels = paste(seq(0, 12000, 2000))),
sub = expression("Cyclone Hours"), main = "IBTrACS Cyclone Hours: 1981 - 2008")
p2 = spplot(hspdf.g, "hours", col = "white", col.regions = cr, sp.layout = list(l2.i),
at = seq(0, 12000, 2000), xlim = bbox(hspdf.i)[1, ], ylim = bbox(hspdf.i)[2,
], colorkey = list(space = "bottom", labels = paste(seq(0, 12000, 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, we map cyclone count hexagons.
cr2 = brewer.pal(6, "Reds")
p3 = spplot(count.pdf, "ib.count", col = "white", col.regions = cr2, sp.layout = list(l2.i),
at = seq(0, 300, 50), xlim = bbox(count.pdf)[1, ], ylim = bbox(count.pdf)[2,
], colorkey = list(space = "bottom", labels = paste(seq(0, 300, 50))),
sub = expression("Per hexagon cyclone count"), main = "IBTrACS Cyclone Counts: 1981 - 2008")
p4 = spplot(countGT.pdf, "gt.count", col = "white", col.regions = cr2, sp.layout = list(l2.i),
at = seq(0, 300, 50), xlim = bbox(count.pdf)[1, ], ylim = bbox(count.pdf)[2,
], colorkey = list(space = "bottom", labels = paste(seq(0, 300, 50))),
sub = expression("Per hexagon cyclone count"), 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)
For further qualitative comparison of the observed and modeled tracks, we create difference maps for cyclone hours and cyclone counts. The values in the hexagons represent observed cyclone hours/counts minus modeled cyclone hours/counts.