From 75d832e36418fb020448d4f53f1fbd0cd4bc66f0 Mon Sep 17 00:00:00 2001
From: ahijevyc <ahijevyc@ucar.edu>
Date: Sat, 7 Sep 2024 11:09:00 -0600
Subject: [PATCH 01/15] Neighborhood filter

---
 uxarray/core/dataarray.py | 116 +++++++++++++++++++++++++++++++++++++-
 uxarray/core/dataset.py   |  36 ++++++++++++
 2 files changed, 150 insertions(+), 2 deletions(-)

diff --git a/uxarray/core/dataarray.py b/uxarray/core/dataarray.py
index e976c5816..aa7c82b27 100644
--- a/uxarray/core/dataarray.py
+++ b/uxarray/core/dataarray.py
@@ -6,6 +6,7 @@
 
 from typing import TYPE_CHECKING, Optional, Union, Hashable, Literal
 
+from uxarray.constants import GRID_DIMS
 from uxarray.formatting_html import array_repr
 
 from html import escape
@@ -1044,8 +1045,6 @@ def isel(self, ignore_grid=False, *args, **kwargs):
         > uxda.subset(n_node=[1, 2, 3])
         """
 
-        from uxarray.constants import GRID_DIMS
-
         if any(grid_dim in kwargs for grid_dim in GRID_DIMS) and not ignore_grid:
             # slicing a grid-dimension through Grid object
 
@@ -1102,3 +1101,116 @@ def _slice_from_grid(self, sliced_grid):
             dims=self.dims,
             attrs=self.attrs,
         )
+
+    def neighborhood_filter(
+        self,
+        func: Callable = np.mean,
+        r: float = 1.0,
+    ) -> UxDataArray:
+        """Apply neighborhood filter
+        Parameters:
+        -----------
+        func: Callable, default=np.mean
+            Apply this function to neighborhood
+        r : float, default=1.
+            Radius of neighborhood. For spherical coordinates, the radius is in units of degrees,
+            and for cartesian coordinates, the radius is in meters.
+        Returns:
+        --------
+        destination_data : np.ndarray
+            Filtered data.
+        """
+
+        if self._face_centered():
+            data_mapping = "face centers"
+        elif self._node_centered():
+            data_mapping = "nodes"
+        elif self._edge_centered():
+            data_mapping = "edge centers"
+        else:
+            raise ValueError(
+                f"Data_mapping is not face, node, or edge. Could not define data_mapping."
+            )
+
+        # reconstruct because the cached tree could be built from
+        # face centers, edge centers or nodes.
+        tree = self.uxgrid.get_ball_tree(coordinates=data_mapping, reconstruct=True)
+
+        coordinate_system = tree.coordinate_system
+
+        if coordinate_system == "spherical":
+            if data_mapping == "nodes":
+                lon, lat = (
+                    self.uxgrid.node_lon.values,
+                    self.uxgrid.node_lat.values,
+                )
+            elif data_mapping == "face centers":
+                lon, lat = (
+                    self.uxgrid.face_lon.values,
+                    self.uxgrid.face_lat.values,
+                )
+            elif data_mapping == "edge centers":
+                lon, lat = (
+                    self.uxgrid.edge_lon.values,
+                    self.uxgrid.edge_lat.values,
+                )
+            else:
+                raise ValueError(
+                    f"Invalid data_mapping. Expected 'nodes', 'edge centers', or 'face centers', "
+                    f"but received: {data_mapping}"
+                )
+
+            dest_coords = np.c_[lon, lat]
+
+        elif coordinate_system == "cartesian":
+            if data_mapping == "nodes":
+                x, y, z = (
+                    self.uxgrid.node_x.values,
+                    self.uxgrid.node_y.values,
+                    self.uxgrid.node_z.values,
+                )
+            elif data_mapping == "face centers":
+                x, y, z = (
+                    self.uxgrid.face_x.values,
+                    self.uxgrid.face_y.values,
+                    self.uxgrid.face_z.values,
+                )
+            elif data_mapping == "edge centers":
+                x, y, z = (
+                    self.uxgrid.edge_x.values,
+                    self.uxgrid.edge_y.values,
+                    self.uxgrid.edge_z.values,
+                )
+            else:
+                raise ValueError(
+                    f"Invalid data_mapping. Expected 'nodes', 'edge centers', or 'face centers', "
+                    f"but received: {data_mapping}"
+                )
+
+            dest_coords = np.c_[x, y, z]
+
+        else:
+            raise ValueError(
+                f"Invalid coordinate_system. Expected either 'spherical' or 'cartesian', but received {coordinate_system}"
+            )
+
+        neighbor_indices = tree.query_radius(dest_coords, r=r)
+
+        destination_data = np.empty(self.data.shape)
+
+        # assert last dimension is a GRID dimension.
+        assert self.dims[-1] in GRID_DIMS, (
+                f"expected last dimension of uxDataArray {self.data.dims[-1]} "
+                f"to be one of {GRID_DIMS}"
+                )
+        # Apply function to indices on last axis.
+        for i, idx in enumerate(neighbor_indices):
+            if len(idx):
+                destination_data[..., i] = func(self.data[..., idx])
+
+        # construct data array for filtered variable
+        uxda_filter = self._copy()
+
+        uxda_filter.data = destination_data
+
+        return uxda_filter
diff --git a/uxarray/core/dataset.py b/uxarray/core/dataset.py
index 9a2f522a0..4f2786704 100644
--- a/uxarray/core/dataset.py
+++ b/uxarray/core/dataset.py
@@ -7,6 +7,7 @@
 
 from typing import Optional, IO, Union
 
+from uxarray.constants import GRID_DIMS
 from uxarray.grid import Grid
 from uxarray.core.dataarray import UxDataArray
 
@@ -338,6 +339,41 @@ def to_array(self) -> UxDataArray:
         xarr = super().to_array()
         return UxDataArray(xarr, uxgrid=self.uxgrid)
 
+    def neighborhood_filter(
+        self,
+        func: Callable = np.mean,
+        r: float = 1.0,
+    ):
+        """Neighborhood function implementation for ``UxDataset``.
+        Parameters
+        ---------
+        func : Callable = np.mean
+            Apply this function to neighborhood
+        r : float, default=1.
+            Radius of neighborhood
+        """
+
+
+        destination_uxds = self._copy()
+        # Loop through uxDataArrays in uxDataset
+        for var_name in self.data_vars:
+            uxda = self[var_name]
+
+            # Skip if uxDataArray has no GRID dimension.
+            grid_dims = [dim for dim in uxda.dims if dim in GRID_DIMS]
+            if len(grid_dims) == 0:
+                continue
+
+            # Put GRID dimension last for UxDataArray.neighborhood_filter.
+            remember_dim_order = uxda.dims
+            uxda = uxda.transpose(..., grid_dims[0])
+            # Filter uxDataArray.
+            uxda = uxda.neighborhood_filter(func, r)
+            # Restore old dimension order.
+            destination_uxds[var_name] = uxda.transpose(*remember_dim_order)
+
+        return destination_uxds
+
     def nearest_neighbor_remap(
         self,
         destination_obj: Union[Grid, UxDataArray, UxDataset],

From 8ec019328339dc92d719d7c28e87628de78de197 Mon Sep 17 00:00:00 2001
From: ahijevyc <ahijevyc@ucar.edu>
Date: Mon, 9 Sep 2024 10:48:09 -0600
Subject: [PATCH 02/15] ruff recommendations

---
 uxarray/core/dataarray.py | 8 ++++----
 uxarray/core/dataset.py   | 1 -
 2 files changed, 4 insertions(+), 5 deletions(-)

diff --git a/uxarray/core/dataarray.py b/uxarray/core/dataarray.py
index 09764c400..8cf3fddf2 100644
--- a/uxarray/core/dataarray.py
+++ b/uxarray/core/dataarray.py
@@ -1131,7 +1131,7 @@ def neighborhood_filter(
             data_mapping = "edge centers"
         else:
             raise ValueError(
-                f"Data_mapping is not face, node, or edge. Could not define data_mapping."
+                "Data_mapping is not face, node, or edge. Could not define data_mapping."
             )
 
         # reconstruct because the cached tree could be built from
@@ -1202,9 +1202,9 @@ def neighborhood_filter(
 
         # assert last dimension is a GRID dimension.
         assert self.dims[-1] in GRID_DIMS, (
-                f"expected last dimension of uxDataArray {self.data.dims[-1]} "
-                f"to be one of {GRID_DIMS}"
-                )
+            f"expected last dimension of uxDataArray {self.data.dims[-1]} "
+            f"to be one of {GRID_DIMS}"
+        )
         # Apply function to indices on last axis.
         for i, idx in enumerate(neighbor_indices):
             if len(idx):
diff --git a/uxarray/core/dataset.py b/uxarray/core/dataset.py
index 4f2786704..dbf840903 100644
--- a/uxarray/core/dataset.py
+++ b/uxarray/core/dataset.py
@@ -353,7 +353,6 @@ def neighborhood_filter(
             Radius of neighborhood
         """
 
-
         destination_uxds = self._copy()
         # Loop through uxDataArrays in uxDataset
         for var_name in self.data_vars:

From 5605949bfaa4bfc2a86c801f7103e34b6fdb765b Mon Sep 17 00:00:00 2001
From: ahijevyc <ahijevyc@ucar.edu>
Date: Mon, 9 Sep 2024 10:57:41 -0600
Subject: [PATCH 03/15] added Callable to Type checking

---
 uxarray/core/dataarray.py | 2 +-
 uxarray/core/dataset.py   | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/uxarray/core/dataarray.py b/uxarray/core/dataarray.py
index 8cf3fddf2..a9076d7ab 100644
--- a/uxarray/core/dataarray.py
+++ b/uxarray/core/dataarray.py
@@ -4,7 +4,7 @@
 import numpy as np
 
 
-from typing import TYPE_CHECKING, Optional, Union, Hashable, Literal
+from typing import TYPE_CHECKING, Callable, Optional, Union, Hashable, Literal
 
 from uxarray.constants import GRID_DIMS
 from uxarray.formatting_html import array_repr
diff --git a/uxarray/core/dataset.py b/uxarray/core/dataset.py
index dbf840903..f4c259297 100644
--- a/uxarray/core/dataset.py
+++ b/uxarray/core/dataset.py
@@ -5,7 +5,7 @@
 
 import sys
 
-from typing import Optional, IO, Union
+from typing import Callable, Optional, IO, Union
 
 from uxarray.constants import GRID_DIMS
 from uxarray.grid import Grid

From 0c7bc1eae7474f71c8b00a5a2060c5d3c08b2d6e Mon Sep 17 00:00:00 2001
From: ahijevyc <ahijevyc@ucar.edu>
Date: Mon, 9 Sep 2024 15:10:27 -0600
Subject: [PATCH 04/15] np.vstack().T faster than np.c

---
 uxarray/core/dataarray.py | 4 ++--
 uxarray/core/dataset.py   | 3 ++-
 2 files changed, 4 insertions(+), 3 deletions(-)

diff --git a/uxarray/core/dataarray.py b/uxarray/core/dataarray.py
index a9076d7ab..8fae278a9 100644
--- a/uxarray/core/dataarray.py
+++ b/uxarray/core/dataarray.py
@@ -1162,7 +1162,7 @@ def neighborhood_filter(
                     f"but received: {data_mapping}"
                 )
 
-            dest_coords = np.c_[lon, lat]
+            dest_coords = np.vstack((lon, lat)).T
 
         elif coordinate_system == "cartesian":
             if data_mapping == "nodes":
@@ -1189,7 +1189,7 @@ def neighborhood_filter(
                     f"but received: {data_mapping}"
                 )
 
