# HW 2 Q 2 JLCY, 9 Nov 2013

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

# part i

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

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

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

# PCA

sstanavg = sstdata

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

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

# Principal Components...
pcs_sst = t(t(zsvd$u) %*% t(sstanavg))

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

# Eigen Spectrum
plot(1:25, lambdas[1:25], type = "l", xlab = "Modes", ylab = "Frac. Var. explained", 
    main = "Egien Variance Spectrum - Winter SST 1900-2012")
points(1:25, lambdas[1:25], col = "red")

par(ask = TRUE)

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

# Plotting

# PC 1

# plot the first spatial component or Eigen Vector pattern..
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

# 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)  #also equal 72*36

# zsvd$u[,1] - first mode
zfull[indexgrid] = zsvd$u[, 1]
zmat = matrix(zfull, nrow = nrows, ncol = ncols)

# Spatial plot
image.plot(xlong, ylat, zmat, ylim = range(-50, 70), main = "E1 Spatial Mode of Variability - Winter SST 1900-2012")
contour(xlong, ylat, (zmat), ylim = range(-50, 70), add = TRUE, nlev = 6, lwd = 2)
# world(add=TRUE,shift=TRUE,ylim=range(-50,70))

library(maps)
library(mapdata)

map("world2Hires", add = T)

par(ask = TRUE)

# Temporal plot
plot(1901:2012, pcs_sst[, 1], type = "l", main = "PC1 Temporal Mode of Variability - Winter SST 1900-2012")

par(ask = TRUE)


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

# check in expert

# SOURCES .KAPLAN .EXTENDED .v2 .ssta T (Dec 1900) (Dec 2012) RANGE T 4
# boxAverage T 12 STEP [X Y][T]svd .Ss ev 1 RANGE /color_smoothing null def
# X Y fig- colors | contours coasts -fig

# SOURCES .KAPLAN .EXTENDED .v2 .ssta T (Dec 1900) (Dec 2012) RANGE T 4
# boxAverage T 12 STEP [X Y][T]svd .Ts ev 1 RANGE

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

# PC 2

# zsvd$u[,2] - second mode
zfull[indexgrid] = zsvd$u[, 2]
zmat = matrix(zfull, nrow = nrows, ncol = ncols)

# Spatial plot
image.plot(xlong, ylat, zmat, ylim = range(-50, 70), main = "E2 Spatial Mode of Variability - Winter SST 1900-2012")
contour(xlong, ylat, (zmat), ylim = range(-50, 70), add = TRUE, nlev = 6, lwd = 2)
map("world2Hires", add = T)

par(ask = TRUE)

# Temporal plot
plot(1901:2012, pcs_sst[, 2], type = "l", main = "PC2 Temporal Mode of Variability - Winter SST 1900-2012")

par(ask = TRUE)

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

# PC 3

# zsvd$u[,3] - third mode
zfull[indexgrid] = zsvd$u[, 3]
zmat = matrix(zfull, nrow = nrows, ncol = ncols)

# Spatial plot
image.plot(xlong, ylat, zmat, ylim = range(-50, 70), main = "E3 Spatial Mode of Variability - Winter SST 1900-2012")
contour(xlong, ylat, (zmat), ylim = range(-50, 70), add = TRUE, nlev = 6, lwd = 2)
map("world2Hires", add = T)

par(ask = TRUE)

# Temporal plot
plot(1901:2012, pcs_sst[, 3], type = "l", main = "PC3 Temporal Mode of Variability - Winter SST 1900-2012")

par(ask = TRUE)

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

# PC 4

# zsvd$u[,4] - fourth mode
zfull[indexgrid] = zsvd$u[, 4]
zmat = matrix(zfull, nrow = nrows, ncol = ncols)

# Spatial plot
image.plot(xlong, ylat, zmat, ylim = range(-50, 70), main = "E4 Spatial Mode of Variability - Winter SST 1900-2012")
contour(xlong, ylat, (zmat), ylim = range(-50, 70), add = TRUE, nlev = 6, lwd = 2)
map("world2Hires", add = T)

par(ask = TRUE)

# Temporal plot
plot(1901:2012, pcs_sst[, 4], type = "l", main = "PC4 Temporal Mode of Variability - Winter SST 1900-2012")

par(ask = TRUE)

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

