Fix HEALPix face area calculation to enforce equal area property (#1379)

* Fix HEALPix face area calculation to enforce equal area property

* o Add area logic to io folder for HEALPix mesh, document notebook and area cal

* o Fix notebook and cleanup

* Add note

* Addressing PR comments

* Update uxarray/grid/grid.py

Co-authored-by: Orhan Eroglu <[email protected]>

* o More cleanup; The compute_face_areas() method is still present in the code but marked as deprecated

* o Add note

* Update area and healpix notebooks

* o Cleanup notebook

* o Add area math

* Small rephrasing on two notebooks

---------

Co-authored-by: Orhan Eroglu <[email protected]>
Co-authored-by: erogluorhan <[email protected]>

Author: rajeeja

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

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

@@ -358,7 +358,7 @@ def test_face_areas_calculate_total_face_area_sphere(self, gridpath, mesh_consta
 def test_face_areas_compute_face_areas_geoflow_small(gridpath):
     """Checks if the GeoFlow Small can generate a face areas output."""
     grid_geoflow = ux.open_grid(gridpath("ugrid", "geoflow-small", "grid.nc"))
-    grid_geoflow.compute_face_areas()
+    grid_geoflow._compute_face_areas()
 
 class TestFaceAreas:
     def test_face_areas_verts_calc_area(self, gridpath, mesh_constants):

test/grid/test_grid.py

@@ -95,6 +95,106 @@ def test_invalid_cells():
     with pytest.raises(ValueError):
         uxda = ux.UxDataset.from_healpix(xrda)
 
+@pytest.mark.parametrize("resolution_level", [0, 2])  # Test lowest and a mid-level
+def test_healpix_equal_area_property(resolution_level):
+    """Test that all HEALPix faces have equal area as per the HEALPix specification.
+
+    HEALPix (Hierarchical Equal Area isoLatitude Pixelization) is designed so that
+    all pixels/faces have exactly the same spherical area.
+    """
+    # Create HEALPix grid with boundaries for area calculation
+    uxgrid = ux.Grid.from_healpix(resolution_level, pixels_only=False)
+
+    # Calculate face areas and expected theoretical area
+    face_areas = uxgrid.face_areas.values
+    nside = hp.order2nside(resolution_level)
+    npix = hp.nside2npix(nside)
+    expected_area_per_pixel = 4 * np.pi / npix
+
+    # All face areas should be equal to the theoretical value
+    # Use a more reasonable tolerance for numerical computations
+    np.testing.assert_allclose(
+        face_areas,
+        expected_area_per_pixel,
+        rtol=1e-12,  # Relaxed from 1e-14 for numerical stability
+        atol=1e-15,  # Added absolute tolerance
+        err_msg=f"HEALPix faces do not have equal areas at resolution {resolution_level}"
+    )
+
+    # Verify total surface area equals 4π
+    total_area = np.sum(face_areas)
+    expected_total_area = 4 * np.pi
+    np.testing.assert_allclose(
+        total_area,
+        expected_total_area,
+        rtol=1e-12,  # Relaxed from 1e-14
+        atol=1e-15,  # Added absolute tolerance
+        err_msg=f"Total HEALPix surface area incorrect at resolution {resolution_level}"
+    )
+
+
+def test_healpix_face_areas_consistency():
+    """Test that HEALPix face areas are consistent across different resolution levels."""
+    resolution_levels = [0, 1]  # Just test basic functionality
+
+    for resolution_level in resolution_levels:
+        uxgrid = ux.Grid.from_healpix(resolution_level, pixels_only=False)
+        face_areas = uxgrid.face_areas.values
+
+        # All faces should have identical areas
+        area_std = np.std(face_areas)
+        area_mean = np.mean(face_areas)
+
+        # Avoid division by zero
+        if area_mean == 0:
+            pytest.fail(f"Mean face area is zero at resolution {resolution_level}")
+
+        relative_std = area_std / area_mean
+
+        # Relative standard deviation should be essentially zero (numerical precision)
+        # Relaxed tolerance for numerical stability
+        assert relative_std < 1e-12, (
+            f"Face areas not equal at resolution {resolution_level}: "
+            f"relative_std={relative_std:.2e}"
+        )
+
+        # Check that face areas match theoretical calculation
+        nside = hp.order2nside(resolution_level)
+        npix = hp.nside2npix(nside)
+        theoretical_area = 4 * np.pi / npix
+
+        # Use consistent tolerance with the parametrized test
+        np.testing.assert_allclose(
+            face_areas,
+            theoretical_area,
+            rtol=1e-12,
+            atol=1e-15,
+            err_msg=f"Face areas incorrect at resolution {resolution_level}"
+        )
+
+
+@pytest.mark.parametrize("resolution_level", [0, 1])  # Just test the scaling relationship
+def test_healpix_area_scaling(resolution_level):
+    """Test that face areas scale correctly with resolution level."""
+    # Create grids at consecutive resolution levels
+    uxgrid_current = ux.Grid.from_healpix(resolution_level, pixels_only=False)
+    uxgrid_next = ux.Grid.from_healpix(resolution_level + 1, pixels_only=False)
+
+    area_current = uxgrid_current.face_areas.values[0]  # All areas are equal
+    area_next = uxgrid_next.face_areas.values[0]
+
+    # Each resolution level increases npix by factor of 4, so area decreases by factor of 4
+    expected_ratio = 4.0
+    actual_ratio = area_current / area_next
+
+    np.testing.assert_allclose(
+        actual_ratio,
+        expected_ratio,
+        rtol=1e-12,
+        err_msg=f"Area scaling incorrect between resolution {resolution_level} and {resolution_level + 1}"
+    )
+
+
 def test_healpix_round_trip_consistency(tmp_path):
     """Test round-trip serialization of HEALPix grid through UGRID and Exodus formats.
 

test/io/test_healpix.py

@@ -60,7 +60,7 @@ def test_calculate_face_area(mesh_constants):
         fill_value=-1,
     )
 
-    area, _ = grid.compute_face_areas()
+    area, _ = grid._compute_face_areas()
     nt.assert_almost_equal(area, mesh_constants['TRI_AREA'], decimal=5)
 
 def test_quadrature():

test/test_helpers.py

@@ -484,9 +484,7 @@ def integrate(
         >>> integral = uxds["psi"].integrate()
         """
         if self.values.shape[-1] == self.uxgrid.n_face:
-            face_areas, face_jacobian = self.uxgrid.compute_face_areas(
-                quadrature_rule, order
-            )
+            face_areas = self.uxgrid.face_areas.values
 
             # perform dot product between face areas and last dimension of data
             integral = np.einsum("i,...i", face_areas, self.values)

uxarray/core/dataarray.py

@@ -586,10 +586,7 @@ def integrate(self, quadrature_rule="triangular", order=4):
 
         integral = 0.0
 
-        # call function to get area of all the faces as a np array
-        face_areas, face_jacobian = self.uxgrid.compute_face_areas(
-            quadrature_rule, order
-        )
+        face_areas = self.uxgrid.face_areas.values
 
         # TODO: Should we fix this requirement? Shouldn't it be applicable to
         # TODO: all variables of dataset or a dataarray instead?