Previous research quantified the sensitivity of per region maximum and limiting intensity to SST for both observed tropical cyclones (TCs) as well as for TCs generated by two low-resolution global climate models (GCMs). While a 1 K increase in sea surface temperature (SST) is associated with a 5–12 m s\(^{-1}\) increase in per region maximum TC intensity over the North Atlantic, the GFDL-HiRAM and FSU/COAPS models produced TCs with sensitivities of less than 2 m s\(^{-1}\) K\(^{-1}\). The research presented here builds on this previous research in two ways: 1) We link the statistically defined sensitivity of observed limiting intensity to SST to the theoretical potential intensity theory, and 2) we examine the ability of a low-resolution GCM to reproduce the observed statistical relationship between potential intensity and SST.

To address these issues, we use NASA’s Modern Era Reanalysis to first estimate the observed sensitivity of per region potential intensity to SST for the 1982–2008 time period. Using a spatial framework, we show that the sensitivity of potential intensity to SST is not statistically different from the sensitivity of observed maximum or limiting intensity to SST. This result therefore links our statistical analysis to the more physically based theory of hurricanes and supports the idea that the maximum wind speeds of the most intense TCs (i.e. those that reach at least 75 % of their potential intensity) should increase as SSTs increase. We similarly calculate the per region potential intensity estimates from the FSU/COAPS GCM. Although the FSU/COAPS model does not capture the observed sensitivity of TC maximum intensity to SST, the model is able to reproduce the observed sensitivity of potential intensity to SST for the 1982–2008 period. This supports the idea that model resolution must be increased in order to accurately represent TCs in climate models. The model generates suitable atmospheric conditions for the development of strong TCs in the correct regions, however, likely as a result of insufficient resolution, strong TCs do not develop.

setwd("/Users/sarahstrazzo/Dropbox/spatialextremes")

require(ggplot2); require(sp); 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)
require(spacetime); 
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")
source("~/Dropbox/spatialextremes/Sarah/Space-Time/arrangeMERRA.R")
source("~/Dropbox/spatialextremes/Sarah/LoI2functions.R")
source("~/Dropbox/spatialextremes/Sarah/pcmin_revised.R")
source("~/Dropbox/spatialextremes/Sarah/arrangeFSU.R")

# function from TszKin Julian on R-bloggers, 2 Sep 2011
rep.row<-function(x,n){      
   matrix(rep(x,each=n),nrow=n)
}

Data

Load best track and GCM track data, subset on tropical cyclones, convert wind speeds to meters per second, and estimate the rotational wind speed by subtracting 60% of the forward motion (maguv)

#load("C:/Users/Sarah/Dropbox/SpatialExtremes/best.use.RData")
load("/Users/sarahstrazzo/Dropbox/spatialextremes/best.use.RData")
begin = 1982; end = 2008
years = begin:end
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
nTC = length(unique(Wind.df$Sid))  # 314

## FSU/COAPS track data
#load("C:/Users/Sarah/Dropbox/SpatialExtremes/fsu_r1i1p1_ATL.RData")
load("/Users/sarahstrazzo/Dropbox/spatialextremes/fsu_r1i1p1_ATL.RData")
fsu.df = fsuTracks.df[fsuTracks.df$Yr >= begin & fsuTracks.df$Yr <= end, ]
fsu.df$WmaxS = fsu.df$WmaxS - fsu.df$maguv * .6
fsu.df = subset(fsu.df, WmaxS > 0)
nTCf = length(unique(fsu.df$Sid))  # 382

Figure 1 highlights the differences in the wind speed distributions for observed (a) and simulated (b) TCs. The distribution of simulated TC wind speeds is truncated at approximately 50 m s\(^{-1}\), while the tail of the observed wind speed distribution extends well beyond 70 m s\(^{-1}\). Table 1 reports basic TC statistics for observed and simulated TCs over the 1982–2008 time period and again illustrates the lack of strong TCs in the simulated data set.

Figure 1: Comparison of wind speed probablity distributions for best-track and GCM-generated TCs (Rotational intensity shown in m/s)

Insert Table 1 Here

Generate a spatial data frame from TC data

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

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

Methods

Define a set of equal area (~407 636 km\(^{2}\)) hexagons that completely cover the set of TC tracks

set.seed(42)
bb = matrix(0, nrow=2, ncol=2)
bb[2,1] = bbox(Wind.sdf)[2,1]*-1.01     
bb[1,1] = bbox(Wind.sdf)[1,1]*1.85
bb[1,2] = bbox(Wind.sdf)[1,2]*1.25
bb[2,2] = bbox(Wind.sdf)[2,2]*1.15
hpt = spsample(Wind.sdf, type="hexagonal", n=200, bb=bb)
hpg = HexPoints2SpatialPolygons(hpt)
#hpg@polygons[[1]]@area * 10^-6   ## 731,935.3 km^2

Overlay TC data onto hexagons

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

############## FSU COAPS ################

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

Create a data frame with the per hexagon maximum TC wind speed (OBS)

windSplit = split(Wind.sdf@data, Wind.sdf$hexid)
maxWinds = matrix(0, nrow=length(windSplit),
                     ncol=dim(Wind.sdf@data)[2])
maxWinds = data.frame(maxWinds)
colnames(maxWinds) = colnames(Wind.sdf@data)
r = 1

for(i in 1:length(windSplit)){
  dat = windSplit[[i]]
  m = max(dat$WmaxS)
  m2 = subset(dat, WmaxS == m)
  maxWinds[r, ] = m2
  r = r+1
}

# dim(maxWinds) = 61, 37

Create a data frame with the per hexagon maximum TC wind speed (GCM)

fsuSplit = split(fsu.sdf@data, fsu.sdf$hexid)
maxFSU = matrix(0, nrow=length(fsuSplit),
                     ncol=dim(fsu.sdf@data)[2])
maxFSU = data.frame(maxFSU)
colnames(maxFSU) = colnames(fsu.sdf@data)
r = 1

