Repository: UXARRAY/uxarray

Issue:
Add `Grid._build_edge_node_connectivity`
Add a function for constructing the connectivity variable `Mesh2_edge_nodes` as described in the UGRID conventions for a Flexible 2D Mesh Topology

Discussion:

Patch:
diff --git a/docs/internal_api/index.rst b/docs/internal_api/index.rst
index 734a59ff..4d585842 100644
--- a/docs/internal_api/index.rst
+++ b/docs/internal_api/index.rst
@@ -19,6 +19,7 @@ Grid Methods
    grid.Grid.__from_ds__
    grid.Grid.__from_vert__
    grid.Grid.__init_grid_var_attrs__
+   grid.Grid._build_edge_node_connectivity
    grid.Grid._populate_cartesian_xyz_coord
    grid.Grid._populate_lonlat_coord
    grid.Grid._build_face_dimension
diff --git a/docs/user_api/index.rst b/docs/user_api/index.rst
index ef2e1d9a..0a0837a2 100644
--- a/docs/user_api/index.rst
+++ b/docs/user_api/index.rst
@@ -40,3 +40,7 @@ Helper Functions
    helpers.get_tri_quadratureDG
    helpers.grid_center_lat_lon
    helpers.parse_grid_type
+   helpers.node_xyz_to_lonlat_rad
+   helpers.node_lonlat_rad_to_xyz
+   helpers.normalize_in_place
+   helpers.close_face_nodes
diff --git a/uxarray/grid.py b/uxarray/grid.py
index 869abe06..50f28710 100644
--- a/uxarray/grid.py
+++ b/uxarray/grid.py
@@ -9,7 +9,7 @@ from ._ugrid import _read_ugrid, _encode_ugrid
 from ._shapefile import _read_shpfile
 from ._scrip import _read_scrip, _encode_scrip
 from ._mpas import _read_mpas
-from .helpers import get_all_face_area_from_coords, parse_grid_type, _convert_node_xyz_to_lonlat_rad, _convert_node_lonlat_rad_to_xyz
+from .helpers import get_all_face_area_from_coords, parse_grid_type, node_xyz_to_lonlat_rad, node_lonlat_rad_to_xyz, close_face_nodes
 from .constants import INT_DTYPE, INT_FILL_VALUE
 
 
@@ -99,6 +99,11 @@ class Grid:
         # initialize convenience attributes
         self.__init_grid_var_attrs__()
 
+        # construct connectivity
+        if self.source_grid != "From vertices":
+            if "Mesh2_edge_nodes" not in self.ds:
+                self._build_edge_node_connectivity()
+
         # build face dimension, possibly safeguard for large datasets
         self._build_face_dimension()
 
@@ -428,6 +433,81 @@ class Grid:
 
         return integral
 
