rm(list=ls()) library(SpatialExtremes) ## First load some data (already formatted) data(rainfall) ## Define the trend surfaces loc.form <- scale.form <- y ~ lon * lat shape.form <- y ~ 1 ## Center and scale the covariates covariates <- scale(coord[,1:2]) M0 <- fitspatgev(rain, covariates, loc.form, scale.form, shape.form) ## Define the grid where we want prediction n.grid <- 250 lon.grid <- seq(645, 770, length = n.grid) lat.grid <- seq(209, 291, length = n.grid) grid <- expand.grid(lon.grid, lat.grid) colnames(grid) <- c("lon", "lat") ## <<-- I told you it's important to ## name columns and rows ;-) ## Since we scale the covariates we need to do the same with our grid grid.scaled <- grid for (i in 1:2) grid.scaled[,i] <- (grid.scaled[,i] - attributes(covariates)\$"scaled:center"[i]) / attributes(covariates)\$"scaled:scale"[i] ## Get the prediction for the 25-year return level predictions <- predict(M0, grid.scaled, ret.per = 25) ## Plot it (yes I know I should write a function that do that ## automagically) col.palette <- colorRampPalette(c("blue", "white", "red")) swiss() image(lon.grid, lat.grid, matrix(predictions[,"Q25"], n.grid), col = col.palette(64), add = TRUE) contour(lon.grid, lat.grid, matrix(predictions[,"Q25"], n.grid), add = TRUE) swiss(add = TRUE)