nrows=72
ncols=36
ntime=1826	#Jan 1856 - Feb 2008
ntime = 1824    #Jan 1856 - Dec 2007
ntime = 1848    #Jan 900 - Dec 2009

ntimep = 1320   # Jan 1900 - Dec 2009
## Start annd end time points for this period
N1 = 529
N2 = 1848   


#ntimep = 1296   # Jan 1900 - Dec 2007
#N1=529
#N2=1824

nyrs = ntime/12	#1856 - 2007 or 1900 - 2009
N = nrows*ncols


### Lat - Long grid..
#missing 1e+30

locs=matrix(scan("sst-lat-long.txt"), ncol=2, byrow=T)
ygrid=seq(-87.5,87.5,by=5)
ny=length(ygrid)

# xgrid=seq(-177.5,177.5,by=5)    #longitude
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){
132
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-1856-Sep2010.r4",what="numeric", n=( nrows * ncols * ntime), size=4,endian="swap")


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


data1=data[,,N1]
	

# 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)
data2=data1[index1]

sstdata=matrix(NA,nrow=ntimep, ncol=nsites)


for(i in N1:N2){
i1=i-N1+1
data1=data[,,i]
index1=index[data1 < 20]
data2=data1[index1]
sstdata[i1,]=data2
}

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

## create annual average data
nyrs1=ntimep/12
sstanavg = matrix(0,nrow=nyrs1, ncol=nsites)

for(i in 1:nsites){
xx=t(matrix(t(sstdata[,i]),nrow=12))
sstanavg[,i]=apply(xx,1,mean)

}

###################### PCA of annual average data

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

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

#Principal Components...
pcs=t(t(zsvd$u) %*% t(sstdata))

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

plot(1:40, lambdas[1:40], type="l", xlab="Modes", ylab="Frac. Var. explained")
points(1:40, lambdas[1:40], col="red")


#plots..
#plot the first spatial component or Eigen Vector pattern..
library(maps)
library(akima)
library(fields)

#zz=interp(xlongs-360,xlats,zsvd$u[,1])
#image(zz)
#map('world',add=TRUE)


#xlongs = locs[,2]
#xlongs[xlongs <0]=xlongs[xlongs < 0]+360

#xlats = locs[,1]
#xx1 = cbind((xlongs-360), xlats)
#zz=Tps(xx1,zsvd$u[,1])
#surface(zz, xlab="Longitude", ylab="Latitude")
#map('world', add=TRUE, fill=TRUE)


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

zfull = rep(NaN,length(index))
zfull[index1]=zsvd$u[,1]
zmat = matrix(zfull,nrow=nrows,ncol=ncols)

image.plot(xlong,ylat,zmat,ylim=range(-40,70))
contour(xlong,ylat,(zmat),ylim=range(-40,70),add=TRUE,nlev=6,lwd=2)
world(add=TRUE,shift=TRUE)


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