Additional Collaborators:
Tim LaRow (FSU-COAPS)
Michael Wehner (NCAR-CAM)
Ming Zhao (GFDL-HiRAM)
Hiroyuki Murakami (MRI-AGCM)

file:///C:/Users/Sarah/Dropbox/SpatialExtremes/Sarah/allModels.html

Introduction

Motivation

We know that tropical cyclones (TCs – including hurricanes, tropical storms, typhoons, and cyclones) can devastate coastal communities:

Given the destruction caused by hurricanes, many people are asking how climate change (e.g., warming sea surface waters) might affect TC activity.

Note that many of these studies use global climate models (GCMs) to make predictions about the effects of climate change on TC activity.

My question:
How well do climate models simulate the statistical relationship between TC intensity and sea surface temperatures (SST)?

If we want to use climate models to understand the future impacts of rising SSTs on TC intensity, we must first make sure that these models are able to reproduce the historical relationship between TC intensity and SST.

Here I examine the relationship between SST and TC intensity using observed and climate model data from the 1979–2010 time period.

The sensitivity of TC limiting intensity to SST: A model intercomparison

Preliminaries

date()
## [1] "Wed Apr 08 09:57:32 2015"

Set the working directory and make the required packages available. I also work from several R scripts to speed things up in certain places. In the future I should make these accessible via my website for transparency and reproducibility purposes!

# Set "eval=TRUE" when using PC
setwd("C:/Users/Sarah/Dropbox/SpatialExtremes")

require(ggplot2); require(rgdal); require(ismev); require(ncdf)
library(mapproj); library(mapdata); library(maptools); library(maps); library(grid); require(plyr)
library(colorRamps); library(RColorBrewer); require(spdep); require(psych); require(boot); require(wesanderson)
source("C:/Users/Sarah/Dropbox/SpatialExtremes/selectData.R")
source("C:/Users/Sarah/Dropbox/SpatialExtremes/ismevplus.R")
source("C:/Users/Sarah/Dropbox/SpatialExtremes/Kerrycanes/sensPoly.R")
source("C:/Users/Sarah/Dropbox/spatialextremes/ohcFunctions.R")
source("C:/Users/Sarah/Dropbox/SpatialExtremes/sensitivityUncertainty.R")

TC Data

Load Observed (best-track) and GCM data. Climate model data come from the following 4 models:

  • FSU-COAPS (100 km, 1982–2008)
  • GFDL-HiRAM (50 km, 1981–2009)
  • NCAR-CAM (25 km, 1979–2005)
  • MRI-AGCM (20 km, 1979–2010)

Subset on tropical storms and hurricanes and convert the wind speeds from knots to m/s. Finally, subtract 60% of the forward motion from the wind speed to approximate the rotational velocity (as per Emanuel (2008))

# Note: End in 2009 rather than 2010 (MRI only goes through 8/2010)
# load("~/Dropbox/SpatialExtremes/best.use.RData")
load("C:/Users/Sarah/Dropbox/SpatialExtremes/best.use.RData")
begin = 1979; end = 2010
Wind.df = best.use[best.use$Yr >= begin & best.use$Yr <= end, ]
Wind.df = subset(Wind.df, Type=="TS" | Type=="HU")
Wind.df$WmaxS = Wind.df$WmaxS * .5144
Wind.df$maguv = Wind.df$maguv * .5144
Wind.df$WmaxS = Wind.df$WmaxS - Wind.df$maguv * .6

#CAM
#load("~/Dropbox/GCMTrackData/CAMTracks_NAtl.RData")
load("C:/Users/Sarah/Dropbox/GCMTrackData/CAMTracks_NAtl.RData")
beginC = 1979; endC=2005
CAM.df = na.df1[na.df1$Yr >= beginC & na.df1$Yr <= endC, ]
CAM.df$WmaxS = CAM.df$WmaxS - CAM.df$maguv * .6

#GFDL-HiRAM (r1)
#load("~/Dropbox/SpatialExtremes/gfdlTracks_r1.df.RData")
load("C:/Users/Sarah/Dropbox/SpatialExtremes/gfdlTracks_r1.df.RData")
beginG = 1981; endG = 2009
GFDL.df = natltracks.df1[natltracks.df1$Yr >= beginG & 
                           natltracks.df1$Yr <= endG, ]
GFDL.df$WmaxS = GFDL.df$WmaxS - GFDL.df$maguv * .6

# FSU-COAPS (r1)
#load("~/Dropbox/SpatialExtremes/fsu_r1i1p1_ATL.RData")
load("C:/Users/Sarah/Dropbox/SpatialExtremes/fsu_r1i1p1_ATL.RData")
beginF = 1982; endF = 2009
FSU.df = fsuTracks.df[fsuTracks.df$Yr >= beginF &
                        fsuTracks.df$Yr <= endF, ]
FSU.df$WmaxS = FSU.df$WmaxS - FSU.df$maguv * .6

# MRI 
#load("~/Dropbox/SpatialExtremes/mriTracks.df.RData")
load("C:/Users/Sarah/Dropbox/SpatialExtremes/mriTracks.df.RData")
MRI.df = mriTracks.df[mriTracks.df$Yr >= begin &
                        mriTracks.df$Yr <= end, ]
MRI.df$WmaxS = MRI.df$WmaxS - MRI.df$maguv * .6

WindC.df = subset(Wind.df, Yr >= beginC & Yr <= endC)  # Obs TCs: 1979--2005
WindG.df = subset(Wind.df, Yr >= beginG & Yr <= endG)  # Obs TCs: 1981--2009
WindF.df = subset(Wind.df, Yr >= beginF & Yr <= endF)  # Obs TCs: 1982--2008

o.TCs = length(unique(Wind.df$Sid))  #372 (Obs)
m.TCs = length(unique(MRI.df$Sid))   #338 (MRI)
c.TCs = length(unique(CAM.df$Sid))   #387 (CAM)
g.TCs = length(unique(GFDL.df$Sid))  #418 (GFDL)
f.TCs = length(unique(FSU.df$Sid))   #392 (FSU)