+    def _build_edge_node_connectivity(self):
+        """Constructs the UGRID connectivity variable (``Mesh2_edge_nodes``)
+        and stores it within the internal (``Grid.ds``) and through the
+        attribute (``Grid.Mesh2_edge_nodes``).
+
+        Additionally, the attributes (``inverse_indices``) and
+        (``fill_value_mask``) are stored for constructing other
+        connectivity variables.
+        """
+        padded_face_nodes = close_face_nodes(self.Mesh2_face_nodes.values,
+                                             self.nMesh2_face,
+                                             self.nMaxMesh2_face_nodes)
+
+        # array of empty edge nodes where each entry is a pair of indices
+        edge_nodes = np.empty((self.nMesh2_face * self.nMaxMesh2_face_nodes, 2),
+                              dtype=INT_DTYPE)
+
+        # first index includes starting node up to non-padded value
+        edge_nodes[:, 0] = padded_face_nodes[:, :-1].ravel()
+
+        # second index includes second node up to padded value
+        edge_nodes[:, 1] = padded_face_nodes[:, 1:].ravel()
+
+        # sorted edge nodes
+        edge_nodes.sort(axis=1)
+
+        # unique edge nodes
+        edge_nodes_unique, inverse_indices = np.unique(edge_nodes,
+                                                       return_inverse=True,
+                                                       axis=0)
+        # find all edge nodes that contain a fill value
+        fill_value_mask = np.logical_or(
+            edge_nodes_unique[:, 0] == INT_FILL_VALUE,
+            edge_nodes_unique[:, 1] == INT_FILL_VALUE)
+
+        # all edge nodes that do not contain a fill value
+        non_fill_value_mask = np.logical_not(fill_value_mask)
+        edge_nodes_unique = edge_nodes_unique[non_fill_value_mask]
+
+        # Update inverse_indices accordingly
+        indices_to_update = np.where(fill_value_mask)[0]
+
+        remove_mask = np.isin(inverse_indices, indices_to_update)
+        inverse_indices[remove_mask] = INT_FILL_VALUE
+
+        # Compute the indices where inverse_indices exceeds the values in indices_to_update
+        indexes = np.searchsorted(indices_to_update,
+                                  inverse_indices,
+                                  side='right')
+        # subtract the corresponding indexes from `inverse_indices`
+        for i in range(len(inverse_indices)):
+            if inverse_indices[i] != INT_FILL_VALUE:
+                inverse_indices[i] -= indexes[i]
+
+        # add Mesh2_edge_nodes to internal dataset
+        self.ds['Mesh2_edge_nodes'] = xr.DataArray(
+            edge_nodes_unique,
+            dims=["nMesh2_edge", "Two"],
+            attrs={
+                "cf_role":
+                    "edge_node_connectivity",
+                "long_name":
+                    "Maps every edge to the two nodes that it connects",
+                "start_index":
+                    INT_DTYPE(0),
+                "inverse_indices":
+                    inverse_indices,
+                "fill_value_mask":
+                    fill_value_mask
+            })
+
+        # set standardized attributes
+        setattr(self, "Mesh2_edge_nodes", self.ds['Mesh2_edge_nodes'])
+        setattr(self, "nMesh2_edge", edge_nodes_unique.shape[0])
+
     def _populate_cartesian_xyz_coord(self):
         """A helper function that populates the xyz attribute in UXarray.ds.
         This function is called when we need to use the cartesian coordinates
@@ -462,7 +542,7 @@ class Grid:
         nodes_lat_rad = np.deg2rad(self.Mesh2_node_y.values)
         nodes_rad = np.stack((nodes_lon_rad, nodes_lat_rad), axis=1)
         nodes_cart = np.asarray(
-            list(map(_convert_node_lonlat_rad_to_xyz, list(nodes_rad))))
+            list(map(node_lonlat_rad_to_xyz, list(nodes_rad))))
 
         self.ds["Mesh2_node_cart_x"] = xr.DataArray(
             data=nodes_cart[:, 0],
@@ -541,7 +621,7 @@ class Grid:
             (self.ds["Mesh2_node_x"].values, self.ds["Mesh2_node_y"].values,
              self.ds["Mesh2_node_z"].values),
             axis=1).tolist()
-        nodes_rad = list(map(_convert_node_xyz_to_lonlat_rad, nodes_cart))
+        nodes_rad = list(map(node_xyz_to_lonlat_rad, nodes_cart))
         nodes_degree = np.rad2deg(nodes_rad)
         self.ds["Mesh2_node_x"] = xr.DataArray(
             data=nodes_degree[:, 0],
diff --git a/uxarray/helpers.py b/uxarray/helpers.py
index 6512e428..aad303f6 100644
--- a/uxarray/helpers.py
+++ b/uxarray/helpers.py
@@ -6,7 +6,6 @@ from numba import njit, config
 import math
 
 from uxarray.constants import INT_DTYPE, INT_FILL_VALUE
-from .constants import INT_DTYPE
 
 config.DISABLE_JIT = False
 
@@ -155,14 +154,14 @@ def calculate_face_area(x,
 
         if (coords_type == "spherical"):
             node1 = np.array(
-                _convert_node_lonlat_rad_to_xyz(
-                    [np.deg2rad(x[0]), np.deg2rad(y[0])]))
+                node_lonlat_rad_to_xyz([np.deg2rad(x[0]),
+                                        np.deg2rad(y[0])]))
             node2 = np.array(
-                _convert_node_lonlat_rad_to_xyz(
+                node_lonlat_rad_to_xyz(
                     [np.deg2rad(x[j + 1]),
                      np.deg2rad(y[j + 1])]))
             node3 = np.array(
-                _convert_node_lonlat_rad_to_xyz(
+                node_lonlat_rad_to_xyz(
                     [np.deg2rad(x[j + 2]),
                      np.deg2rad(y[j + 2])]))
 
@@ -467,7 +466,7 @@ def grid_center_lat_lon(ds):
 
 
 @njit
-def _convert_node_lonlat_rad_to_xyz(node_coord):
+def node_lonlat_rad_to_xyz(node_coord):
     """Helper function to Convert the node coordinate from 2D
     longitude/latitude to normalized 3D xyz.
 
