Zonal mean helpers Fix (#955)

* Initial PR Open

* Finish the change of `pt_within_gcr`

* Finish the change of `_pole_point_inside_polygon`

* Finish the change of `gca_constlat_intersection`

* fix CI fail

* fix facebounds fail

* fix `extreme_gca`

* pre-commit

* Update asv.conf.json

* Update asv.conf.json

* Update uxarray/grid/intersections.py

Co-authored-by: Philip Chmielowiec <[email protected]>

* enable `JIT_CACHE`

* replace using the numba decorate function

* pre-commit

* Update uxarray/grid/integrate.py

* Update uxarray/grid/intersections.py

* Update uxarray/grid/intersections.py

* Update uxarray/grid/integrate.py

* add cross to intersections.py

* add isclose to integrate.py

* fix all() failing with params

* remove local benchmark

* optimize xyz_to_latlon_rad_scalar

* cleanup exception handling

---------

Co-authored-by: Philip Chmielowiec <[email protected]>

Author: hongyuchen1030

benchmarks/asv.conf.json

@@ -94,7 +94,7 @@
         "setuptools_scm": [""],
         "xarray": [""],
         "netcdf4": [""],
-        "pip+pyfma": [""]
+        "pip+git+https://github.com/philipc2/pyfma.git": [""]
     },
 
 

benchmarks/face_bounds.py

@@ -28,3 +28,13 @@ def time_face_bounds(self, grid_path):
     def peakmem_face_bounds(self, grid_path):
         """Peak memory usage obtain ``Grid.face_bounds."""
         face_bounds = self.uxgrid.bounds
+
+class Bounds:
+    def setup(self):
+        self.uxgrid = ux.open_grid(r"C:\Users\chmie\PycharmProjects\ncar-uxarray\uxarray-hongyu\benchmarks\oQU480.grid.nc")
+
+    def teardown(self):
+        del self.uxgrid
+
+    def time_bounds(self):
+        self.uxgrid.bounds

test/test_arcs.py

@@ -30,13 +30,6 @@ class TestIntersectionPoint(TestCase):
 
     def test_pt_within_gcr(self):
         # The GCR that's eexactly 180 degrees will have Value Error raised
-        gcr_180degree_cart = [
-            _lonlat_rad_to_xyz(0.0, 0.0),
-            _lonlat_rad_to_xyz(np.pi, 0.0)
-        ]
-        pt_same_lon_in = _lonlat_rad_to_xyz(0.0, 0.0)
-        with self.assertRaises(ValueError):
-            point_within_gca(np.asarray(pt_same_lon_in), np.asarray(gcr_180degree_cart))
 
         gcr_180degree_cart = [
             _lonlat_rad_to_xyz(0.0, np.pi / 2.0),
@@ -45,27 +38,27 @@ def test_pt_within_gcr(self):
 
         pt_same_lon_in = _lonlat_rad_to_xyz(0.0, 0.0)
         with self.assertRaises(ValueError):
-            point_within_gca(np.asarray(pt_same_lon_in), np.asarray(gcr_180degree_cart))
+            point_within_gca(pt_same_lon_in, gcr_180degree_cart)
 
         # Test when the point and the GCA all have the same longitude
         gcr_same_lon_cart = [
             _lonlat_rad_to_xyz(0.0, 1.5),
             _lonlat_rad_to_xyz(0.0, -1.5)
         ]
         pt_same_lon_in = _lonlat_rad_to_xyz(0.0, 0.0)
-        self.assertTrue(point_within_gca(np.asarray(pt_same_lon_in), np.asarray(gcr_same_lon_cart)))
+        self.assertTrue(point_within_gca(pt_same_lon_in, gcr_same_lon_cart))
 
         pt_same_lon_out = _lonlat_rad_to_xyz(0.0, 1.500000000000001)
-        res = point_within_gca(np.asarray(pt_same_lon_out), np.asarray(gcr_same_lon_cart))
+        res = point_within_gca(pt_same_lon_out, gcr_same_lon_cart)
         self.assertFalse(res)
 
         pt_same_lon_out_2 = _lonlat_rad_to_xyz(0.1, 1.0)
-        res = point_within_gca(np.asarray(pt_same_lon_out_2), np.asarray(gcr_same_lon_cart))
+        res = point_within_gca(pt_same_lon_out_2, gcr_same_lon_cart)
         self.assertFalse(res)
 
         # And if we increase the digital place by one, it should be true again
         pt_same_lon_out_add_one_place = _lonlat_rad_to_xyz(0.0, 1.5000000000000001)
-        res = point_within_gca(np.asarray(pt_same_lon_out_add_one_place), np.asarray(gcr_same_lon_cart))
+        res = point_within_gca(pt_same_lon_out_add_one_place, gcr_same_lon_cart)
         self.assertTrue(res)
 
         # Normal case
@@ -76,7 +69,7 @@ def test_pt_within_gcr(self):
                                                         -0.997]])
         pt_cart_within = np.array(
             [0.25616109352676675, 0.9246590335292105, -0.010021496695000144])
