### Western US Climate Division Data 81- locations
###

wuplocs = matrix(scan("http://civil.colorado.edu/~balajir/CVEN6833/R-sessions/session3/WesternUS-coords.txt"),ncol=5,byrow=T)
nlocs = dim(wuplocs)[1]
xlats = wuplocs[,3]
xlongs = -wuplocs[,4]

nlocs1=nlocs+1
winterprecip = matrix(scan("http://civil.colorado.edu/~balajir/CVEN6833/R-sessions/session3/WesternUS-winterprecip-data.txt"),ncol=nlocs1,byrow=T)
years = winterprecip[,1]
winterprecip = winterprecip[,2:nlocs1]		#first column is year

### Perform PCA

#get variance matrix..
winterprecip1 = scale(winterprecip)
zs=var(winterprecip1)

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

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

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

## Plot for the first 25 fraction Eigen values
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 spatial component or Eigen Vector pattern..
library(maps)
library(akima)
library(fields)


xx1 = cbind(xlongs, xlats)
zz=Tps(xx1,zsvd$u[,1])
surface(zz, xlab="Longitude", ylab="Latitude",xlim=range(-125,-95),ylim=range(30,50))
map('usa', add=TRUE)
### Similarly plot the other two Eigen vectors..


## Plot the First Principal Component 
plot(years, pcs[,1],type="l",xlab="Year",ylab="PC1")

## Similarly plot the first three components