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

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
xlongs = xygrid[,2]
ylats = xygrid[,1]
indextrop = index[data1 < 20 & ylats >= -20 & ylats <= 20]

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

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

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

xygridgp=cbind(ylat,xlong)
xygrids = xygrid
###################### 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]

### SSTdata matrix - rows are years and columns are locations on the globe with data
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

raingrid = xygrid1
xygridr = xygrid

###################  Perform SVD....
C=var(raindata,sstanavgtrop)  #note that the variable with fewer columns should be first
zz=svd(C)

#f1 and g1 are the matrix of time coefficients just like the
#PCS
f1=raindata %*% zz$u
g1=sstanavgtrop %*% zz$v

# f1 %*% t(zz$u) will result in the flows.
# and of course, t(zz$u) %*% zz$u will result in an identfy matrix I
# similarly for g1

#f1[,1] and g1[,1] will have the 'highest' joint covariance and
# so on.. Also means, that f1[,1] is highly related to g1[,1]
# here it is similar to CCA.

#hetrogenous correlation pattern. First time
#coefficient of swe is correlated with the flows data and viceversa

J=min(length(sstanavgtrop[1,]), length(raindata[1,]))

plot(1:J, ((zz$d)^2)/sum((zz$d)^2), xlab="Modes", ylab="Sq. Covariance",
col="red")

# Notice that the first mode captures 'almost' everything  of the joint
# covariance. AS to be expected..

### Spatial map the first hetro corrln map. TC of rain correlated on SST
corrpats1=cor(f1[,1],sstanavgtrop)  

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

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

image.plot(xlong,ylat,zmat,ylim=range(-40,40))
contour(xlong,ylat,(zmat),ylim=range(-40,40),add=TRUE,nlev=6,lwd=2)
world(add=TRUE,shift=TRUE)
title(main="Hetrogeneous Corr. Map - TC1 of Rain with SST")


### Spatial map the first hetro corrln map. TC of SST correlated on Rain
corrpatr1=cor(g1[,1],raindata)


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

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

image.plot(xlong,ylat,zmat)
contour(xlong,ylat,zmat,add=TRUE,nlev=6,lwd=2)
map("world2",add=TRUE)
#world(add=TRUE,shift=TRUE)
title(main="Hetrogeneous Corr. Map - TC1 of SST with Rain")
### Plot the other hetrogeneous correlation maps..