Neighborhood Filter

Neighborhood Filter#

A neighborhood filter replaces the value at each grid element with a reduction of all grid elements whose centers fall within a circular neighborhood of radius r degrees around that element.

Unlike a fixed k-nearest-neighbor average, a radius-based filter imposes a consistent spatial scale across the whole mesh—useful for variable-resolution grids where the number of neighbors varies from region to region.

Supported element types: face-centered, node-centered, and edge-centered data.

API at a glance:

Object	Method
`UxDataArray`	`da.neighborhood(r=5.0).mean()`
`UxDataset`	`ds.neighborhood(r=5.0).mean()`

The returned object is always the same type as the input, with the same grid, dims, coordinates, name, and attributes preserved.

Imports#

from functools import partial

import numpy as np

import uxarray as ux

Visualize the Unfiltered Field#

uxda.plot.polygons(
    cmap="RdBu_r",
    title="Original field (psi)",
    width=700,
    height=400,
)

Other Reductions#

Call the reduction you want as a method. The available ones are mean, sum, min, max, median, ptp, std, var, quantile, and percentile, each running a compiled kernel. Those taking a parameter declare it in their own signature: q for quantile (0–1) and percentile (0–100), ddof for std and var.

# One grouping, reduced four different ways
nb = uxda.neighborhood(r=5.0)

# 90th-percentile filter — highlights local maxima
uxda_p90 = nb.percentile(90)

# Maximum filter
uxda_max = nb.max()

# Median filter — robust to outliers
uxda_med = nb.median()

# Local spread, as a sample standard deviation
uxda_std = nb.std(ddof=1)

print("max filter max   :", uxda_max.values.max())
print("p90 filter max   :", uxda_p90.values.max())
print("median filter max:", uxda_med.values.max())
print("std filter max   :", uxda_std.values.max())

max filter max   : 1.5368961007894795
p90 filter max   : 1.5367148932640757
median filter max: 1.5361502577740587
std filter max   : 0.1639598135773616

(
    uxda_max.plot.polygons(
        cmap="RdBu_r", title="Max filter (r=5°)", width=350, height=250
    )
    + uxda_med.plot.polygons(
        cmap="RdBu_r", title="Median filter (r=5°)", width=350, height=250
    )
).cols(2)

Reductions Without a Method#

If you need something not in that list, hand it to reduce as a callable. It is applied as func(values, axis=-1) to a block whose last axis is the neighborhood, once per grid element, in Python — noticeably slower than the compiled methods, so reach for it only when none of them fits.

# a root-mean-square filter, which has no named equivalent
def rms(values, axis):
    return np.sqrt(np.mean(values**2, axis=axis))


uxda_rms = uxda.neighborhood(r=5.0).reduce(rms)

# `functools.partial` also works, though `.percentile()` is the faster way here
uxda_p90_slow = uxda.neighborhood(r=5.0).reduce(partial(np.percentile, q=90))
print(
    "partial matches the compiled reduction:",
    np.allclose(uxda_p90_slow.values, uxda_p90.values),
)

partial matches the compiled reduction: True

Reusing a Neighborhood Across Reductions#

Each call to neighborhood searches the grid for the neighbors of every element. That search usually costs far more than the reduction itself, so holding onto the object and reducing it several times — as the cell above does — already avoids repeating the expensive part.

Grid.neighborhood goes one step further and shares that search across variables. It is not bound to any data, so its reduction methods take the data as an argument.

nb5 = uxda.uxgrid.neighborhood(r=5.0)
nb5

<Neighborhood on='face centers' r=5.0 n_elements=5400 neighbors_per_element=[9, 12]>

# the search is already done; each of these only runs a reduction
smooth = nb5.mean(uxda)
spread = nb5.std(uxda)
p90 = nb5.percentile(uxda, 90)

print(
    "identical to the data-bound call:",
    np.allclose(smooth.values, uxda.neighborhood(r=5.0).mean().values),
)

identical to the data-bound call: True

n_neighbors reports how many elements fell inside each neighborhood. On a variable-resolution mesh this varies by region, which is worth checking before reading too much into a filtered field.

counts = nb5.n_neighbors
print(
    "neighbors per face: min",
    int(counts.min()),
    " max",
    int(counts.max()),
    " mean",
    float(counts.mean()).__round__(1),
)

neighbors per face: min 9  max 12  mean 9.7

Node- and Edge-Centered Data#

neighborhood works for any data element type. Here we create synthetic node- and edge-centered fields on a HEALPix grid and filter them.

uxgrid = ux.Grid.from_healpix(zoom=3)  # 768 faces, 770 nodes, 1536 edges