Figure 1: Comparison of wind speed densities (top) and cumulative distributions (bottom) for best-track and GCM-generated TCs. Top: Frequency of hourly TC records. Rotational intensity shown in m/s. Bottom: Cumulative Distribution of hourly hurricane records

plot of chunk Fig1

Calculate the mean and skewness for each wind speed distribution.

mC = mean(CAM.df$WmaxS)    # 24.24
mM = mean(MRI.df$WmaxS)    # 27.35
mG = mean(GFDL.df$WmaxS)   # 18.02
mF = mean(FSU.df$WmaxS)    # 19.33
mO = mean(Wind.df$WmaxS)   # 28.83
sC = skew(CAM.df$WmaxS)    # 0.699 
sM = skew(MRI.df$WmaxS)    # 0.479
sG = skew(GFDL.df$WmaxS)   # 0.225
sF = skew(FSU.df$WmaxS)    # -0.0809
sO = skew(Wind.df$WmaxS)   # 1.02

SST Data

Load NOAA’s extended reconstructed SST data set and calculate average SSTs for each of the model time periods. Assign a Lambert conformal conic projection to the SST data (+proj) to create a spatial data frame of SST data.

#nc = open.ncdf("~/Dropbox/SpatialExtremes/sst.mnmean.nc")
nc = open.ncdf("C:/Users/Sarah/Dropbox/SpatialExtremes/sst.mnmean.nc")
avgsstM = reformat.SST(nc, begin, end, 8, 9, 10)
avgsstC = reformat.SST(nc, beginC, endC, 8, 9, 10)
avgsstG = reformat.SST(nc, beginG, endG, 8, 9, 10)
avgsstF = reformat.SST(nc, beginF, endF, 8, 9, 10)
rm(nc)

ll = "+proj=longlat +datum=WGS84"
lcc = "+proj=lcc +lat_1=60 +lat_2=30 +lon_0=-60"

coordinates(avgsstM) = c("lon", "lat")
proj4string(avgsstM) = CRS(ll)
sstM.sdf = spTransform(avgsstM, CRS(lcc))

coordinates(avgsstC) = c("lon", "lat")
proj4string(avgsstC) = CRS(ll)
sstC.sdf = spTransform(avgsstC, CRS(lcc))

coordinates(avgsstG) = c("lon", "lat")
proj4string(avgsstG) = CRS(ll)
sstG.sdf = spTransform(avgsstG, CRS(lcc))

coordinates(avgsstF) = c("lon", "lat")
proj4string(avgsstF) = CRS(ll)
sstF.sdf = spTransform(avgsstF, CRS(lcc))

Hexagons

As with the SST data, assign a Lambert Conformal Conic projection to the TC track data to create spatial data frames.

WindM.sdf = Wind.df
coordinates(WindM.sdf) = c("lon", "lat")
proj4string(WindM.sdf) = CRS(ll)
WindM.sdf = spTransform(WindM.sdf, CRS(lcc))

WindC.sdf = WindC.df
coordinates(WindC.sdf) = c("lon", "lat")
proj4string(WindC.sdf) = CRS(ll)
WindC.sdf = spTransform(WindC.sdf, CRS(lcc))

WindG.sdf = WindG.df
coordinates(WindG.sdf) = c("lon", "lat")
proj4string(WindG.sdf) = CRS(ll)
WindG.sdf = spTransform(WindG.sdf, CRS(lcc))

WindF.sdf = WindF.df
coordinates(WindF.sdf) = c("lon", "lat")
proj4string(WindF.sdf) = CRS(ll)
WindF.sdf = spTransform(WindF.sdf, CRS(lcc))

CAM.sdf = CAM.df
coordinates(CAM.sdf) = c("lon", "lat")
proj4string(CAM.sdf) = CRS(ll)
CAM.sdf = spTransform(CAM.sdf, CRS(lcc))

MRI.sdf = MRI.df
coordinates(MRI.sdf) = c("lon", "lat")
proj4string(MRI.sdf) = CRS(ll)
MRI.sdf = spTransform(MRI.sdf, CRS(lcc))

GFDL.sdf = GFDL.df
coordinates(GFDL.sdf) = c("lon", "lat")
proj4string(GFDL.sdf) = CRS(ll)
GFDL.sdf = spTransform(GFDL.sdf, CRS(lcc))

FSU.sdf = FSU.df
coordinates(FSU.sdf) = c("lon", "lat")
proj4string(FSU.sdf) = CRS(ll)
FSU.sdf = spTransform(FSU.sdf, CRS(lcc))

Define a set of equal area hexagons that cover all sets of track data. Set the random seed for a quasi-random starting location for the tessellation (useful for future testing).

set.seed(1025)
bb = matrix(0, nrow=2, ncol=2)
bb[2,1] = bbox(FSU.sdf)[2,1]*.8     
bb[1,1] = bbox(MRI.sdf)[1,1]*1.35
bb[1,2] = bbox(CAM.sdf)[1,2]*1.08
bb[2,2] = bbox(GFDL.sdf)[2,2]*1.08
hpt = spsample(CAM.sdf, type="hexagonal", n=150, bb=bb)
hpg = HexPoints2SpatialPolygons(hpt)
plot(hpg)
plot(WindM.sdf, add=TRUE, pch=".", col="red")
plot(CAM.sdf, add=TRUE, pch=".", col="blue")
plot(GFDL.sdf, add=TRUE, pch=".", col="green")
plot(FSU.sdf, add=TRUE, pch=".", col="gray")
plot(MRI.sdf, add=TRUE, pch=".", col="magenta")

plot of chunk tessellatebasinOBS

np = length(hpg@polygons)  # 122
area = hpg@polygons[[1]]@area * 10^-6

These hexagons are 816,355.3 km^2 – bigger than Texas but not quite as big as Alaska.

Next, use the “over” function to find out which hexagon each track point occurs in (do this for all five TC data sets). Add a “hexid” column to each of the track spatial data frames.