# 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 = matrix(0, nrow = 81, ncol = 2)

xygrid[, 1] = xlat
xygrid[, 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
}

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

# PCA

pptanavg_m = data_m

pptanavg = scale(pptanavg_m)

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

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

# Principal Components...
pcs_ppt = t(t(zsvd$u) %*% t(pptanavg))

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

# Eigen Spectrum
plot(1:25, lambdas[1:25], type = "l", xlab = "Modes", ylab = "Frac. Var. explained", 
    main = "Egien Variance Spectrum - Winter W US PPT 1900-2012")
points(1:25, lambdas[1:25], col = "red")

par(ask = TRUE)

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

# Plotting

# PC 1

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

# zsvd$u[,1] - first mode
z_mode = zsvd$u[, 1]

# Spatial plot
quilt.plot(lons, lats, z_mode, xlim = range(-125, -100), main = "E1 Spatial Mode of Variability - Winter W US PPT 1900-2012")
US(add = TRUE, col = "grey", lwd = 2)

par(ask = TRUE)

# Temporal plot
plot(1901:2012, pcs_ppt[, 1], type = "l", main = "PC1 Temporal Mode of Variability - Winter SST W US PPT 1900-2012")

par(ask = TRUE)

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

# check from email on 31 Oct 2013

# test=matrix(scan('http://civil.colorado.edu/~balajir/CVEN6833/R-sessions/session2/files-4HW2/western-US-cdiv-winter-1900-2012-precip.txt'),ncol=112,byrow=T)
# test1=scale(t(test)) zz=svd(var(test1))

# plot(1:25,(zz$d/sum(zz$d))[1:25],type='b')

# wcdivlocs =
# matrix(scan('http://civil.colorado.edu/~balajir/CVEN6833/R-sessions/session2/files-4HW2/WesternUS-coords.txt'),ncol=5,byrow=T)
# lats = wcdivlocs[,3] lons = wcdivlocs[,4] lons = -lons #W longitude to be
# negative. quilt.plot(lons, lats, zz$u[,1],xlim=range(-125,-100))
# US(add=TRUE, col='grey', lwd=2,xlim=range(-125,-100))

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

# zsvd$u[,2] - second mode
z_mode = zsvd$u[, 2]

# Spatial plot
quilt.plot(lons, lats, z_mode, xlim = range(-125, -100), main = "E2 Spatial Mode of Variability - Winter W US PPT 1900-2012")
US(add = TRUE, col = "grey", lwd = 2)

par(ask = TRUE)

# Temporal plot
plot(1901:2012, pcs_ppt[, 2], type = "l", main = "PC2 Temporal Mode of Variability - Winter SST W US PPT 1900-2012")

par(ask = TRUE)

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

# zsvd$u[,3] - third mode
z_mode = zsvd$u[, 3]

# Spatial plot
quilt.plot(lons, lats, z_mode, xlim = range(-125, -100), main = "E3 Spatial Mode of Variability - Winter W US PPT 1900-2012")
US(add = TRUE, col = "grey", lwd = 2)

par(ask = TRUE)

# Temporal plot
plot(1901:2012, pcs_ppt[, 3], type = "l", main = "PC3 Temporal Mode of Variability - Winter SST W US PPT 1900-2012")

par(ask = TRUE)

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

# zsvd$u[,4] - fourth mode
z_mode = zsvd$u[, 4]

# Spatial plot
quilt.plot(lons, lats, z_mode, xlim = range(-125, -100), main = "E4 Spatial Mode of Variability - Winter W US PPT 1900-2012")
US(add = TRUE, col = "grey", lwd = 2)

par(ask = TRUE)

# Temporal plot
plot(1901:2012, pcs_ppt[, 4], type = "l", main = "PC4 Temporal Mode of Variability - Winter SST W US PPT 1900-2012")

par(ask = TRUE)

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

# part ii

# plotting 3 x 3 = 9 graphs of the 1st 3 pairs of pcs_sst & pcs_ppt

par(mfrow = c(3, 3))

plot(pcs_sst[, 1], pcs_ppt[, 1], main = "PC1(Winter SST) & PC1(W US Winter PPT)")

plot(pcs_sst[, 2], pcs_ppt[, 1], main = "PC2(Winter SST) & PC1(W US Winter PPT)")