-        self.assertTrue(point_within_gca(np.asarray(pt_cart_within), np.asarray(gcr_cart_2), True))
+        self.assertTrue(point_within_gca(pt_cart_within, gcr_cart_2, True))
 
         # Test other more complicate cases : The anti-meridian case
 
@@ -87,16 +80,16 @@ def test_pt_within_gcr_antimeridian(self):
         gcr_cart = np.array([[0.351, -0.724, 0.593], [0.617, 0.672, 0.410]])
         pt_cart = np.array(
             [0.9438777657502077, 0.1193199333436068, 0.922714737029319])
-        self.assertTrue(point_within_gca(np.asarray(pt_cart), np.asarray(gcr_cart), is_directed=True))
+        self.assertTrue(point_within_gca(pt_cart, gcr_cart, is_directed=True))
         # If we swap the gcr, it should throw a value error since it's larger than 180 degree
         gcr_cart_flip = np.array([[0.617, 0.672, 0.410], [0.351, -0.724,
                                                           0.593]])
         with self.assertRaises(ValueError):
-            point_within_gca(np.asarray(pt_cart), np.asarray(gcr_cart_flip), is_directed=True)
+            point_within_gca(pt_cart, gcr_cart_flip, is_directed=True)
 
         # If we flip the gcr in the undirected mode, it should still work
         self.assertTrue(
-            point_within_gca(np.asarray(pt_cart), np.asarray(gcr_cart_flip), is_directed=False))
+            point_within_gca(pt_cart, gcr_cart_flip, is_directed=False))
 
         # 2nd anti-meridian case
         # GCR vertex0 in radian : [4.104711496596806, 0.5352983676533828],
@@ -107,9 +100,9 @@ def test_pt_within_gcr_antimeridian(self):
         pt_cart_within = np.array(
             [0.6136726305712109, 0.28442243941920053, -0.365605190899831])
         self.assertFalse(
-            point_within_gca(np.asarray(pt_cart_within), np.asarray(gcr_cart_1), is_directed=True))
+            point_within_gca(pt_cart_within, gcr_cart_1, is_directed=True))
         self.assertFalse(
-            point_within_gca(np.asarray(pt_cart_within), np.asarray(gcr_cart_1), is_directed=False))
+            point_within_gca(pt_cart_within, gcr_cart_1, is_directed=False))
 
         # The first case should not work and the second should work
         v1_rad = [0.1, 0.0]
@@ -119,10 +112,10 @@ def test_pt_within_gcr_antimeridian(self):
         gcr_cart = np.array([v1_cart, v2_cart])
         pt_cart = _lonlat_rad_to_xyz(0.01, 0.0)
         with self.assertRaises(ValueError):
-            point_within_gca(np.asarray(pt_cart), np.asarray(gcr_cart), is_directed=True)
+            point_within_gca(pt_cart, gcr_cart, is_directed=True)
         gcr_car_flipped = np.array([v2_cart, v1_cart])
         self.assertTrue(
-            point_within_gca(np.asarray(pt_cart), np.asarray(gcr_car_flipped), is_directed=True))
+            point_within_gca(pt_cart, gcr_car_flipped, is_directed=True))
 
     def test_pt_within_gcr_cross_pole(self):
         gcr_cart = np.array([[0.351, 0.0, 0.3], [-0.351, 0.0, 0.3]])