for(i in 1:length(fsuSplit)){
  dat = fsuSplit[[i]]
  m = max(dat$WmaxS)
  m2 = subset(dat, WmaxS == m)
  maxFSU[r, ] = m2
  r = r+1
}

# dim(maxFSU) = 66, 17

Calculate per hexagon TC counts

stats = ddply(Wind.sdf@data, .(hexid, Sid), summarize,
                 nhr = length(Sid),
                 WmaxS = max(WmaxS))
stats2 = ddply(stats, .(hexid), summarize, 
                 count = length(Sid))
maxWinds$count = stats2$count    # dim = 61 by 38

############## FSU COAPS ################

Fstats = ddply(fsu.sdf@data, .(hexid, Sid), summarize,
                 nhr = length(Sid),
                 WmaxS = max(WmaxS))
Fstats2 = ddply(Fstats, .(hexid), summarize, 
                 count = length(Sid))
maxFSU$count = Fstats2$count    # dim = 66 by 18

Remove grids with too few TCs

maxWinds = subset(maxWinds, WmaxS >= 17)    # dim = 57 by 38
maxWinds = subset(maxWinds, count >= 15)    # dim = 33 by 38

maxFSU = subset(maxFSU, WmaxS >= 17)    # dim = 60 by 5
maxFSU = subset(maxFSU, count >= 15)    # dim = 32 by 5

Define data frame with relevant per hexagon variables

maxWinds = cbind(maxWinds$Yr, maxWinds$Mo, 
                 maxWinds$WmaxS, maxWinds$hexid,
                 maxWinds$count)
colnames(maxWinds) = c("Yr", "Mo", "WmaxS", "hexid", "count")
maxWinds = data.frame(maxWinds)
row.names(maxWinds) = paste("ID", maxWinds$hexid, sep="")
ids = row.names(maxWinds)

############## FSU COAPS ################

maxFSU = cbind(maxFSU$Yr, maxFSU$Mo, 
                 maxFSU$WmaxS, maxFSU$hexid, maxFSU$count)
colnames(maxFSU) = c("Yr", "Mo", "WmaxS", "hexid", "count")
maxFSU = data.frame(maxFSU)
row.names(maxFSU) = paste("ID", maxFSU$hexid, sep="")
FSUids = row.names(maxFSU)

Load and arrange SST data (OBS)

load("~/Dropbox/spatialextremes/Sarah/Space-time/mSST.RData")
A = mSST[,171:197]      
S = mSST[,238:264]
O = mSST[,305:331]
ASOsst = cbind(mSST[,403:404], A, S, O)

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

ASOsst.pdf = over(x=hpg[ids], y=ASOsst.sdf, fn=mean)
if(sum(is.na(ASOsst.pdf) > 0)){
  missValRows = unique(which(is.na(ASOsst.pdf), arr.ind=TRUE)[,1])
  ASOsst.pdf = ASOsst.pdf[-missValRows,]
  windSplit = windSplit[-missValRows]
}

SSTYearMonth = strsplit(substring(colnames(ASOsst.pdf), 2), "M")
SSTYear = as.numeric(sapply(SSTYearMonth, function(x) x[1]))

sst.yr = sapply(unique(SSTYear), 
      function(x) as.vector(rowMeans(ASOsst.pdf[,which(x==SSTYear)])))

sst.yr = data.frame(sst.yr)
rownames(sst.yr) = rownames(maxWinds)
colnames(sst.yr) = paste("Y", 1982:2008, sep="")

rm(mSST)

Arrange SST data (GCM)

FSUsst.pdf = over(x=hpg[FSUids], y=ASOsst.sdf, fn=mean)
if(sum(is.na(FSUsst.pdf) > 0)){
  missValRows = unique(which(is.na(FSUsst.pdf), arr.ind=TRUE)[,1])
  ASOsst.pdf = FSUsst.pdf[-missValRows,]
  fsuSplit = fsuSplit[-missValRows]
}

fsuYearMonth = strsplit(substring(colnames(FSUsst.pdf), 2), "M")
fsuYear = as.numeric(sapply(fsuYearMonth, function(x) x[1]))

fsu.sst.yr = sapply(unique(fsuYear), 
        function(x) as.vector(rowMeans(FSUsst.pdf[,which(x==fsuYear)])))

fsu.sst.yr = data.frame(fsu.sst.yr)

rownames(fsu.sst.yr) = rownames(maxFSU)
colnames(fsu.sst.yr) = paste("Y", 1982:2008, sep="")

Load PI data. PI_12z.RData is calculated from MERRA data over the 1982–2008 time period using 1000–50 hPa at 50 hPa intervals to be consistent with the FSU-COAPS model. PIfsu2.RData omits the same SST grids as are omitted in the MERRA PI calculations (for grids that are to close to land). PI was calculated from monthly means of 12z observations/model points using a version of Kerry Emanuel’s MATLAB code (ftp://texmex.mit.edu/pub/emanuel/TCMAX/) adapted for use in R by Tom Jagger.

load("~/Dropbox/spatialextremes/Sarah/PI_12z.RData")

all.ids = c(ids, FSUids)
all.ids = unique(all.ids)

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

PI.pdf = over(x=hpg[ids], y=PI.sdf, fn=max)
PI.pdf.comp = over(x=hpg[all.ids], y=PI.sdf, fn=mean)
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.pdf.comp) = cn
PI.pdf.comp$PI.merra = rowMeans(PI.pdf.comp)

############## FSU COAPS ################
load("~/Dropbox/spatialextremes/Sarah/PIfsu_12z.RData")

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

fsu.PI.pdf = over(x=hpg[FSUids], y=PIfsu.sdf, fn=max)
fsu.PI.pdf.comp = over(x=hpg[all.ids], y=PIfsu.sdf, fn=mean)
colnames(fsu.PI.pdf) = cn
colnames(fsu.PI.pdf.comp) = cn
fsu.PI.pdf.comp$PI.fsu = rowMeans(fsu.PI.pdf.comp)

Obtain per hexagon max MPI (OBS)

