From c4fd80f0882b1a0ebc9c7d6ed705ffbd3683509d Mon Sep 17 00:00:00 2001 From: Sam Evans <47793072+Sevans711@users.noreply.github.com> Date: Wed, 19 Aug 2026 14:53:32 -0400 Subject: [PATCH 1/4] add Grid.compute_skewness(method="equiangle") --- uxarray/grid/angles.py | 43 ++++++++++++++++++++++++++++++++++++++++-- uxarray/grid/grid.py | 36 ++++++++++++++++++++++++++++++++++- 2 files changed, 76 insertions(+), 3 deletions(-) diff --git a/uxarray/grid/angles.py b/uxarray/grid/angles.py index 678f94cb4..b573409cf 100644 --- a/uxarray/grid/angles.py +++ b/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,43 @@ def _compute_face_node_angles_convex( ) result[i, j] = _small_angle_of_2_vectors(v1, v2) return result + +def _regular_polygon_angle(n, *, degrees=False): + """returns the internal angle for a regular polygon with n sides. + Equivalent to (n-2)*180/n degrees or (n-2)*pi/n radians. + + n can be a scalar or an array of integers (e.g., can use n = n_nodes_per_face) + """ + if degrees: + return (n - 2) * 180 / n + else: + return (n - 2) * np.pi / n + +def _compute_equiangle_skewness(face_node_angles, n_nodes_per_face): + """Returns the equiangle skewness at each face: + max((Amax - Areg) / (pi - Areg), (Amin - Areg) / Areg) + where + Amin, Amax = min, max of the angles at the nodes of the face + Areg = ideal angle for a regular polygon with n_nodes_per_face sides. + + 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 = _regular_polygon_angle(n_nodes_per_face, degrees=False) + term0 = (Amax - Areg) / (np.pi - Areg) + term1 = (Amin - Areg) / 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) diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py index 711b30389..9be9c7a9d 100644 --- a/uxarray/grid/grid.py +++ b/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,37 @@ 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. + + Parameters + ---------- + method: str, defaults to "equiangle" + The method to use for computing skewness. Options are: + - "equiangle": computes the equiangular skewness of each 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, *, From 2127c39afbb9a75c43d6c1e3f657997b22ad6a5f Mon Sep 17 00:00:00 2001 From: Sam Evans <47793072+Sevans711@users.noreply.github.com> Date: Wed, 19 Aug 2026 14:53:53 -0400 Subject: [PATCH 2/4] ruff formatting --- uxarray/grid/angles.py | 6 ++++-- uxarray/grid/grid.py | 6 +----- 2 files changed, 5 insertions(+), 7 deletions(-) diff --git a/uxarray/grid/angles.py b/uxarray/grid/angles.py index b573409cf..94441cd54 100644 --- a/uxarray/grid/angles.py +++ b/uxarray/grid/angles.py @@ -75,6 +75,7 @@ def _compute_face_node_angles_convex( result[i, j] = _small_angle_of_2_vectors(v1, v2) return result + def _regular_polygon_angle(n, *, degrees=False): """returns the internal angle for a regular polygon with n sides. Equivalent to (n-2)*180/n degrees or (n-2)*pi/n radians. @@ -86,6 +87,7 @@ def _regular_polygon_angle(n, *, degrees=False): else: return (n - 2) * np.pi / n + def _compute_equiangle_skewness(face_node_angles, n_nodes_per_face): """Returns the equiangle skewness at each face: max((Amax - Areg) / (pi - Areg), (Amin - Areg) / Areg) @@ -106,8 +108,8 @@ def _compute_equiangle_skewness(face_node_angles, n_nodes_per_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) + Amin = face_node_angles.min("n_max_face_nodes", skipna=True) + Amax = face_node_angles.max("n_max_face_nodes", skipna=True) Areg = _regular_polygon_angle(n_nodes_per_face, degrees=False) term0 = (Amax - Areg) / (np.pi - Areg) term1 = (Amin - Areg) / Areg diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py index 9be9c7a9d..9547f897e 100644 --- a/uxarray/grid/grid.py +++ b/uxarray/grid/grid.py @@ -1976,11 +1976,7 @@ def copy(self): source_dims_dict=self._source_dims_dict, ) - def compute_skewness(self, - *, - method: str = "equiangle", - as_uxarray: bool = False - ): + 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. Parameters From fc897e4334021ca56509b89129fd45cd94d67e46 Mon Sep 17 00:00:00 2001 From: Sam Evans <47793072+Sevans711@users.noreply.github.com> Date: Wed, 19 Aug 2026 15:42:34 -0400 Subject: [PATCH 3/4] fix: equiangle skewness spherical excess previous implementation wasn't accounting for spherical excess! --- uxarray/grid/angles.py | 25 +++++++++++-------------- uxarray/grid/grid.py | 9 ++++++++- 2 files changed, 19 insertions(+), 15 deletions(-) diff --git a/uxarray/grid/angles.py b/uxarray/grid/angles.py index 94441cd54..bbe8d15d2 100644 --- a/uxarray/grid/angles.py +++ b/uxarray/grid/angles.py @@ -76,24 +76,21 @@ def _compute_face_node_angles_convex( return result -def _regular_polygon_angle(n, *, degrees=False): - """returns the internal angle for a regular polygon with n sides. - Equivalent to (n-2)*180/n degrees or (n-2)*pi/n radians. - - n can be a scalar or an array of integers (e.g., can use n = n_nodes_per_face) - """ - if degrees: - return (n - 2) * 180 / n - else: - return (n - 2) * np.pi / n - - def _compute_equiangle_skewness(face_node_angles, n_nodes_per_face): """Returns the equiangle skewness at each face: max((Amax - Areg) / (pi - Areg), (Amin - Areg) / Areg) where Amin, Amax = min, max of the angles at the nodes of the face - Areg = ideal angle for a regular polygon with n_nodes_per_face sides. + 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 ---------- @@ -110,7 +107,7 @@ def _compute_equiangle_skewness(face_node_angles, n_nodes_per_face): """ Amin = face_node_angles.min("n_max_face_nodes", skipna=True) Amax = face_node_angles.max("n_max_face_nodes", skipna=True) - Areg = _regular_polygon_angle(n_nodes_per_face, degrees=False) + Areg = face_node_angles.sum("n_max_face_nodes", skipna=True) / n_nodes_per_face term0 = (Amax - Areg) / (np.pi - Areg) term1 = (Amin - Areg) / Areg # Should just use np.maximum(term0, term1), but that drops UxDataArray type currently, diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py index 9547f897e..d9391618c 100644 --- a/uxarray/grid/grid.py +++ b/uxarray/grid/grid.py @@ -1978,12 +1978,19 @@ def copy(self): 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": computes the equiangular skewness of each face: + equiangle_skewness = max((Amax - Areg) / (pi - Areg), (Amin - Areg) / 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). From 2e085a2199e0ec6350c8d2053df04706da5c9859 Mon Sep 17 00:00:00 2001 From: Sam Evans <47793072+Sevans711@users.noreply.github.com> Date: Thu, 20 Aug 2026 14:37:24 -0400 Subject: [PATCH 4/4] fix skewness formula, add tests, add to api.rst also swaps to allowing "method" kwarg as a positional arg. (formula typo was Amin - Areg, but now fixed to Areg - Amin) --- docs/api.rst | 1 + test/grid/geometry/test_angles.py | 44 +++++++++++++++++++++++++++++++ uxarray/grid/angles.py | 4 +-- uxarray/grid/grid.py | 4 +-- 4 files changed, 49 insertions(+), 4 deletions(-) diff --git a/docs/api.rst b/docs/api.rst index 68b70bf9d..59064fc9c 100644 --- a/docs/api.rst +++ b/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 diff --git a/test/grid/geometry/test_angles.py b/test/grid/geometry/test_angles.py index d42d8f362..5987a2d04 100644 --- a/test/grid/geometry/test_angles.py +++ b/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()) diff --git a/uxarray/grid/angles.py b/uxarray/grid/angles.py index bbe8d15d2..a75fbcdec 100644 --- a/uxarray/grid/angles.py +++ b/uxarray/grid/angles.py @@ -78,7 +78,7 @@ def _compute_face_node_angles_convex( def _compute_equiangle_skewness(face_node_angles, n_nodes_per_face): """Returns the equiangle skewness at each face: - max((Amax - Areg) / (pi - Areg), (Amin - Areg) / Areg) + 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 @@ -109,7 +109,7 @@ def _compute_equiangle_skewness(face_node_angles, n_nodes_per_face): 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 = (Amin - Areg) / 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) diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py index d9391618c..e0eda811b 100644 --- a/uxarray/grid/grid.py +++ b/uxarray/grid/grid.py @@ -1976,7 +1976,7 @@ def copy(self): source_dims_dict=self._source_dims_dict, ) - def compute_skewness(self, *, method: str = "equiangle", as_uxarray: bool = False): + 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, @@ -1987,7 +1987,7 @@ def compute_skewness(self, *, method: str = "equiangle", as_uxarray: bool = Fals 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), (Amin - Areg) / Areg) + 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.