diff --git a/docs/user-guide/from-points.ipynb b/docs/user-guide/from-points.ipynb index 9fa8fce59..62982f5fe 100644 --- a/docs/user-guide/from-points.ipynb +++ b/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, diff --git a/test/test_from_points.py b/test/test_from_points.py index 6cedef164..7bd8d53d6 100644 --- a/test/test_from_points.py +++ b/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) diff --git a/test/test_geometry.py b/test/test_geometry.py index 4d3a30f64..67e9369b9 100644 --- a/test/test_geometry.py +++ b/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) diff --git a/uxarray/grid/geometry.py b/uxarray/grid/geometry.py index 1ae95d99e..ed564aafb 100644 --- a/uxarray/grid/geometry.py +++ b/uxarray/grid/geometry.py @@ -1336,3 +1336,94 @@ def _construct_boundary_edge_indices(edge_face_connectivity): # If an edge only has one face saddling it than the mesh has holes in it edge_with_holes = np.where(edge_face_connectivity[:, 1] == INT_FILL_VALUE)[0] return edge_with_holes + + +def stereographic_projection(lon, lat, central_lon, central_lat): + """Projects a point on the surface of the sphere to a plane using stereographic projection + + Parameters + ---------- + lon: np.ndarray + Longitude coordinates of point + lat: np.ndarray + Latitude coordinate of point + central_lon: np.ndarray + Central longitude of projection + central_lat: np.ndarray + Central latitude of projection + Returns + ------- + x: np.ndarray + 2D x coordinate of projected point + y: np.ndarray + 2D y coordinate of projected point + """ + + # Convert to radians + lon = np.deg2rad(lon) + lat = np.deg2rad(lat) + central_lon = np.deg2rad(central_lon) + central_lat = np.deg2rad(central_lat) + + # Calculate constant used for calculation + k = 2.0 / ( + 1.0 + + np.sin(central_lat) * np.sin(lat) + + np.cos(central_lat) * np.cos(lat) * np.cos(lon - central_lon) + ) + + # Calculate the x and y coordinates + x = k * np.cos(lat) * np.sin(lon - central_lon) + y = k * ( + np.cos(central_lat) * np.sin(lat) + - np.sin(central_lat) * np.cos(lat) * np.cos(lon - central_lon) + ) + + return x, y + + +def inverse_stereographic_projection(x, y, central_lon, central_lat): + """Projects a point on a plane to the surface of the sphere using stereographic projection + + Parameters + ---------- + x: np.ndarray + 2D x coordinates of point + y: np.ndarray + 2D y coordinate of point + central_lon: np.ndarray + Central longitude of projection + central_lat: np.ndarray + Central latitude of projection + Returns + ------- + lon: np.ndarray + Longitude of projected point + lat: np.ndarray + Latitude of projected point + """ + + # If x and y are zero, the lon and lat will also be zero + + if x == 0 and y == 0: + return 0, 0 + + # Convert to radians + central_lat = np.deg2rad(central_lat) + + # Calculate constants used for calculation + p = np.sqrt(x**2 + y**2) + + c = 2 * np.arctan(p / 2) + + # Calculate the lon and lat of the coordinate + lon = central_lon + np.arctan2( + x * np.sin(c), + p * np.cos(central_lat) * np.cos(c) - y * np.sin(central_lat) * np.sin(c), + ) + + lat = np.arcsin( + np.cos(c) * np.sin(central_lat) + ((y * np.sin(c) * central_lat) / p) + ) + + return lon, lat diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py index e7007af26..88a755b0c 100644 --- a/uxarray/grid/grid.py +++ b/uxarray/grid/grid.py @@ -27,7 +27,10 @@ from uxarray.io._icon import _read_icon from uxarray.io._structured import _read_structured_grid from uxarray.io._voronoi import _spherical_voronoi_from_points -from uxarray.io._delaunay import _spherical_delaunay_from_points +from uxarray.io._delaunay import ( + _spherical_delaunay_from_points, + _regional_delaunay_from_points, +) from uxarray.formatting_html import grid_repr @@ -362,6 +365,7 @@ def from_points( - `'spherical_voronoi'`: Constructs a spherical Voronoi diagram. - `'spherical_delaunay'`: Constructs a spherical Delaunay triangulation. + - `'regional_delaunay'`: Constructs a regional Delaunay triangulation. Default is `'spherical_delaunay'`. @@ -387,6 +391,8 @@ def from_points( ds = _spherical_voronoi_from_points(_points, **kwargs) elif method == "spherical_delaunay": ds = _spherical_delaunay_from_points(_points, boundary_points) + elif method == "regional_delaunay": + ds = _regional_delaunay_from_points(_points, boundary_points) else: raise ValueError( f"Unsupported method '{method}'. Expected one of ['spherical_voronoi', 'spherical_delaunay']." diff --git a/uxarray/io/_delaunay.py b/uxarray/io/_delaunay.py index 5522793cc..7c57b5ae3 100644 --- a/uxarray/io/_delaunay.py +++ b/uxarray/io/_delaunay.py @@ -2,6 +2,8 @@ import xarray as xr from scipy.spatial import ConvexHull from uxarray.conventions import ugrid +from uxarray.grid.geometry import stereographic_projection +from scipy.spatial import Delaunay def _spherical_delaunay_from_points(points, boundary_points=None): @@ -61,3 +63,64 @@ def _spherical_delaunay_from_points(points, boundary_points=None): ) return out_ds + + +def _regional_delaunay_from_points(points, boundary_points=None): + """Generates a regional Delaunay triangulation from given points, + excluding triangles where all three nodes are boundary points.""" + out_ds = xr.Dataset() + + # Validate boundary_points if provided + if boundary_points is not None: + boundary_points = np.asarray(boundary_points) + if boundary_points.ndim != 1: + raise ValueError( + "boundary_points must be a 1D array-like of point indices." + ) + if np.any(boundary_points < 0) or np.any(boundary_points >= len(points)): + raise ValueError("boundary_points contain invalid indices.") + + node_x = points[:, 0] + node_y = points[:, 1] + node_z = points[:, 2] + + # Convert Cartesian coordinates to spherical coordinates (longitude and latitude in degrees) + node_lon_rad = np.arctan2(node_y, node_x) + node_lat_rad = np.arcsin(node_z / np.linalg.norm(points, axis=1)) + + node_lon = np.degrees(node_lon_rad) + node_lat = np.degrees(node_lat_rad) + + x_plane, y_plane = stereographic_projection(node_lon, node_lat, 0, 0) + + points = np.column_stack((x_plane, y_plane)) + + delaunay = Delaunay(points) + + # Obtain delaunay triangles + triangles = delaunay.simplices + + if boundary_points is not None: + # Create a boolean mask where True indicates triangles that do not have all nodes as boundary points + mask = ~np.isin(triangles, boundary_points).all(axis=1) + + # Apply the mask to filter out unwanted triangles + triangles = triangles[mask] + + # Assign the filtered triangles to the Dataset + out_ds["face_node_connectivity"] = xr.DataArray( + data=triangles, + dims=ugrid.FACE_NODE_CONNECTIVITY_DIMS, + attrs=ugrid.FACE_NODE_CONNECTIVITY_ATTRS, + ) + + # Assign node coordinates to the Dataset + out_ds["node_lon"] = xr.DataArray( + data=node_lon, dims=ugrid.NODE_DIM, attrs=ugrid.NODE_LON_ATTRS + ) + + out_ds["node_lat"] = xr.DataArray( + data=node_lat, dims=ugrid.NODE_DIM, attrs=ugrid.NODE_LAT_ATTRS + ) + + return out_ds