The influence of model resolution on the simulated sensitivity of North Atlantic tropical cyclone maximum intensity to sea surface temperature

---

---

Introduction

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

Here we 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] "Mon Jan 25 17:56:34 2016"

Set the working directory and make the required packages available. Note: I also work from several R scripts to speed things up in certain places. Please contact Sarah Strazzo (s.e.strazz@gmail.com) for these scripts.

setwd("~/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("~/Dropbox/SpatialExtremes/selectData.R")
source("~/Dropbox/SpatialExtremes/ismevplus.R")
source("~/Dropbox/SpatialExtremes/Kerrycanes/sensPoly.R")
source("~/Dropbox/spatialextremes/ohcFunctions.R")
source("~/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)
• CAM5 0.25 (25 km, 1979–2005)
• CAM5 (1 degree, 1979–2007)
• MRI-AGCM (20 km, 1979–2009)

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")
begin = 1979; end = 2009
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 Quarter Degree
load("~/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

#CAM 1 Degree
load("~/Dropbox/GCMTrackData/CAM1deg.df.RData")
beginC1 = 1979; endC1=2007
CAM1.df = cam1degTracks[cam1degTracks$Yr >= beginC1 & cam1degTracks$Yr <= endC1, ]
CAM1.df$WmaxS = CAM1.df$WmaxS - CAM1.df$maguv * .6

#GFDL-HiRAM (r1)
load("~/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")
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")
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
WindC1.df = subset(Wind.df, Yr >= beginC1 & Yr <= endC1)  
# Obs TCs: 1979--2007
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)
c1d.TCs = length(unique(CAM1.df$Sid))  # 35 (CAM 1 degree)
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

---

---

Calculate the mean and skewness for each wind speed distribution.

mC = mean(CAM.df$WmaxS)    # 24.24
mC1 = mean(CAM1.df$WmaxS) # 17.10
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 
sC1 = skew(CAM1.df$WmaxS) # 0.138
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")
avgsstM = reformat.SST(nc, begin, end, 8, 9, 10)
avgsstC = reformat.SST(nc, beginC, endC, 8, 9, 10)
avgsstC1 = reformat.SST(nc, beginC1, endC1, 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(avgsstC1) = c("lon", "lat")
proj4string(avgsstC1) = CRS(ll)
sstC1.sdf = spTransform(avgsstC1, 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))

WindC1.sdf = WindC1.df
coordinates(WindC1.sdf) = c("lon", "lat")
proj4string(WindC1.sdf) = CRS(ll)
WindC1.sdf = spTransform(WindC1.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))

CAM1.sdf = CAM1.df
coordinates(CAM1.sdf) = c("lon", "lat")
proj4string(CAM1.sdf) = CRS(ll)
CAM1.sdf = spTransform(CAM1.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(CAM1.sdf, add=TRUE, pch=".", col="cyan")
#plot(GFDL.sdf, add=TRUE, pch=".", col="green")
#plot(FSU.sdf, add=TRUE, pch=".", col="gray")
#plot(MRI.sdf, add=TRUE, pch=".", col="magenta")
np = length(hpg@polygons)  # 122
area = hpg@polygons[[1]]@area * 10^-6 
#save(hpg, file="hexagons.RData")

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

hexidOBSc1 = over(x=WindC1.sdf, y=hpg)
WindC1.sdf$hexid = hexidOBSc1

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

hexidCAM1 = over(x=CAM1.sdf, y=hpg)
CAM1.sdf$hexid = hexidCAM1

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)
                  )

OBSstatsC1 = ddply(WindC1.sdf@data, .(hexid, Sid), summarize,
                 nhr = length(Sid),
                 WmaxS = max(WmaxS)
                 )
OBSstatsC12 = ddply(OBSstatsC1, .(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)
                 )

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


CAM1stats = ddply(CAM1.sdf@data, .(hexid, Sid), summarize,
                 nhr = length(Sid),
                 WmaxS = max(WmaxS)
                 )
CAM1stats2 = ddply(CAM1stats, .(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. On second thought, omit 15 TCs threshold given the lack of TC activity in the model. Plus, we’re not actually using LI here, so it doesn’t make sense to limit the number of TCs per hexagon.

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 != 3)
hexDataM.mod = subset(hexDataM.mod, hexid != 4)
hexDataM.mod = subset(hexDataM.mod, hexid != 5)
hexDataM.mod = subset(hexDataM.mod, hexid != 29)
hexDataM.mod = subset(hexDataM.mod, hexid != 43)
hexDataM.mod = subset(hexDataM.mod, hexid != 57)
hexDataM.mod = subset(hexDataM.mod, hexid != 16)
hexDataM.mod = subset(hexDataM.mod, hexid != 17)
hexDataM.mod = subset(hexDataM.mod, hexid != 18)
hexDataM.mod = subset(hexDataM.mod, hexid != 30)
hexDataM.mod = subset(hexDataM.mod, hexid != 31)

hexDataM.obs = subset(hexDataM.obs, hexid != 3)
hexDataM.obs = subset(hexDataM.obs, hexid != 4)
hexDataM.obs = subset(hexDataM.obs, hexid != 5)
hexDataM.obs = subset(hexDataM.obs, hexid != 29)
hexDataM.obs = subset(hexDataM.obs, hexid != 43)
hexDataM.obs = subset(hexDataM.obs, hexid != 57)
hexDataM.obs = subset(hexDataM.obs, hexid != 16)
hexDataM.obs = subset(hexDataM.obs, hexid != 17)
hexDataM.obs = subset(hexDataM.obs, hexid != 18)
hexDataM.obs = subset(hexDataM.obs, hexid != 30)
hexDataM.obs = subset(hexDataM.obs, hexid != 31)

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)     # 54 by 3
hexDataC.obs = subset(hexDataC.obs, countOBS >= 15)  # 32 by 3
hexDataC.mod = subset(CAMstats2, WmaxMOD >= 17)      # 53 by 3
hexDataC.mod = subset(hexDataC.mod, countMOD >= 15)  # 33 by 3

hexDataC.mod = subset(hexDataC.mod, hexid != 3)
hexDataC.mod = subset(hexDataC.mod, hexid != 4)
hexDataC.mod = subset(hexDataC.mod, hexid != 5)
hexDataC.mod = subset(hexDataC.mod, hexid != 29)
hexDataC.mod = subset(hexDataC.mod, hexid != 43)
hexDataC.mod = subset(hexDataC.mod, hexid != 57)
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)

hexDataC.obs = subset(hexDataC.obs, hexid != 3)
hexDataC.obs = subset(hexDataC.obs, hexid != 4)
hexDataC.obs = subset(hexDataC.obs, hexid != 5)
hexDataC.obs = subset(hexDataC.obs, hexid != 29)
hexDataC.obs = subset(hexDataC.obs, hexid != 43)
hexDataC.obs = subset(hexDataC.obs, hexid != 57)
hexDataC.obs = subset(hexDataC.obs, hexid != 16)
hexDataC.obs = subset(hexDataC.obs, hexid != 17)
hexDataC.obs = subset(hexDataC.obs, hexid != 18)
hexDataC.obs = subset(hexDataC.obs, hexid != 30)
hexDataC.obs = subset(hexDataC.obs, hexid != 31)

row.names(hexDataC.obs) = paste("ID", hexDataC.obs$hexid, sep="")
row.names(hexDataC.mod) = paste("ID", hexDataC.mod$hexid, sep="")

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