# Node-centered: a gradient along longitude
node_da = ux.UxDataArray(
    uxgrid.node_lon.values,
    dims=["n_node"],
    uxgrid=uxgrid,
    name="node_lon",
    attrs={"units": "degrees_east"},
)

filtered_node = node_da.neighborhood(r=10.0).mean()
print("node input  dims:", node_da.dims, "  shape:", node_da.shape)
print("node output dims:", filtered_node.dims, "  shape:", filtered_node.shape)
print("attrs preserved:", filtered_node.attrs)

node input  dims: ('n_node',)   shape: (770,)
node output dims: ('n_node',)   shape: (770,)
attrs preserved: {'units': 'degrees_east'}

# Edge-centered: a random field
rng = np.random.default_rng(42)
edge_da = ux.UxDataArray(
    rng.standard_normal(uxgrid.n_edge),
    dims=["n_edge"],
    uxgrid=uxgrid,
    name="edge_noise",
)

filtered_edge = edge_da.neighborhood(r=10.0).mean()
print("edge input  dims:", edge_da.dims, "  shape:", edge_da.shape)
print("edge output dims:", filtered_edge.dims, "  shape:", filtered_edge.shape)

edge input  dims: ('n_edge',)   shape: (1536,)
edge output dims: ('n_edge',)   shape: (1536,)

Multi-Dimensional Data (e.g. Time + Space)#

When a UxDataArray has extra leading dimensions (e.g. time), a neighborhood reduction applies independently at each time step and preserves the full dimension order.

uxds_ts = ux.tutorial.open_dataset("outCSne30-timeseries")
uxda_ts = uxds_ts["psi"]

print("Input  dims:", uxda_ts.dims, " shape:", uxda_ts.shape)

filtered_ts = uxda_ts.neighborhood(r=5.0).mean()

print("Output dims:", filtered_ts.dims, " shape:", filtered_ts.shape)

Input  dims: ('time', 'n_face')  shape: (6, 5400)

Output dims: ('time', 'n_face')  shape: (6, 5400)

The grid and time dimensions are both preserved. Because the filter is applied per time step, memory usage scales with n_time × n_face as expected.

Chaining with xarray Operations#

Because every reduction returns a proper UxDataArray with its uxgrid preserved, you can chain it with any standard xarray operation.

# Apply the filter and then mask values below zero
result = uxda.neighborhood(r=5.0).mean().where(lambda x: x > 0)
print("Masked result type:", type(result).__name__)
print("uxgrid preserved:", result.uxgrid is not None)
print("Positive fraction:", float((result > 0).sum()) / result.size)

Masked result type: UxDataArray
uxgrid preserved: True
Positive fraction: 1.0

# Group by latitude band after smoothing (standard xarray groupby)
import xarray as xr

lat_bins = xr.DataArray(
    np.digitize(uxda.uxgrid.face_lat.values, bins=np.arange(-90, 91, 30)),
    dims=["n_face"],
)

zonal_smooth = uxda.neighborhood(r=5.0).mean().groupby(lat_bins).mean()
print("Grouped result type:", type(zonal_smooth).__name__)
print("Zonal means:", zonal_smooth.values)

Grouped result type: UxDataArray
Zonal means: [1.10201444 1.03766258 1.00581415 0.99418587 0.96233741 0.89798555]

Radius Edge Cases#

Every element is its own neighbor at distance 0, and query_radius rejects a negative radius, so a neighborhood is never empty. r = 0 simply returns the original values, and a radius large enough to span the sphere returns the global reduction everywhere.

The output array is nonetheless allocated with NaN rather than uninitialized memory, so any unexpected gap would show up as an obvious NaN instead of garbage values.

uxgrid_coarse = ux.Grid.from_healpix(zoom=1)  # 48 faces
da_coarse = ux.UxDataArray(
    np.arange(uxgrid_coarse.n_face, dtype=float),
    dims=["n_face"],
    uxgrid=uxgrid_coarse,
)

# r = 0 catches the element itself → output matches the input exactly
filtered_r0 = da_coarse.neighborhood(r=0.0).mean()
print("r = 0:   unchanged?", np.allclose(filtered_r0.values, da_coarse.values))

# r = 360 catches every element → all values equal the global mean
filtered_global = da_coarse.neighborhood(r=360.0).mean()
print(
    "r = 360: all equal global mean?",
    np.allclose(filtered_global.values, da_coarse.values.mean()),
)

r = 0:   unchanged? True
r = 360: all equal global mean? True

Neighborhood Filter

Contents

Neighborhood Filter#

Imports#

Load Sample Data#

Visualize the Unfiltered Field#

Basic Usage: Mean Filter#

Effect of Radius#