@@ -495,7 +494,7 @@ def _convert_node_lonlat_rad_to_xyz(node_coord):
 
 
 @njit
-def _convert_node_xyz_to_lonlat_rad(node_coord):
+def node_xyz_to_lonlat_rad(node_coord):
     """Calculate the latitude and longitude in radiance for a node represented
     in the [x, y, z] 3D Cartesian coordinates.
 
@@ -517,7 +516,7 @@ def _convert_node_xyz_to_lonlat_rad(node_coord):
     if len(node_coord) != 3:
         raise RuntimeError("Input array should have a length of 3: [x, y, z]")
     reference_tolerance = 1.0e-12
-    [dx, dy, dz] = _normalize_in_place(node_coord)
+    [dx, dy, dz] = normalize_in_place(node_coord)
     dx /= np.absolute(dx * dx + dy * dy + dz * dz)
     dy /= np.absolute(dx * dx + dy * dy + dz * dz)
     dz /= np.absolute(dx * dx + dy * dy + dz * dz)
@@ -539,7 +538,7 @@ def _convert_node_xyz_to_lonlat_rad(node_coord):
 
 
 @njit
-def _normalize_in_place(node):
+def normalize_in_place(node):
     """Helper function to project an arbitrary node in 3D coordinates [x, y, z]
     on the unit sphere. It uses the `np.linalg.norm` internally to calculate
     the magnitude.
@@ -562,7 +561,7 @@ def _normalize_in_place(node):
     if len(node) != 3:
         raise RuntimeError("Input array should have a length of 3: [x, y, z]")
 
-    return list(np.array(node) / np.linalg.norm(np.array(node), ord=2))
+    return np.array(node) / np.linalg.norm(np.array(node), ord=2)
 
 
 def _replace_fill_values(grid_var, original_fill, new_fill, new_dtype=None):
@@ -622,3 +621,54 @@ def _replace_fill_values(grid_var, original_fill, new_fill, new_dtype=None):
     grid_var[fill_val_idx] = new_fill
 
     return grid_var
