Finish all necessary changes

Author: hongyuchen1030

test/test_integrate.py

@@ -67,7 +67,7 @@ def test_get_faces_constLat_intersection_info_one_intersection():
     latitude_cart = -0.8660254037844386
     is_latlonface = False
     is_GCA_list = None
-    unique_intersections, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart,
+    unique_intersections,unique_intersections_lonlat, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart,
                                                                                          is_GCA_list, is_latlonface)
     assert len(unique_intersections) == 1
 
@@ -97,7 +97,7 @@ def test_get_faces_constLat_intersection_info_encompass_pole():
 
     is_latlonface = False
     is_GCA_list = None
-    unique_intersections, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart,
+    unique_intersections,unique_intersections_lonlat, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart,
                                                                                          is_GCA_list, is_latlonface)
     assert len(unique_intersections) <= 2 * len(face_edges_cart)
 
@@ -119,7 +119,7 @@ def test_get_faces_constLat_intersection_info_on_pole():
     latitude_cart = -0.9998476951563913
     is_latlonface = False
     is_GCA_list = None
-    unique_intersections, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart,
+    unique_intersections,unique_intersections_lonlat, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart,
                                                                                          is_GCA_list, is_latlonface)
     assert len(unique_intersections) == 2
 
@@ -141,7 +141,7 @@ def test_get_faces_constLat_intersection_info_near_pole():
     latitude_deg = np.rad2deg(latitude_rad)
     is_latlonface = False
     is_GCA_list = None
-    unique_intersections, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart,
+    unique_intersections,unique_intersections_lonlat, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart,
                                                                                          is_GCA_list, is_latlonface)
     assert len(unique_intersections) == 1
 

uxarray/core/zonal.py

@@ -5,7 +5,6 @@
 from uxarray.grid.integrate import _zonal_face_weights, _zonal_face_weights_robust
 
 
-
 def _compute_non_conservative_zonal_mean(uxda, latitudes, use_robust_weights=False):
     """Computes the non-conservative zonal mean across one or more latitudes."""
     uxgrid = uxda.uxgrid
@@ -50,4 +49,3 @@ def _compute_non_conservative_zonal_mean(uxda, latitudes, use_robust_weights=Fal
         )
 
     return result
-

uxarray/grid/geometry.py

@@ -103,7 +103,7 @@ def _unique_points(points, tolerance=ERROR_TOLERANCE):
 
 @njit(cache=True)
 def _pad_closed_face_nodes(
-        face_node_connectivity, n_face, n_max_face_nodes, n_nodes_per_face
+    face_node_connectivity, n_face, n_max_face_nodes, n_nodes_per_face
 ):
     """Pads a closed array of face nodes by inserting the first element at any
     point a fill value is encountered.
@@ -123,14 +123,14 @@ def _pad_closed_face_nodes(
 
 
 def _build_polygon_shells(
-        node_lon,
-        node_lat,
-        face_node_connectivity,
-        n_face,
-        n_max_face_nodes,
-        n_nodes_per_face,
-        projection=None,
-        central_longitude=0.0,
+    node_lon,
+    node_lat,
+    face_node_connectivity,
+    n_face,
+    n_max_face_nodes,
+    n_nodes_per_face,
+    projection=None,
+    central_longitude=0.0,
 ):
     """Builds an array of polygon shells, which can be used with Shapely to
     construct polygons."""
@@ -266,7 +266,7 @@ def _grid_to_polygon_geodataframe(grid, periodic_elements, projection, project,
 
 
 def _build_geodataframe_without_antimeridian(
-        polygon_shells, projected_polygon_shells, antimeridian_face_indices, engine
+    polygon_shells, projected_polygon_shells, antimeridian_face_indices, engine
 ):
     """Builds a ``spatialpandas.GeoDataFrame`` or
     ``geopandas.GeoDataFrame``excluding any faces that cross the