hexidOBSm = over(x=WindM.sdf, y=hpg)
WindM.sdf$hexid = hexidOBSm

hexidOBSc = over(x=WindC.sdf, y=hpg)
WindC.sdf$hexid = hexidOBSc

hexidOBSg = over(x=WindG.sdf, y=hpg)
WindG.sdf$hexid = hexidOBSg

hexidOBSf = over(x=WindF.sdf, y=hpg)
WindF.sdf$hexid = hexidOBSf

hexidCAM = over(x=CAM.sdf, y=hpg)
CAM.sdf$hexid = hexidCAM

hexidMRI = over(x=MRI.sdf, y=hpg)
MRI.sdf$hexid = hexidMRI

hexidGFDL = over(x=GFDL.sdf, y=hpg)
GFDL.sdf$hexid = hexidGFDL

hexidFSU = over(x=FSU.sdf, y=hpg)
FSU.sdf$hexid = hexidFSU

Use the “ddply” function to create data frames for each set of track data that consist of three columns: 1 column indicating the hexagon ID, another column indicating the number of TCs for that hexagon, and a third column indicating the maximum TC wind speed for that hexagon. Thus, the number of rows corresponds to the number of hexagons.

OBSstatsM = ddply(WindM.sdf@data, .(hexid, Sid), summarize,
                 nhr = length(Sid),
                 WmaxS = max(WmaxS)
                 )
OBSstatsM2 = ddply(OBSstatsM, .(hexid), summarize, 
                 countOBS = length(Sid),
                 WmaxOBS = max(WmaxS)
                  )

OBSstatsC = ddply(WindC.sdf@data, .(hexid, Sid), summarize,
                 nhr = length(Sid),
                 WmaxS = max(WmaxS)
                 )
OBSstatsC2 = ddply(OBSstatsC, .(hexid), summarize, 
                 countOBS = length(Sid),
                 WmaxOBS = max(WmaxS)
                  )

OBSstatsG = ddply(WindG.sdf@data, .(hexid, Sid), summarize,
                 nhr = length(Sid),
                 WmaxS = max(WmaxS)
                 )
OBSstatsG2 = ddply(OBSstatsG, .(hexid), summarize, 
                 countOBS = length(Sid),
                 WmaxOBS = max(WmaxS)
                  )


OBSstatsF = ddply(WindF.sdf@data, .(hexid, Sid), summarize,
                 nhr = length(Sid),
                 WmaxS = max(WmaxS)
                 )
OBSstatsF2 = ddply(OBSstatsF, .(hexid), summarize, 
                 countOBS = length(Sid),
                 WmaxOBS = max(WmaxS)
                  )
######################################################

CAMstats = ddply(CAM.sdf@data, .(hexid, Sid), summarize,
                 nhr = length(Sid),
                 WmaxS = max(WmaxS)
                 )
CAMstats2 = ddply(CAMstats, .(hexid), summarize, 
                 countMOD = length(Sid),
                 WmaxMOD = max(WmaxS)
                 )

#########################################################

MRIstats = ddply(MRI.sdf@data, .(hexid, Sid), summarize,
                 nhr = length(Sid),
                 WmaxS = max(WmaxS)
                 )
MRIstats2 = ddply(MRIstats, .(hexid), summarize, 
                 countMOD = length(Sid),
                 WmaxMOD = max(WmaxS)
                 )

#########################################################

GFDLstats = ddply(GFDL.sdf@data, .(hexid, Sid), summarize,
                 nhr = length(Sid),
                 WmaxS = max(WmaxS)
                 )
GFDLstats2 = ddply(GFDLstats, .(hexid), summarize, 
                 countMOD = length(Sid),
                 WmaxMOD = max(WmaxS)
                 )

#########################################################

FSUstats = ddply(FSU.sdf@data, .(hexid, Sid), summarize,
                 nhr = length(Sid),
                 WmaxS = max(WmaxS)
                 )
FSUstats2 = ddply(FSUstats, .(hexid), summarize, 
                 countMOD = length(Sid),
                 WmaxMOD = max(WmaxS)
                 )

Note: I’m not going to look at limiting intensity since it is underestimated over the time period of interest here. Instead, my goal is to understand the relationship between per hexagon maximum intensity and per hexagon average August–October SST.

Set up a data frame with counts, wind speeds, and hexid. Subset on only those hexagons with at least 15 TCs with maximum wind speeds >= 17 m/s. Also, remove any hexagons over land or over the East Pacific.

hexDataM.obs = subset(OBSstatsM2, WmaxOBS >= 17)     # 54 by 3
hexDataM.obs = subset(hexDataM.obs, countOBS >= 15)  # 32 by 3
hexDataM.mod = subset(MRIstats2, WmaxMOD >= 17)      # 53 by 3
hexDataM.mod = subset(hexDataM.mod, countMOD >= 15)  # 33 by 3
hexDataM.mod = subset(hexDataM.mod, hexid != 30)
row.names(hexDataM.obs) = paste("ID", hexDataM.obs$hexid, sep="")
row.names(hexDataM.mod) = paste("ID", hexDataM.mod$hexid, sep="")

##################################################################

hexDataC.obs = subset(OBSstatsC2, WmaxOBS >= 17)     # 52 by 3
hexDataC.obs = subset(hexDataC.obs, countOBS >= 15)  # 32 by 3
hexDataC.mod = subset(CAMstats2, WmaxMOD >= 17)      # 65 by 3
hexDataC.mod = subset(hexDataC.mod, countMOD >= 15)  # 47 by 3
hexDataC.mod = subset(hexDataC.mod, hexid != 16)
hexDataC.mod = subset(hexDataC.mod, hexid != 17)
hexDataC.mod = subset(hexDataC.mod, hexid != 18)
hexDataC.mod = subset(hexDataC.mod, hexid != 30)
hexDataC.mod = subset(hexDataC.mod, hexid != 31)
row.names(hexDataC.obs) = paste("ID", hexDataC.obs$hexid, sep="")
row.names(hexDataC.mod) = paste("ID", hexDataC.mod$hexid, sep="")

