Many analysis techniques for geophysical data require the data be located at regular intervals in space and time. This is so for spectral analysis, digital filtering and wavelet analysis. Other techniques, such as calculating empirical orthogonal functions from a set of time series observed at the same times (though not necessarily at regular intervals, i.e. t ≠ constant) often require that gaps in the data be filled. This is true, for example, in the case of satellite observations that are obscured by cloud or are subject to data drop out from instrument or algorithm failings.

In oceanography, climatologies (e.g. seasonal or monthly means) are typically computed from a compilation of observations made at irregular locations and times, frequently with a sampling distribution that can lead to regional or temporal biases is care is not taken to recognize and address these in the mapping procedure.

These next few lectures will address techniques for producing regularly gridded maps from irregularly sampled data. These techniques and have characteristics of both smoothing or filtering (removing time/space scales) and interpolation (spanning gaps in observations).

We begin by reviewing some basics of matrix and vector algebra, drawing on the “Basic Machinery” described in Chapter 3 of Wunsch (1996). [Wunsch, C., The Ocean Circulation Inverse Problem, Cambridge University Press, 442 pp., 1996.]

Matrix and Vector Algebra and Least squares fitting via the normal equations

Topics covered:

Weighted least squares

Conventional least squares via the normal equations

Least squares solution to data design matrix equation

John’s notes for lecture 1

March 25: Weighted least squares and solutions via Singular Value Decomposition

John’s notes for lecture 2

An example Matlab session to demonstrate least squares fitting using the normal equations, singular value decomposition, and matrix left divide is given in the m-file jwlsqdemo.m in the class matlab directory.

March 29: Local weighted least squares

The quadratic loess smoother

An example of weighted least squares fitting to multivariate data is the quadratic loess smoother.

In one dimension, this can be used to smooth, filter or interpolate a time series of values that may or may not be at regular sampling intervals.

An application in two or more dimensions could be to produce a gridded climatology from a set of data observed at irregularly spaced locations and times, such as a set of shipboard hydrographic observations of temperature, salinity, or other ocean properties.

The flexibility to arbitrarily choose a set of weights means the least squares fitting approach can be tailored to meet a users’ notion of what constitutes a rational weighting scheme based on some a priori knowledge of the processes being observed and modeled.

An example of this is the Climatology of the Australian Regional Seas (CARS):

Ridgway, K., J. R. Dunn and J. L. Wilkin (2002), Ocean interpolation by 4-dimensional weighted least squares: Application to the waters around Australasia, Journal of Atmospheric and Oceanic Technology, 19, 1357-1375. (pdf)

Length and time scales of variability

The quadratic loess smoother requires an a priori choice of the length or time scale to apply in the selection of the smoothing weights. The smoother can be interpreted as a filter, since the linear weighting procedure is effectively implemented as a convolution of the weights with the data. It is more general in the sense that the data do not have to be at regular intervals, because the weights are computed simply as a function of the normalized distance function r. It has been found, empirically, that the effective cutoff frequency of f_c the quadratic loess smoother when it is interpreted as a filter is f_c ≈ L^-1where L is the half width (i.e. the normalization scale) used in the loess smoother.

If the loess smoother is to be used to deliberately remove certain scales of variability (i.e. as a filter), then selection of L is straightforward. However, if the objective is to use the smoother to do the best possible job of interpolating gaps in the data, then the smoothing scale should be adapted to the natural length or time scales of variability in the underlying physical process being observed.

This brings us to method of optimal interpolation (OI), also known as objective mapping or Gauss-Markov smoothing. [See Emery and Thompson, section 4.2.] Optimal interpolation estimates the field being observed at an arbitrary location and time through a linear combination of the available data. The weights used are chosen so that the expected error of the estimate is a minimum in the least squares sense, and the estimate itself is unbiased (i.e. has the same mean as the true field). OI is therefore sometimes referred to as the Best Linear Unbiased Estimator (BLUE) of a field. The natural covariance length and time scales of the data and true field enter into the computation of the linear weights.

Optimal Interpolation

John’s notes for lecture 3

April 1: Optimal interpolation

[Also known as Objective mapping, Objective analysis or Gauss-Markov smoothing. See Emery and Thompson, section 4.2.]

Optimal interpolation estimates the field being observed at an arbitrary location and time through a linear combination of the available data. The weights used are chosen so that the expected error of the estimate is a minimum in the least squares sense, and the estimate itself is unbiased (i.e. has the same mean as the true field). OI is therefore sometimes referred to as the Best Linear Unbiased Estimator (BLUE) of a field. The natural covariance length and time scales of the data and true field enter into the computation of the linear weights.