hexDataC1.obs = subset(OBSstatsC12, WmaxOBS >= 17)     # 54 by 3
hexDataC1.obs = subset(hexDataC1.obs, countOBS >= 15)  # 32 by 3
hexDataC1.mod = subset(CAM1stats2, WmaxMOD >= 17)      # 53 by 3
#hexDataC1.mod = subset(hexDataC1.mod, countMOD >= 15)  # 33 by 3

hexDataC1.mod = subset(hexDataC1.mod, hexid != 3)
hexDataC1.mod = subset(hexDataC1.mod, hexid != 4)
hexDataC1.mod = subset(hexDataC1.mod, hexid != 5)
hexDataC1.mod = subset(hexDataC1.mod, hexid != 29)
hexDataC1.mod = subset(hexDataC1.mod, hexid != 43)
hexDataC1.mod = subset(hexDataC1.mod, hexid != 57)
hexDataC1.mod = subset(hexDataC1.mod, hexid != 16)
hexDataC1.mod = subset(hexDataC1.mod, hexid != 17)
hexDataC1.mod = subset(hexDataC1.mod, hexid != 18)
hexDataC1.mod = subset(hexDataC1.mod, hexid != 30)
hexDataC1.mod = subset(hexDataC1.mod, hexid != 31)

hexDataC1.obs = subset(hexDataC1.obs, hexid != 3)
hexDataC1.obs = subset(hexDataC1.obs, hexid != 4)
hexDataC1.obs = subset(hexDataC1.obs, hexid != 5)
hexDataC1.obs = subset(hexDataC1.obs, hexid != 29)
hexDataC1.obs = subset(hexDataC1.obs, hexid != 43)
hexDataC1.obs = subset(hexDataC1.obs, hexid != 57)
hexDataC1.obs = subset(hexDataC1.obs, hexid != 16)
hexDataC1.obs = subset(hexDataC1.obs, hexid != 17)
hexDataC1.obs = subset(hexDataC1.obs, hexid != 18)
hexDataC1.obs = subset(hexDataC1.obs, hexid != 30)
hexDataC1.obs = subset(hexDataC1.obs, hexid != 31)

row.names(hexDataC1.obs) = paste("ID", hexDataC1.obs$hexid, sep="")
row.names(hexDataC1.mod) = paste("ID", hexDataC1.mod$hexid, sep="")

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

hexDataG.obs = subset(OBSstatsG2, WmaxOBS >= 17)     # 54 by 3
hexDataG.obs = subset(hexDataG.obs, countOBS >= 15)  # 32 by 3
hexDataG.mod = subset(GFDLstats2, WmaxMOD >= 17)      # 53 by 3
hexDataG.mod = subset(hexDataG.mod, countMOD >= 15)  # 33 by 3

hexDataG.mod = subset(hexDataG.mod, hexid != 3)
hexDataG.mod = subset(hexDataG.mod, hexid != 4)
hexDataG.mod = subset(hexDataG.mod, hexid != 5)
hexDataG.mod = subset(hexDataG.mod, hexid != 29)
hexDataG.mod = subset(hexDataG.mod, hexid != 43)
hexDataG.mod = subset(hexDataG.mod, hexid != 57)
hexDataG.mod = subset(hexDataG.mod, hexid != 16)
hexDataG.mod = subset(hexDataG.mod, hexid != 17)
hexDataG.mod = subset(hexDataG.mod, hexid != 18)
hexDataG.mod = subset(hexDataG.mod, hexid != 30)
hexDataG.mod = subset(hexDataG.mod, hexid != 31)

hexDataG.obs = subset(hexDataG.obs, hexid != 3)
hexDataG.obs = subset(hexDataG.obs, hexid != 4)
hexDataG.obs = subset(hexDataG.obs, hexid != 5)
hexDataG.obs = subset(hexDataG.obs, hexid != 29)
hexDataG.obs = subset(hexDataG.obs, hexid != 43)
hexDataG.obs = subset(hexDataG.obs, hexid != 57)
hexDataG.obs = subset(hexDataG.obs, hexid != 16)
hexDataG.obs = subset(hexDataG.obs, hexid != 17)
hexDataG.obs = subset(hexDataG.obs, hexid != 18)
hexDataG.obs = subset(hexDataG.obs, hexid != 30)
hexDataG.obs = subset(hexDataG.obs, 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)     # 54 by 3
hexDataF.obs = subset(hexDataF.obs, countOBS >= 15)  # 32 by 3
hexDataF.mod = subset(FSUstats2, WmaxMOD >= 17)      # 53 by 3
hexDataF.mod = subset(hexDataF.mod, countMOD >= 15)  # 33 by 3

hexDataF.mod = subset(hexDataF.mod, hexid != 3)
hexDataF.mod = subset(hexDataF.mod, hexid != 4)
hexDataF.mod = subset(hexDataF.mod, hexid != 5)
hexDataF.mod = subset(hexDataF.mod, hexid != 29)
hexDataF.mod = subset(hexDataF.mod, hexid != 43)
hexDataF.mod = subset(hexDataF.mod, hexid != 57)
hexDataF.mod = subset(hexDataF.mod, hexid != 16)
hexDataF.mod = subset(hexDataF.mod, hexid != 17)
hexDataF.mod = subset(hexDataF.mod, hexid != 18)
hexDataF.mod = subset(hexDataF.mod, hexid != 30)
hexDataF.mod = subset(hexDataF.mod, hexid != 31)

hexDataF.obs = subset(hexDataF.obs, hexid != 3)
hexDataF.obs = subset(hexDataF.obs, hexid != 4)
hexDataF.obs = subset(hexDataF.obs, hexid != 5)
hexDataF.obs = subset(hexDataF.obs, hexid != 29)
hexDataF.obs = subset(hexDataF.obs, hexid != 43)
hexDataF.obs = subset(hexDataF.obs, hexid != 57)
hexDataF.obs = subset(hexDataF.obs, hexid != 16)
hexDataF.obs = subset(hexDataF.obs, hexid != 17)
hexDataF.obs = subset(hexDataF.obs, hexid != 18)
hexDataF.obs = subset(hexDataF.obs, hexid != 30)
hexDataF.obs = subset(hexDataF.obs, hexid != 31)

row.names(hexDataF.obs) = paste("ID", hexDataF.obs$hexid, sep="")
row.names(hexDataF.mod) = paste("ID", hexDataF.mod$hexid, sep="")

#remove: 58, 59, 44, 31, 32, 19
hexDataM.obs2 = subset(hexDataM.obs, hexid != 58)
hexDataM.obs2 = subset(hexDataM.obs2, hexid != 59)
hexDataM.obs2 = subset(hexDataM.obs2, hexid != 44)
hexDataM.obs2 = subset(hexDataM.obs2, hexid != 32)
hexDataM.obs2 = subset(hexDataM.obs2, hexid != 19)
hspdfM.obs2 = SpatialPolygonsDataFrame(hpg[row.names(hexDataM.obs2)], 
                                      hexDataM.obs2)
avgsstM.obs2 = over(x=hpg[row.names(hexDataM.obs2)], 
                    y=sstM.sdf, fn=mean)
hspdfM.obs2$avgsst = avgsstM.obs2$avgsst

---

---

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

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

hspdfC1.obs = SpatialPolygonsDataFrame(hpg[row.names(hexDataC1.obs)], 
                                      hexDataC1.obs)
hspdfC1.mod = SpatialPolygonsDataFrame(hpg[row.names(hexDataC1.mod)], 
                                      hexDataC1.mod)
avgsstC1.obs = over(x=hpg[row.names(hexDataC1.obs)], 
                    y=sstC1.sdf, fn=mean)