plot(pcs_sst[, 3], pcs_ppt[, 1], main = "PC3(Winter SST) & PC1(W US Winter PPT)")

plot(pcs_sst[, 1], pcs_ppt[, 2], main = "PC1(Winter SST) & PC2(W US Winter PPT)")

plot(pcs_sst[, 2], pcs_ppt[, 2], main = "PC2(Winter SST) & PC2(W US Winter PPT)")

plot(pcs_sst[, 3], pcs_ppt[, 2], main = "PC3(Winter SST) & PC2(W US Winter PPT)")

plot(pcs_sst[, 1], pcs_ppt[, 3], main = "PC1(Winter SST) & PC3(W US Winter PPT)")

plot(pcs_sst[, 2], pcs_ppt[, 3], main = "PC2(Winter SST) & PC3(W US Winter PPT)")

plot(pcs_sst[, 3], pcs_ppt[, 3], main = "PC3(Winter SST) & PC3(W US Winter PPT)")

par(ask = TRUE)

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

# locfit

# load package
library(locfit)

## locfit 1.5-9.1    2013-03-22

# plotting 3 x 3 = 9 graphs of the 1st 3 pairs of pcs_sst & locfit(pcs_ppt)

par(mfrow = c(3, 3))

lf_p1_s1 = locfit(pcs_ppt[, 1] ~ pcs_sst[, 1])
yest_lf_p1_s1 = 1:ntimep
yest_lf_p1_s1 = predict(lf_p1_s1, pcs_sst[, 1], se.fit = T)

plot(pcs_sst[, 1], yest_lf_p1_s1$fit, main = "PC1(Winter SST) & Locfit PC1(W US Winter PPT)")

lf_p1_s2 = locfit(pcs_ppt[, 1] ~ pcs_sst[, 2])
yest_lf_p1_s2 = 1:ntimep
yest_lf_p1_s2 = predict(lf_p1_s2, pcs_sst[, 2], se.fit = T)

plot(pcs_sst[, 2], yest_lf_p1_s2$fit, main = "PC2(Winter SST) & Locfit PC1(W US Winter PPT)")

lf_p1_s3 = locfit(pcs_ppt[, 1] ~ pcs_sst[, 3])
yest_lf_p1_s3 = 1:ntimep
yest_lf_p1_s3 = predict(lf_p1_s3, pcs_sst[, 3], se.fit = T)

plot(pcs_sst[, 3], yest_lf_p1_s3$fit, main = "PC3(Winter SST) & Locfit PC1(W US Winter PPT)")

lf_p2_s1 = locfit(pcs_ppt[, 2] ~ pcs_sst[, 1])
yest_lf_p2_s1 = 1:ntimep
yest_lf_p2_s1 = predict(lf_p2_s1, pcs_sst[, 1], se.fit = T)

plot(pcs_sst[, 1], yest_lf_p2_s1$fit, main = "PC1(Winter SST) & Locfit PC2(W US Winter PPT)")

lf_p2_s2 = locfit(pcs_ppt[, 2] ~ pcs_sst[, 2])
yest_lf_p2_s2 = 1:ntimep
yest_lf_p2_s2 = predict(lf_p2_s2, pcs_sst[, 2], se.fit = T)

plot(pcs_sst[, 2], yest_lf_p2_s2$fit, main = "PC2(Winter SST) & Locfit PC2(W US Winter PPT)")

lf_p2_s3 = locfit(pcs_ppt[, 2] ~ pcs_sst[, 3])
yest_lf_p2_s3 = 1:ntimep
yest_lf_p2_s3 = predict(lf_p2_s3, pcs_sst[, 3], se.fit = T)

plot(pcs_sst[, 3], yest_lf_p2_s3$fit, main = "PC3(Winter SST) & Locfit PC2(W US Winter PPT)")

lf_p3_s1 = locfit(pcs_ppt[, 3] ~ pcs_sst[, 1])
yest_lf_p3_s1 = 1:ntimep
yest_lf_p3_s1 = predict(lf_p3_s1, pcs_sst[, 1], se.fit = T)

plot(pcs_sst[, 1], yest_lf_p3_s1$fit, main = "PC1(Winter SST) & Locfit PC3(W US Winter PPT)")

