Basic Filtering#

The Filter Object#

The core object in GCM-Filters is the gcm_filters.Filter object. Its full documentation below enumerates all possible options.

gcm_filters.Filter(filter_scale: float, dx_min: float, filter_shape: gcm_filters.filter.FilterShape = FilterShape.GAUSSIAN, transition_width: float = 3.141592653589793, ndim: int = 2, n_steps: int = 0, grid_type: gcm_filters.kernels.GridType = GridType.REGULAR, grid_vars: dict = <factory>) None[source]#

A class for applying diffusion-based smoothing filters to gridded data.

Parameters
  • filter_scale (float) – The filter scale, which has different meaning depending on filter shape

  • dx_min (float) – The smallest grid spacing. Should have same units as filter_scale

  • n_steps (int, optional) – Number of total steps in the filter n_steps == 0 means the number of steps is chosen automatically

  • filter_shape (FilterShape) –

    • FilterShape.GAUSSIAN: The target filter has shape \(e^{-(k filter_scale)^2/24}\)

    • FilterShape.TAPER: The target filter has target grid scale Lf. Smaller scales are zeroed out. Scales larger than pi * filter_scale / 2 are left as-is. In between is a smooth transition.

  • transition_width (float, optional) – Width of the transition region in the “Taper” filter. This is a nondimensional parameter. Theoretical minimum is 1; not recommended.

  • ndim (int, optional) – Laplacian is applied on a grid of dimension ndim

  • grid_type (GridType) – what sort of grid we are dealing with

  • grid_vars (dict) – dictionary of extra parameters used to initialize the grid Laplacian

gcm_filters.filter_spec#
Type

FilterSpec

Details related to filter_scale, filter_shape, transition_width, and n_steps can be found in the Filter Theory. The following sections explain the options for grid_type and grid_vars in more detail.

Grid types#

To define a filter, we need to pick a grid and associated Laplacian that matches our data. The currently implemented grid types are:

In [1]: import gcm_filters

In [2]: list(gcm_filters.GridType)
Out[2]: 
[<GridType.REGULAR: 1>,
 <GridType.REGULAR_AREA_WEIGHTED: 2>,
 <GridType.REGULAR_WITH_LAND: 3>,
 <GridType.REGULAR_WITH_LAND_AREA_WEIGHTED: 4>,
 <GridType.IRREGULAR_WITH_LAND: 5>,
 <GridType.MOM5U: 6>,
 <GridType.MOM5T: 7>,
 <GridType.TRIPOLAR_REGULAR_WITH_LAND_AREA_WEIGHTED: 8>,
 <GridType.TRIPOLAR_POP_WITH_LAND: 9>,
 <GridType.VECTOR_C_GRID: 10>]

This list will grow as we implement more Laplacians.

The following table provides an overview of these different grid type options: what grid they are suitable for, whether they handle land (i.e., continental boundaries), what boundary condition the Laplacian operators use, and whether they come with a scalar or vector Laplacian. You can also find links to example usages.

GridType

Grid

Handles land

Boundary condition

Laplacian type

Example

REGULAR

Cartesian grid

no

periodic

Scalar Laplacian

REGULAR_WITH_LAND

Cartesian grid

yes

periodic

Scalar Laplacian

see below

IRREGULAR_WITH_LAND

locally orthogonal grid

yes

periodic

Scalar Laplacian

Example: Different filter types; Example: Filtering on a tripole grid

MOM5U

Velocity-point on Arakawa B-Grid

yes

periodic

Scalar Laplacian

MOM5T

Tracer-point on Arakawa B-Grid

yes

periodic

Scalar Laplacian

TRIPOLAR_POP_WITH_LAND

locally orthogonal grid

yes

tripole

Scalar Laplacian

Example: Filtering on a tripole grid

VECTOR_C_GRID

Arakawa C-Grid

yes

periodic

Vector Laplacian

Example: Viscosity-based filtering with vector Laplacians