Important concepts in optimal interpolation:

OI produces the best linear unbiased estimate of a field from a set of arbitrarily distributed observations.

Central to the estimation procedure is knowledge of the underlying covariance function of the process being observed with the observations (the model-data covariance), and the data being used with themselves (the data-data covariance). This data-data covariance includes an a priori estimate of the uncertainty (error) in the observations themselves.

These covariance patterns should be similar. If the data errors are independent, the error variance adds to the diagonal of data covariance matrix.

Frequently, the covariance is assumed to be homogenous and isotropic, in which case the covariance becomes simply a function of the distance separating the locations of the data and model (grid) points.

If valid, the assumptions of homogeneity and isotropy facilitate the estimation of the shape of the covariance function by taking an ensemble of data covariance binned according to spatial and/or temporal lags.

The OI method produces an objective estimate of the expected error in the result.

The OI technique can be formulated to simultaneously interpolate different data types (e.g. winds and geopotential heights) provided a linear relationship between the model and data (e,g, as in the case of geostrophic winds (currents) computed from geopotential (sea surface) heights).

In the case of geostrophic turbulence, the assumption of isotropy dictates a fixed relationship between the covariance of the individual velocity components and streamfunction.

Simultaneously estimating multiple variables has the advantage that known physical constraints (e.g. continuity, geostrophhy) can be incorporated into the mapping procedure, thereby producing results that are “balanced” kinematically and/or dynamically.

Examples: Combining altimeter sea surface height observations and velocity observations from sequential satellite imagery.

Wilkin, J. L., M. M. Bowen and W. J. Emery (2002), Mapping mesoscale currents by optimal interpolation of satellite radiometer and altimeter data, Ocean Dynamics, 52, 95-103. (pdf).

John’s notes for lecture 4

April 5:

Optimal interpolation example

The matlab scripts cov_mercator.m and oi_mercator.m in the class matlab directory demonstrate fitting a covariance function to a set of synthetic data, and using this function to optimally interpolate to a regular grid.

The data used in this example is ocean temperature taken from the French Mercator http://www.mercator-ocean.fr/en/ operational ocean forecast system for the North Atlantic.

cov_mercator:

The script cov_mercator.m loads the example Mercator snapshot form a mat file, subsamples the data to a small (3%) subset and adds seom normally distributed noise to emulate instrument error (or unresolved high frequency physical variability due e.g. to internal waves in the case of in situ ocean temperature observations).
The lon/lat coordinates of the sub-sampled “data” set are converted to simple x,y, coordinates w.r.t. the southwest corner of the data range, and then the separation distance between all data (r) so that a binned-lagged covariance as a function of r can be computed from the data themselves, i.e. C(r) = <d(x)d(x+r)>.
Two functional forms (Gaussian and Markov) for C(r) are fitted to the estimated covariance using Matlab’s fminsearch function. Note that the covariance at r=0 is not used in the fit because this includes the effect of the independent observational, or error, variance.
The apparent error variance is calculated fomr the difference of the data variance at r=0 (i.e. var(data)) at the estimated signal variance at r=0 that is indicated by the y-intercept of the functional fit, i.e. C(0).

oi_mercator:

Using the fitted Markov covariance function, normalized model-data (C_md) and data-data (C_dd) covariance matrices are computed. The data-data is augmented on the diagonal with the ratio of error to signal variance.
The optimal interpolation fit to the data is computed by direct inversion of C_dd and also by the Matlab matrix-left-divide operation to compute the product C_dd d directly.
Expected errors of the OI are calculated, and qualitatively compared to the actual residuals of the fit to show that, as expected, approximately 65% of the residuals fall within the expected errors.

If the data set (N points) is large, the matrix sizes may get too large for practical handling in Matlab (or any other language) because the OI problem has to solve matrix inversions or simultaneous equation solutions of dimension NxN.

This may be handled by dividing the model grid into subdomains, like tiles, with subsets of the data limited to only those points that fall within the model ‘tile’ plus a ‘halo’ region around the tiles. The halo region should be at least one covariance scale wide to ensure smoothness at the tile boundaries. The data-data matrix Cdd will have to be computed anew for each tile, but computing e.g. order(10) OI operations for order(N/10) data elements may be faster than one OI for order(N) data elements. This is because the computational efforts of the matrix operations scale with N³.