lf_p3_s2 = locfit(pcs_ppt[, 3] ~ pcs_sst[, 2])
yest_lf_p3_s2 = 1:ntimep
yest_lf_p3_s2 = predict(lf_p3_s2, pcs_sst[, 2], se.fit = T)

plot(pcs_sst[, 2], yest_lf_p3_s2$fit, main = "PC2(Winter SST) & Locfit PC3(W US Winter PPT)")

lf_p3_s3 = locfit(pcs_ppt[, 3] ~ pcs_sst[, 3])
yest_lf_p3_s3 = 1:ntimep
yest_lf_p3_s3 = predict(lf_p3_s3, pcs_sst[, 3], se.fit = T)

plot(pcs_sst[, 3], yest_lf_p3_s3$fit, main = "PC3(Winter SST) & Locfit PC3(W US Winter PPT)")

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

# part iii

# plotting 2 x 3 = 6 graphs 1st row: 3 pairs of correlations between 1st 3
# pcs_sst & pptanavg 2st row: 3 pairs of correlations between 1st 3 pcs_ppt
# & sstanavg

par(mfrow = c(2, 3))

# correlate PC1(SST) with PPT
r_pc1sst_ppt = cor(pcs_sst[, 1], pptanavg)

# Spatial plot
quilt.plot(lons, lats, r_pc1sst_ppt[1, ], xlim = range(-125, -100), main = "Corrleation map - PC1(Winter SST) & W US Winter PPT")
US(add = TRUE, col = "grey", lwd = 2)

par(ask = TRUE)

# correlate PC2(SST) with PPT
r_pc2sst_ppt = cor(pcs_sst[, 2], pptanavg)

# Spatial plot
quilt.plot(lons, lats, r_pc2sst_ppt[1, ], xlim = range(-125, -100), main = "Corrleation map - PC2(Winter SST) & W US Winter PPT")
US(add = TRUE, col = "grey", lwd = 2)

par(ask = TRUE)

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

# check with R code

# pc2tpcpcorr=scan('http://iridl.ldeo.columbia.edu/expert/SOURCES/.NOAA/.NCDC/.CIRS/.ClimateDivision/.pcp/T/%28Nov%201900%29%28Dec%202012%29RANGE/T/12/splitstreamgrid/T/%28Nov%29%28Dec%29%28Jan%29%28Feb%29%28Mar%29VALUES%5BT%5Daverage/IDIV/201/202/203/204/205/206/207/401/402/403/404/405/406/407/502/504/505/1001/1002/1003/1004/1005/1006/1007/1008/1009/1010/2401/2402/2403/2404/2405/2406/2407/2601/2602/2603/2604/2901/2902/2904/2905/2906/2907/2908/3501/3502/3503/3504/3505/3506/3507/3508/3509/4201/4202/4203/4204/4205/4206/4207/4501/4502/4503/4504/4505/4506/4507/4508/4509/4510/4801/4802/4803/4804/4805/4806/4807/4808/4809/4810/VALUES/T2//T/renameGRID/SOURCES/.KAPLAN/.EXTENDED/.v2/.ssta/T/%28Nov%201900%29%28Mar%202012%29RANGE/T/12/splitstreamgrid/T/%28Nov%29%28Dec%29%28Jan%29%28Feb%29%28Mar%29VALUES%5BT%5Daverage/T2//T/renameGRID%5BT%5Dstandardize%5BX/Y%5D%5BT%5Dsvd/ev/2/RANGE/.Ts%5BT%5Dcorrelate/data.ch')

# library(fields) wcdivlocs =
# matrix(scan('http://civil.colorado.edu/~balajir/CVEN6833/R-sessions/session2/files-4HW2/WesternUS-coords.txt'),ncol=5,byrow=T)
# lats = wcdivlocs[,3] lons = wcdivlocs[,4] lons = -lons #W longitude to be
# negative. quilt.plot(lons, lats, pc2tpcpcorr,xlim=range(-125,-100))
# US(add=TRUE, col='grey', lwd=2,xlim=range(-125,-100))

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

# correlate PC3(SST) with PPT
r_pc3sst_ppt = cor(pcs_sst[, 3], pptanavg)

