Initial Implementation of Latitude Longitude Face Bounds (#692)

* Initial Draft Algorithm * Revert "Initial Draft Algorithm" This reverts commit 99573b8. * boilerpalte * remove bounding box property for now * Work in Progress * Initial commit * Revert "Initial commit" This reverts commit 1ba6880. * Initial commit * Merge branch 'get-bounding-box' of https://github.com/UXARRAY/uxarray into get-bounding-box merge` * Fix complete * WIP * WIP * WIP * Fix interval looks * fix interval data storage type * fix interval data storage type * Initial draft complete * Support for the `is_latlonface` * Support for the `is_latlonface` and `is_GCA_List` * Finish final draft * Finish * refractor complete * Address PRs * add benchmark to face bounds * update bounds * rename bounds to face_bounds * fix tests * Delete grid_geoflow.exo * rename connectivity * Address PR Review * add warning * add warning --------- Co-authored-by: Hongyu Chen <hyvchen@ucdavis.edu> Co-authored-by: Orhan Eroglu <32553057+erogluorhan@users.noreply.github.com>
UXARRAY · Apr 11, 2024 · 0020265 · 0020265
1 parent b78516b
commit 0020265
Show file tree

Hide file tree

Showing 11 changed files with 1,035 additions and 15 deletions.
diff --git a/benchmarks/face_bounds.py b/benchmarks/face_bounds.py
@@ -0,0 +1,32 @@
+import os
+from pathlib import Path
+
+import uxarray as ux
+
+current_path = Path(os.path.dirname(os.path.realpath(__file__))).parents[0]
+
+grid_quad_hex = current_path / "test" / "meshfiles" / "ugrid" / "quad-hexagon" / "grid.nc"
+grid_geoflow = current_path / "test" / "meshfiles" / "ugrid" / "geoflow-small" / "grid.nc"
+grid_scrip = current_path / "test" / "meshfiles" / "scrip" / "outCSne8" / "outCSne8.nc"
+grid_mpas= current_path / "test" / "meshfiles" / "mpas" / "QU" / "oQU480.231010.nc"
+
+
+
+class FaceBounds:
+
+    params = [grid_quad_hex, grid_geoflow, grid_scrip, grid_mpas]
+
+
+    def setup(self, grid_path):
+        self.uxgrid = ux.open_grid(grid_path)
+
+    def teardown(self, n):
+        del self.uxgrid
+
+    def time_face_bounds(self, grid_path):
+        """Time to obtain ``Grid.face_bounds``"""
+        self.uxgrid.bounds
+
+    def peakmem_face_bounds(self, grid_path):
+        """Peak memory usage obtain ``Grid.face_bounds."""
+        face_bounds = self.uxgrid.bounds
diff --git a/docs/internal_api/index.rst b/docs/internal_api/index.rst
@@ -118,6 +118,7 @@ Connectivity
    grid.connectivity._populate_edge_face_connectivity
    grid.connectivity._populate_n_nodes_per_face
 
+
 Geometry
 --------
 .. autosummary::
@@ -138,6 +139,8 @@ Geometry
    grid.geometry._check_intersection
    grid.geometry._get_latlonbox_width
    grid.geometry._insert_pt_in_latlonbox
+   grid.geometry._populate_face_latlon_bound
+   grid.geometry._populate_bounds
 
 Coordinates
 -----------
@@ -172,6 +175,7 @@ Utils
    grid.utils._newton_raphson_solver_for_gca_constLat
    grid.utils._inv_jacobian
    grid.utils._get_cartesiain_face_edge_nodes
+   grid.utils._get_lonlat_rad_face_edge_nodes
 
 
 

diff --git a/docs/user_api/index.rst b/docs/user_api/index.rst
@@ -275,6 +275,7 @@ Grid Descriptors
 
    Grid.face_areas
    Grid.antimeridian_face_indices
+   Grid.bounds
 
 
 Attributes

diff --git a/test/test_geometry.py b/test/test_geometry.py
diff --git a/test/test_grid.py b/test/test_grid.py
@@ -8,11 +8,13 @@
 
 import uxarray as ux
 
-from uxarray.grid.connectivity import _populate_face_edge_connectivity, _build_edge_face_connectivity
+from uxarray.grid.connectivity import _populate_face_edge_connectivity, _build_edge_face_connectivity, _build_edge_node_connectivity
 
-from uxarray.grid.coordinates import _populate_lonlat_coord
+from uxarray.grid.coordinates import _populate_lonlat_coord, node_lonlat_rad_to_xyz
 
-from uxarray.constants import INT_FILL_VALUE
+from uxarray.constants import INT_FILL_VALUE, ERROR_TOLERANCE
+
+from uxarray.grid.arcs import extreme_gca_latitude
 
 try:
     import constants
@@ -392,7 +394,7 @@ def test_populate_cartesian_xyz_coord(self):
         # These points correspond to the eight vertices of a cube.
         lon_deg = [
             45.0001052295749, 45.0001052295749, 360 - 45.0001052295749,
-            360 - 45.0001052295749
+                                                360 - 45.0001052295749
         ]
         lat_deg = [
             35.2655522903022, -35.2655522903022, 35.2655522903022,
@@ -435,7 +437,7 @@ def test_populate_lonlat_coord(self):
 
         lon_deg = [
             45.0001052295749, 45.0001052295749, 360 - 45.0001052295749,
-            360 - 45.0001052295749
+                                                360 - 45.0001052295749
         ]
         lat_deg = [
             35.2655522903022, -35.2655522903022, 35.2655522903022,
@@ -589,7 +591,7 @@ def _revert_edges_conn_to_face_nodes_conn(
         for face_idx in range(len(face_edges_connectivity)):
             res_face_nodes_connectivity.append(face_nodes_dict[face_idx])
             while len(res_face_nodes_connectivity[face_idx]
-                     ) < original_face_nodes_connectivity.shape[1]:
+                      ) < original_face_nodes_connectivity.shape[1]:
                 res_face_nodes_connectivity[face_idx].append(ux.INT_FILL_VALUE)
 
         return np.array(res_face_nodes_connectivity)
@@ -856,7 +858,7 @@ def test_edge_face_connectivity_sample(self):
                 # shared edge
                 n_shared += 1
             elif edge_face[0] != INT_FILL_VALUE and edge_face[
-                    1] == INT_FILL_VALUE:
+                1] == INT_FILL_VALUE:
                 # edge borders one face
                 n_solo += 1
             else:
@@ -908,3 +910,25 @@ def test_from_face_vertices(self):
         uxgrid = ux.Grid.from_face_vertices(multi_face_latlon, latlon=True)
 
         single_face_cart = [(0.0,)]
+
+
+class TestLatlonBounds(TestCase):
+    def test_populate_bounds_GCA_mix(self):
+        face_1 = [[10.0, 60.0], [10.0, 10.0], [50.0, 10.0], [50.0, 60.0]]
+        face_2 = [[350, 60.0], [350, 10.0], [50.0, 10.0], [50.0, 60.0]]
+        face_3 = [[210.0, 80.0], [350.0, 60.0], [10.0, 60.0], [30.0, 80.0]]
+        face_4 = [[200.0, 80.0], [350.0, 60.0], [10.0, 60.0], [40.0, 80.0]]
+
+        faces = [face_1, face_2, face_3, face_4]
+
+        # Hand calculated bounds for the above faces in radians
+        expected_bounds = [[[0.17453293, 1.07370494],[0.17453293, 0.87266463]],
+                           [[0.17453293, 1.10714872],[6.10865238, 0.87266463]],
+                           [[1.04719755, 1.57079633],[3.66519143, 0.52359878]],
+                           [[1.04719755,1.57079633],[0.,         6.28318531]]]
+
+
+        grid = ux.Grid.from_face_vertices(faces, latlon=True)
+        bounds_xarray = grid.bounds
+        face_bounds = bounds_xarray.values
+        nt.assert_allclose(grid.bounds.values, expected_bounds, atol=ERROR_TOLERANCE)
diff --git a/test/test_helpers.py b/test/test_helpers.py
@@ -14,7 +14,7 @@
 
 from uxarray.grid.coordinates import node_lonlat_rad_to_xyz
 from uxarray.grid.arcs import point_within_gca, _angle_of_2_vectors, in_between
-from uxarray.grid.utils import _get_cartesian_face_edge_nodes
+from uxarray.grid.utils import _get_cartesian_face_edge_nodes, _get_lonlat_rad_face_edge_nodes
 from uxarray.grid.geometry import _pole_point_inside_polygon
 
 try:
@@ -319,5 +319,35 @@ def test_get_cartesian_face_edge_nodes_filled_value(self):
             'North', face_edges_connectivity_cartesian)
 
         # Assert that the result is True
+        self.assertTrue(result)
+
+    def test_get_lonlat_face_edge_nodes_pipeline(self):
+        # Create the vertices for the grid, based around the North Pole
+
+        vertices = [[0.5, 0.5, 0.5], [-0.5, 0.5, 0.5], [-0.5, -0.5, 0.5],
+                    [0.5, -0.5, 0.5]]
+
+        #Normalize the vertices
+        vertices = [x / np.linalg.norm(x) for x in vertices]
+
+        # Construct the grid from the vertices
+        grid = ux.Grid.from_face_vertices(vertices, latlon=False)
+        face_edges_connectivity_lonlat = _get_lonlat_rad_face_edge_nodes(
+            grid.face_node_connectivity.values[0],
+            grid.face_edge_connectivity.values[0],
+            grid.edge_node_connectivity.values, grid.node_lon.values,
+            grid.node_lat.values)
+
+        # Convert all the values into cartesian coordinates
+        face_edges_connectivity_cartesian = []
+        for i in range(len(face_edges_connectivity_lonlat)):
+            edge = face_edges_connectivity_lonlat[i]
+            edge_cart =   [node_lonlat_rad_to_xyz(node) for node in edge]
+            face_edges_connectivity_cartesian.append(edge_cart)
 
+        # Check that the face_edges_connectivity_cartesian works as an input to _pole_point_inside_polygon
+        result = ux.grid.geometry._pole_point_inside_polygon(
+            'North', np.array(face_edges_connectivity_cartesian))
+
+        # Assert that the result is True
         self.assertTrue(result)
diff --git a/uxarray/constants.py b/uxarray/constants.py
@@ -4,7 +4,11 @@
 INT_DTYPE = np.intp
 INT_FILL_VALUE = np.iinfo(INT_DTYPE).min
 
-ERROR_TOLERANCE = 1.0e-15
+# According to Ogita, Takeshi & Rump, Siegfried & Oishi, Shin’ichi. (2005). Accurate Sum and Dot Product.
+#     SIAM J. Scientific Computing. 26. 1955-1988. 10.1137/030601818.
+# half of the working precision is already the most optimal value for the error tolerance,
+# more tailored values will be used in the future.
+ERROR_TOLERANCE = 1.0e-8
 
 ENABLE_JIT_CACHE = True
 ENABLE_JIT = True

diff --git a/uxarray/grid/arcs.py b/uxarray/grid/arcs.py
@@ -93,6 +93,9 @@ def point_within_gca(pt, gca_cart, is_directed=False):
     if np.isclose(
         abs(GCRv1_lonlat[0] - GCRv0_lonlat[0]), np.pi, rtol=0, atol=ERROR_TOLERANCE
     ):
+        # Special case, if the pt is on the pole point, then set its longitude to the GCRv0_lonlat[0]
+        if np.isclose(abs(pt_lonlat[1]), np.pi / 2, rtol=0, atol=ERROR_TOLERANCE):
+            pt_lonlat[0] = GCRv0_lonlat[0]
         if not np.isclose(
             GCRv0_lonlat[0], pt_lonlat[0], rtol=0, atol=ERROR_TOLERANCE
         ) and not np.isclose(