Grid types with scalar Laplacians can be used for filtering scalar fields (such as temperature), and grid types with vector Laplacians can be used for filtering vector fields (such as velocity).

Grid types for simple fixed factor filtering#

The remaining grid types are for a special type of filtering: simple fixed factor filtering to achieve a fixed coarsening factor (see also the Filter Theory). If you specify one of the following grid types for your data, gcm_filters will internally transform your original (locally orthogonal) grid to a uniform Cartesian grid with dx = dy = 1, and perform fixed factor filtering on the uniform grid. After this is done, gcm_filters transforms the filtered field back to your original grid. In practice, this coordinate transformation is achieved by area weighting and deweighting (see Filter Theory). This is why the following grid types have the suffix AREA_WEIGHTED.

GridType

Grid

Handles land

Boundary condition

Laplacian type

Example

REGULAR_AREA_WEIGHTED

locally orthogonal grid

no

periodic

Scalar Laplacian

REGULAR_WITH_LAND_AREA_WEIGHTED

locally orthogonal grid

yes

periodic

Scalar Laplacian

Example: Different filter types; Example: Filtering on a tripole grid

TRIPOLAR_REGULAR_WITH_LAND_AREA_WEIGHTED

locally orthogonal grid

yes

tripole

Scalar Laplacian

Example: Filtering on a tripole grid

Grid variables#

Each grid type from the above two tables has different grid variables that must be provided as Xarray DataArrays. For example, let’s assume we are on a Cartesian grid (with uniform grid spacing equal to 1), and we want to use the grid type REGULAR_WITH_LAND. To find out what the required grid variables for this grid type are, we can use this utility function.

In [3]: gcm_filters.required_grid_vars(gcm_filters.GridType.REGULAR_WITH_LAND)
Out[3]: ['wet_mask']

wet_mask is a binary array representing the topography on our grid. Here the convention is that the array is 1 in the ocean (“wet points”) and 0 on land (“dry points”).

In [4]: import numpy as np

In [5]: import xarray as xr

In [6]: ny, nx = (128, 256)

In [7]: mask_data = np.ones((ny, nx))