# Spatial plot
quilt.plot(lons, lats, r_pc3sst_ppt[1, ], xlim = range(-125, -100), main = "Corrleation map - PC3(Winter SST) & W US Winter PPT")
US(add = TRUE, col = "grey", lwd = 2)

par(ask = TRUE)

# correlate PC1(PPT) with SST
r_pc1ppt_sst = cor(pcs_ppt[, 1], sstanavg)

zfull = rep(NaN, nglobe)  #also equal 72*36
zfull[indexgrid] = r_pc1ppt_sst[1, ]
zmat = matrix(zfull, nrow = nrows, ncol = ncols)

image.plot(xlong, ylat, zmat, ylim = range(-50, 70), main = "Corrleation map - PC1(US W Winter PPT) & Winter SST")
contour(xlong, ylat, (zmat), ylim = range(-50, 70), add = TRUE, nlev = 6, lwd = 2)
map("world2Hires", add = T)

par(ask = TRUE)

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

# check in expert

# expert SOURCES .KAPLAN .EXTENDED .v2 .ssta T (Nov 1900) (Mar 2012) RANGE T
# 12 splitstreamgrid T (Nov) (Dec) (Jan) (Feb) (Mar) VALUES [T]average T2 /T
# renameGRID SOURCES .NOAA .NCDC .CIRS .ClimateDivision .pcp T (Nov 1900)
# (Dec 2012) RANGE T 12 splitstreamgrid T (Nov) (Dec) (Jan) (Feb) (Mar)
# VALUES [T]average T2 /T renameGRID IDIV 201 202 203 204 205 206 207 401
# 402 403 404 405 406 407 502 504 505 1001 1002 1003 1004 1005 1006 1007
# 1008 1009 1010 2401 2402 2403 2404 2405 2406 2407 2601 2602 2603 2604 2901
# 2902 2904 2905 2906 2907 2908 3501 3502 3503 3504 3505 3506 3507 3508 3509
# 4201 4202 4203 4204 4205 4206 4207 4501 4502 4503 4504 4505 4506 4507 4508
# 4509 4510 4801 4802 4803 4804 4805 4806 4807 4808 4809 4810 VALUES
# [T]standardize [IDIV][T]svd ev 1 RANGE .Ts T 2 index .T replaceGRID
# [T]correlate /color_smoothing null def X Y fig- colors | contours coasts
# -fig

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

# correlate PC2(PPT) with SST
r_pc2ppt_sst = cor(pcs_ppt[, 2], sstanavg)

zfull = rep(NaN, nglobe)  #also equal 72*36
zfull[indexgrid] = r_pc2ppt_sst[1, ]
zmat = matrix(zfull, nrow = nrows, ncol = ncols)

image.plot(xlong, ylat, zmat, ylim = range(-50, 70), main = "Corrleation map - PC2(US W Winter PPT) & Winter SST")
contour(xlong, ylat, (zmat), ylim = range(-50, 70), add = TRUE, nlev = 6, lwd = 2)
map("world2Hires", add = T)

par(ask = TRUE)

# correlate PC3(PPT) with SST
r_pc3ppt_sst = cor(pcs_ppt[, 3], sstanavg)

zfull = rep(NaN, nglobe)  #also equal 72*36
zfull[indexgrid] = r_pc3ppt_sst[1, ]
zmat = matrix(zfull, nrow = nrows, ncol = ncols)

image.plot(xlong, ylat, zmat, ylim = range(-50, 70), main = "Corrleation map - PC3(US W Winter PPT) & Winter SST")
contour(xlong, ylat, (zmat), ylim = range(-50, 70), add = TRUE, nlev = 6, lwd = 2)
map("world2Hires", add = T)

par(ask = TRUE)

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

# part iv

# SST

nrows = 72
ncols = 36

ntimep = 37  # Dec 1975 to Mar 1976 - Dec 2011 to Dec 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('http://iridl.ldeo.columbia.edu/expert/SOURCES/.KAPLAN/.EXTENDED/.v2/.ssta/T/%28Dec%201975%29%28Dec%202012%29RANGE/T/4/boxAverage/T/12/STEP/data.r4',what='numeric',
# n=( nrows * ncols * ntime), size=4, endian='swap')

data = readBin("C:\\Users\\jfklee90\\Desktop\\CVEN 6833\\HW 2\\Kaplan-SST-DJFM-1975-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
}

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

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

