grace | bundle

Tools for manipulating and analysing GRACE level-2 data

balaji devaraju
devaraju at ife uni-hannover de

Filtering GRACE data

Filtering is an envitable step in using GRACE data. Here we will use first the Gaussian filter to filter the data, then we will use the destriping filter and then finally we will use an unconventional filter, the spectral cosine filter

Gaussian filter

% Script for filtering GRACE data

% Read GRACE data
gsm = readgrace('/path/to/GRACE/data');

% Compute the GRACE mean
gsm_mean = getgracemean(gsm, [2005 1], [2010 12]);

% Remove mean
sc = cs2sc(gsm{10, 9} - gsm_mean);

% Maximum degree of spherical harmonic expansion
lmax = gsm{10, 8};       

% SH synthesis
qty = 'smd';    % Surface Mass Density
grd = 'block';  % Grid type
N   = 360;      % Grid size 0.5 degree grid
nrm = false;    % Remove normal field

% Unfiltered field
f = gshs_(sc, 'quant', qty, ...
              'grid', grd, ...
              'gridsize', N, ...
              'sub_WGS84', nrm);

% Plot the field
% Indicate whether the longitude of the field is from
% 0 to 360 (1) or -180 to 180 (0)
ldir = 1; 
h = mapfield(f, grd, ldir);

% Gaussian filter radius in radians
psi  = 4*pi/180; 

% Compute the Gaussian filter
Bg = gaussfltr(psi, lmax);

% Filter the field with the Gaussian filter
sc_fil = bsxfun(@times, sc, Bg);

% Filtered field
f_fil = gshs_(sc_fil, 'quant', qty, ...
                      'grid', grd, ...
                      'gridsize', N, ...
                      'sub_WGS84', nrm);

h_fil = mapfield(f_fil, grd, ldir);

Destriping filter

Now let us filter the GRACE data using the destriping filter

% Settings for the destriping filter
poly_degree = 2;
start_degree = 10;
start_order = 10;
end_degree = lmax;
end_degree = lmax;

% Destriping the GRACE data
sc_ds = dstrpngmtrx(sc, lmax, poly_degree, ...
                                start_order, end_order, ...
                                start_degree, end_degree);

Since destriping filter is a cascade filter you have to filter it again with Gaussian filter.

sc_ds_fil = bsxfun(@times, sc_ds, Bg);

% Destriped field
f_ds_fil = gshs_(sc_ds_fil, 'quant', qty, ...
                              'grid', grd, ...
                              'gridsize', N, ...
                              'sub_WGS84', nrm);

h_ds_fil = mapfield(f_ds_fil, grd, ldir);

Spectral cosine filter

Now let us use a filter that is seldom used in the GRACE community, the spectral cosine filter

% Parameters for the spectral cosine filter
lc = lmax;  % Degree after which you want to omit the signal
k  = 4;     % Cosine function raised to the power of k
lstart = 6; % We will not fit the first six degrees

% Compute the filter
Bc = spkcosine(lc, k, lstart, lmax);

% Filter the destriped coefficients
sc_ds_cos = bsxfun(@times, sc_ds, Bc);

% Destriped and spectral cosine filtered field
f_ds_cos = gshs_(sc_ds_cos, 'quant', qty, ...
                              'grid', grd, ...
                              'gridsize', N, ...
                              'sub_WGS84', nrm);

h_ds_cos = mapfield(f_ds_cos, grd, ldir);