Spatial filters

Oceanostics provides spatial filters that operate directly on Oceananigans fields. All filters are built on top of Oceananigans' KernelFunctionOperation, so they compose with the rest of the Oceananigans ecosystem (outputs, reductions, other operations, etc.) and run on both CPU and GPU.

Box filter

The BoxFilter computes a local running-mean (box average) of a field over one or more grid directions. The stencil size is controlled by the N keyword: each filtered direction uses an N-point symmetric average centred on the current cell — i.e. N cells contribute to one filtered output value. N must be an odd integer ≥ 3, and refers to the size of the filter stencil, not the grid.

Multi-direction filters are fused into a single kernel at compile time, so a 3D box filter performs one pass over the data, not three.

Basic usage

Build a reusable filter once with the keyword-only constructor, then apply it to any field by calling it — bf(ψ) returns a KernelFunctionOperation that you wrap in a Field:

julia> using Oceananigans, Oceanostics

julia> grid = RectilinearGrid(size=(32, 32), x=(0, 1), z=(0, 1),
                              topology=(Periodic, Flat, Bounded));

julia> c = CenterField(grid);

julia> set!(c, (x, z) -> sin(2π * x) * z);

julia> bf = BoxFilter(; dims=(1, 3), N=5)
BoxFilter(dims=(1, 3), N=5, boundary=:shrink)

julia> c̄ = Field(bf(c));

julia> size(c̄)
(32, 1, 32)

The same bf can be applied to as many fields as you like. A one-step shortcut BoxFilter(ψ; dims, N, boundary), which builds and applies the filter in a single call, is also accepted throughout.

Boundary handling

For Periodic directions, stencil offsets are always wrapped — the boundary keyword is silently ignored. For Bounded directions the boundary keyword selects how out-of-bounds offsets are treated:

`boundary`	Behaviour
`:shrink` (default)	Drop out-of-bounds offsets from sum and count (honest local average; stencil shrinks near walls).
`:edge`	Replicate the nearest interior cell (`ψ[1]` or `ψ[N]`).
`(left=a, right=b)`	Pad with constant `a` on the low side and `b` on the high side.

A single spec applies to every filtered dimension, or a tuple gives per-dimension control:

julia> bf_edge  = BoxFilter(; dims=(1, 3), N=3, boundary=:edge);

julia> bf_mixed = BoxFilter(; dims=(1, 3), N=3, boundary=(:shrink, (left=0.0, right=0.0)));

julia> size(Field(bf_edge(c))) == size(Field(bf_mixed(c))) == (32, 1, 32)
true

API reference

Oceanostics.SpatialFilters.BoxFilter — Type

BoxFilter(ψ; dims, N, boundary)

One-step form: apply a box filter to ψ directly, returning the KernelFunctionOperation that computes its local box-average. Equivalent to BoxFilter(; dims, N, boundary)(ψ).

Refer to BoxFilter(; dims, N, boundary) for the full description of the keyword arguments and boundary handling.

source

BoxFilter(; dims, N, boundary)

Build a reusable, field-less box filter that computes a local box-average over the directions listed in dims. The returned object is callable: applying it to a field, bf(ψ), returns a KernelFunctionOperation. Build the filter once and reuse it across many fields, or pass it to other diagnostics that accept a filter.

dims is a tuple of distinct integers drawn from (1, 2, 3) (where 1, 2, 3 correspond to x, y, z). A single integer (e.g. dims=2) is accepted as shorthand for filtering along that one direction.

N is the total number of grid points used by the filter stencil along each filtered direction — i.e. how many cells are averaged together to produce one filtered output value (e.g. N=3 is a 3-point running mean, N=5 is a 5-point running mean). N must be an odd integer ≥ 3 so the stencil is symmetric around the current cell. (This is the size of the filter stencil — not the size of the grid.)

A multi-directional filter is assembled as a single KernelFunctionOperation whose kernel function is a 1D BoxFilterKernel{d₁}, with the next dimension's BoxFilterKernel{d₂} (and so on) threaded into the argument list. The box average is separable, so when the operation is the operand of a Field (the standard Field(bf(ψ)) / compute! path) it is evaluated as d sequential 1D passes through intermediate fields — d × N reads per output cell instead of Nᵈ. If the filtered field is composed into another AbstractOperation (e.g. 2 * bf(c)) the original fused single-kernel evaluation runs instead.

Boundary handling

Stencil offsets that leave the interior 1:Nd_grid of a direction (where Nd_grid is the number of cells along that direction) are handled per-direction. For Periodic directions offsets are always wrapped periodically, independent of the boundary keyword. For Bounded directions the boundary keyword picks the policy (default: :shrink):

:shrink — drop out-of-bounds offsets from both the sum and the count, so the filter is an honest local average whose effective stencil shrinks near a wall. This is the default for Bounded directions.
:edge — replicate the boundary-cell value (reads ψ[1] or ψ[Nd_grid] for offsets past either end).
(left=a, right=b) — pad with constant a on the low-index side and b on the high-index side (a and b are promoted to a common type).

boundary may be a single spec applied to every filtered dim, or a tuple with one entry per dim in dims (in the order the user passed them):

