bat365 Homepage

`filt2` documentation

The filt2 function performs a highpass, lowpass, bandpass, or bandstop 2D gaussian filter on gridded data such as topographic, atmospheric, oceanographic, or any kind of geospatial data. This function is designed to make it easy to remove features longer or shorter than a given characteristic wavelength. The input grid can contain NaNs!

Syntax

Zf = filt2(Z,res,lambda,filtertype)

Description

Zf = filt2(Z,res,lambda,filtertype) filters 2D dataset Z that has resolution res, to an approximate wavelength lambda. If the filtertype is 'lp' or 'hp' for lowpass or highpass, lambda must be a scalar value. If the filtertype is 'bp' or 'bs' for bandpass or bandstop, lambda must be a two-element array of the two cutoff wavelengths.

Explanation of this type of filter

There are many ways to filter a cat. The approach implemented by filt2 is less complex than a 2D FFT, but slightly more nuanced than a simple 2D moving average.

If you have a data matrix Z and you want to create a 2D lowpass-filtered matrix Zlow, a common way to do it is to fill each pixel in Zlow with the average value of all the nearby pixels in Z. You might specify a box size, say, 25 pixels by 25 pixels, and go down each row and column of Zlow, filling it with the average value of all the 25 by 25 surrounding pixels from Z. That approach effectively gets rid of the high-frequency noise, but it's a bit inelegant because a pixel 12 pixels away from the center of the box contributes just as much to the filtered center pixel value as the center pixel itself. Further, a moving square has contributions from corners, which are 1.4 times farther away from the center than the top or side pixels. It often makes more sense to use a weighting system that is symmetric about the center pixel, where pixels close to the center contribute more than the faraway pixels.

The 2D moving gaussian window implemented by filt2 is similar to a moving 2D square window, but the contribution of each pixel is weighted by a gaussian curve which is symmetric about the center pixel. So where the weighting of a moving square window would look like a box in 3D space, a gaussian window looks more like a rounded cone shape. Here's the weighting curve generated by filt2 if you give it a cutoff wavelength of 10 km and a for a 200 m resolution dataset:

lambda = 10; % cutoff "wavelength" in kilometers
res = 0.2;   % 200 m resolution of dataset

% Define standard deviation of gaussian curve:
sigma = (lambda/res) /(2*pi);

% Create the 2D gaussian curve:
f = fspecial('gaussian',2*ceil(2.6*sigma)+1,sigma);

surf(f)
axis tight
xlabel 'pixels'
ylabel 'pixels'
zlabel 'weight in moving average'

The value sigma here is about 8 pixels, at which point the weighting is about 0.6 times the amplitude of the center.

Example

Consider a 100 km by 100 km elevation dataset that has a resolution of 200 m. It has some long 25 km wavelength features aligned with north/south direction, some short ~5 km features oriented diagonally, and a considerable amount of random noise. There's also a block of missing data. Here's what your datset looks like:

res = 0.2; % 200 m resolution

x = 0:res:100; % eastings from 0 to 100 km
y = 0:res:100; % northings from 0 to 100 km
[X,Y] = meshgrid(x,y);

% Z contains 25 km features, ~5 km diagonal features, and noise:
Z = cos(2*pi*X/25)+cos(2*pi*(X+Y)/7)+randn(size(X));

% Z also has some missing data:
Z(100:120,100:120) = nan;

imagesc(x,y,Z);
axis xy image
caxis([-1 1])
title ' original Z matrix '
xlabel ' eastings (km) '
ylabel ' northings (km) '

Filter out all features that are shorter than about 15 km:

Zlow = filt2(Z,res,15,'lp');

imagesc(x,y,Zlow);
axis xy image
caxis([-1 1])
title ' 15 km lowpass filtered Z matrix '
xlabel ' eastings (km) '
ylabel ' northings (km) '

Similarly, filter out all features that are longer than 15 km:

Zhi = filt2(Z,res,15,'hp');
imagesc(x,y,Zhi);
axis xy image
caxis([-1 1])
title ' 15 km highpass filtered Z matrix '
xlabel ' eastings (km) '
ylabel ' northings (km) '

Or we can keep only the diagonal lines by removing the long 25 km wavelength features and removing the short-wavelength speckle noise: This is where you'll see it becomes a balancing act--tinker with the wavelengths and you'll realize the wider the bandpass region, the more noise gets through, but the more signal gets through too. Tighten up the bandpass region and you'll also be filtering out some of the 5 km wavelengths we want to keep because this filter does not have a sharp cutoff frequency:

Zbp = filt2(Z,res,[4 7],'bp');
imagesc(x,y,Zbp);
axis xy image
caxis([-1 1])
title ' 4 to 7 km bandpass filtered Z matrix '
xlabel ' eastings (km) '
ylabel ' northings (km) '

Similar to the bandpass filter, we can use a bandstop filter to get rid of the 5 km wavelength diagonal features:

Zbs = filt2(Z,res,[3 12],'bs');
imagesc(x,y,Zbs);
axis xy image
caxis([-1 1])
title ' 3 to 12 km bandstop filtered Z matrix '
xlabel ' eastings (km) '
ylabel ' northings (km) '

Special note on georeferenced data.

This function is great at filtering 2D spatial data, but do not use it on datasets in which latitudes and longitudes are equally spaced! That's because filtering to N degrees does not make any sense. Instead, this function is for use with data in projected coordinates, where grid dimensions are meter or kilometers, or even feet if you're so inclined. But do not use it on the type of grids that are created by cdtgrid.

Author Info

This function was written by Chad A. Greene of the University of Texas Institute for Geophysics in November 2016; however, all I did was repackage Carlos Adrian Vargas Aguilera's superb ndnanfilter function to make it a little easier to use for this specific application. Many thanks to Carlos for his well-thought-out code and clear documentation.

filt2 documentation