%% |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!  
% 
% See also: <filt1_documentation.html |filt1|>. 
% 
% <CDT_Contents.html Back to Climate Data Tools Contents>.
%% 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_documentation.html |cdtgrid|>. 
% 
%% Author Info
% This function was written by <http://www.chadagreene.com 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 <http://www.mathworks.com/matlabcentral/fileexchange/20417 |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.  