If there are too many data within a few covariance scales to practically invert Cdd, then it is likely that there are more data than necessary to resolve the mapped field. This indicates that a shorter covariance scale can probably be used. Alternatively, it is probably safe to decimate the data (just use less of it, thereby making N smaller) or average the observations in small bins. For independent errors, the binning step will reduce the expected error (the noise variance) of the binned values, and this information can be carried through the analysis.

Expected errors

A posteriori testing of error estimates can be done to see whether the proportion of residuals within the expected errors is statistically consistent. More in depth tests would compare the results to a set of independent data, such as from another instrument, or data withheld from the OI itself. See Walker and Wilkin (1998) for an example of checking the validity of error estimates.

Related topics: Kriging. See references and web resources below.

No notes for lecture 5. Run the cov_mercator and oi_mercator Matlab scripts.

April 8: Empirical Orthogonal Functions

The terminology Principal Components (PC) and Empirical Orthogonal Functions (EOF) are generally used interchangeably in the earth sciences. Oceanographers almost always refer to EOFs, whereas meteorologists and climatologists mix both terms. Some people make a distinction between PC and EOF on the basis of whether the

Eigenvector theory

In its simplest formulation, EOFs are eigenvectors of the data covariance matrix. The eigenvectors are commonly referred to as modes or loadings.

The eigenvectors of the ‘time’ (column dimension) and ‘space’ (row dimension) coordinates are complimentary. They have the same eigenvalues. It can be shown that the eigenvalues and both sets of eigenvectors are inter-related through the singular value decomposition of the data matrix.

EOFs can be obtained by solving for the eigenvalues of the covariance matrix, or by Singular Value Decomposition of the data matrix thereby effectively by-passing computation of the covariance matrix. The choice is an issue of computational convenience. Issues are the size of the matrix (M>N or M<N) and the difficulties associated with dealing with missing data.

John’s notes for lecture 6

The 2-dimensional eigenvalue problem solved in class can be viewed by running twodimeig.m in Matlab.

April 12: Empirical Orthogonal Functions – solution by singular value decomposition

John’s notes for lecture 7

Matlab script cacsst.m demonstrates calculating EOFs from the data matrix by singular value decomposition.

The data used are CAC monthly observed SST obtained from the National Virtual Ocean Data System (NVODS) Live Access Server at NOAA PMEL. The following link goes straight to the data set selection and subset utility:

http://ferret.pmel.noaa.gov/NVODS/servlets/dataset?catitem=5890

April 15: Empirical Orthogonal Functions

Expressing EOFs with dimensional units. Skill scores.

Estimating significance levels for retaining EOFs

Qualitative: the “scree” plot method
Quantitative: Use Monte Carlo method to estimate the 95% confidence limits at which the fraction of variance of the data in each mode is significant.
See Matlab script cacsstsignificance.m
See the discussion in Jollife (2002) chapter 6 for more details.
North’s Rule of Thumb

John’s notes for lecture 8

April 19:

Processes that are not standing wave in nature

If lag correlation of amplitude time series shows phase quadrature in time then there are possibly propagating waves. In this situation, conventional EOF analysis tedns to partition the propagating wave into two standing waves of similar amplitude. When examining fractional variance eigenvalues the possibility of this occurrence should be noted. Extensions to the EOF method designed to identify propagating waves are the Hilbert EOF (often called Complex EOF) method.

How to deal with gaps?

Ignore and just compute covariance matrix elements based on the common points in each time series.
Introduce randomly generated data to fill gaps and test for M realizations
Fill the gaps in each data series using e.g. optimal interpolation or weighted regression.

Rotated EOFs

Example of varimax rotation of Pacififc SST EOFs. (Richman)

See example script rotationexample.m and the varimax algorithm it uses varimax_kaplan.m obtained from http://erizo.ucdavis.edu/~dmk/notes/EOFs/EOFs.html

John’s notes for lecture 9

Topics I never got to:

Robustness of the analysis – things to look for:

Are the patterns sensitive to the choice of domain? If you partition that domain into separate regions, do you get the same loadings in each? Are the patterns sensitive to the sample of data used? If you split the time series in two, do you get the same spatial pattern for the first and seconds halves of the data set?

Putting it together: using OI and EOFs together to map 3D datasets

