Re: more on eof calculus and missing values

[Billy Kessler <kessler@xxxxxxxxxxxxx>]

Just to add briefly to this discussion.

The SVD code I use comes from Numerical Recipes (section 2.9) 
which has a good discussion of this. It is definitely a faster 
computation than the "traditional" eigenvalue decomposition 
used for the presently-available routine. It has a significant
disadvantage, however, which has delayed its implementation as
a Ferret routine: the input data matrix must be larger in the
time dimension than in the space dimension. (As Mick pointed
out, all the spatial dimensions are collapsed into one, so the
input to these is a (space,time) array). This restriction means
that the routine is not generally applicable. It can be worked 
around by transposing the array and retransposing after, but I
haven't had time to do implement that.

I don't know any way to do the SVD algorithm with gappy time
series. Perhaps someone else does.

Re Mick's comment about the objective estimate of the time function
with missing data. The estimate is best in the sense of making the
least change to the overall variance of the time series of each
EOF mode.

The algorithm is from Chelton and Davis (1982) JPO, 12, 757-784.
See the appendix to that paper, which is an objective technique
for low-pass filtering in the presence of gaps. That is essentially
what is needed to find the time functions, which are the projections
of each data point onto each eigenvector, and are thus a sum over 
ALL spatial locations for each time realization. If ANY gridpoint 
has a missing value at time t_1, then the sum cannot be found, and 
the gap appears in EVERY time series at t_1. Thus a few scattered 
gaps decimate the entire calculation. What the algorithm does is 
to estimate the best "low-pass filter" of the product eigenvector*
(data value) across the width of the gap, and estimate the error 
involved. The error depends on the length of the gap and the level 
of high-frequency variability of the time series that is unsampled 
due to the gap. The calculation is about NT*NX^3 long. (Of course 
if there are no gaps, that loop is skipped). 

Billy K