This function sets up the grid for the SST data. It finds the undefined indices and then organizes the SST data for use into a matrix that is nrows x ncols where the number of rows is the number of years and the number of columns are the number of SST grid points. It outputs the SSTdata and the xygrid to be used throughout the homework.
SSTxygrid <- function(data) {
ygrid <- seq(-87.5, 87.5, by = 5)
xgrid <- seq(27.5, 382.5, by = 5)
ny <- length(ygrid) # length of ygrid = 36
nx <- length(xgrid) # length of 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]
}
}
data <- array(data = data, dim = c(Nrows, Ncols, Ntime))
xygrid1 <- xygrid[index1, ]
Nsites <- length(index1)
SSTdata = matrix(NA, nrow = Ntime, ncol = Nsites)
for (i in 1:Ntime) {
data1 <- data[, , i]
data2 <- data1[index1]
SSTdata[i, ] <- data2
}
indexgrid <- index1
rm("data")
results <- list(SSTdata = SSTdata, grid = xygrid)
return(results)
}
This function performs a full PCA analysis. You can specify the number of lambdas to use in the Eigen value spectrum. You can also decide whether or not you want to plot the temporal and spatial PC patterns. It is set-up to automatically give the first four in both temporal and spatial PCs unless otherwise stated. It is also coded to specify if your spatial pattern is 'SST' or 'WUSP' (western united states precipitation). Besides the plots, it outputs values used throughout this homework assignment (i.e. $u, $v, $d, $pcs, $lambdas, etc…)
FullSVD <- function(data, Nlamb, temporal, spatial, xygrid) {
### Perform SVD analysis
zs <- var(data)
zsvd <- svd(zs)
pcs <- t(t(zsvd$u) %*% t(data))
lambdas <- (zsvd$d/sum(zsvd$d))
### Plot Eigen Value Spectrum
par(mfrow = c(1, 1))
plot(1:Nlamb, lambdas[1:Nlamb], type = "l", xlab = "Modes/PCSs", ylab = "Fraction of Variance Explained",
main = "EVS")
points(1:Nlamb, lambdas[1:Nlamb], col = "red")
### Temporal PC Plots
if (temporal == "TRUE") {
par(mfrow = c(2, 2))
years <- dim(data)[1]
for (i in 1:4) {
plot(1:years, pcs[, i], type = "l", ylab = paste("PC", i))
}
par(oma = c(2, 2, 3, 2))
title("First Four Temporal PCs", outer = TRUE)
}
### Spatial PC Plots Plot Spatial SST PCs
if (spatial == "SST") {
setwd("/Users/emilygill/Documents/University of Colorado/Fall 2013/CVEN6833/HW2")
# index1 <- scan('index1.txt')
Nglobe <- 72 * 36
# set up xygrid
ygrid <- seq(-87.5, 87.5, by = 5)
ny <- length(ygrid) # length of ygrid = 36
xgrid <- seq(27.5, 382.5, by = 5)
nx <- length(xgrid) # length of 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]
}
}
# prepare to plot
par(mfrow = c(2, 2))
titles <- c("1st PC", "2nd PC", "3rd PC", "4th PC")
for (i in 1:4) {
xlong <- sort(unique(xygrid[, 2]))
ylat <- sort(unique(xygrid[, 1]))
zfull <- rep(NaN, Nglobe)
zfull[index1] <- zsvd$u[, i]
zmat <- matrix(zfull, nrow = 72, ncol = 36)
image.plot(xlong, ylat, zmat, ylim = range(-40, 70), main = titles[i])
contour(xlong, ylat, (zmat), ylim = range(-40, 70), add = TRUE,
nlev = 6, lwd = 2)
map(database = "world2", add = TRUE)
}
par(oma = c(2, 2, 3, 2))
title("First Four Spatial PCs for SST", outer = TRUE)
}
# Plot Spatial W. US Precip PCs
if (spatial == "WUSP") {
par(mfrow = c(2, 2))
titles <- c("1st PC", "2nd PC", "3rd PC", "4th PC")
for (i in 1:4) {
# quilt.plot(lons, lats, zsvd.P$u[,i],xlim=range(-125,-100),main=titles[i])
# US(add=TRUE, col='grey', lwd=2,xlim=range(-125,-100)) fit <-
# Tps(latlonUS,zsvd.P$u[,i])
fit <- interp(lons, lats, zsvd$u[, i])
surface(fit, main = titles[i])
map("world", add = TRUE)
}
par(oma = c(2, 2, 3, 2))
title("First Four Spatial PCs for Precipitation", outer = TRUE)
}
# # Plot Spatial AZ PCs if (spatial == 'AZP'){ par(mfrow=c(2,2)) titles <-
# c('1st PC','2nd PC','3rd PC','4th PC')
# }
results <- list(u = zsvd$u, v = zsvd$v, d = zsvd$d, pcs = pcs, lambdas = lambdas)
}