Assessing uncetainty arising from hexagon configurations

How much uncertainty in the sensitivity of LI to SST can be attributed to the specific spatial configuration of the hexagon grid? I'll test this by generating 500 different hexagon grids, each slightly (and randomly) offset from the others in space, and then calculating the sensitivity of LI to SST for each grid.

Set the working directory

setwd("C:/Users/Sarah/Dropbox/SpatialExtremes")

Acquire the necessary packages

require(ggplot2)
require(rgdal)
require(ismev)
require(ncdf)
library(mapproj)
library(mapdata)
library(maptools)
library(maps)
library(grid)
library(colorRamps)
library(RColorBrewer)
require(spdep)
require(psych)
require(boot)
require(plyr)

Load sensitivityUncertainty.R, which automates the procedure for estimating the sensitivity of LI to SST. Also load R files for estimating the LI

source("selectData.R")
source("ismevplus.R")
source("sensitivityUncertainty.R")

Load observations

load("best.use.RData")
begin = 1981
end = 2010
Wind.df = best.use[best.use$Yr >= begin & best.use$Yr <= end, ]

Wind.df = subset(Wind.df, Type == "*")
Wind.df$WmaxS = Wind.df$WmaxS * 0.5144
quantile(Wind.df$WmaxS, prob = seq(0.6, 0.8, 0.01))  #33 m/s is the 69th percentile

##   60%   61%   62%   63%   64%   65%   66%   67%   68%   69%   70%   71% 
## 28.44 28.87 29.43 29.95 30.48 30.87 31.24 31.80 32.41 32.99 33.34 33.49 
##   72%   73%   74%   75%   76%   77%   78%   79%   80% 
## 33.76 34.29 34.86 35.53 36.01 36.48 37.21 37.90 38.50

Wind.df$maguv = Wind.df$maguv * 0.5144
Wind.df = subset(Wind.df, WmaxS >= 33)
Wind.df$WmaxS = Wind.df$WmaxS - Wind.df$maguv * 0.6
length(unique(Wind.df$Sid))  # 188 tropical cyclones (1981-2010), 268 (1967-2010)

## [1] 188

Load and reorganize SST

nc = open.ncdf("sst.mnmean.nc")
sstvar = nc$var$sst
ncdata = get.var.ncdf(nc, sstvar)
dim(ncdata)

## [1]  180   89 1901

object.size(ncdata)

## 243632368 bytes

vals = lapply(nc$var$sst$dim, function(x) as.vector(x$vals))
vals[[1]] = (vals[[1]] + 180)%%360 - 180
vals[[2]] = rev(vals[[2]])
timedate = as.POSIXlt(86400 * vals[[3]], origin = ISOdatetime(1800, 1, 1, 0, 
    0, 0, tz = "GMT"), tz = "GMT")
timecolumn = paste("Y", 1900 + timedate$year, "M", formatC(as.integer(timedate$mo + 
    1), 1, flag = "0"), sep = "")
names(vals) = sapply(nc$var$sst$dim, "[", "name")
vals = vals[1:2]
ncdata1 = ncdata[, (dim(ncdata)[2]:1), ]
dims = dim(ncdata1)
dim(ncdata1) = c(dims[1] * dims[2], dims[3])
colnames(ncdata1) = timecolumn
ncdata1 = as.data.frame(ncdata1)
ncdataframe = cbind(expand.grid(vals), ncdata1)
misbyrow = apply(ncdataframe, 1, function(x) sum(is.na(x)))
ncdataframe = ncdataframe[misbyrow == 0, ]

ncdataframe = subset(ncdataframe, lon <= 10 & lon >= -100 & lat >= 0 & lat <= 
    70)
ll = "+proj=longlat +datum=WGS84"
lcc = "+proj=lcc +lat_1=60 +lat_2=30 +lon_0=-60"
years = begin:end
sehur = paste("Y", years, sep = "")
sstAug = ncdataframe[paste(sehur, "M", formatC(8, 1, flag = "0"), sep = "")]
sstSep = ncdataframe[paste(sehur, "M", formatC(9, 1, flag = "0"), sep = "")]
sstOct = ncdataframe[paste(sehur, "M", formatC(10), sep = "")]
sstAugmean = rowMeans(sstAug)
sstSepmean = rowMeans(sstSep)
sstOctmean = rowMeans(sstOct)
ASOsst = cbind(sstAugmean, sstSepmean, sstOctmean)
avgsst = rowMeans(ASOsst)
avgsst = as.data.frame(avgsst)
avgsst$lon = ncdataframe$lon
avgsst$lat = ncdataframe$lat
coordinates(avgsst) = c("lon", "lat")
proj4string(avgsst) = CRS(ll)
SST.sdf = spTransform(avgsst, CRS(lcc))

