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

### Perform SVD - i.e. joint analysis


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

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


## set the data grid and index
raingrid = xygrid1
xygridrain = xygrid
indexrain = index1
##########################################################
### (i)
###################  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:25, (((zz$d)^2)/sum((zz$d)^2))[1:25], xlab="Modes", ylab="Sq. Covariance",
col="red")

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


#### (ii) 
### 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(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]=corrpats1
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)
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(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]=corrpatr1      ###### 
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)


title(main="Hetrogeneous Corr. Map - TC1 of SST with Rain")
### Plot the other hetrogeneous correlation maps..


################
#### (iii)

plot(1951:2009, f1[,1],xlab="Year",ylab="Rainfall Time Component 1", typee="b")

plot(1951:2009, g1[,1],xlab="Year",ylab="SST Time Component 1", typee="b")

## similarly the other three Time components


##### (iv)
plot(g1[,1],sstpcs[,1],xlab="SST Time component 1", ylab="SST PC1")

plot(f1[,1],rainpcs[,1],xlab="Rainfall Time component 1", ylab="Rainfall PC1")

### similarly other three components and