From GHER

(Difference between revisions)

Revision as of 09:40, 7 June 2013

divand: n-dimensional variational data analysis

divand performs an n-dimensional variational analysis of arbitrarily located observations. Observation will be interpolated on a curvilinear grid in 2, 3 or more dimensions.

Requirements

Your need GNU Octave (version 3.4 or later) or MATLAB to run divand.

Installing

Download the source code from http://modb.oce.ulg.ac.be/mediawiki/upload/Alex/divand/divand-1.0.tar.gz.
Extract the package
In MATLAB add the directory divand-x.y to the search path. In octave add the command to your .octaverc file (for testing it is sufficient if you move your current working directory to the divand-x.y directory):

addpath('/path/to/divand-x.y.');

Testing

A test script test_divand is included to verify the correct functioning of the toolbox. The script should be run in a Octave or MATLAB session. All tests should pass without warning or error.

run test_interp_1d: (max difference=2.22045e-16)   OK  
run test_interp_2d: (max difference=2.22045e-16)   OK  
run test_interp_regular: (max difference=2.22045e-16)   OK  
run test_sparse_diff: (max difference=0)   OK  
run test_1dvar:   OK  
run test_2dvar: (max difference=4.28659e-16)   OK  
run test_2dvar_check: (max rms=2.39788e-16)   OK  
run test_2dvar_adv:   OK  
run test_2dvar_rellen:   OK  
run test_2dvar_lenxy:   OK  
run test_2dvar_check_correrr: (max rms=7.45461e-16)   OK  
run test_2dvar_constrain: (max rms=9.41909e-10)   OK  
run test_2dvar_cyclic: (max rms=2.71352e-14)   OK  
run test_2dvar_eof_var:   OK  
run test_2dvar_eof_check:   OK  
run test_3dvar: (max difference=0.0182573)   OK  
run test_3dvar_large_stacked: (max difference=2.40796e-15)   OK  
run test_4dvar: (max difference=0.0111807)   OK

Documentation

Issue in MATLAB/Octave the command help divand:

 [fi,err,s] = divand(mask,pmn,xi,x,f,len,lambda,...); 
 [fi,err] = divand(mask,pmn,xi,x,f,len,lambda,...); 
 [fi,s] = divand(mask,pmn,xi,x,f,len,lambda,...); 
 
 Perform an n-dimentional variational analysis of ariabitrarly located observations.

Input

mask: binary mask delimiting the domain. 1 is inside and 0 outside. For oceanographic application, this is the land-sea mask.
pmn: scale factor of the grid. pmn is a cell array with n elements. Every element represents the scale factor of the corresponding dimension. Its inverse is the local resolution of the grid in a particular dimension.
xi: cell array with n elements. Every element represents a coordinate of the final grid on which the observations are interpolated
x: cell array with n elements. Every element represents a coordinate of the observations
f: value of the observations *minus* the background estimate (m-by-1 array). (see note)
len: correlation length

if len is scalar, then the correlation length is the same in all directions and everywhere in the domain
if len is cell array of n scalars, then the correlation length can depend on the direction but for a given direction the correlation length is the same everywhere in the domain. The option is for example used when analysing 3D domain (longitude, latitude and depth) since the vertical scale is much smaller than the horizontal scale.
if len is cell array of n arrays, then the correlation length can depend on the direction and on the location. The individual arrays should have the same size as the domain.

lambda: signal-to-noise ratio of observations (if lambda is a scalar). The larger this value is, the closer is the field fi to the observation.

If lambda is a scalar, then R = 1/lambda I, where R is the observation error covariance matrix,
If lambda is a vector, then R = diag(lambda)
If lambda is a matrix (or a CovarParam object), then R = lambda

Optional input arguments

Optional input arguments specified as pairs of keyword and values:

'velocity', vel: velocity of advection constraint. The default is no-advection constraint

'alpha': alpha is vector of coefficients multiplying various terms in the cost function.

The first element multiplies the norm. The other i-th element of alpha multiplies the (i+1)-th derivative. Per default, the highest derivative is m = ceil(1+n/2) where n is the dimension of the problem. The values of alpha are per default the (m+1)th row of the Pascal triangle:

           m=0         1
           m=1       1   1
           m=1     1   2   1     (n=1,2)
           m=2   1   3   3   1   (n=3,4)
           ...

'diagnostics': 0 or 1 turns diagnostic and debugging information on (1) or off (0, default). If on, they will be returned as the last output argument

'constraint': a structure with user specified constraint

'moddim': modulo for cyclic dimension (vector with n elements). Zero is used for non-cyclic dimensions. Halo points should not be included for cyclic dimensions. For example if the first dimension is cyclic, then the grid point corresponding to mask(1,j) should be between mask(end,1) (left neighbord) and mask(2,j) (right neighbord)

'fracdim': fractional indices (n-by-m array). If this array is specified, then x and xi are not used.

Note: 'velocity' and 'constraint' may appear multiple times

Output

fi: the analysed field
err: error variance of the analysis field relative to the error variance of the background
s: structure with some diagnostics mainly for debugging purposes.

Note

If zero is not a valid first guess for your variable (as it is the case for e.g. ocean temperature), you have to subtract the first guess from the observations before calling divand and then add the first guess back in.

Example

see divand_simple_example.m

% A simple example of divand is 2 dimensions
% with observations from an analytical function


% observations
x = rand(75,1);
y = rand(75,1);
f = sin(x*6) .* cos(y*6);

% final grid
[xi,yi] = ndgrid(linspace(0,1,30));

% reference field
fref = sin(xi*6) .* cos(yi*6);

% all points are valid points
mask = ones(size(xi));

% this problem has a simple cartesian metric
% pm is the inverse of the resolution along the 1st dimension
% pn is the inverse of the resolution along the 2nd dimension

pm = ones(size(xi)) / (xi(2,1)-xi(1,1));
pn = ones(size(xi)) / (yi(1,2)-yi(1,1));
  
% correlation length
len = 0.1;

% signal-to-noise ratio
lambda = 20;

% fi is the interpolated field
fi = divand(mask,{pm,pn},{xi,yi},{x,y},f,len,lambda);

% plotting of results
subplot(1,2,1); 
pcolor(xi,yi,fref);
shading flat,colorbar
ca = caxis;
set(findobj(gcf,'FontSize',10),'FontSize',16)
title('Reference field and observation loc.');
hold on
plot(x,y,'k.');

subplot(1,2,2); 
hold on
pcolor(xi,yi,fi);
shading flat,colorbar
caxis(ca);
set(findobj(gcf,'FontSize',10),'FontSize',16)
title('Interpolated field');

set(gcf,'PaperUnits','centimeters','PaperPosition',[0 0 25 12 ]);
print -dpng divand_simple_example.png

Fun

A educational web application has been developed to reconstruct a field based on point "observations". The user must choose in an optimal way the location of 10 observation such that the analysed field obtained by divand based on these observation is as close as possible to the original field.

Troubleshooting

Cholesky factorization fails with the following error:

??? Error using ==> chol
CHOLMOD internal error.

Error in ==> CovarIS.factorize at 5
  [self.RP, self.q, self.PP] = chol (self.IS);

This might happen if:

memory is insufficient,
matrix contains a NaN introduced through metric scale factor equal to zero or a user constraint.

@@ Line 174: / Line 174: @@
 = Troubleshooting =
-Choleky factorization fails with the following error:
+Cholesky factorization fails with the following error:
   ??? Error using ==> chol

Divand