Read data

time_elapsed=0
mainDir=getwd()
knitr::opts_knit$set(root.dir = mainDir)
setwd(mainDir)
#load the libary file
suppressPackageStartupMessages(source("lib.R"))
source("Lib.R")

#Max summer precipitation
data=read.delim("Max_Prec.txt", header = TRUE, sep = "\t", dec = ".")#in mmm
#Total summer precipitation
data1=read.delim("Total_Prec.txt", header = TRUE, sep = "\t", dec = ".")#in mmm
#Stations information (lat,long, elev, etc...)
info=read.delim("station_location_swus.txt", header = TRUE, sep = "\t", dec = ".")
#load covariates
#MEI values
ENSO=read.delim("meiv.txt", header = TRUE, sep = "\t", dec = ".")
ENSO_S=cbind(ENSO[(ENSO$YEAR>=min(data$YEAR) & ENSO$YEAR<=max(data$YEAR)),1],apply(ENSO[(ENSO$YEAR>=min(data$YEAR) & ENSO$YEAR<=max(data$YEAR)),7:11], 1, FUN = mean))
colnames(ENSO_S)=c("year","val")
#PDO values
PDO=read.delim("PDO.txt", header = TRUE, sep = "\t", dec = ".")
PDO_S=cbind(PDO[(PDO$YEAR>=min(data$YEAR) & PDO$YEAR<=max(data$YEAR)),1],apply(PDO[(PDO$YEAR>=min(data$YEAR) & PDO$YEAR<=max(data$YEAR)),7:10], 1, FUN = mean))
colnames(PDO_S)=c("year","val")

elev_grid=read.delim("Elevation_grid_1deg.txt", header = TRUE, sep = "\t", dec = ".")

getting the standarized total precip and the covariates data frame

Pmean=apply(data1[,2:ncol(data1)],1,mean)
Pmean=(Pmean-mean(Pmean))/sd(Pmean)#nomalized Pemean
datafit=as.data.frame(cbind(ENSO_S,PDO_S[,2],Pmean))
colnames(datafit)=c("t","ENSO","PDO","PT")

Fit a non stationary GEV for each location considering only the location parameter as nonstationary

\(Y(s_i,t)\sim GEV(\mu(s_i,t),\sigma(s_i),\xi(s_i))\)

where

\(\mu(s_i,t)=\alpha_{\mu_0}(s_i)+\alpha_{\mu_1}ENSO(t)+\alpha_{\mu_2}PDO(t)+\alpha_{\mu_2}PT(t)\)

\(log(\sigma(s_i))=\alpha_{\sigma_0}(s_i)\)

ptm = proc.time()
name_par=c("Station","mu_0","mu_1","mu_2","mu_3","sigma","xi")
params=data.frame(matrix(data=-NA,nrow = dim(data)[2]-1,ncol=length(name_par)))#to save the GEV parameter of each station
params[,1]=info$Station
colnames(params)=name_par
gev_ob=list(desc="")
for (i in 2:dim(data)[2]) {
  # print(i-1)
  
  fit.gev= fevd(data[,i],datafit,location.fun=~ENSO+PDO+PT,use.phi = TRUE)
  gev_ob[[paste("Est_",i-1,sep = "")]]=fit.gev
  params[i-1,2:dim(params)[2]]=fit.gev$results$par
  # cov_mat[[paste("st_",i-1,sep = "")]]=parcov.fevd(fit.gev)#solve(fit.gev$results$hessian)
}
time_elapsed=time_elapsed+(proc.time() - ptm)[3]
ptm = proc.time()
X=info[,c(4,3,5)]# long, lat, and elev
X1=elev_grid #long, lat, and elev of 1deg grid
n.samples=1000#the number of MCMC iterations
burn=0.5#numner of samples discarded, warmup
thin1=1#sample thinning factor, if thin = 10 then 1 in 10 samples are considered between start and end
displ=TRUE
thin1=1
plots=list(desc="plots")#to save posterior distribution plots
par=list(desc="parameter values")#to save the postrior vaues of the GEV coefficients over the grid
acep_rate=c()#to save the aceptance_rate for each coefficient simulated
for (i in 2:dim(params)[2]) {
  Y=params[,i]
  result=get.samples_1(Y,X,X1,colnames(params)[i],n.samples,burn,thin1,displ=FALSE)
  plots[[paste("Dist",(i-1),sep = "")]]=result$p
  par[[paste("Median_par",(i-1),sep = "")]]=result$Pred
  par[[paste("par",(i-1),sep = "")]]=result$y
  acep_rate=c(acep_rate,result$acep_r[2])
}
time_elapsed=time_elapsed+(proc.time() - ptm)[3]
sprintf("Aceptance rate for each coefficients are (percent): %s ",acep_rate)
## [1] "Aceptance rate for each coefficients are (percent): 35.6 "
## [2] "Aceptance rate for each coefficients are (percent): 37.4 "
## [3] "Aceptance rate for each coefficients are (percent): 37.2 "
## [4] "Aceptance rate for each coefficients are (percent): 39.2 "
## [5] "Aceptance rate for each coefficients are (percent): 47 "  
## [6] "Aceptance rate for each coefficients are (percent): 46.2 "
sprintf("elapsed time: %s minutes",round((proc.time() - ptm)[3]/60,2))
## [1] "elapsed time: 0.09 minutes"

