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 678f94cb4..a75fbcdec 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,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) diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py index 711b30389..e0eda811b 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,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, *,