library(evgam)

data(USgusts)

# the following three lines reduce USgusts to only
# those rows corresponding to annual maxima
USgusts.gev <- split(USgusts, USgusts$WMO_id)
USgusts.gev <- lapply(USgusts.gev, function(x) split(x, x$year))
USgusts.gev <- do.call(rbind, lapply(USgusts.gev, function(x) do.call(rbind, lapply(x, function(y) y[which.max(y$gust),]))))

# here we use a thin-plate spline to capture spatial variation in the
# GEV's location, scale and shape parameters and allow the splines
# to have different basis dimensions
fmla.gev <- list(gust ~ s(lon, lat, k=50), ~ s(lon, lat, k=40),  ~ s(lon, lat, k=30))

# fit the model; family="gev" is default
fit.gev <- fevgam(fmla.gev, USgusts.gev)

# this can help identify redundant smooths or model misspecification
summary(fit.gev)

# plot smooths
plot(fit.gev)

# map 100-year return level estimate
library(maps)
state <- map("state")
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)
over.land <- !is.na(map.where("state", df.plot$lon, df.plot$lat))
df.plot <- subset(df.plot, over.land)
return.period <- 100
df.plot <- cbind(df.plot, predict(fit.gev, df.plot, type="quantile", probs=1-1/return.period))

plot.mat <- matrix(NA, length(lon.plot), length(lat.plot))
plot.mat[over.land] <- df.plot[,"q:0.99"]
image(lon.plot, lat.plot, plot.mat, xlab="longitude", ylab="latitude")
map("state", add=TRUE)