avgsstC1.mod = over(x=hpg[row.names(hexDataC1.mod)], 
                    y=sstC1.sdf, fn=mean)
hspdfC1.obs$avgsst = avgsstC1.obs$avgsst
hspdfC1.mod$avgsst = avgsstC1.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).

#centers1 = coordinates(hspdfM.obs)
#centers2 = coordinates(hspdf.sst)
#l2 = list("sp.text", centers1, hspdfM.obs$hexid, col="red", cex=1.1)
#l3 = list("sp.text", centers2, hspdf.sst$hexid, col="black", cex=1.1)

#crCount = brewer.pal(8, "Blues")     # Count
crCount = blue2yellow(8)
rngCount = seq(0, 120, 15)          # Count

crWind = wes_palette(name = "Zissou", n=7, type = "continuous") 
rngWind = seq(17, 87, 10)           # Wmax

crSST = brewer.pal(5, "YlOrRd")     # SST
rngSST = seq(25, 30, 1)             # SST


p0 = spplot(hspdfM.obs, "countOBS", col="white", col.regions=crCount,
            sp.layout=list(l1), xlim=bb[1, ],
            ylim=c(bb[2, 1]*-1.1, bb[2, 2]),
            at=rngCount, 
            colorkey=list(space="bottom", 
            labels=c(rngCount[1:length(rngCount)-1], '>105')),
            sub=expression("Observed TC Frequency"))

p1 = spplot(hspdfM.obs, "WmaxOBS", col="white", col.regions=crWind,
            sp.layout=list(l1), xlim=bb[1, ],
            ylim=c(bb[2, 1]*-1.1, bb[2, 2]),
            at=rngWind, 
            colorkey=list(space="bottom", labels=paste(rngWind)),
            sub=expression(paste("Observed Maximum Intensity (m"~s^-1~")")))

p2 = spplot(hspdf.sst, "avgsst", col="white", col.regions=crSST,
            sp.layout=list(l1), xlim=bb[1, ],
            ylim=c(bb[2, 1]*-1.1, bb[2, 2]),
            at=rngSST, 
            colorkey=list(space="bottom", labels=paste(rngSST)),
            sub=expression(paste("Mean August--October SST (",degree,"C)")))


p0 = update(p0, main = textGrob("(a)", x = unit(0.05, "npc")),
            par.settings = list(fontsize = list(text = 14)))
p1 = update(p1, main = textGrob("(b)", x = unit(0.05, "npc")), 
            par.settings = list(fontsize = list(text = 14)))
p2 = update(p2, main = textGrob("(c)", x = unit(0.05, "npc")), 
            par.settings = list(fontsize = list(text = 14)))

plot(p0, split=c(1, 1, 1, 3), more=TRUE)
plot(p1, split=c(1, 2, 1, 3), more=TRUE)
plot(p2, split=c(1, 3, 1, 3), more=FALSE)

---

---

Figure 3: Frequency maps. Red text indicates per hexagon observed TC counts

#centers.M = coordinates(hspdfM.mod)
#centers.C = coordinates(hspdfC.mod)
#centers.C1 = coordinates(hspdfC1.mod)
#centers.G = coordinates(hspdfG.mod)
#centers.F = coordinates(hspdfF.mod)

#l.M = list("sp.text", centers.M, hspdfM.mod$hexid, col="red", cex=1.1)
#l.C = list("sp.text", centers.C, hspdfC.mod$hexid, col="red", cex=1.1)
#l.C1 = list("sp.text", centers.C1, hspdfC1.mod$hexid, col="red", cex=1.1)
#l.G = list("sp.text", centers.G, hspdfG.mod$hexid, col="red", cex=1.1)
#l.F = list("sp.text", centers.F, hspdfF.mod$hexid, col="red", cex=1.1)

hspdfplot = hspdfF.mod
hspdfplot$countMOD[hspdfplot$countMOD >= 120] = 115

p1 = spplot(hspdfM.mod, "countMOD", col="white", col.regions=crCount,
            sp.layout=list(l1), xlim=bb[1, ],
            ylim=c(bb[2, 1]*-1.1, bb[2, 2]),
            at=rngCount, 
            colorkey=list(space="bottom", 
            labels=c(rngCount[1:length(rngCount)-1], '> 105')),
            sub=expression("TC Frequency: 1979-2009"))

p2 = spplot(hspdfC.mod, "countMOD", col="white", col.regions=crCount,
            sp.layout=list(l1), xlim=bb[1, ],
            ylim=c(bb[2, 1]*-1.1, bb[2, 2]),
            at=rngCount, 
            colorkey=list(space="bottom", 
            labels=c(rngCount[1:length(rngCount)-1], '> 105')),
            sub=expression("TC Frequency: 1979-2005"))

p3 = spplot(hspdfG.mod, "countMOD", col="white", col.regions=crCount,
            sp.layout=list(l1), xlim=bb[1, ],
            ylim=c(bb[2, 1]*-1.1, bb[2, 2]),
            at=rngCount, 
            colorkey=list(space="bottom", 
            labels=c(rngCount[1:length(rngCount)-1], '> 105')),
            sub=expression("TC Frequency: 1981-2009"))

p4 = spplot(hspdfplot, "countMOD", col="white", col.regions=crCount,
            sp.layout=list(l1), xlim=bb[1, ],
            ylim=c(bb[2, 1]*-1.1, bb[2, 2]),
            at=rngCount, 
            colorkey=list(space="bottom", 
            labels=c(rngCount[1:length(rngCount)-1], '> 105')),
            sub=expression("TC Frequency: 1982-2008"))

p5 = spplot(hspdfC1.mod, "countMOD", col="white", col.regions=crCount,
            sp.layout=list(l1), xlim=bb[1, ],
            ylim=c(bb[2, 1]*-1.1, bb[2, 2]),
            at=rngCount, 
            colorkey=list(space="bottom", 
            labels=c(rngCount[1:length(rngCount)-1], '> 105')),
            sub=expression("TC Frequency: 1979-2008"))


p1 = update(p1, main = textGrob("(a) MRI", x = unit(0.05, "npc")),
            par.settings = list(fontsize = list(text = 14)))
p2 = update(p2, main = textGrob("(b) CAM5 0.25", x = unit(0.05, "npc")), 
            par.settings = list(fontsize = list(text = 14)))
p3 = update(p3, main = textGrob("(c) GFDL", x = unit(0.05, "npc")), 
            par.settings = list(fontsize = list(text = 14)))
p4 = update(p4, main = textGrob("(d) FSU", x = unit(0.05, "npc")), 
            par.settings = list(fontsize = list(text = 14)))
p5 = update(p5, main = textGrob("(e) CAM5 1", x = unit(0.05, "npc")), 
            par.settings = list(fontsize = list(text = 14)))

plot(p1, split=c(1, 1, 2, 3), more=TRUE)
plot(p2, split=c(2, 1, 2, 3), more=TRUE)
plot(p3, split=c(1, 2, 2, 3), more=TRUE)
plot(p4, split=c(2, 2, 2, 3), more=TRUE)
plot(p5, split=c(1, 3, 2, 3), more=FALSE)

rm(hspdfplot)

---

---

Figure 4: Per hexagon maximum intensity maps. Blue text indicates per hexagon observed maximum TC intensity

p1 = spplot(hspdfM.mod, "WmaxMOD", col="white", col.regions=crWind,
            sp.layout=list(l1), xlim=bb[1, ],
            ylim=c(bb[2, 1]*-1.1, bb[2, 2]),
            at=rngWind, 
            colorkey=list(space="bottom", labels=paste(rngWind)),
            sub=expression(paste("Maximum Intensity (m "~s^-1~")")))