@@ -295,10 +295,10 @@ def _build_geodataframe_without_antimeridian(
 
 
 def _build_geodataframe_with_antimeridian(
-        polygon_shells,
-        projected_polygon_shells,
-        antimeridian_face_indices,
-        engine,
+    polygon_shells,
+    projected_polygon_shells,
+    antimeridian_face_indices,
+    engine,
 ):
     """Builds a ``spatialpandas.GeoDataFrame`` or ``geopandas.GeoDataFrame``
     including any faces that cross the antimeridian."""
@@ -317,9 +317,9 @@ def _build_geodataframe_with_antimeridian(
 
 
 def _build_corrected_shapely_polygons(
-        polygon_shells,
-        projected_polygon_shells,
-        antimeridian_face_indices,
+    polygon_shells,
+    projected_polygon_shells,
+    antimeridian_face_indices,
 ):
     if projected_polygon_shells is not None:
         # use projected shells if a projection is applied
@@ -438,7 +438,7 @@ def _build_corrected_polygon_shells(polygon_shells):
 
 
 def _grid_to_matplotlib_polycollection(
-        grid, periodic_elements, projection=None, **kwargs
+    grid, periodic_elements, projection=None, **kwargs
 ):
     """Constructs and returns a ``matplotlib.collections.PolyCollection``"""
 
@@ -642,7 +642,7 @@ def _get_polygons(grid, periodic_elements, projection=None, apply_projection=Tru
 
 
 def _grid_to_matplotlib_linecollection(
-        grid, periodic_elements, projection=None, **kwargs
+    grid, periodic_elements, projection=None, **kwargs
 ):
     """Constructs and returns a ``matplotlib.collections.LineCollection``"""
 
@@ -688,6 +688,7 @@ def _convert_shells_to_polygons(shells):
     return polygons
 
 
+# TODO: If we don't actually use this one, remove it
 def _pole_point_inside_polygon_cartesian(pole, face_edges_xyz):
     if isinstance(pole, str):
         pole = POLE_NAME_TO_INT[pole]
@@ -913,7 +914,7 @@ def _check_intersection(ref_edge_xyz, edges_xyz):
         for i in range(n_edges):
             edge_xyz = edges_xyz[i]
             if allclose(
-                    intersection_point, edge_xyz[0], atol=ERROR_TOLERANCE
+                intersection_point, edge_xyz[0], atol=ERROR_TOLERANCE
             ) or allclose(intersection_point, edge_xyz[1], atol=ERROR_TOLERANCE):
                 return 0
 
@@ -1140,10 +1141,10 @@ def insert_pt_in_latlonbox(old_box, new_pt, is_lon_periodic=True):
 
 @njit(cache=True)
 def _populate_face_latlon_bound(
-        face_edges_xyz,
-        face_edges_lonlat,
-        is_latlonface=False,
-        is_GCA_list=None,
+    face_edges_xyz,
+    face_edges_lonlat,
+    is_latlonface=False,
+    is_GCA_list=None,
 ):
     # Check if face_edges contains pole points
     has_north_pole = pole_point_inside_polygon(1, face_edges_xyz, face_edges_lonlat)
@@ -1198,9 +1199,9 @@ def _populate_face_latlon_bound(
             # Check if the node matches the pole point or if the pole point is within the edge
             max_abs_diff = np.max(np.abs(n1_cart - pole_point_xyz))
             if max_abs_diff <= ERROR_TOLERANCE or point_within_gca(
-                    pole_point_xyz,
-                    n1_cart,
-                    n2_cart,
+                pole_point_xyz,
+                n1_cart,
+                n2_cart,
             ):
                 is_center_pole = False
                 face_latlon_array = insert_pt_in_latlonbox(
@@ -1292,15 +1293,15 @@ def _populate_face_latlon_bound(
 
             # Insert extreme latitude points into the latlonbox
             if (
-                    abs(node1_lat_rad - lat_max) > ERROR_TOLERANCE
-                    and abs(node2_lat_rad - lat_max) > ERROR_TOLERANCE
+                abs(node1_lat_rad - lat_max) > ERROR_TOLERANCE
+                and abs(node2_lat_rad - lat_max) > ERROR_TOLERANCE
             ):
                 face_latlon_array = insert_pt_in_latlonbox(
                     face_latlon_array, np.array([lat_max, node1_lon_rad])
                 )
             elif (
-                    abs(node1_lat_rad - lat_min) > ERROR_TOLERANCE
-                    and abs(node2_lat_rad - lat_min) > ERROR_TOLERANCE
+                abs(node1_lat_rad - lat_min) > ERROR_TOLERANCE
+                and abs(node2_lat_rad - lat_min) > ERROR_TOLERANCE
             ):
                 face_latlon_array = insert_pt_in_latlonbox(
                     face_latlon_array, np.array([lat_min, node1_lon_rad])
@@ -1315,23 +1316,23 @@ def _populate_face_latlon_bound(
 
 @njit(cache=True, parallel=True)
 def compute_temp_latlon_array(
-        face_node_connectivity,
-        faces_edges_cartesian,
-        faces_edges_lonlat_rad,
-        n_nodes_per_face,
-        is_latlonface,
-        is_face_GCA_list,
-        INT_FILL_VALUE,
+    face_node_connectivity,
+    faces_edges_cartesian,
+    faces_edges_lonlat_rad,
+    n_nodes_per_face,
+    is_latlonface,
+    is_face_GCA_list,
+    INT_FILL_VALUE,
 ):
     n_face = face_node_connectivity.shape[0]
     temp_latlon_array = np.full((n_face, 2, 2), INT_FILL_VALUE, dtype=np.float64)
     for face_idx in prange(n_face):
         cur_face_edges_cartesian = faces_edges_cartesian[
-                                   face_idx, 0: n_nodes_per_face[face_idx]
-                                   ]
+            face_idx, 0 : n_nodes_per_face[face_idx]
+        ]
         cur_faces_edges_lonlat_rad = faces_edges_lonlat_rad[
-                                     face_idx, 0: n_nodes_per_face[face_idx]
-                                     ]
+            face_idx, 0 : n_nodes_per_face[face_idx]
+        ]
         if is_face_GCA_list is not None:
             is_GCA_list = is_face_GCA_list[face_idx]
         else:
@@ -1347,9 +1348,14 @@ def compute_temp_latlon_array(
 
 
 def _populate_faces_edges_cartesian(grid, return_array=False):
-    faces_edges_cartesian = _get_cartesian_faces_edge_nodes(grid.face_node_connectivity.values, grid.n_face,
-                                                            grid.n_max_face_edges, grid.node_x.values,
-                                                            grid.node_y.values, grid.node_z.values)
+    faces_edges_cartesian = _get_cartesian_faces_edge_nodes(
+        grid.face_node_connectivity.values,
+        grid.n_face,
+        grid.n_max_face_edges,
+        grid.node_x.values,
+        grid.node_y.values,
+        grid.node_z.values,
+    )
 
     faces_edges_cartesian_xarray = xr.DataArray(
         faces_edges_cartesian,
@@ -1358,7 +1364,7 @@ def _populate_faces_edges_cartesian(grid, return_array=False):
             "cf_role": "faces_edges_cartesian",
             "_FillValue": INT_FILL_VALUE,
             "long_name": "Provide the Cartesian coordinates of the edge for each face",
-            "start_index": INT_DTYPE(0)
+            "start_index": INT_DTYPE(0),
         },
     )
 
@@ -1369,9 +1375,13 @@ def _populate_faces_edges_cartesian(grid, return_array=False):
 
 
 def _populate_faces_edges_spherical(grid, return_array=False):
-    faces_edges_lonlat_rad = _get_lonlat_rad_faces_edge_nodes(grid.face_node_connectivity.values, grid.n_face,
-                                                              grid.n_max_face_edges, grid.node_lon.values,
-                                                              grid.node_lat.values)
+    faces_edges_lonlat_rad = _get_lonlat_rad_faces_edge_nodes(
+        grid.face_node_connectivity.values,
+        grid.n_face,
+        grid.n_max_face_edges,
+        grid.node_lon.values,
+        grid.node_lat.values,
+    )
 
     faces_edges_lonlat_rad_xarray = xr.DataArray(
         faces_edges_lonlat_rad,
@@ -1380,7 +1390,7 @@ def _populate_faces_edges_spherical(grid, return_array=False):
             "cf_role": "faces_edges_spherical",
             "_FillValue": INT_FILL_VALUE,
             "long_name": "Provide the Spherical coordinates (lon lat) in radians of the edge for each face",
-            "start_index": INT_DTYPE(0)
+            "start_index": INT_DTYPE(0),
         },
     )
 
@@ -1391,10 +1401,7 @@ def _populate_faces_edges_spherical(grid, return_array=False):
 
 
 def _populate_bounds(
-        grid,
-        is_latlonface: bool = False,
-        is_face_GCA_list=None,
-        return_array=False
+    grid, is_latlonface: bool = False, is_face_GCA_list=None, return_array=False
 ):
     """Populates the bounds of the grid based on the geometry of its faces,
     taking into account special conditions such as faces crossing the
@@ -1408,17 +1415,6 @@ def _populate_bounds(
         The grid object whose internal representation will be updated with the computed
         face bounds.
 
-    faces_edges_cartesian : np.ndarray, optional
-        An array of shape (n_face, n_max_face_edges, 2, 3) containing the Cartesian coordinates
-        of each face's edges. This array may include dummy values for grids with holes.
-        Default is None.
-
-    faces_edges_lonlat_rad : np.ndarray, optional
-        An array of shape (n_face, n_max_face_edges, 2, 2) containing the longitude and latitude
-        (in radians) coordinates of each face's edges. This array may include dummy values for
-        grids with holes.
-        Default is None.
-
     is_latlonface : bool, optional
         A global flag that indicates if faces are latlon faces. If True, all faces
         are treated as latlon faces, meaning that all edges are either longitude or
@@ -1563,16 +1559,16 @@ def stereographic_projection(lon, lat, central_lon, 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)
+        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)
+        np.cos(central_lat) * np.sin(lat)
+        - np.sin(central_lat) * np.cos(lat) * np.cos(lon - central_lon)
     )
 
     return x, y
@@ -1608,7 +1604,7 @@ def inverse_stereographic_projection(x, y, central_lon, central_lat):
     central_lat = np.deg2rad(central_lat)
 
     # Calculate constants used for calculation
-    p = np.sqrt(x ** 2 + y ** 2)
+    p = np.sqrt(x**2 + y**2)
 
     c = 2 * np.arctan(p / 2)
 
@@ -1627,9 +1623,9 @@ def inverse_stereographic_projection(x, y, central_lon, central_lat):
 
 @njit(cache=True)
 def point_in_face(
-        edges_xyz,
-        point_xyz,
-        inclusive=True,
+    edges_xyz,
+    point_xyz,
+    inclusive=True,
 ):
     """Determines if a point lies inside a face.
 
@@ -1679,9 +1675,9 @@ def point_in_face(
     for ind in range(len(edges_xyz)):
         # If the point lies on an edge, return True if inclusive
         if point_within_gca(
-                point_xyz,
-                edges_xyz[ind][0],
-                edges_xyz[ind][1],
+            point_xyz,
+            edges_xyz[ind][0],
+            edges_xyz[ind][1],
         ):
             if inclusive:
                 return True
@@ -1774,7 +1770,7 @@ def _populate_max_face_radius(self):
 
 @njit(cache=True)
 def calculate_max_face_radius(
-        face_node_connectivity, node_lats_rad, node_lons_rad, face_lats_rad, face_lons_rad
+    face_node_connectivity, node_lats_rad, node_lons_rad, face_lats_rad, face_lons_rad
 ):
     """Finds the max face radius in the mesh.
     Parameters
@@ -1848,7 +1844,7 @@ def haversine_distance(lon_a, lat_a, lon_b, lat_b):
 
     # Haversine formula
     equation_in_sqrt = (np.sin(dlat / 2) ** 2) + np.cos(lat_a) * np.cos(lat_b) * (
-            np.sin(dlon / 2) ** 2
+        np.sin(dlon / 2) ** 2
     )
     distance = 2 * np.arcsin(np.sqrt(equation_in_sqrt))