# HW 2 Q 3 JLCY, 9 Nov 2013

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

# SST

nrows = 72
ncols = 36

ntimep = 112  # Dec 1900 to Mar 1901 - Dec 2011 to Mar 2012

ntime = ntimep

nglobe = nrows * ncols

# Lat - Long grid..  create grid below:

ygrid = seq(-87.5, 87.5, by = 5)
ny = length(ygrid)  # = 36

xgrid = seq(27.5, 382.5, by = 5)
xgrid[xgrid > 180] = xgrid[xgrid > 180]
nx = length(xgrid)  # = 72

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 -- already annual averaged

data = readBin("C:\\Users\\jfklee90\\Desktop\\CVEN 6833\\HW 2\\Kaplan-SST-DJFM-1900-2012.r4", 
    what = "numeric", n = (nrows * ncols * ntime), size = 4, endian = "swap")

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

# read data for the 1st month
data1 = data[, , 1]

# the lat -long data grid..

# all the grid points on the globe
index = 1:(nx * ny)

# only non-missing data, missing data = 1.E30 include only grid points on
# the ocean sort out NaN values
index1 = index[data1 != "NaN"]

# choose only ocean's grid points for xygrid
xygrid1 = xygrid[index1, ]
x1 = xygrid1[, 2]

nsites = length(index1)
data2 = data1[index1]

### SSTdata matrix - rows are months annd columns are locations
sstdata = matrix(0, nrow = ntimep, ncol = nsites)

for (i in 1:ntimep) {
    data1 = data[, , i]
    index1 = index[data1 != "NaN"]
    data2 = data1[index1]
    sstdata[i, ] = data2
}

sstanavg = scale(sstdata)

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

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

# PPT

ntimep = 112  # Dec 1900 to Mar 1901 - Dec 2011 to Mar 2012

ntime = ntimep

nlocation = 81  # 81 stations in W US

# Lat - Long grid..  create grid below:

# read coord online
# locs=matrix(scan('http://civil.colorado.edu/~balajir/CVEN6833/R-sessions/session2/files-4HW2/WesternUS-coords.txt'),
# ncol=5, byrow=T)

# read coord from pc
locs = matrix(scan("C:\\Users\\jfklee90\\Desktop\\CVEN 6833\\HW 2\\WesternUS-coords.txt"), 
    ncol = 5, byrow = T)

xlat = locs[, 3]
ylon = locs[, 4]

xygrid_ppt = matrix(0, nrow = 81, ncol = 2)

xygrid_ppt[, 1] = xlat
xygrid_ppt[, 2] = ylon

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

# Read W US PPT data -- already annual averaged

# data=scan('http://civil.colorado.edu/~balajir/CVEN6833/R-sessions/session2/files-4HW2/western-US-cdiv-winter-1900-2012-precip.txt')

data = scan("C:\\Users\\jfklee90\\Desktop\\CVEN 6833\\HW 2\\W-US-cdiv-ppt-DJFM-1900-2012.txt")

data_m <- array(dim = c(ntime, nlocation))

# set data matrix: 112 years x 81 locations
i = 1
j = 1
for (i in 1:81) {
    data_m[, i] = data[j:(j + 111)]
    i = i + 1
    j = j + 112
}

pptanavg_m = data_m

pptanavg = scale(pptanavg_m)

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

# Perform SVD

C = var(pptanavg, sstanavg)  #note that the variable with fewer columns should be first
zz = svd(C)

# f1_ppt and g1_sst are the matrix of time coefficients just like the PCS
f1_ppt = pptanavg %*% zz$u
g1_sst = sstanavg %*% 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 vice versa

# Singular value spectrum: J matrix
J = min(length(sstanavg[1, ]), length(pptanavg[1, ]))

# zz$d - eigenvalues
plot(1:J, ((zz$d)^2)/sum((zz$d)^2), xlab = "Modes", ylab = "Sq. Covariance", 
    col = "red", main = "Singular Value Spectrum - Winter SST & W US PPT 1900-2012")

par(ask = TRUE)

# Notice that the first mode captures 'almost' everything of the joint
# covariance, as expected

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

# Heterogeneous Correlation Map & Time Coefficient Plot

# Heterogeneous Correlation Map: Time Coefficent of PPT correlate on SST

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

## Loading required package: spam
## Loading required package: grid
## Spam version 0.30-1 (2013-08-25) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction 
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
## 
## Attaching package: 'spam'
## 
## The following object is masked from 'package:base':
## 
##     backsolve, forwardsolve

library(graphics)

# 1st Time Coefficent of PPT correlate on SST
corr_ppt1_sst = cor(f1_ppt[, 1], sstanavg)

# 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(xygrid[, 2]))
ylat = sort(unique(xygrid[, 1]))
zfull = rep(NaN, nglobe)

zfull[indexgrid] = corr_ppt1_sst[1, ]
zmat1 = matrix(zfull, nrow = nrows, ncol = ncols)

