Implements Grid.compute_skewness(..., method="equiangle")

docs/api.rst

@@ -195,6 +195,7 @@ Methods
    Grid.calculate_total_face_area
    Grid.compute_face_areas
    Grid.compute_face_node_angles
+   Grid.compute_skewness
    Grid.construct_face_centers
    Grid.get_ball_tree
    Grid.get_kd_tree

test/grid/geometry/test_angles.py

@@ -65,3 +65,47 @@ def test_face_node_angles_hexagons_and_pentagons():
     spherical_excess = angles.sum('n_max_face_nodes') - (grid.n_nodes_per_face - 2) * np.pi
     face_areas = grid.compute_face_areas()
     assert np.allclose(spherical_excess, face_areas, atol=0, rtol=1e-9)
+
+def test_equiangle_skewness():
+    """just some tests about Grid.compute_skewness(method="equiangle")..."""
+
+    ## spot check with some triangles
+    # make a tiny 90,60,30 degrees triangle:
+    # (n1)
+    #  |   %%
+    #  |       %%
+    # (n0) ------ (n2)
+    node_lon = [0, 0, np.sqrt(3)]
+    node_lat = [0, 1, 0]
+    face_node_connectivity = [[0, 1, 2]]
+    grid = ux.Grid.from_topology(node_lon, node_lat, face_node_connectivity)
+    # expect skewness = max((Amax - Areg) / (180 degrees - Areg), (Areg - Amin) / Areg)
+    # here, Areg ~= 60 degrees (it's a tiny triangle), Amax = 90 degrees, Amin = 30 degrees
+    #   --> max((90 - 60) / (180 - 60), (60 - 30) / 60) --> max(0.25, 0.5)
+    skewness = grid.compute_skewness(method="equiangle")
+    assert np.isclose(skewness, 0.5, atol=1e-4, rtol=0)
+
+    ## much larger "would-be 90,60,30" triangle, to ensure accounting for spherical geometry.
+    # It's actually 90, 68.6, 36.2 degrees, due to spherical geometry.
+    node_lon = [0, 0, 30*np.sqrt(3)]
+    node_lat = [0, 30, 0]
+    face_node_connectivity = [[0, 1, 2]]
+    grid = ux.Grid.from_topology(node_lon, node_lat, face_node_connectivity)
+    skewness = grid.compute_skewness(method="equiangle")
+    assert not np.isclose(skewness, 0.5, atol=1e-2, rtol=0)
+    assert 0.42 < skewness < 0.46  # hard-coding "expected" answer.
+    # If using correct angles but not spherical geometry in Areg, would assume Areg = 60
+    #   --> skewness = max((90 - 68.6) / (180 - 68.6), (68.6 - 36.2) / 68.6) --> max(0.192, 0.472)
+    # Therefore, the test above is indeed sensitive to
+    #   "does the skewness formula actually use the proper Areg based on spherical geometry?"
+
+    ## spot check of grid with "not all faces have same number of nodes"
+    grid = ux.tutorial.open_grid('mpas-QU-480')
+    skewness = grid.compute_skewness(method="equiangle")
+    assert isinstance(skewness, xr.DataArray)
+    # there happens to be at least one very non-skew shape in this grid
+    assert 0 < skewness.min() < 1e-8
+    # there aren't any very skew shapes in this grid
+    assert skewness.max() < 0.15
+    # skewness is well defined for all faces
+    assert not np.any(skewness.isnull())

uxarray/grid/angles.py

@@ -16,8 +16,7 @@ def _compute_face_node_angles_convex(
     face_node_connectivity,
     n_nodes_per_face,
 ):
-    """
-    Calculate the angles at each node for each face, assuming convex faces
+    """Returns angles [in radians] at each node for each face, assuming convex faces
     and a spherical geometry (these assumptions occur throughout uxarray).
 
     Parameters
@@ -75,3 +74,42 @@ def _compute_face_node_angles_convex(
             )
             result[i, j] = _small_angle_of_2_vectors(v1, v2)
     return result
