Merge branch 'main' into rajeeja/zonal-mean-analytic-band

Author: rajeeja

.github/workflows/ci.yml

@@ -19,7 +19,7 @@ jobs:
       fail-fast: false
       matrix:
         os: [ "ubuntu-latest", "macos-latest", "windows-latest"]
-        python-version: [ "3.10", "3.11", "3.12", "3.13"]
+        python-version: [ "3.11", "3.12", "3.13"]
     steps:
       - name: Cancel previous runs
         uses: styfle/[email protected]

.pre-commit-config.yaml

@@ -14,7 +14,7 @@ repos:
 
 - repo: https://github.com/astral-sh/ruff-pre-commit
   # Ruff version.
-  rev: v0.13.0
+  rev: v0.13.1
   hooks:
     - id: ruff
       name: lint with ruff

docs/api.rst

@@ -148,7 +148,6 @@ Methods
    Grid.copy
    Grid.chunk
    Grid.validate
-   Grid.compute_face_areas
    Grid.calculate_total_face_area
    Grid.normalize_cartesian_coordinates
    Grid.construct_face_centers
@@ -453,6 +452,17 @@ on each face.
    UxDataArray.topological_all
    UxDataArray.topological_any
 
+Azimuthal
+~~~~~~~~~
+
+Azimuthal aggregations apply an aggregation (i.e. averaging) along circles of constant great-circle distance from a specified point on the sphere.
+
+
+.. autosummary::
+   :toctree: generated/
+
+   UxDataArray.azimuthal_mean
+
 Zonal Average
 ~~~~~~~~~~~~~
 .. autosummary::

docs/user-guide/area_calc.ipynb

@@ -109,7 +109,9 @@
    "cell_type": "markdown",
    "id": "b464bd06cf9dbf3b",
    "metadata": {},