# Heterogeneous Correlation Map
image.plot(xlong, ylat, zmat1, ylim = range(-50, 70), main = "1st Mode of W US Winter PPT correlate on Winter SST - 1900-2012")
contour(xlong, ylat, (zmat1), ylim = range(-50, 70), add = TRUE, nlev = 10, 
    lwd = 2)
# world(add=TRUE,shift=TRUE,ylim=range(-50,70))

library(maps)
library(mapdata)

map("world2Hires", add = T)

par(ask = TRUE)

# Time Coefficient Plot
plot(1901:2012, f1_ppt[, 1], type = "l", main = "1st Time Coefficient - W US Winter PPT 1900-2012")

par(ask = TRUE)

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

# 2nd Time Coefficent of PPT correlate on SST
corr_ppt2_sst = cor(f1_ppt[, 2], sstanavg)

zfull[indexgrid] = corr_ppt2_sst[1, ]
zmat2 = matrix(zfull, nrow = nrows, ncol = ncols)

# Heterogeneous Correlation Map
image.plot(xlong, ylat, zmat2, ylim = range(-50, 70), main = "2nd Mode of W US Winter PPT correlate on Winter SST - 1900-2012")
contour(xlong, ylat, (zmat2), ylim = range(-50, 70), add = TRUE, nlev = 10, 
    lwd = 2)
map("world2Hires", add = T)

par(ask = TRUE)

# Time Coefficient Plot
plot(1901:2012, f1_ppt[, 2], type = "l", main = "2nd Time Coefficient - W US Winter PPT 1900-2012")

par(ask = TRUE)

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

# 3rd Time Coefficent of PPT correlate on SST
corr_ppt3_sst = cor(f1_ppt[, 3], sstanavg)

zfull[indexgrid] = corr_ppt3_sst[1, ]
zmat3 = matrix(zfull, nrow = nrows, ncol = ncols)

# Heterogeneous Correlation Map
image.plot(xlong, ylat, zmat3, ylim = range(-50, 70), main = "3rd Mode of W US Winter PPT correlate on Winter SST - 1900-2012")
contour(xlong, ylat, (zmat3), ylim = range(-50, 70), add = TRUE, nlev = 10, 
    lwd = 2)
map("world2Hires", add = T)

par(ask = TRUE)

# Time Coefficient Plot
plot(1901:2012, f1_ppt[, 3], type = "l", main = "3rd Time Coefficient - W US Winter PPT 1900-2012")

par(ask = TRUE)

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

# 4th Time Coefficent of PPT correlate on SST
corr_ppt4_sst = cor(f1_ppt[, 4], sstanavg)

zfull[indexgrid] = corr_ppt4_sst[1, ]
zmat4 = matrix(zfull, nrow = nrows, ncol = ncols)

# Heterogeneous Correlation Map
image.plot(xlong, ylat, zmat4, ylim = range(-50, 70), main = "4th Mode of W US Winter PPT correlate on Winter SST - 1900-2012")
contour(xlong, ylat, (zmat4), ylim = range(-50, 70), add = TRUE, nlev = 10, 
    lwd = 2)
map("world2Hires", add = T)

par(ask = TRUE)

# Time Coefficient Plot
plot(1901:2012, f1_ppt[, 4], type = "l", main = "4th Time Coefficient - W US Winter PPT 1900-2012")

par(ask = TRUE)

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

# Heterogeneous Correlation Map: Time Coefficent of SST correlate on PPT

# 1st Time Coefficent of SST correlate on PPT
corr_sst1_ppt = cor(g1_sst[, 1], pptanavg)

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

# the data is on a grid so fill the entire globaal grid with NaN and then
# populate the ocean grids with the Eigen vector

lats = locs[, 3]
lons = locs[, 4]
lons = -lons  #W longitude to be negative.

z_mode1 = corr_sst1_ppt[1, ]

# Heterogeneous Correlation Map
quilt.plot(lons, lats, z_mode1, xlim = range(-125, -100), main = "1st Mode of Winter SST correlate on W US Winter PPT - 1900-2012")
US(add = TRUE, col = "grey", lwd = 2)

par(ask = TRUE)

# Time Coefficient Plot
plot(1901:2012, g1_sst[, 1], type = "l", main = "1st Time Coefficient - Winter SST 1900-2012")

par(ask = TRUE)

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

# 2nd Time Coefficent of SST correlate on PPT
corr_sst2_ppt = cor(g1_sst[, 2], pptanavg)

z_mode2 = corr_sst2_ppt[1, ]

# Heterogeneous Correlation Map
quilt.plot(lons, lats, z_mode2, xlim = range(-125, -100), main = "2nd Mode of Winter SST correlate on W US Winter PPT - 1900-2012")
US(add = TRUE, col = "grey", lwd = 2)

par(ask = TRUE)

# Time Coefficient Plot
plot(1901:2012, g1_sst[, 2], type = "l", main = "2nd Time Coefficient - Winter SST 1900-2012")

par(ask = TRUE)

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

# 3rd Time Coefficent of SST correlate on PPT
corr_sst3_ppt = cor(g1_sst[, 3], pptanavg)