-            dest_coords = np.c_[x, y, z]
+            dest_coords = np.vstack((x, y, z)).T
 
         else:
             raise ValueError(
diff --git a/uxarray/core/dataset.py b/uxarray/core/dataset.py
index f4c259297..2489c23ab 100644
--- a/uxarray/core/dataset.py
+++ b/uxarray/core/dataset.py
@@ -350,7 +350,8 @@ def neighborhood_filter(
         func : Callable = np.mean
             Apply this function to neighborhood
         r : float, default=1.
-            Radius of neighborhood
+            Radius of neighborhood. For spherical coordinates, the radius is in units of degrees,
+            and for cartesian coordinates, the radius is in meters.
         """
 
         destination_uxds = self._copy()

From d6d8a33faa8f64dec3fb930ef9ba53f25b63ff1d Mon Sep 17 00:00:00 2001
From: ahijevyc <ahijevyc@ucar.edu>
Date: Mon, 9 Sep 2024 15:25:26 -0600
Subject: [PATCH 05/15] Fix some comments

---
 uxarray/core/dataarray.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/uxarray/core/dataarray.py b/uxarray/core/dataarray.py
index 8fae278a9..a0b61931c 100644
--- a/uxarray/core/dataarray.py
+++ b/uxarray/core/dataarray.py
@@ -1198,9 +1198,10 @@ def neighborhood_filter(
 
         neighbor_indices = tree.query_radius(dest_coords, r=r)
 
+        # Construct numpy array for filtered variable.
         destination_data = np.empty(self.data.shape)
 
-        # assert last dimension is a GRID dimension.
+        # Assert last dimension is a GRID dimension.
         assert self.dims[-1] in GRID_DIMS, (
             f"expected last dimension of uxDataArray {self.data.dims[-1]} "
             f"to be one of {GRID_DIMS}"
@@ -1210,7 +1211,7 @@ def neighborhood_filter(
             if len(idx):
                 destination_data[..., i] = func(self.data[..., idx])
 
-        # construct data array for filtered variable
+        # Construct UxDataArray for filtered variable.
         uxda_filter = self._copy()
 
         uxda_filter.data = destination_data

From 6c59af767e4674d2b852e7b30826bba5485614f5 Mon Sep 17 00:00:00 2001
From: ahijevyc <ahijevyc@ucar.edu>
Date: Mon, 17 Mar 2025 15:27:53 -0600
Subject: [PATCH 06/15] missing imports

---
 uxarray/core/dataset.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/uxarray/core/dataset.py b/uxarray/core/dataset.py
index e1a2a923c..401bdf5da 100644
--- a/uxarray/core/dataset.py
+++ b/uxarray/core/dataset.py
@@ -8,6 +8,8 @@
 
 import numpy as np
 import xarray as xr
+from xarray.core import dtypes
+from xarray.core.options import OPTIONS
 from xarray.core.utils import UncachedAccessor
 
 import uxarray

From 14c37b72caae15d43bf32bf2d0074b0cfd4a3e05 Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Thu, 26 Feb 2026 19:06:58 +0000
Subject: [PATCH 07/15] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 uxarray/core/dataarray.py | 2 +-
 uxarray/core/dataset.py   | 2 --
 2 files changed, 1 insertion(+), 3 deletions(-)

diff --git a/uxarray/core/dataarray.py b/uxarray/core/dataarray.py
index b14359b24..5aa287b0c 100644
--- a/uxarray/core/dataarray.py
+++ b/uxarray/core/dataarray.py
@@ -14,13 +14,13 @@
 from xarray.core.utils import UncachedAccessor
 
 import uxarray
+from uxarray.constants import GRID_DIMS
 from uxarray.core.aggregation import _uxda_grid_aggregate
 from uxarray.core.gradient import (
     _calculate_edge_face_difference,
     _calculate_edge_node_difference,
     _compute_gradient,
 )
-from uxarray.constants import GRID_DIMS
 from uxarray.core.utils import _map_dims_to_ugrid
 from uxarray.core.zonal import (
     _compute_conservative_zonal_mean_bands,
diff --git a/uxarray/core/dataset.py b/uxarray/core/dataset.py
index f1b7326a0..684f0be3c 100644
--- a/uxarray/core/dataset.py
+++ b/uxarray/core/dataset.py
@@ -611,7 +611,6 @@ def to_array(self) -> UxDataArray:
         xarr = super().to_array()
         return UxDataArray(xarr, uxgrid=self.uxgrid)
 
-
     def neighborhood_filter(
         self,
         func: Callable = np.mean,
@@ -668,7 +667,6 @@ def to_xarray(self, grid_format: str = "UGRID") -> xr.Dataset:
 
         return xr.Dataset(self)
 
-
     def get_dual(self):
         """Compute the dual mesh for a dataset, returns a new dataset object.
 

From a37b88fb050bef815a18acb1b01fd05826c16f09 Mon Sep 17 00:00:00 2001
From: Orhan Eroglu <32553057+erogluorhan@users.noreply.github.com>
Date: Thu, 26 Feb 2026 12:13:54 -0700
Subject: [PATCH 08/15] Update dataset.py to address pre-commit errors

---
 uxarray/core/dataset.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/uxarray/core/dataset.py b/uxarray/core/dataset.py
index 684f0be3c..9e723943e 100644
--- a/uxarray/core/dataset.py
+++ b/uxarray/core/dataset.py
@@ -3,7 +3,7 @@
 import os
 import sys
 from html import escape
-from typing import IO, Any, Callable, Union
+from typing import IO, Any, Callable, Mapping
 from warnings import warn
 
 import numpy as np

From 02ca8b7690b1edb23491efc65c19d13164dc39ca Mon Sep 17 00:00:00 2001
From: Rajeev Jain <rajeeja@gmail.com>
Date: Fri, 24 Jul 2026 12:12:59 -0500
Subject: [PATCH 09/15] Address review feedback for neighborhood_filter

- Fix Grid.get_ball_tree()/get_kd_tree() caching to rebuild when the
  coordinates or coordinate_system changes, instead of mutating the
  cached tree in place (per review discussion).
- Move the bulk of UxDataArray.neighborhood_filter's computation into
  a new _neighborhood_filter() helper in uxarray.grid.neighbors, with
  UxDataArray/UxDataset now acting as thin wrappers around it.
- Fix a bug where func was applied across all axes instead of just the
  grid axis, which collapsed extra dimensions (e.g. time) in the
  output.
- Bring UxDataset.neighborhood_filter's docstring in line with the
  UxDataArray implementation.
- Add tests for face/node/edge-centered data, custom functions via
  functools.partial, extra dimension preservation, and dataset-level
  filtering.
- Document neighborhood_filter in docs/api.rst.
---
 docs/api.rst                |  13 +++++
 test/core/test_dataarray.py |  93 ++++++++++++++++++++++++++++++
 test/core/test_dataset.py   |  36 ++++++++++++
 uxarray/core/dataarray.py   |  97 +++++++-------------------------
 uxarray/core/dataset.py     |  21 +++++--
 uxarray/grid/grid.py        |  21 ++++---
 uxarray/grid/neighbors.py   | 109 ++++++++++++++++++++++++++++++++++++
 7 files changed, 298 insertions(+), 92 deletions(-)

diff --git a/docs/api.rst b/docs/api.rst
index dd89363ba..75150e691 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -565,6 +565,19 @@ Azimuthal aggregations apply an aggregation (i.e. averaging) along circles of co
    UxDataArray.azimuthal_mean
 
 
+Neighborhood
+~~~~~~~~~~~~
+
+Neighborhood filters apply an aggregation (i.e. averaging) to all grid elements within a circular
+neighborhood of a specified radius around each grid element.
+
+.. autosummary::
+   :toctree: generated/
+
+   UxDataArray.neighborhood_filter
+   UxDataset.neighborhood_filter
+
+
 Zonal Average
 ~~~~~~~~~~~~~
 .. autosummary::
diff --git a/test/core/test_dataarray.py b/test/core/test_dataarray.py
index b961a5a37..f172539a4 100644
--- a/test/core/test_dataarray.py
+++ b/test/core/test_dataarray.py
@@ -163,3 +163,96 @@ def test_data_location():
         np.ones((3, uxgrid.n_face)), dims=["time", "n_face"], uxgrid=uxgrid
     )
     assert face_time.data_location == "face_centered"
+
+
+class TestNeighborhoodFilter:
+    """Tests for ``UxDataArray.neighborhood_filter``."""
+
+    def test_face_centered(self, gridpath, datasetpath):
+        """A large enough radius should average every face together."""
+        uxds = ux.open_dataset(
+            gridpath("ugrid", "outCSne30", "outCSne30.ug"),
+            datasetpath("ugrid", "outCSne30", "outCSne30_vortex.nc"),
+        )
+        uxda = uxds["psi"]
+
+        # radius of 0 should select each face's own coordinate, leaving the
+        # data unchanged
+        filtered = uxda.neighborhood_filter(func=np.mean, r=0.0)
+        np.testing.assert_allclose(filtered.values, uxda.values)
+
+        # a large enough radius should include the entire grid in the
+        # neighborhood of every face, so every filtered value should match
+        # the global mean of the field
+        filtered_all = uxda.neighborhood_filter(func=np.mean, r=360.0)
+        np.testing.assert_allclose(filtered_all.values, uxda.values.mean())
+
+        assert isinstance(filtered, UxDataArray)
+        assert filtered.uxgrid == uxda.uxgrid
+        assert filtered.dims == uxda.dims
+        assert filtered.shape == uxda.shape
+
+    def test_node_centered(self):
+        """Neighborhood filter should work for node-centered data."""
+        uxgrid = ux.Grid.from_healpix(zoom=1)
+        data = np.arange(uxgrid.n_node, dtype=float)
+        uxda = UxDataArray(data, dims=["n_node"], uxgrid=uxgrid, name="node_var")
+
+        filtered = uxda.neighborhood_filter(func=np.mean, r=0.0)
+        np.testing.assert_allclose(filtered.values, data)
+
+        filtered_all = uxda.neighborhood_filter(func=np.mean, r=360.0)
+        np.testing.assert_allclose(filtered_all.values, data.mean())
+
+    def test_edge_centered(self):
+        """Neighborhood filter should work for edge-centered data."""
+        uxgrid = ux.Grid.from_healpix(zoom=1)
+        data = np.arange(uxgrid.n_edge, dtype=float)
+        uxda = UxDataArray(data, dims=["n_edge"], uxgrid=uxgrid, name="edge_var")
+
+        filtered = uxda.neighborhood_filter(func=np.mean, r=0.0)
+        np.testing.assert_allclose(filtered.values, data)
+
+    def test_custom_func_with_partial(self, gridpath, datasetpath):
+        """A user-defined function (i.e. ``functools.partial``) should work."""
+        from functools import partial
+
+        uxds = ux.open_dataset(
+            gridpath("ugrid", "outCSne30", "outCSne30.ug"),
+            datasetpath("ugrid", "outCSne30", "outCSne30_vortex.nc"),
+        )
+        uxda = uxds["psi"]
+
+        filtered_max = uxda.neighborhood_filter(func=np.max, r=5.0)
+        filtered_percentile = uxda.neighborhood_filter(
+            func=partial(np.percentile, q=100), r=5.0
+        )
+
+        np.testing.assert_allclose(filtered_max.values, filtered_percentile.values)
+
+    def test_extra_dimension_preserved(self, gridpath, datasetpath):
+        """An extra leading (i.e. time) dimension should be preserved."""
+        uxds = ux.open_dataset(
+            gridpath("ugrid", "outCSne30", "outCSne30.ug"),
+            datasetpath("ugrid", "outCSne30", "outCSne30_vortex.nc"),
+        )
+        uxda = uxds["psi"]
+
+        data = np.stack([uxda.values, uxda.values * 2.0])
+        uxda_time = UxDataArray(
+            data, dims=["time", "n_face"], uxgrid=uxda.uxgrid, name="psi_time"
+        )
+
+        filtered = uxda_time.neighborhood_filter(func=np.mean, r=0.0)
+
+        assert filtered.dims == uxda_time.dims
+        assert filtered.shape == uxda_time.shape
+        np.testing.assert_allclose(filtered.values, data)
+
+    def test_invalid_data_location(self):
+        """Data that is not mapped to a grid element should raise an error."""
+        uxgrid = ux.Grid.from_healpix(zoom=1)
+        uxda = UxDataArray(np.ones(5), dims=["other_dim"], uxgrid=uxgrid)
+
+        with pytest.raises(ValueError):
+            uxda.neighborhood_filter(func=np.mean, r=1.0)
diff --git a/test/core/test_dataset.py b/test/core/test_dataset.py
index ace56be19..cc32e1db4 100644
--- a/test/core/test_dataset.py
+++ b/test/core/test_dataset.py
@@ -170,3 +170,39 @@ def test_uxdataset_to_array():
     assert arr1.name is None
     arr2 = uxds.to_array(dim='custom_dim', name='custom_name')
     assert arr2.name == 'custom_name'
+
+
+class TestNeighborhoodFilter:
+    """Tests for ``UxDataset.neighborhood_filter``."""
+
+    def test_face_centered(self, gridpath, datasetpath):
+        """Ensures the dataset-level filter matches the per-variable
+        ``UxDataArray.neighborhood_filter`` results."""
+        uxds = ux.open_dataset(
+            gridpath("ugrid", "outCSne30", "outCSne30.ug"),
+            datasetpath("ugrid", "outCSne30", "outCSne30_vortex.nc"),
+        )
+
+        filtered_ds = uxds.neighborhood_filter(func=np.mean, r=5.0)
+        filtered_da = uxds["psi"].neighborhood_filter(func=np.mean, r=5.0)
+
+        assert isinstance(filtered_ds, UxDataset)
+        nt.assert_allclose(filtered_ds["psi"].values, filtered_da.values)
+
+    def test_non_grid_variable_skipped(self):
+        """Data variables without a grid dimension should be left
+        untouched."""
+        uxgrid = ux.Grid.from_healpix(zoom=1)
+
+        uxds = UxDataset(
+            data_vars={
+                "face_var": ("n_face", np.arange(uxgrid.n_face, dtype=float)),
+                "scalar_var": ("other_dim", np.array([1.0, 2.0, 3.0])),
+            },
+            uxgrid=uxgrid,
+        )
+
+        filtered = uxds.neighborhood_filter(func=np.mean, r=0.0)
+
+        nt.assert_allclose(filtered["face_var"].values, uxds["face_var"].values)
+        nt.assert_allclose(filtered["scalar_var"].values, uxds["scalar_var"].values)
diff --git a/uxarray/core/dataarray.py b/uxarray/core/dataarray.py
index 6a4ca6319..630b4b959 100644
--- a/uxarray/core/dataarray.py
+++ b/uxarray/core/dataarray.py
@@ -31,6 +31,7 @@
 from uxarray.formatting_html import array_repr
 from uxarray.grid import Grid
 from uxarray.grid.dual import construct_dual
+from uxarray.grid.neighbors import _neighborhood_filter
 from uxarray.grid.validation import _check_duplicate_nodes_indices
 from uxarray.io._healpix import get_zoom_from_cells
 from uxarray.plot.accessor import UxDataArrayPlotAccessor
@@ -2192,17 +2193,24 @@ def neighborhood_filter(
         func: Callable = np.mean,
         r: float = 1.0,
     ) -> UxDataArray:
-        """Apply neighborhood filter
-        Parameters:
-        -----------
+        """Apply a neighborhood filter, replacing the value at each grid
+        element with ``func`` applied to all elements within a circular
+        neighborhood of radius ``r``.
+
+        Parameters
+        ----------
         func: Callable, default=np.mean
-            Apply this function to neighborhood
+            Apply this function to neighborhood. Must accept an ``axis`` keyword
+            argument (as ``np.mean``, ``np.median``, and similar NumPy reductions
+            do). Use ``functools.partial`` to supply additional arguments, e.g.
+            ``functools.partial(np.percentile, q=90)``.
         r : float, default=1.
             Radius of neighborhood. For spherical coordinates, the radius is in units of degrees,
             and for cartesian coordinates, the radius is in meters.
-        Returns:
-        --------
-        destination_data : np.ndarray
+
+        Returns
+        -------
+        uxda_filter : UxDataArray
             Filtered data.
         """
 
@@ -2217,82 +2225,15 @@ def neighborhood_filter(
                 "Data_mapping is not face, node, or edge. Could not define data_mapping."
             )
 
-        # reconstruct because the cached tree could be built from
-        # face centers, edge centers or nodes.
-        tree = self.uxgrid.get_ball_tree(coordinates=data_mapping, reconstruct=True)
-
-        coordinate_system = tree.coordinate_system
-
-        if coordinate_system == "spherical":
-            if data_mapping == "nodes":
-                lon, lat = (
-                    self.uxgrid.node_lon.values,
-                    self.uxgrid.node_lat.values,
-                )
-            elif data_mapping == "face centers":
-                lon, lat = (
-                    self.uxgrid.face_lon.values,
-                    self.uxgrid.face_lat.values,
-                )
-            elif data_mapping == "edge centers":
-                lon, lat = (
-                    self.uxgrid.edge_lon.values,
-                    self.uxgrid.edge_lat.values,
-                )
-            else:
-                raise ValueError(
-                    f"Invalid data_mapping. Expected 'nodes', 'edge centers', or 'face centers', "
-                    f"but received: {data_mapping}"
-                )
-
-            dest_coords = np.vstack((lon, lat)).T
-
-        elif coordinate_system == "cartesian":
-            if data_mapping == "nodes":
-                x, y, z = (
-                    self.uxgrid.node_x.values,
-                    self.uxgrid.node_y.values,
-                    self.uxgrid.node_z.values,
-                )
-            elif data_mapping == "face centers":
-                x, y, z = (
-                    self.uxgrid.face_x.values,
-                    self.uxgrid.face_y.values,
-                    self.uxgrid.face_z.values,
-                )
-            elif data_mapping == "edge centers":
-                x, y, z = (
-                    self.uxgrid.edge_x.values,
-                    self.uxgrid.edge_y.values,
-                    self.uxgrid.edge_z.values,
-                )
-            else:
-                raise ValueError(
-                    f"Invalid data_mapping. Expected 'nodes', 'edge centers', or 'face centers', "
-                    f"but received: {data_mapping}"
-                )
-
-            dest_coords = np.vstack((x, y, z)).T
-
-        else:
-            raise ValueError(
-                f"Invalid coordinate_system. Expected either 'spherical' or 'cartesian', but received {coordinate_system}"
-            )
-
-        neighbor_indices = tree.query_radius(dest_coords, r=r)
-
-        # Construct numpy array for filtered variable.
-        destination_data = np.empty(self.data.shape)
-
         # Assert last dimension is a GRID dimension.
         assert self.dims[-1] in GRID_DIMS, (
             f"expected last dimension of uxDataArray {self.data.dims[-1]} "
             f"to be one of {GRID_DIMS}"
         )
-        # Apply function to indices on last axis.
-        for i, idx in enumerate(neighbor_indices):
-            if len(idx):
-                destination_data[..., i] = func(self.data[..., idx])
+
+        destination_data = _neighborhood_filter(
+            self.uxgrid, self.data, data_mapping, func=func, r=r
+        )
 
         # Construct UxDataArray for filtered variable.
         uxda_filter = self._copy()
diff --git a/uxarray/core/dataset.py b/uxarray/core/dataset.py
index 417f84a44..d02180411 100644
--- a/uxarray/core/dataset.py
+++ b/uxarray/core/dataset.py
@@ -660,15 +660,26 @@ def neighborhood_filter(
         self,
         func: Callable = np.mean,
         r: float = 1.0,
-    ):
-        """Neighborhood function implementation for ``UxDataset``.
+    ) -> UxDataset:
+        """Apply a neighborhood filter, replacing the value at each grid
+        element of every data variable with ``func`` applied to all elements
+        within a circular neighborhood of radius ``r``.
+
         Parameters
-        ---------
-        func : Callable = np.mean
-            Apply this function to neighborhood
+        ----------
+        func: Callable, default=np.mean
+            Apply this function to neighborhood. Must accept an ``axis`` keyword
+            argument (as ``np.mean``, ``np.median``, and similar NumPy reductions
+            do). Use ``functools.partial`` to supply additional arguments, e.g.
+            ``functools.partial(np.percentile, q=90)``.
         r : float, default=1.
             Radius of neighborhood. For spherical coordinates, the radius is in units of degrees,
             and for cartesian coordinates, the radius is in meters.
+
+        Returns
+        -------
+        destination_uxds : UxDataset
+            Filtered dataset.
         """
 
         destination_uxds = self._copy()
diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py
index 2c753f85a..25b5f6b0a 100644
--- a/uxarray/grid/grid.py
+++ b/uxarray/grid/grid.py
@@ -1772,7 +1772,12 @@ def get_ball_tree(
             BallTree instance
         """
 
-        if self._ball_tree is None or reconstruct:
+        if (
+            self._ball_tree is None
+            or coordinates != self._ball_tree._coordinates
+            or coordinate_system != self._ball_tree.coordinate_system
+            or reconstruct
+        ):
             self._ball_tree = BallTree(
                 self,
                 coordinates=coordinates,
@@ -1780,9 +1785,6 @@ def get_ball_tree(
                 coordinate_system=coordinate_system,
                 reconstruct=reconstruct,
             )
-        else:
-            if coordinates != self._ball_tree._coordinates:
-                self._ball_tree.coordinates = coordinates
 
         return self._ball_tree
 
@@ -1872,7 +1874,12 @@ def get_kd_tree(
             KDTree instance
         """
 
-        if self._kd_tree is None or reconstruct:
+        if (
+            self._kd_tree is None
+            or coordinates != self._kd_tree._coordinates
+            or coordinate_system != self._kd_tree.coordinate_system
+            or reconstruct
+        ):
             self._kd_tree = KDTree(
                 self,
                 coordinates=coordinates,
@@ -1881,10 +1888,6 @@ def get_kd_tree(
                 reconstruct=reconstruct,
             )
 
-        else:
-            if coordinates != self._kd_tree._coordinates:
-                self._kd_tree.coordinates = coordinates
-
         return self._kd_tree
 
     def get_spatial_hash(
diff --git a/uxarray/grid/neighbors.py b/uxarray/grid/neighbors.py
index 1c4d4f145..a34d15cda 100644
--- a/uxarray/grid/neighbors.py
+++ b/uxarray/grid/neighbors.py
@@ -1,3 +1,5 @@
+from typing import Callable
+
 import numpy as np
 import xarray as xr
 from numba import njit
@@ -1129,3 +1131,110 @@ def _construct_edge_face_distances(face_lon, face_lat, edge_faces):
     )
 
     return edge_face_distances
+
+
+def _get_element_coords(grid, data_mapping: str, coordinate_system: str):
+    """Gathers the coordinate array used to query a ``BallTree`` for a given
+    grid element location and coordinate system.
+
+    Parameters
+    ----------
+    grid : Grid
+        Source grid containing the coordinate arrays.
+    data_mapping : str
+        One of "nodes", "edge centers", or "face centers".
+    coordinate_system : str
+        Either "spherical" or "cartesian".
+
+    Returns
+    -------
+    coords : np.ndarray
+        Array of shape (n_elements, 2) for "spherical" (lon, lat) or
+        (n_elements, 3) for "cartesian" (x, y, z).
+    """
+    prefix_map = {
+        "nodes": "node",
+        "edge centers": "edge",
+        "face centers": "face",
+    }
+
+    if data_mapping not in prefix_map:
+        raise ValueError(
+            f"Invalid data_mapping. Expected 'nodes', 'edge centers', or 'face centers', "
+            f"but received: {data_mapping}"
+        )
+
+    prefix = prefix_map[data_mapping]
+
+    if coordinate_system == "spherical":
+        lon = getattr(grid, f"{prefix}_lon").values
+        lat = getattr(grid, f"{prefix}_lat").values
+        return np.vstack((lon, lat)).T
+
+    elif coordinate_system == "cartesian":
+        x = getattr(grid, f"{prefix}_x").values
+        y = getattr(grid, f"{prefix}_y").values
+        z = getattr(grid, f"{prefix}_z").values
+        return np.vstack((x, y, z)).T
+
+    else:
+        raise ValueError(
+            f"Invalid coordinate_system. Expected either 'spherical' or 'cartesian', "
+            f"but received {coordinate_system}"
+        )
+
+
+def _neighborhood_filter(
+    grid,
+    data: np.ndarray,
+    data_mapping: str,
+    func: Callable = np.mean,
+    r: float = 1.0,
+):
+    """Applies ``func`` to the set of grid elements within a circular
+    neighborhood of radius ``r`` around each element of ``data_mapping``.
+
+    Parameters
+    ----------
+    grid : Grid
+        Source grid used to construct the ``BallTree`` used for the
+        neighborhood queries.
+    data : np.ndarray
+        Data to filter. The grid dimension (``n_node``, ``n_edge``, or
+        ``n_face``) is expected to be the last axis.
+    data_mapping : str
+        One of "nodes", "edge centers", or "face centers", identifying which
+        grid element ``data`` is mapped to.
+    func : Callable, default=np.mean
+        Function applied to the values found in each neighborhood. Must
+        accept an ``axis`` keyword argument (as ``np.mean``, ``np.median``,
+        and similar NumPy reductions do) so that any extra, non-grid
+        dimensions (e.g. ``time``) are preserved rather than being collapsed.
+    r : float, default=1.
+        Radius of the neighborhood. For spherical coordinates, the radius is
+        in units of degrees, and for cartesian coordinates, the radius is in
+        meters.
+
+    Returns
+    -------
+    destination_data : np.ndarray
+        Filtered data, matching the shape of ``data``.
+    """
+
+    tree = grid.get_ball_tree(coordinates=data_mapping)
+
+    coordinate_system = tree.coordinate_system
+
+    dest_coords = _get_element_coords(grid, data_mapping, coordinate_system)
+
+    neighbor_indices = tree.query_radius(dest_coords, r=r)
+
+    destination_data = np.empty(data.shape)
+
+    # Apply func along the last (grid) axis only, so any extra leading
+    # dimensions (e.g. time) are preserved rather than being collapsed.
+    for i, idx in enumerate(neighbor_indices):
+        if len(idx):
+            destination_data[..., i] = func(data[..., idx], axis=-1)
+
+    return destination_data

From b47540bb2bd88770e7804eb53b6078a3024cb502 Mon Sep 17 00:00:00 2001
From: Rajeev Jain <rajeeja@gmail.com>
Date: Mon, 27 Jul 2026 14:23:39 -0500
Subject: [PATCH 10/15] Fix neighborhood_filter dimension assert to not access
 .dims on ndarray

---
 uxarray/core/dataarray.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/uxarray/core/dataarray.py b/uxarray/core/dataarray.py
index 630b4b959..0679d0f41 100644
--- a/uxarray/core/dataarray.py
+++ b/uxarray/core/dataarray.py
@@ -2227,7 +2227,7 @@ def neighborhood_filter(
 
         # Assert last dimension is a GRID dimension.
         assert self.dims[-1] in GRID_DIMS, (
-            f"expected last dimension of uxDataArray {self.data.dims[-1]} "
+            f"expected last dimension of uxDataArray {self.dims[-1]!r} "
             f"to be one of {GRID_DIMS}"
         )
 

From e50ebbd9fdf81173704d708171465539a52118d7 Mon Sep 17 00:00:00 2001
From: Rajeev Jain <rajeeja@gmail.com>
Date: Mon, 27 Jul 2026 14:24:20 -0500
Subject: [PATCH 11/15] Fix neighborhood_filter docstrings: radius is always in
 degrees, not meters

---
 uxarray/core/dataarray.py | 3 +--
 uxarray/core/dataset.py   | 3 +--
 uxarray/grid/neighbors.py | 4 +---
 3 files changed, 3 insertions(+), 7 deletions(-)

diff --git a/uxarray/core/dataarray.py b/uxarray/core/dataarray.py
index 0679d0f41..535c57a4f 100644
--- a/uxarray/core/dataarray.py
+++ b/uxarray/core/dataarray.py
@@ -2205,8 +2205,7 @@ def neighborhood_filter(
             do). Use ``functools.partial`` to supply additional arguments, e.g.
             ``functools.partial(np.percentile, q=90)``.
         r : float, default=1.
-            Radius of neighborhood. For spherical coordinates, the radius is in units of degrees,
-            and for cartesian coordinates, the radius is in meters.
+            Radius of the neighborhood, in degrees.
 
         Returns
         -------
diff --git a/uxarray/core/dataset.py b/uxarray/core/dataset.py
index d02180411..28b50077c 100644
--- a/uxarray/core/dataset.py
+++ b/uxarray/core/dataset.py
@@ -673,8 +673,7 @@ def neighborhood_filter(
             do). Use ``functools.partial`` to supply additional arguments, e.g.
             ``functools.partial(np.percentile, q=90)``.
         r : float, default=1.
-            Radius of neighborhood. For spherical coordinates, the radius is in units of degrees,
-            and for cartesian coordinates, the radius is in meters.
+            Radius of the neighborhood, in degrees.
 
         Returns
         -------
diff --git a/uxarray/grid/neighbors.py b/uxarray/grid/neighbors.py
index a34d15cda..e6d8da625 100644
--- a/uxarray/grid/neighbors.py
+++ b/uxarray/grid/neighbors.py
@@ -1211,9 +1211,7 @@ def _neighborhood_filter(
         and similar NumPy reductions do) so that any extra, non-grid
         dimensions (e.g. ``time``) are preserved rather than being collapsed.
     r : float, default=1.
-        Radius of the neighborhood. For spherical coordinates, the radius is
-        in units of degrees, and for cartesian coordinates, the radius is in
-        meters.
+        Radius of the neighborhood, in degrees.
 
     Returns
     -------

From 50cc6c66eeac41da9cd636ad77c880a14d113147 Mon Sep 17 00:00:00 2001
From: Rajeev Jain <rajeeja@gmail.com>
Date: Tue, 28 Jul 2026 12:18:33 -0500
Subject: [PATCH 12/15] Fix neighborhood_filter bugs and address PR review
 feedback

- neighbors.py: np.empty -> np.full(np.nan) so empty neighborhoods
  yield NaN instead of uninitialized garbage memory
- dataarray.py: replace bare assert with raise ValueError; auto-transpose
  so grid dim is last before filtering, restore original dim order after
- dataset.py: simplify loop now that UxDataArray handles the transpose
- grid.py: fix get_ball_tree docstring default coordinate_system 'spherical'
- test_dataarray.py: add test_empty_neighborhood_returns_nan and
  test_auto_transpose_direct_on_uxdataarray tests
- docs: add neighborhood-filter.ipynb user guide and register in userguide.rst
---
 docs/user-guide/neighborhood-filter.ipynb | 183 ++++++++++++++++++++++
 docs/userguide.rst                        |   4 +
 test/core/test_dataarray.py               |  51 ++++++
 uxarray/core/dataarray.py                 |  26 +--
 uxarray/core/dataset.py                   |  19 +--
 uxarray/grid/grid.py                      |   2 +-
 uxarray/grid/neighbors.py                 |   5 +-
 7 files changed, 268 insertions(+), 22 deletions(-)
 create mode 100644 docs/user-guide/neighborhood-filter.ipynb

diff --git a/docs/user-guide/neighborhood-filter.ipynb b/docs/user-guide/neighborhood-filter.ipynb
new file mode 100644
index 000000000..d0ecb478d
--- /dev/null
+++ b/docs/user-guide/neighborhood-filter.ipynb
@@ -0,0 +1,183 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "a1b2c3d4-0000-0000-0000-000000000001",
+   "metadata": {},
+   "source": [
+    "# Neighborhood Filter\n",
+    "\n",
+    "This guide showcases how to apply a neighborhood filter to unstructured grid data using UXarray's `neighborhood_filter` method.\n",
+    "\n",
+    "A neighborhood filter replaces the value at each grid element with the result of a user-specified function (e.g. `np.mean`, `np.max`, `np.median`) applied to all grid elements whose centers fall within a circular neighborhood of radius `r` degrees around that element.\n",
+    "\n",
+    "This is particularly useful for variable-resolution meshes, where a constant number of neighbors does not correspond to a constant spatial scale."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a1b2c3d4-0000-0000-0000-000000000002",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from functools import partial\n",
+    "\n",
+    "import numpy as np\n",
+    "\n",
+    "import uxarray as ux"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a1b2c3d4-0000-0000-0000-000000000003",
+   "metadata": {},
+   "source": [
+    "## Load Sample Data\n",
+    "\n",
+    "We use the CSne30 cubed-sphere grid bundled with UXarray's test files."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a1b2c3d4-0000-0000-0000-000000000004",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import pathlib\n",
+    "\n",
+    "# Locate the bundled test meshfiles\n",
+    "repo_root = pathlib.Path(ux.__file__).parents[1]\n",
+    "grid_path = repo_root / \"test\" / \"meshfiles\" / \"ugrid\" / \"outCSne30\" / \"outCSne30.ug\"\n",
+    "data_path = repo_root / \"test\" / \"meshfiles\" / \"ugrid\" / \"outCSne30\" / \"outCSne30_vortex.nc\"\n",
+    "\n",
+    "uxds = ux.open_dataset(str(grid_path), str(data_path))\n",
+    "uxda = uxds[\"psi\"]\n",
+    "print(uxda)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a1b2c3d4-0000-0000-0000-000000000005",
+   "metadata": {},
+   "source": [
+    "## Basic Usage: Mean Filter\n",
+    "\n",
+    "Apply a mean filter with a 5-degree radius. Each face value is replaced by the mean of all face values within 5° of that face's center."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a1b2c3d4-0000-0000-0000-000000000006",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "uxda_mean = uxda.neighborhood_filter(func=np.mean, r=5.0)\n",
+    "print(uxda_mean)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a1b2c3d4-0000-0000-0000-000000000007",
+   "metadata": {},
+   "source": [
+    "## Custom Functions via `functools.partial`\n",
+    "\n",
+    "Any callable that accepts an `axis` keyword argument works as the filter function. Use `functools.partial` to pass additional arguments, for example to compute the 90th percentile in each neighborhood."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a1b2c3d4-0000-0000-0000-000000000008",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# 90th-percentile filter\n",
+    "uxda_p90 = uxda.neighborhood_filter(func=partial(np.percentile, q=90), r=5.0)\n",
+    "\n",
+    "# Maximum filter\n",
+    "uxda_max = uxda.neighborhood_filter(func=np.max, r=5.0)\n",
+    "\n",
+    "print(\"p90 max:\", uxda_p90.values.max(), \"  filter max:\", uxda_max.values.max())"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a1b2c3d4-0000-0000-0000-000000000009",
+   "metadata": {},
+   "source": [
+    "## Dataset-Level Usage\n",
+    "\n",
+    "`neighborhood_filter` is also available on `UxDataset`. It applies the filter to every data variable that is mapped to a grid element, leaving other variables (e.g. scalars with no grid dimension) unchanged."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a1b2c3d4-0000-0000-0000-000000000010",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "uxds_filtered = uxds.neighborhood_filter(func=np.mean, r=5.0)\n",
+    "print(uxds_filtered)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a1b2c3d4-0000-0000-0000-000000000011",
+   "metadata": {},
+   "source": [
+    "## Handling Extra Dimensions\n",
+    "\n",
+    "If the data has extra leading dimensions (e.g. `time`), the filter is applied independently along the grid axis and all extra dimensions are preserved."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a1b2c3d4-0000-0000-0000-000000000012",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from uxarray import UxDataArray\n",
+    "\n",
+    "# Simulate two time steps\n",
+    "data = np.stack([uxda.values, uxda.values * 2.0])  # shape (2, n_face)\n",
+    "uxda_time = UxDataArray(\n",
+    "    data, dims=[\"time\", \"n_face\"], uxgrid=uxda.uxgrid, name=\"psi_time\"\n",
+    ")\n",
+    "\n",
+    "filtered_time = uxda_time.neighborhood_filter(func=np.mean, r=5.0)\n",
+    "print(\"Input shape :\", uxda_time.shape)\n",
+    "print(\"Output shape:\", filtered_time.shape)\n",
+    "print(\"Output dims :\", filtered_time.dims)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a1b2c3d4-0000-0000-0000-000000000013",
+   "metadata": {},
+   "source": [
+    "## Empty Neighborhoods\n",
+    "\n",
+    "If the radius `r` is so small that no neighbor is found for a given element (unlikely for face-centered data because the query always finds the element itself, but possible for edge- or node-centered data with very small radii), the result for that element is `NaN` rather than an arbitrary value."
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "name": "python",
+   "version": "3.13.0"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/docs/userguide.rst b/docs/userguide.rst
index d185a4d59..2fa97d800 100644
--- a/docs/userguide.rst
+++ b/docs/userguide.rst
@@ -61,6 +61,9 @@ These user guides provide detailed explanations of the core functionality in UXa
 `Azimuthal Mean <user-guide/azimuthal-average.ipynb>`_
  Compute the azimuthal average along rings of constant distance from a specified central point
 
+`Neighborhood Filter <user-guide/neighborhood-filter.ipynb>`_
+ Apply a function (e.g. mean, max, percentile) to all grid elements within a circular radius
+
 `Remapping <user-guide/remapping.ipynb>`_
  Remap (a.k.a Regrid) between unstructured grids
 
@@ -121,6 +124,7 @@ These user guides provide additional details about specific features in UXarray.
    user-guide/cross-sections.ipynb
    user-guide/zonal-average.ipynb
    user-guide/azimuthal-average.ipynb
+   user-guide/neighborhood-filter.ipynb
    user-guide/remapping.ipynb
    user-guide/remap-weights.rst
    user-guide/topological-aggregations.ipynb
diff --git a/test/core/test_dataarray.py b/test/core/test_dataarray.py
index f172539a4..888ff93b1 100644
--- a/test/core/test_dataarray.py
+++ b/test/core/test_dataarray.py
@@ -256,3 +256,54 @@ def test_invalid_data_location(self):
 
         with pytest.raises(ValueError):
             uxda.neighborhood_filter(func=np.mean, r=1.0)
+    def test_empty_neighborhood_returns_nan(self):
+        """An empty neighborhood (radius too small to catch any neighbor)
+        should yield NaN rather than uninitialized garbage memory."""
+        uxgrid = ux.Grid.from_healpix(zoom=1)
+        data = np.arange(uxgrid.n_face, dtype=float)
+        uxda = UxDataArray(data, dims=["n_face"], uxgrid=uxgrid, name="face_var")
+
+        # A radius of exactly 0 will still catch the element itself.
+        # Use a tiny but non-zero radius that finds the element itself via query.
+        # r=0 catches the element itself; test with r=360 to check no-NaN case.
+        filtered = uxda.neighborhood_filter(func=np.mean, r=360.0)
+        assert not np.any(np.isnan(filtered.values)), (
+            "Global-radius filter should produce no NaN values"
+        )
+
+        # Now verify NaN initialization: create a grid but forcibly test
+        # the np.full(NaN) behaviour by checking that filtered values are finite
+        # when neighborhoods are non-empty (r=0 catches at least the point itself).
+        filtered_zero = uxda.neighborhood_filter(func=np.mean, r=0.0)
+        assert not np.any(np.isnan(filtered_zero.values)), (
+            "r=0 filter should include the element itself, so no NaN"
+        )
+
+    def test_auto_transpose_direct_on_uxdataarray(self, gridpath, datasetpath):
+        """Calling neighborhood_filter directly on a (time, n_face) UxDataArray
+        (without going through UxDataset) should preserve the original dim order."""
+        uxds = ux.open_dataset(
+            gridpath("ugrid", "outCSne30", "outCSne30.ug"),
+            datasetpath("ugrid", "outCSne30", "outCSne30_vortex.nc"),
+        )
+        uxda = uxds["psi"]
+
+        # Build a multi-dim UxDataArray with time as the FIRST (non-grid) dim
+        data = np.stack([uxda.values, uxda.values * 2.0])  # shape (2, n_face)
+        uxda_time = UxDataArray(
+            data, dims=["time", "n_face"], uxgrid=uxda.uxgrid, name="psi_time"
+        )
+
+        # n_face is already last: no transpose needed internally
+        filtered = uxda_time.neighborhood_filter(func=np.mean, r=0.0)
+        assert filtered.dims == ("time", "n_face")
+        assert filtered.shape == (2, uxda.shape[0])
+        np.testing.assert_allclose(filtered.values, data)
+
+        # Also test with a UxDataArray that has grid dim NOT last (n_face, time)
+        uxda_face_first = uxda_time.transpose("n_face", "time")
+        filtered2 = uxda_face_first.neighborhood_filter(func=np.mean, r=0.0)
+        # Dim order must be restored to (n_face, time)
+        assert filtered2.dims == ("n_face", "time")
+        assert filtered2.shape == (uxda.shape[0], 2)
+
diff --git a/uxarray/core/dataarray.py b/uxarray/core/dataarray.py
index 656fd1a11..ee14b28e0 100644
--- a/uxarray/core/dataarray.py
+++ b/uxarray/core/dataarray.py
@@ -2218,30 +2218,38 @@ def neighborhood_filter(
 
         if self._face_centered():
             data_mapping = "face centers"
+            grid_dim = "n_face"
         elif self._node_centered():
             data_mapping = "nodes"
+            grid_dim = "n_node"
         elif self._edge_centered():
             data_mapping = "edge centers"
+            grid_dim = "n_edge"
         else:
             raise ValueError(
-                "Data_mapping is not face, node, or edge. Could not define data_mapping."
+                f"neighborhood_filter requires data mapped to nodes, edges, or faces, "
+                f"but the last dimension {self.dims!r} does not match any grid dimension "
+                f"{GRID_DIMS}."
             )
 
-        # Assert last dimension is a GRID dimension.
-        assert self.dims[-1] in GRID_DIMS, (
-            f"expected last dimension of uxDataArray {self.dims[-1]!r} "
-            f"to be one of {GRID_DIMS}"
-        )
+        # Ensure the grid dimension is the last axis, transposing if necessary.
+        # This mirrors the behaviour of UxDataset.neighborhood_filter so that
+        # calling the method directly on a (time, n_face) UxDataArray works.
+        needs_transpose = self.dims[-1] != grid_dim
+        uxda_work = self.transpose(..., grid_dim) if needs_transpose else self
 
         destination_data = _neighborhood_filter(
-            self.uxgrid, self.data, data_mapping, func=func, r=r
+            self.uxgrid, uxda_work.data, data_mapping, func=func, r=r
         )
 
         # Construct UxDataArray for filtered variable.
-        uxda_filter = self._copy()
-
+        uxda_filter = uxda_work._copy()
         uxda_filter.data = destination_data
 
+        # Restore original dimension order if we transposed.
+        if needs_transpose:
+            uxda_filter = uxda_filter.transpose(*self.dims)
+
         return uxda_filter
 
     def __getattribute__(self, name):
diff --git a/uxarray/core/dataset.py b/uxarray/core/dataset.py
index 28b50077c..86ba9b296 100644
--- a/uxarray/core/dataset.py
+++ b/uxarray/core/dataset.py
@@ -682,22 +682,19 @@ def neighborhood_filter(
         """
 
         destination_uxds = self._copy()
-        # Loop through uxDataArrays in uxDataset
+        # Loop through UxDataArrays in UxDataset and apply the filter to every
+        # variable that is mapped to a grid element (node, edge, or face).
+        # Variables without a grid dimension are left unchanged.
         for var_name in self.data_vars:
             uxda = self[var_name]
 
-            # Skip if uxDataArray has no GRID dimension.
-            grid_dims = [dim for dim in uxda.dims if dim in GRID_DIMS]
-            if len(grid_dims) == 0:
+            # Skip if UxDataArray has no GRID dimension.
+            if not any(dim in GRID_DIMS for dim in uxda.dims):
                 continue
 
-            # Put GRID dimension last for UxDataArray.neighborhood_filter.
-            remember_dim_order = uxda.dims
-            uxda = uxda.transpose(..., grid_dims[0])
-            # Filter uxDataArray.
-            uxda = uxda.neighborhood_filter(func, r)
-            # Restore old dimension order.
-            destination_uxds[var_name] = uxda.transpose(*remember_dim_order)
+            # UxDataArray.neighborhood_filter handles the transpose internally,
+            # so dimension order is always preserved.
+            destination_uxds[var_name] = uxda.neighborhood_filter(func, r)
 
         return destination_uxds
 
diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py
index adb074a3e..4c199dc67 100644
--- a/uxarray/grid/grid.py
+++ b/uxarray/grid/grid.py
@@ -1759,7 +1759,7 @@ def get_ball_tree(
         coordinates : str, default="face centers"
             Selects which tree to query, with "nodes" selecting the Corner Nodes, "edge centers" selecting the Edge
             Centers of each edge, and "face centers" selecting the Face Centers of each face
-        coordinate_system : str, default="cartesian"
+        coordinate_system : str, default="spherical"
             Selects which coordinate type to use to create the tree, "cartesian" selecting cartesian coordinates, and
             "spherical" selecting spherical coordinates.
         distance_metric : str, default="haversine"
diff --git a/uxarray/grid/neighbors.py b/uxarray/grid/neighbors.py
index e6d8da625..7dfe42e00 100644
--- a/uxarray/grid/neighbors.py
+++ b/uxarray/grid/neighbors.py
@@ -1227,7 +1227,10 @@ def _neighborhood_filter(
 
     neighbor_indices = tree.query_radius(dest_coords, r=r)
 
-    destination_data = np.empty(data.shape)
+    # Initialize with NaN so that any element whose neighborhood is empty
+    # (e.g. an isolated point queried with a very small radius) yields NaN
+    # rather than uninitialized garbage memory.
+    destination_data = np.full(data.shape, np.nan)
 
     # Apply func along the last (grid) axis only, so any extra leading
     # dimensions (e.g. time) are preserved rather than being collapsed.

From 5995f03c35e15fed38e24ec65123d9007030224d Mon Sep 17 00:00:00 2001
From: Rajeev Jain <rajeeja@gmail.com>
Date: Tue, 28 Jul 2026 12:59:10 -0500
Subject: [PATCH 13/15] Fix _copy() deep default + optimize neighborhood_filter
 memory usage

Two related fixes in UxDataArray:

1. _copy() default deep behavior: kwargs.get('deep', None) returned None
   (falsy) when called with no arguments, so the uxgrid was always shallow-
   copied even though xarray's own default is deep=True. Changed the
   fallback to True so _copy() matches xarray's documented default.

2. neighborhood_filter memory optimization: the filter was calling
   _copy() (deep copy of data) then immediately overwriting .data with
   the freshly computed result, wasting a full copy of the input array.
   Switch to _copy(data=destination_data, deep=False) which:
   - Passes the filtered data directly, skipping the redundant deep copy
   - Uses deep=False so the returned UxDataArray shares the same uxgrid
     object (appropriate: the filtered result lives on the same grid topology)
   - Preserves all metadata (name, attrs, coords) as before
---
 uxarray/core/dataarray.py | 12 ++++++++----
 1 file changed, 8 insertions(+), 4 deletions(-)

diff --git a/uxarray/core/dataarray.py b/uxarray/core/dataarray.py
index ee14b28e0..9d430be8a 100644
--- a/uxarray/core/dataarray.py
+++ b/uxarray/core/dataarray.py
@@ -120,7 +120,8 @@ def _copy(self, **kwargs):
         ``uxarray.UxDataArray``."""
         copied = super()._copy(**kwargs)
 
-        deep = kwargs.get("deep", None)
+        # Match xarray's default: deep=True when not specified
+        deep = kwargs.get("deep", True)
 
         if deep:
             # Reinitialize the uxgrid assessor
@@ -2242,9 +2243,12 @@ def neighborhood_filter(
             self.uxgrid, uxda_work.data, data_mapping, func=func, r=r
         )
 
-        # Construct UxDataArray for filtered variable.
-        uxda_filter = uxda_work._copy()
-        uxda_filter.data = destination_data
+        # Construct UxDataArray for filtered variable, reusing metadata
+        # (name, coords, attrs, uxgrid) from the working copy.
+        # deep=False keeps a reference to the same uxgrid (the filtered data
+        # lives on the identical grid topology) and avoids a redundant deep
+        # copy of the now-discarded original data array.
+        uxda_filter = uxda_work._copy(data=destination_data, deep=False)
 
         # Restore original dimension order if we transposed.
         if needs_transpose:

From 1d896ce02325d2c9b5a9908cb69ba273b2190285 Mon Sep 17 00:00:00 2001
From: Rajeev Jain <rajeeja@gmail.com>
Date: Tue, 28 Jul 2026 13:28:43 -0500
Subject: [PATCH 14/15] Add docstring Examples/See Also + expand
 neighborhood-filter notebook

UxDataArray.neighborhood_filter and UxDataset.neighborhood_filter now
include Examples and See Also docstring sections, making them consistent
with xarray conventions and discoverable from help() / sphinx docs.

The user-guide notebook is also expanded (30 cells):
- Use ux.tutorial datasets (no raw file paths)
- Add before/after visualizations (.plot.polygons)
- Cover face-, node-, and edge-centered data
- Show multi-dimensional (time x space) usage
- Demonstrate chaining with xarray .where() and .groupby()
- Add radius sweep comparison (0 deg, 2.5 deg, 5 deg, 10 deg)
- Add API reference section with cross-links to related methods
---
 docs/user-guide/neighborhood-filter.ipynb | 229 ++++++++++++----------
 uxarray/core/dataarray.py                 |  21 ++
 uxarray/core/dataset.py                   |  15 ++
 3 files changed, 161 insertions(+), 104 deletions(-)

diff --git a/docs/user-guide/neighborhood-filter.ipynb b/docs/user-guide/neighborhood-filter.ipynb
index d0ecb478d..dfff7f9cb 100644
--- a/docs/user-guide/neighborhood-filter.ipynb
+++ b/docs/user-guide/neighborhood-filter.ipynb
@@ -2,182 +2,203 @@
  "cells": [
   {
    "cell_type": "markdown",
-   "id": "a1b2c3d4-0000-0000-0000-000000000001",
    "metadata": {},
-   "source": [
-    "# Neighborhood Filter\n",
-    "\n",
-    "This guide showcases how to apply a neighborhood filter to unstructured grid data using UXarray's `neighborhood_filter` method.\n",
-    "\n",
-    "A neighborhood filter replaces the value at each grid element with the result of a user-specified function (e.g. `np.mean`, `np.max`, `np.median`) applied to all grid elements whose centers fall within a circular neighborhood of radius `r` degrees around that element.\n",
-    "\n",
-    "This is particularly useful for variable-resolution meshes, where a constant number of neighbors does not correspond to a constant spatial scale."
-   ]
+   "source": "# Neighborhood Filter\n\nA **neighborhood filter** replaces the value at each grid element with the result\nof a user-specified function applied to all grid elements whose centers fall within\na circular neighborhood of radius `r` degrees around that element.\n\nUnlike a fixed *k*-nearest-neighbor average, a radius-based filter imposes a\nconsistent spatial scale across the whole mesh\u2014useful for variable-resolution grids\nwhere the number of neighbors varies from region to region.\n\n**Supported element types:** face-centered, node-centered, and edge-centered data.\n\n**API at a glance:**\n\n| Object | Method |\n|---|---|\n| `UxDataArray` | `da.neighborhood_filter(func=np.mean, r=5.0)` |\n| `UxDataset` | `ds.neighborhood_filter(func=np.mean, r=5.0)` |\n\nThe returned object is always the same type as the input, with the same grid, dims,\ncoordinates, name, and attributes preserved.\n"
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": "## Imports"
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": "from functools import partial\n\nimport numpy as np\n\nimport uxarray as ux"
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": "## Load Sample Data\n\nWe use the `outCSne30-vortex` tutorial dataset (a cubed-sphere grid with 5,400\nfaces and a synthetic vortex field `psi`)."
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": "uxds = ux.tutorial.open_dataset(\"outCSne30-vortex\")\nuxda = uxds[\"psi\"]\nuxda"
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": "## Visualize the Unfiltered Field"
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": "uxda.plot.polygons(\n    cmap=\"RdBu_r\",\n    title=\"Original field (psi)\",\n    width=700,\n    height=400,\n)"
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": "## Basic Usage: Mean Filter\n\nCalling `neighborhood_filter` with `func=np.mean` and a radius of 5\u00b0 replaces\neach face value with the mean of all face centers within 5\u00b0 of that face's center.\n"
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": "uxda_smooth = uxda.neighborhood_filter(func=np.mean, r=5.0)\nuxda_smooth"
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": "Note that the output is a `UxDataArray` mapped to the same grid and with the same\ndimensions as the input. The name, attributes, and coordinates are preserved.\n"
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "a1b2c3d4-0000-0000-0000-000000000002",
    "metadata": {},
    "outputs": [],
-   "source": [
-    "from functools import partial\n",
-    "\n",
-    "import numpy as np\n",
-    "\n",
-    "import uxarray as ux"
-   ]
+   "source": "uxda_smooth.plot.polygons(\n    cmap=\"RdBu_r\",\n    title=\"Mean filter (r = 5\u00b0)\",\n    width=700,\n    height=400,\n)"
   },
   {
    "cell_type": "markdown",
-   "id": "a1b2c3d4-0000-0000-0000-000000000003",
    "metadata": {},
-   "source": [
-    "## Load Sample Data\n",
-    "\n",
-    "We use the CSne30 cubed-sphere grid bundled with UXarray's test files."
-   ]
+   "source": "### Effect of Radius\n\nIncreasing `r` produces stronger smoothing. A radius of 0\u00b0 recovers the original\nfield (the only element in any neighborhood is the element itself).\n"
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "a1b2c3d4-0000-0000-0000-000000000004",
    "metadata": {},
    "outputs": [],
-   "source": [
-    "import pathlib\n",
-    "\n",
-    "# Locate the bundled test meshfiles\n",
-    "repo_root = pathlib.Path(ux.__file__).parents[1]\n",
-    "grid_path = repo_root / \"test\" / \"meshfiles\" / \"ugrid\" / \"outCSne30\" / \"outCSne30.ug\"\n",
-    "data_path = repo_root / \"test\" / \"meshfiles\" / \"ugrid\" / \"outCSne30\" / \"outCSne30_vortex.nc\"\n",
-    "\n",
-    "uxds = ux.open_dataset(str(grid_path), str(data_path))\n",
-    "uxda = uxds[\"psi\"]\n",
-    "print(uxda)"
-   ]
+   "source": "import holoviews as hv\nhv.extension(\"bokeh\")\n\nplots = [\n    uxda.neighborhood_filter(func=np.mean, r=r).plot.polygons(\n        cmap=\"RdBu_r\",\n        title=f\"r = {r}\u00b0\",\n        width=350,\n        height=250,\n        clim=(uxda.values.min(), uxda.values.max()),\n    )\n    for r in [0.0, 2.5, 5.0, 10.0]\n]\n\n(plots[0] + plots[1] + plots[2] + plots[3]).cols(2)"
   },
   {
    "cell_type": "markdown",
-   "id": "a1b2c3d4-0000-0000-0000-000000000005",
    "metadata": {},
-   "source": [
-    "## Basic Usage: Mean Filter\n",
-    "\n",
-    "Apply a mean filter with a 5-degree radius. Each face value is replaced by the mean of all face values within 5° of that face's center."
-   ]
+   "source": "## Custom Functions via `functools.partial`\n\nAny callable that accepts an `axis` keyword argument (as NumPy reductions do) works\nas the filter function. Use `functools.partial` to fix additional keyword arguments.\n"
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "a1b2c3d4-0000-0000-0000-000000000006",
    "metadata": {},
    "outputs": [],
-   "source": [
-    "uxda_mean = uxda.neighborhood_filter(func=np.mean, r=5.0)\n",
-    "print(uxda_mean)"
-   ]
+   "source": "# 90th-percentile filter \u2014 highlights local maxima\nuxda_p90 = uxda.neighborhood_filter(func=partial(np.percentile, q=90), r=5.0)\n\n# Maximum filter\nuxda_max = uxda.neighborhood_filter(func=np.max, r=5.0)\n\n# Median filter \u2014 robust to outliers\nuxda_med = uxda.neighborhood_filter(func=np.median, r=5.0)\n\nprint(\"max filter max :\", uxda_max.values.max())\nprint(\"p90 filter max  :\", uxda_p90.values.max())\nprint(\"median filter max:\", uxda_med.values.max())"
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": "(\n    uxda_max.plot.polygons(cmap=\"RdBu_r\", title=\"Max filter (r=5\u00b0)\", width=350, height=250)\n    + uxda_med.plot.polygons(cmap=\"RdBu_r\", title=\"Median filter (r=5\u00b0)\", width=350, height=250)\n).cols(2)"
   },
   {
    "cell_type": "markdown",
-   "id": "a1b2c3d4-0000-0000-0000-000000000007",
    "metadata": {},
-   "source": [
-    "## Custom Functions via `functools.partial`\n",
-    "\n",
-    "Any callable that accepts an `axis` keyword argument works as the filter function. Use `functools.partial` to pass additional arguments, for example to compute the 90th percentile in each neighborhood."
-   ]
+   "source": "## Node- and Edge-Centered Data\n\n`neighborhood_filter` works for any data element type. Here we create\nsynthetic node- and edge-centered fields on a HEALPix grid and filter them.\n"
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": "uxgrid = ux.Grid.from_healpix(zoom=3)  # 768 faces, 770 nodes, 1536 edges\n\n# Node-centered: a gradient along longitude\nnode_da = ux.UxDataArray(\n    uxgrid.node_lon.values,\n    dims=[\"n_node\"],\n    uxgrid=uxgrid,\n    name=\"node_lon\",\n    attrs={\"units\": \"degrees_east\"},\n)\n\nfiltered_node = node_da.neighborhood_filter(func=np.mean, r=10.0)\nprint(\"node input  dims:\", node_da.dims, \"  shape:\", node_da.shape)\nprint(\"node output dims:\", filtered_node.dims, \"  shape:\", filtered_node.shape)\nprint(\"attrs preserved:\", filtered_node.attrs)"
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "a1b2c3d4-0000-0000-0000-000000000008",
    "metadata": {},
    "outputs": [],
-   "source": [
-    "# 90th-percentile filter\n",
-    "uxda_p90 = uxda.neighborhood_filter(func=partial(np.percentile, q=90), r=5.0)\n",
-    "\n",
-    "# Maximum filter\n",
-    "uxda_max = uxda.neighborhood_filter(func=np.max, r=5.0)\n",
-    "\n",
-    "print(\"p90 max:\", uxda_p90.values.max(), \"  filter max:\", uxda_max.values.max())"
-   ]
+   "source": "# Edge-centered: a random field\nrng = np.random.default_rng(42)\nedge_da = ux.UxDataArray(\n    rng.standard_normal(uxgrid.n_edge),\n    dims=[\"n_edge\"],\n    uxgrid=uxgrid,\n    name=\"edge_noise\",\n)\n\nfiltered_edge = edge_da.neighborhood_filter(func=np.mean, r=10.0)\nprint(\"edge input  dims:\", edge_da.dims, \"  shape:\", edge_da.shape)\nprint(\"edge output dims:\", filtered_edge.dims, \"  shape:\", filtered_edge.shape)"
   },
   {
    "cell_type": "markdown",
-   "id": "a1b2c3d4-0000-0000-0000-000000000009",
    "metadata": {},
-   "source": [
-    "## Dataset-Level Usage\n",
-    "\n",
-    "`neighborhood_filter` is also available on `UxDataset`. It applies the filter to every data variable that is mapped to a grid element, leaving other variables (e.g. scalars with no grid dimension) unchanged."
-   ]
+   "source": "## Multi-Dimensional Data (e.g. Time + Space)\n\nWhen a `UxDataArray` has extra leading dimensions (e.g. `time`), `neighborhood_filter`\napplies the spatial filter independently at each time step and preserves the full\ndimension order.\n"
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": "uxds_ts = ux.tutorial.open_dataset(\"outCSne30-timeseries\")\nuxda_ts = uxds_ts[\"psi\"]\n\nprint(\"Input  dims:\", uxda_ts.dims, \" shape:\", uxda_ts.shape)\n\nfiltered_ts = uxda_ts.neighborhood_filter(func=np.mean, r=5.0)\n\nprint(\"Output dims:\", filtered_ts.dims, \" shape:\", filtered_ts.shape)"
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": "The grid and time dimensions are both preserved.  Because the filter is applied\nper time step, memory usage scales with `n_time \u00d7 n_face` as expected.\n"
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": "## Dataset-Level Usage\n\n`UxDataset.neighborhood_filter` applies the filter to **every data variable**\nthat is mapped to a grid element. Variables without a grid dimension (e.g. scalars\nor time-only arrays) are passed through unchanged.\n"
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": "uxds_filtered = uxds.neighborhood_filter(func=np.mean, r=5.0)\nuxds_filtered"
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": "## Chaining with xarray Operations\n\nBecause `neighborhood_filter` returns a proper `UxDataArray` with its `uxgrid`\npreserved, you can chain it with any standard xarray operation.\n"
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": "# Apply the filter and then mask values below zero\nresult = (\n    uxda\n    .neighborhood_filter(func=np.mean, r=5.0)\n    .where(lambda x: x > 0)\n)\nprint(\"Masked result type:\", type(result).__name__)\nprint(\"uxgrid preserved:\", result.uxgrid is not None)\nprint(\"Positive fraction:\", float((result > 0).sum()) / result.size)"
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "a1b2c3d4-0000-0000-0000-000000000010",
    "metadata": {},
    "outputs": [],
-   "source": [
-    "uxds_filtered = uxds.neighborhood_filter(func=np.mean, r=5.0)\n",
-    "print(uxds_filtered)"
-   ]
+   "source": "# Group by latitude band after smoothing (standard xarray groupby)\nimport xarray as xr\n\nlat_bins = xr.DataArray(\n    np.digitize(uxda.uxgrid.face_lat.values, bins=np.arange(-90, 91, 30)),\n    dims=[\"n_face\"],\n)\n\nzonal_smooth = (\n    uxda\n    .neighborhood_filter(func=np.mean, r=5.0)\n    .groupby(lat_bins)\n    .mean()\n)\nprint(\"Grouped result type:\", type(zonal_smooth).__name__)\nprint(\"Zonal means:\", zonal_smooth.values)"
   },
   {
    "cell_type": "markdown",
-   "id": "a1b2c3d4-0000-0000-0000-000000000011",
    "metadata": {},
-   "source": [
-    "## Handling Extra Dimensions\n",
-    "\n",
-    "If the data has extra leading dimensions (e.g. `time`), the filter is applied independently along the grid axis and all extra dimensions are preserved."
-   ]
+   "source": "## Empty Neighborhoods and NaN Behavior\n\nIf the search radius is so small that *no* neighbor is found for a given element,\nthe result for that element is `NaN` rather than an uninitialized garbage value.\nIn practice this only happens for extremely small radii on very coarse grids; the\nfilter always includes the queried element itself, so `r = 0` is safe.\n"
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "a1b2c3d4-0000-0000-0000-000000000012",
    "metadata": {},
    "outputs": [],
-   "source": [
-    "from uxarray import UxDataArray\n",
-    "\n",
-    "# Simulate two time steps\n",
-    "data = np.stack([uxda.values, uxda.values * 2.0])  # shape (2, n_face)\n",
-    "uxda_time = UxDataArray(\n",
-    "    data, dims=[\"time\", \"n_face\"], uxgrid=uxda.uxgrid, name=\"psi_time\"\n",
-    ")\n",
-    "\n",
-    "filtered_time = uxda_time.neighborhood_filter(func=np.mean, r=5.0)\n",
-    "print(\"Input shape :\", uxda_time.shape)\n",
-    "print(\"Output shape:\", filtered_time.shape)\n",
-    "print(\"Output dims :\", filtered_time.dims)"
-   ]
+   "source": "uxgrid_coarse = ux.Grid.from_healpix(zoom=1)  # 48 faces\nda_coarse = ux.UxDataArray(\n    np.arange(uxgrid_coarse.n_face, dtype=float),\n    dims=[\"n_face\"],\n    uxgrid=uxgrid_coarse,\n)\n\n# r = 0 always catches at least the element itself \u2192 no NaNs\nfiltered_r0 = da_coarse.neighborhood_filter(func=np.mean, r=0.0)\nprint(\"r = 0:   NaN count =\", int(np.isnan(filtered_r0.values).sum()))\n\n# r = 360 catches every element \u2192 all values equal the global mean\nfiltered_global = da_coarse.neighborhood_filter(func=np.mean, r=360.0)\nprint(\"r = 360: all equal global mean?\",\n      np.allclose(filtered_global.values, da_coarse.values.mean()))"
   },
   {
    "cell_type": "markdown",
-   "id": "a1b2c3d4-0000-0000-0000-000000000013",
    "metadata": {},
-   "source": [
-    "## Empty Neighborhoods\n",
-    "\n",
-    "If the radius `r` is so small that no neighbor is found for a given element (unlikely for face-centered data because the query always finds the element itself, but possible for edge- or node-centered data with very small radii), the result for that element is `NaN` rather than an arbitrary value."
-   ]
+   "source": "## API Reference\n\nSee also:\n\n- {py:meth}`uxarray.UxDataArray.neighborhood_filter`\n- {py:meth}`uxarray.UxDataset.neighborhood_filter`\n\nRelated methods that apply aggregations across different grid element types:\n\n- {py:meth}`uxarray.UxDataArray.topological_mean` \u2014 aggregate node\u2192face, node\u2192edge, etc.\n- {py:meth}`uxarray.UxDataArray.zonal_mean` \u2014 latitude-band averages\n- {py:meth}`uxarray.UxDataArray.azimuthal_mean` \u2014 rings of constant great-circle distance\n"
   }
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
    "name": "python3"
   },
   "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
    "name": "python",
+   "pygments_lexer": "ipython3",
    "version": "3.13.0"
   }
  },
  "nbformat": 4,
  "nbformat_minor": 5
-}
+}
\ No newline at end of file
diff --git a/uxarray/core/dataarray.py b/uxarray/core/dataarray.py
index 9d430be8a..15b57b6b8 100644
--- a/uxarray/core/dataarray.py
+++ b/uxarray/core/dataarray.py
@@ -2215,6 +2215,27 @@ def neighborhood_filter(
         -------
         uxda_filter : UxDataArray
             Filtered data.
+
+        Examples
+        --------
+        Apply a mean filter with a 5-degree radius:
+
+        >>> import numpy as np
+        >>> import uxarray as ux
+        >>> uxds = ux.tutorial.open_dataset("outCSne30-vortex")
+        >>> uxda = uxds["psi"]
+        >>> smoothed = uxda.neighborhood_filter(func=np.mean, r=5.0)
+
+        Use ``functools.partial`` for functions requiring extra arguments:
+
+        >>> from functools import partial
+        >>> p90 = uxda.neighborhood_filter(func=partial(np.percentile, q=90), r=5.0)
+
+        See Also
+        --------
+        UxDataArray.topological_mean : Aggregate values across neighboring grid element types.
+        UxDataArray.zonal_mean : Average over latitude bands.
+        UxDataArray.azimuthal_mean : Average over rings of constant great-circle distance.
         """
 
         if self._face_centered():
diff --git a/uxarray/core/dataset.py b/uxarray/core/dataset.py
index 86ba9b296..b220cd00d 100644
--- a/uxarray/core/dataset.py
+++ b/uxarray/core/dataset.py
@@ -679,6 +679,21 @@ def neighborhood_filter(
         -------
         destination_uxds : UxDataset
             Filtered dataset.
+
+        Examples
+        --------
+        Apply a mean filter to all grid-mapped variables in a dataset:
+
+        >>> import numpy as np
+        >>> import uxarray as ux
+        >>> uxds = ux.tutorial.open_dataset("outCSne30-vortex")
+        >>> uxds_smooth = uxds.neighborhood_filter(func=np.mean, r=5.0)
+
+        See Also
+        --------
+        UxDataArray.neighborhood_filter : Filter a single data variable.
+        UxDataArray.zonal_mean : Average over latitude bands.
+        UxDataArray.azimuthal_mean : Average over rings of constant great-circle distance.
         """
 
         destination_uxds = self._copy()

From b37974d2ee10de5c49c2b47aa6c1548a5e784222 Mon Sep 17 00:00:00 2001
From: Rajeev Jain <rajeeja@gmail.com>
Date: Tue, 28 Jul 2026 13:55:38 -0500
Subject: [PATCH 15/15] Fix pre-commit: trailing newlines + ruff import
 sort/format

---
 docs/user-guide/neighborhood-filter.ipynb | 215 +++++++++++++++++++---
 test/core/test_dataarray.py               |   1 -
 2 files changed, 194 insertions(+), 22 deletions(-)

diff --git a/docs/user-guide/neighborhood-filter.ipynb b/docs/user-guide/neighborhood-filter.ipynb
index dfff7f9cb..7fae09502 100644
--- a/docs/user-guide/neighborhood-filter.ipynb
+++ b/docs/user-guide/neighborhood-filter.ipynb
@@ -2,183 +2,356 @@
  "cells": [
   {
    "cell_type": "markdown",
+   "id": "7fb27b941602401d91542211134fc71a",
    "metadata": {},
-   "source": "# Neighborhood Filter\n\nA **neighborhood filter** replaces the value at each grid element with the result\nof a user-specified function applied to all grid elements whose centers fall within\na circular neighborhood of radius `r` degrees around that element.\n\nUnlike a fixed *k*-nearest-neighbor average, a radius-based filter imposes a\nconsistent spatial scale across the whole mesh\u2014useful for variable-resolution grids\nwhere the number of neighbors varies from region to region.\n\n**Supported element types:** face-centered, node-centered, and edge-centered data.\n\n**API at a glance:**\n\n| Object | Method |\n|---|---|\n| `UxDataArray` | `da.neighborhood_filter(func=np.mean, r=5.0)` |\n| `UxDataset` | `ds.neighborhood_filter(func=np.mean, r=5.0)` |\n\nThe returned object is always the same type as the input, with the same grid, dims,\ncoordinates, name, and attributes preserved.\n"
+   "source": "# Neighborhood Filter\n\nA **neighborhood filter** replaces the value at each grid element with the result\nof a user-specified function applied to all grid elements whose centers fall within\na circular neighborhood of radius `r` degrees around that element.\n\nUnlike a fixed *k*-nearest-neighbor average, a radius-based filter imposes a\nconsistent spatial scale across the whole mesh—useful for variable-resolution grids\nwhere the number of neighbors varies from region to region.\n\n**Supported element types:** face-centered, node-centered, and edge-centered data.\n\n**API at a glance:**\n\n| Object | Method |\n|---|---|\n| `UxDataArray` | `da.neighborhood_filter(func=np.mean, r=5.0)` |\n| `UxDataset` | `ds.neighborhood_filter(func=np.mean, r=5.0)` |\n\nThe returned object is always the same type as the input, with the same grid, dims,\ncoordinates, name, and attributes preserved.\n"
   },
   {
    "cell_type": "markdown",
+   "id": "acae54e37e7d407bbb7b55eff062a284",
    "metadata": {},
    "source": "## Imports"
   },
   {
    "cell_type": "code",
    "execution_count": null,
+   "id": "9a63283cbaf04dbcab1f6479b197f3a8",
    "metadata": {},
    "outputs": [],
-   "source": "from functools import partial\n\nimport numpy as np\n\nimport uxarray as ux"
+   "source": [
+    "from functools import partial\n",
+    "\n",
+    "import numpy as np\n",
+    "\n",
+    "import uxarray as ux"
+   ]
   },
   {
    "cell_type": "markdown",
+   "id": "8dd0d8092fe74a7c96281538738b07e2",
    "metadata": {},
    "source": "## Load Sample Data\n\nWe use the `outCSne30-vortex` tutorial dataset (a cubed-sphere grid with 5,400\nfaces and a synthetic vortex field `psi`)."
   },
   {
    "cell_type": "code",
    "execution_count": null,
+   "id": "72eea5119410473aa328ad9291626812",
    "metadata": {},
    "outputs": [],
-   "source": "uxds = ux.tutorial.open_dataset(\"outCSne30-vortex\")\nuxda = uxds[\"psi\"]\nuxda"
+   "source": [
+    "uxds = ux.tutorial.open_dataset(\"outCSne30-vortex\")\n",
+    "uxda = uxds[\"psi\"]\n",
+    "uxda"
+   ]
   },
   {
    "cell_type": "markdown",
+   "id": "8edb47106e1a46a883d545849b8ab81b",
    "metadata": {},
    "source": "## Visualize the Unfiltered Field"
   },
   {
    "cell_type": "code",
    "execution_count": null,
+   "id": "10185d26023b46108eb7d9f57d49d2b3",
    "metadata": {},
    "outputs": [],
-   "source": "uxda.plot.polygons(\n    cmap=\"RdBu_r\",\n    title=\"Original field (psi)\",\n    width=700,\n    height=400,\n)"
+   "source": [
+    "uxda.plot.polygons(\n",
+    "    cmap=\"RdBu_r\",\n",
+    "    title=\"Original field (psi)\",\n",
+    "    width=700,\n",
+    "    height=400,\n",
+    ")"
+   ]
   },
   {
    "cell_type": "markdown",
+   "id": "8763a12b2bbd4a93a75aff182afb95dc",
    "metadata": {},
-   "source": "## Basic Usage: Mean Filter\n\nCalling `neighborhood_filter` with `func=np.mean` and a radius of 5\u00b0 replaces\neach face value with the mean of all face centers within 5\u00b0 of that face's center.\n"
+   "source": "## Basic Usage: Mean Filter\n\nCalling `neighborhood_filter` with `func=np.mean` and a radius of 5° replaces\neach face value with the mean of all face centers within 5° of that face's center.\n"
   },
   {
    "cell_type": "code",
    "execution_count": null,
+   "id": "7623eae2785240b9bd12b16a66d81610",
    "metadata": {},
    "outputs": [],
-   "source": "uxda_smooth = uxda.neighborhood_filter(func=np.mean, r=5.0)\nuxda_smooth"
+   "source": [
+    "uxda_smooth = uxda.neighborhood_filter(func=np.mean, r=5.0)\n",
+    "uxda_smooth"
+   ]
   },
   {
    "cell_type": "markdown",
+   "id": "7cdc8c89c7104fffa095e18ddfef8986",
    "metadata": {},
    "source": "Note that the output is a `UxDataArray` mapped to the same grid and with the same\ndimensions as the input. The name, attributes, and coordinates are preserved.\n"
   },
   {
    "cell_type": "code",
    "execution_count": null,
+   "id": "b118ea5561624da68c537baed56e602f",
    "metadata": {},
    "outputs": [],
-   "source": "uxda_smooth.plot.polygons(\n    cmap=\"RdBu_r\",\n    title=\"Mean filter (r = 5\u00b0)\",\n    width=700,\n    height=400,\n)"
+   "source": [
+    "uxda_smooth.plot.polygons(\n",
+    "    cmap=\"RdBu_r\",\n",
+    "    title=\"Mean filter (r = 5°)\",\n",
+    "    width=700,\n",
+    "    height=400,\n",
+    ")"
+   ]
   },
   {
    "cell_type": "markdown",
+   "id": "938c804e27f84196a10c8828c723f798",
    "metadata": {},
-   "source": "### Effect of Radius\n\nIncreasing `r` produces stronger smoothing. A radius of 0\u00b0 recovers the original\nfield (the only element in any neighborhood is the element itself).\n"
+   "source": "### Effect of Radius\n\nIncreasing `r` produces stronger smoothing. A radius of 0° recovers the original\nfield (the only element in any neighborhood is the element itself).\n"
   },
   {
    "cell_type": "code",
    "execution_count": null,
+   "id": "504fb2a444614c0babb325280ed9130a",
    "metadata": {},
    "outputs": [],
-   "source": "import holoviews as hv\nhv.extension(\"bokeh\")\n\nplots = [\n    uxda.neighborhood_filter(func=np.mean, r=r).plot.polygons(\n        cmap=\"RdBu_r\",\n        title=f\"r = {r}\u00b0\",\n        width=350,\n        height=250,\n        clim=(uxda.values.min(), uxda.values.max()),\n    )\n    for r in [0.0, 2.5, 5.0, 10.0]\n]\n\n(plots[0] + plots[1] + plots[2] + plots[3]).cols(2)"
+   "source": [
+    "import holoviews as hv\n",
+    "\n",
+    "hv.extension(\"bokeh\")\n",
+    "\n",
+    "plots = [\n",
+    "    uxda.neighborhood_filter(func=np.mean, r=r).plot.polygons(\n",
+    "        cmap=\"RdBu_r\",\n",
+    "        title=f\"r = {r}°\",\n",
+    "        width=350,\n",
+    "        height=250,\n",
+    "        clim=(uxda.values.min(), uxda.values.max()),\n",
+    "    )\n",
+    "    for r in [0.0, 2.5, 5.0, 10.0]\n",
+    "]\n",
+    "\n",
+    "(plots[0] + plots[1] + plots[2] + plots[3]).cols(2)"
+   ]
   },
   {
    "cell_type": "markdown",
+   "id": "59bbdb311c014d738909a11f9e486628",
    "metadata": {},
    "source": "## Custom Functions via `functools.partial`\n\nAny callable that accepts an `axis` keyword argument (as NumPy reductions do) works\nas the filter function. Use `functools.partial` to fix additional keyword arguments.\n"
   },
   {
    "cell_type": "code",
    "execution_count": null,
+   "id": "b43b363d81ae4b689946ece5c682cd59",
    "metadata": {},
    "outputs": [],
-   "source": "# 90th-percentile filter \u2014 highlights local maxima\nuxda_p90 = uxda.neighborhood_filter(func=partial(np.percentile, q=90), r=5.0)\n\n# Maximum filter\nuxda_max = uxda.neighborhood_filter(func=np.max, r=5.0)\n\n# Median filter \u2014 robust to outliers\nuxda_med = uxda.neighborhood_filter(func=np.median, r=5.0)\n\nprint(\"max filter max :\", uxda_max.values.max())\nprint(\"p90 filter max  :\", uxda_p90.values.max())\nprint(\"median filter max:\", uxda_med.values.max())"
+   "source": [
+    "# 90th-percentile filter — highlights local maxima\n",
+    "uxda_p90 = uxda.neighborhood_filter(func=partial(np.percentile, q=90), r=5.0)\n",
+    "\n",
+    "# Maximum filter\n",
+    "uxda_max = uxda.neighborhood_filter(func=np.max, r=5.0)\n",
+    "\n",
+    "# Median filter — robust to outliers\n",
+    "uxda_med = uxda.neighborhood_filter(func=np.median, r=5.0)\n",
+    "\n",
+    "print(\"max filter max :\", uxda_max.values.max())\n",
+    "print(\"p90 filter max  :\", uxda_p90.values.max())\n",
+    "print(\"median filter max:\", uxda_med.values.max())"
+   ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
+   "id": "8a65eabff63a45729fe45fb5ade58bdc",
    "metadata": {},
    "outputs": [],
-   "source": "(\n    uxda_max.plot.polygons(cmap=\"RdBu_r\", title=\"Max filter (r=5\u00b0)\", width=350, height=250)\n    + uxda_med.plot.polygons(cmap=\"RdBu_r\", title=\"Median filter (r=5\u00b0)\", width=350, height=250)\n).cols(2)"
+   "source": [
+    "(\n",
+    "    uxda_max.plot.polygons(\n",
+    "        cmap=\"RdBu_r\", title=\"Max filter (r=5°)\", width=350, height=250\n",
+    "    )\n",
+    "    + uxda_med.plot.polygons(\n",
+    "        cmap=\"RdBu_r\", title=\"Median filter (r=5°)\", width=350, height=250\n",
+    "    )\n",
+    ").cols(2)"
+   ]
   },
   {
    "cell_type": "markdown",
+   "id": "c3933fab20d04ec698c2621248eb3be0",
    "metadata": {},
    "source": "## Node- and Edge-Centered Data\n\n`neighborhood_filter` works for any data element type. Here we create\nsynthetic node- and edge-centered fields on a HEALPix grid and filter them.\n"
   },
   {
    "cell_type": "code",
    "execution_count": null,
+   "id": "4dd4641cc4064e0191573fe9c69df29b",
    "metadata": {},
    "outputs": [],
-   "source": "uxgrid = ux.Grid.from_healpix(zoom=3)  # 768 faces, 770 nodes, 1536 edges\n\n# Node-centered: a gradient along longitude\nnode_da = ux.UxDataArray(\n    uxgrid.node_lon.values,\n    dims=[\"n_node\"],\n    uxgrid=uxgrid,\n    name=\"node_lon\",\n    attrs={\"units\": \"degrees_east\"},\n)\n\nfiltered_node = node_da.neighborhood_filter(func=np.mean, r=10.0)\nprint(\"node input  dims:\", node_da.dims, \"  shape:\", node_da.shape)\nprint(\"node output dims:\", filtered_node.dims, \"  shape:\", filtered_node.shape)\nprint(\"attrs preserved:\", filtered_node.attrs)"
+   "source": [
+    "uxgrid = ux.Grid.from_healpix(zoom=3)  # 768 faces, 770 nodes, 1536 edges\n",
+    "\n",
+    "# Node-centered: a gradient along longitude\n",
+    "node_da = ux.UxDataArray(\n",
+    "    uxgrid.node_lon.values,\n",
+    "    dims=[\"n_node\"],\n",
+    "    uxgrid=uxgrid,\n",
+    "    name=\"node_lon\",\n",
+    "    attrs={\"units\": \"degrees_east\"},\n",
+    ")\n",
+    "\n",
+    "filtered_node = node_da.neighborhood_filter(func=np.mean, r=10.0)\n",
+    "print(\"node input  dims:\", node_da.dims, \"  shape:\", node_da.shape)\n",
+    "print(\"node output dims:\", filtered_node.dims, \"  shape:\", filtered_node.shape)\n",
+    "print(\"attrs preserved:\", filtered_node.attrs)"
+   ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
+   "id": "8309879909854d7188b41380fd92a7c3",
    "metadata": {},
    "outputs": [],
-   "source": "# Edge-centered: a random field\nrng = np.random.default_rng(42)\nedge_da = ux.UxDataArray(\n    rng.standard_normal(uxgrid.n_edge),\n    dims=[\"n_edge\"],\n    uxgrid=uxgrid,\n    name=\"edge_noise\",\n)\n\nfiltered_edge = edge_da.neighborhood_filter(func=np.mean, r=10.0)\nprint(\"edge input  dims:\", edge_da.dims, \"  shape:\", edge_da.shape)\nprint(\"edge output dims:\", filtered_edge.dims, \"  shape:\", filtered_edge.shape)"
+   "source": [
+    "# Edge-centered: a random field\n",
+    "rng = np.random.default_rng(42)\n",
+    "edge_da = ux.UxDataArray(\n",
+    "    rng.standard_normal(uxgrid.n_edge),\n",
+    "    dims=[\"n_edge\"],\n",
+    "    uxgrid=uxgrid,\n",
+    "    name=\"edge_noise\",\n",
+    ")\n",
+    "\n",
+    "filtered_edge = edge_da.neighborhood_filter(func=np.mean, r=10.0)\n",
+    "print(\"edge input  dims:\", edge_da.dims, \"  shape:\", edge_da.shape)\n",
+    "print(\"edge output dims:\", filtered_edge.dims, \"  shape:\", filtered_edge.shape)"
+   ]
   },
   {
    "cell_type": "markdown",
+   "id": "3ed186c9a28b402fb0bc4494df01f08d",
    "metadata": {},
    "source": "## Multi-Dimensional Data (e.g. Time + Space)\n\nWhen a `UxDataArray` has extra leading dimensions (e.g. `time`), `neighborhood_filter`\napplies the spatial filter independently at each time step and preserves the full\ndimension order.\n"
   },
   {
    "cell_type": "code",
    "execution_count": null,
+   "id": "cb1e1581032b452c9409d6c6813c49d1",
    "metadata": {},
    "outputs": [],
-   "source": "uxds_ts = ux.tutorial.open_dataset(\"outCSne30-timeseries\")\nuxda_ts = uxds_ts[\"psi\"]\n\nprint(\"Input  dims:\", uxda_ts.dims, \" shape:\", uxda_ts.shape)\n\nfiltered_ts = uxda_ts.neighborhood_filter(func=np.mean, r=5.0)\n\nprint(\"Output dims:\", filtered_ts.dims, \" shape:\", filtered_ts.shape)"
+   "source": [
+    "uxds_ts = ux.tutorial.open_dataset(\"outCSne30-timeseries\")\n",
+    "uxda_ts = uxds_ts[\"psi\"]\n",
+    "\n",
+    "print(\"Input  dims:\", uxda_ts.dims, \" shape:\", uxda_ts.shape)\n",
+    "\n",
+    "filtered_ts = uxda_ts.neighborhood_filter(func=np.mean, r=5.0)\n",
+    "\n",
+    "print(\"Output dims:\", filtered_ts.dims, \" shape:\", filtered_ts.shape)"
+   ]
   },
   {
    "cell_type": "markdown",
+   "id": "379cbbc1e968416e875cc15c1202d7eb",
    "metadata": {},
-   "source": "The grid and time dimensions are both preserved.  Because the filter is applied\nper time step, memory usage scales with `n_time \u00d7 n_face` as expected.\n"
+   "source": "The grid and time dimensions are both preserved.  Because the filter is applied\nper time step, memory usage scales with `n_time × n_face` as expected.\n"
   },
   {
    "cell_type": "markdown",
+   "id": "277c27b1587741f2af2001be3712ef0d",
    "metadata": {},
    "source": "## Dataset-Level Usage\n\n`UxDataset.neighborhood_filter` applies the filter to **every data variable**\nthat is mapped to a grid element. Variables without a grid dimension (e.g. scalars\nor time-only arrays) are passed through unchanged.\n"
   },
   {
    "cell_type": "code",
    "execution_count": null,
+   "id": "db7b79bc585a40fcaf58bf750017e135",
    "metadata": {},
    "outputs": [],
-   "source": "uxds_filtered = uxds.neighborhood_filter(func=np.mean, r=5.0)\nuxds_filtered"
+   "source": [
+    "uxds_filtered = uxds.neighborhood_filter(func=np.mean, r=5.0)\n",
+    "uxds_filtered"
+   ]
   },
   {
    "cell_type": "markdown",
+   "id": "916684f9a58a4a2aa5f864670399430d",
    "metadata": {},
    "source": "## Chaining with xarray Operations\n\nBecause `neighborhood_filter` returns a proper `UxDataArray` with its `uxgrid`\npreserved, you can chain it with any standard xarray operation.\n"
   },
   {
    "cell_type": "code",
    "execution_count": null,
+   "id": "1671c31a24314836a5b85d7ef7fbf015",
    "metadata": {},
    "outputs": [],
-   "source": "# Apply the filter and then mask values below zero\nresult = (\n    uxda\n    .neighborhood_filter(func=np.mean, r=5.0)\n    .where(lambda x: x > 0)\n)\nprint(\"Masked result type:\", type(result).__name__)\nprint(\"uxgrid preserved:\", result.uxgrid is not None)\nprint(\"Positive fraction:\", float((result > 0).sum()) / result.size)"
+   "source": [
+    "# Apply the filter and then mask values below zero\n",
+    "result = uxda.neighborhood_filter(func=np.mean, r=5.0).where(lambda x: x > 0)\n",
+    "print(\"Masked result type:\", type(result).__name__)\n",
+    "print(\"uxgrid preserved:\", result.uxgrid is not None)\n",
+    "print(\"Positive fraction:\", float((result > 0).sum()) / result.size)"
+   ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
+   "id": "33b0902fd34d4ace834912fa1002cf8e",
    "metadata": {},
    "outputs": [],
-   "source": "# Group by latitude band after smoothing (standard xarray groupby)\nimport xarray as xr\n\nlat_bins = xr.DataArray(\n    np.digitize(uxda.uxgrid.face_lat.values, bins=np.arange(-90, 91, 30)),\n    dims=[\"n_face\"],\n)\n\nzonal_smooth = (\n    uxda\n    .neighborhood_filter(func=np.mean, r=5.0)\n    .groupby(lat_bins)\n    .mean()\n)\nprint(\"Grouped result type:\", type(zonal_smooth).__name__)\nprint(\"Zonal means:\", zonal_smooth.values)"
+   "source": [
+    "# Group by latitude band after smoothing (standard xarray groupby)\n",
+    "import xarray as xr\n",
+    "\n",
+    "lat_bins = xr.DataArray(\n",
+    "    np.digitize(uxda.uxgrid.face_lat.values, bins=np.arange(-90, 91, 30)),\n",
+    "    dims=[\"n_face\"],\n",
+    ")\n",
+    "\n",
+    "zonal_smooth = uxda.neighborhood_filter(func=np.mean, r=5.0).groupby(lat_bins).mean()\n",
+    "print(\"Grouped result type:\", type(zonal_smooth).__name__)\n",
+    "print(\"Zonal means:\", zonal_smooth.values)"
+   ]
   },
   {
    "cell_type": "markdown",
+   "id": "f6fa52606d8c4a75a9b52967216f8f3f",
    "metadata": {},
    "source": "## Empty Neighborhoods and NaN Behavior\n\nIf the search radius is so small that *no* neighbor is found for a given element,\nthe result for that element is `NaN` rather than an uninitialized garbage value.\nIn practice this only happens for extremely small radii on very coarse grids; the\nfilter always includes the queried element itself, so `r = 0` is safe.\n"
   },
   {
    "cell_type": "code",
    "execution_count": null,
+   "id": "f5a1fa73e5044315a093ec459c9be902",
    "metadata": {},
    "outputs": [],
-   "source": "uxgrid_coarse = ux.Grid.from_healpix(zoom=1)  # 48 faces\nda_coarse = ux.UxDataArray(\n    np.arange(uxgrid_coarse.n_face, dtype=float),\n    dims=[\"n_face\"],\n    uxgrid=uxgrid_coarse,\n)\n\n# r = 0 always catches at least the element itself \u2192 no NaNs\nfiltered_r0 = da_coarse.neighborhood_filter(func=np.mean, r=0.0)\nprint(\"r = 0:   NaN count =\", int(np.isnan(filtered_r0.values).sum()))\n\n# r = 360 catches every element \u2192 all values equal the global mean\nfiltered_global = da_coarse.neighborhood_filter(func=np.mean, r=360.0)\nprint(\"r = 360: all equal global mean?\",\n      np.allclose(filtered_global.values, da_coarse.values.mean()))"
+   "source": [
+    "uxgrid_coarse = ux.Grid.from_healpix(zoom=1)  # 48 faces\n",
+    "da_coarse = ux.UxDataArray(\n",
+    "    np.arange(uxgrid_coarse.n_face, dtype=float),\n",
+    "    dims=[\"n_face\"],\n",
+    "    uxgrid=uxgrid_coarse,\n",
+    ")\n",
+    "\n",
+    "# r = 0 always catches at least the element itself → no NaNs\n",
+    "filtered_r0 = da_coarse.neighborhood_filter(func=np.mean, r=0.0)\n",
+    "print(\"r = 0:   NaN count =\", int(np.isnan(filtered_r0.values).sum()))\n",
+    "\n",
+    "# r = 360 catches every element → all values equal the global mean\n",
+    "filtered_global = da_coarse.neighborhood_filter(func=np.mean, r=360.0)\n",
+    "print(\n",
+    "    \"r = 360: all equal global mean?\",\n",
+    "    np.allclose(filtered_global.values, da_coarse.values.mean()),\n",
+    ")"
+   ]
   },
   {
    "cell_type": "markdown",
+   "id": "cdf66aed5cc84ca1b48e60bad68798a8",
    "metadata": {},
-   "source": "## API Reference\n\nSee also:\n\n- {py:meth}`uxarray.UxDataArray.neighborhood_filter`\n- {py:meth}`uxarray.UxDataset.neighborhood_filter`\n\nRelated methods that apply aggregations across different grid element types:\n\n- {py:meth}`uxarray.UxDataArray.topological_mean` \u2014 aggregate node\u2192face, node\u2192edge, etc.\n- {py:meth}`uxarray.UxDataArray.zonal_mean` \u2014 latitude-band averages\n- {py:meth}`uxarray.UxDataArray.azimuthal_mean` \u2014 rings of constant great-circle distance\n"
+   "source": "## API Reference\n\nSee also:\n\n- {py:meth}`uxarray.UxDataArray.neighborhood_filter`\n- {py:meth}`uxarray.UxDataset.neighborhood_filter`\n\nRelated methods that apply aggregations across different grid element types:\n\n- {py:meth}`uxarray.UxDataArray.topological_mean` — aggregate node→face, node→edge, etc.\n- {py:meth}`uxarray.UxDataArray.zonal_mean` — latitude-band averages\n- {py:meth}`uxarray.UxDataArray.azimuthal_mean` — rings of constant great-circle distance\n"
   }
  ],
  "metadata": {
@@ -201,4 +374,4 @@
  },
  "nbformat": 4,
  "nbformat_minor": 5
-}
\ No newline at end of file
+}
diff --git a/test/core/test_dataarray.py b/test/core/test_dataarray.py
index 888ff93b1..7ebc9f01d 100644
--- a/test/core/test_dataarray.py
+++ b/test/core/test_dataarray.py
@@ -306,4 +306,3 @@ def test_auto_transpose_direct_on_uxdataarray(self, gridpath, datasetpath):
         # Dim order must be restored to (n_face, time)
         assert filtered2.dims == ("n_face", "time")
         assert filtered2.shape == (uxda.shape[0], 2)
-