p2 = spplot(hspdfC.mod, "WmaxMOD", col="white", col.regions=crWind,
            sp.layout=list(l1), xlim=bb[1, ],
            ylim=c(bb[2, 1]*-1.1, bb[2, 2]),
            at=rngWind, 
            colorkey=list(space="bottom", labels=paste(rngWind)),
            sub=expression(paste("Maximum Intensity (m "~s^-1~")")))

p3 = spplot(hspdfG.mod, "WmaxMOD", col="white", col.regions=crWind,
            sp.layout=list(l1), xlim=bb[1, ],
            ylim=c(bb[2, 1]*-1.1, bb[2, 2]),
            at=rngWind, 
            colorkey=list(space="bottom", labels=paste(rngWind)),
            sub=expression(paste("Maximum Intensity (m "~s^-1~")")))

p4 = spplot(hspdfF.mod, "WmaxMOD", col="white", col.regions=crWind,
            sp.layout=list(l1), xlim=bb[1, ],
            ylim=c(bb[2, 1]*-1.1, bb[2, 2]), 
            at=rngWind, 
            colorkey=list(space="bottom", labels=paste(rngWind)),
            sub=expression(paste("Maximum Intensity (m "~s^-1~")")))

p5 = spplot(hspdfC1.mod, "WmaxMOD", col="white", col.regions=crWind,
            sp.layout=list(l1), xlim=bb[1, ],
            ylim=c(bb[2, 1]*-1.1, bb[2, 2]),
            at=rngWind, 
            colorkey=list(space="bottom", labels=paste(rngWind)),
            sub=expression(paste("Maximum Intensity (m "~s^-1~")")))

p1 = update(p1, main = textGrob("(a) MRI", x = unit(0.05, "npc")),
            par.settings = list(fontsize = list(text = 14)))
p2 = update(p2, main = textGrob("(b) CAM5 0.25", x = unit(0.05, "npc")), 
            par.settings = list(fontsize = list(text = 14)))
p3 = update(p3, main = textGrob("(c) GFDL", x = unit(0.05, "npc")), 
            par.settings = list(fontsize = list(text = 14)))
p4 = update(p4, main = textGrob("(d) FSU", x = unit(0.05, "npc")), 
            par.settings = list(fontsize = list(text = 14)))
p5 = update(p5, main = textGrob("(e) CAM5 1", x = unit(0.05, "npc")), 
            par.settings = list(fontsize = list(text = 14)))

plot(p1, split=c(1, 1, 2, 3), more=TRUE)
plot(p2, split=c(2, 1, 2, 3), more=TRUE)
plot(p3, split=c(1, 2, 2, 3), more=TRUE)
plot(p4, split=c(2, 2, 2, 3), more=TRUE)
plot(p5, split=c(1, 3, 2, 3), more=FALSE)

---

---

---

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.31, 8.64
## 6.47 +/- 1.05  
## t-val = 6.17
## R^2 = 0.6136

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

lmM.mod = lm(WmaxMOD ~ avgsst, data=datM.mod)  
resM.mod = residuals(lmM.mod)
ciM.mod = confint(lmM.mod)    #-2.67, 2.66
## -7.94e-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.54, 9.06
## 6.80 +/- 1.10  
## t-val = 6.21
## R^2 = 0.616

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 CAM 1 degree:

datC1.obs = as.data.frame(hspdfC1.obs)
datC1.obs = subset(datC1.obs, avgsst >= 25)
datC1.mod = as.data.frame(hspdfC1.mod)
datC1.mod = subset(datC1.mod, avgsst >= 25)

lmC1.obs = lm(WmaxOBS ~ avgsst, data=datC1.obs)   
resC1.obs = residuals(lmC1.obs)
ciC1.obs = confint(lmC1.obs)   # 4.54, 9.02
## 6.8 +/- 1.09  
## t-val = 6.25
## R^2 = 0.6196

lmC1.mod = lm(WmaxMOD ~ avgsst, data=datC1.mod)
resC1.mod = residuals(lmC1.mod)
ciC1.mod = confint(lmC1.mod)   # -3.79, 4.71
## 0.463 +/- 2.09  
## t-val = 0.221
## R^2 = 0.00144

r2.C1.obs = summary(lmC1.obs)$r.squared
r2.C1.mod = summary(lmC1.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.63, 8.30
## 6.0 +/- 1.13  
## t-val = 5.27
## R^2 = 0.537

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.64, 8.34
## 6.0 +/- 1.13  
## t-val = 5.64
## R^2 = 0.534

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.

pdatMRI = get.pdat(hspdfM.obs, hspdfM.mod)
pdatMRI$lab = rep("(a) MRI", dim(pdatMRI)[1])

pdatCAM = get.pdat(hspdfC.obs, hspdfC.mod)
pdatCAM$lab = rep("(b) CAM5 0.25", dim(pdatCAM)[1])

pdatGFDL = get.pdat(hspdfG.obs, hspdfG.mod)
pdatGFDL$lab = rep("(c) GFDL", dim(pdatGFDL)[1])

pdatFSU = get.pdat(hspdfF.obs, hspdfF.mod)
pdatFSU$lab = rep("(d) FSU", dim(pdatFSU)[1])

pdatCAM1 = get.pdat(hspdfC1.obs, hspdfC1.mod)
pdatCAM1$lab = rep("(e) CAM5 1", dim(pdatCAM1)[1])

pdat = rbind(pdatCAM, pdatMRI, pdatGFDL, pdatFSU, pdatCAM1)
xlabel = expression(paste("Sea surface temperature [", degree, "C]"))
ylabel = expression(paste("Maximum Intensity [m ", s^-1, "]"))

p25 = ggplot(pdat, aes(x=avgsst, y=Wmax, color=type)) + geom_point()  
p25 = p25 + facet_grid(~lab)
p25 = p25 + theme_bw() + xlab(xlabel) + ylab(ylabel)
p25 = p25 + geom_smooth(method="lm", 
                        aes(weight=count, color=type, fill=type))
p25 = p25 + scale_color_manual( values = c("blue","red"))
p25 = p25 + scale_fill_manual( values = c("blue","red"))
p25 = p25 + theme(axis.text=element_text(size=14), 
                    axis.title=element_text(size=16,face="bold"))
p25 = p25 + theme(legend.text=element_text(size=14),
                    legend.title=element_text(size=14, face="bold"))
p25 = p25 + theme(strip.text.x = element_text(size = 12, face="bold"))
p25 = p25 + theme(legend.title=element_blank())
p25

Subsetting on same set of hexagons for an MRI - OBS comparison

#Plot all hexagons -- the gulf of mexico
p1 = spplot(hspdfM.obs2, "WmaxOBS", col="white", col.regions=crWind,
            sp.layout=list(l1), xlim=bb[1, ],
            ylim=c(bb[2, 1]*-1.1, bb[2, 2]),
            at=rngWind, 
            colorkey=list(space="bottom", labels=paste(rngWind)),
            sub=expression(paste("Observed Maximum Intensity (m "~s^-1~")")))

p2 = spplot(hspdfM.mod, "WmaxMOD", col="white", col.regions=crWind,
            sp.layout=list(l1), xlim=bb[1, ],
            ylim=c(bb[2, 1]*-1.1, bb[2, 2]),
            at=rngWind, 
            colorkey=list(space="bottom", labels=paste(rngWind)),
            sub=expression(paste("MRI Maximum Intensity (m "~s^-1~")")))

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)

Comparing two slopes (Observed sensitivity with GoM vs. without GoM regions)

