library(ggplot2)
library(geoR) #for geodata
library(spBayes) #for splm
library(fields)


library(leaps)
library(MPV)          # for calculating PRESS
library(fields)       # for surface to plot surface plot
library(ggplot2)
library(rgdal)        # for converting projection of coordinates
library(reshape2)     # for using the melt function (Convert an object into a molten data frame)
library(scales) 
library(akima)
library(corrplot)
library(gridExtra)
library(glmnet)
library(ggmap)
library(sm) 
library(locfit)
library(prodlim)
library(maps)         # useful for graphs function map

set.seed(1)

(i) Perform a PCA on the winter global SST anomalies -Plot the Eigen variance spectrum for the first 25 modes -Plot the leading 4 spatial (Eigen vectors) and temporal (PCs) modes of variability (ii) Perform a rotated PCA (rotate the first 6 PCs) and plot the leading 4 spatial and temporal modes of variability. Compare the results with (i) above. (iv) Correlate the first four PCs of rainfall with the SSTs and vice-versa and, show correlation maps.

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)

sstscale = scale(sstannavg)
zs=var(sstscale)

#do an Eigen decomposition..
zsvd=svd(zs)

#Principal Components...
pcs=sstscale %*%zsvd$u
#Eigen Values.. - fraction variance 
lambdas=(zsvd$d/sum(zsvd$d))
#Plot the eigen variance spectrum for the first 25 modes
#Eigen spectrum
ggplot() +
 geom_line(mapping = aes(c(1:25), lambdas[1:25])) +
 geom_point(mapping = aes(c(1:25), lambdas[1:25]), shape = 21, size = 3, color = "gray30", fill ="cadetblue1") +
 labs(title = "Eigen Spectrum",x = "Modes",y = "Frac. Var. explained")+
  theme(plot.title = element_text(size = 13, face = "bold", hjust = 0.5))

sprintf("First PC explained %s percent of the variance",round(sum(lambdas[1])*100))
sprintf("First 2 PC explained %s percent of the variance",round(sum(lambdas[1:2])*100))
sprintf("First 3 PCs explained %s percent of the variance",round(sum(lambdas[1:3])*100))
sprintf("First 4 PCs explained %s percent of the variance",round(sum(lambdas[1:4])*100))

'First PC explained 21 percent of the variance'

'First 2 PC explained 34 percent of the variance'

'First 3 PCs explained 41 percent of the variance'

'First 4 PCs explained 47 percent of the variance'

PCA

xlong = sort(unique(xygrid[,2]))
ylat = sort(unique(xygrid[,1]))

par(mfrow = c(4, 2))
par(mar = c(3, 4, 2, 1))
for(i in 1:4){
  zfull = rep(NaN,nglobe)   #also equal 72*36
  zfull[indexgrid]=zsvd$u[,i]
  zmat = matrix(zfull,nrow=nx,ncol=ny)
  image.plot(xlong,ylat,zmat,ylim=range(-40,70),ann=FALSE)
  mtext(paste("EOF no.",i,sep=""), side = 3, line = 0.2, cex = 0.8)
  contour(xlong,ylat,(zmat),ylim=range(-40,70),add=TRUE,nlev=6,lwd=2)
  map("world2",add=T)
  box()
  
  plot( scale(pcs[,i]),type="l",xlab="Year",axes=FALSE,ann=FALSE)
  axis(2)
  axis(1)
  mtext(paste("PC no.",i,sep=""), side = 3, line = 0.2, cex = 0.8)
}

Rotated PCA

zrot = varimax(zsvd$u[,1:6],normalize=FALSE)
rotpcs=sstannavg %*% zrot$loadings
xlong = sort(unique(xygrid[,2]))
ylat = sort(unique(xygrid[,1]))
par(mfrow = c(4, 2))
par(mar = c(3, 4, 2, 1))

for(i in 1:4){
  zfull = rep(NaN,nglobe)   #also equal 72*36
  zfull[indexgrid]=zrot$loadings[,i]
  zmat = matrix(zfull,nrow=nx,ncol=ny)
  image.plot(xlong,ylat,zmat,ylim=range(-50,70),ann=FALSE)
  mtext(paste("Rotated EOF no.",i,sep=""), side = 3, line = 0.2, cex = 0.8)
  contour(xlong,ylat,(zmat),ylim=range(-50,70),add=TRUE,nlev=6,lwd=2)
  map("world2",add=T)
  box()
  
  plot(scale(rotpcs[,i]),type="l",xlab="Year",axes=FALSE,ann=FALSE)
  axis(2)
  axis(1)
  mtext(paste("Rotated PC no.",i,sep=""), side = 3, line = 0.2, cex = 0.8)
}

N = 49  #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

### Perform PCA