##################################################################

hexDataG.obs = subset(OBSstatsG2, WmaxOBS >= 17)     # 52 by 3
hexDataG.obs = subset(hexDataG.obs, countOBS >= 15)  # 32 by 3
hexDataG.mod = subset(GFDLstats2, WmaxMOD >= 17)      # 51 by 3
hexDataG.mod = subset(hexDataG.mod, countMOD >= 15)  # 39 by 3
hexDataG.mod = subset(hexDataG.mod, hexid != 17)
hexDataG.mod = subset(hexDataG.mod, hexid != 30)
hexDataG.mod = subset(hexDataG.mod, hexid != 31)
row.names(hexDataG.obs) = paste("ID", hexDataG.obs$hexid, sep="")
row.names(hexDataG.mod) = paste("ID", hexDataG.mod$hexid, sep="")

##################################################################

hexDataF.obs = subset(OBSstatsF2, WmaxOBS >= 17)     # 52 by 3
hexDataF.obs = subset(hexDataF.obs, countOBS >= 15)  # 32 by 3
hexDataF.mod = subset(FSUstats2, WmaxMOD >= 17)      # 55 by 3
hexDataF.mod = subset(hexDataF.mod, countMOD >= 15)  # 31 by 3
row.names(hexDataF.obs) = paste("ID", hexDataF.obs$hexid, sep="")
row.names(hexDataF.mod) = paste("ID", hexDataF.mod$hexid, sep="")

Make spatial polygons data frames for plotting purposes. Also, calculate the per hexagon average August-October SST.

hspdfM.obs = SpatialPolygonsDataFrame(hpg[row.names(hexDataM.obs)], 
                                      hexDataM.obs)
hspdfM.mod = SpatialPolygonsDataFrame(hpg[row.names(hexDataM.mod)], 
                                      hexDataM.mod)
avgsstM.obs = over(x=hpg[row.names(hexDataM.obs)], y=sstM.sdf, fn=mean)
avgsstM.mod = over(x=hpg[row.names(hexDataM.mod)], y=sstM.sdf, fn=mean)
hspdfM.obs$avgsst = avgsstM.obs$avgsst
hspdfM.mod$avgsst = avgsstM.mod$avgsst

hspdf.sst = subset(hspdfM.obs@data, avgsst >= 25)
hspdf.sst = SpatialPolygonsDataFrame(hpg[row.names(hspdf.sst)], hspdf.sst)

##################################################################

hspdfC.obs = SpatialPolygonsDataFrame(hpg[row.names(hexDataC.obs)], 
                                      hexDataC.obs)
hspdfC.mod = SpatialPolygonsDataFrame(hpg[row.names(hexDataC.mod)], 
                                      hexDataC.mod)
avgsstC.obs = over(x=hpg[row.names(hexDataC.obs)], y=sstC.sdf, fn=mean)
avgsstC.mod = over(x=hpg[row.names(hexDataC.mod)], y=sstC.sdf, fn=mean)
hspdfC.obs$avgsst = avgsstC.obs$avgsst
hspdfC.mod$avgsst = avgsstC.mod$avgsst

##################################################################

hspdfG.obs = SpatialPolygonsDataFrame(hpg[row.names(hexDataG.obs)], 
                                      hexDataG.obs)
hspdfG.mod = SpatialPolygonsDataFrame(hpg[row.names(hexDataG.mod)], 
                                      hexDataG.mod)
avgsstG.obs = over(x=hpg[row.names(hexDataG.obs)], y=sstG.sdf, fn=mean)
avgsstG.mod = over(x=hpg[row.names(hexDataG.mod)], y=sstG.sdf, fn=mean)
hspdfG.obs$avgsst = avgsstG.obs$avgsst
hspdfG.mod$avgsst = avgsstG.mod$avgsst

##################################################################

hspdfF.obs = SpatialPolygonsDataFrame(hpg[row.names(hexDataF.obs)], 
                                      hexDataF.obs)
hspdfF.mod = SpatialPolygonsDataFrame(hpg[row.names(hexDataF.mod)], 
                                      hexDataF.mod)
avgsstF.obs = over(x=hpg[row.names(hexDataF.obs)], y=sstF.sdf, fn=mean)
avgsstF.mod = over(x=hpg[row.names(hexDataF.mod)], y=sstF.sdf, fn=mean)
hspdfF.obs$avgsst = avgsstF.obs$avgsst
hspdfF.mod$avgsst = avgsstF.mod$avgsst

Set up map borders

cl = map("worldHires", xlim=c(-120, 20), ylim=c(-10, 70), plot=TRUE, resolution=.7)
clp = map2SpatialLines(cl, proj4string=CRS(ll))
clp = spTransform(clp, CRS(lcc))
l1 = list("sp.lines", clp, col="gray")

Figure 2: Per hexagon frequency and intensity maps for observed hexagons over the 1979-2010 (MRI) period).

plot of chunk Fig2

Figure 3: Frequency maps for grids with at least 15 observed AND modeled TCs. Red text indicates per hexagon observed TC counts

plot of chunk Fig3

Figure 4: Per hexagon maximum intensity maps for grids with at least 15 modeled TCs. Blue text indicates per hexagon observed maximum TC intensity

plot of chunk Fig4

Estimating the senstivity of maximum intensity to SST

Subset on only those hexagons with SST >= 25CC and regress per hexagon maximum intensity onto SST to obtain an estimate of the sensitivity of TC maximum intensity to SST (the slope coefficient). This first is done for the MRI model:

datM.obs = as.data.frame(hspdfM.obs)
datM.obs = subset(datM.obs, avgsst >= 25)
#datM.obs2 = as.data.frame(hspdfM.obs2)
#datM.obs2 = subset(datM.obs2, avgsst >= 25)
datM.mod = as.data.frame(hspdfM.mod)
datM.mod = subset(datM.mod, avgsst >= 25)