C = hspdfM.obs@data
C = subset(C, avgsst >= 25)
C$type = rep("All", dim(C)[1])
D = hspdfM.obs2@data
D = subset(D, avgsst >= 25)
D$type = rep("No Gulf of Mexico", dim(D)[1])

pdat = rbind(C, D)
xlabel = expression(paste("Sea surface temperature [", degree, "C]"))
ylabel = expression(paste("Maximum Intensity [m ", s^-1, "]"))

p25 = ggplot(pdat, aes(x=avgsst, y=WmaxOBS, color=type)) + geom_point()  
p25 = p25 + theme_bw() + xlab(xlabel) + ylab(ylabel)
p25 = p25 + geom_smooth(method="lm", aes(color=type, fill=type))
p25 = p25 + scale_color_manual( values = c("blue","red"))
p25 = p25 + scale_fill_manual( values = c("blue","red"))
p25 = p25 + theme(axis.text=element_text(size=14), 
                    axis.title=element_text(size=16,face="bold"))
p25 = p25 + theme(legend.text=element_text(size=14),
                    legend.title=element_text(size=14, face="bold"))
p25 = p25 + theme(legend.title=element_blank())
p25

---

---

---

---

---

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")

---

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

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

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(CAM1.sdf, type="hexagonal", n=150, bb=bb)
  hpg = HexPoints2SpatialPolygons(hpts)
  allSensCAM1[i,] = get.sensitivity3(WindC1.sdf, 
                                     CAM1.sdf, hpg, 
                                     sstC1.sdf, 1979, 2007)
}

colnames(allSensCAM1) = c("OBS", "CAM1")
range(allSensCAM1$OBS)  ## 3.18, 10.39
range(allSensCAM1$CAM1)  ## -2.38, 0.751
median(allSensCAM1$OBS) ##  7.31
median(allSensCAM1$CAM1) ##  -0.986

save(allSensCAM1, file="allSensCAM1.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/allSensCAM1.RData")
load("~/Dropbox/SpatialExtremes/allSensMRI.RData")
load("~/Dropbox/SpatialExtremes/allSensGFDL.RData")
load("~/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) CAM5 0.25", 200)
SC.df = data.frame(SC.df)