Using EOFs to map 3-D datasets (see papers I. Fukumori and C. Wunsch, N. Bindoff, and H. Mercier). Steps are compute vertical modes of variability as the loadings of a sample of hydrographic station data, and then optimal interpolation of the map of spatial loadings to a regular grid. Reconstruct the full climatology by using the vertical modes to project the horizontally mapped amplitudes to all depths.

Inferring patterns from correlations (subsurface projection).

Subsurface ocean density anomalies can be inferred by noting that the amplitudes of the gravest EOFs of vertical anomalies correlate with direct observations of SSH from satellites. A regression of satellite observed SSH anomalies on in situ dynamics height anomalies from CTD data gives the factor by which to multiply vertical EOF structures to get the subsurface anomaly estimate.

Complex EOFs (for wave phenomena)

Cyclostationary EOF (for periodic but not simple harmonic processes e.g. Quasi-Biennial Oscillation)

Canonical Correlation Analysis

If variables are expected to be closely related in variability, but with different spatial patterns (e.g. ocean temperature and wind speed), they cannot be combined into am multi-variate EOF. CCA offers ways to identify linked variation in time that has distinct spatial patterns

Wind data to complement the example SST data set:

FSU Pacific Monthly Pseudostress Averages

http://ferret.pmel.noaa.gov/NVODS/servlets/constrain?var=18352

Schedule of student presentations

April 22:

Thomas Atkins
Chaochao Gao
Alex Kahl
Natasa Skific

April 26:

Stephanie Weber
Gordon Zhang
Keleigh McAllister
Mike Previdi

April 29:

Meenal Gogate
Mike Foster
Haibin Li
Yonghua Chen

May 3:

Jaci Secora
John Krasting
Sindia Sosdian
Gregg Foti
Jung Hoon Kim

References and study resources

References

Daley, R., 1991, Atmospheric Data Analysis, Cambridge University Press, 457 pp.

Denman, K.L. and H.J. Freeland. 1985. Correlation scales, objective mapping and a statistical test of geostrophy over the continental shelf. J. Mar. Res. 43(3), 517-539.

Holbrook, N. J. and N. L. Bindoff, 2000:A statistically efficient mapping technique for four-dimensional ocean temperature data. J. Atmos. Ocean. Technol., 17, 831-846.
[pdf]

J.D. Horel. 1984. Complex principal component analysis: theory and examples. Journal of Climate and Applied Meteorology , vol. 23,no. 12 , pp. 1660-73 , Dec. 1984.

Jollife, I. T., Principal Component Analysis, Springer Series in Statistics, Springer-Verlag, New York, 487 pp., 2^nd ed. 2002.

Jolliffe, I. T. and M.B. Richman. 1987. Rotation of principal components: some comments (with reply). Journal of Climatology, vol. 7, no. 5, pp. 507-20.

Kelly, K., 1985, The influence of winds and topography on he sea surface temperature patterns over the northern California slope, J. Geophys. Res., 90, 11,783-11,798.

Merrifield. M. and R.T. Guza, 1990. Detecting propagating signals with complex empirical orthogonal functions: A cautionary note. Journal of Physical Oceanography, vol. 20, no. 10 , pp. 1628-33, Oct. 1990.

Mestas-Nunez, A.M, 2000: Orthogonality properties of rotated empirical modes. International Journal of Climatology, 20:1509-1516. [pdf]

Richman, M. B., 1986: Rotation of principal components. Journal of Climatology, vol. 6, no. 3, pp. 293-335.

Ridgway, K., J. R. Dunn and J. L. Wilkin, 2002, Ocean interpolation by 4-dimensional weighted least squares: Application to the waters around Australasia, Journal of Atmospheric and Oceanic Technology, 19, 1357-1375. (pdf)

Strang, G., 1980, Linear Algebra and Its Applications, Academic Press.

Walker, A.E. and J.L. Wilkin, 1998, Optimal averaging of NOAA/NASA Pathfinder satellite sea surface temperature data, Journal of Geophysical Research, 103, 12,869-12,883. (pdf)

Wilkin, J. L., M. M. Bowen and W. J. Emery, 2002, Mapping mesoscale currents by optimal interpolation of satellite radiometer and altimeter data, Ocean Dynamics, 52, 95-103. (pdf).

Wunsch, C., The Ocean Circulation Inverse Problem, Cambridge University Press, 442 pp., 1996.

http://www.atmos.washington.edu/~dennis/552_Notes_4.pdf

Varimax rotation algorithm for Matlab:

http://erizo.ucdavis.edu/~dmk/notes/EOFs/EOFs.html