Two-dimensional spatial filters, developed for application to ocean and atmosphere model output and gridded observations
The functions in this repository apply various two-dimensional spatial filters to gridded data. These data can be output from a model simulation, or observations in a gridded layout (e.g., from satellite-based products, reanalyses, Argo). There is no requirement that the grid points be evenly spaced or aligned with longitude or latitude lines; however, the filters will work better if any curvature in the grid occurs on scales larger than the filter threshold scale. For input fields being filtered, non-finite (NaN, -Inf/Inf) or zero values are assumed to be land areas, and buffering and/or adjustments are implemented to reduce (but not entirely eliminate) variance loss at the land boundaries. Currently these functions are available in Matlab (should also work in GNU Octave) and Python.
The most widely applicable filtering functions in this repository implement 2-D low-pass filters. The residual of applying these functions can also be used to isolate smaller scales, and these filters can also be combined to effectively "bandpass" for a specific range of spatial scales. The core filter in these functions is applied as a window in 2-D wavenumber space, after a planar fit to the current tile has been removed and the tile is nested in a larger 2x2 tile of zeros to buffer the edges. The transfer function for the 2-D spectral coefficients is
where deg_threshold), steepness_factor representing how abrupt the cutoff is, and
The low-pass filters here should yield broadly similar results to those implemented in the Python package gcm-filters. The spatial-filters-2d low-pass filters are most similar to the TAPER filters in gcm-filters. Note that the definition of wavenumber/wavelength cutoff in gcm-filters is different from the half-power thresholds used as inputs in the functions here. Another key difference from gcm-filters is that spatial-filters-2d functions take longitude and latitude 2-D arrays as inputs and compute the relevant grid spacings automatically. However, if conservation of area- or transect-integrated quantities is desired in the low-pass filtered field, then it is necessary for the user to weight the input field with the relevant grid cell area or face length.
The low-pass filters in spatial-filters-2d take as input half-power wavelength thresholds defined in the deg_threshold vector (in units of degrees latitude) with a corresponding lat_deg_threshold vector indicating the latitude of each entry in deg_threshold. The threshold cutoff wavelengths are linearly interpolated between the latitudes defined in lat_deg_threshold. For example, if lat_deg_threshold = [-90 0 90] and deg_threshold = [2 20 2], then the half-power cutoff wavelength is 2 deg lat
Other aspects of the low-pass filtering process are described here:
In order to apply the filter at different latitudes with different threshold scales, the domain is first divided into a set of overlapping squares or tiles. There are two parameters ratio_side_deg_threshold and ratio_side_spacing that set the size of the tile squares and the spacing between them, respectively. ratio_side_deg_threshold has a recommended range between 2 and 4, with 2.5 or 3 usually an optimal choice. If ratio_side_deg_threshold = 2.5 and the local threshold wavelength is 10 degrees latitude, then the side of each square is 25 degrees latitude wide (and the equivalent distance in longitude). It is important that this value be somewhat larger than 1 (so that the threshold scale can be resolved) but not too large (otherwise you can get strange artifacts from the filtering especially near coastlines). ratio_side_spacing has a recommended range between 2 and 5, with 3 often an optimal choice. If ratio_side_spacing = 3 with squares of 25 degrees latitude on a side, then the centers of adjacent squares are separated by
To improve the filtered fields near coastlines (as determined by the presence of NaN, -Inf/Inf, or zero values), a coastal "buffer" is prepared. Before the final application of the filter, the field minus planar fit is filtered a number of times as specified by n_iter_coastbuffer, and after each application, the unfiltered values (minus planar fit) are restored to water grid cells. The effect is that after this process, the water grid cells are unchanged, but the land cells near the coastlines are buffered, so that the transitions in the low-pass filtered field are not abrupt near the coastlines, and attenuation of the filtered field near land areas is much reduced (though not entirely eliminated). Once the planar fit is restored, a new planar fit is computed and removed again prior to applying the filter one final time, with the filtered values over water grid cells retained. The fully low-passed field then has the planar fit restored, the original tile is recovered from the larger nested tile, and land areas are re-masked.
For the functions specifically intended to handle horizontal velocities (which contain vel in the name), an additional adjustment is applied. Namely, the low-pass filtered "wet" mask (with values of 1 for water or air grid cells and 0 for land) is used to compute a local angle for the low-passed coastline, and the coast-perpendicular component of the velocity field is set to zero before the final application of the filter. Note that the objective of this adjustment is not to eliminate low-passed flow across the real coastline, which would be impossible while retaining the scale properties of the filtered field. This correction merely reduces the low-passed cross-shore flow while retaining the low-passed alongshore component, a useful adjustment for improving large-scale transport representation in the presence of strong near-coastal flows (e.g., western boundary currents).
spat_filter_2d_lp.m: For filtering of most scalar quantities (including vertical velocity) in a non-periodic, locally Cartesian grid. This is best used for regional subsets of model output or data, though it could be applied to more complex global grids if the user manually stitches the connections of the grid together before calling the function. It can also be used to filter horizontal velocities if the user does not want the near-coastal velocity field adjustment described above.
spat_filter_vel_2d_lp.m: For filtering of horizontal velocity components
spatial_filters_2d_lp.py: A Python module containing the function spat_filter_2d_lp, equivalent to the Matlab function above. Can be imported into Python code or a Jupyter notebook by calling from spatial_filters_2d_lp import *.