SC1 = c(allSensCAM1$OBS, allSensCAM1$CAM1)
Label = c(rep("Best-track", 100), rep("GCM", 100))
SC1.df = data.frame(Label = Label, Sens = SC1)
SC1.df$type = rep("(e) CAM5 1", 200)
SC1.df = data.frame(SC1.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, SC1.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

**** Comparison with the sensitivity of PI to SST ****

load("~/Dropbox/spatialextremes/Sarah/PI_12z.RData")
load("~/Dropbox/spatialextremes/Sarah/PIfsu_12z.RData")
load("~/Dropbox/spatialextremes/Sarah/PIcam1.RData")
load("~/Dropbox/spatialextremes/Sarah/PIcamQ.RData")
load("~/Dropbox/spatialextremes/Sarah/PIgfdl.RData")
load("~/Dropbox/spatialextremes/Sarah/PImri.RData")

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

PI.comp.pdf = over(x=hpg[rownames(hexDataF.obs)], y=PI.sdf, fn=max)
PI.pdf = over(x=hpg[rownames(datF.obs)], y=PI.sdf, fn=max)
months = rep(c("M08", "M09", "M10"), 27)
yrs = rep(1982:2008, each=3)
cn = paste("Y", yrs, months, sep="")
colnames(PI.pdf) = cn
colnames(PI.comp.pdf) = cn
datF.obs$PI = do.call(pmax, PI.pdf)
hspdfF.obs$PIobs = do.call(pmax, PI.comp.pdf)

##################  FSU-COAPS  #######################

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

PIfsu.comp.pdf = over(x=hpg[rownames(hexDataF.obs)], y=PIfsu.sdf, fn=max)
PIfsu.pdf = over(x=hpg[rownames(datF.mod)], y=PIfsu.sdf, fn=max)
colnames(PIfsu.pdf) = cn
colnames(PIfsu.comp.pdf) = cn
datF.mod$PI = do.call(pmax, PIfsu.pdf)
hspdfF.obs$PIfsu = do.call(pmax, PIfsu.comp.pdf)

##################  CAM 1deg  #######################

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

PIcam1.comp.pdf = over(x=hpg[rownames(hexDataF.obs)], y=PIcam1.sdf, fn=max)
PIcam1.pdf = over(x=hpg[rownames(datC1.mod)], y=PIcam1.sdf, fn=max)
colnames(PIcam1.pdf) = cn
colnames(PIcam1.comp.pdf) = cn
datC1.mod$PI = do.call(pmax, PIcam1.pdf)
hspdfF.obs$PIcam1 = do.call(pmax, PIcam1.comp.pdf)

################  CAM quarter  ####################

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

PIcamQ.comp.pdf = over(x=hpg[rownames(hexDataF.obs)], y=PIcamQ.sdf, fn=max)
PIcamQ.pdf = over(x=hpg[rownames(datC.mod)], y=PIcamQ.sdf, fn=max)
colnames(PIcamQ.pdf) = cn
colnames(PIcamQ.pdf) = cn
datC.mod$PI = do.call(pmax, PIcamQ.pdf)
hspdfF.obs$PIcamQ = do.call(pmax, PIcamQ.comp.pdf)


##################  GFDL  #######################

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

PIgfdl.comp.pdf = over(x=hpg[rownames(hexDataF.obs)], y=PIgfdl.sdf, fn=max)
PIgfdl.pdf = over(x=hpg[rownames(datG.mod)], y=PIgfdl.sdf, fn=max)
colnames(PIgfdl.pdf) = cn
colnames(PIgfdl.comp.pdf) = cn
datG.mod$PI = do.call(pmax, PIgfdl.pdf)
hspdfF.obs$PIgfdl = do.call(pmax, PIgfdl.comp.pdf)


##################  MRI  #######################

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

PImri.comp.pdf = over(x=hpg[rownames(hexDataF.obs)], y=PImri.sdf, fn=max)
PImri.pdf = over(x=hpg[rownames(datM.mod)], y=PImri.sdf, fn=max)
colnames(PImri.pdf) = cn
colnames(PImri.comp.pdf) = cn
datM.mod$PI = do.call(pmax, PImri.pdf)
hspdfF.obs$PImri = do.call(pmax, PImri.comp.pdf)

save as 10 by 12

p1 = spplot(hspdfF.obs, "PIobs", col="white", col.regions=crWind,
            sp.layout=list(l1), xlim=bb[1, ],
            ylim=c(bb[2, 1]*-1.1, bb[2, 2]),
            at=rngWind, 
            colorkey=list(space="bottom", labels=paste(rngWind)),
            sub=expression(paste("Potential Intensity (m "~s^-1~")")))

p2 = spplot(hspdfF.obs, "PImri", col="white", col.regions=crWind,
            sp.layout=list(l1), xlim=bb[1, ],
            ylim=c(bb[2, 1]*-1.1, bb[2, 2]),
            at=rngWind, 
            colorkey=list(space="bottom", labels=paste(rngWind)),
            sub=expression(paste("Potential Intensity (m "~s^-1~")")))

p3 = spplot(hspdfF.obs, "PIcamQ", col="white", col.regions=crWind,
            sp.layout=list(l1), xlim=bb[1, ],
            ylim=c(bb[2, 1]*-1.1, bb[2, 2]),
            at=rngWind, 
            colorkey=list(space="bottom", labels=paste(rngWind)),
            sub=expression(paste("Potential Intensity (m "~s^-1~")")))

p4 = spplot(hspdfF.obs, "PIgfdl", col="white", col.regions=crWind,
            sp.layout=list(l1), xlim=bb[1, ],
            ylim=c(bb[2, 1]*-1.1, bb[2, 2]),
            at=rngWind, 
            colorkey=list(space="bottom", labels=paste(rngWind)),
            sub=expression(paste("Potential Intensity (m "~s^-1~")")))

p5 = spplot(hspdfF.obs, "PIfsu", col="white", col.regions=crWind,
            sp.layout=list(l1), xlim=bb[1, ],
            ylim=c(bb[2, 1]*-1.1, bb[2, 2]),
            at=rngWind, 
            colorkey=list(space="bottom", labels=paste(rngWind)),
            sub=expression(paste("Potential Intensity (m "~s^-1~")")))

p6 = spplot(hspdfF.obs, "PIcam1", col="white", col.regions=crWind,
            sp.layout=list(l1), xlim=bb[1, ],
            ylim=c(bb[2, 1]*-1.1, bb[2, 2]),
            at=rngWind, 
            colorkey=list(space="bottom", labels=paste(rngWind)),
            sub=expression(paste("Potential Intensity (m "~s^-1~")")))

p1 = update(p1, main = textGrob("(a) MERRA", x = unit(0.05, "npc")),
            par.settings = list(fontsize = list(text = 14)))
p2 = update(p2, main = textGrob("(b) MRI", x = unit(0.05, "npc")), 
            par.settings = list(fontsize = list(text = 14)))
p3 = update(p3, main = textGrob("(c) CAM5 0.25", x = unit(0.05, "npc")), 
            par.settings = list(fontsize = list(text = 14)))
p4 = update(p4, main = textGrob("(d) GFDL", x = unit(0.05, "npc")), 
            par.settings = list(fontsize = list(text = 14)))
p5 = update(p5, main = textGrob("(e) FSU", x = unit(0.05, "npc")), 
            par.settings = list(fontsize = list(text = 14)))
p6 = update(p6, main = textGrob("(f) CAM5 1", x = unit(0.05, "npc")), 
            par.settings = list(fontsize = list(text = 14)))

plot(p1, split=c(1, 1, 2, 3), more=TRUE)
plot(p2, split=c(2, 1, 2, 3), more=TRUE)
plot(p3, split=c(1, 2, 2, 3), more=TRUE)
plot(p4, split=c(2, 2, 2, 3), more=TRUE)
plot(p5, split=c(1, 3, 2, 3), more=TRUE)
plot(p6, split=c(2, 3, 2, 3), more=FALSE)

hspdfF.obs$Fdiff = hspdfF.obs$PIobs - hspdfF.obs$PIfsu
hspdfF.obs$Gdiff = hspdfF.obs$PIobs - hspdfF.obs$PIgfdl
hspdfF.obs$Mdiff = hspdfF.obs$PIobs - hspdfF.obs$PImri
hspdfF.obs$CQdiff = hspdfF.obs$PIobs - hspdfF.obs$PIcamQ
hspdfF.obs$C1diff = hspdfF.obs$PIobs - hspdfF.obs$PIcam1

Difference Maps

crPD = rev(brewer.pal(10, "RdBu"))
rngPD = seq(-10, 10, 2)

p1 = spplot(hspdfF.obs, "Mdiff", col="white", col.regions=crPD,
            sp.layout=list(l1), xlim=bb[1, ],
            ylim=c(bb[2, 1]*-1.1, bb[2, 2]),
            at=rngPD, 
            colorkey=list(space="bottom", labels=paste(rngPD)),
            sub=expression(paste("MERRA - Model PI (m "~s^-1~")")))

p2 = spplot(hspdfF.obs, "CQdiff", col="white", col.regions=crPD,
            sp.layout=list(l1), xlim=bb[1, ],
            ylim=c(bb[2, 1]*-1.1, bb[2, 2]),
            at=rngPD, 
            colorkey=list(space="bottom", labels=paste(rngPD)),
            sub=expression(paste("MERRA - Model PI (m "~s^-1~")")))

p3 = spplot(hspdfF.obs, "Gdiff", col="white", col.regions=crPD,
            sp.layout=list(l1), xlim=bb[1, ],
            ylim=c(bb[2, 1]*-1.1, bb[2, 2]),
            at=rngPD, 
            colorkey=list(space="bottom", labels=paste(rngPD)),
            sub=expression(paste("MERRA - Model PI (m "~s^-1~")")))

p4 = spplot(hspdfF.obs, "Fdiff", col="white", col.regions=crPD,
            sp.layout=list(l1), xlim=bb[1, ],
            ylim=c(bb[2, 1]*-1.1, bb[2, 2]),
            at=rngPD, 
            colorkey=list(space="bottom", labels=paste(rngPD)),
            sub=expression(paste("MERRA - Model PI (m "~s^-1~")")))

p5 = spplot(hspdfF.obs, "C1diff", col="white", col.regions=crPD,
            sp.layout=list(l1), xlim=bb[1, ],
            ylim=c(bb[2, 1]*-1.1, bb[2, 2]),
            at=rngPD, 
            colorkey=list(space="bottom", labels=paste(rngPD)),
            sub=expression(paste("MERRA - Model PI (m "~s^-1~")")))

p1 = update(p1, main = textGrob("(a) MRI", x = unit(0.05, "npc")),
            par.settings = list(fontsize = list(text = 14)))
p2 = update(p2, main = textGrob("(b) CAM5 0.25", x = unit(0.05, "npc")), 
            par.settings = list(fontsize = list(text = 14)))
p3 = update(p3, main = textGrob("(c) GFDL", x = unit(0.05, "npc")), 
            par.settings = list(fontsize = list(text = 14)))
p4 = update(p4, main = textGrob("(d) FSU", x = unit(0.05, "npc")), 
            par.settings = list(fontsize = list(text = 14)))
p5 = update(p5, main = textGrob("(e) CAM5 1", x = unit(0.05, "npc")), 
            par.settings = list(fontsize = list(text = 14)))

plot(p1, split=c(1, 1, 2, 3), more=TRUE)
plot(p2, split=c(2, 1, 2, 3), more=TRUE)
plot(p3, split=c(1, 2, 2, 3), more=TRUE)
plot(p4, split=c(2, 2, 2, 3), more=TRUE)
plot(p5, split=c(1, 3, 2, 3), more=FALSE)

Need wind, hexid, count, SST

o1 = cbind("SST" = datF.obs$avgsst, "Wind" = datF.obs$WmaxOBS)
o2 = cbind("SST"= datF.obs$avgsst, "Wind" = datF.obs$PI)
o.all = data.frame(rbind(o1, o2))
o.all$Type = factor(c(rep("Wmax", dim(o1)[1]), rep("PI", dim(o2)[1])))
o.all$Data = rep("(a) MERRA", dim(o.all)[1])

m1 = cbind("SST" = datM.mod$avgsst, "Wind" = datM.mod$WmaxMOD)
m2 = cbind("SST"= datM.mod$avgsst, "Wind" = datM.mod$PI)
m.all = data.frame(rbind(m1, m2))
m.all$Type = factor(c(rep("Wmax", dim(m1)[1]), rep("PI", dim(m2)[1])))
m.all$Data = rep("(b) MRI", dim(m.all)[1])

cq1 = cbind("SST" = datC.mod$avgsst, "Wind" = datC.mod$WmaxMOD)
cq2 = cbind("SST"= datC.mod$avgsst, "Wind" = datC.mod$PI)
cq.all = data.frame(rbind(cq1, cq2))
cq.all$Type = factor(c(rep("Wmax", dim(cq1)[1]), rep("PI", dim(cq2)[1])))
cq.all$Data = rep("(c) CAM5 0.25", dim(cq.all)[1])

g1 = cbind("SST" = datG.mod$avgsst, "Wind" = datG.mod$WmaxMOD)
g2 = cbind("SST"= datG.mod$avgsst, "Wind" = datG.mod$PI)
g.all = data.frame(rbind(g1, g2))
g.all$Type = factor(c(rep("Wmax", dim(g1)[1]), rep("PI", dim(g2)[1])))
g.all$Data = rep("(d) GFDL", dim(g.all)[1])

f1 = cbind("SST" = datF.mod$avgsst, "Wind" = datF.mod$WmaxMOD)
f2 = cbind("SST"= datF.mod$avgsst, "Wind" = datF.mod$PI)
f.all = data.frame(rbind(f1, f2))
f.all$Type = factor(c(rep("Wmax", dim(f1)[1]), rep("PI", dim(f2)[1])))
f.all$Data = rep("(e) FSU", dim(f.all)[1])

c11 = cbind("SST" = datC1.mod$avgsst, "Wind" = datC1.mod$WmaxMOD)
c12 = cbind("SST"= datC1.mod$avgsst, "Wind" = datC1.mod$PI)
c1.all = data.frame(rbind(c11, c12))
c1.all$Type = factor(c(rep("Wmax", dim(c11)[1]), rep("PI", dim(c12)[1])))
c1.all$Data = rep("(f) CAM5 1", dim(c1.all)[1])

all = rbind(o.all, m.all, cq.all, g.all, f.all, c1.all)

xlabel = expression(paste("Sea surface temperature [", degree, "C]"))
ylabel = expression(paste("Intensity [m ", s^-1, "]"))

pS = ggplot(all, aes(x=SST, y=Wind, color=Type)) 
pS = pS + geom_point() + facet_grid(~Data)
pS = pS + theme_bw() + xlab(xlabel) + ylab(ylabel)
pS = pS + geom_smooth(method="lm", aes(color=Type, fill=Type))
pS = pS + theme(axis.text=element_text(size=14), 
                    axis.title=element_text(size=16,face="bold"))
pS = pS + theme(legend.text=element_text(size=14),
                    legend.title=element_text(size=14, face="bold"))
pS = pS + theme(strip.text=element_text(size=14),
                    strip.text=element_text(size=14, face="bold"))
pS = pS + scale_color_manual(values = c("black", "purple"))
pS = pS + scale_fill_manual(values = c("black", "purple"))
pS = pS + theme(legend.title=element_blank())
pS

lmPI.obs = lm(PI ~ avgsst, data=datF.obs)   
resPI.obs = residuals(lmPI.obs)  
ciPI.obs = confint(lmPI.obs)   # 2.94, 4.35
# Slope = 3.64
# S.E. = 0.341
# T-val = 10.7
# R^2 = 0.8268

lmPI.mri = lm(PI ~ avgsst, data=datM.mod)   
resPI.mri = residuals(lmPI.mri)  
ciPI.mri = confint(lmPI.mri)   # 2.72, 4.15
# Slope = 3.43
# S.E. = 0.3473
# T-val = 9.89
# R^2 = 0.810

lmPI.camQ = lm(PI ~ avgsst, data=datC.mod)   
resPI.camQ = residuals(lmPI.camQ)  
ciPI.camQ = confint(lmPI.camQ)   # 3.43, 4.94
# Slope = 4.18
# S.E. = 0.368
# T-val = 11.38
# R^2 = 0.817

lmPI.gfdl = lm(PI ~ avgsst, data=datG.mod)   
resPI.gfdl = residuals(lmPI.gfdl)  
ciPI.gfdl = confint(lmPI.gfdl)   # 3.34, 4.91
# Slope = 4.12
# S.E. = 0.383
# T-val = 10.78
# R^2 = 0.811

lmPI.fsu = lm(PI ~ avgsst, data=datF.mod)   
resPI.fsu = residuals(lmPI.fsu)  
ciPI.fsu = confint(lmPI.fsu)   # 2.31, 4.28
# Slope = 3.30
# S.E. = 0.477
# T-val = 6.92
# R^2 = 0.675

lmPI.cam1 = lm(PI ~ avgsst, data=datC1.mod)   
resPI.cam1 = residuals(lmPI.cam1)  
ciPI.cam1 = confint(lmPI.cam1)   # 2.31, 4.06
# Slope = 3.19
# S.E. = 0.405
# T-val = 7.87
# R^2 = 0.827

Make wind shear comparison with FSU-COAPS model (OBS, FSU-COAPS, DIFF)

load("~/Dropbox/SpatialExtremes/Sarah/Space-time/vwsData.RData")
load("~/Dropbox/SpatialExtremes/Sarah/vwsFSU.RData")

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

vws.pdf = over(x=hpg[rownames(hexDataF.obs)], y=vwsData.sdf, fn=mean)
vwsYearMonth = strsplit(substring(colnames(vws.pdf), 5), "M")
vwsYear = as.numeric(sapply(vwsYearMonth, function(x) x[1]))

vws.yr = sapply(unique(vwsYear), 
    function(x) as.vector(rowMeans(vws.pdf[,which(x==vwsYear)])))

dimnames(vws.yr) = list(id=rownames(vws.pdf), 
                   Year=paste("Y",unique(vwsYear), sep="")) 

vws.yr = data.frame(vws.yr)
vws.yr = vws.yr[,4:30]

be = rep("vwsY", 81)
ye = rep(1982:2008, each=3)
m = rep(c("M08", "M09", "M10"), 27)
colnames(vwsFSU) = c("lon", "lat", paste(be, ye, m, sep=""))

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

vws.fsu.pdf = over(x=hpg[rownames(hexDataF.obs)], y=vwsFSU.sdf, fn=mean)
vws.fsu.pdf2 = over(x=hpg[rownames(hexDataF.mod)], y=vwsFSU.sdf, fn=mean)
vwsYearMonth = strsplit(substring(colnames(vws.fsu.pdf), 5), "M")
vwsYear = as.numeric(sapply(vwsYearMonth, function(x) x[1]))

vws.fsu.yr = sapply(unique(vwsYear), 
    function(x) as.vector(rowMeans(vws.fsu.pdf[,which(x==vwsYear)])))
vws.fsu.yr2 = sapply(unique(vwsYear), 
    function(x) as.vector(rowMeans(vws.fsu.pdf2[,which(x==vwsYear)])))

dimnames(vws.fsu.yr) = list(id=rownames(vws.fsu.pdf), 
                   Year=paste("Y",unique(vwsYear), sep=""))
dimnames(vws.fsu.yr2) = list(id=rownames(vws.fsu.pdf2), 
                   Year=paste("Y",unique(vwsYear), sep="")) 

vws.fsu.yr = data.frame(vws.fsu.yr)
vws.fsu.yr2 = data.frame(vws.fsu.yr2)

hspdfF.obs$vwsOBS = rowMeans(vws.yr)
hspdfF.obs$vwsFSU = rowMeans(vws.fsu.yr)
hspdfF.obs$vwsDIFF = hspdfF.obs$vwsOBS - hspdfF.obs$vwsFSU

hspdfF.mod$vwsFSU = rowMeans(vws.fsu.yr2)

dat.1 = as.data.frame(hspdfF.obs)
dat.2 = as.data.frame(hspdfF.mod)


lm.1 = lm(WmaxOBS ~ vwsOBS, data = dat.1)  
### slope = -1.76, ste = 0.5961, t-val = -2.96, R^2 = 0.2381
lm.1.1 = lm(WmaxOBS ~ avgsst + vwsOBS, data=dat.1)
### slopeVWS = 0.591, ste = 0.6688, t-val = 0.884

lm.2 = lm(WmaxMOD ~ vwsFSU, data = dat.2)
### slope = -0.248, ste = 0.2207, t-val = -1.12
lm.2.1 = lm(WmaxMOD ~ avgsst + vwsFSU, data = dat.2)
###

crPD = rev(brewer.pal(8, "RdBu"))
rngPD = seq(-8, 8, 2)

spplot(hspdfF.obs, "vwsDIFF", col="white", col.regions=crPD,
            sp.layout=list(l1), xlim=bb[1, ],
            ylim=c(bb[2, 1]*-1.1, bb[2, 2]),
            at=rngPD, 
            colorkey=list(space="bottom", labels=paste(rngPD)),
            sub=expression(paste("MERRA - FSU vertical shear (m "~s^-1~")")))

Add pmin to the mix and calculate the sensitivity of pmin to SST

load("~/Dropbox/GCMTrackData/CAM.pmin.RData")
load("~/Dropbox/GCMTrackData/CAM1.pmin.RData")
load("~/Dropbox/GCMTrackData/mri.pmin.RData")
load("~/Dropbox/GCMTrackData/fsuPmin_r1.RData")
load("~/Dropbox/GCMTrackData/gfdlPmin_r1.RData")

CAM.pmin = data.frame(cbind("lon" = CAM.pmin$lon, 
                 "lat" = CAM.pmin$lat, 
                 "pmin" = CAM.pmin$pmin))
CAM.pmin = subset(CAM.pmin, lon >= min(CAM.df$lon) & lon <= max(CAM.df$lon))

CAM1.pmin = data.frame(cbind("lon" = CAM1.pmin$lon, 
                 "lat" = CAM1.pmin$lat, 
                 "pmin" = CAM1.pmin$pmin))
CAM1.pmin = subset(CAM1.pmin, lon >= min(CAM1.df$lon) & lon <= max(CAM1.df$lon))

MRI.pmin = data.frame(cbind("lon" = mri.pmin$lon, 
                 "lat" = mri.pmin$lat, 
                 "pmin" = mri.pmin$pmin))
MRI.pmin = subset(MRI.pmin, lon >= min(MRI.df$lon) & lon <= max(MRI.df$lon))

FSU.pmin = data.frame(cbind("lon" = fsuPmin_r1$lon, 
                 "lat" = fsuPmin_r1$lat, 
                 "pmin" = fsuPmin_r1$pmin))
FSU.pmin = subset(FSU.pmin, lon >= min(FSU.df$lon) & lon <= max(FSU.df$lon))

GFDL.pmin = data.frame(cbind("lon" = gfdlPmin_r1$lon, 
                 "lat" = gfdlPmin_r1$lat, 
                 "pmin" = gfdlPmin_r1$pmin))
GFDL.pmin = subset(GFDL.pmin, lon >= min(GFDL.df$lon) & lon <= max(GFDL.df$lon))

OBS.pmin = data.frame(cbind("lon" = Wind.df$lon, 
                            "lat" = Wind.df$lat,
                            "pmin" = Wind.df$pmin))

Overlay pmin onto hexagons

coordinates(CAM.pmin) = c("lon", "lat")
proj4string(CAM.pmin) = CRS(ll)
CAM.pmin.sdf = spTransform(CAM.pmin, CRS(lcc))
CAM.pmin.pdf = over(x=hpg[rownames(datC.mod)], 
                    y=CAM.pmin.sdf, fn=min)
datC.mod$pmin = CAM.pmin.pdf$pmin


coordinates(CAM1.pmin) = c("lon", "lat")
proj4string(CAM1.pmin) = CRS(ll)
CAM1.pmin.sdf = spTransform(CAM1.pmin, CRS(lcc))
CAM1.pmin.pdf = over(x=hpg[rownames(datC1.mod)], 
                    y=CAM1.pmin.sdf, fn=min)
datC1.mod$pmin = CAM1.pmin.pdf$pmin


coordinates(FSU.pmin) = c("lon", "lat")
proj4string(FSU.pmin) = CRS(ll)
FSU.pmin.sdf = spTransform(FSU.pmin, CRS(lcc))
FSU.pmin.pdf = over(x=hpg[rownames(datF.mod)], 
                    y=FSU.pmin.sdf, fn=min)
datF.mod$pmin = FSU.pmin.pdf$pmin


coordinates(GFDL.pmin) = c("lon", "lat")
proj4string(GFDL.pmin) = CRS(ll)
GFDL.pmin.sdf = spTransform(GFDL.pmin, CRS(lcc))
GFDL.pmin.pdf = over(x=hpg[rownames(datG.mod)], 
                    y=GFDL.pmin.sdf, fn=min)
datG.mod$pmin = GFDL.pmin.pdf$pmin


coordinates(MRI.pmin) = c("lon", "lat")
proj4string(MRI.pmin) = CRS(ll)
MRI.pmin.sdf = spTransform(MRI.pmin, CRS(lcc))
MRI.pmin.pdf = over(x=hpg[rownames(datM.mod)], 
                    y=MRI.pmin.sdf, fn=min)
datM.mod$pmin = MRI.pmin.pdf$pmin


coordinates(OBS.pmin) = c("lon", "lat")
proj4string(OBS.pmin) = CRS(ll)
OBS.pmin.sdf = spTransform(OBS.pmin, CRS(lcc))
OBS.pmin.pdf = over(x=hpg[rownames(datM.obs)], 
                    y=OBS.pmin.sdf, fn=min)
datM.obs$pmin = OBS.pmin.pdf$pmin

Calculate sensitivity of pmin to SST

lmPmin.obs = lm(pmin ~ avgsst, data=datM.obs)   
ciPmin.obs = confint(lmPmin.obs)   # -18.22, -7.909
# Slope = -13.1
# S.E. = 2.491
# T-val = -5.24 ***
# R^2 = 0.5445 

lmPmin.mri = lm(pmin ~ avgsst, data=datM.mod)   
ciPmin.mri = confint(lmPmin.mri)   # -4.55, 8.37
# Slope = 1.91
# S.E. = 3.123
# T-val = 0.611
# R^2 = 0.016

lmPmin.camQ = lm(pmin ~ avgsst, data=datC.mod)   
ciPmin.camQ = confint(lmPmin.camQ)   # -7.561, 5.564
# Slope = -0.9982
# S.E. = 3.209
# T-val = -0.311
# R^2 = 0.00333

lmPmin.gfdl = lm(pmin ~ avgsst, data=datG.mod)   
ciPmin.gfdl = confint(lmPmin.gfdl)   # -3.72, 6.16
# Slope = 1.22
# S.E. = 2.41
# T-val = 0.508
# R^2 = 0.00948

lmPmin.fsu = lm(pmin ~ avgsst, data=datF.mod)   
ciPmin.fsu = confint(lmPmin.fsu)   # -0.75, 12
# Slope = 5.59
# S.E. = 3.07
# T-val = 1.82
# R^2 = 0.675

lmPmin.cam1 = lm(pmin ~ avgsst, data=datC1.mod)   
ciPmin.cam1 = confint(lmPmin.cam1)   # 1.62, 9.43
# Slope = 5.53
# S.E. = 1.81
# T-val = 3.05
# R^2 = 0

————————

---