# PCA

sstanavg = sstdata

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

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

# Principal Components...
pcs_sst = t(t(zsvd$u) %*% t(sstanavg))

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

par(mfrow = c(1, 1))

# Eigen Spectrum
plot(1:25, lambdas[1:25], type = "l", xlab = "Modes", ylab = "Frac. Var. explained", 
    main = "Egien Variance Spectrum - Winter SST 1975-2012")
points(1:25, lambdas[1:25], col = "red")

par(ask = TRUE)

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

# check in expert

# SOURCES .KAPLAN .EXTENDED .v2 .ssta T (Dec 1975) (Dec 2012) RANGE T 4
# boxAverage T 12 STEP [T]standardize [X Y][T]svd .evaln ev 1 25 RANGE

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

# Plotting

# PC 1

# plot the first spatial component or Eigen Vector pattern..

# 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)  #also equal 72*36

# zsvd$u[,1] - first mode
zfull[indexgrid] = zsvd$u[, 1]
zmat = matrix(zfull, nrow = nrows, ncol = ncols)

# Spatial plot
image.plot(xlong, ylat, zmat, ylim = range(-50, 70), main = "E1 Spatial Mode of Variability - Winter SST 1975-2012")
contour(xlong, ylat, (zmat), ylim = range(-50, 70), add = TRUE, nlev = 6, lwd = 2)
map("world2Hires", add = T)

par(ask = TRUE)

# Temporal plot
plot(1976:2012, pcs_sst[, 1], type = "l", main = "PC1 Temporal Mode of Variability - Winter SST 1975-2012")

par(ask = TRUE)

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

# check in expert

# SOURCES .KAPLAN .EXTENDED .v2 .ssta T (Dec 1975) (Dec 2012) RANGE T 4
# boxAverage T 12 STEP [T]standardize [X Y][T]svd .Ss ev 1 RANGE
# /color_smoothing null def X Y fig- colors | contours coasts -fig

# SOURCES .KAPLAN .EXTENDED .v2 .ssta T (Dec 1975) (Dec 2012) RANGE T 4
# boxAverage T 12 STEP [T]standardize [X Y][T]svd .Ts ev 1 RANGE

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

# PC 2

# zsvd$u[,2] - second mode
zfull[indexgrid] = zsvd$u[, 2]
zmat = matrix(zfull, nrow = nrows, ncol = ncols)

# Spatial plot
image.plot(xlong, ylat, zmat, ylim = range(-50, 70), main = "E2 Spatial Mode of Variability - Winter SST 1975-2012")
contour(xlong, ylat, (zmat), ylim = range(-50, 70), add = TRUE, nlev = 6, lwd = 2)
# world(add=TRUE,shift=TRUE,ylim=range(-50,70))
map("world2Hires", add = T)

par(ask = TRUE)

# Temporal plot
plot(1976:2012, pcs_sst[, 2], type = "l", main = "PC2 Temporal Mode of Variability - Winter SST 1975-2012")

par(ask = TRUE)

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

# PC 3

# zsvd$u[,3] - third mode
zfull[indexgrid] = zsvd$u[, 3]
zmat = matrix(zfull, nrow = nrows, ncol = ncols)

# Spatial plot
image.plot(xlong, ylat, zmat, ylim = range(-50, 70), main = "E3 Spatial Mode of Variability - Winter SST 1975-2012")
contour(xlong, ylat, (zmat), ylim = range(-50, 70), add = TRUE, nlev = 6, lwd = 2)
# world(add=TRUE,shift=TRUE,ylim=range(-50,70))
map("world2Hires", add = T)

par(ask = TRUE)

# Temporal plot
plot(1976:2012, pcs_sst[, 3], type = "l", main = "PC3 Temporal Mode of Variability - Winter SST 1975-2012")

par(ask = TRUE)

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

# PC 4

# zsvd$u[,4] - fourth mode
zfull[indexgrid] = zsvd$u[, 4]
zmat = matrix(zfull, nrow = nrows, ncol = ncols)

# Spatial plot
image.plot(xlong, ylat, zmat, ylim = range(-50, 70), main = "E4 Spatial Mode of Variability - Winter SST 1975-2012")
contour(xlong, ylat, (zmat), ylim = range(-50, 70), add = TRUE, nlev = 6, lwd = 2)
# world(add=TRUE,shift=TRUE,ylim=range(-50,70))
map("world2Hires", add = T)