ym = character()
PI.pdf$MPI = do.call(pmax, PI.pdf)
for(i in 1:dim(PI.pdf)[1]){
  m = which(PI.pdf[i,1:81] == PI.pdf$MPI[i])
  ym[i] = colnames(PI.pdf[m])
}

PIYearMonth = strsplit(ym, "M")
PIYear = sapply(PIYearMonth, function(x) x[1])
PIYear = strsplit(PIYear, "Y")
PIYear = as.numeric(sapply(PIYear, function(x) x[2]))
PIMonth = as.numeric(sapply(PIYearMonth, function(x) x[2]))


maxWinds$MPI = PI.pdf$MPI
maxWinds$YrMPI = PIYear
maxWinds$MoMPI = PIMonth
maxWinds$sst = rowMeans(sst.yr)

Obtain per hexagon max MPI (GCM)

ym = character()
fsu.PI.pdf$MPI = do.call(pmax, fsu.PI.pdf)
for(i in 1:dim(fsu.PI.pdf)[1]){
  m = which(fsu.PI.pdf[i,1:81] == fsu.PI.pdf$MPI[i])
  ym[i] = colnames(fsu.PI.pdf[m])
}

fsuPIYearMonth = strsplit(ym, "M")
fsuPIYear = sapply(fsuPIYearMonth, function(x) x[1])
fsuPIYear = strsplit(fsuPIYear, "Y")
fsuPIYear = as.numeric(sapply(fsuPIYear, function(x) x[2]))
fsuPIMonth = as.numeric(sapply(fsuPIYearMonth, function(x) x[2]))


maxFSU$MPI = fsu.PI.pdf$MPI
maxFSU$YrMPI = fsuPIYear
maxFSU$MoMPI = fsuPIMonth
maxFSU$sst = rowMeans(fsu.sst.yr)

Create a spatial polygons data frame (for plotting purposes)

maxWinds$W.scale = maxWinds$WmaxS/mean(maxWinds$WmaxS)
rownames(maxWinds) = paste("ID", maxWinds$hexid, sep="")
wind.pdf = SpatialPolygonsDataFrame(hpg[rownames(maxWinds)], maxWinds)

maxFSU$W.scale = maxFSU$WmaxS/mean(maxFSU$WmaxS)
rownames(maxFSU) = paste("ID", maxFSU$hexid, sep="")
fsu.pdf = SpatialPolygonsDataFrame(hpg[rownames(maxFSU)], maxFSU)

Set up map

We begin by subdividing the North Atlantic basin into equal area hexagon cells. For consistency, these individual hexagons will be referred to as “regions” throughout the remainder of this document. Each region has an area of 731,935.3 km\(^{2}\), which is slightly larger than the U.S. state of Texas. Following this step, we overlay TC track data onto the hexagonal regions, as illustrated in Fig. 2a for best-track (BT) TC track data and in Fig. 2b for simulated track data. Both best-track and simulated data are overlaid onto the same grid to facilitate comparison. Figure 2c(d) depicts the per region maximum best-track(simulated) TC intensity over the 1982–2008 period. The strongest best-track TCs tend to occur over the Gulf of Mexico and reach intensities exceeding 70 m \(s^{-1}\), while the strongest simulated TCs occur in the central Atlantic and do not exceed 50 m \(s^{-1}\).

Figure 2 also displays the total number of unique TCs that passed over each region (i.e. the per region TC count) for best-track (e) and simulated (f) data. These figures clearly illustrate the preferred formation region for FSU/COAPS TCs over the southern-central portion of the North Atlantic basin. Here the model generates far more TCs than are observed. Conversly, the model generates and tracks far fewer TCs than are observed over the Gulf of Mexico. For the following analyses, we subset on those regions with TC counts of at least 15, where only those TCs with maximum wind speeds of at least 17 m s\(^{-1}\) are included in the count.

Fig. 2 The top row depicts all observed (a) and modeled (b) TC tracks from the 1982–2008 period overlaid onto the same grid of equal area hexagons. Shading is done in proportion to TC category on the Saffir-Simpson scale with blues representing tropical storms and reds representing strong hurricanes. The bottom row displays per region (i.e. hexagon) TC counts for observed (c) and modeled (d) TCs. Only hexagons with at least 15 observed TCs (a, c) or 15 modeled TCs (b, d) are displayed.

In addition to track data, we similarly we overlay the set of gridded ASO SSTs onto the same set of regions and compute the per region average ASO SST over the 1982–2008 time period (Fig. 3).

Fig. 3 Per region average August–October SST over the 1982–2008 time period. Units are given in \(^{\deg}C\) Only regions with at least 15 observed TCs are displayed.

We next calculate limiting and potential intensity. Limiting intensity is a statistical upper limit on per region TC maximum intensity that is estimated from the climatology of TC maximum intensities.

Estimate LI for observed and simulated TCs using a combined GPD-Poisson model.

oDat = stats
oDat = merge(oDat, maxWinds, by="hexid")
oDat = split(oDat, oDat$hexid)
obs.param = add.LI(oDat, 0.52, 1982, 2008)
maxWinds = cbind(maxWinds, "thresh"=obs.param$thresh,
                 "rate"=obs.param$rate, "lambda"=obs.param$lambda,
                 "sigma"=obs.param$sigma, "xi"=obs.param$xi,
                 "theo"=obs.param$theo)

fDat = Fstats
fDat = merge(fDat, maxFSU, by="hexid")
fDat = split(fDat, fDat$hexid)
fsu.param = add.LI(fDat, 0.52, 1982, 2008)
maxFSU = cbind(maxFSU, "thresh"=fsu.param$thresh,
               "rate"=fsu.param$rate, "lambda"=fsu.param$lambda,
               "sigma"=fsu.param$sigma, "xi"=fsu.param$xi,
               "theo"=fsu.param$theo)

maxWinds.all = maxWinds
maxFSU.all = maxFSU

Organize data to display using ggplot

maxWinds$lat = coordinates(wind.pdf)[,2]
maxFSU$lat = coordinates(fsu.pdf)[,2]

