[Thread Prev][Thread Next][Index]
Re: more on eof calculus and missing values
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
[Thread Prev][Thread Next][Index]
Dept of Commerce /
NOAA /
OAR /
PMEL /
TMAP
Contact Us | Privacy Policy | Disclaimer | Accessibility Statement