# HW 2 Q 1 JLCY, 9 Nov 2013

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

# part i

nrows = 72
ncols = 36

# ntime = 1341 #Jan 1900 - Sep 2011
ntime = 1365  #Jan 1900 - Sep 2013

# ntimep = 1332 # Jan 1900 - Dec 2010
ntimep = 1356  # Jan 1900 - Dec 2012

nglobe = nrows * ncols
# N = nrows*ncols

nyrs = ntimep/12  #1900 - 2012 

# Lat - Long grid..

# read data online
# locs=matrix(scan('http://civil.colorado.edu/~balajir/CVEN6833/R-sessions/session2/files-4HW2/sst-lat-long.txt'),
# ncol=2, byrow=T)

# read data from pc locs=matrix(scan(file =
# 'C:\\Users\\jfklee90\\Desktop\\CVEN 6833\\HW
# 2\\sst-lat-long.txt'), ncol=2, byrow=T)

# 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]-360 #longitude on 0-360 grid if
# needed
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..

# data=readBin('http://civil.colorado.edu/~balajir/CVEN6833/R-sessions/session2/files-4HW2/Kaplan-SST-Jan1900-Sep2011.r4',what='numeric',
# n=( nrows * ncols * ntime), size=4,endian='swap')
# data=readBin('Kaplan-SST-Jan1900-Sep2011.r4',what='numeric', n=( nrows *
# ncols * ntime), size=4,endian='swap')
# data=readBin('C:\\Users\\jfklee90\\Desktop\\CVEN 6833\\R session
# 2\\Kaplan-SST-Jan1900-Sep2011.r4',what='numeric', n=( nrows * ncols *
# ntime), size=4,endian='swap')

