[R-sig-Geo] spacetime EOF implementation - interpretation question

[Andrew Vitale]

Hello,

I'm performing an empirical orthogonal function analysis on sea level
pressure using the spacetime package, and I'm having trouble determining
which elements of the function output I should use to interpret the results.

It seems that most examples in the literature calculate empirical
orthogonal functions (EOFs) and their associated principle component (PC)
time series (e.g. http://www.met.reading.ac.uk/~han/Monitor/eofprimer.pdf ,
pp 11-12, Fig. 2, 3).  As far as I can tell, I should be able to extract
all of the relevant information using only the spatial mode in the
spacetime:::EOF function.  However, the ambiguous terminology in the
literature and the availability of the "temporal" mode in spacetime:::EOF
have me uncertain of my interpretation.

I have been treating the output of spacetime:::EOF with
returnPredictions=TRUE as the EOFs, and I have been treating the
spacetime:::EOF $rotation element with returnPredictions=FALSE as the PCs.

Am I extracting the proper information from the spacetime:::EOF results?
Are the results of the spatial and temporal modes of spacetime:::EOF both
required to perform an EOF analysis of a sea level pressure field, or are
these two separate analyses?

Below is a small example of how I am currently extracting the spatial EOFs
and their associated PCs.

Thanks for any comments,
Andrew Vitale


library(raster)
library(spacetime)

## Create a raster object that mimics my data
## which is actually a raster stack of 20,718
## daily sea level pressure anomolies
a = array(rnorm(9*9*100), c(9,9,100))
b = brick(a)
b = setValues(b, a)
z = seq(as.Date('1980-01-01'), as.Date('1980-04-09'),
        by = 'day')
b = setZ(b, z)

## coerce the raster stack to an STFDF object
st = as(b, 'STFDF')

## calculate EOFs in spatial mode using spacetime
seof = EOF(st, how='spatial', returnPredictions=FALSE, scale.=TRUE)
seof_preds = EOF(st, how='spatial', returnPredictions=TRUE, scale.=TRUE)

## Create a raster stack of the predictions, which seem to be
## generally referred to as EOFs in the literature
EOFs = stack(seof_preds)
## extract and standardize the rotation element of the
## seof object, which seems to be the eigenvectors, which seem
## to also be referred to as PCs in the literature
PCs = apply(seof$rotation, 2, function(x) (x - mean(x)) / sd(x))

## Plot the EOFs as maps (seof_preds)
plot(EOFs)

## Plot the first PC (seof$rotation [standardized])
x11(height=7, width=10)
plot(x=as.Date(row.names(PCs), format='X%Y.%m.%d'),
     y=PCs[ , 1], type='l', xlab='Date', ylab='Standardized PC1')


-- 
*Andrew P. Vitale*
Masters Student
Department of Geography
University of Nevada, Reno
vitale232 at gmail.com

	[[alternative HTML version deleted]]