# Remove hexagons with xi > 0
b1 = which(maxWinds$hexid==77)
b2 = which(maxWinds$hexid==110)
b3 = which(maxWinds$xi > 0)
bad = c(b1, b2, b3)
maxWinds = maxWinds[!maxWinds$hexid==77,]
maxWinds = maxWinds[!maxWinds$hexid==110,]
maxWinds = subset(maxWinds, xi < 0)

obDat = cbind(maxWinds$Yr, maxWinds$WmaxS, maxWinds$hexid, 
              maxWinds$count, maxWinds$sst)
colnames(obDat) = c("Year", "Wind", "Hexid", "Count", "SST")
obDat = data.frame(obDat)
obDat$type = as.factor(rep("W_max", dim(obDat)[1]))
obDat$id = rep(1:dim(obDat)[1])

piDat = cbind(maxWinds$YrMPI, maxWinds$MPI, maxWinds$hexid, 
              maxWinds$count, maxWinds$sst)
colnames(piDat) = c("Year", "Wind", "Hexid", "Count", "SST")
piDat = data.frame(piDat)
piDat$type = as.factor(rep("PI", dim(piDat)[1]))
piDat$id = obDat$id
piDat$PIrank = rank(-piDat$Wind)

liDat = cbind(maxWinds$Yr, maxWinds$theo, maxWinds$hexid,
              maxWinds$count, maxWinds$sst)
colnames(liDat) = c("Year", "Wind", "Hexid", "Count", "SST")
liDat = data.frame(liDat)
liDat$type = as.factor(rep("LI", dim(liDat)[1]))
liDat$id = obDat$id
liDat$PIrank = piDat$PIrank
obDat$PIrank = piDat$PIrank

allDat = rbind(obDat, piDat, liDat)
allDat$data = rep("Best-track (Observed)", dim(allDat)[1])

Organize data to display using ggplot

maxFSU = maxFSU[!maxFSU$hexid==111,]

moDat = cbind(maxFSU$Yr, maxFSU$WmaxS, maxFSU$hexid, 
              maxFSU$count, maxFSU$sst)
colnames(moDat) = c("Year", "Wind", "Hexid", "Count", "SST")
moDat = data.frame(moDat)
moDat$type = as.factor(rep("W_max", dim(moDat)[1]))
moDat$id = rep(1:dim(moDat)[1])

m.piDat = cbind(maxFSU$YrMPI, maxFSU$MPI, maxFSU$hexid, 
                maxFSU$count, maxFSU$sst)
colnames(m.piDat) = c("Year", "Wind", "Hexid", "Count", "SST")
m.piDat = data.frame(m.piDat)
m.piDat$type = as.factor(rep("PI", dim(m.piDat)[1]))
m.piDat$id = moDat$id
m.piDat$PIrank = rank(-m.piDat$Wind)

m.liDat = cbind(maxFSU$Yr, maxFSU$theo, maxFSU$hexid, 
                maxFSU$count, maxFSU$sst)
colnames(m.liDat) = c("Year", "Wind", "Hexid", "Count", "SST")
m.liDat = data.frame(m.liDat)
m.liDat$type = as.factor(rep("LI", dim(m.liDat)[1]))
m.liDat$id = moDat$id
m.liDat$PIrank = m.piDat$PIrank
moDat$PIrank = m.piDat$PIrank

m.allDat = rbind(moDat, m.piDat, m.liDat)
m.allDat$data = rep("FSU-COAPS", dim(m.allDat)[1])

Create another spatial polygons data frame for plotting purposes

maxWinds$LI.scale = maxWinds$theo/mean(maxWinds$theo)
rownames(maxWinds) = paste("ID", maxWinds$hexid, sep="")
wind.pdf = SpatialPolygonsDataFrame(hpg[rownames(maxWinds)], maxWinds)

maxFSU$LI.scale = maxFSU$theo/mean(maxFSU$theo)
rownames(maxFSU) = paste("ID", maxFSU$hexid, sep="")
fsu.pdf = SpatialPolygonsDataFrame(hpg[rownames(maxFSU)], maxFSU)

Maximum, limiting, and potential intensity: FSU/COAPS vs. Reanalysis and observations

Fig. 4 Per region limiting intensity (LI) estimated from observed (a) and model-generated (b) TC data. The bottom panels depict per region limiting intensity scaled by the mean per region limiting intensity. Scaled values below 1 indicate a limiting intensity that below the basin wide average value, while values above 1 indicate a limiting intensity that is above the basin wide average value. Limiting intensity is epxressed in units of m/s.

Load SLP data for plotting purposes

load("~/Dropbox/SpatialExtremes/Sarah/slpData.RData")
load("~/Dropbox/SpatialExtremes/Sarah/slpFSU.RData")

slpData$mSLP = rowMeans(slpData[,3:83])
slpData$mSLP = slpData$mSLP/100
slpFSU$mSLP = rowMeans(slpFSU[,3:83])

lon = seq(-110.625, 0.625, 1.25)
lat = seq(-3.125, 70.625, 1.25)
slp = slpData$mSLP
slp.mat = matrix(slp, ncol=60)

lonF = unique(slpFSU$lon)
latF = unique(slpFSU$lat)
slpF = slpFSU$mSLP
slp.matF = matrix(slpF, ncol=75)

Fig. 5 Raw sea level pressure data from MERRA (left) and the FSU/COAP model (right). Contouring is shown at 2 hPa intervals.’

rm(slpData)
rm(slpFSU)

Fig. 6 Per region PI estimated from MERRA (a) and simulated (b) data. The bottom panels depict per region PI scaled by the mean regional PI over the entire North Atlantic. Scaled values below 1 indicate an PI that below the basin wide average PI, while values above 1 indicate an PI that is above the basin wide average PI. PI is epxressed in units of m/s

Fig. 7 Bar plots comparing per hexagon maximum observed intensity (black), maximum potential intensity (red), and limiting intensity (blue). Each triad of black, red, and blue bars refers to a different region. Regions are ordered from highest to lowest PI