Histogram of the GEV coeficients

MLE

name=colnames(params)[2:dim(params)[2]]
p1=ggplot() +
  geom_histogram(aes(as.vector(params[,2]),..density..),fill="gray50",col='black') +
  labs(title = TeX(paste('Histogram of $\\',name[1],'$',sep = "")),
       x = "Predictions")+
  #theme_light()+
  theme(plot.title = element_text(size = 13, face = "bold", hjust = 0.5))

p2=ggplot() +
  geom_histogram(aes(as.vector(params[,3]),..density..),fill="gray50",col='black') +
  labs(title = TeX(paste('Histogram of $\\',name[2],'$',sep = "")),
       x = "Predictions")+
  #theme_light()+
  theme(plot.title = element_text(size = 13, face = "bold", hjust = 0.5))

p3=ggplot() +
  geom_histogram(aes(as.vector(params[,4]),..density..),fill="gray50",col='black') +
  labs(title = TeX(paste('Histogram of $\\',name[3],'$',sep = "")),
       x = "Predictions")+
  #theme_light()+
  theme(plot.title = element_text(size = 13, face = "bold", hjust = 0.5))

p4=ggplot() +
  geom_histogram(aes(as.vector(params[,5]),..density..),fill="gray50",col='black') +
  labs(title = TeX(paste('Histogram of $\\',name[4],'$',sep = "")),
       x = "Predictions")+
  #theme_light()+
  theme(plot.title = element_text(size = 13, face = "bold", hjust = 0.5))

p5=ggplot() +
  geom_histogram(aes(as.vector(params[,6]),..density..),fill="gray50",col='black') +
  labs(title = TeX(paste('Histogram of $\\',name[5],'$',sep = "")),
       x = "Predictions")+
# theme_light()+
theme(plot.title = element_text(size = 13, face = "bold", hjust = 0.5))

p6=ggplot() +
  geom_histogram(aes(as.vector(params[,7]),..density..),fill="gray50",col='black') +
  labs(title = TeX(paste('Histogram of $\\',name[6],'$',sep = "")),
       x = "Predictions")+
# theme_light()+
theme(plot.title = element_text(size = 13, face = "bold", hjust = 0.5))


multiplot(p1,p2,p3,p4,p5,p6, cols = 3)

** Spatial Bayesian model**

multiplot(plots$Dist1,plots$Dist2,plots$Dist3,plots$Dist4,plots$Dist5,plots$Dist6, cols = 3)

Spatial maps of posterior parameters, GEV fitted parameters, and standard error

to get the locations of the grid cells in a rectangular grid

ygrid=seq(31.5,42,by=1)
# ygrid_1=seq(range(ygrid)[1],range(ygrid)[2],by=res_grid)
ny=length(ygrid) # number of latitudinal locations

xgrid=seq(-114.5,-102,by=1)
# xgrid_1=seq(range(xgrid)[1],range(xgrid)[2],by=res_grid)
nx=length(xgrid) # number of longitudinal locations

nglobe = nx*ny
# creation of spatial grid
xygrid=matrix(0,nrow=nglobe,ncol=2)
i=0
for(iy in 1:ny){
  for(ix in 1:nx){
    i=i+1
    xygrid[i,2]=ygrid[iy]
    xygrid[i,1]=xgrid[ix]
  }
}

#find the locations in the rectangular grid
index1=c()
for (i in 1:nrow(X1)) {
  tmp=which(X1[i,1]==xygrid[,1] & X1[i,2]==xygrid[,2])
  index1=c(index1,tmp)
}

Spatial Bayesian model 

r_pal=rev(brewer.pal(10,"RdBu"))
MainStates <- map_data("state")
xlim1=c(-115,-101.3)
ylim1=c(29.5,42.5)
namex=seq(from=-114,to=-102,by=2)
labx=paste(abs(namex[1]),"\u00B0W",sep = "")
for (i in 2:length(namex)) {
  labx=c(labx,paste(abs(namex[i]),"\u00B0W",sep = ""))
}
namey=seq(from=30,to=42,by=2)
laby=paste(abs(namey[1]),"\u00B0N",sep = "")
for (i in 2:length(namey)) {
  laby=c(laby,paste(abs(namey[i]),"\u00B0N",sep = ""))
}