#get variance matrix..
winterprecip1 = scale(winterprecip)


zs=var(winterprecip1)

#do an Eigen decomposition..
zsvd=svd(zs)
#Eigen vectors
Preceof=zsvd$u
#Principal Components...
precippcs=winterprecip1 %*% zsvd$u

#Eigen Values.. - fraction variance 
lambdas=(zsvd$d/sum(zsvd$d))

#Plot the eigen variance spectrum for the first 25 modes
#Eigen spectrum
ggplot() +
 geom_line(mapping = aes(c(1:25), lambdas[1:25])) +
 geom_point(mapping = aes(c(1:25), lambdas[1:25]), shape = 21, size = 3, color = "gray30", fill ="cadetblue1") +
 labs(title = "Eigen Spectrum",x = "Modes",y = "Frac. Var. explained")+
  theme(plot.title = element_text(size = 13, face = "bold", hjust = 0.5))

sprintf("First 3 PCs explained %s percent of the variance",round(sum(lambdas[1:3])*100))
sprintf("First 4 PCs explained %s percent of the variance",round(sum(lambdas[1:4])*100))

'First 3 PCs explained 38 percent of the variance'

'First 4 PCs explained 44 percent of the variance'

par(mfrow = c(4, 2))
par(mar = c(3, 4, 2, 1))
for(i in 1:4){
  quilt.plot(xlongs, xlats, zsvd$u[,i],xlim=range(-114,-103),ann=FALSE)
  mtext(paste("EOF no.",i,sep=""), side = 3, line = 0.2, cex = 0.8)
  grid(col="gray70",lty=2)
  US(add=TRUE, col="gray40", lwd=1,xlim=range(-114,-103))
  box()
  
  plot(1964:2018, precippcs[,i],type="l",axes=FALSE,ann=FALSE)
  axis(2)
  axis(1)
  mtext(paste("PC no.",i,sep=""), side = 3, line = 0.2, cex = 0.8)
}

zrot = varimax(zsvd$u[,1:6],normalize=FALSE)
rotpcs=t(t(zrot$loadings) %*% t(winterprecip1))
## plot the first 4 rotated spatial mode and the rotated PCs
par(mfrow = c(4, 2))
par(mar = c(3, 4, 2, 1))
for(i in 1:4){
  quilt.plot(xlongs, xlats, zrot$loadings[,i],xlim=range(-114,-103),ann=FALSE)
  mtext(paste("Rotated EOF no.",i,sep=""), side = 3, line = 0.2, cex = 0.8)
  grid(col="gray70",lty=2)
  US(add=TRUE, col="gray40", lwd=1,xlim=range(-114,-103))
  box()
  
  plot(1964:2018, scale(rotpcs[,i]),type="l",axes=FALSE,ann=FALSE)
  axis(2)
  axis(1)
  mtext(paste("Rotated PC no.",i,sep=""), side = 3, line = 0.2, cex = 0.8)
}

Pgrid=cbind(xlongs, xlats)

cor_pc=cor(precippcs[,1:4], pcs[,1:4])# correlation between PCS
par(mfrow = c(4, 4))
par( oma=c(2,0,4,0)) # save some room for the legend
par(mar = c(4, 4, 2, 1))
for (i in 1:4) {
  for (j in 1:4) {
    if(i==1 & j==1){
      plot(pcs[,j],precippcs[,i],type="p",xlab=paste("R2=",round(cor_pc[i,j]*1000)/1000,sep = ""),ylab=as.character(i),xlim=c(-20,20),ylim=c(-10,10),main=as.character(j),axes=FALSE,font.lab=2)
      axis(2)
      axis(1)
    }else if(i==1){
      plot(pcs[,j],precippcs[,i],type="p",xlab=paste("R2=",round(cor_pc[i,j]*1000)/1000,sep = ""),ylab="",xlim=c(-20,20),ylim=c(-10,10),main=as.character(j),axes=FALSE,font.lab=2 )
      axis(2)
      axis(1)
    }else if(j==1){
      plot(pcs[,j],precippcs[,i],type="p",xlab=paste("R2=",round(cor_pc[i,j]*1000)/1000,sep = ""),ylab=as.character(i),xlim=c(-20,20),ylim=c(-10,10),axes=FALSE,font.lab=2 )
      axis(2)
      axis(1)
    }
    else {
      plot(pcs[,j],precippcs[,i],type="p",xlab=paste("R2=",round(cor_pc[i,j]*1000)/1000,sep = ""),ylab="",xlim=c(-20,20),ylim=c(-10,10),axes=FALSE,font.lab=2 )
      axis(2)
      axis(1)
    }}}
mtext("Correlation between four PCs of winter SSTs with that of the western US precipitation", outer = TRUE, side = 3, cex = 1.2, line = 1)