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


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


### Read in JJAS SST average and JJAS AISMR

nrows=72
ncols=36
ntime = 59    #JJAS 1951 - JJAS 2009
ntimep = 59   # JJAS 1951 - JJAS 2009
N = nrows*ncols


### Lat - Long grid..
ygrid=seq(-87.5,87.5,by=5)
ny=length(ygrid)

xgrid=seq(27.5,382.5,by=5)
#xgrid[xgrid > 180]=xgrid[xgrid > 180]-360	#longitude on 0-360 grid if needed
xgrid[xgrid > 180]=xgrid[xgrid > 180]
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 Kaplan SST data.. for 1951 - 2009

data=readBin("Kaplan-SST-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, put it to a large number..
data1[data1 == "NaN"]=1e+30	

# the lat -long data grid..

index=1:(nx*ny)

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

#x1[x1 < 0]= x1[x1 < 0] + 360
#xygrid1[,2]=x1

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

### SSTdata matrix - rows are years and columns are locations on the globe with data
sstdata=matrix(NA,nrow=ntimep, ncol=nsites)


for(i in 1:ntimep){
data1=data[,,i]
data1[data1 == "NaN"]=1e+30	
index1=index[data1 < 20]
data2=data1[index1]
sstdata[i,]=data2
}

## Index of locations corresponding to Global Tropics
### grids with non-missing data and between 25N-25S
xlongs = xygrid[,2]
ylats = xygrid[,1]
indextrop = index[data1 < 20 & ylats >= -25 & ylats <= 25]

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

### global Tropics
xlongs = xygrid1[,2]
ylats = xygrid1[,1]
xlong = xlongs[ylats >= -25 & ylats <= 25]
ylat = ylats[ylats >= -25 & ylats <= 25]
index1=1:length(xlongs)
index = index1[ylats >= -25 & ylats <= 25]

### Tropical seasonal average.. - the data already is seasonal average
sstanavgtrop = sstdata[,index]

xygridgp=cbind(ylat,xlong)

## write out the grid locations.. if needed
#write(t(xygridp),file="kaplan-sst-pac-locs.txt",ncol=2)

### set the data grid and the index
xygridtemp = xygrid
indextemp = indextrop   ##replace with index if you want global
###################### 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
####################################################
########## sstanavgtrop and raindata are the matrices of 
#######    tropical SST anomalies and grid rainfall, respectively.
##########  both for the summer season Jun-Sep average

####### indextrop and indexrain are the grid indices with data
########   

##########################################################
### (i)
###################### PCA 
## PCA on the global tropical summer SST 

#get variance matrix..
zs=var(sstanavgtrop)

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

#Principal Components...
pcs=sstanavgtrop %*%  zsvd$u 
sstpcs = 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")


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


# 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(xygridtemp[,2]))
ylat = sort(unique(xygridtemp[,1]))
nrows=length(xlong)
ncols=length(ylat)
nglobe = nrows*ncols   # also equal to 72*36

zfull = rep(NaN,nglobe)
zfull[indextrop]=zsvd$u[,2]      ###### 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(-40,40),xlab="longitude",ylab="latitude")
contour(xlong,ylat,(zmat),ylim=range(-40,40),add=TRUE,nlev=6,lwd=2)
map('world2',add=TRUE,ylim=range(-40,40))
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 SSTs - Global tropic")

### Similarly plot the other three Eigen vectors..

##########################################################

### (ii)

################# 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

# 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[,1]      ###### 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))
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")

### Similarly plot the other three Eigen vectors..

##################

### (iii)
########### scatterplot the rainfall PCs against SST PCs
#### and plot a smooth line through them

plot(sstpcs[,1],rainpcs[,1],xlab="SST PC1", ylab="Rain PC1")

plot(sstpcs[,1],rainpcs[,2],xlab="SST PC1", ylab="Rain PC2")

### similarly for PCs 2 through 4


### #########
#### (iv) correlate rainfall PC with SST data and vice versa
############# 

zcorr = cor(rainpcs[,1],sstanavgtrop)


xlong = sort(unique(xygridtemp[,2]))
ylat = sort(unique(xygridtemp[,1]))
nrows=72
ncols=36
nglobe = length(xlong)*length(ylat)   # also equal to 72*36

zfull = rep(NaN,nglobe)
zfull[indextrop]=zcorr      
zmat = matrix(zfull,nrow=nrows,ncol=ncols)

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


### similarly correlate the next three PCs and plot the correlation maps

### likewise correlate SST PC with rainfall
zcorr = cor(sstpcs[,1],raindata1)

xlong = sort(unique(xygridrain[,2]))
ylat = sort(unique(xygridrain[,1]))
nrows=35
ncols=33
nglobe = length(xlong)*length(ylat)   # also equal to 35*33

zfull = rep(NaN,nglobe)
zfull[indexrain]=zcorr     
zmat = matrix(zfull,nrow=nrows,ncol=ncols)


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


#### similarly correlate the next three PCCS andd plot the correlation maps