As we strive to better understand the relationship between tropical cyclones (TCs) and climate, global climate models (GCMs) have become a valuable tool and resource. However, despite notable progress in physics and resolution, GCMs are still unable to realistically simulate the inner core structure of intense TCs. Given this limitation, we must exercise caution when using GCMs to predict and understand how climate change may affect future TC activity, particularly with respect to intensity. Although climate models do not yet realistically simulate TC intensity and structure, they do a relatively better job of capturing the large scale environmental fields that influence how intense a TC may become (cite). Therefore, if we can use observational data to understand how various environmental variables (e.g., sea surface temperature, vertical wind shear, mid-level relative humidity) influence the lifetime maximum intensity of TCs, we may be able to utilize GCM simulations of relevant environmental variables to infer information about TC intensity in the future.
setwd("C:/Users/Sarah/Dropbox/SpatialExtremes")
require(ggplot2)
require(rgdal)
require(ismev)
require(ncdf)
library(mapproj)
library(mapdata)
library(maptools)
library(maps)
library(grid)
require(plyr)
library(colorRamps)
library(RColorBrewer)
require(spdep)
require(psych)
require(boot)
source("C:/Users/Sarah/Dropbox/SpatialExtremes/selectData.R")
source("C:/Users/Sarah/Dropbox/SpatialExtremes/ismevplus.R")
source("C:/Users/Sarah/Dropbox/SpatialExtremes/Kerrycanes/sensPoly.R")
source("C:/Users/Sarah/Dropbox/SpatialExtremes/ohcFunctions.R")
source("C:/Users/Sarah/Dropbox/SpatialExtremes/sensitivityUncertainty.R")
source("C:/Users/Sarah/Dropbox/SpatialExtremes/sensTrends.R")
source("C:/Users/Sarah/Dropbox/SpatialExtremes/Sarah/Space-Time/arrangeMERRA.R")
Load best track data, subset on tropical cycloens, 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/Research/Data/best.use2.RData")
begin = 1979
end = 2012
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 * 0.5144
quantile(Wind.df$WmaxS, prob = seq(0.6, 0.8, 0.01))
## 60% 61% 62% 63% 64% 65% 66% 67% 68% 69% 70% 71%
## 32.97 33.27 33.43 33.52 33.72 34.08 34.53 34.95 35.46 35.90 36.21 36.68
## 72% 73% 74% 75% 76% 77% 78% 79% 80%
## 37.24 37.79 38.33 38.66 39.10 39.79 40.51 41.08 41.60
Wind.df$maguv = Wind.df$maguv * 0.5144
Wind.df$WmaxS = Wind.df$WmaxS - Wind.df$maguv * 0.6
length(unique(Wind.df$Sid)) # 410
## [1] 410
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))
Define a set of equal area (~543 515 km2) hexagons that completely cover the set of TC tracks
set.seed(1025)
bb = matrix(0, nrow = 2, ncol = 2)
bb[2, 1] = bbox(Wind.sdf)[2, 1] * 0.8
bb[1, 1] = bbox(Wind.sdf)[1, 1] * 1.35
bb[1, 2] = bbox(Wind.sdf)[1, 2] * 1.08
bb[2, 2] = bbox(Wind.sdf)[2, 2] * 1.08
hpt = spsample(Wind.sdf, type = "hexagonal", n = 100, bb = bb)
hpg = HexPoints2SpatialPolygons(hpt)
# plot(hpg) plot(Wind.sdf, add=TRUE, pch='.', col='red')
np = length(hpg@polygons) # 77
area = hpg@polygons[[1]]@area * 10^-6
# 815271.9 km^2
Overlay TC data onto hexagons
hexid = over(x = Wind.sdf, y = hpg)
Wind.sdf$hexid = hexid
Split the TC data by hexagon ID, subset on only those hexagons with at least 15 TCs, find the per hexagon wind max, and finally define a matrix that gives the number of TCs per year for each hexagon (each row corresponds to a specific hexagon and each column specifies a year).
years = 1979:2012
Wind.split = split(Wind.sdf@data, Wind.sdf$hexid)
less = which(sapply(Wind.split, function(x) length(unique(x$Sid)) < 15))
Wind.split = Wind.split[-less]
Wind.max = lapply(Wind.split, function(x) get.max(x, maxfield = "WmaxS"))
Wind.max = do.call(rbind.data.frame, Wind.max) # dim = 1730 X 37
Wind.max = subset(Wind.max, WmaxS >= 17) # dim = 1410 X 37
Wind.max = split(Wind.max, Wind.max$hexid)
Wind.count = t(sapply(Wind.max, function(x) table(factor(x$Yr, level = years))))
Define a data frame similar to the per hexagon per year counts matrix defined above, but this time, record the per hexagon per year maximum TC wind speed
Wind.max.yr = matrix(0, nrow = length(Wind.max), ncol = length(1979:2012))
Wind.max.yr = data.frame(Wind.max.yr)
# Is there a less loopy way to do this? Look into it...though loops do
# the trick
for (i in 1:length(Wind.max)) {
y = 1979
for (j in 1:length(1979:2012)) {
if (Wind.count[i, j] > 0) {
all = Wind.max[[i]]
all = subset(all, Yr == y)
maxw = max(all$WmaxS)
Wind.max.yr[i, j] = maxw
} else {
Wind.max.yr[i, j] = 0
}
y = y + 1
}
}
Add row and column names to the count and wind speed data frames. Also, remove hexagons over land (47/48)
colnames(Wind.count) = paste("Y", colnames(Wind.count), sep = "")
colnames(Wind.max.yr) = colnames(Wind.count)
rownames(Wind.count) = paste("ID", rownames(Wind.count), sep = "")
rownames(Wind.max.yr) = rownames(Wind.count)
Wind.count = data.frame(Wind.count)
Wind.count = Wind.count[rownames(Wind.count) != "ID47", ]
Wind.count = Wind.count[rownames(Wind.count) != "ID48", ]
Wind.max.yr = Wind.max.yr[rownames(Wind.max.yr) != "ID47", ]
Wind.max.yr = Wind.max.yr[rownames(Wind.max.yr) != "ID48", ]
Create spatial polygons data frames for the count and wind speed data frames
Wind.count.pdf = SpatialPolygonsDataFrame(hpg[rownames(Wind.count)], Wind.count)
Wind.count.pdf$all = rowSums(Wind.count.pdf@data)
Wind.count.pdf$allRate = Wind.count.pdf$all/length(years)
Wind.max.yr.pdf = SpatialPolygonsDataFrame(hpg[rownames(Wind.max.yr)], Wind.max.yr)
Set up map borders
cl = map("world", xlim = c(-115, 15), ylim = c(-5, 70), resolution = 100)
clp = map2SpatialLines(cl, proj4string = CRS(ll))
clp = spTransform(clp, CRS(lcc))
l1 = list("sp.lines", clp, col = "black", lwd = 0.25)
Utilize the “spacetime” package to better visualize the spatial/temporal structure of the data. First, we create space-time data frames for the count/wind speed data
require(spacetime)
## Loading required package: spacetime
## Warning: package 'spacetime' was built under R version 3.0.3
# Create a long table of data...with space varying faster?
hex = hpg[rownames(Wind.max.yr)]
dates = paste("01jan", years, sep = "")
dates = as.Date(dates, format = "%d%B%Y")
wlong = reshape(Wind.max.yr, v.names = "Wmax", times = dates, varying = colnames(Wind.max.yr),
ids = rownames(Wind.max.yr), direction = "long")
wData = data.frame(values = wlong$Wmax, ID = wlong$id)
wstfdf = STFDF(hex, dates, wData)
# Do the same for the count data
flong = reshape(Wind.count, v.names = "count", times = dates, varying = colnames(Wind.count),
ids = rownames(Wind.count), direction = "long")
fData = data.frame(values = flong$count, ID = flong$id)
fstfdf = STFDF(hex, dates, fData, endTime = as.POSIXct(dates))
Create a multi-panel plot of per hexagon max TC wind speed for each year. NOTE, this may take several minutes to render
cr = brewer.pal(7, "Reds") # Max Wind Speed
rng = seq(17, 87, 10) # Max Wind Speed
stplot(wstfdf, col.regions = cr, sp.layout = list(l1), at = rng, colorkey = list(space = "bottom",
labels = paste(rng)), col = "white", names.attr = 1979:2012, sub = expression("Per grid max wind speed (m/s)"))
Create a multi-panel plot of per hexagon TC count for each year. NOTE, this may take several minutes to render
blue.pal = colorRampPalette(c(rgb(158, 202, 225, maxColorValue = 255), rgb(66,
146, 198, maxColorValue = 255), rgb(33, 113, 181, maxColorValue = 255),
rgb(8, 69, 148, maxColorValue = 255)))
rng2 = seq(1, 9, 2) # Count
stplot(fstfdf, col.regions = blue.pal(4), sp.layout = list(l1), at = rng2, colorkey = list(space = "bottom",
labels = paste(rng2)), col = "white", names.attr = 1979:2012, sub = expression("Per grid TC frequency"))
FOR PRELIMINARY EXPLORATORY PURPOSES ONLY! First, divide the hexagons into three subregions: Gulf of Mexico/Carribbean (roughly), Northern basin, and central basin.
GM = c(35, 36, 37, 25, 26, 27, 14, 15, 16, 6)
N = c(61, 62, 63, 49, 50, 51, 52, 53, 38, 39, 40, 41)
C = c(28, 29, 30, 31, 32, 17, 18, 19, 20, 21, 7, 8, 9)
explore.max = Wind.max.yr
hID = lapply(Wind.split, function(x) unique(x$hexid))
hID = as.vector(unlist(hID))
bad = which(hID == 48)
bad2 = which(hID == 47)
hID = hID[-bad]
hID = hID[-bad2]
explore.max$hID = hID
GoM = subset(explore.max, hID %in% GM)
GoM$SR = factor(rep("GoM", dim(GoM)[1]))
NAt = subset(explore.max, hID %in% N)
NAt$SR = factor(rep("NAt", dim(NAt)[1]))
CAt = subset(explore.max, hID %in% C)
CAt$SR = factor(rep("CAt", dim(CAt)[1]))
explore.max = rbind(GoM, NAt, CAt)
Create a plot depicting these three subregions
explore.max.pdf = SpatialPolygonsDataFrame(hpg[rownames(explore.max)], explore.max)
crq = brewer.pal(3, "Set1")
centers = coordinates(explore.max.pdf)
l2 = list("sp.text", centers, explore.max.pdf$hID, col = "yellow", cex = 1.1)
spplot(explore.max.pdf, "SR", col = "white", sp.layout = list(l1, l2), col.regions = brewer.pal(3,
"Set1"), colorkey = list(space = "bottom"), sub = expression("Subregions"))
Create space-time data frames for each subregion
# Create a long table of data...with space varying faster?
hexG = hpg[rownames(GoM)]
GoM = GoM[, 1:34]
GoMlong = reshape(GoM, v.names = "Wmax", times = dates, varying = colnames(GoM),
ids = rownames(GoM), direction = "long")
GoMData = data.frame(values = GoMlong$Wmax, ID = GoMlong$id)
GoMstfdf = STFDF(hexG, dates, GoMData, endTime = as.POSIXct(dates))
#### Do the same for the Northern potion of the basin
hexN = hpg[rownames(NAt)]
NAt = NAt[, 1:34]
NAtlong = reshape(NAt, v.names = "Wmax", times = dates, varying = colnames(NAt),
ids = rownames(NAt), direction = "long")
NAtData = data.frame(values = NAtlong$Wmax, ID = NAtlong$id)
NAtstfdf = STFDF(hexN, dates, NAtData, endTime = as.POSIXct(dates))
#### Do the same for the Central potion of the basin
hexC = hpg[rownames(CAt)]
CAt = CAt[, 1:34]
CAtlong = reshape(CAt, v.names = "Wmax", times = dates, varying = colnames(CAt),
ids = rownames(CAt), direction = "long")
CAtData = data.frame(values = CAtlong$Wmax, ID = CAtlong$id)
CAtstfdf = STFDF(hexC, dates, CAtData, endTime = as.POSIXct(dates))
Generate Hovmoller type plots for per hexagon maximum wind speed. Hexagon ID is given on the x axis, while time is shown vertically. Show one hovmoller for each of the three subregions
crG = brewer.pal(7, "Reds")
rng = seq(17, 87, 10)
crN = brewer.pal(7, "Blues")
crC = brewer.pal(7, "Greens")
p1 = stplot(GoMstfdf, mode = "xt", col.regions = crG, at = rng, xlab = NULL,
colorkey = list(space = "bottom"), main = expression("Gulf of Mexico"),
sub = expression("Per grid max wind speed"))
p2 = stplot(NAtstfdf, mode = "xt", col.regions = crN, at = rng, xlab = NULL,
ylab = NULL, colorkey = list(space = "bottom"), main = expression("Northern Atlantic"),
sub = expression("Per grid max wind speed (m/s)"))
p3 = stplot(CAtstfdf, mode = "xt", col.regions = crC, at = rng, xlab = NULL,
ylab = NULL, colorkey = list(space = "bottom"), main = expression("Central Altantic"),
sub = expression("Per grid max Wind speed"))
plot(p1, split = c(1, 1, 3, 1), more = TRUE)
plot(p2, split = c(2, 1, 3, 1), more = TRUE)
plot(p3, split = c(3, 1, 3, 1), more = FALSE)
Now that the TC data are formatted, we shift our focus to the environmental data (model covariates). For the intensity model, I'd like to use ASO SST, ASO 800-250 hPa vertical shear, ASO average 400-700 hPa relative humidity, the southern oscillation index (SOI), and the difference between ASO SST and ASO 100 mb temperature.
For the count model, I will be using ASO SST, MJ NAO, ASO SOI, and S sunspots (as in Hodges et al. 2014).
Load covariate data
load("C:/Users/Sarah/Dropbox/SpatialExtremes/Sarah/Space-time/mSST.RData")
load("C:/Users/Sarah/Dropbox/SpatialExtremes/Sarah/Space-time/ttData.RData")
load("C:/Users/Sarah/Dropbox/SpatialExtremes/Sarah/Space-time/rhData.RData")
load("C:/Users/Sarah/Dropbox/SpatialExtremes/Sarah/Space-time/vwsData.RData")
Format the SST data similarly to the TC data, with each row representing a specific region/hexagon, and each column representing a specific year
ASOsst = mSST[, 69:170] # ASO
ASOsst = cbind(mSST[, 205:206], ASOsst)
coordinates(ASOsst) = c("lon", "lat")
proj4string(ASOsst) = CRS(ll)
ASOsst.sdf = spTransform(ASOsst, CRS(lcc))
ASOsst.pdf = over(x = hpg[row.names(Wind.count)], y = ASOsst.sdf, fn = mean)
SSTYearMonth = strsplit(substring(colnames(ASOsst.pdf), 2), "M")
SSTYear = as.numeric(sapply(SSTYearMonth, function(x) x[1]))
ASOsst.mean.year = sapply(unique(SSTYear), function(x) as.vector(rowMeans(ASOsst.pdf[,
which(x == SSTYear)])))
dimnames(ASOsst.mean.year) = list(id = rownames(ASOsst.pdf), Year = paste("Y",
unique(SSTYear), sep = ""))
SST.ASO.anom = data.frame(t(sapply(1:nrow(ASOsst.mean.year), function(i) ASOsst.mean.year[i,
] - mean(ASOsst.mean.year[i, ]))))
dimnames(SST.ASO.anom) = list(id = rownames(ASOsst.mean.year), Year = paste("Y",
unique(SSTYear), sep = ""))
Format the 850-200 hPa vertical wind shear (vws) data similarly to the TC and SST data, with each row representing a specific region/hexagon, and each column representing a specific year Note: hexagons 25 and 35 contain grid points with missing values. Hexagon 25 is only missing 1 grid point of 38), while hexagon 35 is missing 16 points (of 42).
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[row.names(Wind.count)], y = vwsData.sdf, fn = mean)
vwsYearMonth = strsplit(substring(colnames(vws.pdf), 5), "M")
vwsYear = as.numeric(sapply(vwsYearMonth, function(x) x[1]))
vws.mean.year = sapply(unique(vwsYear), function(x) as.vector(rowMeans(vws.pdf[,
which(x == vwsYear)])))
dimnames(vws.mean.year) = list(id = rownames(vws.pdf), Year = paste("Y", unique(vwsYear),
sep = ""))
VWS.ASO.anom = data.frame(t(sapply(1:nrow(vws.mean.year), function(i) vws.mean.year[i,
] - mean(vws.mean.year[i, ]))))
dimnames(VWS.ASO.anom) = list(id = rownames(vws.mean.year), Year = paste("Y",
unique(vwsYear), sep = ""))
Format the average 400-700 hPa relative humidity (rh) data similarly to the TC and SST data, with each row representing a specific region/hexagon, and each column representing a specific year Note: No missing values for the relative humidity data, oddly.
coordinates(rhData) = c("lon", "lat")
proj4string(rhData) = CRS(ll)
rhData.sdf = spTransform(rhData, CRS(lcc))
rh.pdf = over(x = hpg[row.names(Wind.count)], y = rhData.sdf, fn = mean)
rhYearMonth = strsplit(substring(colnames(rh.pdf), 4), "M")
rhYear = as.numeric(sapply(rhYearMonth, function(x) x[1]))
rh.mean.year = sapply(unique(rhYear), function(x) as.vector(rowMeans(rh.pdf[,
which(x == rhYear)])))
dimnames(rh.mean.year) = list(id = rownames(rh.pdf), Year = paste("Y", unique(rhYear),
sep = ""))
RH.ASO.anom = data.frame(t(sapply(1:nrow(rh.mean.year), function(i) rh.mean.year[i,
] - mean(rh.mean.year[i, ]))))
dimnames(RH.ASO.anom) = list(id = rownames(rh.mean.year), Year = paste("Y",
unique(rhYear), sep = ""))
RH.ASO.anom = RH.ASO.anom * 100 # as percentage?
Format the average difference between SST and 100 mb temp similarly to the other covariates, with each row representing a specific region/hexagon, and each column representing a specific year Note: No missing values for the relative humidity data, oddly. Also note: I may not be able to use this given the high correlation to SST
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[row.names(Wind.count)], y = ttData.sdf, fn = mean)
ttYearMonth = strsplit(substring(colnames(tt.pdf), 6), "M")
ttYear = as.numeric(sapply(ttYearMonth, function(x) x[1]))
tt.mean.year = sapply(unique(ttYear), function(x) as.vector(rowMeans(tt.pdf[,
which(x == ttYear)])))
dimnames(tt.mean.year) = list(id = rownames(tt.pdf), Year = paste("Y", unique(ttYear),
sep = ""))
sst.mean.year = ASOsst.mean.year + 273.15
TDiff.mean.year = sst.mean.year - tt.mean.year
TDiff.ASO.anom = data.frame(t(sapply(1:nrow(TDiff.mean.year), function(i) TDiff.mean.year[i,
] - mean(TDiff.mean.year[i, ]))))
dimnames(TDiff.ASO.anom) = list(id = rownames(TDiff.mean.year), Year = paste("Y",
unique(ttYear), sep = ""))
Read and format SOI data. We'll use ASO SOI values
soi = read.table("C:/Users/Sarah/Dropbox/SpatialExtremes/Sarah/Space-time/SOI.txt",
header = T)
soi = cbind(soi$YEAR, soi$AUG, soi$SEP, soi$OCT)
ASO = rowMeans(soi[, 2:4])
ASO = as.vector(ASO)
soi = cbind(soi, ASO)
colnames(soi) = c("Year", "AUG", "SEP", "OCT", "ASO")
soi = data.frame(soi)
soi = subset(soi, Year > 1978 & Year < 2013)
As with the TC data, set up space-time data frames for the covariate data
# SST Data
SSTlong = reshape(SST.ASO.anom, v.names = "SST", times = dates, varying = colnames(SST.ASO.anom),
ids = rownames(SST.ASO.anom), direction = "long")
sstData = data.frame(values = SSTlong$SST, ID = SSTlong$id)
SSTstfdf = STFDF(hex, dates, sstData)
# VWS Data
VWSlong = reshape(VWS.ASO.anom, v.names = "VWS", times = dates, varying = colnames(VWS.ASO.anom),
ids = rownames(VWS.ASO.anom), direction = "long")
vwsData = data.frame(values = VWSlong$VWS, ID = VWSlong$id)
VWSstfdf = STFDF(hex, dates, vwsData, endTime = as.POSIXct(dates))
# RH Data
RHlong = reshape(RH.ASO.anom, v.names = "RH", times = dates, varying = colnames(RH.ASO.anom),
ids = rownames(RH.ASO.anom), direction = "long")
rhData = data.frame(values = RHlong$RH, ID = RHlong$id)
RHstfdf = STFDF(hex, dates, rhData, endTime = as.POSIXct(dates))
# TDiff data
TDlong = reshape(TDiff.ASO.anom, v.names = "TDiff", times = dates, varying = colnames(TDiff.ASO.anom),
ids = rownames(TDiff.ASO.anom), direction = "long")
tdData = data.frame(values = TDlong$TDiff, ID = TDlong$id)
TDstfdf = STFDF(hex, dates, tdData, endTime = as.POSIXct(dates))
Create a multi-panel plot of per hexagon ASO SST for each year. NOTE, this may take several minutes to render
cr = rev(brewer.pal(10, "RdBu"))
rng = seq(-2.5, 2.5, 0.5)
stplot(SSTstfdf, col.regions = cr, sp.layout = list(l1), at = rng, colorkey = list(space = "bottom",
labels = paste(rng)), col = "white", names.attr = 1979:2012, sub = expression("Per grid ASO SST anomaly"))
Notice that the pattern is dominated by more negative anomalies early in the period transitioning to more positive anomalies later (trend, decadal/multidecadal variability, both?).
Create a multi-panel plot of per hexagon ASO VWS for each year. NOTE, this may take several minutes to render
cr = rev(brewer.pal(8, "RdBu"))
rng = seq(-8, 8, 2)
stplot(VWSstfdf, col.regions = cr, sp.layout = list(l1), at = rng, colorkey = list(space = "bottom",
labels = paste(rng)), col = "white", names.attr = 1979:2012, sub = expression("Per grid ASO 850-200 hPa vertical shear anomaly"))
I'm still not sure whether ASO average is most suitable; however, for a first try, I'll use it and see what happens. It seems possible that using a more temporally precise VWS (same month as Wmax, or even same or previous day) might be problematic since the presence of the TC should affect the area averaged shear.
Create a multi-panel plot of per hexagon ASO RH for each year. NOTE, this may take several minutes to render
cr = rev(brewer.pal(8, "RdBu"))
rng = seq(-12, 12, 3)
stplot(RHstfdf, col.regions = cr, sp.layout = list(l1), at = rng, colorkey = list(space = "bottom",
labels = paste(rng)), col = "white", names.attr = 1979:2012, sub = expression("Per grid ASO 400-700 hPa RH anomaly (%)"))
Again, I'm not sure whether the ASO average is the most approriate, but I'll give it a try.
Create a multi-panel plot of per hexagon ASO SST - 100hPa T for each year. NOTE, this may take several minutes to render
cr = rev(brewer.pal(8, "RdBu"))
rng = seq(-4, 4, 1)
stplot(TDstfdf, col.regions = cr, sp.layout = list(l1), at = rng, colorkey = list(space = "bottom",
labels = paste(rng)), col = "white", names.attr = 1979:2012, sub = expression("Per grid ASO SST - 100 hPa Temp anomaly"))
As with SST, we see a transition from negative to positive anomalies through time. This pattern is likely accentuated here because of the combined effect of increasing SSTs and decreasing tropopause temps (in MERRA) over this time period (see Emanuel et al. 2013 and Vecchi et al. 2013).
Create side-by-side space-time (Hovmoller) plots for Wmax and SST
crWmax = brewer.pal(7, "Reds")
rngWmax = seq(17, 87, 10)
crSST = rev(brewer.pal(10, "RdBu"))
rngSST = seq(-2.5, 2.5, 0.5)
p4 = stplot(wstfdf, mode = "xt", col.regions = crWmax, at = rngWmax, xlab = NULL,
ylab = NULL, colorkey = list(space = "bottom"), sub = expression("Per grid max Wind speed"))
p5 = stplot(SSTstfdf, mode = "xt", col.regions = crSST, at = rngSST, xlab = NULL,
ylab = NULL, colorkey = list(space = "bottom"), sub = expression("Per grid ASO SST Anomaly"))
plot(p4, split = c(1, 1, 2, 1), more = TRUE)
plot(p5, split = c(2, 1, 2, 1), more = FALSE)
These plots seem like a convenient means of (preliminarily) exploring relationships between per hex Wmax and each of the covariates. The purpose of these is simply initial visualizatiion. Obviously we can't conclude anything on the basis of these plots.
Create side-by-side space-time (Hovmoller) plots for Wmax and wind shear
crVWS = rev(brewer.pal(8, "RdBu"))
rngVWS = seq(-8, 8, 2)
p6 = stplot(wstfdf, mode = "xt", col.regions = crWmax, at = rngWmax, xlab = NULL,
ylab = NULL, colorkey = list(space = "bottom"), sub = expression("Per grid max Wind speed"))
p7 = stplot(VWSstfdf, mode = "xt", col.regions = crVWS, at = rngVWS, xlab = NULL,
ylab = NULL, colorkey = list(space = "bottom"), sub = expression("Per grid ASO 850-200 hPa vertical shear Anomaly"))
plot(p6, split = c(1, 1, 2, 1), more = TRUE)
plot(p7, split = c(2, 1, 2, 1), more = FALSE)
Are higher (lower) ASO wind shear anomalies associated with lower (higher) Wmax? This plot seems to suggest that wind shear may be a useful predictor.
Create side-by-side space-time (Hovmoller) plots for Wmax and RH
crRH = rev(brewer.pal(8, "RdBu"))
rngRH = seq(-12, 12, 3)
p8 = stplot(wstfdf, mode = "xt", col.regions = crWmax, at = rngWmax, xlab = NULL,
ylab = NULL, colorkey = list(space = "bottom"), sub = expression("Per grid max Wind speed"))
p9 = stplot(RHstfdf, mode = "xt", col.regions = crRH, at = rngRH, xlab = NULL,
ylab = NULL, colorkey = list(space = "bottom"), sub = expression("Per grid ASO 400-700 hPa RH Anomaly"))
plot(p8, split = c(1, 1, 2, 1), more = TRUE)
plot(p9, split = c(2, 1, 2, 1), more = FALSE)
Create side-by-side space-time (Hovmoller) plots for Wmax and SST-100 mb T
crTD = rev(brewer.pal(8, "RdBu"))
rngTD = seq(-4, 4, 1)
p10 = stplot(wstfdf, mode = "xt", col.regions = crWmax, at = rngWmax, xlab = NULL,
ylab = NULL, colorkey = list(space = "bottom"), sub = expression("Per grid max Wind speed"))
p11 = stplot(TDstfdf, mode = "xt", col.regions = crTD, at = rngTD, xlab = NULL,
ylab = NULL, colorkey = list(space = "bottom"), sub = expression("Per grid ASO SST-100 hPa Temp Anomaly"))
plot(p10, split = c(1, 1, 2, 1), more = TRUE)
plot(p11, split = c(2, 1, 2, 1), more = FALSE)
First, set up the frequency model similarly to Hodges et al. 2014.
Create a data frame for the per TC max winds that includes hexid for space-time model
Wind.max = do.call(rbind.data.frame, Wind.max)
Wind.max = cbind(Wind.max$name, Wind.max$Yr, Wind.max$Mo, Wind.max$Wmax, Wind.max$hexid)
colnames(Wind.max) = c("name", "Yr", "Mo", "WmaxS", "hexid")
Wind.max = data.frame(Wind.max)
Acquire tropopause temperature data from MERRA files. Note that I save 100, 150, and 200mb as potential tropopause height candidates. Based on what I've examined so far, 100 mb seems to be the best candidate in terms of minimum temperature. This is also what Emanuel Solomon et al. used in their 2013 paper.
source("C:/Users/Sarah/Dropbox/SpatialExtremes/Sarah/Space-time/arrangeMERRA.R")
ttData = get.tt(1979)
for (i in 1980:2012) {
yt = get.tt(i)
ttData = cbind(ttData, yt[, 3:11])
}
save(ttData, file = "C:/Users/Sarah/Dropbox/SpatialExtremes/Sarah/Space-time/ttData.RData")
Calculate monthly 850-200 mb wind shear values from MERRA u and v wind fields.
source("C:/Users/Sarah/Dropbox/SpatialExtremes/Sarah/Space-time/arrangeMERRA.R")
vwsData = get.vws(1979)
for (i in 1980:2012) {
yv = get.vws(i)
vwsData = cbind(vwsData, yv[, 3:5])
}
save(vwsData, file = "C:/Users/Sarah/Dropbox/SpatialExtremes/Sarah/Space-time/vwsData.RData")
source("C:/Users/Sarah/Dropbox/SpatialExtremes/Sarah/Space-time/arrangeMERRA.R")
rhData = get.rh(1979)
for (i in 1980:2012) {
yr = get.rh(i)
rhData = cbind(rhData, yr[, 3:5])
}
save(rhData, file = "C:/Users/Sarah/Dropbox/SpatialExtremes/Sarah/Space-time/rhData.RData")
Load SST data and reformat it (using reformat.SST.temp) to generate a data frame of monthly SSTs (JJASO) for the years 1979–2010. Use relative.SST to generate a data frame of monthly tropical average (30N to 30S) SSTs for the same period. Assign a lambert conformal conic projection to each of the TC and SST data frames
nc = open.ncdf("H:/Data/sst.mnmean.nc")
mSST = reformat.SST.temp(nc, 1979, 2012) ## (monthly SST)
save(mSST, file = "C:/Users/Sarah/Dropbox/SpatialExtremes/Sarah/Space-time/mSST.RData")
rm(nc)