+
+
+def _compute_equiangle_skewness(face_node_angles, n_nodes_per_face):
+    """Returns the equiangle skewness at each face:
+        max((Amax - Areg) / (pi - Areg), (Areg - Amin) / Areg)
+    where
+        Amin, Amax = min, max of the angles at the nodes of the face
+        Areg = internal angle at all nodes for a regular polygon with
+            the same number of sides and covering the same area as this face.
+
+    In a flat geometry, the sum of angles in a polygon with n sides is (n-2)*pi.
+    Splitting the angles equally to form a regular polygon yields Areg_flat = (n-2)*pi/n.
+    However, for a spherical geometry, the sum of angles depends on face area:
+        sum(angles) = (n-2)*pi + face_area / sphere_radius^2
+    Areg should be based on a regular polygon with same area as the corresponding face,
+    so, splitting the angles equally to form a regular polygon yields simply:
+        Areg = sum(angles)/n.
+
+    Parameters
+    ----------
+    face_node_angles : xr.DataArray or UxDataArray with dims 'n_face', 'n_max_face_nodes'
+        Angles [in radians] at each node of each face.
+    n_nodes_per_face : xr.DataArray or UxDataArray with dims 'n_face'
+        Number of nodes for each face.
+
+    Returns
+    -------
+    xr.DataArray or UxDataArray with dims 'n_face'
+        Equiangle skewness for each face.
+        Type matches the input type (xr.DataArray or UxDataArray).
+    """
+    Amin = face_node_angles.min("n_max_face_nodes", skipna=True)
+    Amax = face_node_angles.max("n_max_face_nodes", skipna=True)
+    Areg = face_node_angles.sum("n_max_face_nodes", skipna=True) / n_nodes_per_face
+    term0 = (Amax - Areg) / (np.pi - Areg)
+    term1 = (Areg - Amin) / Areg
+    # Should just use np.maximum(term0, term1), but that drops UxDataArray type currently,
+    # so use where as a workaround for now. TODO: swap to np.maximum after fixing issue #1685.
+    return term0.where(term0 > term1, term1)

uxarray/grid/grid.py

@@ -19,7 +19,10 @@
 from uxarray.cross_sections import GridCrossSectionAccessor
 from uxarray.errors import DataCenteringError, DimensionError, GridInvalidError
 from uxarray.formatting_html import grid_repr
-from uxarray.grid.angles import _compute_face_node_angles_convex
+from uxarray.grid.angles import (
+    _compute_equiangle_skewness,
+    _compute_face_node_angles_convex,
+)
 from uxarray.grid.area import _get_all_face_area_from_coords
 from uxarray.grid.bounds import _populate_face_bounds
 from uxarray.grid.connectivity import (
@@ -1973,6 +1976,40 @@ def copy(self):
             source_dims_dict=self._source_dims_dict,
         )
 
+    def compute_skewness(self, method: str = "equiangle", *, as_uxarray: bool = False):
+        """Returns the skewness of each face in the grid, computed using the specified method.
+        Skewness is a measure of how much a face deviates from being regular,
+        e.g. having equal angles at all nodes. Values close to 0 indicate a regular face,
+        while values close to 1 indicate a highly skewed / nearly degenerate face.
+
+        Parameters
+        ----------
+        method: str, defaults to "equiangle"
+            The method to use for computing skewness. Options are:
+            - "equiangle": computes the equiangular skewness of each face:
+                equiangle_skewness = max((Amax - Areg) / (pi - Areg), (Areg - Amin) / Areg)
+                where Amin, Amax = min, max of the angles at the nodes of the face,
+                and Areg = internal angle at all nodes for a regular polygon with
+                the same number of sides and covering the same area as this face.
+            - (other options not yet implemented)
+        as_uxarray: bool, defaults to False
+            Whether to return a uxarray.DataArray (if True) or an xarray.DataArray (if False).
+            If True, equivalent to uxarray.DataArray(self.compute_skewness(..., as_uxarray=False), uxgrid=self).
+
+        Returns
+        -------
+        skewness : xr.DataArray or uxarray.UxDataArray (if as_uxarray=True)
+            The skewness of each face in the grid.
+            Has 'n_face' dimension, with same size as in self.
+        """
+        if method == "equiangle":
+            face_node_angles = self.compute_face_node_angles(as_uxarray=as_uxarray)
+            return _compute_equiangle_skewness(face_node_angles, self.n_nodes_per_face)
+        else:
+            raise NotImplementedError(
+                f"Skewness computation method '{method}' is not implemented."
+            )
+
     def compute_face_node_angles(
         self,
         *,