lmM.obs = lm(WmaxOBS ~ avgsst, data=datM.obs)
resM.obs = residuals(lmM.obs)
ciM.obs = confint(lmM.obs)   # 4.54, 8.69
## 6.6 +/- 1.01  
## t-val = 6.57
## R^2 = 0.6329

#lmM.obs2 = lm(WmaxOBS ~ avgsst, data=datM.obs2)
#ciM.obs2 = confint(lmM.obs2)   # 3.51, 8.74
## 6.1 +/- 1.25  
## t-val = 4.91
## R^2 = 0.5591

lmM.mod = lm(WmaxMOD ~ avgsst, data=datM.mod)  
resM.mod = residuals(lmM.mod)
ciM.mod = confint(lmM.mod)    #-2.67, 2.66
## -7.9e-4 +/- 1.29  
## t-val = -0.001
## R^2 = 1.651e-8

r2.M.obs = summary(lmM.obs)$r.squared
r2.M.mod = summary(lmM.mod)$r.squared

Same as previus chunk, but this time for the CAM:

datC.obs = as.data.frame(hspdfC.obs)
datC.obs = subset(datC.obs, avgsst >= 25)
datC.mod = as.data.frame(hspdfC.mod)
datC.mod = subset(datC.mod, avgsst >= 25)

lmC.obs = lm(WmaxOBS ~ avgsst, data=datC.obs)   
resC.obs = residuals(lmC.obs)
ciC.obs = confint(lmC.obs)   # 4.78, 9.12
## 6.9 +/- 1.05  
## t-val = 6.60
## R^2 = 0.6353

lmC.mod = lm(WmaxMOD ~ avgsst, data=datC.mod)
resC.mod = residuals(lmC.mod)
ciC.mod = confint(lmC.mod)   # -0.79, 5.07
## 2.1 +/- 1.43  
## t-val = 1.49
## R^2 = 0.07134

r2.C.obs = summary(lmC.obs)$r.squared
r2.C.mod = summary(lmC.mod)$r.squared

Same as previus chunk, but this time for the GFDL:

datG.obs = as.data.frame(hspdfG.obs)
datG.obs = subset(datG.obs, avgsst >= 25)
datG.mod = as.data.frame(hspdfG.mod)
datG.mod = subset(datG.mod, avgsst >= 25)

lmG.obs = lm(WmaxOBS ~ avgsst, data=datG.obs) 
resG.obs = residuals(lmG.obs)
ciG.obs = confint(lmG.obs)     #3.94, 8.45
## 6.2 +/- 1.10  
## t-val = 5.65
## R^2 = 0.5611

lmG.mod = lm(WmaxMOD ~ avgsst, data=datG.mod)
resG.mod = residuals(lmG.mod)
ciG.mod = confint(lmG.mod)   # -0.635, 2.21
## 7.9e-1 +/- 6.93e-1  
## t-val = 1.14
## R^2 = 0.04555

r2.G.obs = summary(lmG.obs)$r.squared
r2.G.mod = summary(lmG.mod)$r.squared

Same as previus chunk, but this time for the FSU:

datF.obs = as.data.frame(hspdfF.obs)
datF.obs = subset(datF.obs, avgsst >= 25)
datF.mod = as.data.frame(hspdfF.mod)
datF.mod = subset(datF.mod, avgsst >= 25)

lmF.obs = lm(WmaxOBS ~ avgsst, data=datF.obs)  
resF.obs = residuals(lmF.obs)
ciF.obs = confint(lmF.obs)   # 3.95, 8.49
## 6.2 +/- 1.10  
## t-val = 5.64
## R^2 = 0.5601

lmF.mod = lm(WmaxMOD ~ avgsst, data=datF.mod) 
resF.mod = residuals(lmF.mod)
ciF.mod = confint(lmF.mod)    # -3.01, 1.68
## -6.63e-1 +/- 1.13  
## t-val = -0.585
## R^2 = 0.01464

r2.F.obs = summary(lmF.obs)$r.squared
r2.F.mod = summary(lmF.mod)$r.squared

Figure 5: Comparing the sensitivity of LI to SST.

plot of chunk Fig5

Obtain the regression residuals

MRI.25.obs = subset(hspdfM.obs@data, avgsst >= 25)
MRI.25.mod = subset(hspdfM.mod@data, avgsst >= 25)
MRI.25.obs$res = resM.obs
MRI.25.mod$res = resM.mod
MRI.25.obs = SpatialPolygonsDataFrame(hpg[rownames(MRI.25.obs)], MRI.25.obs)
MRI.25.mod = SpatialPolygonsDataFrame(hpg[rownames(MRI.25.mod)], MRI.25.mod)

CAM.25.obs = subset(hspdfC.obs@data, avgsst >= 25)
CAM.25.mod = subset(hspdfC.mod@data, avgsst >= 25)
CAM.25.obs$res = resC.obs
CAM.25.mod$res = resC.mod
CAM.25.obs = SpatialPolygonsDataFrame(hpg[rownames(CAM.25.obs)], CAM.25.obs)
CAM.25.mod = SpatialPolygonsDataFrame(hpg[rownames(CAM.25.mod)], CAM.25.mod)

GFD.25.obs = subset(hspdfG.obs@data, avgsst >= 25)
GFD.25.mod = subset(hspdfG.mod@data, avgsst >= 25)
GFD.25.obs$res = resG.obs
GFD.25.mod$res = resG.mod
GFD.25.obs = SpatialPolygonsDataFrame(hpg[rownames(GFD.25.obs)], GFD.25.obs)
GFD.25.mod = SpatialPolygonsDataFrame(hpg[rownames(GFD.25.mod)], GFD.25.mod)

FSU.25.obs = subset(hspdfF.obs@data, avgsst >= 25)
FSU.25.mod = subset(hspdfF.mod@data, avgsst >= 25)
FSU.25.obs$res = resF.obs
FSU.25.mod$res = resF.mod
FSU.25.obs = SpatialPolygonsDataFrame(hpg[rownames(FSU.25.obs)], FSU.25.obs)
FSU.25.mod = SpatialPolygonsDataFrame(hpg[rownames(FSU.25.mod)], FSU.25.mod)

