##Define the coordinate of each location
n.site <- 30
locations <- matrix(runif(2*n.site, 0, 10), ncol = 2)
colnames(locations) <- c("lon", "lat")
##Simulate a max-stable process - with unit Frechet margins
data <- rmaxstab(40, locations, cov.mod = "whitmat", sill = 1, range = 30,
smooth = 0.5)
##Now define the spatial model for the GEV parameters
param.loc <- -10 - 4 * locations[,1] + locations[,2]^2
param.scale <- 5 + locations[,1] + locations[,2]^2 / 10
param.shape <- rep(.2, n.site)
##Transform the unit Frechet margins to GEV
for (i in 1:n.site)
data[,i] <- frech2gev(data[,i], param.loc[i], param.scale[i],
param.shape[i])
##Define a model for the GEV margins to be fitted
##shape ~ 1 stands for the GEV shape parameter is constant
##over the region
loc.form <- loc ~ lon + I(lat^2)
scale.form <- scale ~ lon + I(lat^2)
shape.form <- shape ~ lat + lon
## 1- Fit a max-stable process
schlather <- fitmaxstab(data, locations, "whitmat", loc.form, scale.form,
shape.form)
## 2- Produce a map of the location, scale and shape parameters
map(schlather, "loc", col = rainbow(80))
title("Location")
map(schlather, "scale", col = heat.colors(80))
title("Scale")
map(schlather, "shape", col = topo.colors(100))
title("Shape")
## 3- Produce a map for the 50 years return level
new.ranges <- cbind(c(0, 15), c(0, 15))
colnames(new.ranges) <- c("lon", "lat")
map(schlather, "quant", ret.per = 50 , ranges = new.ranges)
title("50-year return level")Run the code above in your browser using DataLab