data = readBin("C:\\Users\\jfklee90\\Desktop\\CVEN 6833\\HW 2\\Kaplan-SST-Jan1900-Sep2013.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
index1 = index[data1 < 20]

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

# x1[x1 < 0]= x1[x1 < 0] + 360 xygrid1[,2]=x1

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

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

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

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

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

## create annual average data

nyrs1 = ntimep/12  # = nyrs
sstanavg = matrix(0, nrow = nyrs1, ncol = nsites)

for (i in 1:nsites) {
    xx = t(matrix(t(sstdata[, i]), nrow = 12))  # set 12 rows at a time
    sstanavg[, i] = apply(xx, 1, mean)  # calculate xx's row mean
}

## write out the grid locations..
write(t(xygrid1), file = "C:\\Users\\jfklee90\\Desktop\\CVEN 6833\\HW 2\\kaplan-sst-locs.txt", 
    ncol = 2)

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

# PCA

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

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

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

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

plot(1:25, lambdas[1:25], type = "l", xlab = "Modes", ylab = "Frac. Var. explained", 
    main = "Egien Variance Spectrum")
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")
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(1900:2012, pcs[, 1], type = "l", main = "PC1 Temporal Mode of Variability")

par(ask = TRUE)

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

# 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")
contour(xlong, ylat, (zmat), ylim = range(-50, 70), add = TRUE, nlev = 6, lwd = 2)
map("world2Hires", add = T)

par(ask = TRUE)

# Temporal plot
plot(1900:2012, pcs[, 2], type = "l", main = "PC2 Temporal Mode of Variability")

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")
contour(xlong, ylat, (zmat), ylim = range(-50, 70), add = TRUE, nlev = 6, lwd = 2)
map("world2Hires", add = T)

par(ask = TRUE)

# Temporal plot
plot(1900:2012, pcs[, 3], type = "l", main = "PC3 Temporal Mode of Variability")

par(ask = TRUE)

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

# PC 4

# zsvd$u[,4] - third 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")
contour(xlong, ylat, (zmat), ylim = range(-50, 70), add = TRUE, nlev = 6, lwd = 2)
map("world2Hires", add = T)

par(ask = TRUE)

# Temporal plot
plot(1900:2012, pcs[, 4], type = "l", main = "PC4 Temporal Mode of Variability")

par(ask = TRUE)

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

# part ii

# remove ESNO mode: PC 2

# number of modes
nmodes = length(zsvd$u[, 2])

# modes to move: removing second mode
nrm = c(2)

# empty matrix
E = matrix(0, nrow = nmodes, ncol = nmodes)

E[, nrm] = zsvd$u[, nrm]
sstan_rm = pcs %*% t(E)

# new variable without second mode
sstanavg_rm = sstanavg - sstan_rm

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

# PCA after removed PC 2

# get variance matrix..
zs_rm = var(scale(sstanavg_rm))

# do an Eigen(Singular Value) decomposition..
zsvd_rm = svd(zs_rm)

# Principal Components...
pcs_rm = t(t(zsvd_rm$u) %*% t(sstanavg_rm))

# Eigen Values.. - fraction variance
lambdas_rm = (zsvd_rm$d/sum(zsvd_rm$d))

plot(1:25, lambdas_rm[1:25], type = "l", xlab = "Modes", ylab = "Frac. Var. explained", 
    main = "Egien Variance Spectrum - ENSO Removed")
points(1:25, lambdas_rm[1:25], col = "red")

par(ask = TRUE)

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

# PC 1

# zsvd_rm$u[,1] - first mode
zfull[indexgrid] = zsvd_rm$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 - ENSO Removed")
contour(xlong, ylat, (zmat), ylim = range(-50, 70), add = TRUE, nlev = 6, lwd = 2)
map("world2Hires", add = T)

par(ask = TRUE)

# Temporal plot
plot(1900:2012, pcs_rm[, 1], type = "l", main = "PC1 Temporal Mode of Variability - ENSO Removed")

par(ask = TRUE)

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

# PC 2

# zsvd_rm$u[,2] - second mode
zfull[indexgrid] = zsvd_rm$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 - ENSO Removed")
contour(xlong, ylat, (zmat), ylim = range(-50, 70), add = TRUE, nlev = 6, lwd = 2)
map("world2Hires", add = T)

par(ask = TRUE)

# Temporal plot
plot(1900:2012, pcs_rm[, 2], type = "l", main = "PC2 Temporal Mode of Variability - ENSO Removed")

par(ask = TRUE)

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

# PC 3

# zsvd_rm$u[,3] - third mode
zfull[indexgrid] = zsvd_rm$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 - ENSO Removed")
contour(xlong, ylat, (zmat), ylim = range(-50, 70), add = TRUE, nlev = 6, lwd = 2)
map("world2Hires", add = T)

par(ask = TRUE)

# Temporal plot
plot(1900:2012, pcs_rm[, 3], type = "l", main = "PC3 Temporal Mode of Variability - ENSO Removed")

par(ask = TRUE)

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

# PC 4

# zsvd_rm$u[,4] - third mode
zfull[indexgrid] = zsvd_rm$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 - ENSO Removed")
contour(xlong, ylat, (zmat), ylim = range(-50, 70), add = TRUE, nlev = 6, lwd = 2)
map("world2Hires", add = T)

par(ask = TRUE)

# Temporal plot
plot(1900:2012, pcs_rm[, 4], type = "l", main = "PC4 Temporal Mode of Variability - ENSO Removed")

par(ask = TRUE)

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

# part iii

# 1900 - 1935

nrows = 72
ncols = 36

ntimep = 36  # 1900 - 1935

ntime = ntimep

nglobe = nrows * ncols
# N = 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-Jan1900-Dec1935.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 a = is.nan(data1)
index1 = index[data1 != "NaN"]

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

# x1[x1 < 0]= x1[x1 < 0] + 360 xygrid1[,2]=x1

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 = t(t(zsvd$u) %*% t(sstanavg))

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

plot(1:25, lambdas[1:25], type = "l", xlab = "Modes", ylab = "Frac. Var. explained", 
    main = "Egien Variance Spectrum 1900-1935")
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
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 1900-1935")
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(1900:1935, pcs[, 1], type = "l", main = "PC1 Temporal Mode of Variability 1900-1935")

par(ask = TRUE)

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

# 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 1900-1935")
contour(xlong, ylat, (zmat), ylim = range(-50, 70), add = TRUE, nlev = 6, lwd = 2)
map("world2Hires", add = T)

par(ask = TRUE)

# Temporal plot
plot(1900:1935, pcs[, 2], type = "l", main = "PC2 Temporal Mode of Variability 1900-1935")

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 1900-1935")
contour(xlong, ylat, (zmat), ylim = range(-50, 70), add = TRUE, nlev = 6, lwd = 2)
map("world2Hires", add = T)

par(ask = TRUE)

# Temporal plot
plot(1900:1935, pcs[, 3], type = "l", main = "PC3 Temporal Mode of Variability 1900-1935")

par(ask = TRUE)

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

# PC 4

# zsvd$u[,4] - third 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 1900-1935")
contour(xlong, ylat, (zmat), ylim = range(-50, 70), add = TRUE, nlev = 6, lwd = 2)
map("world2Hires", add = T)

par(ask = TRUE)

# Temporal plot
plot(1900:1935, pcs[, 4], type = "l", main = "PC4 Temporal Mode of Variability 1900-1935")

par(ask = TRUE)

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

# 1935 - 1975

nrows = 72
ncols = 36

ntimep = 41  # 1935 - 1975

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-Jan1935-Dec1975.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 = t(t(zsvd$u) %*% t(sstanavg))

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

plot(1:25, lambdas[1:25], type = "l", xlab = "Modes", ylab = "Frac. Var. explained", 
    main = "Egien Variance Spectrum 1935-1975")
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
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 1935-1975")
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(1935:1975, pcs[, 1], type = "l", main = "PC1 Temporal Mode of Variability 1935-1975")

par(ask = TRUE)

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

# 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 1935-1975")
contour(xlong, ylat, (zmat), ylim = range(-50, 70), add = TRUE, nlev = 6, lwd = 2)
map("world2Hires", add = T)

par(ask = TRUE)

# Temporal plot
plot(1935:1975, pcs[, 2], type = "l", main = "PC2 Temporal Mode of Variability 1935-1975")

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 1935-1975")
contour(xlong, ylat, (zmat), ylim = range(-50, 70), add = TRUE, nlev = 6, lwd = 2)
map("world2Hires", add = T)

par(ask = TRUE)

# Temporal plot
plot(1935:1975, pcs[, 3], type = "l", main = "PC3 Temporal Mode of Variability 1935-1975")

par(ask = TRUE)

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

# PC 4

# zsvd$u[,4] - third 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 1935-1975")
contour(xlong, ylat, (zmat), ylim = range(-50, 70), add = TRUE, nlev = 6, lwd = 2)
map("world2Hires", add = T)

par(ask = TRUE)

# Temporal plot
plot(1935:1975, pcs[, 4], type = "l", main = "PC4 Temporal Mode of Variability 1935-1975")

par(ask = TRUE)

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

# 1975 - 2012

nrows = 72
ncols = 36

ntimep = 38  # 1975 - 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-Jan1975-Dec2012.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 = t(t(zsvd$u) %*% t(sstanavg))

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

plot(1:25, lambdas[1:25], type = "l", xlab = "Modes", ylab = "Frac. Var. explained", 
    main = "Egien Variance Spectrum 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
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 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))

library(maps)
library(mapdata)

map("world2Hires", add = T)

par(ask = TRUE)

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

par(ask = TRUE)

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

# 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 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(1975:2012, pcs[, 2], type = "l", main = "PC2 Temporal Mode of Variability 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 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(1975:2012, pcs[, 3], type = "l", main = "PC3 Temporal Mode of Variability 1975-2012")

par(ask = TRUE)

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

# PC 4

# zsvd$u[,4] - third 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 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(1975:2012, pcs[, 4], type = "l", main = "PC4 Temporal Mode of Variability 1975-2012")

par(ask = TRUE)

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