test/test_geometry.py

@@ -27,6 +27,13 @@
 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"
+
+
+# List of grid files to test
+grid_files_latlonBound = [grid_quad_hex, grid_geoflow, gridfile_CSne8, grid_mpas]
 
 class TestAntimeridian(TestCase):
 
@@ -56,6 +63,8 @@ def test_linecollection_execution(self):
         lines = uxgrid.to_linecollection()
 
 
+
+
 class TestPredicate(TestCase):
 
     def test_pole_point_inside_polygon_from_vertice_north(self):
@@ -366,6 +375,19 @@ def test_extreme_gca_latitude_max(self):
                                expected_max_latitude,
                                delta=ERROR_TOLERANCE)
 
+    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 ]])
+
+        # Calculate the maximum latitude
+        max_latitude = ux.grid.arcs.extreme_gca_latitude(gca_cart, 'max')
+
+        # Check if the maximum latitude is correct
+        expected_max_latitude = self._max_latitude_rad_iterative(gca_cart)
+        self.assertAlmostEqual(max_latitude,
+                               expected_max_latitude,
+                               delta=ERROR_TOLERANCE)
+
     def test_extreme_gca_latitude_min(self):
         # Define a great circle arc that is symmetrical around 0 degrees longitude
         gca_cart = np.array([
@@ -621,6 +643,103 @@ def test_populate_bounds_antimeridian(self):
         bounds = _populate_face_latlon_bound(face_edges_connectivity_cartesian, face_edges_connectivity_lonlat)
         nt.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE)
 
+    def test_populate_bounds_equator(self):
+        # the face is touching the equator
+        face_edges_cart = np.array([
+            [[0.99726469, -0.05226443, -0.05226443], [0.99862953, 0.0, -0.05233596]],
+            [[0.99862953, 0.0, -0.05233596], [1.0, 0.0, 0.0]],
+            [[1.0, 0.0, 0.0], [0.99862953, -0.05233596, 0.0]],
+            [[0.99862953, -0.05233596, 0.0], [0.99726469, -0.05226443, -0.05226443]]
+        ]
+        )
+        # 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([[-0.05235988, 0], [6.23082543, 0]])
+        nt.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE)
+
+    def test_populate_bounds_southSphere(self):
+        # The face is near the south pole but doesn't contains the pole
+        face_edges_cart = np.array([
+            [[-1.04386773e-01, -5.20500333e-02, -9.93173799e-01], [-1.04528463e-01, -1.28010448e-17, -9.94521895e-01]],
+            [[-1.04528463e-01, -1.28010448e-17, -9.94521895e-01], [-5.23359562e-02, -6.40930613e-18, -9.98629535e-01]],
+            [[-5.23359562e-02, -6.40930613e-18, -9.98629535e-01], [-5.22644277e-02, -5.22644277e-02, -9.97264689e-01]],
+            [[-5.22644277e-02, -5.22644277e-02, -9.97264689e-01], [-1.04386773e-01, -5.20500333e-02, -9.93173799e-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.51843645, -1.45388627], [3.14159265, 3.92699082]])
+        nt.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE)
+
+    def test_populate_bounds_near_pole(self):
+        # The face is near the south pole but doesn't contains the pole
+        face_edges_cart = np.array([
+            [[3.58367950e-01, 0.00000000e+00, -9.33580426e-01], [3.57939780e-01, 4.88684203e-02, -9.32465008e-01]],
+            [[3.57939780e-01, 4.88684203e-02, -9.32465008e-01], [4.06271283e-01, 4.78221112e-02, -9.12500241e-01]],
+            [[4.06271283e-01, 4.78221112e-02, -9.12500241e-01], [4.06736643e-01, 2.01762691e-16, -9.13545458e-01]],
+            [[4.06736643e-01, 2.01762691e-16, -9.13545458e-01], [3.58367950e-01, 0.00000000e+00, -9.33580426e-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], [0,0.13568803]])
+        nt.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE)
+
+    def test_populate_bounds_near_pole2(self):
+        # The face is near the south pole but doesn't contains the pole
+        face_edges_cart = np.array([
+            [[3.57939780e-01, -4.88684203e-02, -9.32465008e-01], [3.58367950e-01, 0.00000000e+00, -9.33580426e-01]],
+            [[3.58367950e-01, 0.00000000e+00, -9.33580426e-01], [4.06736643e-01, 2.01762691e-16, -9.13545458e-01]],
+            [[4.06736643e-01, 2.01762691e-16, -9.13545458e-01], [4.06271283e-01, -4.78221112e-02, -9.12500241e-01]],
+            [[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]])
+        nt.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE)
+
+    def test_populate_bounds_long_face(self):
+        """Test case where one of the face edges is a longitude GCA."""
+        face_edges_cart = np.array([
+            [[9.9999946355819702e-01, -6.7040475551038980e-04, 8.0396590055897832e-04],
+             [9.9999439716339111e-01, -3.2541253603994846e-03, -8.0110825365409255e-04]],
+            [[9.9999439716339111e-01, -3.2541253603994846e-03, -8.0110825365409255e-04],
+             [9.9998968839645386e-01, -3.1763643492013216e-03, -3.2474612817168236e-03]],
+            [[9.9998968839645386e-01, -3.1763643492013216e-03, -3.2474612817168236e-03],
+             [9.9998861551284790e-01, -8.2993711112067103e-04, -4.7004125081002712e-03]],
+            [[9.9998861551284790e-01, -8.2993711112067103e-04, -4.7004125081002712e-03],
+             [9.9999368190765381e-01, 1.7522916896268725e-03, -3.0944822356104851e-03]],
+            [[9.9999368190765381e-01, 1.7522916896268725e-03, -3.0944822356104851e-03],
+             [9.9999833106994629e-01, 1.6786820488050580e-03, -6.4892979571595788e-04]],
+            [[9.9999833106994629e-01, 1.6786820488050580e-03, -6.4892979571595788e-04],
+             [9.9999946355819702e-01, -6.7040475551038980e-04, 8.0396590055897832e-04]]
+        ])
+
+        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)
+
+        # The expected bounds should not contains the south pole [0,-0.5*np.pi]
+        self.assertTrue(bounds[1][0] !=  0.0)
+
+
+
+
     def test_populate_bounds_node_on_pole(self):
         # Generate a normal face that is crossing the antimeridian
         vertices_lonlat = [[10.0, 90.0], [10.0, 10.0], [50.0, 10.0], [50.0, 60.0]]