Calculate the pattern correlation of the residuals for the MRI vs. observed regressions:

mri.25.o = cbind(MRI.25.obs$hexid, 
               MRI.25.obs$res)
colnames(mri.25.o) = c("hexid", "res.o")
mri.25.o = data.frame(mri.25.o)

mri.25.m = cbind(MRI.25.mod$hexid, 
               MRI.25.mod$res)
colnames(mri.25.m) = c("hexid", "res.m")
mri.25.m = data.frame(mri.25.m)

mri.25.all = merge(mri.25.o, mri.25.m, by="hexid")
cor.test(mri.25.all$res.o, mri.25.all$res.m)
## 
##  Pearson's product-moment correlation
## 
## data:  mri.25.all$res.o and mri.25.all$res.m
## t = 3.454, df = 19, p-value = 0.00266
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2587 0.8302
## sample estimates:
##   cor 
## 0.621
# 0.259, 0.830 (0.621)

Same as previous chunk, but for the CAM:

cam.25.o = cbind(CAM.25.obs$hexid, 
               CAM.25.obs$res)
colnames(cam.25.o) = c("hexid", "res.o")
cam.25.o = data.frame(cam.25.o)

cam.25.m = cbind(CAM.25.mod$hexid, 
               CAM.25.mod$res)
colnames(cam.25.m) = c("hexid", "res.m")
cam.25.m = data.frame(cam.25.m)

cam.25.all = merge(cam.25.o, cam.25.m, by="hexid")
cor.test(cam.25.all$res.o, cam.25.all$res.m)
## 
##  Pearson's product-moment correlation
## 
## data:  cam.25.all$res.o and cam.25.all$res.m
## t = 1.982, df = 24, p-value = 0.059
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.01434  0.66572
## sample estimates:
##    cor 
## 0.3751
# -0.0143, 0.666 (0.375)

Same as previous chunk, but for the GFDL:

gfd.25.o = cbind(GFD.25.obs$hexid, 
               GFD.25.obs$res)
colnames(gfd.25.o) = c("hexid", "res.o")
gfd.25.o = data.frame(gfd.25.o)

gfd.25.m = cbind(GFD.25.mod$hexid, 
               GFD.25.mod$res)
colnames(gfd.25.m) = c("hexid", "res.m")
gfd.25.m = data.frame(gfd.25.m)

gfd.25.all = merge(gfd.25.o, gfd.25.m, by="hexid")
cor.test(gfd.25.all$res.o, gfd.25.all$res.m)
## 
##  Pearson's product-moment correlation
## 
## data:  gfd.25.all$res.o and gfd.25.all$res.m
## t = 2.817, df = 23, p-value = 0.009781
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1392 0.7513
## sample estimates:
##    cor 
## 0.5065
# 0.139, 0.751 (0.506)

Same as previous chunk, but for the FSU:

fsu.25.o = cbind(FSU.25.obs$hexid, 
               FSU.25.obs$res)
colnames(fsu.25.o) = c("hexid", "res.o")
fsu.25.o = data.frame(fsu.25.o)

fsu.25.m = cbind(FSU.25.mod$hexid, 
               FSU.25.mod$res)
colnames(fsu.25.m) = c("hexid", "res.m")
fsu.25.m = data.frame(fsu.25.m)

fsu.25.all = merge(fsu.25.o, fsu.25.m, by="hexid")
cor.test(fsu.25.all$res.o, fsu.25.all$res.m)
## 
##  Pearson's product-moment correlation
## 
## data:  fsu.25.all$res.o and fsu.25.all$res.m
## t = 3.116, df = 20, p-value = 0.005447
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1977 0.8004
## sample estimates:
##    cor 
## 0.5716
# 0.198, 0.800 (0.572)

Fig. 6 The spatial distribution of the residuals for the regression of PI onto SSTs > 25. The top panel shows residuals for MERRA data while the bottom panel shows residuals for FSU/COAPS data

plot of chunk Fig6

Relationship to latitude

Does latitude describe any of the variability in per hexagon maximum TC intensity?

Add hexagon center latitudes to data frames

datM.obs$lat = coordinates(MRI.25.obs)[,2]
datM.mod$lat = coordinates(MRI.25.mod)[,2]

datC.obs$lat = coordinates(CAM.25.obs)[,2]
datC.mod$lat = coordinates(CAM.25.mod)[,2]

datG.obs$lat = coordinates(GFD.25.obs)[,2]
datG.mod$lat = coordinates(GFD.25.mod)[,2]

datF.obs$lat = coordinates(FSU.25.obs)[,2]
datF.mod$lat = coordinates(FSU.25.mod)[,2]

Regress Wmax ~ avgsst + lat (MRI)

cor.test(datM.obs$avgsst, datM.obs$lat)
## 
##  Pearson's product-moment correlation
## 
## data:  datM.obs$avgsst and datM.obs$lat
## t = -3.406, df = 25, p-value = 0.002231
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.7768 -0.2328
## sample estimates:
##    cor 
## -0.563
# -0.23, -0.78
lmM.lat.obs = lm(WmaxOBS ~ avgsst + lat, data=datM.obs)
# avgsst: 7.4 +/1 1.21   ***    t-val = 6.11
# lat: 1.9e-6 +/-1.68e-6        t-val = 1.15
# R^2 = 0.652

cor.test(datM.mod$avgsst, datM.mod$lat)
## 
##  Pearson's product-moment correlation
## 
## data:  datM.mod$avgsst and datM.mod$lat
## t = -5.393, df = 23, p-value = 1.77e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.8820 -0.4997
## sample estimates:
##     cor 
## -0.7473
# -0.50, -0.88
lmM.lat.mod = lm(WmaxMOD ~ avgsst + lat, data=datM.mod)
# avgsst: 3.9 +/- 1.64     *    t-val = 2.37
# lat: 7.6e-6 +/-2.43e-6   **   t-val = 3.17
# R^2 = 0.3132