-   "source": "However, even if you load in a HEALPix grid specifying that you do not want the connectivity upfront, they can still be constructed when desired because of UXarray's internal design."
+   "source": [
+    "However, even if you load in a HEALPix grid specifying that you do not want the connectivity upfront, they can still be constructed when desired because of UXarray's internal design."
+   ]
   },
   {
    "cell_type": "code",
@@ -191,7 +193,9 @@
    "cell_type": "markdown",
    "id": "e1ef5671b3006b0c",
    "metadata": {},
-   "source": "Before remapping, we can plot our Source and Destination grids."
+   "source": [
+    "Before remapping, we can plot our Source and Destination grids."
+   ]
   },
   {
    "cell_type": "code",
@@ -209,7 +213,9 @@
    "cell_type": "markdown",
    "id": "d9f677166dd1b9f",
    "metadata": {},
-   "source": "We can now perform our remapping. In this case, we apply a simple nearest neighbor remapping."
+   "source": [
+    "We can now perform our remapping. In this case, we apply a simple nearest neighbor remapping."
+   ]
   },
   {
    "cell_type": "code",
@@ -226,7 +232,9 @@
    "cell_type": "markdown",
    "id": "f829232283eeb4",
    "metadata": {},
-   "source": "Our original data variable now resides on our HEALPix grid."
+   "source": [
+    "Our original data variable now resides on our HEALPix grid."
+   ]
   },
   {
    "cell_type": "code",
@@ -259,6 +267,94 @@
     "psi_hp.to_dataset().to_xarray(grid_format=\"HEALPix\").to_netcdf(\"psi_healpix.nc\")"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "healpix_equal_area_section",
+   "metadata": {},
+   "source": [
+    "## HEALPix Equal Area Property\n",
+    "\n",
+    "HEALPix's fundamental property is that **all pixels at a given resolution have exactly the same spherical area**. This \"equal area\" property is crucial for scientific applications such as global averaging, zonal means, and conservation properties in regridding operations.\n",
+    "\n",
+    "The theoretical area of each HEALPix pixel is given by:\n",
+    "\n",
+    "$$\n",
+    "A_{\\text{pixel}} = \\frac{4\\pi}{N_{\\text{pix}}} = \\frac{4\\pi}{12 \\cdot 4^{\\text{resolution}}}\n",
+    "$$\n",
+    "\n",
+    "where  \n",
+    "$N_{\\text{pix}} = 12 \\cdot 4^{\\text{resolution}}$ is the total number of pixels at a given resolution level, and  \n",
+    "$4\\pi$ is the total surface area of the unit sphere in steradians. \n",
+    "\n",
+    "### Geometric Representation Considerations\n",
+    "\n",
+    "UXarray represents HEALPix grids using Great Circle Arcs (GCAs) to define pixel boundaries, following UGRID conventions. However, this geometric representation can introduce systematic errors when computing areas numerically, potentially violating HEALPix's equal-area property. \n",
+    "\n",
+    "```{note}\n",
+    "To alleviate the impacts of this systematic differences between UXarray and HEALPix, we adjust our `Grid.face_areas` property to fulfill the HEALPix equal area property, making sure that all the faces in a HEALPix mesh have the same theoretical HEALPix area.\n",
+    "```\n",
+    "\n",
+    "Let's demonstrate this with HEALPix's 12 base pixels:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "healpix_area_demo",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "\n",
+    "# Create HEALPix grid at resolution 0 (12 base pixels)\n",
+    "grid = ux.Grid.from_healpix(0, pixels_only=False)\n",
+    "\n",
+    "# Using face_areas property ensures equal areas for HEALPix\n",
+    "hp_areas = grid.face_areas.values\n",
+    "print(f\"Standard deviation: {np.std(hp_areas):.2e}\")\n",
+    "print(f\"All pixels have area: {hp_areas[0]:.6f} steradians\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "6e8b7185-c3ae-4f50-939b-4cb28a952e6f",
+   "metadata": {},
+   "source": [
+    "### Still need geometric face area calculations for a HEALPix mesh?\n",
+    "\n",
+    "For most use cases, the `Grid.face_areas` property provides the recommended approach for accessing face areas. However, if you specifically need geometric calculations of individual HEALPix faces as they are represented in UXarray (rather than the theoretical equal areas), you may want to access the internal computation method, `Grid._compute_face_areas()`. Note that this approach may not preserve HEALPix's equal-area property due to geometric representation differences. \n",
+    "\n",
+    "Look at the following example:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "dfd494b5-0b3e-472f-8896-bd3566747bdb",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# For advanced use cases: access geometric calculations (may not preserve equal-area property)\n",
+    "hp_geometric_areas, face_jacobians = grid._compute_face_areas(\n",
+    "    quadrature_rule=\"triangular\", order=4\n",
+    ")\n",
+    "print(\n",
+    "    f\"Geometric std deviation: {hp_geometric_areas.std():.2e} (vs theoretical equal areas)\"\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "healpix_error_analysis",
+   "metadata": {},
+   "source": [
+    "If you are interested in futher details of the systematic errors due to geometric calculation with a clear spatial pattern, our analysis across resolution levels 0-7 shows:\n",
+    "\n",
+    "- **Equatorial pixels** (lat≈0°): +20.5% area error\n",
+    "- **Mid-latitude pixels** (lat≈±42°): -9.9% area error  \n",
+    "- **Maximum errors**: Persist at ~10% even at fine resolutions\n"
+   ]
+  },
   {
    "cell_type": "markdown",
    "id": "7988c071af403bf2",
@@ -293,7 +389,9 @@
    "cell_type": "markdown",
    "id": "11ff4e1ad036ea53",
    "metadata": {},
-   "source": "The interface above looks almost identical to what you would see if you loaded in the file directly with Xarray."
+   "source": [
+    "The interface above looks almost identical to what you would see if you loaded in the file directly with Xarray."
+   ]
   },
   {
    "cell_type": "code",
@@ -359,7 +457,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.12.3"
+   "version": "3.12.6"
   }
  },
  "nbformat": 4,

docs/user-guide/healpix.ipynb

@@ -138,10 +138,10 @@ def test_resample_reduces_time_dimension(gridpath):
     uxds = ux.UxDataset(ds, uxgrid=ux.open_grid(gridpath("ugrid", "outCSne30", "outCSne30.ug")))
 
     # Test monthly resampling reduces from 365 days to 12 months
-    monthly = uxds.resample(time="1M").mean()
+    monthly = uxds.resample(time="1ME").mean()
     assert "time" in monthly.dims, "time dimension missing after resample"
-    assert monthly.dims["time"] < uxds.dims["time"], "time dimension not reduced"
-    assert monthly.dims["time"] <= 12, "monthly resampling should give 12 or fewer time points"
+    assert monthly.sizes["time"] < uxds.sizes["time"], "time dimension not reduced"
+    assert monthly.sizes["time"] <= 12, "monthly resampling should give 12 or fewer time points"
 
 
 def test_resample_with_cftime(gridpath):
@@ -166,10 +166,10 @@ def test_resample_with_cftime(gridpath):
     uxds = ux.UxDataset(ds, uxgrid=ux.open_grid(gridpath("ugrid", "outCSne30", "outCSne30.ug")))
 
     # Test that quarterly resampling works with cftime
-    quarterly = uxds.resample(time="Q").mean()
+    quarterly = uxds.resample(time="QE").mean()
     assert hasattr(quarterly, "uxgrid"), "uxgrid not preserved with cftime resampling"
     assert "time" in quarterly.dims, "time dimension missing after cftime resample"
-    assert quarterly.dims["time"] < uxds.dims["time"], "time dimension not reduced with cftime"
+    assert quarterly.sizes["time"] < uxds.sizes["time"], "time dimension not reduced with cftime"
 
 
 def test_rolling_preserves_uxgrid(gridpath):
@@ -237,7 +237,7 @@ def test_coarsen_preserves_uxgrid(gridpath):
 
     # Test that coarsen reduces dimension correctly
     assert len(da_result.time) == 8, "coarsen by 3 should reduce 24 points to 8"
-    assert ds_result.dims["time"] == 8, "coarsen should reduce time dimension"
+    assert ds_result.sizes["time"] == 8, "coarsen should reduce time dimension"
 
 
 def test_weighted_preserves_uxgrid(gridpath, datasetpath):
@@ -251,7 +251,7 @@ def test_weighted_preserves_uxgrid(gridpath, datasetpath):
         gridpath("ugrid", "outCSne30", "outCSne30.ug"),
         datasetpath("ugrid", "outCSne30", "outCSne30_var2.nc")
     )
