library(chron)
library(fields)
library(evgam)

data(USgusts)
USgusts$date <- chron(paste(USgusts$month, USgusts$day, USgusts$year, sep="/"))
USgusts$cyc <- as.numeric(strftime(USgusts$date, "%j")) / 366
hurricane.states <- unique(subset(USgusts, lon > -95 & lat < 40)$state)[-c(10, 11, 18, 19)]
USgusts <- subset(USgusts, state %in% hurricane.states)

# here we use a tensor product between cyclic cubic and 
# thin-plate splines to capture an interaction between 
# within-year and spatial variation in the quantile,
# i.e. the asymmetric Laplace distribution's location
# and scale parameters

fmla.ald <- gust ~ te(lon, lat, cyc, d=c(2, 1), k=c(30, 8), bs=c('tp', 'cc'))
# fit the quantile regression model using 97th percentile
# might take a minute or two
fit.ald <- fevgam(fmla.ald, USgusts, family="ald", tau=.97, outer="nlm", trace=2)

USgusts$threshold <- predict(fit.ald)$location
USgusts$excess <- USgusts$gust - USgusts$threshold

# create a data frame with just threshold exceedances
USgusts.gpd <- subset(USgusts, excess > 0)

fmla.gpd <- excess ~ te(lon, lat, cyc, d=c(2, 1), k=c(30, 8), bs=c('tp', 'cc'))
fit.gpd <- fevgam(fmla.gpd, USgusts.gpd, family="gpd", trace=2)
USgusts.gpd <- cbind(USgusts.gpd, predict(fit.gpd, type="response"))

summary(fit.gpd)

## 100-year return level estimation 

state <- map("state", region=hurricane.states)
lon.plot <- seq(min(state$x, na.rm=TRUE), max(state$x, na.rm=TRUE), l=100)
lat.plot <- seq(min(state$y, na.rm=TRUE), max(state$y, na.rm=TRUE), l=100)
df.plot <- expand.grid(lon=lon.plot, lat=lat.plot)
is.hurricane <- sapply(strsplit(map.where("state", df.plot$lon, df.plot$lat), ":"), function(x) x[1]) %in% hurricane.states
df.plot <- subset(df.plot, is.hurricane)
return.period <- 100 # no. days
df.plot$rl1 <- NA

## without extremal index (not advisable)

for (i in seq_len(nrow(df.plot))) {
df.i <- as.list(df.plot[i, c("lon", "lat")])
df.i$cyc <- ppoints(20)
df.i <- expand.grid(df.i)
df.i <- cbind(df.i, predict(fit.gpd, df.i, type="response"))
df.i$threshold <- predict(fit.ald, df.i, type="response")$location
df.plot$rl1[i] <- qblockmax(df.i$threshold, df.i$scale, df.i$shape, fit.ald$tau, 365.25, 1, 1, p=1-1/return.period)
}

plot.mat1 <- matrix(NA, length(lon.plot), length(lat.plot))
plot.mat1[is.hurricane] <- df.plot$rl1
image.plot(lon.plot, lat.plot, plot.mat1, xlab="longitude", ylab="latitude")
lines(state$x, state$y)

## with extremal index

# the following does a probability integral transform to convert to
# censored unit Frechet margins
USgusts$exi <- USgusts$excess > 0
USgusts$frech <- -1/log(fit.ald$tau)
USgusts$frech[USgusts$exi] <- -1/log(1 - (1 - fit.ald$tau) * (1 - pgpd(USgusts.gpd$excess, 0, USgusts.gpd$scale, USgusts.gpd$shape)))

# this creates an integer time variable
# which is needed for calculating consecutive days'
# maxima
USgusts$chron <- as.numeric(chron(paste(USgusts$year, USgusts$month, USgusts$day, sep=""), format="ymd"))

USgusts.exi <- split(USgusts, USgusts$WMO_id)

USgusts.exi <- do.call(rbind, lapply(USgusts.exi, give2DayMax, cons="chron", ynm="frech"))

# estimate extremal index with 
# thin-plate spline for spatial variation
# and probit link
# smooth parameter estimate tends not to converge too well
# when there's only one smoothing parameter
fit.exi <- fevgam(frech ~  te(lon, lat, cyc, d=c(2, 1), k=c(30, 8), bs=c('tp', 'cc')), USgusts.exi, family="exi", exilink="probit", outer="nlm", trace=2)
# plot estimated spline for extremal index (excludes constant)

for (i in seq_len(nrow(df.plot))) {
df.i <- as.list(df.plot[i, c("lon", "lat")])
df.i$cyc <- ppoints(20)
df.i <- expand.grid(df.i)
df.i <- cbind(df.i, predict(fit.gpd, df.i, type="response"))
df.i$threshold <- predict(fit.ald, df.i, type="response")$location
df.i$exi <- predict(fit.exi, df.i, type="response")$location
df.plot$rl2[i] <- qblockmax(df.i$threshold, df.i$scale, df.i$shape, fit.ald$tau, 365.25, 1, df.i$exi, p=1-1/return.period)
}

plot.mat2 <- matrix(NA, length(lon.plot), length(lat.plot))
plot.mat2[is.hurricane] <- df.plot$rl2
image.plot(lon.plot, lat.plot, plot.mat2, xlab="longitude", ylab="latitude")
lines(state$x, state$y)

# compare 100-year return level estimates with and without extremal index
image.plot(lon.plot, lat.plot, plot.mat2 / plot.mat1, xlab="longitude", ylab="latitude")
lines(state$x, state$y)