@@ -1116,6 +1235,8 @@ 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]]
@@ -1173,7 +1294,7 @@ def test_populate_bounds_node_on_pole(self):
         nt.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE)
 
     def test_populate_bounds_edge_over_pole(self):
-        # Generate a normal face that is crossing the antimeridian
+        # Generate a normal face that is around the north pole
         vertices_lonlat = [[210.0, 80.0], [350.0, 60.0], [10.0, 60.0], [30.0, 80.0]]
         vertices_lonlat = np.array(vertices_lonlat)
 
@@ -1288,6 +1409,31 @@ def test_populate_bounds_GCAList_mix(self):
             nt.assert_allclose(face_bounds[i], expected_bounds[i], atol=ERROR_TOLERANCE)
 
 
+# Test class
+class TestLatlonBoundsFiles:
+
+    def test_face_bounds(self):
+        """Test to ensure ``Grid.face_bounds`` works correctly for all grid
+        files."""
+        for grid_path in grid_files_latlonBound:
+            try:
+                # Open the grid file
+                self.uxgrid = ux.open_grid(grid_path)
+
+                # Test: Ensure the bounds are obtained
+                bounds = self.uxgrid.bounds
+                assert bounds is not None, f"Grid.face_bounds should not be None for {grid_path}"
+
+            except Exception as e:
+                # Print the failing grid file and re-raise the exception
+                print(f"Test failed for grid file: {grid_path}")
+                raise e
+
+            finally:
+                # Clean up the grid object
+                del self.uxgrid
+
+
 class TestGeoDataFrame(TestCase):
 
     def test_to_gdf(self):