r2.latM.obs = summary(lmM.lat.obs)$r.squared
r2.latM.mod = summary(lmM.lat.mod)$r.squared

Regress Wmax ~ avgsst + lat (CAM)

cor.test(datC.obs$avgsst, datC.obs$lat)
# -0.23, -0.78
lmC.lat.obs = lm(WmaxOBS ~ avgsst + lat, data=datC.obs)
# avgsst: 8.0 +/- 1.23      ***    t-val = 6.51
# lat: 2.7e-6 +/- 1.70e-6          t-val = 1.58
# R^2 = 0.6695

cor.test(datC.mod$avgsst, datC.mod$lat)
# -0.25, -0.76
lmC.lat.mod = lm(WmaxMOD ~ avgsst + lat, data=datC.mod)
# avgsst: 4.2 +/- 1.62      *    t-val = 2.58
# lat: 4.9e-6 +/- 2.19e-6   *    t-val = 2.25
# R^2 = 0.2131

r2.latC.obs = summary(lmC.lat.obs)$r.squared
r2.latC.mod = summary(lmC.lat.mod)$r.squared

Regress Wmax ~ avgsst + lat (GFDL)

cor.test(datG.obs$avgsst, datG.obs$lat)
# -0.23, -0.78
lmG.lat.obs = lm(WmaxOBS ~ avgsst + lat, data=datG.obs)
# avgsst: 7.1 +/- 1.31      ***    t-val = 5.43
# lat: 2.2e-6 +/- 1.81e-6          t-val = 1.23
# R^2 = 0.5872

cor.test(datG.mod$avgsst, datG.mod$lat)
# -0.26, -0.78
lmG.lat.mod = lm(WmaxMOD ~ avgsst + lat, data=datG.mod)
# avgsst: 2.3 +/- 6.92e-1   **     t-val = 3.32
# lat: 3.8e-6 +/- 9.93e-7   ***    t-val = 3.80
# R^2 = 0.3862

r2.latG.obs = summary(lmG.lat.obs)$r.squared
r2.latG.mod = summary(lmG.lat.mod)$r.squared

Regress Wmax ~ avgsst + lat (FSU)

cor.test(datF.obs$avgsst, datF.obs$lat)
# -0.23, -0.77
lmF.lat.obs = lm(WmaxOBS ~ avgsst + lat, data=datF.obs)
# avgsst: 7.1 +/- 1.32      ***    t-val = 5.41
# lat: 2.2e-6 +/- 1.81e-6          t-val = 1.23
# R^2 = 0.5860

cor.test(datF.mod$avgsst, datF.mod$lat)
# -0.35, -0.83
lmF.lat.mod = lm(WmaxMOD ~ avgsst + lat, data=datF.mod)
# avgsst: 2.8 +/- 1.01      *      t-val = 2.80
# lat: 6.0e-6 +/- 1.13e-6   ***    t-val = 5.30
# R^2 = 0.5673

r2.latF.obs = summary(lmF.lat.obs)$r.squared
r2.latF.mod = summary(lmF.lat.mod)$r.squared

Improvements in explained variability by adding latitude as a predictor:

(R2Full - R2SSTonly)/R2SSTonly * 100%

## MRI
MRI.p.obs = r2.latM.obs*100 - r2.M.obs*100           # 1.91 %
MRI.p.mod = r2.latM.mod*100 - r2.M.mod*100           # 31.3 %

## CAM
CAM.p.obs = r2.latC.obs*100 - r2.C.obs*100           # 3.42 %
CAM.p.mod = r2.latC.mod*100 - r2.C.mod*100           # 14.2 % 

## GFDL
GFD.p.obs = r2.latG.obs*100 - r2.G.obs*100           # 2.61 %
GFD.p.mod = r2.latG.mod*100 - r2.G.mod*100           # 34.1 %

## FSU
FSU.p.obs = r2.latF.obs*100 - r2.F.obs*100           # 2.59 %
FSU.p.mod = r2.latF.mod*100 - r2.F.mod*100           # 55.3 %

Quantifying uncertainty associated with the spatial lattice

How does the specific starting point of the hexagon structure affect the results?

Generate 100 different tesselations, each slightly offset in space, and calculate the sensitivity for each. First, the CAM:

allSensCAM = matrix(0, ncol=2, nrow=100)
allSensCAM = data.frame(allSensCAM)

set.seed(1025)
bb = matrix(0, nrow=2, ncol=2)
bb[2,1] = bbox(FSU.sdf)[2,1]*.8     
bb[1,1] = bbox(MRI.sdf)[1,1]*1.35
bb[1,2] = bbox(CAM.sdf)[1,2]*1.08
bb[2,2] = bbox(GFDL.sdf)[2,2]*1.08

for(i in 1:100){
  print(i)
  hpts = spsample(CAM.sdf, type="hexagonal", n=150, bb=bb)
  hpg = HexPoints2SpatialPolygons(hpts)
  allSensCAM[i,] = get.sensitivity3(WindC.sdf, CAM.sdf, hpg, sstC.sdf, 1979, 2005)
}

colnames(allSensCAM) = c("OBS", "CAM")
range(allSensCAM$OBS)  ## 5.38, 11.3
range(allSensCAM$CAM)  ## -0.131, 5.04
median(allSensCAM$OBS) ##  7.88
median(allSensCAM$CAM) ##  1.59

save(allSensCAM, file="allSensCAM.RData")

Same as previous chunk, but for the MRI:

allSensMRI = matrix(0, ncol=2, nrow=100)
allSensMRI = data.frame(allSensMRI)

set.seed(1025)
bb = matrix(0, nrow=2, ncol=2)
bb[2,1] = bbox(FSU.sdf)[2,1]*.8     
bb[1,1] = bbox(MRI.sdf)[1,1]*1.35
bb[1,2] = bbox(CAM.sdf)[1,2]*1.08
bb[2,2] = bbox(GFDL.sdf)[2,2]*1.08