source("lib.R")
for (i in 2:dim(params)[2]) {
  Y=params[,i]
  plot_sp_grid(par[[paste("Median_par",(i-1),sep = "")]],Y,X[,1:2],xgrid,ygrid,nx*ny,index1,colnames(params)[i],xlim1,ylim1,namex,namey,labx,laby,r_pal)
}

Estimate 2-year and 100-year return level of max seasonal precipitation for a wet year

finding the wet year based on the average spatial 3-day precipitation

ind_pc_max=which(datafit$PT==max(datafit$PT))#pick a wet year
sprintf("The wet year is: %s",data$YEAR[ind_pc_max])
## [1] "The wet year is: 1984"

get Simulations for the wet year

ptm = proc.time()
period=c(2,100)# return periods considered
set.seed(3000)
# xygrid=elev_grid$elev_data[,1:2]
nsim=ncol(par$par1)#number of samples
estimate=list(T2=matrix(NA,nrow(par$par1),ncol(par$par1)),T100=matrix(NA,nrow(par$par1),ncol(par$par1)))

for (sim in 1:nsim) {
  par1=data.frame(matrix(data=NA,nrow = dim(X1)[1],ncol=dim(params)[2]-1))#to save the GEV coefficients of the grid
  for (i in 1:(dim(params)[2]-1)) {
    par1[,i]=par[[paste("par",(i),sep = "")]][,sim]
  }
  
  for (i in 1:nrow(X1)) {
    loc=as.matrix(datafit[ind_pc_max,2:4]) %*% as.numeric(par1[i,2:4])+as.numeric(par1[i,1])
    loc=as.numeric(loc)
    scal=exp(as.numeric(par1[i,5]))
    shap=as.numeric(par1[i,6])
    tmp1=as.numeric(qgev(1-1/period, xi = shap, mu = loc, beta = scal, lower.tail = TRUE))#column 1: 2yr, column 2:100yr
    estimate$T2[i,sim]=tmp1[1]
    estimate$T100[i,sim]=tmp1[2]
    
  }
}
time_elapsed=time_elapsed+(proc.time() - ptm)[3]
sprintf("elapsed time: %s minutes",round((proc.time() - ptm)[3]/60,2))
## [1] "elapsed time: 0.75 minutes"

plot of the median for T=2yr and T=100yr

med_predT2=apply(estimate$T2, 1, median)
med_predT100=apply(estimate$T100, 1, median)
Y=as.numeric(data[ind_pc_max,2:ncol(data)])
par( oma=c(0,0,3,1)) # save some room for the legend
par(mfrow = c(1, 3))
par(mar = c(4.5, 4.7, 2, 1.5))

# plot of the observed precip for the wet year
quilt.plot(X[,1],X[,2],Y,axes=F,xlab="Longitude",ylab="Latitude",main="Observed",ylim=ylim1,xlim=xlim1,col=r_pal,cex.main=1.4,cex.lab=1.3)
US(add=TRUE, col="gray30", lwd=2)
box()
axis(1, at=namex, labels=labx)
axis(2, at=namey, labels=laby)

zfull = rep(NaN,nglobe)
zfull[index1]=med_predT2
zmat = matrix(zfull,nrow=length(xgrid),ncol=length(ygrid))
image.plot(xgrid,ygrid,zmat,ylim=ylim1,xlim=xlim1,main ="Median T=2 yrs estimated using spatial Bayesian model",axes=F,xlab="Longitude",ylab="Latitude",cex.main=1.4,cex.lab=1.3,col=r_pal)
US(add=TRUE, col="gray30", lwd=2)
box()
axis(1, at=namex, labels=labx,cex.axis=1)
axis(2, at=namey, labels=laby,cex.axis=1)

zfull = rep(NaN,nglobe)
zfull[index1]=med_predT100
zmat = matrix(zfull,nrow=length(xgrid),ncol=length(ygrid))
image.plot(xgrid,ygrid,zmat,ylim=ylim1,xlim=xlim1,main ="Median T=100 yrs estimated using spatial Bayesian model",axes=F,xlab="Longitude",ylab="Latitude",cex.main=1.4,cex.lab=1.3,col=r_pal)
# contour(xgrid,ygrid,zmat,ylim=ylim1,xlim=xlim1,add=TRUE,nlev=4,lwd=1)
US(add=TRUE, col="gray30", lwd=2)
# grid(col="gray30",lty=2)
box()
axis(1, at=namex, labels=labx,cex.axis=1)
axis(2, at=namey, labels=laby,cex.axis=1)

mtext(paste("Summer Season 3-day Maximum Precipitation for the Wettest Year (",data$YEAR[ind_pc_max],"), (mm)",sep = ""), outer = TRUE, side = 3, cex = 1.6, line = 1)