Merge pull request #1052 from UXARRAY/zedwick/to_planar

Regional Delaunay

Author: rajeeja

docs/user-guide/from-points.ipynb

@@ -80,7 +80,7 @@
     "uxgrid_global = ux.open_grid(\"../../test/meshfiles/ugrid/outCSne30/outCSne30.ug\")\n",
     "uxgrid_global_ocean = ux.open_grid(\"../../test/meshfiles/mpas/QU/oQU480.231010.nc\")\n",
     "uxgrid_global_ocean.normalize_cartesian_coordinates()\n",
-    "uxgrid_regional = uxgrid_global.subset.nearest_neighbor((0.0, 0.0), k=9)\n",
+    "uxgrid_regional = uxgrid_global.subset.nearest_neighbor((0.0, 0.0), k=50)\n",
     "\n",
     "(\n",
     "    uxgrid_global.plot.face_centers(\n",
@@ -723,9 +723,105 @@
    "source": [
     "## Regional Data\n",
     "\n",
-    "```{warning}\n",
-    "Constructing a grid from regional point data is not yet supported.\n",
-    "```"
+    "The regional delaunay method can be used to construct a grid from points in a regional area.\n",
+    "\n",
+    "### How It Works\n",
+    "\n",
+    "1. **Input Points on the Sphere**:\n",
+    "   - The method accepts input points defined in spherical coordinates (e.g., latitude and longitude) or Cartesian coordinates (x, y, z) that lie on the surface of the sphere. They are internally converted to normalized Cartesian coordinates.\n",
+    "\n",
+    "2. **Computing the Regional Delaunay Diagram**:\n",
+    "   - The method projects the points to a 2D plane using stereographic projection, followed by SciPy's `Delaunay` triangulation method to construct the grid.\n",
+    "\n",
+    "3. **Extracting Triangles**:\n",
+    "   - Once the triangles of 2D points are determined, the connectivity of the triangular faces are extracted. These triangles represent the Delaunay triangulation on the sphere's surface, ensuring that no point is inside the circumcircle of any triangle, which is a key property of Delaunay triangulations."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5e15b602-70e1-434d-a0a4-51d5913abb6a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "grid_r = ux.Grid.from_points(points_regional, method=\"regional_delaunay\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6e4d838a-4f7f-474f-a2f9-47569df02a91",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "grid_r.plot(\n",
+    "    linewidth=0.5,\n",
+    "    periodic_elements=\"exclude\",\n",
+    "    height=500,\n",
+    "    width=1000,\n",
+    "    title=\"Regional Delaunay Regions\",\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "03dbafcc-a603-466d-b886-adb15764bd4c",
+   "metadata": {},
+   "source": [
+    "### Antimerdian\n",
+    "\n",
+    "This also works on regions wrapping the antimerdian"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "95e37aa6-6591-47b0-8a64-9d38f9dec664",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "antimerdian_region = uxgrid_global.subset.bounding_circle((-180, 0), 10)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2d630797-b2c8-46dc-bda2-60a69ed266d6",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "x_antimerdian_region, y_antimerdian_region, z_antimerdian_region = (\n",
+    "    antimerdian_region.face_x.values,\n",
+    "    antimerdian_region.face_y.values,\n",
+    "    antimerdian_region.face_z.values,\n",
+    ")\n",
+    "antimerdian_region = (x_antimerdian_region, y_antimerdian_region, z_antimerdian_region)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8f5e70c5-c31e-4183-af2a-24ed1dcaba94",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "grid_r = ux.Grid.from_points(antimerdian_region, method=\"regional_delaunay\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a83c11a8-fbc9-4dc9-ba10-8c3babdae83f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "grid_r.plot(\n",
+    "    linewidth=0.5,\n",
+    "    periodic_elements=\"exclude\",\n",
+    "    height=500,\n",
+    "    width=1000,\n",
+    "    title=\"Regional Delaunay Regions\",\n",
+    ")"
    ]
   }
  ],
@@ -745,7 +841,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.12.0"
+   "version": "3.12.7"
   }
  },
  "nbformat": 4,

test/test_from_points.py

@@ -25,6 +25,23 @@ def test_spherical_delaunay():
     assert uxgrid_dt_xyz.triangular
     assert uxgrid_dt_latlon.triangular
 
+
+def test_regional_delaunay():
+    uxgrid = ux.open_grid(grid_path)
+
+    uxgrid_regional = uxgrid.subset.nearest_neighbor((0.0, 0.0), k=50)
+
+    points_xyz = (uxgrid_regional.node_x.values, uxgrid_regional.node_y.values, uxgrid_regional.node_z.values)
+    points_latlon = (uxgrid_regional.node_lon.values, uxgrid_regional.node_lat.values)
+
+    uxgrid_dt_xyz = ux.Grid.from_points(points_xyz, method='regional_delaunay')
+    uxgrid_dt_latlon = ux.Grid.from_points(points_latlon, method='regional_delaunay')
+
+    assert uxgrid_dt_xyz.n_node == uxgrid_dt_latlon.n_node == len(points_xyz[0])
+    assert uxgrid_dt_xyz.triangular
+    assert uxgrid_dt_latlon.triangular
+
+
 def test_spherical_voronoi():
     uxgrid = ux.open_grid(grid_path)
     points_xyz = (uxgrid.node_x.values, uxgrid.node_y.values, uxgrid.node_z.values)