+
+
+def close_face_nodes(Mesh2_face_nodes, nMesh2_face, nMaxMesh2_face_nodes):
+    """Closes (``Mesh2_face_nodes``) by inserting the first node index after
+    the last non-fill-value node.
+
+    Parameters
+    ----------
+    Mesh2_face_nodes : np.ndarray
+        Connectivity array for constructing a face from its nodes
+    nMesh2_face : constant
+        Number of faces
+    nMaxMesh2_face_nodes : constant
+        Max number of nodes that compose a face
+
+    Returns
+    ----------
+    closed : ndarray
+        Closed (padded) Mesh2_face_nodes
+
+    Example
+    ----------
+    Given face nodes with shape [2 x 5]
+        [0, 1, 2, 3, FILL_VALUE]
+        [4, 5, 6, 7, 8]
+    Pads them to the following with shape [2 x 6]
+        [0, 1, 2, 3, 0, FILL_VALUE]
+        [4, 5, 6, 7, 8, 4]
+    """
+
+    # padding to shape [nMesh2_face, nMaxMesh2_face_nodes + 1]
+    closed = np.ones((nMesh2_face, nMaxMesh2_face_nodes + 1),
+                     dtype=INT_DTYPE) * INT_FILL_VALUE
+
+    # set all non-paded values to original face nodee values
+    closed[:, :-1] = Mesh2_face_nodes.copy()
+
+    # instance of first fill value
+    first_fv_idx_2d = np.argmax(closed == INT_FILL_VALUE, axis=1)
+
+    # 2d to 1d index for np.put()
+    first_fv_idx_1d = first_fv_idx_2d + (
+        (nMaxMesh2_face_nodes + 1) * np.arange(0, nMesh2_face))
+
+    # column of first node values
+    first_node_value = Mesh2_face_nodes[:, 0].copy()
+
+    # insert first node column at occurrence of first fill value
+    np.put(closed.ravel(), first_fv_idx_1d, first_node_value)
+
+    return closed


Test patch:
diff --git a/test/test_grid.py b/test/test_grid.py
index 80cab545..5b3b7fd9 100644
--- a/test/test_grid.py
+++ b/test/test_grid.py
@@ -271,12 +271,6 @@ class TestGrid(TestCase):
         self.assertEqual(n_faces, grid.nMesh2_face)
         self.assertEqual(n_face_nodes, grid.nMaxMesh2_face_nodes)
 
-    # def test_init_dimension_attrs(self):
-
-    # TODO: Move to test_shpfile/scrip when implemented
-    # use external package to read?
-    # https://gis.stackexchange.com/questions/113799/how-to-read-a-shapefile-in-python
-
     def test_read_shpfile(self):
         """Reads a shape file and write ugrid file."""
         with self.assertRaises(RuntimeError):
@@ -486,3 +480,55 @@ class TestPopulateCoordinates(TestCase):
             nt.assert_almost_equal(vgrid.ds["Mesh2_node_y"].values[i],
                                    lat_deg[i],
                                    decimal=12)