Loading 100 hPa temperature from MERRA and FSU/COAPS

load("~/Dropbox/SpatialExtremes/Sarah/tFSU.RData")
load("~/Dropbox/SpatialExtremes/Sarah/Space-time/ttData.RData")

PI.comp = data.frame(cbind(PI.pdf.comp$PI.merra, fsu.PI.pdf.comp$PI.fsu))
colnames(PI.comp) = c("PI.merra", "PI.fsu")
rownames(PI.comp) = rownames(PI.pdf.comp)

Format MERRA 100 hPa temperatures, with each row representing a specific region/hexagon, and each column representing a specific year. Note: I may not be able to use this given the high correlation to SST. I initially used 70 hPa temps as this level was the coldest; however, because the FSU/COAPS model does not include a 70 hPa level, I will use MERRA 100 hPa temperatures for a more fair comparison with the model.

tt = ttData[,3:308]
tt100col = seq(3, 306, 3)
tt100 = tt[,tt100col]
tt100 = tt100[,10:90]
ttData = cbind(ttData[,1:2], tt100)

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

tt.pdf = over(x=hpg[all.ids], y=ttData.sdf, fn=mean)
ttYearMonth = strsplit(substring(colnames(tt.pdf), 6), "M")
ttYear = as.numeric(sapply(ttYearMonth, function(x) x[1]))

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

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

tt.yr = data.frame(tt.yr)
PI.comp$TT.merra = rowMeans(tt.yr)

Format FSU/COAPS 100 hPa temperatures

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

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

ttFSU.pdf = over(x=hpg[all.ids], y=ttFSU.sdf, fn=mean)
ttYearMonth = strsplit(substring(colnames(ttFSU.pdf), 6), "M")
ttYear = as.numeric(sapply(ttYearMonth, function(x) x[1]))

tt.fsu.yr = sapply(unique(ttYear), 
    function(x) as.vector(rowMeans(ttFSU.pdf[,which(x==ttYear)])))

dimnames(tt.fsu.yr) = list(id=rownames(ttFSU.pdf), 
          Year=paste("Y",unique(ttYear), sep="")) 

tt.fsu.yr = data.frame(tt.fsu.yr)
PI.comp$TT.fsu = rowMeans(tt.fsu.yr)

Create spatial polygons data frame for plotting purposes

PI.comp$PI.diff = PI.comp$PI.merra - PI.comp$PI.fsu
PI.comp$TT.diff = PI.comp$TT.merra - PI.comp$TT.fsu

PIcomp.pdf = SpatialPolygonsDataFrame(hpg[all.ids], PI.comp)

Fig. 8 Per region MERRA - FSU/COAPS mean PI in units of m/s (a), and per region MERRA - FSU/COAPS mean 100 hPa temperature in units of K (b). Both temperature and PI means are taken over the entire 1982–2008 time period

The sensitivity of potential intensity to SST

Calculate sensitivity of LI, Wmax, and PI for observed TCs:

maxWinds25 = subset(maxWinds, sst >= 25)
maxFSU25 = subset(maxFSU, sst >= 25)

lm.w = summary(lm(WmaxS ~ sst, data=maxWinds))
lm.w.25 = summary(lm(WmaxS ~ sst, data=maxWinds25))
Wmax.res = residuals(lm.w.25)

lm.pi = summary(lm(MPI ~ sst, data=maxWinds))
lm.pi.25 = summary(lm(MPI ~ sst, data=maxWinds25))
MPI.res = residuals(lm.pi.25)

lm.li = summary(lm(theo ~ sst, data=maxWinds))
lm.li.25 = summary(lm(theo ~ sst, data=maxWinds25))
theo.res = residuals(lm.li.25)

Calculate the sensitivity of Wmax, LI, and PI to SST for FSU/COAPS TCs:

lm.w.fsu = summary(lm(WmaxS ~ sst, data=maxFSU))
lm.w.25.fsu = summary(lm(WmaxS ~ sst, data=maxFSU25))
Wmax.res.fsu = residuals(lm.w.25.fsu)

lm.pi.fsu = summary(lm(MPI ~ sst, data=maxFSU))
lm.pi.25.fsu = summary(lm(MPI ~ sst, data=maxFSU25))
MPI.res.fsu = residuals(lm.pi.25.fsu)

lm.li.fsu = summary(lm(theo ~ sst, data=maxFSU))
lm.li.25.fsu = summary(lm(theo ~ sst, data=maxFSU25))
theo.res.fsu = residuals(lm.li.25.fsu)

Fig. 9 The sensitivity of per hexagon maximum, maximum potential, and limiting intensity to SST for observed (BT) and simulated (FSU/COAPS) TCs. The shaded area gives the 95% confidence interval on the regressions

maxWinds25 = subset(maxWinds, sst >= 25)
maxFSU25 = subset(maxFSU, sst >= 25)

maxWinds25$Wmax.res = Wmax.res
maxWinds25$MPI.res = MPI.res
maxWinds25$theo.res = theo.res

maxFSU25$Wmax.res = Wmax.res.fsu
maxFSU25$MPI.res = MPI.res.fsu
maxFSU25$theo.res = theo.res.fsu

wind25.pdf = SpatialPolygonsDataFrame(hpg[rownames(maxWinds25)], 
                                      maxWinds25)
FSU25.pdf = SpatialPolygonsDataFrame(hpg[rownames(maxFSU25)], 
                                      maxFSU25)

Fig. 10 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

The sensitivity of maximum intensity to atmospheric variables

load("~/Dropbox/SpatialExtremes/Sarah/tFSU.RData")
load("~/Dropbox/SpatialExtremes/Sarah/vwsFSU.RData")
load("~/Dropbox/SpatialExtremes/Sarah/rhFSU.RData")