for(i in 1:100){
  print(i)
  hpts = spsample(MRI.sdf, type="hexagonal", n=150, bb=bb)
  hpg = HexPoints2SpatialPolygons(hpts)
  allSensMRI[i,] = get.sensitivity3(WindM.sdf, MRI.sdf, hpg, sstM.sdf, 1979, 2010)
}

colnames(allSensMRI) = c("OBS", "MRI")
range(allSensMRI$OBS)   # 5.30, 10.43
range(allSensMRI$MRI)   # -1.43, 3.30
median(allSensMRI$OBS)  # 7.41
median(allSensMRI$MRI)  # 0.223

save(allSensMRI, file="allSensMRI.RData")

Same as previous chunk, but for the GFDL:

allSensGFDL = matrix(0, ncol=2, nrow=100)
allSensGFDL = data.frame(allSensGFDL)

set.seed(1025)
bb = matrix(0, nrow=2, ncol=2)
bb[2,1] = bbox(FSU.sdf)[2,1]*.8     
bb[1,1] = bbox(MRI.sdf)[1,1]*1.35
bb[1,2] = bbox(CAM.sdf)[1,2]*1.08
bb[2,2] = bbox(GFDL.sdf)[2,2]*1.08

for(i in 1:100){
  print(i)
  hpts = spsample(GFDL.sdf, type="hexagonal", n=150, bb=bb)
  hpg = HexPoints2SpatialPolygons(hpts)
  allSensGFDL[i,] = get.sensitivity3(WindG.sdf, GFDL.sdf, hpg, sstG.sdf, 1981, 2009)
}

colnames(allSensGFDL) = c("OBS", "GFDL")
range(allSensGFDL$OBS)   ## 5.22, 10.44
range(allSensGFDL$GFDL)  ## 0.191, 2.75
median(allSensGFDL$OBS)  # 7.07
median(allSensGFDL$GFDL) # 1.06

save(allSensGFDL, file="allSensGFDL.RData")

Same as previous chunk, but for the FSU:

allSensFSU = matrix(0, ncol=2, nrow=100)
allSensFSU = data.frame(allSensFSU)

set.seed(1025)
bb = matrix(0, nrow=2, ncol=2)
bb[2,1] = bbox(FSU.sdf)[2,1]*.8     
bb[1,1] = bbox(MRI.sdf)[1,1]*1.35
bb[1,2] = bbox(CAM.sdf)[1,2]*1.08
bb[2,2] = bbox(GFDL.sdf)[2,2]*1.08

for(i in 1:100){
  print(i)
  hpts = spsample(FSU.sdf, type="hexagonal", n=150, bb=bb)
  hpg = HexPoints2SpatialPolygons(hpts)
  allSensFSU[i,] = get.sensitivity3(WindF.sdf, FSU.sdf, hpg, sstF.sdf, 1982, 2009)
}

colnames(allSensFSU) = c("OBS", "FSU")
range(allSensFSU$OBS)   # 5.12, 10.46
range(allSensFSU$FSU)   # -2.00, 0.590
median(allSensFSU$OBS)  # 7.10
median(allSensFSU$FSU)  # -0.632

save(allSensFSU, file="allSensFSU.RData")

Fig. 8: Estimating uncertainty on sensitivity resulting from the specific configuration of the hexagon grid. We use per hexagon maximum intensity in place of the theoretical LI

#load("~/Dropbox/SpatialExtremes/allSensCAM.RData")
#load("~/Dropbox/SpatialExtremes/allSensMRI.RData")
#load("~/Dropbox/SpatialExtremes/allSensGFDL.RData")
#load("~/Dropbox/SpatialExtremes/allSensFSU.RData")

load("C:/Users/Sarah/Dropbox/SpatialExtremes/allSensCAM.RData")
load("C:/Users/Sarah/Dropbox/SpatialExtremes/allSensMRI.RData")
load("C:/Users/Sarah/Dropbox/SpatialExtremes/allSensGFDL.RData")
load("C:/Users/Sarah/Dropbox/SpatialExtremes/allSensFSU.RData")


SC = c(allSensCAM$OBS, allSensCAM$CAM)
Label = c(rep("Best-track", 100), rep("GCM", 100))
SC.df = data.frame(Label = Label, Sens = SC)
SC.df$type = rep("(b) CAM", 200)
SC.df = data.frame(SC.df)

SM = c(allSensMRI$OBS, allSensMRI$MRI)
Label = c(rep("Best-track", 100), rep("GCM", 100))
SM.df = data.frame(Label = Label, Sens = SM)
SM.df$type = rep("(a) MRI", 200)
SM.df = data.frame(SM.df)

SG = c(allSensGFDL$OBS, allSensGFDL$GFDL)
Label = c(rep("Best-track", 100), rep("GCM", 100))
SG.df = data.frame(Label = Label, Sens = SG)
SG.df$type = rep("(c) GFDL", 200)
SG.df = data.frame(SG.df)

SF = c(allSensFSU$OBS, allSensFSU$FSU)
Label = c(rep("Best-track", 100), rep("GCM", 100))
SF.df = data.frame(Label = Label, Sens = SF)
SF.df$type = rep("(d) FSU", 200)
SF.df = data.frame(SF.df)

sens.df = rbind(SC.df, SM.df, SG.df, SF.df)
sens.df = data.frame(sens.df)

p = ggplot(sens.df, aes(Sens, fill = Label)) + theme_bw() 
p = p + geom_density(alpha=0.2) + facet_grid(~type)
p = p + xlab(expression(paste("Sensitivity of maximum intensity to SST [m s",{}^-1, " K", {}^-1,"]"))) + 
    ylab("Relative Frequency")
p = p + ggtitle('a') + theme(plot.title=element_text(hjust=2))
p = p + scale_fill_manual( values = c("blue","red"))
p = p + theme(axis.text=element_text(size=14), 
                    axis.title=element_text(size=16,face="bold"))
p = p + theme(legend.text=element_text(size=14),
                    legend.title=element_text(size=14, face="bold"))
p = p + theme(strip.text.x = element_text(size = 12, face="bold"))
p = p + theme(legend.title=element_blank())
p

plot of chunk Fig8