library(RColorBrewer) # Color schemes for maps
library(fields) # Color schemes for maps
# information obtained from: http://iridl.ldeo.columbia.edu/SOURCES/.KAPLAN/.EXTENDED/.v2/
ntime = 55 #Nov-Mar 1901 - Nov-Mar 2018
### Lat - Long grid..
ygrid=seq(-87.5,87.5,by=5)
ny=length(ygrid) # number of latitudinal locations
xgrid=seq(27.5,382.5,by=5)
nx=length(xgrid) # number of longitudinal locations
nglobe = nx*ny
# creation of spatial grid
xygrid=matrix(0,nrow=nx*ny,ncol=2)
i=0
for(iy in 1:ny){
for(ix in 1:nx){
i=i+1
xygrid[i,1]=ygrid[iy]
xygrid[i,2]=xgrid[ix]
}
}
nyrs = 55 #Nov-Mar 1901 - Nov-Mar 2018
data=readBin("Kaplan-SST-DJFM1900-DJFM2018.r4",what="numeric", n=( nx * ny * nyrs), size=4,endian="swap")
data = array(data = data, dim=c( nx, ny, nyrs ) )
data1=data[,,1]
# the lat -long data grid..
index=1:(nx*ny)
index1=index[data1 < 20 & data1 != "NaN"] # only non-missing data.
xygrid1=xygrid[index1,]
nsites=length(index1)
data2=data1[index1]
### SSTdata matrix - rows are seasonal (i.e. one value per year)
## and columns are locations
sstannavg=matrix(NA,nrow=nyrs, ncol=nsites)
for(i in 1:nyrs){
data1=data[,,i]
index1=index[data1 < 20 & data1 != "NaN"]
data2=data1[index1]
sstannavg[i,]=data2
}
indexgrid = index1
rm("data") #remove the object data to clear up space
## write out the grid locations..
write(t(xygrid1),file="kaplan-sst-locs.txt",ncol=2)
N = 55 #Nov-Mar 1964 - Nov-Mar 2013
wuplocs = matrix(scan("Precipitaton_data - Copy.txt"),ncol =4,byrow=T)
#wuplocs
nlocs = dim(wuplocs)[1]
xlats = wuplocs[,2]
xlongs = wuplocs[,1]
#xlongs
nlocs1=nlocs+1
winterprecip = matrix(scan("Max_Winter_Seas_Prec - Copy.txt"),ncol=nlocs1,byrow=T)
years = winterprecip[,1]
winterprecip = winterprecip[,2:nlocs1] #first column is year
nlevel=11
#coltab=rev(brewer.pal(nlevel-1,))
zr=range(winterprecip[N,])
break_a=round(seq(zr[1],zr[2],by=((zr[2]-zr[1])/(nlevel-1))))
quilt.plot(xlongs, xlats,winterprecip[N,],main = "2013 winter precipitation Western-US",zlim=zr)
grid(col="gray70",lty=2)
US(add=TRUE, col="gray50", lwd=2)
box()
sstannavg2=sstannavg[1:N,]
sstannavg1 = scale(sstannavg2)
zs=var(sstannavg1)
#do an Eigen decomposition..
zsvd=svd(zs)
ssteof = zsvd$u
#Principal Components...
pcs= sstannavg1 %*% zsvd$u
#Eigen Values.. - fraction variance
lambdas=(zsvd$d/sum(zsvd$d))
npc =4 # number of PCs to use for CCA
sstPC = pcs[,1:npc]
sstpcfull = pcs
winterprecip1 = scale(winterprecip)
precSd = apply(winterprecip, 2, sd)
precMean = apply(winterprecip, 2, mean)
zs=var(winterprecip1)
#do an Eigen decomposition..
zsvd=svd(zs)
#Eigen vectors
Preceof=zsvd$u
#Principal Components...
precippcs=winterprecip1 %*% zsvd$u
precPC = precippcs[,1:npc]
M=dim(sstPC)[2]
J=dim(precPC)[2]
J=min(M,J)
N = length(precPC[,1])
Qx1 = qr.Q(qr(sstPC))
Qy1 = qr.Q(qr(precPC))
T11 = qr.R(qr(sstPC))
T22 = qr.R(qr(precPC))
VV=t(Qx1) %*% Qy1
# BB = svd(VV)$v
# AA = svd(VV)$u
BB = solve(T22) %*% svd(VV)$v * sqrt(N-1)
wm1 = precPC %*% BB
AA = solve(T11) %*% svd(VV)$u * sqrt(N-1)
vm1 = sstPC %*% AA
M=dim(sstPC)[2]
J=dim(precPC)[2]
J=min(M,J)
N = length(precPC[,1])
Qx1 = qr.Q(qr(sstPC))
Qy1 = qr.Q(qr(precPC))
T11 = qr.R(qr(sstPC))
T22 = qr.R(qr(precPC))
VV=t(Qx1) %*% Qy1
# BB = svd(VV)$v
# AA = svd(VV)$u
BB = solve(T22) %*% svd(VV)$v * sqrt(N-1)
wm1 = precPC %*% BB
AA = solve(T11) %*% svd(VV)$u * sqrt(N-1)
vm1 = sstPC %*% AA
### Predict the winter precipitation PCS
betahat = solve(t(AA) %*% t(sstPC)%*% sstPC %*% AA) %*% t(AA) %*% t(sstPC) %*% precPC
ypred=sstPC %*% AA %*% betahat
### first npc PCs from the PC forecast above and the remaining PCs are set to
## their means - i.e., 0
N1 = dim(winterprecip1)[2]-(npc)
precPCpred = cbind(ypred,matrix(rep(0,N),ncol=N1,nrow=N))
## back transform to get the winter precipitation field
#precpred = rpredcPCpred %*% preceof
#precPCpred = cbind(ypred,rainpcfull[,npc1:dim(precip)[2]])
### Keep only the first npc Eigen Vectors and set rest to zero
E = matrix(0,nrow=dim(winterprecip1)[2],ncol=dim(winterprecip1)[2])
E[,1:npc]=Preceof[,1:npc]
## back transform to get the winter precipitation field
precpred = precPCpred %*% t(E)
### rescale the winter precipitation
precpred=t(t(precpred)*precSd + precMean)
### correlate the forecasted PCs with the historical PCs
pccor = diag(cor(ypred,precPC))
print("the correlation between observed and forecasted PCs is:")
print(round(pccor*1000)/1000)
[1] "the correlation between observed and forecasted PCs is:" [1] 0.177 0.354 0.151 0.215
preccor = diag(cor(precpred,winterprecip))
#plot of R2 between the observed and predicted winter precipitation at each grid point
zr=range(preccor)
quilt.plot(xlongs, xlats,preccor ,nx=30, ny=30,main = bquote(~ R^2 ~ "between the observed and predicted winter precipitation"),zlim=zr)
grid(col="gray70",lty=2)
US(add=TRUE, col="gray50", lwd=2,)
box()
N = length(precPC[,1])
ypred=matrix(0,nrow=N,ncol=J)
index=1:N
for(i in 1:N){
#drop a point..
index1=index[index != i]
xp=precPC[i,] #prediction point
#get the rest of the data..
X=precPC[index1,]
swemeans=colMeans(X)
swesds=apply(X,2,sd)
Y=sstPC[index1,]
flowmeans=colMeans(Y)
flowsds=apply(Y,2,sd)
## scale the data to be predicted..
xp = (xp-swemeans)/swesds
# get the first 4 PCS of SWE
X=scale(X)
Sxx=var(X)
xeof=svd(Sxx)
xpcs=X %*% xeof$u
xpcs=xpcs
X=xpcs
### Perform CCA on the 4 PCS
X1 = scale(X)
Y1=scale(Y)
Qx1 = qr.Q(qr(X1))
Qy1 = qr.Q(qr(Y1))
VV=t(Qx1) %*% Qy1
T11 = qr.R(qr(X1))
T22 = qr.R(qr(Y1))
BB = solve(T22) %*% svd(VV)$v * sqrt(length(Y1[,1]-1))
BB = svd(VV)$v
wm1=Y1 %*% BB
Fyy = var(Y1) %*% BB
AA = solve(T11) %*% svd(VV)$u * sqrt(length(Y1[,1]-1))
AA = svd(VV)$u
vm1 = X1 %*% AA
cancorln = svd(VV)$d[1:J]
#### get the regression coefficient..
betahat = solve(t(AA) %*% t(X1)%*% X1 %*% AA) %*% t(AA) %*% t(X1) %*% Y1
## get the first 4 Pc values of xp
vmp = (xp %*% xeof$u)
vmp = vmp %*% AA
### predict
yp = vmp %*% betahat
### scale back
yp = yp*flowsds + flowmeans
ypred[i,]=yp
}
pccor = diag(cor(ypred,precPC))
print("the correlation between observed and forecasted PCs is:")
print(round(pccor*1000)/1000)
[1] "the correlation between observed and forecasted PCs is:" [1] 0.498 0.754 -0.008 0.072