+
+
+class TestConnectivity(TestCase):
+    mpas_filepath = current_path / "meshfiles" / "mpas" / "QU" / "mesh.QU.1920km.151026.nc"
+    exodus_filepath = current_path / "meshfiles" / "exodus" / "outCSne8" / "outCSne8.g"
+    ugrid_filepath_01 = current_path / "meshfiles" / "ugrid" / "outCSne30" / "outCSne30.ug"
+    ugrid_filepath_02 = current_path / "meshfiles" / "ugrid" / "outRLL1deg" / "outRLL1deg.ug"
+    ugrid_filepath_03 = current_path / "meshfiles" / "ugrid" / "ov_RLL10deg_CSne4" / "ov_RLL10deg_CSne4.ug"
+
+    def test_build_edge_nodes(self):
+        """Tests the construction of (``Mesh2_edge_nodes``) on an MPAS grid
+        with known edge nodes."""
+
+        # grid with known edge node connectivity
+        mpas_grid_xr = xr.open_dataset(self.mpas_filepath)
+        mpas_grid_ux = ux.Grid(mpas_grid_xr)
+        edge_nodes_expected = mpas_grid_ux.ds['Mesh2_edge_nodes'].values
+
+        # arrange edge nodes in the same manner as Grid._build_edge_node_connectivity
+        edge_nodes_expected.sort(axis=1)
+        edge_nodes_expected = np.unique(edge_nodes_expected, axis=0)
+
+        # construct edge nodes
+        mpas_grid_ux._build_edge_node_connectivity()
+        edge_nodes_output = mpas_grid_ux.ds['Mesh2_edge_nodes'].values
+
+        assert np.array_equal(edge_nodes_expected, edge_nodes_output)
+
+        # euler's formula (n_face = n_edges - n_nodes + 2)
+        n_face = mpas_grid_ux.nMesh2_node
+        n_node = mpas_grid_ux.nMesh2_face
+        n_edge = edge_nodes_output.shape[0]
+
+        assert (n_face == n_edge - n_node + 2)
+
+    def test_edge_nodes_euler(self):
+        """Verifies that (``nMesh2_edge``) follows euler's formula."""
+        grid_paths = [
+            self.exodus_filepath, self.ugrid_filepath_01,
+            self.ugrid_filepath_02, self.ugrid_filepath_03
+        ]
+
+        for grid_path in grid_paths:
+            grid_xr = xr.open_dataset(grid_path)
+            grid_ux = ux.Grid(grid_xr)
+
+            n_face = grid_ux.nMesh2_node
+            n_node = grid_ux.nMesh2_face
+            n_edge = grid_ux.nMesh2_edge
+
+            # euler's formula (n_face = n_edges - n_nodes + 2)
+            assert (n_face == n_edge - n_node + 2)
diff --git a/test/test_helpers.py b/test/test_helpers.py
index 905d2d4a..fd6a748c 100644
--- a/test/test_helpers.py
+++ b/test/test_helpers.py
@@ -95,33 +95,31 @@ class TestGridCenter(TestCase):
 class TestCoordinatesConversion(TestCase):
 
     def test_normalize_in_place(self):
-        [x, y, z] = ux.helpers._normalize_in_place(
+        [x, y, z] = ux.helpers.normalize_in_place(
             [random.random(), random.random(),
              random.random()])
         self.assertLessEqual(np.absolute(np.sqrt(x * x + y * y + z * z) - 1),
                              err_tolerance)
 
-    def test_convert_node_xyz_to_lonlat_rad(self):
-        [x, y, z] = ux.helpers._normalize_in_place([
+    def test_node_xyz_to_lonlat_rad(self):
+        [x, y, z] = ux.helpers.normalize_in_place([
             random.uniform(-1, 1),
             random.uniform(-1, 1),
             random.uniform(-1, 1)
         ])
-        [lon, lat] = ux.helpers._convert_node_xyz_to_lonlat_rad([x, y, z])
-        [new_x, new_y,
-         new_z] = ux.helpers._convert_node_lonlat_rad_to_xyz([lon, lat])
+        [lon, lat] = ux.helpers.node_xyz_to_lonlat_rad([x, y, z])
+        [new_x, new_y, new_z] = ux.helpers.node_lonlat_rad_to_xyz([lon, lat])
         self.assertLessEqual(np.absolute(new_x - x), err_tolerance)
         self.assertLessEqual(np.absolute(new_y - y), err_tolerance)
         self.assertLessEqual(np.absolute(new_z - z), err_tolerance)
 
-    def test_convert_node_latlon_rad_to_xyz(self):
+    def test_node_latlon_rad_to_xyz(self):
         [lon, lat] = [
             random.uniform(0, 2 * np.pi),
             random.uniform(-0.5 * np.pi, 0.5 * np.pi)
         ]
-        [x, y, z] = ux.helpers._convert_node_lonlat_rad_to_xyz([lon, lat])
-        [new_lon,
-         new_lat] = ux.helpers._convert_node_xyz_to_lonlat_rad([x, y, z])
+        [x, y, z] = ux.helpers.node_lonlat_rad_to_xyz([lon, lat])
+        [new_lon, new_lat] = ux.helpers.node_xyz_to_lonlat_rad([x, y, z])
         self.assertLessEqual(np.absolute(new_lon - lon), err_tolerance)
         self.assertLessEqual(np.absolute(new_lat - lat), err_tolerance)