test/test_geometry.py

@@ -12,7 +12,8 @@
 from uxarray.grid.coordinates import _populate_node_latlon, _lonlat_rad_to_xyz, _normalize_xyz, _xyz_to_lonlat_rad
 from uxarray.grid.arcs import extreme_gca_latitude
 from uxarray.grid.utils import _get_cartesian_face_edge_nodes, _get_lonlat_rad_face_edge_nodes
-from uxarray.grid.geometry import _populate_face_latlon_bound, _populate_bounds
+from uxarray.grid.geometry import _populate_face_latlon_bound, _populate_bounds, stereographic_projection, \
+    inverse_stereographic_projection
 
 from spatialpandas.geometry import MultiPolygon
 
@@ -27,14 +28,14 @@
 grid_files = [gridfile_CSne8, gridfile_geoflow]
 data_files = [datafile_CSne30, datafile_geoflow]
 
-grid_quad_hex = current_path/ "meshfiles" / "ugrid" / "quad-hexagon" / "grid.nc"
-grid_geoflow = current_path/ "meshfiles" / "ugrid" / "geoflow-small" / "grid.nc"
-grid_mpas = current_path/ "meshfiles" / "mpas" / "QU" / "oQU480.231010.nc"
-
+grid_quad_hex = current_path / "meshfiles" / "ugrid" / "quad-hexagon" / "grid.nc"
+grid_geoflow = current_path / "meshfiles" / "ugrid" / "geoflow-small" / "grid.nc"
+grid_mpas = current_path / "meshfiles" / "mpas" / "QU" / "oQU480.231010.nc"
 
 # List of grid files to test
 grid_files_latlonBound = [grid_quad_hex, grid_geoflow, gridfile_CSne8, grid_mpas]
 
+
 class TestAntimeridian(TestCase):
 
     def test_crossing(self):
@@ -63,8 +64,6 @@ def test_linecollection_execution(self):
         lines = uxgrid.to_linecollection()
 
 
-
-
 class TestPredicate(TestCase):
 
     def test_pole_point_inside_polygon_from_vertice_north(self):
@@ -377,7 +376,7 @@ def test_extreme_gca_latitude_max(self):
 
     def test_extreme_gca_latitude_max_short(self):
         # Define a great circle arc in 3D space that has a small span
-        gca_cart = np.array( [[ 0.65465367, -0.37796447, -0.65465367], [ 0.6652466,  -0.33896007, -0.6652466 ]])
+        gca_cart = np.array([[0.65465367, -0.37796447, -0.65465367], [0.6652466, -0.33896007, -0.6652466]])
 
         # Calculate the maximum latitude
         max_latitude = ux.grid.arcs.extreme_gca_latitude(gca_cart, 'max')
@@ -463,7 +462,7 @@ def test_insert_pt_in_empty_state(self):
 class TestLatlonBoundsGCA(TestCase):
 
     def _get_cartesian_face_edge_nodes_testcase_helper(
-            self,face_nodes_ind, face_edges_ind, edge_nodes_grid, node_x, node_y, node_z
+            self, face_nodes_ind, face_edges_ind, edge_nodes_grid, node_x, node_y, node_z
     ):
         """This function is only used to help generating the testcase and
         should not be used in the actual implementation. Construct an array to
@@ -522,7 +521,7 @@ def _get_cartesian_face_edge_nodes_testcase_helper(
         return cartesian_coordinates
 
     def _get_lonlat_rad_face_edge_nodes_testcase_helper(
-            self,face_nodes_ind, face_edges_ind, edge_nodes_grid, node_lon, node_lat
+            self, face_nodes_ind, face_edges_ind, edge_nodes_grid, node_lon, node_lat
     ):
         """This function is only used to help generating the testcase and
         should not be used in the actual implementation. Construct an array to
@@ -691,7 +690,7 @@ def test_populate_bounds_near_pole(self):
             [[_xyz_to_lonlat_rad(*edge[0]), _xyz_to_lonlat_rad(*edge[1])] for edge in face_edges_cart])
 
         bounds = _populate_face_latlon_bound(face_edges_cart, face_edges_lonlat)
-        expected_bounds = np.array([[-1.20427718, -1.14935491], [0,0.13568803]])
+        expected_bounds = np.array([[-1.20427718, -1.14935491], [0, 0.13568803]])
         nt.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE)
 
     def test_populate_bounds_near_pole2(self):