-    n_face = uxds_base.dims["n_face"]
+    n_face = uxds_base.sizes["n_face"]
 
     # Create data with time and face dimensions
     temp_data = np.random.rand(10, n_face)

test/core/test_accessors.py

@@ -0,0 +1,68 @@
+import pytest
+import uxarray as ux
+import numpy as np
+
+
+
+def test_gaussian(gridpath, datasetpath):
+    uxds = ux.open_dataset(
+        gridpath("mpas", "dyamond-30km", "gradient_grid_subset.nc"),
+        datasetpath("mpas", "dyamond-30km", "gradient_data_subset.nc")
+    )
+
+    res = uxds['gaussian'].azimuthal_mean(center_coord=(45, 0), outer_radius=2, radius_step=0.5)
+
+    # Expects decreasing values from center
+    valid_vals = res[1:]
+
+    np.testing.assert_array_less(
+        valid_vals.diff("radius").values, 1e-12
+    )
+
+
+
+def test_inverse_gaussian(gridpath, datasetpath):
+
+    uxds =  ux.open_dataset(
+        gridpath("mpas", "dyamond-30km", "gradient_grid_subset.nc"),
+        datasetpath("mpas", "dyamond-30km", "gradient_data_subset.nc")
+    )
+
+    res = uxds['inverse_gaussian'].azimuthal_mean(center_coord=(45, 0), outer_radius=2, radius_step=0.5)
+
+    # Expects increasing values from center
+    atol = 1e-12
+    diffs = res[1:].diff("radius").values
+    diffs = diffs[np.isfinite(diffs)]
+    np.testing.assert_array_less(-atol, diffs)
+
+def test_non_zero_hit_counts(gridpath, datasetpath):
+    uxds =  ux.open_dataset(
+        gridpath("mpas", "dyamond-30km", "gradient_grid_subset.nc"),
+        datasetpath("mpas", "dyamond-30km", "gradient_data_subset.nc")
+    )
+
+    res, hit_counts = uxds['inverse_gaussian'].azimuthal_mean(center_coord=(45, 0), outer_radius=2, radius_step=0.5, return_hit_counts=True)
+
+    # At least one hit after the first circle
+    assert np.all(hit_counts[1:] > 1)
+
+    assert 'radius' in hit_counts.dims
+    assert hit_counts.sizes['radius'] == res.sizes['radius']
+
+def test_zero_hit_counts(gridpath, datasetpath):
+    uxds =  ux.open_dataset(
+        gridpath("mpas", "dyamond-30km", "gradient_grid_subset.nc"),
+        datasetpath("mpas", "dyamond-30km", "gradient_data_subset.nc")
+    )
+
+    # Outside of grid domain
+    res, hit_counts = uxds['inverse_gaussian'].azimuthal_mean(center_coord=(-45, 0), outer_radius=2, radius_step=0.5, return_hit_counts=True)
+
+    assert 'radius' in hit_counts.dims
+    assert hit_counts.sizes['radius'] == res.sizes['radius']
+
+    # No hits
+    assert np.all(hit_counts == 0)
+
+    print(res)

