#### rotation Example
### Obtain summer (Jun-Sep, JJAS) average rainfall over India
#### and average sea surface temperature (SST) over the tropics 25N - 25S


###################### REAd in RAjeevan JJAS AISMR

### Read in Rajeevan Data (summer JJAS average) 1951 - 2009 and perform PCA

nrows=35
ncols=33
ntime = 59    #JJAS 1951 - JJAS 2009
ntimep = 59   # JJAS 1951 - JJAS 2000
N = nrows*ncols


### Lat - Long grid..

ygrid=seq(6.5,38.5,by=1)
ny=length(ygrid)

xgrid=seq(66.5,100.5,by=1)
nx=length(xgrid)

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]
}

}

# REad Rajeevan seasonal average data..

data=readBin("Rajeevan-AISMR-JJAS1951-JJAS2009.r4",what="numeric", n=( nrows * ncols * ntime), size=4,endian="swap")


data <- array(data = data, dim=c( nrows, ncols, ntime ) )


data1=data[,,1]
# Missing value is NaN, change it to a negative number.
data1[data1 == "NaN"]=-99999.
# the lat -long data grid..

index=1:(nx*ny)

index1=index[data1 >= 0]	# only non-missing data.
xygrid1=xygrid[index1,]
x1=xygrid1[,2]

nsites=length(index1)	# locations with data -i.e. global locations
data2=data1[index1]

### Rain data matrix - rows are years and columns are locations on the grid over India
raindata=matrix(NA,nrow=ntimep, ncol=nsites)


for(i in 1:ntimep){
data1=data[,,i]
data1[data1 == "NaN"]=-99999.
index1=index[data1 >= 0]
data2=data1[index1]
raindata[i,]=data2
}

indexgrid = index1

rm("data")	#remove the object data to clear up space

## set the data grid and index
raingrid = xygrid1
xygridrain = xygrid
indexrain = index1

################# Now perform the PCA on the summer rainfall

#get variance matrix..
raindata1=scale(raindata)
zs=var(raindata1)

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

#Principal Components...
pcs=raindata1 %*% zsvd$u
rainpcs = pcs

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

### Plot the first 25 modes
plot(1:25, lambdas[1:25], type="l", xlab="Modes", ylab="Frac. Var. explained")
points(1:25, lambdas[1:25], col="red")

e first four spatial component or Eigen Vector pattern.. and the
## time seris of the first four PCs

library(maps)
library(akima)
library(fields)


# the data is on a grid so fill the entire globaal grid with NaN and then populate the ocean grids with
# the Eigen vector
xlong = sort(unique(xygridrain[,2]))
ylat = sort(unique(xygridrain[,1]))
nrows=35
ncols=33
nglobe = nrows*ncols   # also equal to 35*33

zfull = rep(NaN,nglobe)
zfull[indexrain]=zsvd$u[,3]      ###### replace this with zsvd$u[,2] for second Eigen Vector etc.
zmat = matrix(zfull,nrow=nrows,ncol=ncols)

image.plot(xlong,ylat,zmat,ylim=range(2,40),zlim=range(-0.2,0.2))
contour(xlong,ylat,(zmat),ylim=range(2,40),add=TRUE,nlev=6,lwd=2)
map('world2',add=TRUE)
grid(col="black",lty=1)

## Plot PC1. Similarly, plot PCs 2 through 4
plot(1951:2009, pcs[,1],xlab="Year",ylab="PC1",type="b")
title(main="PC1 of summer rainfall over India")

###  Rotate first six PCS
##################
zrot = varimax(zsvd$u[,1:6],normalize=FALSE)
zrot = promax(zsvd$u[,1:6],m=2)
xlong = sort(unique(xygridrain[,2]))
ylat = sort(unique(xygridrain[,1]))
nrows=35
ncols=33
nglobe = nrows*ncols   # also equal to 35*33

zfull = rep(NaN,nglobe)
zfull[indexrain]=zrot$loadings[,3]      ###### replace this with zsvd$u[,2] for second Eigen Vector etc.
zmat = matrix(zfull,nrow=nrows,ncol=ncols)

image.plot(xlong,ylat,zmat,ylim=range(2,40),zlim=range(-0.2,0.2))
contour(xlong,ylat,(zmat),ylim=range(2,40),add=TRUE,nlev=6,lwd=2)
map('world2',add=TRUE)
grid(col="black",lty=1)