load("~/Dropbox/SpatialExtremes/Sarah/Space-time/ttData.RData")
load("~/Dropbox/SpatialExtremes/Sarah/Space-time/rhData.RData")
load("~/Dropbox/SpatialExtremes/Sarah/Space-time/vwsData.RData")
if(sum(is.na(vwsData)) > 0){
missingVals = which(is.na(vwsData), arr.ind=TRUE)   
badRows = unique(missingVals[,1])
vwsData = vwsData[-badRows,]
}

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

vws.pdf = over(x=hpg[rownames(maxWinds)], y=vwsData.sdf, fn=mean)
vws.pdf.a = over(x=hpg[rownames(maxWinds.all)], 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)])))
vws.yr.a = sapply(unique(vwsYear), 
    function(x) as.vector(rowMeans(vws.pdf.a[,which(x==vwsYear)])))

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

vws.yr = data.frame(vws.yr)
vws.yr.a = data.frame(vws.yr.a)
vws.yr = vws.yr[,4:30]
vws.yr.a = vws.yr.a[,4:30]
coordinates(rhData) = c("lon", "lat")
proj4string(rhData) = CRS(ll)
rhData.sdf = spTransform(rhData, CRS(lcc))

rh.pdf = over(x=hpg[rownames(maxWinds)], y=rhData.sdf, fn=mean)
rh.pdf.a = over(x=hpg[rownames(maxWinds.all)], y=rhData.sdf, fn=mean)
rhYearMonth = strsplit(substring(colnames(rh.pdf), 4), "M")
rhYear = as.numeric(sapply(rhYearMonth, function(x) x[1]))

rh.yr = sapply(unique(rhYear), 
      function(x) as.vector(rowMeans(rh.pdf[,which(x==rhYear)])))
rh.yr.a = sapply(unique(rhYear), 
      function(x) as.vector(rowMeans(rh.pdf.a[,which(x==rhYear)])))

dimnames(rh.yr) = list(id=rownames(rh.pdf), 
                  Year=paste("Y",unique(rhYear), sep=""))
dimnames(rh.yr.a) = list(id=rownames(rh.pdf.a), 
                  Year=paste("Y",unique(rhYear), sep=""))

rh.yr = data.frame(rh.yr)
rh.yr.a = data.frame(rh.yr.a)
rh.yr = rh.yr[,4:30]
rh.yr.a = rh.yr.a[,4:30]
tt = ttData[,3:308]
tt100col = seq(3, 306, 3)
tt100 = tt[,tt100col]
ttData = cbind(ttData[,1:2], tt100)

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

tt.pdf = over(x=hpg[rownames(maxWinds)], y=ttData.sdf, fn=mean)
tt.pdf.a = over(x=hpg[rownames(maxWinds.all)], y=ttData.sdf, fn=mean)
ttYearMonth = strsplit(substring(colnames(tt.pdf), 6), "M")
ttYear = as.numeric(sapply(ttYearMonth, function(x) x[1]))

tt.yr = sapply(unique(ttYear), 
    function(x) as.vector(rowMeans(tt.pdf[,which(x==ttYear)])))
tt.yr.a = sapply(unique(ttYear), 
    function(x) as.vector(rowMeans(tt.pdf.a[,which(x==ttYear)])))

dimnames(tt.yr) = list(id=rownames(tt.pdf), 
          Year=paste("Y",unique(ttYear), sep=""))
dimnames(tt.yr.a) = list(id=rownames(tt.pdf.a), 
          Year=paste("Y",unique(ttYear), sep=""))

tt.yr = data.frame(tt.yr)
tt.yr.a = data.frame(tt.yr.a)
tt.yr = tt.yr[,4:30]
tt.yr.a = tt.yr.a[,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))

vwsFSU.pdf.a = over(x=hpg[rownames(maxWinds.all)], y=vwsFSU.sdf, fn=mean)
vwsFSU.pdf = over(x=hpg[rownames(maxFSU)], y=vwsFSU.sdf, fn=mean)
vwsYearMonth = strsplit(substring(colnames(vwsFSU.pdf), 5), "M")
vwsYear = as.numeric(sapply(vwsYearMonth, function(x) x[1]))

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

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

vws.fsu.yr = data.frame(vws.fsu.yr)
vws.fsu.yr.a = data.frame(vws.fsu.yr.a)
rhFSU = cbind(tFSU[,1:2], rhFSU)
be = rep("rhY", 81)
ye = rep(1982:2008, each=3)
m = rep(c("M08", "M09", "M10"), 27)
colnames(rhFSU) = c("lon", "lat", paste(be, ye, m, sep=""))

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

rhFSU.pdf = over(x=hpg[rownames(maxFSU)], y=rhFSU.sdf, fn=mean)
rhFSU.pdf.a = over(x=hpg[rownames(maxWinds.all)], y=rhFSU.sdf, fn=mean)
rhYearMonth = strsplit(substring(colnames(rhFSU.pdf), 4), "M")
rhYear = as.numeric(sapply(rhYearMonth, function(x) x[1]))

rh.fsu.yr = sapply(unique(rhYear), 
    function(x) as.vector(rowMeans(rhFSU.pdf[,which(x==rhYear)])))
rh.fsu.yr.a = sapply(unique(rhYear), 
    function(x) as.vector(rowMeans(rhFSU.pdf.a[,which(x==rhYear)])))

dimnames(rh.fsu.yr) = list(id=rownames(rhFSU.pdf), 
                  Year=paste("Y",unique(rhYear), sep="")) 
dimnames(rh.fsu.yr.a) = list(id=rownames(rhFSU.pdf.a), 
                  Year=paste("Y",unique(rhYear), sep="")) 

rh.fsu.yr = data.frame(rh.fsu.yr)
rh.fsu.yr.a = data.frame(rh.fsu.yr.a)
be = rep("t100Y", 81)
ye = rep(1982:2008, each=3)
m = rep(c("M08", "M09", "M10"), 27)
colnames(tFSU) = c("lon", "lat", paste(be, ye, m, sep=""))

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