Create spatial data frame from Winds.df

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

Test grid uncertainty

Generate 500 sensitivies for different hexagon grids. Here I estimate LI using the GPD/Poisson model.

allSens = rep(0, 500)
set.seed(7777)
for (i in 1:500) {
    hpts = spsample(Wind.sdf, type = "hexagonal", n = 120, bb = bbox(Wind.sdf) * 
        1.08)
    hpg = HexPoints2SpatialPolygons(hpts)
    allSens[i] = get.sensitivity(Wind.sdf, hpg, SST.sdf)
}

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

Examine range/median/etc. of sensitivity values

load("allSens.RData")
range(allSens)

## [1] -8.524 67.227

median(allSens)

## [1] 5.949

quantile(allSens, 0.025)

##  2.5% 
## 1.537

quantile(allSens, 0.975)

## 97.5% 
## 8.491

allSens = as.data.frame(allSens)
colnames(allSens) = c("Sens")

Ignore low and high values (note: need to come up with a better way to omit grids with xi > 0)

testing = subset(allSens, Sens > 2 & Sens < 9.5)
range(testing$Sens)

## [1] 2.050 9.489

median(testing$Sens)

## [1] 5.993

Plot a histogram of the sensitivities:

p = ggplot(allSens, aes(Sens)) + geom_histogram(fill = "grey", binwidth = 0.25) + 
    theme_bw()
p = p + xlab(expression(paste("Sensitivity [m s", {
}^-1, " K", {
}^-1, "]")))
p = p + ylab("Frequency")
p

Plot a histogram of sensitivities, this time leaving out extremes:

p = ggplot(testing, aes(Sens)) + geom_histogram(fill = "grey", binwidth = 0.25) + 
    theme_bw()
p = p + xlab(expression(paste("Sensitivity [m s", {
}^-1, " K", {
}^-1, "]")))
p = p + ylab("Frequency")
p

Generate 500 sensitivies for different hexagon grids. Here I estimate LI using the GPD/Poisson model, but increase the hexagon area (n = 40)

allSensBIG = rep(0, 500)
set.seed(7777)
for (i in 1:500) {
    hpts = spsample(Wind.sdf, type = "hexagonal", n = 40, bb = bbox(Wind.sdf) * 
        1.08)
    hpg = HexPoints2SpatialPolygons(hpts)
    allSensBIG[i] = get.sensitivity(Wind.sdf, hpg, SST.sdf)
}

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

load("allSensBIG.RData")
range(allSensBIG)

## [1]  3.211 13.909

median(allSensBIG)

## [1] 7.024

quantile(allSensBIG, 0.025)

##  2.5% 
## 4.065

quantile(allSensBIG, 0.975)

## 97.5% 
## 10.62

allSensBIG = as.data.frame(allSensBIG)
colnames(allSensBIG) = c("Sens")

Plot histogram for sensitivities from BIG hexagons

p = ggplot(allSensBIG, aes(Sens)) + geom_histogram(fill = "grey", binwidth = 0.25) + 
    theme_bw()
p = p + xlab(expression(paste("Sensitivity [m s", {
}^-1, " K", {
}^-1, "]")))
p = p + ylab("Frequency")
p

Generate 500 sensitivies for different hexagon grids. Here I use the per hexagon maximum wind speed:

allSens2 = rep(0, 500)
set.seed(7777)
for (i in 1:500) {
    hpts = spsample(Wind.sdf, type = "hexagonal", n = 120, bb = bbox(Wind.sdf) * 
        1.08)
    hpg = HexPoints2SpatialPolygons(hpts)
    allSens2[i] = get.sensitivity2(Wind.sdf, hpg, SST.sdf)
}

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

Examine range/median etc. for allSens2

load("allSens2.RData")
range(allSens2)

## [1] 4.247 8.274

median(allSens2)

## [1] 6.536

quantile(allSens2, 0.025)

##  2.5% 
## 4.863

quantile(allSens2, 0.975)

## 97.5% 
## 7.968

allSens2 = as.data.frame(allSens2)
colnames(allSens2) = c("Sens")

Plot histogram for sensitivities from per hexagon max intensity

p = ggplot(allSens2, aes(Sens)) + geom_histogram(fill = "grey", binwidth = 0.25) + 
    theme_bw()
p = p + xlab(expression(paste("Sensitivity [m s", {
}^-1, " K", {
}^-1, "]")))
p = p + ylab("Frequency")
p