test/core/test_azimuthal.py

@@ -16,7 +16,7 @@ def ds():
 
     uxds["fc"] = uxds["fc"].assign_coords(face_id=("n_face", np.arange(uxgrid.n_face)))
     uxds["nc"] = uxds["nc"].assign_coords(node_id=("n_node", np.arange(uxgrid.n_node)))
-    uxds["t"] = uxds["t"].assign_coords(time_id=("time", np.arange(uxds.dims["time"])))
+    uxds["t"] = uxds["t"].assign_coords(time_id=("time", np.arange(uxds.sizes["time"])))
 
     return uxds
 
@@ -109,7 +109,7 @@ def test_groupby(self, ds):
         assert hasattr(grouped, 'uxgrid')
 
     def test_assign_coords(self, ds):
-        n = ds.dims['n_face']
+        n = ds.sizes['n_face']
         new_coord = xr.DataArray(np.arange(n) * 10, dims=['n_face'])
         out = ds.assign_coords(scaled_id=new_coord)
         assert isinstance(out, ux.UxDataset)
@@ -121,7 +121,7 @@ def test_expand_dims(self, ds):
         out = ds.expand_dims({'member': 1})
         assert isinstance(out, ux.UxDataset)
         assert hasattr(out, 'uxgrid') and out.uxgrid is ds.uxgrid
-        assert 'member' in out.dims and out.dims['member'] == 1
+        assert 'member' in out.dims and out.sizes['member'] == 1
         assert 'member' in out['t'].dims and out['t'].sizes['member'] == 1
 
     def test_method_chaining(self, ds):
@@ -147,7 +147,7 @@ def test_stack_unstack(self, ds):
         assert unstacked['fc'].shape == ds_fc['fc'].shape
 
     def test_sortby(self, ds):
-        n = ds.dims['n_face']
+        n = ds.sizes['n_face']
         ds_fc = ds[['fc']].assign_coords(reverse_id=('n_face', np.arange(n)[::-1]))
         out = ds_fc.sortby('reverse_id')
         assert isinstance(out, ux.UxDataset)

test/core/test_inheritance.py

@@ -0,0 +1 @@
+"""Integration test module."""