ttFSU.pdf = over(x=hpg[rownames(maxFSU)], y=ttFSU.sdf, fn=mean)
ttFSU.pdf.a = over(x=hpg[rownames(maxWinds.all)], y=ttFSU.sdf, fn=mean)
ttYearMonth = strsplit(substring(colnames(ttFSU.pdf), 6), "M")
ttYear = as.numeric(sapply(ttYearMonth, function(x) x[1]))

tt.fsu.yr = sapply(unique(ttYear), 
    function(x) as.vector(rowMeans(ttFSU.pdf[,which(x==ttYear)])))
tt.fsu.yr.a = sapply(unique(ttYear), 
    function(x) as.vector(rowMeans(ttFSU.pdf.a[,which(x==ttYear)])))

dimnames(tt.fsu.yr) = list(id=rownames(ttFSU.pdf), 
          Year=paste("Y",unique(ttYear), sep="")) 
dimnames(tt.fsu.yr.a) = list(id=rownames(ttFSU.pdf.a), 
          Year=paste("Y",unique(ttYear), sep="")) 

tt.fsu.yr = data.frame(tt.fsu.yr)
tt.fsu.yr.a = data.frame(tt.fsu.yr.a)
vws.o = rowMeans(vws.yr.a)
vws.f = rowMeans(vws.fsu.yr.a)
rh.o = rowMeans(rh.yr.a)*100
rh.f = rowMeans(rh.fsu.yr.a)*100
tt.o = rowMeans(tt.yr.a)
tt.f = rowMeans(tt.fsu.yr.a)

atm.data = data.frame(cbind(vws.o, rh.o, tt.o, 
                            vws.f, rh.f, tt.f))
atm.data$vws.diff = atm.data$vws.o - atm.data$vws.f
atm.data$rh.diff = atm.data$rh.o - atm.data$rh.f
atm.data$tt.diff = atm.data$tt.o - atm.data$tt.f

atm.data.pdf = SpatialPolygonsDataFrame(hpg[rownames(maxWinds.all)],
                                    atm.data)

Set up atmospheric variables for observed TCs

covar = matrix(0, nrow=dim(maxWinds)[1], ncol=3)
covar = data.frame(covar)

for(i in 1:dim(maxWinds)[1]){
  y = which(colnames(vws.yr) == paste("Y", maxWinds$Yr[i], sep=""),
            arr.ind=TRUE)
  covar[i,1] = vws.yr[i, y]
  covar[i,2] = rh.yr[i, y]*100
  covar[i,3] = tt.yr[i, y]
}

colnames(covar) = c("shear", "rh", "Ttemp")

maxWinds = cbind(maxWinds, covar)

maxWinds.pdf = SpatialPolygonsDataFrame(hpg[rownames(maxWinds)],
                                    maxWinds)

Set up atmospheric variables for FSU/COAPS TCs

covar = matrix(0, nrow=dim(maxWinds)[1], ncol=3)
covar = data.frame(covar)

for(i in 1:dim(maxFSU)[1]){
  y = which(colnames(vws.fsu.yr) == paste("Y", maxFSU$Yr[i], sep=""),
            arr.ind=TRUE)
  covar[i,1] = vws.fsu.yr[i, y]
  covar[i,2] = rh.fsu.yr[i, y]*100
  covar[i,3] = tt.fsu.yr[i, y]
}

colnames(covar) = c("shear", "rh", "Ttemp")

maxFSU = cbind(maxFSU, covar)
maxFSU.pdf = SpatialPolygonsDataFrame(hpg[rownames(maxFSU)],
                                    maxFSU)

Fig. 11 Spatial patterns of atmospheric variables

Fig. 12 Difference Maps

require(corrplot)
## Loading required package: corrplot
dat = subset(maxWinds, sst > 25)
dat2 = subset(maxFSU, sst > 25)

cor.test(dat$rh, dat$sst)     # 0.542 (0.196, 0.768)
## 
##  Pearson's product-moment correlation
## 
## data:  dat$rh and dat$sst
## t = 3.1616, df = 24, p-value = 0.004214
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1960751 0.7682359
## sample estimates:
##       cor 
## 0.5422436
cor.test(dat2$rh, dat2$sst)   # 0.178 (-0.217, 0.523)
## 
##  Pearson's product-moment correlation
## 
## data:  dat2$rh and dat2$sst
## t = 0.9045, df = 25, p-value = 0.3743
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.2166532  0.5226728
## sample estimates:
##       cor 
## 0.1780172
cor.test(dat$shear, dat$sst)     # -0.503 (-0.144, -0.745)
## 
##  Pearson's product-moment correlation
## 
## data:  dat$shear and dat$sst
## t = -2.854, df = 24, p-value = 0.008755
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.7453926 -0.1441338
## sample estimates:
##        cor 
## -0.5033817
cor.test(dat2$shear, dat2$sst)   # -0.487 (-0.131, -0.731)
## 
##  Pearson's product-moment correlation
## 
## data:  dat2$shear and dat2$sst
## t = -2.786, df = 25, p-value = 0.01003
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.7314603 -0.1309526
## sample estimates:
##        cor 
## -0.4867443
cor.test(dat$Ttemp, dat$sst)     # -0.767 (-0.538, -0.889)
## 
##  Pearson's product-moment correlation
## 
## data:  dat$Ttemp and dat$sst
## t = -5.8326, df = 24, p-value = 5.149e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.8893088 -0.5379381
## sample estimates:
##        cor 
## -0.7657332
cor.test(dat2$Ttemp, dat2$sst)   # -0.844 (-0.683, -0.927)
## 
##  Pearson's product-moment correlation
## 
## data:  dat2$Ttemp and dat2$sst
## t = -7.863, df = 25, p-value = 3.215e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.9266926 -0.6827767
## sample estimates:
##        cor 
## -0.8438432
cor.test(dat$rh, dat$WmaxS)      #  0.476 (0.109, 0.728)
## 
##  Pearson's product-moment correlation
## 
## data:  dat$rh and dat$WmaxS
## t = 2.6512, df = 24, p-value = 0.01398
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1086177 0.7289155
## sample estimates:
##       cor 
## 0.4759451
cor.test(dat2$rh, dat2$WmaxS)    # -0.705 (-0.856, -0.443)
## 
##  Pearson's product-moment correlation
## 
## data:  dat2$rh and dat2$WmaxS
## t = -4.9644, df = 25, p-value = 4.087e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.8555240 -0.4432373
## sample estimates:
##        cor 
## -0.7045752
cor.test(dat$shear, dat$WmaxS)   # -0.547 (-0.202, -0.771)
## 
##  Pearson's product-moment correlation
## 
## data:  dat$shear and dat$WmaxS
## t = -3.1988, df = 24, p-value = 0.003852
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.7708355 -0.2021967
## sample estimates:
##        cor 
## -0.5467282
cor.test(dat2$shear, dat2$WmaxS) # -0.102 (-0.464, 0.289)
## 
##  Pearson's product-moment correlation
## 
## data:  dat2$shear and dat2$WmaxS
## t = -0.5115, df = 25, p-value = 0.6135
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.4638438  0.2894393
## sample estimates:
##        cor 
## -0.1017682
cor.test(dat$Ttemp, dat$WmaxS)   #  -0.4881 (-0.115, -0.732)
## 
##  Pearson's product-moment correlation
## 
## data:  dat$Ttemp and dat$WmaxS
## t = -2.6902, df = 24, p-value = 0.01279
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.7321807 -0.1155322
## sample estimates:
##        cor 
## -0.4813432
cor.test(dat2$Ttemp, dat2$WmaxS)    #  0.199 (-0.196, 0.538)
## 
##  Pearson's product-moment correlation
## 
## data:  dat2$Ttemp and dat2$WmaxS
## t = 1.013, df = 25, p-value = 0.3208
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1962560  0.5379856
## sample estimates:
##       cor 
## 0.1985672
cmat.o = cbind(dat$sst, dat$shear, dat$rh, dat$Ttemp)
colnames(cmat.o) = c("SST", "Shear", "RH", "100 hPa T")
cmat.m = cbind(dat2$sst, dat2$shear, dat2$rh, dat2$Ttemp)
colnames(cmat.m) = c("SST", "Shear", "RH", "100 hPa T")