BoxFilter(; dims=(1, 2), N=7, boundary=:edge)
BoxFilter(; dims=(1, 2), N=7, boundary=(:shrink, :edge))
BoxFilter(; dims=(1,),   N=7, boundary=(left=0.0, right=1.0))

Because every policy wraps, clamps, or skips indices up front, halo_size(grid) does not constrain N: a small halo on a bounded direction is fine. The output location matches the location of the field the filter is applied to, which can be any input that supports the standard Oceananigans ψ[i, j, k] indexing contract (e.g. a Field or any AbstractOperation).

For Periodic directions the stencil must span at most one period: N ≤ 2*Nd_grid + 1, where Nd_grid is the number of cells along that direction. This is enforced when the filter is applied.

Examples

Build a box filter and apply it on a 2D (xz) grid that is periodic in x and bounded in z. The boundary=:edge keyword explicitly overrides the default :shrink policy for the bounded z-direction; x is Periodic so the boundary spec is silently ignored for that direction regardless.

julia> using Oceananigans, Oceanostics

julia> grid = RectilinearGrid(size=(8, 8), x=(0, 1), z=(0, 1),
                              topology=(Periodic, Flat, Bounded));

julia> c = CenterField(grid);

julia> bf = BoxFilter(; dims=(1, 3), N=5, boundary=:edge)
BoxFilter(dims=(1, 3), N=5, boundary=:edge)

julia> bf(c) isa KernelFunctionOperation
true

A one-step shortcut BoxFilter(ψ; dims, N, boundary) is also accepted, which applies the filter to ψ immediately (equivalent to BoxFilter(; dims, N, boundary)(ψ)).

source

Gaussian filter

The GaussianFilter computes a Gaussian-weighted local average of a field over one or more grid directions. Along a single filtered direction the filtered value at cell $i$ is

\[\bar{\psi}_i = \frac{\sum_m \Delta_m \, e^{-(x_m - x_i)^2 / (2\sigma^2)} \, \psi_m} {\sum_m \Delta_m \, e^{-(x_m - x_i)^2 / (2\sigma^2)}} ,\]

where $m$ runs over the cells of the stencil centred on $i$; $x_m$ and $\Delta_m$ are the coordinate and width of cell $m$ along the filtered direction, $x_i$ is the current cell's coordinate, and $\sigma$ is the kernel's standard deviation in physical units (the same as the grid spacing). In words: each cell's Gaussian weight $e^{-(x_m - x_i)^2 / (2\sigma^2)}$ is multiplied by that cell's own width $\Delta_m$, and the same $\Delta_m$-weighted sum appears in the normalizing denominator. The normalization returns a constant field unchanged and keeps every boundary policy consistent, while the $\Delta_m$ factor is the quadrature weight (the $dx'$ of $\int G_\sigma(x-x')\,\psi(x')\,dx' \big/ \int G_\sigma(x-x')\,dx'$) that makes the discrete sum approximate the continuous Gaussian convolution. Multi-direction filters apply this to each filtered direction in turn.

σ is the only required parameter beyond dims — N is inferred per-direction from σ and the minimum cell spacing in that direction so the stencil extends at least 2σ on each side of the current cell. To override, pass N explicitly: a single odd integer (≥ 3) applies to every filtered dim, or a tuple of odd integers sets one count per dim.

dims and boundary work identically to BoxFilter.

For a worked end-to-end example — coarse-graining a turbulent flow, a filter-width sweep, and a subfilter tracer flux — see the Spatial filtering example.

Variably spaced (stretched) grids

GaussianFilter works on stretched grids, but it is around 4 times slower than on regular grids:

Along a uniformly spaced direction every $\Delta_m$ is equal, so it factors out of the sum and cancels between numerator and denominator. The weights then reduce to a plain $e^{-(x_m - x_i)^2 / (2\sigma^2)}$ that is identical for every cell, so they are precomputed once and looked up — the fast path used on regular grids.
Along a variably spaced direction $x_m$ and $\Delta_m$ differ from cell to cell, so the weights cannot be precomputed and are evaluated on the fly. Keeping the per-cell $\Delta_m$ factor is what stops the average from being biased toward finely resolved regions; it makes the filter preserve constants exactly and linear fields to quadrature accuracy.

Directions are handled independently, so a grid that is uniform in x and y but stretched in z uses the fast path for the horizontal directions and the node-distance path for the vertical:

julia> stretched_grid = RectilinearGrid(size=(16, 16),
                                        x=(0, 1),
                                        z=k -> -1 + (k-1)/16 + 0.05sin(2π*(k-1)/16),
                                        topology=(Periodic, Flat, Bounded));

julia> cz = CenterField(stretched_grid); set!(cz, (x, z) -> sin(2π*x) * z);

julia> gf = GaussianFilter(; dims=(1, 3), σ=0.1)
GaussianFilter(dims=(1, 3), σ=0.1, N=nothing, boundary=:shrink)

julia> c̄z = Field(gf(cz));

julia> c̄z isa Field
true

Because the filter stores only its parameters (not a grid or field), the same gf applies to a field on any grid — here the uniform grid from the box-filter example above:

julia> c̄_gauss = Field(gf(c));

julia> c̄_gauss isa Field
true

Pass N to override the default stencil:

julia> gf_wide   = GaussianFilter(; dims=(1,), σ=0.05, N=11);

julia> gf_perdim = GaussianFilter(; dims=(1, 3), σ=0.05, N=(7, 11));

julia> (Field(gf_wide(c)) isa Field, Field(gf_perdim(c)) isa Field)
(true, true)

API reference

Oceanostics.SpatialFilters.GaussianFilter — Type

GaussianFilter(ψ; dims, σ, N, boundary)

One-step form: apply a Gaussian filter to ψ directly, returning the KernelFunctionOperation that computes its Gaussian-weighted local average. Equivalent to GaussianFilter(; dims, σ, N, boundary)(ψ).

Refer to GaussianFilter(; dims, σ, N, boundary) for the full description of the keyword arguments, the Gaussian weighting, stretched-grid handling, and boundary handling.

source

GaussianFilter(; dims, σ, N, boundary)

Build a reusable, field-less Gaussian filter that computes a Gaussian-weighted local average over the directions listed in dims. The returned object is callable: applying it to a field, gf(ψ), returns a KernelFunctionOperation. Build the filter once and reuse it across many fields, or pass it to other diagnostics that accept a filter.

σ is the standard deviation of the Gaussian kernel in physical units (the same units as the grid spacing). The filter approximates the continuous Gaussian convolution ∫ G_σ(x-x') ψ(x') dx' / ∫ G_σ(x-x') dx': a stencil cell of width Δₘ whose centre sits a physical distance r from the current cell centre contributes weight Δₘ · exp(-r² / 2σ²), and the result is normalized by the sum of the surviving weights, so all boundary policies behave consistently. On a uniform direction the Δₘ factor is constant and cancels, recovering the plain exp(-r² / 2σ²) weighting.