@@ -703,13 +702,12 @@ def test_populate_bounds_near_pole2(self):
             [[4.06271283e-01, -4.78221112e-02, -9.12500241e-01], [3.57939780e-01, -4.88684203e-02, -9.32465008e-01]]
         ])
 
-
         # Apply the inverse transformation to get the lat lon coordinates
         face_edges_lonlat = np.array(
             [[_xyz_to_lonlat_rad(*edge[0]), _xyz_to_lonlat_rad(*edge[1])] for edge in face_edges_cart])
 
         bounds = _populate_face_latlon_bound(face_edges_cart, face_edges_lonlat)
-        expected_bounds = np.array([[-1.20427718, -1.14935491], [6.147497,4.960524e-16]])
+        expected_bounds = np.array([[-1.20427718, -1.14935491], [6.147497, 4.960524e-16]])
         nt.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE)
 
     def test_populate_bounds_long_face(self):
@@ -735,10 +733,7 @@ def test_populate_bounds_long_face(self):
         bounds = _populate_face_latlon_bound(face_edges_cart, face_edges_lonlat)
 
         # The expected bounds should not contains the south pole [0,-0.5*np.pi]
-        self.assertTrue(bounds[1][0] !=  0.0)
-
-
-
+        self.assertTrue(bounds[1][0] != 0.0)
 
     def test_populate_bounds_node_on_pole(self):
         # Generate a normal face that is crossing the antimeridian
@@ -828,7 +823,7 @@ def test_populate_bounds_pole_inside(self):
 class TestLatlonBoundsLatLonFace(TestCase):
 
     def _get_cartesian_face_edge_nodes_testcase_helper(
-            self,face_nodes_ind, face_edges_ind, edge_nodes_grid, node_x, node_y, node_z
+            self, face_nodes_ind, face_edges_ind, edge_nodes_grid, node_x, node_y, node_z
     ):
         """This function is only used to help generating the testcase and
         should not be used in the actual implementation. Construct an array to
@@ -887,7 +882,7 @@ def _get_cartesian_face_edge_nodes_testcase_helper(
         return cartesian_coordinates
 
     def _get_lonlat_rad_face_edge_nodes_testcase_helper(
-            self,face_nodes_ind, face_edges_ind, edge_nodes_grid, node_lon, node_lat
+            self, face_nodes_ind, face_edges_ind, edge_nodes_grid, node_lon, node_lat
     ):
         """This function is only used to help generating the testcase and
         should not be used in the actual implementation. Construct an array to
@@ -1088,7 +1083,7 @@ def test_populate_bounds_pole_inside(self):
 
 class TestLatlonBoundsGCAList(TestCase):
     def _get_cartesian_face_edge_nodes_testcase_helper(
-            self,face_nodes_ind, face_edges_ind, edge_nodes_grid, node_x, node_y, node_z
+            self, face_nodes_ind, face_edges_ind, edge_nodes_grid, node_x, node_y, node_z
     ):
         """This function is only used to help generating the testcase and
         should not be used in the actual implementation. Construct an array to
@@ -1147,7 +1142,7 @@ def _get_cartesian_face_edge_nodes_testcase_helper(
         return cartesian_coordinates
 
     def _get_lonlat_rad_face_edge_nodes_testcase_helper(
-            self,face_nodes_ind, face_edges_ind, edge_nodes_grid, node_lon, node_lat
+            self, face_nodes_ind, face_edges_ind, edge_nodes_grid, node_lon, node_lat
     ):
         """This function is only used to help generating the testcase and
         should not be used in the actual implementation. Construct an array to
@@ -1235,8 +1230,6 @@ def test_populate_bounds_normal(self):
                                              is_GCA_list=[True, False, True, False])
         nt.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE)
 
-
-
     def test_populate_bounds_antimeridian(self):
         # Generate a normal face that is crossing the antimeridian
         vertices_lonlat = [[350, 60.0], [350, 10.0], [50.0, 10.0], [50.0, 60.0]]
@@ -1436,7 +1429,6 @@ def test_face_bounds(self):
 
 class TestGeoDataFrame(TestCase):
 
-
     def test_engine(self):
         uxgrid = ux.open_grid(gridfile_geoflow)
         for engine in ['geopandas', 'spatialpandas']:
@@ -1483,3 +1475,20 @@ def test_cache_and_override(self):
         gdf_f = uxgrid.to_geodataframe(exclude_antimeridian=True)
 
         assert gdf_f is not gdf_e
+
+
+class TestStereographicProjection(TestCase):
+    def test_stereographic_projection(self):
+        lon = np.array(0)
+        lat = np.array(0)
+
+        central_lon = np.array(0)
+        central_lat = np.array(0)
+
+        x, y = stereographic_projection(lon, lat, central_lon, central_lat)
+
+        new_lon, new_lat = inverse_stereographic_projection(x, y, central_lon, central_lat)
+
+        self.assertTrue(lon == new_lon)
+        self.assertTrue(lat == new_lat)
+        self.assertTrue(x == y == 0)