cr2 = rev(brewer.pal(10, "RdBu"))

M1 = cor(cmat.o)
M2 = cor(cmat.m)
par(mfrow=c(1,2))
corrplot.mixed(M1, cl.ratio=0.3, col=cr2, 
         tl.pos="d",
         tl.cex=0.9,
         tl.col="black")
corrplot.mixed(M2, cl.ratio=0.3, col=cr2, 
         tl.pos="d", 
         tl.cex=0.9,
         tl.col="black")

Multiple Regressions: We first using atmospheric covariate data corresponding the ASO average from the year during which the per region maximum TC intensity occurred.

summary(lm(WmaxS ~ sst + rh + shear, data=dat))
# Coeff = 6.6  | -0.0601 | -0.6679  
# S.E. = 1.24  |  0.187  | 0.507
# t-val = 5.33 | -0.322  | -1.32
# M R^2 = 0.6702  (adjusted)

summary(lm(MPI ~ sst + rh + shear, data=dat))
# Coeff = 4.0  | -0.0520 | -0.0372  
# S.E. = 0.424 |  0.0640 |  0.174
# t-val = 9.34 | -0.812  | -0.215
# M R^2 = 0.8288


summary(lm(WmaxS ~ rh, data=dat))
# Coeff = 0.59,  S.E. = 0.221,  t-val = 2.65, R^2 = 0.2265
summary(lm(WmaxS ~ shear, data=dat))
# Coeff = -1.88,  S.E. = 0.587,  t-val = -3.20, R^2 = 0.2989
summary(lm(WmaxS ~ Ttemp, data=dat))
# Coeff = -1.43,  S.E. = 0.531,  t-val = -2.69, R^2 = 0.2317

summary(lm(MPI ~ rh, data=dat))
# Coeff = 0.26,  S.E. = 0.107,  t-val = 2.42, R^2 = 0.1961
summary(lm(MPI ~ shear, data=dat))
# Coeff = -0.74,  S.E. = 0.297,  t-val = -2.48, R^2 = 0.2040
summary(lm(MPI ~ Ttemp, data=dat))
# Coeff = -0.744,  S.E. = 0.245,  t-val = -3.04, R^2 = 0.2783

Same as above, but for the FSU/COAPS model.

summary(lm(WmaxS ~ sst + rh + shear, data=dat2))
# Coeff = 0.27  | -0.37   | -0.11  
# S.E. =  0.813 |  0.0771 |  0.229
# t-val = 0.331 | -4.80   | -0.479
# M R^2 = 0.5103  (adjusted)

summary(lm(MPI ~ sst + rh + shear, data=dat2))
# Coeff = 4.4  | -0.0518 |  0.41  
# S.E. = 0.431 |  0.0409 |  0.122
# t-val = 10.2 | -1.27   |  3.38
# M R^2 = 0.8261


summary(lm(WmaxS ~ rh, data=dat2))
# Coeff = -0.36,  S.E. = 0.073,  t-val = -4.96, R^2 = 0.496
summary(lm(WmaxS ~ shear, data=dat2))
# Coeff = -0.14,  S.E. = 0.272,  t-val = -0.511, R^2 = 0.0104
summary(lm(WmaxS ~ Ttemp, data=dat2))
# Coeff =  0.42,  S.E. = 0.412,  t-val = 1.01, R^2 = 0.0394

summary(lm(MPI ~ rh, data=dat2))
# Coeff = 0.032,  S.E. = 0.0919,  t-val = 0.35, R^2 = 4.89e-3
summary(lm(MPI ~ shear, data=dat2))
# Coeff = -0.20,  S.E. = 0.240,  t-val = -0.851, R^2 = 0.0282
summary(lm(MPI ~ Ttemp, data=dat2))
# Coeff = -1.1,  S.E. = 0.306,  t-val = -3.54, R^2 = 0.334