GaussianFilter supports both uniformly and variably spaced (stretched) directions, choosing the implementation per direction so the regular-grid case keeps its original speed:

Uniform direction — the weights are identical for every cell, so they are precomputed once in cell-offset units (σ_cells = σ / Δ, weight exp(-Δi² / 2σ_cells²) at cell offset Δi) and looked up in the hot loop. This is the fast path used on regular grids.
Stretched direction — the physical footprint of a fixed cell offset varies from cell to cell, so the weights cannot be precomputed; each is evaluated on the fly from the node coordinates and widths (Δₘ · exp(-(xₘ - xᵢ)² / 2σ²)). The cell-width factor Δₘ is the quadrature weight of the continuous convolution; it stops the average from being biased toward finely resolved regions and keeps constants (and, to quadrature accuracy, linear fields) preserved. The two paths agree exactly where they overlap (a uniform direction has xₘ - xᵢ = Δm·Δ and constant Δₘ). Directions are decided independently, so a grid that is uniform in x but stretched in z uses the fast path for x and the node-distance path for z.

N is the total number of grid points used by the filter stencil along each filtered direction — i.e. how many cells contribute to a single filtered output value. N must be an odd integer ≥ 3 so the stencil is symmetric around the current cell. (This is the size of the filter stencil — not the size of the grid.) If unspecified, N is inferred per-direction from σ and the smallest spacing along that direction so that the stencil extends at least 2σ to each side of the current cell everywhere (on a stretched direction it therefore reaches farther than 2σ where cells are large, which is harmless since the Gaussian weights there are tiny). To override, pass either a single odd integer (applied to every filtered dim) or a tuple with one odd-integer count per dim in dims (in the order the user passed them). For Periodic directions the stencil must span at most one period (N ≤ 2*Nd_grid + 1, where Nd_grid is the number of cells along that direction); this is enforced when the filter is applied.

See BoxFilter for the dims and boundary keyword documentation.

Performance notes

A multi-direction Gaussian filter is mathematically separable. Applying the filter returns a single composable KernelFunctionOperation, but when that operation is the operand of a Field (the standard Field(gf(ψ)) / compute! path), the implementation evaluates the filter as a sequence of 1D passes through intermediate fields. This reduces the per-output read count from N^d to d × N, which is the main reason multi-direction filters with wide stencils are competitive on GPUs. Mixed-spacing filters stage the same way — each direction's 1D pass uses its own (uniform or stretched) kernel. If the filtered field is composed into another AbstractOperation (e.g. 2 * gf(c)) it falls back to the fused, single-kernel evaluation.

Examples

julia> using Oceananigans, Oceanostics

julia> grid = RectilinearGrid(size=(8, 8), x=(0, 1), z=(0, 1),
                              topology=(Periodic, Flat, Bounded));

julia> c = CenterField(grid);

julia> gf = GaussianFilter(; dims=(1, 3), σ=0.1)
GaussianFilter(dims=(1, 3), σ=0.1, N=nothing, boundary=:shrink)

julia> gf(c) isa KernelFunctionOperation
true

A one-step shortcut GaussianFilter(ψ; dims, σ, N, boundary) is also accepted, which applies the filter to ψ immediately (equivalent to GaussianFilter(; dims, σ, N, boundary)(ψ)).

source