From 5be721e4312ee247fa833f6fbacb28b016563af0 Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Date: Tue, 6 Aug 2024 11:11:16 -0500 Subject: [PATCH 1/9] add check for normalization --- uxarray/grid/grid.py | 34 ++++++++++++++++++++++++++++++++++ uxarray/grid/validation.py | 32 ++++++++++++++++++++++++-------- 2 files changed, 58 insertions(+), 8 deletions(-) diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py index 5b92c1026..841661e10 100644 --- a/uxarray/grid/grid.py +++ b/uxarray/grid/grid.py @@ -37,6 +37,7 @@ _set_desired_longitude_range, _populate_node_latlon, _populate_node_xyz, + _normalize_xyz, ) from uxarray.grid.connectivity import ( _populate_edge_node_connectivity, @@ -70,6 +71,7 @@ _check_connectivity, _check_duplicate_nodes, _check_area, + _check_normalization, ) from xarray.core.utils import UncachedAccessor @@ -1158,6 +1160,38 @@ def compute_face_areas( return self._face_areas, self._face_jacobian + def normalize_cartesian_coordinates(self): + """Normalizes Cartesian coordinates.""" + + if _check_normalization(self): + # check if coordinates are already normalized + return + + if "node_x" in self._ds: + # normalize node coordinates + node_x, node_y, node_z = _normalize_xyz( + self.node_x.values, self.node_y.values, self.node_z.values + ) + self.node_x.data = node_x + self.node_y.data = node_y + self.node_z.data = node_z + if "edge_x" in self._ds: + # normalize edge coordinates + edge_x, edge_y, edge_z = _normalize_xyz( + self.edge_x.values, self.edge_y.values, self.edge_z.values + ) + self.edge_x.data = edge_x + self.edge_y.data = edge_y + self.edge_z.data = edge_z + if "face_x" in self._ds: + # normalize face coordinates + face_x, face_y, face_z = _normalize_xyz( + self.face_x.values, self.face_y.values, self.face_z.values + ) + self.face_x.data = face_x + self.face_y.data = face_y + self.face_z.data = face_z + def to_xarray(self, grid_format: Optional[str] = "ugrid"): """Returns a xarray Dataset representation in a specific grid format from the Grid object. diff --git a/uxarray/grid/validation.py b/uxarray/grid/validation.py index 1a4852711..f4c212d2c 100644 --- a/uxarray/grid/validation.py +++ b/uxarray/grid/validation.py @@ -5,7 +5,7 @@ # validation helper functions -def _check_connectivity(self): +def _check_connectivity(grid): """Check if all nodes are referenced by at least one element. If not, the mesh may have hanging nodes and may not a valid UGRID @@ -14,28 +14,28 @@ def _check_connectivity(self): # Check if all nodes are referenced by at least one element # get unique nodes in connectivity - nodes_in_conn = np.unique(self.face_node_connectivity.values.flatten()) + nodes_in_conn = np.unique(grid.face_node_connectivity.values.flatten()) # remove negative indices/fill values from the list nodes_in_conn = nodes_in_conn[nodes_in_conn >= 0] # check if the size of unique nodes in connectivity is equal to the number of nodes - if nodes_in_conn.size == self.n_node: + if nodes_in_conn.size == grid.n_node: print("-All nodes are referenced by at least one element.") return True else: warn( "Some nodes may not be referenced by any element. {0} and {1}".format( - nodes_in_conn.size, self.n_node + nodes_in_conn.size, grid.n_node ), RuntimeWarning, ) return False -def _check_duplicate_nodes(self): +def _check_duplicate_nodes(grid): """Check if there are duplicate nodes in the mesh.""" - coords1 = np.column_stack((np.vstack(self.node_lon), np.vstack(self.node_lat))) + coords1 = np.column_stack((np.vstack(grid.node_lon), np.vstack(grid.node_lat))) unique_nodes, indices = np.unique(coords1, axis=0, return_index=True) duplicate_indices = np.setdiff1d(np.arange(len(coords1)), indices) @@ -52,9 +52,9 @@ def _check_duplicate_nodes(self): return True -def _check_area(self): +def _check_area(grid): """Check if each face area is greater than our constant ERROR_TOLERANCE.""" - areas = self.face_areas + areas = grid.face_areas # Check if area of any face is close to zero if np.any(np.isclose(areas, 0, atol=ERROR_TOLERANCE)): warn( @@ -65,3 +65,19 @@ def _check_area(self): else: print("-No face area is close to zero.") return True + + +def _check_normalization(grid): + """Checks whether all the cartesiain coordinates are normalized.""" + + if "node_x" in grid._ds: + if not (grid.node_x**2 + grid.node_y**2 + grid.node_z**2).all(): + return False + if "edge_x" in grid._ds: + if not (grid.edge_x**2 + grid.edge_y**2 + grid.edge_z**2).all(): + return False + if "face_x" in grid._ds: + if not (grid.face_x**2 + grid.face_y**2 + grid.face_z**2).all(): + return False + + return True From 03336919ba8e908b7c61e98c60cb3c47cd1a1638 Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Date: Fri, 6 Sep 2024 14:00:52 -0500 Subject: [PATCH 2/9] add tests, benchmark --- benchmarks/mpas_ocean.py | 15 +++++++++++++++ test/test_grid.py | 10 ++++++++++ 2 files changed, 25 insertions(+) diff --git a/benchmarks/mpas_ocean.py b/benchmarks/mpas_ocean.py index d2baaf86b..eb256303c 100644 --- a/benchmarks/mpas_ocean.py +++ b/benchmarks/mpas_ocean.py @@ -164,3 +164,18 @@ def time_nearest_neighbor_remapping(self): def time_inverse_distance_weighted_remapping(self): self.uxds_480["bottomDepth"].remap.inverse_distance_weighted(self.uxds_120.uxgrid) + + +class ConstructTreeStructures: + param_names = ['resolution'] + params = ['480km', '120km'] + + def setup(self, resolution): + self.uxgrid = ux.open_grid(file_path_dict[resolution][0]) + + def teardown(self, resolution): + del self.uxgrid + + def time_check_norm(self, resolution): + from uxarray.grid.validation import _check_normalization + _check_normalization(self.uxgrid) diff --git a/test/test_grid.py b/test/test_grid.py index 10faba696..448aca831 100644 --- a/test/test_grid.py +++ b/test/test_grid.py @@ -957,3 +957,13 @@ def test_populate_bounds_MPAS(self): uxgrid = ux.Grid.from_dataset(xrds, use_dual=True) bounds_xarray = uxgrid.bounds pass + + +class TestNormalizeExistingCoordinates(TestCase): + gridfile_mpas = current_path / "meshfiles" / "mpas" / "QU" / "mesh.QU.1920km.151026.nc" + + def test_mpas_norm(self): + from uxarray.grid.validation import _check_normalization + uxgrid = ux.open_grid(self.gridfile_mpas) + + _check_normalization(uxgrid) From b2a1094aaedc3bc7949904f364afec7e11c047da Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Date: Fri, 6 Sep 2024 14:21:47 -0500 Subject: [PATCH 3/9] update norm check, update tests --- test/test_grid.py | 24 ++++++++++++++++++++++-- uxarray/grid/validation.py | 12 +++++++++--- 2 files changed, 31 insertions(+), 5 deletions(-) diff --git a/test/test_grid.py b/test/test_grid.py index 448aca831..ae1ff11fb 100644 --- a/test/test_grid.py +++ b/test/test_grid.py @@ -961,9 +961,29 @@ def test_populate_bounds_MPAS(self): class TestNormalizeExistingCoordinates(TestCase): gridfile_mpas = current_path / "meshfiles" / "mpas" / "QU" / "mesh.QU.1920km.151026.nc" + gridfile_CSne30 = current_path / "meshfiles" / "ugrid" / "outCSne30" / "outCSne30.ug" - def test_mpas_norm(self): + def test_non_norm_initial(self): + """Check the normalization of coordinates that were initially parsed as + non-normalized.""" from uxarray.grid.validation import _check_normalization uxgrid = ux.open_grid(self.gridfile_mpas) - _check_normalization(uxgrid) + # Make the coordinates not normalized + uxgrid.node_x.data = 5 * uxgrid.node_x.data + uxgrid.node_y.data = 5 * uxgrid.node_y.data + uxgrid.node_z.data = 5 * uxgrid.node_z.data + assert not _check_normalization(uxgrid) + + uxgrid.normalize_cartesian_coordinates() + + assert _check_normalization(uxgrid) + + def test_norm_initial(self): + """Coordinates should be normalized for grids that we construct + them.""" + from uxarray.grid.validation import _check_normalization + uxgrid = ux.open_grid(self.gridfile_CSne30) + + assert _check_normalization(uxgrid) + pass diff --git a/uxarray/grid/validation.py b/uxarray/grid/validation.py index f4c212d2c..f1a26f21f 100644 --- a/uxarray/grid/validation.py +++ b/uxarray/grid/validation.py @@ -71,13 +71,19 @@ def _check_normalization(grid): """Checks whether all the cartesiain coordinates are normalized.""" if "node_x" in grid._ds: - if not (grid.node_x**2 + grid.node_y**2 + grid.node_z**2).all(): + if not ( + np.isclose((grid.node_x**2 + grid.node_y**2 + grid.node_z**2), 1.0) + ).all(): return False if "edge_x" in grid._ds: - if not (grid.edge_x**2 + grid.edge_y**2 + grid.edge_z**2).all(): + if not ( + np.isclose((grid.node_x**2 + grid.node_y**2 + grid.node_z**2), 1.0) + ).all(): return False if "face_x" in grid._ds: - if not (grid.face_x**2 + grid.face_y**2 + grid.face_z**2).all(): + if not ( + np.isclose((grid.node_x**2 + grid.node_y**2 + grid.node_z**2), 1.0) + ).all(): return False return True From f293e6793569bb0e81386c875b7119aaecf79e79 Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Date: Fri, 6 Sep 2024 16:11:10 -0500 Subject: [PATCH 4/9] add flag to keep track of norm --- uxarray/grid/grid.py | 3 +++ uxarray/grid/validation.py | 9 +++++++++ 2 files changed, 12 insertions(+) diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py index c3722ed0c..ed74f5aa2 100644 --- a/uxarray/grid/grid.py +++ b/uxarray/grid/grid.py @@ -176,6 +176,9 @@ def __init__( self._ball_tree = None self._kd_tree = None + # flag to track if coordinates are normalized + self._normalized = None + # set desired longitude range to [-180, 180] _set_desired_longitude_range(self._ds) diff --git a/uxarray/grid/validation.py b/uxarray/grid/validation.py index f1a26f21f..d24bee138 100644 --- a/uxarray/grid/validation.py +++ b/uxarray/grid/validation.py @@ -70,20 +70,29 @@ def _check_area(grid): def _check_normalization(grid): """Checks whether all the cartesiain coordinates are normalized.""" + if grid._normalized is not None: + return grid._normalized + if "node_x" in grid._ds: if not ( np.isclose((grid.node_x**2 + grid.node_y**2 + grid.node_z**2), 1.0) ).all(): + grid._normalized = False return False if "edge_x" in grid._ds: if not ( np.isclose((grid.node_x**2 + grid.node_y**2 + grid.node_z**2), 1.0) ).all(): + grid._normalized = False return False if "face_x" in grid._ds: if not ( np.isclose((grid.node_x**2 + grid.node_y**2 + grid.node_z**2), 1.0) ).all(): + grid._normalized = False return False + # set the grid as normalized + grid._normalized = True + return True From caf0f04cf87cb071bcb51c65fa739fa5bdb9f5da Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Date: Fri, 6 Sep 2024 16:12:19 -0500 Subject: [PATCH 5/9] fix logic and add comment --- uxarray/grid/validation.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/uxarray/grid/validation.py b/uxarray/grid/validation.py index d24bee138..7e1bc91aa 100644 --- a/uxarray/grid/validation.py +++ b/uxarray/grid/validation.py @@ -70,7 +70,8 @@ def _check_area(grid): def _check_normalization(grid): """Checks whether all the cartesiain coordinates are normalized.""" - if grid._normalized is not None: + if grid._normalized is True: + # grid is already normalized, no need to run extra checks return grid._normalized if "node_x" in grid._ds: From d208cf66978062a1f9fc55e101293a176858944a Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Date: Fri, 6 Sep 2024 16:13:37 -0500 Subject: [PATCH 6/9] add error tollerance --- uxarray/grid/validation.py | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/uxarray/grid/validation.py b/uxarray/grid/validation.py index 7e1bc91aa..b42b8cbb3 100644 --- a/uxarray/grid/validation.py +++ b/uxarray/grid/validation.py @@ -76,19 +76,31 @@ def _check_normalization(grid): if "node_x" in grid._ds: if not ( - np.isclose((grid.node_x**2 + grid.node_y**2 + grid.node_z**2), 1.0) + np.isclose( + (grid.node_x**2 + grid.node_y**2 + grid.node_z**2), + 1.0, + atol=ERROR_TOLERANCE, + ) ).all(): grid._normalized = False return False if "edge_x" in grid._ds: if not ( - np.isclose((grid.node_x**2 + grid.node_y**2 + grid.node_z**2), 1.0) + np.isclose( + (grid.node_x**2 + grid.node_y**2 + grid.node_z**2), + 1.0, + atol=ERROR_TOLERANCE, + ) ).all(): grid._normalized = False return False if "face_x" in grid._ds: if not ( - np.isclose((grid.node_x**2 + grid.node_y**2 + grid.node_z**2), 1.0) + np.isclose( + (grid.node_x**2 + grid.node_y**2 + grid.node_z**2), + 1.0, + atol=ERROR_TOLERANCE, + ) ).all(): grid._normalized = False return False From 7f033ecbd003b2e1e2f408247a05e6099b88b3c2 Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Date: Mon, 16 Sep 2024 20:59:55 -0500 Subject: [PATCH 7/9] Update test_grid.py --- test/test_grid.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/test/test_grid.py b/test/test_grid.py index 0ed631e82..044dcf464 100644 --- a/test/test_grid.py +++ b/test/test_grid.py @@ -281,7 +281,6 @@ def test_read_scrip(self): # Test read from scrip and from ugrid for grid class grid_CSne8 = ux.open_grid(gridfile_CSne8) # tests from scrip - pass class TestOperators(TestCase): @@ -926,7 +925,6 @@ def test_from_dataset(self): xrds = xr.open_dataset(self.gridfile_scrip) uxgrid = ux.Grid.from_dataset(xrds) - pass def test_from_face_vertices(self): single_face_latlon = [(0.0, 90.0), (-180, 0.0), (0.0, -90)] @@ -965,7 +963,6 @@ def test_populate_bounds_MPAS(self): xrds = xr.open_dataset(self.gridfile_mpas) uxgrid = ux.Grid.from_dataset(xrds, use_dual=True) bounds_xarray = uxgrid.bounds - pass class TestNormalizeExistingCoordinates(TestCase): From a39fba31f2e25a245861bb955ea028d40a09a5ea Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Date: Mon, 16 Sep 2024 21:00:36 -0500 Subject: [PATCH 8/9] remove pass --- test/test_grid.py | 1 - 1 file changed, 1 deletion(-) diff --git a/test/test_grid.py b/test/test_grid.py index 044dcf464..57499248e 100644 --- a/test/test_grid.py +++ b/test/test_grid.py @@ -992,4 +992,3 @@ def test_norm_initial(self): uxgrid = ux.open_grid(self.gridfile_CSne30) assert _check_normalization(uxgrid) - pass From 77ebb94bb8b723b2835ce34803223bd5c4ca4aa8 Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Date: Mon, 16 Sep 2024 21:09:00 -0500 Subject: [PATCH 9/9] update test --- test/test_grid.py | 3 +-- uxarray/grid/validation.py | 1 - 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/test/test_grid.py b/test/test_grid.py index 57499248e..a31223ffe 100644 --- a/test/test_grid.py +++ b/test/test_grid.py @@ -960,8 +960,7 @@ def test_populate_bounds_GCA_mix(self): nt.assert_allclose(grid.bounds.values, expected_bounds, atol=ERROR_TOLERANCE) def test_populate_bounds_MPAS(self): - xrds = xr.open_dataset(self.gridfile_mpas) - uxgrid = ux.Grid.from_dataset(xrds, use_dual=True) + uxgrid = ux.open_grid(self.gridfile_mpas) bounds_xarray = uxgrid.bounds diff --git a/uxarray/grid/validation.py b/uxarray/grid/validation.py index eefa7561f..b42b8cbb3 100644 --- a/uxarray/grid/validation.py +++ b/uxarray/grid/validation.py @@ -1,7 +1,6 @@ import numpy as np from warnings import warn - from uxarray.constants import ERROR_TOLERANCE