z_mode3 = corr_sst3_ppt[1, ]

# Heterogeneous Correlation Map
quilt.plot(lons, lats, z_mode3, xlim = range(-125, -100), main = "3rd Mode of Winter SST correlate on W US Winter PPT - 1900-2012")
US(add = TRUE, col = "grey", lwd = 2)

par(ask = TRUE)

# Time Coefficient Plot
plot(1901:2012, g1_sst[, 3], type = "l", main = "3rd Time Coefficient - Winter SST 1900-2012")

par(ask = TRUE)

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

# 4th Time Coefficent of SST correlate on PPT
corr_sst4_ppt = cor(g1_sst[, 4], pptanavg)

z_mode4 = corr_sst4_ppt[1, ]

# Heterogeneous Correlation Map
quilt.plot(lons, lats, z_mode4, xlim = range(-125, -100), main = "4th Mode of Winter SST correlate on W US Winter PPT - 1900-2012")
US(add = TRUE, col = "grey", lwd = 2)

par(ask = TRUE)

# Time Coefficient Plot
plot(1901:2012, g1_sst[, 4], type = "l", main = "4th Time Coefficient - Winter SST 1900-2012")

par(ask = TRUE)

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

# group pars of spatial maps & temporal plots for better visualization

par(mfrow = c(2, 1))

# Heterogeneous Correlation Map

# 1st mode

image.plot(xlong, ylat, zmat1, ylim = range(-50, 70), main = "1st Mode of W US Winter PPT correlate on Winter SST - 1900-2012")
contour(xlong, ylat, (zmat1), ylim = range(-50, 70), add = TRUE, nlev = 10, 
    lwd = 2)
map("world2Hires", add = T)

quilt.plot(lons, lats, z_mode1, xlim = range(-125, -100), main = "1st Mode of Winter SST correlate on W US Winter PPT - 1900-2012")
US(add = TRUE, col = "grey", lwd = 2)

par(ask = TRUE)

# 2nd mode

image.plot(xlong, ylat, zmat2, ylim = range(-50, 70), main = "2nd Mode of W US Winter PPT correlate on Winter SST - 1900-2012")
contour(xlong, ylat, (zmat2), ylim = range(-50, 70), add = TRUE, nlev = 10, 
    lwd = 2)
map("world2Hires", add = T)

quilt.plot(lons, lats, z_mode2, xlim = range(-125, -100), main = "2nd Mode of Winter SST correlate on W US Winter PPT - 1900-2012")
US(add = TRUE, col = "grey", lwd = 2)

par(ask = TRUE)

# 3rd mode

image.plot(xlong, ylat, zmat3, ylim = range(-50, 70), main = "3rd Mode of W US Winter PPT correlate on Winter SST - 1900-2012")
contour(xlong, ylat, (zmat3), ylim = range(-50, 70), add = TRUE, nlev = 10, 
    lwd = 2)
map("world2Hires", add = T)

quilt.plot(lons, lats, z_mode3, xlim = range(-125, -100), main = "3rd Mode of Winter SST correlate on W US Winter PPT - 1900-2012")
US(add = TRUE, col = "grey", lwd = 2)

par(ask = TRUE)

# 4th mode

image.plot(xlong, ylat, zmat4, ylim = range(-50, 70), main = "4th Mode of W US Winter PPT correlate on Winter SST - 1900-2012")
contour(xlong, ylat, (zmat4), ylim = range(-50, 70), add = TRUE, nlev = 10, 
    lwd = 2)
map("world2Hires", add = T)

quilt.plot(lons, lats, z_mode4, xlim = range(-125, -100), main = "4th Mode of Winter SST correlate on W US Winter PPT - 1900-2012")
US(add = TRUE, col = "grey", lwd = 2)

par(ask = TRUE)

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

# Time Coefficient Plot

# 1st mode

plot(1901:2012, g1_sst[, 1], type = "l", main = "1st Time Coefficient - Winter SST 1900-2012")

plot(1901:2012, f1_ppt[, 1], type = "l", main = "1st Time Coefficient - W US Winter PPT 1900-2012")

par(ask = TRUE)

# 2nd mode

plot(1901:2012, g1_sst[, 2], type = "l", main = "2nd Time Coefficient - Winter SST 1900-2012")

plot(1901:2012, f1_ppt[, 2], type = "l", main = "2nd Time Coefficient - W US Winter PPT 1900-2012")

par(ask = TRUE)

# 3rd mode

plot(1901:2012, g1_sst[, 3], type = "l", main = "3rd Time Coefficient - Winter SST 1900-2012")

plot(1901:2012, f1_ppt[, 3], type = "l", main = "3rd Time Coefficient - W US Winter PPT 1900-2012")

par(ask = TRUE)

# 4th mode

plot(1901:2012, g1_sst[, 4], type = "l", main = "4th Time Coefficient - Winter SST 1900-2012")

plot(1901:2012, f1_ppt[, 4], type = "l", main = "4th Time Coefficient - W US Winter PPT 1900-2012")

par(ask = TRUE)

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