In [8]: mask_data[(ny // 4):(3 * ny // 4), (nx // 4):(3 * nx // 4)] = 0

In [9]: wet_mask = xr.DataArray(mask_data, dims=['y', 'x'])
In [10]: wet_mask.plot()
Out[10]: <matplotlib.collections.QuadMesh at 0x7f1791a1f0a0>
_images/wet_mask.png

We have made a big island.

Note

Some more complicated grid types require more grid variables. The units for these variables should be consistent, but no specific system of units is required. For example, if grid cell edge lengths are defined using kilometers, then the filter scale and dx_min should also be defined using kilometers, and the grid cell areas should be defined in square kilometers.

Creating the Filter Object#

We create a filter object as follows.

In [11]: filter = gcm_filters.Filter(
   ....:     filter_scale=4,
   ....:     dx_min=1,
   ....:     filter_shape=gcm_filters.FilterShape.TAPER,
   ....:     grid_type=gcm_filters.GridType.REGULAR_WITH_LAND,
   ....:     grid_vars={'wet_mask': wet_mask}
   ....: )
   ....: 

In [12]: filter
Out[12]: Filter(filter_scale=4, dx_min=1, filter_shape=<FilterShape.TAPER: 2>, transition_width=3.141592653589793, ndim=2, n_steps=16, grid_type=<GridType.REGULAR_WITH_LAND: 3>)

The string representation for the filter object in the last line includes some of the parameters it was initiliazed with, to help us keep track of what we are doing. We have created a Taper filter that will filter our data by a fixed factor of 4.

Applying the Filter#

We can now apply the filter object that we created above to some data. Let’s create a random 3D cube of data that matches our grid.

In [13]: nt = 10

In [14]: data = np.random.rand(nt, ny, nx)

In [15]: da = xr.DataArray(data, dims=['time', 'y', 'x'])

In [16]: da
Out[16]: 
<xarray.DataArray (time: 10, y: 128, x: 256)>
array([[[6.79549819e-01, 2.56747257e-01, 1.90721202e-01, ...,
         4.23222276e-01, 8.80245446e-01, 6.70726214e-01],
        [2.66900790e-01, 1.53166817e-01, 3.32715806e-01, ...,
         6.76637783e-01, 4.19362420e-01, 7.51364233e-01],
        [3.44317934e-01, 2.10520784e-01, 1.58466537e-01, ...,
         3.89028642e-01, 9.70924763e-01, 5.49384944e-01],
        ...,
        [7.39874157e-01, 5.26131831e-01, 5.42453249e-01, ...,
         6.13632715e-01, 3.92399522e-01, 6.81270717e-01],
        [1.63759161e-01, 1.36273087e-01, 3.44976852e-01, ...,
         2.00889166e-01, 5.41555982e-01, 8.57385411e-01],
        [3.49960990e-01, 2.52985802e-01, 5.78762990e-01, ...,
         8.51936868e-01, 7.00822717e-01, 6.13251856e-01]],

       [[1.59983362e-01, 1.15682581e-01, 4.08507300e-01, ...,
         1.98882122e-01, 1.49109785e-01, 2.51243938e-01],
        [9.49797992e-01, 1.74470075e-01, 4.57082532e-01, ...,
         7.76905072e-02, 5.70286134e-01, 1.65427115e-01],
        [1.61425635e-01, 4.95416678e-01, 4.08463855e-01, ...,
         7.80769068e-01, 5.19746017e-01, 9.74830416e-01],
...
        [8.34564643e-01, 9.57975921e-01, 5.45621494e-01, ...,
         6.90327446e-01, 1.71680422e-01, 8.75053635e-01],
        [8.59290463e-01, 7.97925211e-02, 2.03877795e-01, ...,
         6.72224069e-01, 1.89766030e-01, 6.23632708e-01],
        [9.21238407e-01, 5.94241891e-01, 7.95765022e-01, ...,
         7.24370070e-01, 6.24748656e-01, 4.26456127e-01]],

       [[2.16166333e-01, 7.37701912e-01, 7.96036394e-01, ...,
         2.09145433e-01, 5.11066799e-01, 9.87273989e-01],
        [2.60280974e-02, 9.84357153e-01, 1.73164069e-01, ...,
         4.42358290e-01, 4.28379257e-01, 7.43537312e-02],
        [3.58534288e-01, 4.98661804e-01, 6.78755499e-01, ...,
         1.27274231e-01, 1.00086670e-01, 1.44241363e-01],
        ...,
        [5.48673465e-01, 7.73927589e-01, 5.55437022e-01, ...,
         3.63146238e-01, 1.55075435e-01, 8.48222897e-02],
        [6.63697602e-01, 6.22412802e-03, 5.55260112e-01, ...,
         3.03529269e-01, 2.70588937e-03, 8.54114031e-01],
        [6.36990023e-01, 9.90165830e-01, 7.93079363e-01, ...,
         4.24575998e-01, 1.25311318e-01, 5.41165209e-01]]])
Dimensions without coordinates: time, y, x

We now mask our data with the wet_mask.

In [17]: da_masked = da.where(wet_mask)
In [18]: da_masked.isel(time=0).plot()
Out[18]: <matplotlib.collections.QuadMesh at 0x7f17918feca0>
_images/data.png

Now that we have some data, we can apply our filter. We need to specify which dimension names to apply the filter over. In this case, it is y, x.

Warning

The dimension order matters! Since some filters deal with anisotropic grids, the latitude / y dimension must appear first in order to obtain the correct result. That is not an issue for this simple (isotropic) toy example but needs to be kept in mind for applications on real GCM grids.

In [19]: %time da_filtered = filter.apply(da_masked, dims=['y', 'x'])
CPU times: user 104 ms, sys: 20 ms, total: 124 ms
Wall time: 124 ms
In [20]: da_filtered
Out[20]: 
<xarray.DataArray (time: 10, y: 128, x: 256)>
array([[[0.43413226, 0.33930291, 0.3401849 , ..., 0.65326858,
         0.64184421, 0.56156196],
        [0.41568458, 0.32833427, 0.34727813, ..., 0.66847332,
         0.64937917, 0.55156423],
        [0.40112062, 0.3370507 , 0.36826005, ..., 0.62983306,
         0.6094977 , 0.51796927],
        ...,
        [0.52411866, 0.49092206, 0.4834373 , ..., 0.52316595,
         0.54588984, 0.55287247],
        [0.48920299, 0.43349416, 0.42068515, ..., 0.54024529,
         0.55381381, 0.54288957],
        [0.45715121, 0.37535553, 0.36418936, ..., 0.59629554,
         0.59711076, 0.55120981]],

       [[0.35450707, 0.34523498, 0.35364483, ..., 0.31850903,
         0.34463275, 0.36117314],
        [0.40291918, 0.38665239, 0.38660331, ..., 0.37358883,
         0.40202735, 0.41493939],
        [0.48096979, 0.45777828, 0.45487027, ..., 0.4571414 ,
         0.49517675, 0.50290719],
...
        [0.52931764, 0.52521728, 0.53425387, ..., 0.56520506,
         0.56429544, 0.54668307],
        [0.53115583, 0.53547974, 0.55154455, ..., 0.53008455,
         0.54181649, 0.53710865],
        [0.47967983, 0.49013693, 0.52328617, ..., 0.53048715,
         0.51789828, 0.4932214 ]],

       [[0.50523046, 0.60479949, 0.66656192, ..., 0.40545442,
         0.38875615, 0.42108813],
        [0.43001846, 0.55591293, 0.65567449, ..., 0.39281355,
         0.34169526, 0.34738571],
        [0.38960964, 0.54221809, 0.67817483, ..., 0.36449229,
         0.30169467, 0.29967612],
        ...,
        [0.4388221 , 0.50009338, 0.54172735, ..., 0.31925286,
         0.33280214, 0.37580614],
        [0.50142614, 0.57938142, 0.62230636, ..., 0.32142382,
         0.34947691, 0.41463169],
        [0.53877662, 0.62465022, 0.66808914, ..., 0.37214112,
         0.38872119, 0.44823388]]])
Dimensions without coordinates: time, y, x

Let’s visualize what the filter did.

In [21]: da_filtered.isel(time=0).plot()
Out[21]: <matplotlib.collections.QuadMesh at 0x7f179187a2b0>
_images/data_filtered.png

Using Dask#

Up to now, we have filtered eagerly; when we called .apply, the results were computed immediately and stored in memory. GCM-Filters is also designed to work seamlessly with Dask array inputs. With dask, we can filter lazily, deferring the filter computations and possibly executing them in parallel. We can do this with our synthetic data by converting them to dask.

In [22]: da_dask = da_masked.chunk({'time': 2})

In [23]: da_dask
Out[23]: 
<xarray.DataArray (time: 10, y: 128, x: 256)>
dask.array<xarray-<this-array>, shape=(10, 128, 256), dtype=float64, chunksize=(2, 128, 256), chunktype=numpy.ndarray>
Dimensions without coordinates: time, y, x

We now filter our data lazily.

In [24]: da_filtered_lazy = filter.apply(da_dask, dims=['y', 'x'])

In [25]: da_filtered_lazy
Out[25]: 
<xarray.DataArray (time: 10, y: 128, x: 256)>
dask.array<transpose, shape=(10, 128, 256), dtype=float64, chunksize=(2, 128, 256), chunktype=numpy.ndarray>
Dimensions without coordinates: time, y, x

Nothing has actually been computed yet. We can trigger computation as follows:

In [26]: %time da_filtered_computed = da_filtered_lazy.compute()
CPU times: user 203 ms, sys: 11.9 ms, total: 215 ms
Wall time: 128 ms

Here we got only a very modest speedup because our example data are too small. For bigger data, the performance benefit will be more evident.