par(ask = TRUE)

# Temporal plot
plot(1976:2012, pcs_sst[, 4], type = "l", main = "PC4 Temporal Mode of Variability - Winter SST 1975-2012")

par(ask = TRUE)

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

# PPT

ntimep = 37  # Dec 1975 to Mar 1976 - 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 = matrix(0, nrow = 81, ncol = 2)

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

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

# Read W US PPT data -- already annual averaged

# data=scan('http://iridl.ldeo.columbia.edu/expert/SOURCES/.NOAA/.NCDC/.CIRS/.ClimateDivision/IDIV/201/202/203/204/205/206/207/401/402/403/404/405/406/407/502/504/505/1001/1002/1003/1004/1005/1006/1007/1008/1009/1010/2401/2402/2403/2404/2405/2406/2407/2601/2602/2603/2604/2901/2902/2904/2905/2906/2907/2908/3501/3502/3503/3504/3505/3506/3507/3508/3509/4201/4202/4203/4204/4205/4206/4207/4501/4502/4503/4504/4505/4506/4507/4508/4509/4510/4801/4802/4803/4804/4805/4806/4807/4808/4809/4810/VALUES/.pcp/T/%28Dec%201975%29%28Dec%202012%29RANGE/T/4/boxAverage/T/12/STEP/data.ch')

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

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

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

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

# PCA

pptanavg_m = data_m

pptanavg = scale(pptanavg_m)

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

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

# Principal Components...
pcs_ppt = t(t(zsvd$u) %*% t(pptanavg))

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

# Eigen Spectrum
plot(1:25, lambdas[1:25], type = "l", xlab = "Modes", ylab = "Frac. Var. explained", 
    main = "Egien Variance Spectrum - Winter W US PPT 1975-2012")
points(1:25, lambdas[1:25], col = "red")

par(ask = TRUE)

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

# Plotting

# PC 1

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

# zsvd$u[,1] - first mode
z_mode = zsvd$u[, 1]

# Spatial plot
quilt.plot(lons, lats, z_mode, xlim = range(-125, -100), main = "E1 Spatial Mode of Variability - Winter W US PPT 1975-2012")
US(add = TRUE, col = "grey", lwd = 2)

par(ask = TRUE)

# Temporal plot
plot(1976:2012, pcs_ppt[, 1], type = "l", main = "PC1 Temporal Mode of Variability - Winter SST W US PPT 1975-2012")

par(ask = TRUE)

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

# zsvd$u[,2] - second mode
z_mode = zsvd$u[, 2]

# Spatial plot
quilt.plot(lons, lats, z_mode, xlim = range(-125, -100), main = "E2 Spatial Mode of Variability - Winter W US PPT 1975-2012")
US(add = TRUE, col = "grey", lwd = 2)

par(ask = TRUE)

# Temporal plot
plot(1976:2012, pcs_ppt[, 2], type = "l", main = "PC2 Temporal Mode of Variability - Winter SST W US PPT 1975-2012")

par(ask = TRUE)

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

# zsvd$u[,3] - third mode
z_mode = zsvd$u[, 3]

# Spatial plot
quilt.plot(lons, lats, z_mode, xlim = range(-125, -100), main = "E3 Spatial Mode of Variability - Winter W US PPT 1975-2012")
US(add = TRUE, col = "grey", lwd = 2)

par(ask = TRUE)

# Temporal plot
plot(1976:2012, pcs_ppt[, 3], type = "l", main = "PC3 Temporal Mode of Variability - Winter SST W US PPT 1975-2012")

par(ask = TRUE)

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

# zsvd$u[,4] - fourth mode
z_mode = zsvd$u[, 4]

# Spatial plot
quilt.plot(lons, lats, z_mode, xlim = range(-125, -100), main = "E4 Spatial Mode of Variability - Winter W US PPT 1975-2012")
US(add = TRUE, col = "grey", lwd = 2)

par(ask = TRUE)

# Temporal plot
plot(1976:2012, pcs_ppt[, 4], type = "l", main = "PC4 Temporal Mode of Variability - Winter SST W US PPT 1975-2012")

par(ask = TRUE)

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