DRAFT: Build_edge_from_nodes method

docs/user_api/index.rst

@@ -167,6 +167,7 @@ Connectivity
    :toctree: _autosummary
 
    grid.connectivity.close_face_nodes
+   grid.connectivity.build_edge_from_nodes
 
 Coordinates
 -----------

test/test_grid.py

@@ -8,7 +8,8 @@
 
 import uxarray as ux
 
-from uxarray.grid.connectivity import _build_edge_node_connectivity, _build_face_edges_connectivity
+from uxarray.grid.connectivity import _build_edge_node_connectivity, _build_face_edges_connectivity, \
+    build_edge_from_nodes
 
 from uxarray.grid.coordinates import _populate_cartesian_xyz_coord, _populate_lonlat_coord
 
@@ -33,7 +34,6 @@
 
 
 class TestGrid(TestCase):
-
     grid_CSne30 = ux.open_grid(gridfile_CSne30)
     grid_RLL1deg = ux.open_grid(gridfile_RLL1deg)
     grid_RLL10deg_CSne4 = ux.open_grid(gridfile_RLL10deg_CSne4)
@@ -288,6 +288,16 @@ 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
 
+    def test_build_edge_from_nodes(self):
+        ugrid = current_path / "meshfiles" / "ugrid" / "ov_RLL10deg_CSne4" / "ov_RLL10deg_CSne4.ug"
+        mpas = current_path / "meshfiles" / "mpas" / "QU" / "mesh.QU.1920km.151026.nc"
+        exodus = current_path / "meshfiles" / "exodus" / "outCSne8" / "outCSne8.g"
+        scrip = current_path / "meshfiles" / "scrip" / "outCSne8" / "outCSne8.nc"
+        face_nodes_connectivity = np.array(
+            [np.array([[0, 0], [0, 10], [10, 0], [10, 10]])])
+        gridOpen = ux.open_grid(face_nodes_connectivity)
+        build_edge_from_nodes(gridOpen, grid_type='delaunay', plot=False)
+
 
 class TestOperators(TestCase):
     grid_CSne30_01 = ux.open_grid(gridfile_CSne30)
@@ -304,7 +314,6 @@ def test_ne(self):
 
 
 class TestFaceAreas(TestCase):
-
     grid_CSne30 = ux.open_grid(gridfile_CSne30)
 
     def test_calculate_total_face_area_triangle(self):
@@ -319,7 +328,7 @@ def test_calculate_total_face_area_triangle(self):
                                   islatlon=False,
                                   isconcave=False)
 
-        #calculate area
+        # calculate area
         area_gaussian = grid_verts.calculate_total_face_area(
             quadrature_rule="gaussian", order=5)
         nt.assert_almost_equal(area_gaussian, constants.TRI_AREA, decimal=3)

uxarray/grid/connectivity.py

@@ -1,5 +1,6 @@
 import numpy as np
 import xarray as xr
+from scipy.spatial import Delaunay, Voronoi
 
 from uxarray.constants import INT_DTYPE, INT_FILL_VALUE
 
@@ -230,3 +231,61 @@ def _build_face_edges_connectivity(grid):
             "start_index": INT_DTYPE(0),
             "long_name": "Maps every edge to the two nodes that it connects",
         })
+
+
+def build_edge_from_nodes(self, grid_type='delaunay', plot=False):
+    """Constructs edges from given node points using either delaunay
+    triangulation or voronoi diagram."""
+    x = self.Mesh2_node_x.data
+    y = self.Mesh2_node_y.data
+    points = np.column_stack((x, y))
+
+    if "Mesh2_face_nodes" in self._ds.variables:
+        self._ds = self._ds.drop_vars("Mesh2_face_nodes")
+
+    if grid_type == 'delaunay':
+        # Perform Delaunay triangulation
+        grid = Delaunay(points)
+
+        face_nodes = []
+        for simplex in grid.simplices:
+            face_nodes.append(simplex)
+        self._ds["Mesh2_face_nodes"] = (["nMesh2_face",
+                                         "nMaxMesh2_face_nodes"], face_nodes)
+    elif grid_type == 'voronoi':
+        # Perform Voronoi diagram construction
+        grid = Voronoi(points)
+
+        face_nodes = []
+        for ridge_vertices in grid.ridge_vertices:
+            if -1 not in ridge_vertices:  # Ignore infinite regions
+                face_nodes.append(ridge_vertices)
+        self._ds["Mesh2_face_nodes"] = (["nMesh2_face",
+                                         "nMaxMesh2_face_nodes"], face_nodes)
+    else:
+        raise ValueError("Invalid grid_type. Use 'delaunay' or 'voronoi'.")
+
+    # if plot:
+    # # Example data
+    # x = self.Mesh2_node_x.data
+    # y = self.Mesh2_node_y.data
+
+    # # Array of connections
+    # connections = self._ds["Mesh2_face_nodes"]
+
+    # # Plot the lines connecting the points
+    # for connection in connections:
+    # connection = np.append(connection, connection[0])
+    # x_points = x[connection]
+    # y_points = y[connection]
+    # plt.plot(x_points, y_points, color='blue')
+
+    # # Plot the individual points
+    # plt.scatter(x, y, color='red', marker='o', label='Points')
+
+    # # Add labels and title
+    # plt.xlabel('X-axis')
+    # plt.ylabel('Y-axis')
+    # plt.title('Lines Connecting Points')
+    # plt.legend()
+    # plt.show()