Repository: UXARRAY/uxarray

Issue:
Implement `extreme_gca_lat` using strictly cartesian coordinates

Hints from issue comments:

Solution patch:
diff --git a/uxarray/grid/arcs.py b/uxarray/grid/arcs.py
index f01418986..dc53a2bd7 100644
--- a/uxarray/grid/arcs.py
+++ b/uxarray/grid/arcs.py
@@ -2,8 +2,6 @@
 import math
 
 
-from uxarray.grid.coordinates import _xyz_to_lonlat_rad_scalar
-
 from uxarray.grid.coordinates import (
     _normalize_xyz_scalar,
 )
@@ -184,18 +182,6 @@ def clip_scalar(a, a_min, a_max):
         return a
 
 
-def _extreme_gca_latitude_cartesian(gca_cart, extreme_type):
-    # should be shape [2, 2]
-    gca_lonlat = np.array(
-        [
-            _xyz_to_lonlat_rad_scalar(gca_cart[0, 0], gca_cart[0, 1], gca_cart[0, 2]),
-            _xyz_to_lonlat_rad_scalar(gca_cart[1, 0], gca_cart[1, 1], gca_cart[1, 2]),
-        ]
-    )
-
-    return extreme_gca_latitude(gca_cart, gca_lonlat, extreme_type)
-
-
 @njit(cache=True)
 def extreme_gca_latitude(gca_cart, gca_lonlat, extreme_type):
     """
@@ -277,3 +263,109 @@ def extreme_gca_latitude(gca_cart, gca_lonlat, extreme_type):
         return max(lat_n1, lat_n2)
     else:
         return min(lat_n1, lat_n2)
+
+
+@njit(cache=True)
+def extreme_gca_z(gca_cart, extreme_type):
+    """
+    Calculate the maximum or minimum latitude of a great circle arc defined
+    by two 3D points.
+
+    Parameters
+    ----------
+    gca_cart : numpy.ndarray
+        An array containing two 3D vectors that define a great circle arc.
+
+    extreme_type : str
+        The type of extreme latitude to calculate. Must be either 'max' or 'min'.
+
+    Returns
+    -------
+    float
+        The maximum or minimum z of the great circle arc
+
+    Raises
+    ------
+    ValueError
+        If `extreme_type` is not 'max' or 'min'.
+    """
+
+    # Validate extreme_type
+    if (extreme_type != "max") and (extreme_type != "min"):
+        raise ValueError("extreme_type must be either 'max' or 'min'")
+
+    # Extract the two points
+    n1 = gca_cart[0]
+    n2 = gca_cart[1]
+
+    # Compute dot product
+    dot_n1_n2 = dot(n1, n2)
+
+    # Compute denominator
+    denom = (n1[2] + n2[2]) * (dot_n1_n2 - 1.0)
+
+    # (z) coordinate
+    z_n1 = gca_cart[0][2]
+    z_n2 = gca_cart[1][2]
+
+    # Check if the denominator is zero to avoid division by zero
+    if denom != 0.0:
+        d_a_max = (n1[2] * dot_n1_n2 - n2[2]) / denom
+
+        # Handle cases where d_a_max is very close to 0 or 1
+        if isclose(d_a_max, 0.0, atol=ERROR_TOLERANCE) or isclose(
+            d_a_max, 1.0, atol=ERROR_TOLERANCE
+        ):
+            d_a_max = clip_scalar(d_a_max, 0.0, 1.0)
+
+        # Check if d_a_max is within the valid range
+        if (d_a_max > 0.0) and (d_a_max < 1.0):
+            # Compute the intermediate point on the GCA
+            node3 = (1.0 - d_a_max) * n1 + d_a_max * n2
+
+            # Normalize the intermediate point
+            x, y, z = _normalize_xyz_scalar(node3[0], node3[1], node3[2])
+            node3_normalized = np.empty(3)
+            node3_normalized[0] = x
+            node3_normalized[1] = y
+            node3_normalized[2] = z
+
+            d_z = clip_scalar(node3_normalized[2], -1.0, 1.0)
+
+            if extreme_type == "max":
+                return max3(d_z, z_n1, z_n2)
+            else:
+                return min3(d_z, z_n1, z_n2)
+
+    # If denom is zero or d_a_max is not in (0,1), return max or min of lat_n1 and lat_n2
+    if extreme_type == "max":
+        return max(z_n1, z_n2)
+    else:
+        return min(z_n1, z_n2)
+
+
+@njit(cache=True)
+def compute_arc_length(pt_a, pt_b):
+    """
+    Compute the great circle arc length between two points on a unit sphere at constant latitude.
+
+    Parameters
+    ----------
+    pt_a : tuple or array-like
+        First point coordinates (x, y, z) on unit sphere
+    pt_b : tuple or array-like
+        Second point coordinates (x, y, z) on unit sphere
+
+    Returns
+    -------
+    float
+        Arc length between the points at their constant latitude
+    """
+    x1, y1, z1 = pt_a
+    x2, y2, z2 = pt_b
+    rho = np.sqrt(1.0 - z1 * z2)
+    cross_2d = x1 * y2 - y1 * x2
+    dot_2d = x1 * x2 + y1 * y2
+    delta_theta = np.atan2(cross_2d, dot_2d)
+
+    return rho * abs(delta_theta)
diff --git a/uxarray/grid/integrate.py b/uxarray/grid/integrate.py
index cb2d9d1d2..d0beade57 100644
--- a/uxarray/grid/integrate.py
+++ b/uxarray/grid/integrate.py
@@ -1,100 +1,115 @@
 import numpy as np
+import polars as pl
+from numba import njit, prange, types
+from numba.typed import List
+
+from uxarray.grid.intersections import (
+    gca_const_lat_intersection,
+    get_number_of_intersections,
+)
 from uxarray.constants import ERROR_TOLERANCE, INT_FILL_VALUE
-from uxarray.grid.intersections import gca_const_lat_intersection
 from uxarray.grid.coordinates import _xyz_to_lonlat_rad
-import pandas as pd
-
 from uxarray.utils.computing import isclose
+from uxarray.grid.arcs import compute_arc_length
 
 DUMMY_EDGE_VALUE = [INT_FILL_VALUE, INT_FILL_VALUE, INT_FILL_VALUE]
 
+point_type = types.UniTuple(types.float64, 3)
+edge_type = types.UniTuple(types.int64, 2)
 
-def _get_zonal_faces_weight_at_constLat(
-    faces_edges_cart_candidate,
-    latitude_cart,
-    face_latlon_bound_candidate,
-    is_latlonface=False,
-    is_face_GCA_list=None,
-):
-    """Utilize the sweep line algorithm to calculate the weight of each face at
-    a constant latitude.
 
-     Parameters
+def _is_edge_gca(is_GCA_list, is_latlonface, edges_z):
+    """Determine if each edge is a Great Circle Arc (GCA) or a constant
+    latitude line in a vectorized manner.
+
+    Parameters:
     ----------
-    faces_edges_cart_candidate : np.ndarray
-        A list of the candidate face polygon represented by edges in Cartesian coordinates.
-        The faces must be selected in the previous step such that they will be intersected by the constant latitude.
-        It should have the same shape as the face_latlon_bound_candidate.
-        The input should not contain any 'INT_FILL_VALUE'. Shape: (n_faces(candidate), n_edges, 2, 3)
+    is_GCA_list : np.ndarray or None
+        An array indicating whether each edge is a GCA (True) or a constant latitude line (False).
+        Shape: (n_edges). If None, edge types are determined based on `is_latlonface` and the z-coordinates.
+    is_latlonface : bool
+        Flag indicating if all edges should be considered as lat-lon faces, which implies all edges
+        are either constant latitude or longitude lines.
+    edges_z : np.ndarray
+        Array containing the z-coordinates for each vertex of the edges. This is used to determine
+        whether edges are on the equator or if they are aligned in latitude when `is_GCA_list` is None.
+        Shape should be (n_edges, 2).
 
-    latitude_cart : float
-        The latitude in Cartesian coordinates (The normalized z coordinate)
+    Returns:
+    -------
+    np.ndarray
+        A boolean array where each element indicates whether the corresponding edge is considered a GCA.
+        True for GCA, False for constant latitude line.
+    """
+    if is_GCA_list is not None:
+        return is_GCA_list
+    if is_latlonface:
+        return ~np.isclose(edges_z[:, 0], edges_z[:, 1], atol=ERROR_TOLERANCE)
+    return ~(
+        np.isclose(edges_z[:, 0], 0, atol=ERROR_TOLERANCE)
+        & np.isclose(edges_z[:, 1], 0, atol=ERROR_TOLERANCE)
+    )
 
-    face_latlon_bound_candidate : np.ndarray
-        The list of latitude and longitude bounds of candidate faces.
-        It should have the same shape as the face_edges_cart_candidate.
-        Shape: (n_faces(candidate),,2, 2),
-        [...,[lat_min, lat_max], [lon_min, lon_max],...]
 
-    is_latlonface : bool, optional, default=False
-        A global flag to indicate if faces are latlon face. If True, then treat all faces as latlon faces. Latlon face means
-        That all edge is either a longitude or constant latitude line. If False, then all edges are GCA.
-         Default is False.
+def _zonal_face_weights_robust(
+    faces_edges_cart_candidate: np.ndarray,
+    latitude_cart: float,
+    face_latlon_bound_candidate: np.ndarray,
+    is_latlonface: bool = False,
+    is_face_GCA_list: np.ndarray | None = None,
+) -> pl.DataFrame:
+    """
+    Utilize the sweep line algorithm to calculate the weight of each face at
+    a constant latitude, returning a Polars DataFrame.
 
-    is_face_GCA_list : np.ndarray, optional (default=None)
-        A list of boolean values that indicates if the edge in that face is a GCA. Shape: (n_faces,n_edges).
-        True means edge face is a GCA.
-        False mean this edge is a constant latitude.
-        If None, all edges are considered as GCA. This attribute will overwrite the is_latlonface attribute.
+    Parameters
+    ----------
+    faces_edges_cart_candidate : np.ndarray
+        A list of the candidate face polygon represented by edges in Cartesian coordinates.
+        Shape: (n_faces(candidate), n_edges, 2, 3)
+    latitude_cart : float
+        The latitude in Cartesian coordinates (the normalized z coordinate).
+    face_latlon_bound_candidate : np.ndarray
+        An array with shape (n_faces, 2, 2), each face entry like [[lat_min, lat_max],[lon_min, lon_max]].
+    is_latlonface : bool, default=False
+        Global flag indicating if faces are lat-lon faces (edges are constant lat or long).
+    is_face_GCA_list : np.ndarray | None, default=None
+        Boolean array (n_faces, n_edges) indicating which edges are GCAs (True) or constant-lat (False).
+        If None, all edges are considered GCA.
 
     Returns
     -------
-    weights_df : pandas.DataFrame
-        A DataFrame with the calculated weights of each face. The DataFrame has two columns:
-        - 'face_index': The index of the face (integer).
-        - 'weight': The calculated weight of the face in radian (float).
-        The DataFrame is indexed by the face indices, providing a mapping from each face to its corresponding weight.
-
-    Notes
-    -----
-    Special handling is implemented for the cases when the latitude_cart is close to 1 or -1,
-    which corresponds to the poles (90 and -90 degrees). In these cases, if a pole point is
-    inside a face, that face's value is the only value that should be considered. If the pole
-    point is not inside any face, it lies on the boundary of surrounding faces, and their weights
-    are considered evenly since they only contain points rather than intervals.
-    This treatment is hard-coded in the function and should be tested with appropriate test cases.
+    weights_df : pl.DataFrame
+        DataFrame with columns ["face_index", "weight"], containing the per-face weights
+        (as a fraction of the total length of intersection).
     """
 
-    # Special case if the latitude_cart is 1 or -1, meaning right at the pole
-    # If the latitude_cart is close to 1 or -1 (indicating the pole), handle it separately.
-    # The -90 and 90 treatment is hard-coded in the function, based on:
-    # If a pole point is inside a face, then this face's value is the only value that should be considered.
-    # If the pole point is not inside any face, then it's on the boundary of faces around it, so their weights are even.
+    # Special case: latitude_cart close to +1 or -1 (near poles)
     if isclose(latitude_cart, 1, atol=ERROR_TOLERANCE) or isclose(
         latitude_cart, -1, atol=ERROR_TOLERANCE
     ):
-        # Now all candidate faces( the faces around the pole) are considered as the same weight
-        # If the face encompases the pole, then the weight is 1
-        weights = {
-            face_index: 1 / len(faces_edges_cart_candidate)
-            for face_index in range(len(faces_edges_cart_candidate))
-        }
-        # Convert weights to DataFrame
-        weights_df = pd.DataFrame(
-            list(weights.items()), columns=["face_index", "weight"]
+        # Evenly distribute weight among candidate faces
+        n_faces = len(faces_edges_cart_candidate)
+        weights = {face_index: 1.0 / n_faces for face_index in range(n_faces)}
+        # Convert dict to Polars DataFrame
+        return pl.DataFrame(
+            list(weights.items()), schema=["face_index", "weight"], orient="row"
         )
-        return weights_df
 
     intervals_list = []
 
-    # Iterate through all faces and their edges
+    # Iterate over faces
     for face_index, face_edges in enumerate(faces_edges_cart_candidate):
-        # Remove the Int_fill_value from the face_edges
+        # Remove edges that contain INT_FILL_VALUE
         face_edges = face_edges[np.all(face_edges != INT_FILL_VALUE, axis=(1, 2))]
+
+        # Which edges are GCA vs constant-lat?
         if is_face_GCA_list is not None:
             is_GCA_list = is_face_GCA_list[face_index]
         else:
             is_GCA_list = None
+
+        # Retrieve intervals for the current face
         face_interval_df = _get_zonal_face_interval(
             face_edges,
             latitude_cart,
@@ -102,83 +117,240 @@ def _get_zonal_faces_weight_at_constLat(
             is_latlonface=is_latlonface,
             is_GCA_list=is_GCA_list,
         )
-        # If any end of the interval is NaN
-        if face_interval_df.isnull().values.any():
-            # Skip this face as it is just being touched by the constant latitude
+
+        # Check if there are any null values in face_interval_df
+        has_null = face_interval_df.select(pl.col("*").is_null().any()).row(0)[0]
+        if has_null:
+            # Skip this face (only "touched" by the latitude)
             continue
-        # Check if the DataFrame is empty (start and end are both 0)
-        if (face_interval_df["start"] == 0).all() and (
-            face_interval_df["end"] == 0
-        ).all():
-            # Skip this face as it is just being touched by the constant latitude
+
+        # Check if all start == 0 and all end == 0
+        all_start_zero = face_interval_df.select((pl.col("start") == 0).all()).row(0)[0]
+        all_end_zero = face_interval_df.select((pl.col("end") == 0).all()).row(0)[0]
+        if all_start_zero and all_end_zero:
+            # Skip face being merely touched
             continue
-        else:
-            for _, row in face_interval_df.iterrows():
-                intervals_list.append(
-                    {"start": row["start"], "end": row["end"], "face_index": face_index}
-                )
 
-    intervals_df = pd.DataFrame(intervals_list)
+        # Add each interval row to intervals_list
+        for row in face_interval_df.iter_rows(named=True):
+            intervals_list.append(
+                {
+                    "start": row["start"],
+                    "end": row["end"],
+                    "face_index": face_index,
+                }
+            )
+
+    # Build a Polars DataFrame from intervals
+    intervals_df = pl.DataFrame(intervals_list)
+
+    # Process intervals to get overlap contributions
     try:
         overlap_contributions, total_length = _process_overlapped_intervals(
             intervals_df
         )
 
-        # Calculate weights for each face
-        weights = {
-            face_index: overlap_contributions.get(face_index, 0.0) / total_length
-            for face_index in range(len(faces_edges_cart_candidate))
-        }
+        # Build final weights dict
+        weights = {}
+        n_faces = len(faces_edges_cart_candidate)
+        for face_index in range(n_faces):
+            # fraction of total for this face
+            weights[face_index] = (
+                overlap_contributions.get(face_index, 0.0) / total_length
+            )
 
-        # Convert weights to DataFrame
-        weights_df = pd.DataFrame(
-            list(weights.items()), columns=["face_index", "weight"]
+        # Return as Polars DataFrame
+        weights_df = pl.DataFrame(
+            list(weights.items()), schema=["face_index", "weight"], orient="row"
         )
         return weights_df
 
     except ValueError:
-        print(f"Face index: {face_index}")
-        print(f"Face edges information: {face_edges}")
-        print(f"Constant z0: {latitude_cart}")
-        print(
-            f"Face latlon bound information: {face_latlon_bound_candidate[face_index]}"
-        )
-        print(f"Face interval information: {face_interval_df}")
-        # Handle the exception or propagate it further if necessary
+        # If an exception occurs, you can print debug info here if needed
         raise
 
 
-def _is_edge_gca(is_GCA_list, is_latlonface, edges_z):
-    """Determine if each edge is a Great Circle Arc (GCA) or a constant
-    latitude line in a vectorized manner.
+def _get_zonal_face_interval(
+    face_edges_cart: np.ndarray,
+    latitude_cart: float,
+    face_latlon_bound: np.ndarray,
+    is_latlonface: bool = False,
+    is_GCA_list: np.ndarray | None = None,
+) -> pl.DataFrame:
+    """
+    Processes a face polygon represented by edges in Cartesian coordinates
+    to find intervals where the face intersects with a given latitude. This
+    function handles directed and undirected Great Circle Arcs (GCAs) and edges
+    at constant latitude, returning a Polars DataFrame with columns ["start", "end"].
 
-    Parameters:
+    Requires the face edges to be sorted in counter-clockwise order, and the span of the
+    face in longitude should be less than pi. Also, all arcs/edges length should be within pi.
+
+    Users can specify which edges are GCAs and which are constant latitude using `is_GCA_list`.
+    However, edges on the equator are always treated as constant latitude edges regardless of
+    `is_GCA_list`.
+
+    Parameters
     ----------
-    is_GCA_list : np.ndarray or None
-        An array indicating whether each edge is a GCA (True) or a constant latitude line (False).
-        Shape: (n_edges). If None, edge types are determined based on `is_latlonface` and the z-coordinates.
-    is_latlonface : bool
-        Flag indicating if all edges should be considered as lat-lon faces, which implies all edges
-        are either constant latitude or longitude lines.
-    edges_z : np.ndarray
-        Array containing the z-coordinates for each vertex of the edges. This is used to determine
-        whether edges are on the equator or if they are aligned in latitude when `is_GCA_list` is None.
-        Shape should be (n_edges, 2).
+    face_edges_cart : np.ndarray
+        A face polygon represented by edges in Cartesian coordinates. Shape: (n_edges, 2, 3)
+    latitude_cart : float
+        The latitude in cartesian, the normalized Z coordinates.
+    face_latlon_bound : np.ndarray
+        The latitude and longitude bounds of the face. Shape: (2, 2), [[lat_min, lat_max], [lon_min, lon_max]]
+    is_latlonface : bool, optional, default=False
+        A global flag to indicate if faces are latlon face. If True, then treat all faces as latlon faces. Latlon face means
+        That all edge is either a longitude or constant latitude line. If False, then all edges are GCA.
+         Default is False. This attribute will overwrite the is_latlonface attribute.
+    is_GCA_list : np.ndarray, optional, default=False
+        An array indicating if each edge is a GCA (True) or a constant latitude (False). Shape: (n_edges,).
+        If None, all edges are considered as GCAs. Default is None.
 
-    Returns:
+    Returns
     -------
-    np.ndarray
-        A boolean array where each element indicates whether the corresponding edge is considered a GCA.
-        True for GCA, False for constant latitude line.
+    Intervals_df : pl.DataFrame
+        a Polars DataFrame with columns ["start", "end"]
     """
-    if is_GCA_list is not None:
-        return is_GCA_list
-    if is_latlonface:
-        return ~np.isclose(edges_z[:, 0], edges_z[:, 1], atol=ERROR_TOLERANCE)
-    return ~(
-        np.isclose(edges_z[:, 0], 0, atol=ERROR_TOLERANCE)
-        & np.isclose(edges_z[:, 1], 0, atol=ERROR_TOLERANCE)
-    )
+
+    face_lon_bound_left, face_lon_bound_right = face_latlon_bound[1]
+
+    try:
+        unique_intersections, pt_lon_min, pt_lon_max = (
+            _get_faces_constLat_intersection_info(
+                face_edges_cart, latitude_cart, is_GCA_list, is_latlonface
+            )
+        )
+
+        # If there's exactly one intersection, the face is only "touched"
+        if len(unique_intersections) == 1:
+            return pl.DataFrame({"start": [0.0], "end": [0.0]})
+
+        # Convert intersection points to (lon, lat) in radians
+        longitudes = np.array(
+            [_xyz_to_lonlat_rad(*pt.tolist())[0] for pt in unique_intersections]
+        )
+
+        # Handle special wrap-around cases (crossing anti-meridian, etc.)
+        if face_lon_bound_left >= face_lon_bound_right or (
+            face_lon_bound_left == 0 and face_lon_bound_right == 2 * np.pi
+        ):
+            if not (
+                (pt_lon_max >= np.pi and pt_lon_min >= np.pi)
+                or (0 <= pt_lon_max <= np.pi and 0 <= pt_lon_min <= np.pi)
+            ):
+                if pt_lon_max != 2 * np.pi and pt_lon_min != 0:
+                    # Add wrap-around points
+                    longitudes = np.append(longitudes, [0.0, 2 * np.pi])
+                elif pt_lon_max >= np.pi and pt_lon_min == 0:
+                    # If min is 0, but we really need 2*pi
+                    longitudes[longitudes == 0] = 2.0 * np.pi
+
+        # Sort unique longitudes
+        longitudes = np.unique(longitudes)
+        longitudes.sort()
+
+        # Pair sorted longitudes into intervals
+        starts = longitudes[::2]
+        ends = longitudes[1::2]
+
+        # Create Polars DataFrame
+        intervals_df = pl.DataFrame({"start": starts, "end": ends})
+        intervals_df_sorted = intervals_df.sort("start")
+
+        return intervals_df_sorted
+
+    except ValueError as e:
+        default_print_options = np.get_printoptions()
+        if str(e) == (
+            "No intersections are found for the face, please make sure the "
+            "build_latlon_box generates the correct results"
+        ):
+            np.set_printoptions(precision=16, suppress=False)
+            print(
+                "ValueError: No intersections are found for the face, make sure build_latlon_box is correct"
+            )
+            print(f"Face edges info:\n{face_edges_cart}")
+            print(f"Constant z_0: {latitude_cart}")
+            print(f"Face latlon bound:\n{face_latlon_bound}")
+            np.set_printoptions(**default_print_options)
+            raise
+        else:
+            np.set_printoptions(precision=17, suppress=False)
+            print(f"Face edges info:\n{face_edges_cart}")
+            print(f"Constant z_0: {latitude_cart}")
+            print(f"Face latlon bound:\n{face_latlon_bound}")
+            np.set_printoptions(**default_print_options)
+            raise
+
+
+def _process_overlapped_intervals(intervals_df: pl.DataFrame):
+    """Process the overlapped intervals using the sweep line algorithm.
+
+    This function processes multiple intervals per face using a sweep line algorithm,
+    calculating both individual face contributions and total length while handling
+    overlaps. The algorithm moves through sorted interval events (starts and ends),
+    maintaining a set of active faces and distributing overlap lengths equally.
+
+    Parameters
+    ----------
+    intervals_df : pl.DataFrame
+        A Polars DataFrame containing the intervals and corresponding face indices.
+        Required columns:
+            - start : numeric
+                Starting position of each interval
+            - end : numeric
+                Ending position of each interval
+            - face_index : int or str
+                Identifier for the face associated with each interval
+
+    Returns
+    -------
+    tuple[dict, float]
+        A tuple containing:
+        - dict: Maps face indices to their contributions to the total length,
+               where overlapping segments are weighted equally among active faces
+        - float: The total length of all intervals considering their overlaps
+    """
+
+    events = []
+    # Iterate Polars rows as dictionaries
+    for row in intervals_df.iter_rows(named=True):
+        events.append((row["start"], "start", row["face_index"]))
+        events.append((row["end"], "end", row["face_index"]))
+
+    # Sort the events by (position, event_type)
+    # so that 'start' comes before 'end' if position ties
+    events.sort(key=lambda x: (x[0], x[1]))
+
+    active_faces = set()
+    last_position = None
+    total_length = 0.0
+    overlap_contributions = {}
+
+    for position, event_type, face_idx in events:
+        if last_position is not None and active_faces:
+            segment_length = position - last_position
+            # Each face gets an equal share of this segment
+            segment_weight = segment_length / len(active_faces)
+            for active_face in active_faces:
+                overlap_contributions[active_face] = (
+                    overlap_contributions.get(active_face, 0.0) + segment_weight
+                )
+            total_length += segment_length
+
+        if event_type == "start":
+            active_faces.add(face_idx)
+        elif event_type == "end":
+            if face_idx in active_faces:
+                active_faces.remove(face_idx)
+            else:
+                raise ValueError(
+                    f"Error: Trying to remove face_idx {face_idx} not in active_faces"
+                )
+
+        last_position = position
+
+    return overlap_contributions, total_length
 
 
 def _get_faces_constLat_intersection_info(
@@ -237,10 +409,10 @@ def _get_faces_constLat_intersection_info(
         for idx, edge in enumerate(valid_edges):
             if is_GCA[idx]:
                 intersections = gca_const_lat_intersection(edge, latitude_cart)
-
-                if intersections.size == 0:
+                n_intersections = get_number_of_intersections(intersections)
+                if n_intersections == 0:
                     continue
-                elif intersections.shape[0] == 2:
+                elif n_intersections == 2:
                     intersections_pts_list_cart.extend(intersections)
                 else:
                     intersections_pts_list_cart.append(intersections[0])
@@ -249,7 +421,6 @@ def _get_faces_constLat_intersection_info(
     unique_intersections = np.unique(intersections_pts_list_cart, axis=0)
 
     if len(unique_intersections) == 2:
-        # TODO: vectorize?
         unique_intersection_lonlat = np.array(
             [_xyz_to_lonlat_rad(pt[0], pt[1], pt[2]) for pt in unique_intersections]
         )
@@ -285,212 +456,221 @@ def _get_faces_constLat_intersection_info(
         )
 
 
-def _get_zonal_face_interval(
-    face_edges_cart,
-    latitude_cart,
-    face_latlon_bound,
-    is_latlonface=False,
-    is_GCA_list=None,
-):
-    """Processes a face polygon represented by edges in Cartesian coordinates
-    to find intervals where the face intersects with a given latitude. This
-    function handles directed and undirected Great Circle Arcs (GCAs) and edges
-    at constant latitude.
+@njit(cache=True)
+def _add_edge(edges_list, i0, i1):
+    """Insert an edge into a list of edges in sorted order, ensuring no duplicates.
 
-    Requires the face edges to be sorted in counter-clockwise order, and the span of the
-    face in longitude should be less than pi. Also, all arcs/edges length should be within pi.
+    Parameters
+    ----------
+    edges_list : List
+        List of edge tuples (i0, i1) where each i represents a point index
+    i0 : int
+        First point index of the edge
+    i1 : int
+        Second point index of the edge
+    """
+    if i1 < i0:
+        i0, i1 = i1, i0
+
+    # Linear search for duplicates
+    for k in range(len(edges_list)):
+        e0, e1 = edges_list[k]
+        if e0 == i0 and e1 == i1:
+            return
+
+    edges_list.append((i0, i1))
 
-    Users can specify which edges are GCAs and which are constant latitude using `is_GCA_list`.
-    However, edges on the equator are always treated as constant latitude edges regardless of
-    `is_GCA_list`.
+
+@njit(cache=True)
+def _get_point_index(points_list, px, py, pz, tol=1e-12):
+    """
+    Find or create an index for a 3D point, checking for existing points within tolerance.
 
     Parameters
     ----------
-    face_edges_cart : np.ndarray
-        A face polygon represented by edges in Cartesian coordinates. Shape: (n_edges, 2, 3)
-    latitude_cart : float
-        The latitude in cartesian, the normalized Z coordinates.
-    face_latlon_bound : np.ndarray
-        The latitude and longitude bounds of the face. Shape: (2, 2), [[lat_min, lat_max], [lon_min, lon_max]]
-    is_latlonface : bool, optional, default=False
-        A global flag to indicate if faces are latlon face. If True, then treat all faces as latlon faces. Latlon face means
-        That all edge is either a longitude or constant latitude line. If False, then all edges are GCA.
-         Default is False. This attribute will overwrite the is_latlonface attribute.
-    is_GCA_list : np.ndarray, optional, default=False
-        An array indicating if each edge is a GCA (True) or a constant latitude (False). Shape: (n_edges,).
-        If None, all edges are considered as GCAs. Default is None.
+    points_list : List
+        List of point tuples (x, y, z)
+    px : float
+        X-coordinate of the point
+    py : float
+        Y-coordinate of the point
+    pz : float
+        Z-coordinate of the point
+    tol : float, optional
+        Tolerance for considering points as identical, default 1e-12
 
     Returns
     -------
-    Intervals_df : pd.DataFrame
-        A DataFrame containing the intervals stored as pandas Intervals for the face.
-        The columns of the DataFrame are: ['start', 'end']
+    int
+        Index of the matching point if found, or index of newly added point
     """
-    face_lon_bound_left, face_lon_bound_right = face_latlon_bound[1]
+    for i in range(len(points_list)):
+        (ex, ey, ez) = points_list[i]
+        if abs(ex - px) < tol and abs(ey - py) < tol and abs(ez - pz) < tol:
+            return i
 
-    # Call the vectorized function to process all edges
-    try:
-        # Call the vectorized function to process all edges
-        unique_intersections, pt_lon_min, pt_lon_max = (
-            _get_faces_constLat_intersection_info(
-                face_edges_cart, latitude_cart, is_GCA_list, is_latlonface
-            )
-        )
+    # If no match, add a new point
+    idx_new = len(points_list)
+    points_list.append((px, py, pz))
+    return idx_new
 
-        # Handle the special case where the unique_intersections is 1, which means the face is just being touched
-        if len(unique_intersections) == 1:
-            # If the face is just being touched, then just return the empty DataFrame
-            return pd.DataFrame({"start": [0.0], "end": [0.0]}, index=[0])
 
-        # Convert intersection points to longitude-latitude
-        longitudes = np.array(
-            [_xyz_to_lonlat_rad(*pt.tolist())[0] for pt in unique_intersections]
-        )
-
-        # Handle special wrap-around cases by checking the face bounds
-        if face_lon_bound_left >= face_lon_bound_right or (
-            face_lon_bound_left == 0 and face_lon_bound_right == 2 * np.pi
-        ):
-            if not (
-                (pt_lon_max >= np.pi and pt_lon_min >= np.pi)
-                or (0 <= pt_lon_max <= np.pi and 0 <= pt_lon_min <= np.pi)
-            ):
-                # If the anti-meridian is crossed, instead of just being touched,add the wrap-around points
-                if pt_lon_max != 2 * np.pi and pt_lon_min != 0:
-                    # They're at different sides of the 0-lon, adding wrap-around points
-                    longitudes = np.append(longitudes, [0.0, 2 * np.pi])
-                elif pt_lon_max >= np.pi and pt_lon_min == 0:
-                    # That means the face is actually from pt_lon_max to 2*pi.
-                    # Replace the 0 in longnitude with 2*pi
-                    longitudes[longitudes == 0] = 2 * np.pi
+@njit(cache=True)
+def _compute_face_arc_length(face_edges_xyz, z):
+    """
+    Compute the total arc length of a face's intersection with a line of constant latitude.
 
-        # Ensure longitudes are sorted
-        longitudes = np.unique(longitudes)
-        longitudes.sort()
+    Parameters
+    ----------
+    face_edges_xyz : np.ndarray
+        Array of shape (n_edges, 2, 3) containing the xyz coordinates of face edges
+    z : float
+        Z-coordinate of the constant latitude line
 
-        # Split the sorted longitudes into start and end points of intervals
-        starts = longitudes[::2]  # Start points
-        ends = longitudes[1::2]  # End points
+    Returns
+    -------
+    float
+        Total arc length of all intersections between the face and the constant latitude line
+    """
+    n_edges = face_edges_xyz.shape[0]
+
+    # 1) Typed lists for points and edges
+    points_list = List.empty_list(point_type)
+    edges_list = List.empty_list(edge_type)
+    singles_list = List.empty_list(types.int64)
+
+    # 2) Gather intersections from each edge
+    for e in range(n_edges):
+        edge = face_edges_xyz[e]  # shape (2,3)
+        intersections = gca_const_lat_intersection(edge, z)  # shape (2,3)
+        n_int = get_number_of_intersections(intersections)
+
+        if n_int == 1:
+            px0, py0, pz0 = intersections[0]
+            idx0 = _get_point_index(points_list, px0, py0, pz0)
+            singles_list.append(idx0)
+
+        elif n_int == 2:
+            px0, py0, pz0 = intersections[0]
+            px1, py1, pz1 = intersections[1]
+            idx0 = _get_point_index(points_list, px0, py0, pz0)
+            idx1 = _get_point_index(points_list, px1, py1, pz1)
+            _add_edge(edges_list, idx0, idx1)
+
+    # 3) Convert points_list to a (N,3) NumPy array
+    n_points = len(points_list)
+    points_array = np.empty((n_points, 3), dtype=np.float64)
+    for i in range(n_points):
+        (xx, yy, zz) = points_list[i]
+        points_array[i, 0] = xx
+        points_array[i, 1] = yy
+        points_array[i, 2] = zz
+
+    # 4) Finalize singles: connect each single to nearest neighbor
+    for s_idx in singles_list:
+        sx, sy, sz = points_array[s_idx]
+        best_i = -1
+        best_dist_sq = np.inf
+        for i in range(n_points):
+            if i == s_idx:
+                continue
+            dx = points_array[i, 0] - sx
+            dy = points_array[i, 1] - sy
+            dz = points_array[i, 2] - sz
+            dist_sq = dx * dx + dy * dy + dz * dz
+            if dist_sq < best_dist_sq:
+                best_dist_sq = dist_sq
+                best_i = i
+
+        # If we found a neighbor, add the edge in canonical form
+        if best_i >= 0 and best_i != s_idx:
+            _add_edge(edges_list, s_idx, best_i)
+
+    # 5) Sum arc lengths for all edges
+    total_length = 0.0
+    for k in range(len(edges_list)):
+        i0, i1 = edges_list[k]
+        total_length += compute_arc_length(points_array[i0], points_array[i1])
 
-        # Create the intervals DataFrame
-        Intervals_df = pd.DataFrame({"start": starts, "end": ends})
-        # For consistency, sort the intervals by start
-        interval_df_sorted = Intervals_df.sort_values(by="start").reset_index(drop=True)
-        return interval_df_sorted
+    return total_length
 
-    except ValueError as e:
-        default_print_options = np.get_printoptions()
-        if (
-            str(e)
-            == "No intersections are found for the face, please make sure the build_latlon_box generates the correct results"
-        ):
-            # Set print options for full precision
-            np.set_printoptions(precision=16, suppress=False)
 
-            print(
-                "ValueError: No intersections are found for the face, please make sure the build_latlon_box generates the correct results"
-            )
-            print(f"Face edges information: {face_edges_cart}")
-            print(f"Constant z_0: {latitude_cart}")
-            print(f"Face latlon bound information: {face_latlon_bound}")
+@njit(cache=True, parallel=True)
+def _zonal_face_weights_util_numba(
+    face_edges_xyz: np.ndarray,
+    n_edges_per_face: np.ndarray,
+    z: float,
+) -> np.ndarray:
+    """
+    Calculate normalized weights for faces intersecting a constant latitude using Numba
 
-            # Reset print options to default
-            np.set_printoptions(**default_print_options)
+    Parameters
+    ----------
+    face_edges_xyz : np.ndarray
+        Array of shape (n_face, max_edges, 2, 3) containing face edge coordinates
+    n_edges_per_face : np.ndarray
+        Array of shape (n_face,) containing the number of edges for each face
+    z : float
+        Z-coordinate of the constant latitude line
 
-            raise
-        else:
-            # Set print options for full precision
-            np.set_printoptions(precision=17, suppress=False)
+    Returns
+    -------
+    np.ndarray
+        Array of shape (n_face,) containing normalized weights for each face
+    """
+    n_face = face_edges_xyz.shape[0]
+    arc_lengths = np.zeros(n_face, dtype=np.float64)
 
-            print(f"Face edges information: {face_edges_cart}")
-            print(f"Constant z_0: {latitude_cart}")
-            print(f"Face latlon bound information: {face_latlon_bound}")
+    # 1) Pole Case: evenly distribute weights
+    if np.isclose(z, 1.0, atol=ERROR_TOLERANCE) or np.isclose(
+        z, -1.0, atol=ERROR_TOLERANCE
+    ):
+        return np.ones(n_face, dtype=np.float64) / n_face
 
-            # Reset print options to default
-            np.set_printoptions(**default_print_options)
+    # 2) Regular Case
+    for face_idx in prange(n_face):
+        n_edge = n_edges_per_face[face_idx]
+        face_data = face_edges_xyz[face_idx, :n_edge]  # shape (n_e, 2, 3)
+        arc_lengths[face_idx] = _compute_face_arc_length(face_data, z)
 
-            raise  # Re-raise the exception if it's not the expected ValueError
+    total_arc = np.sum(arc_lengths)
+    return arc_lengths / total_arc
 
 
-def _process_overlapped_intervals(intervals_df):
-    """Process the overlapped intervals using the sweep line algorithm,
-    considering multiple intervals per face.
+def _zonal_face_weights(
+    face_edges_xyz: np.ndarray,
+    face_bounds: np.ndarray,
+    n_edges_per_face: np.ndarray,
+    z: float,
+    check_equator: bool = False,
+) -> np.ndarray:
+    """
+    Calculate weights for faces intersecting a line of constant latitude, used for non-conservative zonal averaging.
 
     Parameters
     ----------
-    intervals_df : pd.DataFrame
-        The DataFrame that contains the intervals and the corresponding face index.
-        The columns of the DataFrame are: ['start', 'end', 'face_index']
+    face_edges_xyz : np.ndarray
+        Array of shape (n_face, max_edges, 2, 3) containing face edge coordinates
+    face_bounds : np.ndarray
+        Array containing bounds for each face
+    n_edges_per_face : np.ndarray
+        Array of shape (n_face,) containing the number of edges for each face
+    z : float
+        Z-coordinate of the constant latitude line
+    check_equator : bool
+        Whether to use a more precise weighting scheme near the equator
 
     Returns
     -------
-    dict
-        A dictionary where keys are face indices and values are their respective contributions to the total length.
-    float
-        The total length of all intervals considering their overlaps.
-
-    Example
-    -------
-    >>> intervals_data = [
-    ...     {"start": 0.0, "end": 100.0, "face_index": 0},
-    ...     {"start": 50.0, "end": 150.0, "face_index": 1},
-    ...     {"start": 140.0, "end": 150.0, "face_index": 2},
-    ... ]
-    >>> intervals_df = pd.DataFrame(intervals_data)
-    >>> overlap_contributions, total_length = _process_overlapped_intervals(
-    ...     intervals_df
-    ... )
-    >>> print(overlap_contributions)
-    >>> print(total_length)
+    np.ndarray
+        Array of weights for each face intersecting the latitude line
     """
-    events = []
-    for idx, row in intervals_df.iterrows():
-        events.append((row["start"], "start", row["face_index"]))
-        events.append((row["end"], "end", row["face_index"]))
-
-    events.sort(key=lambda x: (x[0], x[1]))
 
-    active_faces = set()
-    last_position = None
-    total_length = 0.0
-    overlap_contributions = {}
+    if check_equator:
+        # If near equator, use original approach
+        if np.isclose(z, 0.0, atol=ERROR_TOLERANCE):
+            overlap_result = _zonal_face_weights_robust(face_edges_xyz, z, face_bounds)
+            return overlap_result["weight"].to_numpy()
 
-    for position, event_type, face_idx in events:
-        if face_idx == 51:
-            pass
-        if last_position is not None and active_faces:
-            segment_length = position - last_position
-            segment_weight = segment_length / len(active_faces) if active_faces else 0
-            for active_face in active_faces:
-                overlap_contributions[active_face] = (
-                    overlap_contributions.get(active_face, 0) + segment_weight
-                )
-            total_length += segment_length
-
-        if event_type == "start":
-            # use try catch to handle the case where the face_idx is not be able to be added
-            try:
-                active_faces.add(face_idx)
-            except Exception as e:
-                print(f"An error occurred: {e}")
-                print(f"Face index: {face_idx}")
-                print(f"Position: {position}")
-                print(f"Event type: {event_type}")
-                print(f"Active faces: {active_faces}")
-                print(f"Last position: {last_position}")
-                print(f"Total length: {total_length}")
-                print(f"Overlap contributions: {overlap_contributions}")
-                print(f"Intervals data: {intervals_df}")
-                raise
-
-        elif event_type == "end":
-            if face_idx in active_faces:
-                active_faces.remove(face_idx)
-            else:
-                raise ValueError(
-                    f"Error: Trying to remove face_idx {face_idx} which is not in active_faces"
-                )
-
-        last_position = position
-
-    return overlap_contributions, total_length
+    # Otherwise, use the Numba approach
+    return _zonal_face_weights_util_numba(face_edges_xyz, n_edges_per_face, z)
diff --git a/uxarray/grid/intersections.py b/uxarray/grid/intersections.py
index cd9b702ed..0eabb7b26 100644
--- a/uxarray/grid/intersections.py
+++ b/uxarray/grid/intersections.py
@@ -1,17 +1,14 @@
 import numpy as np
 from uxarray.constants import MACHINE_EPSILON, ERROR_TOLERANCE, INT_DTYPE
 from uxarray.grid.utils import (
-    _newton_raphson_solver_for_gca_constLat,
     _angle_of_2_vectors,
 )
 from uxarray.grid.arcs import (
     in_between,
-    _extreme_gca_latitude_cartesian,
+    extreme_gca_z,
     point_within_gca,
 )
-import platform
-import warnings
-from uxarray.utils.computing import cross_fma, allclose, cross, norm
+from uxarray.utils.computing import allclose, cross, norm
 
 
 from numba import njit, prange
@@ -365,137 +362,93 @@ def gca_gca_intersection(gca_a_xyz, gca_b_xyz):
     return res[:count, :]
 
 
-def gca_const_lat_intersection(gca_cart, constZ, fma_disabled=True, verbose=False):
+@njit(cache=True)
+def gca_const_lat_intersection(gca_cart, const_z):
     """Calculate the intersection point(s) of a Great Circle Arc (GCA) and a
     constant latitude line in a Cartesian coordinate system.
 
-    To reduce relative errors, the Fused Multiply-Add (FMA) operation is utilized.
-    A warning is raised if the given coordinates are not in the cartesian coordinates, or
-    they cannot be accurately handled using floating-point arithmetic.
-
     Parameters
     ----------
     gca_cart : [2, 3] np.ndarray Cartesian coordinates of the two end points GCA.
-    constZ : float
+    const_z : float
         The constant latitude represented in cartesian of the latitude line.
-    fma_disabled : bool, optional (default=True)
-        If True, the FMA operation is disabled. And a naive `np.cross` is used instead.
-    verbose : bool, optional (default=False)
-        If True, the function prints out the intermediate results.
 
     Returns
     -------
     np.ndarray
-        Cartesian coordinates of the intersection point(s) the shape is [n_intersext_pts, 3].
-
-    Raises
-    ------
-    ValueError
-        If the input GCA is not in the cartesian [x, y, z] format.
+        Cartesian coordinates of the intersection point(s) the shape is [2, 3]. If no intersections are found,
+        all values a `nan`. If one intersection is found, the first column represent the intersection point, and
+        if two intersections are found, each column represents a point.
 
-    Warning
-    -------
-        If running on the Windows system with fma_disabled=False since the C/C++ implementation of FMA in MS Windows
-        is fundamentally broken. (bug report: https://bugs.python.org/msg312480)
-
-        If the intersection point cannot be converged using the Newton-Raphson method, the initial guess intersection
-        point is used instead, proceed with caution.
     """
+    res = np.empty((2, 3))
+    res.fill(np.nan)
+
     x1, x2 = gca_cart
 
     # Check if the constant latitude has the same latitude as the GCA endpoints
-    # We are using the relative tolerance and ERROR_TOLERANCE since the constZ is calculated from np.sin, which
-    # may have some floating-point error.
-    res = None
-    if np.isclose(x1[2], constZ, rtol=ERROR_TOLERANCE, atol=ERROR_TOLERANCE):
-        res = np.array([x1]) if res is None else np.vstack((res, x1))
+    x1_at_const_z = np.isclose(
+        x1[2], const_z, rtol=ERROR_TOLERANCE, atol=ERROR_TOLERANCE
+    )
+    x2_at_const_z = np.isclose(
+        x2[2], const_z, rtol=ERROR_TOLERANCE, atol=ERROR_TOLERANCE
+    )
+
+    if x1_at_const_z and x2_at_const_z:
+        res[0] = x1
+        res[1] = x2
+        return res
+    elif x1_at_const_z:
+        res[0] = x1
+        return res
+    elif x2_at_const_z:
+        res[0] = x2
+        return res
 
-    if np.isclose(x2[2], constZ, rtol=ERROR_TOLERANCE, atol=ERROR_TOLERANCE):
-        res = np.array([x2]) if res is None else np.vstack((res, x2))
+    # If the constant latitude is not the same as the GCA endpoints, calculate the intersection point
+    z_min = extreme_gca_z(gca_cart, extreme_type="min")
+    z_max = extreme_gca_z(gca_cart, extreme_type="max")
 
-    if res is not None:
+    # Check if the constant latitude is within the GCA range
+    if not in_between(z_min, const_z, z_max):
         return res
 
-    # If the constant latitude is not the same as the GCA endpoints, calculate the intersection point
-    lat_min = _extreme_gca_latitude_cartesian(gca_cart, extreme_type="min")
-    lat_max = _extreme_gca_latitude_cartesian(gca_cart, extreme_type="max")
+    n = cross(x1, x2)
 
-    constLat_rad = np.arcsin(constZ)
+    nx, ny, nz = n
 
-    # Check if the constant latitude is within the GCA range
-    # Because the constant latitude is calculated from np.sin, which may have some floating-point error,
-    if not in_between(lat_min, constLat_rad, lat_max):
-        pass
-        return np.array([])
+    s_tilde = np.sqrt(nx**2 + ny**2 - (nx**2 + ny**2 + nz**2) * const_z**2)
+    p1_x = -(1.0 / (nx**2 + ny**2)) * (const_z * nx * nz + s_tilde * ny)
+    p2_x = -(1.0 / (nx**2 + ny**2)) * (const_z * nx * nz - s_tilde * ny)
+    p1_y = -(1.0 / (nx**2 + ny**2)) * (const_z * ny * nz - s_tilde * nx)
+    p2_y = -(1.0 / (nx**2 + ny**2)) * (const_z * ny * nz + s_tilde * nx)
 
-    if fma_disabled:
-        n = cross(x1, x2)
+    p1 = np.array([p1_x, p1_y, const_z])
+    p2 = np.array([p2_x, p2_y, const_z])
 
-    else:
-        # Raise a warning for Windows users
-        if platform.system() == "Windows":
-            warnings.warn(
-                "The C/C++ implementation of FMA in MS Windows is reportedly broken. Use with care. (bug report: "
-                "https://bugs.python.org/msg312480)"
-                "The single rounding cannot be guaranteed, hence the relative error bound of 3u cannot be guaranteed."
-            )
-        n = cross_fma(x1, x2)
+    p1_intersects_gca = point_within_gca(p1, gca_cart[0], gca_cart[1])
+    p2_intersects_gca = point_within_gca(p2, gca_cart[0], gca_cart[1])
 
-    nx, ny, nz = n
+    if p1_intersects_gca and p2_intersects_gca:
+        res[0] = p1
+        res[1] = p2
+    elif p1_intersects_gca:
+        res[0] = p1
+    elif p2_intersects_gca:
+        res[0] = p2
 
-    s_tilde = np.sqrt(nx**2 + ny**2 - (nx**2 + ny**2 + nz**2) * constZ**2)
-    p1_x = -(1.0 / (nx**2 + ny**2)) * (constZ * nx * nz + s_tilde * ny)
-    p2_x = -(1.0 / (nx**2 + ny**2)) * (constZ * nx * nz - s_tilde * ny)
-    p1_y = -(1.0 / (nx**2 + ny**2)) * (constZ * ny * nz - s_tilde * nx)
-    p2_y = -(1.0 / (nx**2 + ny**2)) * (constZ * ny * nz + s_tilde * nx)
-
-    p1 = np.array([p1_x, p1_y, constZ])
-    p2 = np.array([p2_x, p2_y, constZ])
-
-    res = None
-
-    # Now test which intersection point is within the GCA range
-    if point_within_gca(p1, gca_cart[0], gca_cart[1]):
-        try:
-            converged_pt = _newton_raphson_solver_for_gca_constLat(
-                p1, gca_cart, verbose=verbose
-            )
-
-            if converged_pt is None:
-                # The point is not be able to be converged using the jacobi method, raise a warning and continue with p2
-                warnings.warn(
-                    "The intersection point cannot be converged using the Newton-Raphson method. "
-                    "The initial guess intersection point is used instead, procced with caution."
-                )
-                res = np.array([p1]) if res is None else np.vstack((res, p1))
-            else:
-                res = (
-                    np.array([converged_pt])
-                    if res is None
-                    else np.vstack((res, converged_pt))
-                )
-        except RuntimeError:
-            raise RuntimeError(f"Error encountered with initial guess: {p1}")
-
-    if point_within_gca(p2, gca_cart[0], gca_cart[1]):
-        try:
-            converged_pt = _newton_raphson_solver_for_gca_constLat(
-                p2, gca_cart, verbose=verbose
-            )
-            if converged_pt is None:
-                # The point is not be able to be converged using the jacobi method, raise a warning and continue with p2
-                warnings.warn(
-                    "The intersection point cannot be converged using the Newton-Raphson method. "
-                    "The initial guess intersection point is used instead, procced with caution."
-                )
-                res = np.array([p2]) if res is None else np.vstack((res, p2))
-            else:
-                res = (
-                    np.array([converged_pt])
-                    if res is None
-                    else np.vstack((res, converged_pt))
-                )
-        except RuntimeError:
-            raise RuntimeError(f"Error encountered with initial guess: {p2}")
-
-    return res if res is not None else np.array([])
+    return res
+
+
+@njit(cache=True)
+def get_number_of_intersections(arr):
+    """Returns the number of intersection points for the output of the gca-const-lat intersection."""
+    row1_is_nan = np.all(np.isnan(arr[0]))
+    row2_is_nan = np.all(np.isnan(arr[1]))
+
+    if row1_is_nan and row2_is_nan:
+        return 0
+    elif row2_is_nan:
+        return 1
+    else:
+        return 2
diff --git a/uxarray/grid/utils.py b/uxarray/grid/utils.py
index 228fd96fb..d9a6621c9 100644
--- a/uxarray/grid/utils.py
+++ b/uxarray/grid/utils.py
@@ -1,7 +1,5 @@
 import numpy as np
-from uxarray.constants import INT_FILL_VALUE, MACHINE_EPSILON
-import warnings
-import uxarray.utils.computing as ac_utils
+from uxarray.constants import INT_FILL_VALUE
 
 from numba import njit
 
@@ -80,142 +78,6 @@ def _angle_of_2_vectors(u, v):
         return 2 * np.pi - angle_u_v_rad
 
 
-def _inv_jacobian(x0, x1, y0, y1, z0, z1, x_i_old, y_i_old):
-    """Calculate the inverse Jacobian matrix for a given set of parameters.
-
-    Parameters
-    ----------
-    x0 : float
-        Description of x0.
-    x1 : float
-        Description of x1.
-    y0 : float
-        Description of y0.
-    y1 : float
-        Description of y1.
-    z0 : float
-        Description of z0.
-    z1 : float
-        Description of z1.
-    x_i_old : float
-        Description of x_i_old.
-    y_i_old : float
-        Description of y_i_old.
-
-    Returns
-    -------
-    numpy.ndarray or None
-        The inverse Jacobian matrix if it is non-singular, or None if a singular matrix is encountered.
-
-    Notes
-    -----
-    This function calculates the inverse Jacobian matrix based on the provided parameters. If the Jacobian matrix
-    is singular, a warning is printed, and None is returned.
-    """
-
-    # d_dx = (x0 * x_i_old - x1 * x_i_old * z0 + y0 * y_i_old * z1 - y1 * y_i_old * z0 - y1 * y_i_old * z0)
-    # d_dy = 2 * (x0 * x_i_old * z1 - x1 * x_i_old * z0 + y0 * y_i_old * z1 - y1 * y_i_old * z0)
-    #
-    # # row 1
-    # J[0, 0] = y_i_old / d_dx
-    # J[0, 1] = (x0 * z1 - z0 * x1) / d_dy
-    # # row 2
-    # J[1, 0] = x_i_old / d_dx
-    # J[1, 1] = (y0 * z1 - z0 * y1) / d_dy
-
-    # The Jacobian Matrix
-    jacobian = [
-        [ac_utils._fmms(y0, z1, z0, y1), ac_utils._fmms(x0, z1, z0, x1)],
-        [2 * x_i_old, 2 * y_i_old],
-    ]
-
-    # First check if the Jacobian matrix is singular
-    if np.linalg.matrix_rank(jacobian) < 2:
-        warnings.warn("The Jacobian matrix is singular.")
-        return None
-
-    try:
-        inverse_jacobian = np.linalg.inv(jacobian)
-    except np.linalg.LinAlgError as e:
-        # Print out the error message
-
-        cond_number = np.linalg.cond(jacobian)
-        print(f"Condition number: {cond_number}")
-        print(f"Jacobian matrix:\n{jacobian}")
-        print(f"An error occurred: {e}")
-        raise
-
-    return inverse_jacobian
-
-
-def _newton_raphson_solver_for_gca_constLat(
-    init_cart, gca_cart, max_iter=1000, verbose=False
-):
-    """Solve for the intersection point between a great circle arc and a
-    constant latitude.
-
-    Args:
-        init_cart (np.ndarray): Initial guess for the intersection point.
-        w0_cart (np.ndarray): First vector defining the great circle arc.
-        w1_cart (np.ndarray): Second vector defining the great circle arc.
-        max_iter (int, optional): Maximum number of iterations. Defaults to 1000.
-        verbose (bool, optional): Whether to print verbose output. Defaults to False.
-
-    Returns:
-        np.ndarray or None: The intersection point or None if the solver fails to converge.
-    """
-    tolerance = MACHINE_EPSILON * 100
-    w0_cart, w1_cart = gca_cart
-    error = float("inf")
-    constZ = init_cart[2]
-    y_guess = np.array(init_cart[0:2])
-    y_new = y_guess
-
-    _iter = 0
-
-    while error > tolerance and _iter < max_iter:
-        f_vector = np.array(
-            [
-                np.dot(
-                    np.cross(w0_cart, w1_cart),
-                    np.array([y_guess[0], y_guess[1], constZ]),
-                ),
-                y_guess[0] * y_guess[0]
-                + y_guess[1] * y_guess[1]
-                + constZ * constZ
-                - 1.0,
-            ]
-        )
-
-        try:
-            j_inv = _inv_jacobian(
-                w0_cart[0],
-                w1_cart[0],
-                w0_cart[1],
-                w1_cart[1],
-                w0_cart[2],
-                w1_cart[2],
-                y_guess[0],
-                y_guess[1],
-            )
-
-            if j_inv is None:
-                return None
-        except RuntimeError as e:
-            print(f"Encountered an error: {e}")
-            raise
-
-        y_new = y_guess - np.matmul(j_inv, f_vector)
-        error = np.max(np.abs(y_guess - y_new))
-        y_guess = y_new
-
-        if verbose:
-            print(f"Newton method iter: {_iter}, error: {error}")
-        _iter += 1
-
-    return np.append(y_new, constZ)
-
-
 def _swap_first_fill_value_with_last(arr):
     """Swap the first occurrence of INT_FILL_VALUE in each sub-array with the
     last value in the sub-array.


Test patch:
diff --git a/test/test_geometry.py b/test/test_geometry.py
index 7bf21cf6b..7de043479 100644
--- a/test/test_geometry.py
+++ b/test/test_geometry.py
@@ -10,7 +10,7 @@
 from uxarray.constants import ERROR_TOLERANCE, INT_FILL_VALUE
 import uxarray.utils.computing as ac_utils
 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, _extreme_gca_latitude_cartesian
+from uxarray.grid.arcs import extreme_gca_latitude, extreme_gca_z
 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, _pole_point_inside_polygon_cartesian, stereographic_projection, inverse_stereographic_projection
 
@@ -308,13 +308,13 @@ def test_extreme_gca_latitude_max():
         _normalize_xyz(*[-0.5, 0.5, 0.5])
     ])
 
-    max_latitude = _extreme_gca_latitude_cartesian(gca_cart, 'max')
-    expected_max_latitude = _max_latitude_rad_iterative(gca_cart)
+    max_latitude = extreme_gca_z(gca_cart, 'max')
+    expected_max_latitude = np.cos(_max_latitude_rad_iterative(gca_cart))
     assert np.isclose(max_latitude, expected_max_latitude, atol=ERROR_TOLERANCE)
 
     gca_cart = np.array([[0.0, 0.0, 1.0], [1.0, 0.0, 0.0]])
-    max_latitude = _extreme_gca_latitude_cartesian(gca_cart, 'max')
-    expected_max_latitude = np.pi / 2  # 90 degrees in radians
+    max_latitude = extreme_gca_z(gca_cart, 'max')
+    expected_max_latitude = 1.0
     assert np.isclose(max_latitude, expected_max_latitude, atol=ERROR_TOLERANCE)
 
 def test_extreme_gca_latitude_max_short():
@@ -322,13 +322,13 @@ def test_extreme_gca_latitude_max_short():
     gca_cart = np.array([[0.65465367, -0.37796447, -0.65465367], [0.6652466, -0.33896007, -0.6652466]])
 
     # Calculate the maximum latitude
-    max_latitude = _extreme_gca_latitude_cartesian(gca_cart, 'max')
+    max_latitude = np.asin(extreme_gca_z(gca_cart, 'max'))
 
     # Check if the maximum latitude is correct
     expected_max_latitude = _max_latitude_rad_iterative(gca_cart)
     assert np.isclose(max_latitude,
-                            expected_max_latitude,
-                            atol=ERROR_TOLERANCE)
+                      expected_max_latitude,
+                      atol=ERROR_TOLERANCE)
 
 
 def test_extreme_gca_latitude_min():
@@ -337,13 +337,13 @@ def test_extreme_gca_latitude_min():
         _normalize_xyz(*[-0.5, 0.5, -0.5])
     ])
 
-    min_latitude = _extreme_gca_latitude_cartesian(gca_cart, 'min')
+    min_latitude = np.asin(extreme_gca_z(gca_cart, 'min'))
     expected_min_latitude = _min_latitude_rad_iterative(gca_cart)
     assert np.isclose(min_latitude, expected_min_latitude, atol=ERROR_TOLERANCE)
 
     gca_cart = np.array([[0.0, 0.0, -1.0], [1.0, 0.0, 0.0]])
-    min_latitude = _extreme_gca_latitude_cartesian(gca_cart, 'min')
-    expected_min_latitude = -np.pi / 2  # 90 degrees in radians
+    min_latitude = np.asin(extreme_gca_z(gca_cart, 'min'))
+    expected_min_latitude = -np.pi / 2
     assert np.isclose(min_latitude, expected_min_latitude, atol=ERROR_TOLERANCE)
 
 def test_get_latlonbox_width():
@@ -440,9 +440,9 @@ def test_populate_bounds_normal_latlon_bounds_gca():
     vertices_rad = np.radians(vertices_lonlat)
     vertices_cart = np.vstack([_lonlat_rad_to_xyz(vertices_rad[:, 0], vertices_rad[:, 1])]).T
     lat_max = max(np.deg2rad(60.0),
-                  _extreme_gca_latitude_cartesian(np.array([vertices_cart[0], vertices_cart[3]]), extreme_type="max"))
+                  np.asin(extreme_gca_z(np.array([vertices_cart[0], vertices_cart[3]]), extreme_type="max")))
     lat_min = min(np.deg2rad(10.0),
-                  _extreme_gca_latitude_cartesian(np.array([vertices_cart[1], vertices_cart[2]]), extreme_type="min"))
+                  np.asin(extreme_gca_z(np.array([vertices_cart[1], vertices_cart[2]]), extreme_type="min")))
     lon_min = np.deg2rad(10.0)
     lon_max = np.deg2rad(50.0)
     grid = ux.Grid.from_face_vertices(vertices_lonlat, latlon=True)
@@ -467,9 +467,9 @@ def test_populate_bounds_antimeridian_latlon_bounds_gca():
     vertices_rad = np.radians(vertices_lonlat)
     vertices_cart = np.vstack([_lonlat_rad_to_xyz(vertices_rad[:, 0], vertices_rad[:, 1])]).T
     lat_max = max(np.deg2rad(60.0),
-                  _extreme_gca_latitude_cartesian(np.array([vertices_cart[0], vertices_cart[3]]), extreme_type="max"))
+                  np.asin(extreme_gca_z(np.array([vertices_cart[0], vertices_cart[3]]), extreme_type="max")))
     lat_min = min(np.deg2rad(10.0),
-                  _extreme_gca_latitude_cartesian(np.array([vertices_cart[1], vertices_cart[2]]), extreme_type="min"))
+                  np.asin(extreme_gca_z(np.array([vertices_cart[1], vertices_cart[2]]), extreme_type="min")))
     lon_min = np.deg2rad(350.0)
     lon_max = np.deg2rad(50.0)
     grid = ux.Grid.from_face_vertices(vertices_lonlat, latlon=True)
@@ -578,7 +578,7 @@ def test_populate_bounds_node_on_pole_latlon_bounds_gca():
     vertices_cart = np.vstack([_lonlat_rad_to_xyz(vertices_rad[:, 0], vertices_rad[:, 1])]).T
     lat_max = np.pi / 2
     lat_min = min(np.deg2rad(10.0),
-                  _extreme_gca_latitude_cartesian(np.array([vertices_cart[1], vertices_cart[2]]), extreme_type="min"))
+                  np.asin(extreme_gca_z(np.array([vertices_cart[1], vertices_cart[2]]), extreme_type="min")))
     lon_min = np.deg2rad(10.0)
     lon_max = np.deg2rad(50.0)
     grid = ux.Grid.from_face_vertices(vertices_lonlat, latlon=True)
@@ -604,7 +604,7 @@ def test_populate_bounds_edge_over_pole_latlon_bounds_gca():
     vertices_cart = np.vstack([_lonlat_rad_to_xyz(vertices_rad[:, 0], vertices_rad[:, 1])]).T
     lat_max = np.pi / 2
     lat_min = min(np.deg2rad(60.0),
-                  _extreme_gca_latitude_cartesian(np.array([vertices_cart[1], vertices_cart[2]]), extreme_type="min"))
+                  np.asin(extreme_gca_z(np.array([vertices_cart[1], vertices_cart[2]]), extreme_type="min")))
     lon_min = np.deg2rad(210.0)
     lon_max = np.deg2rad(30.0)
     grid = ux.Grid.from_face_vertices(vertices_lonlat, latlon=True)
@@ -630,7 +630,7 @@ def test_populate_bounds_pole_inside_latlon_bounds_gca():
     vertices_cart = np.vstack([_lonlat_rad_to_xyz(vertices_rad[:, 0], vertices_rad[:, 1])]).T
     lat_max = np.pi / 2
     lat_min = min(np.deg2rad(60.0),
-                  _extreme_gca_latitude_cartesian(np.array([vertices_cart[1], vertices_cart[2]]), extreme_type="min"))
+                  np.asin(extreme_gca_z(np.array([vertices_cart[1], vertices_cart[2]]), extreme_type="min")))
     lon_min = 0
     lon_max = 2 * np.pi
     grid = ux.Grid.from_face_vertices(vertices_lonlat, latlon=True)
diff --git a/test/test_integrate.py b/test/test_integrate.py
index f2566189f..ad9881c04 100644
--- a/test/test_integrate.py
+++ b/test/test_integrate.py
@@ -1,22 +1,32 @@
-import uxarray as ux
 import os
 from pathlib import Path
+import os
+from pathlib import Path
+
 import numpy as np
+import numpy.testing as nt
 import pandas as pd
+import polars as pl
+import pytest
+from polars.testing import assert_frame_equal
 
-import numpy.testing as nt
+from uxarray.constants import ERROR_TOLERANCE
 
 import uxarray as ux
-import pytest
-from uxarray.constants import ERROR_TOLERANCE, INT_FILL_VALUE
+from uxarray.constants import INT_FILL_VALUE
 from uxarray.grid.coordinates import _lonlat_rad_to_xyz
-from uxarray.grid.integrate import _get_zonal_face_interval, _process_overlapped_intervals, _get_zonal_faces_weight_at_constLat,_get_faces_constLat_intersection_info
+from uxarray.grid.integrate import _get_zonal_face_interval, _process_overlapped_intervals, \
+    _get_faces_constLat_intersection_info, _zonal_face_weights, \
+    _zonal_face_weights_robust
+
+from uxarray.grid.utils import _get_cartesian_face_edge_nodes
 
 current_path = Path(os.path.dirname(os.path.realpath(__file__)))
 
 gridfile_ne30 = current_path / "meshfiles" / "ugrid" / "outCSne30" / "outCSne30.ug"
 dsfile_var2_ne30 = current_path / "meshfiles" / "ugrid" / "outCSne30" / "outCSne30_var2.nc"
 
+
 def test_single_dim():
     """Integral with 1D data mapped to each face."""
     uxgrid = ux.open_grid(gridfile_ne30)
@@ -27,6 +37,7 @@ def test_single_dim():
     assert integral.ndim == len(dims) - 1
     nt.assert_almost_equal(integral, 4 * np.pi)
 
+
 def test_multi_dim():
     """Integral with 3D data mapped to each face."""
     uxgrid = ux.open_grid(gridfile_ne30)
@@ -37,6 +48,7 @@ def test_multi_dim():
     assert integral.ndim == len(dims) - 1
     nt.assert_almost_equal(integral, np.ones((5, 5)) * 4 * np.pi)
 
+
 def test_get_faces_constLat_intersection_info_one_intersection():
     face_edges_cart = np.array([
         [[-5.4411371445381629e-01, -4.3910468172333759e-02, -8.3786164521844386e-01],
@@ -55,9 +67,11 @@ 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, is_GCA_list, is_latlonface)
+    unique_intersections, 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
 
+
 def test_get_faces_constLat_intersection_info_encompass_pole():
     face_edges_cart = np.array([
         [[0.03982285692494229, 0.00351700770436231, 0.9992005658140627],
@@ -83,9 +97,11 @@ 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, is_GCA_list, is_latlonface)
+    unique_intersections, 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)
 
+
 def test_get_faces_constLat_intersection_info_on_pole():
     face_edges_cart = np.array([
         [[-5.2264427688714095e-02, -5.2264427688714102e-02, -9.9726468863423734e-01],
@@ -103,9 +119,11 @@ 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, is_GCA_list, is_latlonface)
+    unique_intersections, 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
 
+
 def test_get_faces_constLat_intersection_info_near_pole():
     face_edges_cart = np.array([
         [[-5.1693346290592648e-02, 1.5622531297347531e-01, -9.8636780641686628e-01],
@@ -123,9 +141,11 @@ 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, is_GCA_list, is_latlonface)
+    unique_intersections, 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
 
+
 def test_get_zonal_face_interval():
     """Test the _get_zonal_face_interval function for correct interval computation."""
     vertices_lonlat = [[1.6 * np.pi, 0.25 * np.pi],
@@ -154,6 +174,7 @@ def test_get_zonal_face_interval():
 
     nt.assert_array_almost_equal(actual_values_sorted, expected_values_sorted, decimal=13)
 
+
 def test_get_zonal_face_interval_empty_interval():
     """Test the _get_zonal_face_interval function for cases where the interval is empty."""
     face_edges_cart = np.array([
@@ -177,8 +198,9 @@ def test_get_zonal_face_interval_empty_interval():
     ])
 
     res = _get_zonal_face_interval(face_edges_cart, latitude_cart, face_latlon_bounds)
-    expected_res = pd.DataFrame({"start": [0.0], "end": [0.0]})
-    pd.testing.assert_frame_equal(res, expected_res)
+    expected_res = pl.DataFrame({"start": [0.0], "end": [0.0]})
+    assert_frame_equal(res, expected_res)
+
 
 def test_get_zonal_face_interval_encompass_pole():
     """Test the _get_zonal_face_interval function for cases where the face encompasses the pole inside."""
@@ -205,22 +227,23 @@ def test_get_zonal_face_interval_encompass_pole():
         [np.arcsin(0.9992005658145248), 0.5 * np.pi],
         [0, 2 * np.pi]
     ])
-    expected_df = pd.DataFrame({
+    expected_df = pl.DataFrame({
         'start': [0.000000, 1.101091, 2.357728, 3.614365, 4.871002, 6.127640],
         'end': [0.331721, 1.588358, 2.844995, 4.101632, 5.358270, 6.283185]
     })
 
     res = _get_zonal_face_interval(face_edges_cart, latitude_cart, face_latlon_bounds)
 
-    pd.testing.assert_frame_equal(res, expected_df)
+    assert_frame_equal(res, expected_df)
+
 
 def test_get_zonal_face_interval_FILL_VALUE():
     """Test the _get_zonal_face_interval function for cases where there are dummy nodes."""
     dummy_node = [INT_FILL_VALUE, INT_FILL_VALUE, INT_FILL_VALUE]
     vertices_lonlat = [[1.6 * np.pi, 0.25 * np.pi],
-                         [1.6 * np.pi, -0.25 * np.pi],
-                         [0.4 * np.pi, -0.25 * np.pi],
-                         [0.4 * np.pi, 0.25 * np.pi]]
+                       [1.6 * np.pi, -0.25 * np.pi],
+                       [0.4 * np.pi, -0.25 * np.pi],
+                       [0.4 * np.pi, 0.25 * np.pi]]
     vertices = [_lonlat_rad_to_xyz(*v) for v in vertices_lonlat]
 
     face_edge_nodes = np.array([[vertices[0], vertices[1]],
@@ -244,6 +267,7 @@ def test_get_zonal_face_interval_FILL_VALUE():
 
     nt.assert_array_almost_equal(actual_values_sorted, expected_values_sorted, decimal=13)
 
+
 def test_get_zonal_face_interval_GCA_constLat():
     vertices_lonlat = [[-0.4 * np.pi, 0.25 * np.pi],
                        [-0.4 * np.pi, -0.25 * np.pi],
@@ -273,11 +297,12 @@ def test_get_zonal_face_interval_GCA_constLat():
 
     nt.assert_array_almost_equal(actual_values_sorted, expected_values_sorted, decimal=13)
 
+
 def test_get_zonal_face_interval_equator():
     """Test that the face interval is correctly computed when the latitude
             is at the equator."""
     vertices_lonlat = [[-0.4 * np.pi, 0.25 * np.pi], [-0.4 * np.pi, 0.0],
-                        [0.4 * np.pi, 0.0], [0.4 * np.pi, 0.25 * np.pi]]
+                       [0.4 * np.pi, 0.0], [0.4 * np.pi, 0.25 * np.pi]]
 
     vertices = [_lonlat_rad_to_xyz(*v) for v in vertices_lonlat]
 
@@ -287,9 +312,9 @@ def test_get_zonal_face_interval_equator():
                                 [vertices[3], vertices[0]]])
 
     interval_df = _get_zonal_face_interval(face_edge_nodes, 0.0,
-                                            np.array([[-0.25 * np.pi, 0.25 * np.pi], [1.6 * np.pi,
-                                                                                        0.4 * np.pi]]),
-                                            is_GCA_list=np.array([True, True, True, True]))
+                                           np.array([[-0.25 * np.pi, 0.25 * np.pi], [1.6 * np.pi,
+                                                                                     0.4 * np.pi]]),
+                                           is_GCA_list=np.array([True, True, True, True]))
     expected_interval_df = pd.DataFrame({
         'start': [1.6 * np.pi, 0.0],
         'end': [2.0 * np.pi, 00.4 * np.pi]
@@ -306,9 +331,9 @@ def test_get_zonal_face_interval_equator():
 
     # Even if we change the is_GCA_list to False, the result should be the same
     interval_df = _get_zonal_face_interval(face_edge_nodes, 0.0,
-                                            np.array([[-0.25 * np.pi, 0.25 * np.pi], [1.6 * np.pi,
-                                                                                        0.4 * np.pi]]),
-                                            is_GCA_list=np.array([True, False, True, False]))
+                                           np.array([[-0.25 * np.pi, 0.25 * np.pi], [1.6 * np.pi,
+                                                                                     0.4 * np.pi]]),
+                                           is_GCA_list=np.array([True, False, True, False]))
     expected_interval_df = pd.DataFrame({
         'start': [1.6 * np.pi, 0.0],
         'end': [2.0 * np.pi, 00.4 * np.pi]
@@ -324,6 +349,7 @@ def test_get_zonal_face_interval_equator():
     # Asserting almost equal arrays
     nt.assert_array_almost_equal(actual_values_sorted, expected_values_sorted, decimal=10)
 
+
 def test_process_overlapped_intervals_overlap_and_gap():
     intervals_data = [
         {
@@ -353,22 +379,43 @@ def test_process_overlapped_intervals_overlap_and_gap():
         },
     ]
 
-    df = pd.DataFrame(intervals_data)
-    df['interval'] = df.apply(lambda row: pd.Interval(
-        left=row['start'], right=row['end'], closed='both'),
-                              axis=1)
-    df['interval'] = df['interval'].astype('interval[float64]')
-
-    expected_overlap_contributions = np.array({
+    # Create Polars DataFrame with explicit types
+    df = pl.DataFrame(
+        {
+            'start': pl.Series([x['start'] for x in intervals_data], dtype=pl.Float64),
+            'end': pl.Series([x['end'] for x in intervals_data], dtype=pl.Float64),
+            'face_index': pl.Series([x['face_index'] for x in intervals_data], dtype=pl.Int64)
+        }
+    )
+
+    # Expected results
+    expected_overlap_contributions = {
         0: 75.0,
         1: 70.0,
         2: 5.0,
         3: 100.0,
         4: 90.0
-    })
+    }
+
+    # Process intervals
     overlap_contributions, total_length = _process_overlapped_intervals(df)
-    assert total_length == 340.0
-    nt.assert_array_equal(overlap_contributions, expected_overlap_contributions)
+
+    # Assertions
+    assert abs(total_length - 340.0) < 1e-10  # Using small epsilon for float comparison
+
+    # Check each contribution matches expected value
+    for face_idx, expected_value in expected_overlap_contributions.items():
+        assert abs(overlap_contributions[face_idx] - expected_value) < 1e-10, \
+            f"Mismatch for face_index {face_idx}: expected {expected_value}, got {overlap_contributions[face_idx]}"
+
+    # Check that we have all expected face indices
+    assert set(overlap_contributions.keys()) == set(expected_overlap_contributions.keys()), \
+        "Mismatch in face indices"
+
+    # Check total contributions sum matches total length
+    assert abs(sum(overlap_contributions.values()) - total_length) < 1e-10, \
+        "Sum of contributions doesn't match total length"
+
 
 def test_process_overlapped_intervals_antimeridian():
     intervals_data = [
@@ -399,94 +446,126 @@ def test_process_overlapped_intervals_antimeridian():
         },
     ]
 
-    df = pd.DataFrame(intervals_data)
-    df['interval'] = df.apply(lambda row: pd.Interval(
-        left=row['start'], right=row['end'], closed='both'),
-                              axis=1)
-    df['interval'] = df['interval'].astype('interval[float64]')
-
-    expected_overlap_contributions = np.array({
+    # Create Polars DataFrame with explicit types
+    df = pl.DataFrame(
+        {
+            'start': pl.Series([x['start'] for x in intervals_data], dtype=pl.Float64),
+            'end': pl.Series([x['end'] for x in intervals_data], dtype=pl.Float64),
+            'face_index': pl.Series([x['face_index'] for x in intervals_data], dtype=pl.Int64)
+        }
+    )
+
+    # Expected results for antimeridian case
+    expected_overlap_contributions = {
         0: 105.0,
         1: 25.0,
         2: 175.0,
         3: 45.0
-    })
+    }
+
+    # Process intervals
     overlap_contributions, total_length = _process_overlapped_intervals(df)
-    assert total_length == 350.0
-    nt.assert_array_equal(overlap_contributions, expected_overlap_contributions)
+
+    # Assert total length
+    assert abs(total_length - 350.0) < 1e-10, \
+        f"Expected total length 350.0, got {total_length}"
+
+    # Check each contribution matches expected value
+    for face_idx, expected_value in expected_overlap_contributions.items():
+        assert abs(overlap_contributions[face_idx] - expected_value) < 1e-10, \
+            f"Mismatch for face_index {face_idx}: expected {expected_value}, got {overlap_contributions[face_idx]}"
+
+    # Verify all expected face indices are present
+    assert set(overlap_contributions.keys()) == set(expected_overlap_contributions.keys()), \
+        "Mismatch in face indices"
+
+    # Verify total contributions sum matches total length
+    sum_contributions = sum(overlap_contributions.values())
+    assert abs(sum_contributions - total_length) < 1e-10, \
+        f"Sum of contributions ({sum_contributions}) doesn't match total length ({total_length})"
+
 
 def test_get_zonal_faces_weight_at_constLat_equator():
-    face_0 = [[np.deg2rad(350), np.deg2rad(40)], [np.deg2rad(350), np.deg2rad(20)],
-                  [np.deg2rad(10), np.deg2rad(20)], [np.deg2rad(10), np.deg2rad(40)]]
-    face_1 = [[np.deg2rad(5), np.deg2rad(20)], [np.deg2rad(5), np.deg2rad(10)],
-                [np.deg2rad(25), np.deg2rad(10)], [np.deg2rad(25), np.deg2rad(20)]]
-    face_2 = [[np.deg2rad(30), np.deg2rad(40)], [np.deg2rad(30), np.deg2rad(20)],
-                [np.deg2rad(40), np.deg2rad(20)], [np.deg2rad(40), np.deg2rad(40)]]
+    face_0 = [[1.7 * np.pi, 0.25 * np.pi], [1.7 * np.pi, 0.0],
+              [0.3 * np.pi, 0.0], [0.3 * np.pi, 0.25 * np.pi]]
+    face_1 = [[0.3 * np.pi, 0.0], [0.3 * np.pi, -0.25 * np.pi],
+              [0.6 * np.pi, -0.25 * np.pi], [0.6 * np.pi, 0.0]]
+    face_2 = [[0.3 * np.pi, 0.25 * np.pi], [0.3 * np.pi, 0.0], [np.pi, 0.0],
+              [np.pi, 0.25 * np.pi]]
+    face_3 = [[0.7 * np.pi, 0.0], [0.7 * np.pi, -0.25 * np.pi],
+              [np.pi, -0.25 * np.pi], [np.pi, 0.0]]
 
     # Convert the face vertices to xyz coordinates
     face_0 = [_lonlat_rad_to_xyz(*v) for v in face_0]
     face_1 = [_lonlat_rad_to_xyz(*v) for v in face_1]
     face_2 = [_lonlat_rad_to_xyz(*v) for v in face_2]
-
+    face_3 = [_lonlat_rad_to_xyz(*v) for v in face_3]
 
     face_0_edge_nodes = np.array([[face_0[0], face_0[1]],
-                                    [face_0[1], face_0[2]],
-                                    [face_0[2], face_0[3]],
-                                    [face_0[3], face_0[0]]])
+                                  [face_0[1], face_0[2]],
+                                  [face_0[2], face_0[3]],
+                                  [face_0[3], face_0[0]]])
     face_1_edge_nodes = np.array([[face_1[0], face_1[1]],
-                                    [face_1[1], face_1[2]],
-                                    [face_1[2], face_1[3]],
-                                    [face_1[3], face_1[0]]])
+                                  [face_1[1], face_1[2]],
+                                  [face_1[2], face_1[3]],
+                                  [face_1[3], face_1[0]]])
     face_2_edge_nodes = np.array([[face_2[0], face_2[1]],
-                                    [face_2[1], face_2[2]],
-                                    [face_2[2], face_2[3]],
-                                    [face_2[3], face_2[0]]])
-
-    face_0_latlon_bound = np.array([[np.deg2rad(20), np.deg2rad(40)],
-                                    [np.deg2rad(350), np.deg2rad(10)]])
-    face_1_latlon_bound = np.array([[np.deg2rad(10), np.deg2rad(20)],
-                                    [np.deg2rad(5), np.deg2rad(25)]])
-    face_2_latlon_bound = np.array([[np.deg2rad(20), np.deg2rad(40)],
-                                    [np.deg2rad(30), np.deg2rad(40)]])
+                                  [face_2[1], face_2[2]],
+                                  [face_2[2], face_2[3]],
+                                  [face_2[3], face_2[0]]])
+    face_3_edge_nodes = np.array([[face_3[0], face_3[1]],
+                                  [face_3[1], face_3[2]],
+                                  [face_3[2], face_3[3]],
+                                  [face_3[3], face_3[0]]])
 
+    face_0_latlon_bound = np.array([[0.0, 0.25 * np.pi],
+                                    [1.7 * np.pi, 0.3 * np.pi]])
+    face_1_latlon_bound = np.array([[-0.25 * np.pi, 0.0],
+                                    [0.3 * np.pi, 0.6 * np.pi]])
+    face_2_latlon_bound = np.array([[0.0, 0.25 * np.pi],
+                                    [0.3 * np.pi, np.pi]])
+    face_3_latlon_bound = np.array([[-0.25 * np.pi, 0.0],
+                                    [0.7 * np.pi, np.pi]])
 
     latlon_bounds = np.array([
-        face_0_latlon_bound, face_1_latlon_bound, face_2_latlon_bound
+        face_0_latlon_bound, face_1_latlon_bound, face_2_latlon_bound,
+        face_3_latlon_bound
     ])
 
-    sum = 17.5 + 17.5 + 10
-    expected_weight_df = pd.DataFrame({
-        'face_index': [0, 1, 2],
-        'weight': [17.5 / sum, 17.5/sum, 10/sum]
-    })
-
-    # Assert the results is the same to the 3 decimal places
-    weight_df = _get_zonal_faces_weight_at_constLat(np.array([
-        face_0_edge_nodes, face_1_edge_nodes, face_2_edge_nodes
-    ]), np.sin(np.deg2rad(20)), latlon_bounds, is_latlonface=True)
+    face_edges_cart = np.array([
+        face_0_edge_nodes, face_1_edge_nodes, face_2_edge_nodes,
+        face_3_edge_nodes
+    ])
 
+    constLat_cart = 0.0
 
-    nt.assert_array_almost_equal(weight_df, expected_weight_df, decimal=3)
+    weights = _zonal_face_weights(face_edges_cart,
+                                  latlon_bounds,
+                                  np.array([4, 4, 4, 4]),
+                                  z=constLat_cart,
+                                  check_equator=True)
 
+    expected_weights = np.array([0.46153, 0.11538, 0.30769, 0.11538])
 
+    nt.assert_array_almost_equal(weights, expected_weights, decimal=3)
 
     # A error will be raise if we don't set is_latlonface=True since the face_2 will be concave if
     # It's edges are all GCA
     with pytest.raises(ValueError):
-        _get_zonal_faces_weight_at_constLat(np.array([
+        _zonal_face_weights_robust(np.array([
             face_0_edge_nodes, face_1_edge_nodes, face_2_edge_nodes
         ]), np.deg2rad(20), latlon_bounds)
 
 
 def test_get_zonal_faces_weight_at_constLat_regular():
     face_0 = [[1.7 * np.pi, 0.25 * np.pi], [1.7 * np.pi, 0.0],
-                [0.3 * np.pi, 0.0], [0.3 * np.pi, 0.25 * np.pi]]
+              [0.3 * np.pi, 0.0], [0.3 * np.pi, 0.25 * np.pi]]
     face_1 = [[0.4 * np.pi, 0.3 * np.pi], [0.4 * np.pi, 0.0],
-                [0.5 * np.pi, 0.0], [0.5 * np.pi, 0.3 * np.pi]]
+              [0.5 * np.pi, 0.0], [0.5 * np.pi, 0.3 * np.pi]]
     face_2 = [[0.5 * np.pi, 0.25 * np.pi], [0.5 * np.pi, 0.0], [np.pi, 0.0],
-                [np.pi, 0.25 * np.pi]]
+              [np.pi, 0.25 * np.pi]]
     face_3 = [[1.2 * np.pi, 0.25 * np.pi], [1.2 * np.pi, 0.0],
-                [1.6 * np.pi, -0.01 * np.pi], [1.6 * np.pi, 0.25 * np.pi]]
+              [1.6 * np.pi, -0.01 * np.pi], [1.6 * np.pi, 0.25 * np.pi]]
 
     # Convert the face vertices to xyz coordinates
     face_0 = [_lonlat_rad_to_xyz(*v) for v in face_0]
@@ -495,21 +574,21 @@ def test_get_zonal_faces_weight_at_constLat_regular():
     face_3 = [_lonlat_rad_to_xyz(*v) for v in face_3]
 
     face_0_edge_nodes = np.array([[face_0[0], face_0[1]],
-                                    [face_0[1], face_0[2]],
-                                    [face_0[2], face_0[3]],
-                                    [face_0[3], face_0[0]]])
+                                  [face_0[1], face_0[2]],
+                                  [face_0[2], face_0[3]],
+                                  [face_0[3], face_0[0]]])
     face_1_edge_nodes = np.array([[face_1[0], face_1[1]],
-                                    [face_1[1], face_1[2]],
-                                    [face_1[2], face_1[3]],
-                                    [face_1[3], face_1[0]]])
+                                  [face_1[1], face_1[2]],
+                                  [face_1[2], face_1[3]],
+                                  [face_1[3], face_1[0]]])
     face_2_edge_nodes = np.array([[face_2[0], face_2[1]],
-                                    [face_2[1], face_2[2]],
-                                    [face_2[2], face_2[3]],
-                                    [face_2[3], face_2[0]]])
+                                  [face_2[1], face_2[2]],
+                                  [face_2[2], face_2[3]],
+                                  [face_2[3], face_2[0]]])
     face_3_edge_nodes = np.array([[face_3[0], face_3[1]],
-                                    [face_3[1], face_3[2]],
-                                    [face_3[2], face_3[3]],
-                                    [face_3[3], face_3[0]]])
+                                  [face_3[1], face_3[2]],
+                                  [face_3[2], face_3[3]],
+                                  [face_3[3], face_3[0]]])
 
     face_0_latlon_bound = np.array([[0.0, 0.25 * np.pi],
                                     [1.7 * np.pi, 0.3 * np.pi]])
@@ -525,33 +604,258 @@ def test_get_zonal_faces_weight_at_constLat_regular():
         face_3_latlon_bound
     ])
 
-    expected_weight_df = pd.DataFrame({
-        'face_index': [0, 1, 2, 3],
-        'weight': [0.375, 0.0625, 0.3125, 0.25]
-    })
-
-    # Assert the results is the same to the 3 decimal places
-    weight_df = _get_zonal_faces_weight_at_constLat(np.array([
+    face_edges_cart = np.array([
         face_0_edge_nodes, face_1_edge_nodes, face_2_edge_nodes,
         face_3_edge_nodes
-    ]), np.sin(0.1 * np.pi), latlon_bounds)
+    ])
+
+    constLat_cart = np.sin(0.1 * np.pi)
+
+    weights = _zonal_face_weights(face_edges_cart,
+                                  latlon_bounds,
+                                  np.array([4, 4, 4, 4]),
+                                  z=constLat_cart)
+
+    expected_weights = np.array([0.375, 0.0625, 0.3125, 0.25])
+
+    nt.assert_array_almost_equal(weights, expected_weights)
+
+    # expected_weight_df = pd.DataFrame({
+    #     'face_index': [0, 1, 2, 3],
+    #     'weight': [0.375, 0.0625, 0.3125, 0.25]
+    # })
+    #
+    # # Assert the results is the same to the 3 decimal places
+    # weight_df = _get_zonal_faces_weight_at_constLat(np.array([
+    #     face_0_edge_nodes, face_1_edge_nodes, face_2_edge_nodes,
+    #     face_3_edge_nodes
+    # ]), np.sin(0.1 * np.pi), latlon_bounds)
+    #
+    # nt.assert_array_almost_equal(weight_df, expected_weight_df, decimal=3)
 
-    nt.assert_array_almost_equal(weight_df, expected_weight_df, decimal=3)
 
 def test_get_zonal_faces_weight_at_constLat_on_pole_one_face():
-    #The face is touching the pole， so the weight should be 1.0 since there's only 1 face
+    # The face is touching the pole, so the weight should be 1.0 since there's only 1 face
     face_edges_cart = np.array([[
         [[-5.22644277e-02, -5.22644277e-02, -9.97264689e-01],
-            [-5.23359562e-02, -6.40930613e-18, -9.98629535e-01]],
-
+         [-5.23359562e-02, -6.40930613e-18, -9.98629535e-01]],
         [[-5.23359562e-02, -6.40930613e-18, -9.98629535e-01],
-            [6.12323400e-17, 0.00000000e+00, -1.00000000e+00]],
+         [6.12323400e-17, 0.00000000e+00, -1.00000000e+00]],
+        [[6.12323400e-17, 0.00000000e+00, -1.00000000e+00],
+         [3.20465306e-18, -5.23359562e-02, -9.98629535e-01]],
+        [[3.20465306e-18, -5.23359562e-02, -9.98629535e-01],
+         [-5.22644277e-02, -5.22644277e-02, -9.97264689e-01]]
+    ]])
+
+    # Corrected face_bounds
+    face_bounds = np.array([
+        [-1.57079633, -1.4968158],
+        [3.14159265, 0.]
+    ])
+    constLat_cart = -1
+
+    weights = _zonal_face_weights(face_edges_cart,
+                                  np.array([4, 4, 4, 4]),
+                                  z=constLat_cart)
+
+    expected_weights = np.array([1.0])
+
+    nt.assert_array_equal(weights, expected_weights)
+
+    # weight_df = _get_zonal_faces_weight_at_constLat(face_edges_cart, constLat_cart, face_bounds)
+    #
+    # # Create expected Polars DataFrame
+    # expected_weight_df = pl.DataFrame(
+    #     {
+    #         "face_index": pl.Series([0], dtype=pl.Int64),
+    #         "weight": pl.Series([1.0], dtype=pl.Float64)
+    #     }
+    # )
+    #
+    # # Assert equality using Polars
+    # assert_frame_equal(weight_df, expected_weight_df), \
+    #     f"Expected:\n{expected_weight_df}\nGot:\n{weight_df}"
+
+
+def test_get_zonal_faces_weight_at_constLat_on_pole_faces():
+    # There will be 4 faces touching the pole, so the weight should be 0.25 for each face
+    face_edges_cart = np.array([
+        [
+            [[5.22644277e-02, -5.22644277e-02, 9.97264689e-01], [5.23359562e-02, 0.00000000e+00, 9.98629535e-01]],
+            [[5.23359562e-02, 0.00000000e+00, 9.98629535e-01], [6.12323400e-17, 0.00000000e+00, 1.00000000e+00]],
+            [[6.12323400e-17, 0.00000000e+00, 1.00000000e+00], [3.20465306e-18, -5.23359562e-02, 9.98629535e-01]],
+            [[3.20465306e-18, -5.23359562e-02, 9.98629535e-01], [5.22644277e-02, -5.22644277e-02, 9.97264689e-01]]
+        ],
+        [
+            [[5.23359562e-02, 0.00000000e+00, 9.98629535e-01], [5.22644277e-02, 5.22644277e-02, 9.97264689e-01]],
+            [[5.22644277e-02, 5.22644277e-02, 9.97264689e-01], [3.20465306e-18, 5.23359562e-02, 9.98629535e-01]],
+            [[3.20465306e-18, 5.23359562e-02, 9.98629535e-01], [6.12323400e-17, 0.00000000e+00, 1.00000000e+00]],
+            [[6.12323400e-17, 0.00000000e+00, 1.00000000e+00], [5.23359562e-02, 0.00000000e+00, 9.98629535e-01]]
+        ],
+        [
+            [[3.20465306e-18, -5.23359562e-02, 9.98629535e-01], [6.12323400e-17, 0.00000000e+00, 1.00000000e+00]],
+            [[6.12323400e-17, 0.00000000e+00, 1.00000000e+00], [-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], [3.20465306e-18, -5.23359562e-02, 9.98629535e-01]]
+        ],
+        [
+            [[6.12323400e-17, 0.00000000e+00, 1.00000000e+00], [3.20465306e-18, 5.23359562e-02, 9.98629535e-01]],
+            [[3.20465306e-18, 5.23359562e-02, 9.98629535e-01], [-5.22644277e-02, 5.22644277e-02, 9.97264689e-01]],
+            [[-5.22644277e-02, 5.22644277e-02, 9.97264689e-01], [-5.23359562e-02, -6.40930613e-18, 9.98629535e-01]],
+            [[-5.23359562e-02, -6.40930613e-18, 9.98629535e-01], [6.12323400e-17, 0.00000000e+00, 1.00000000e+00]]
+        ]
+    ])
+
+    face_bounds = np.array([
+        [[1.4968158, 1.57079633], [4.71238898, 0.0]],
+        [[1.4968158, 1.57079633], [0.0, 1.57079633]],
+        [[1.4968158, 1.57079633], [3.14159265, 0.0]],
+        [[1.4968158, 1.57079633], [0.0, 3.14159265]]
+    ])
+
+    constLat_cart = 1.0
+
+    weights = _zonal_face_weights(face_edges_cart,
+                                  np.array([4, 4, 4, 4]),
+                                  z=constLat_cart)
 
+    expected_weights = np.array([0.25, 0.25, 0.25, 0.25])
+
+    nt.assert_array_equal(weights, expected_weights)
+
+    #
+    # weight_df = _get_zonal_faces_weight_at_constLat(face_edges_cart, constLat_cart, face_bounds)
+    #
+    # # Create expected Polars DataFrame
+    # expected_weight_df = pl.DataFrame(
+    #     {
+    #         'face_index': pl.Series([0, 1, 2, 3], dtype=pl.Int64),
+    #         'weight': pl.Series([0.25, 0.25, 0.25, 0.25], dtype=pl.Float64)
+    #     }
+    # )
+    #
+    # # Assert equality using Polars
+    # assert_frame_equal(weight_df, expected_weight_df), \
+    #     f"Expected:\n{expected_weight_df}\nGot:\n{weight_df}"
+
+
+def test_get_zonal_faces_weight_at_constLat_on_pole_one_face():
+    # The face is touching the pole, so the weight should be 1.0 since there's only 1 face
+    face_edges_cart = np.array([[
+        [[-5.22644277e-02, -5.22644277e-02, -9.97264689e-01],
+         [-5.23359562e-02, -6.40930613e-18, -9.98629535e-01]],
+        [[-5.23359562e-02, -6.40930613e-18, -9.98629535e-01],
+         [6.12323400e-17, 0.00000000e+00, -1.00000000e+00]],
         [[6.12323400e-17, 0.00000000e+00, -1.00000000e+00],
-            [3.20465306e-18, -5.23359562e-02, -9.98629535e-01]],
+         [3.20465306e-18, -5.23359562e-02, -9.98629535e-01]],
+        [[3.20465306e-18, -5.23359562e-02, -9.98629535e-01],
+         [-5.22644277e-02, -5.22644277e-02, -9.97264689e-01]]
+    ]])
 
+    # Corrected face_bounds
+    face_bounds = np.array([
+        [-1.57079633, -1.4968158],
+        [3.14159265, 0.]
+    ])
+    constLat_cart = -1
+
+    weights = _zonal_face_weights(face_edges_cart,
+                                  np.array([4, 4, 4, 4]),
+                                  z=constLat_cart)
+
+    expected_weights = np.array([0.25, 0.25, 0.25, 0.25])
+
+    nt.assert_array_equal(weights, expected_weights)
+
+    # weight_df = _get_zonal_faces_weight_at_constLat(face_edges_cart, constLat_cart, face_bounds)
+    #
+    # # Create expected Polars DataFrame
+    # expected_weight_df = pl.DataFrame(
+    #     {
+    #         "face_index": pl.Series([0], dtype=pl.Int64),
+    #         "weight": pl.Series([1.0], dtype=pl.Float64)
+    #     }
+    # )
+    #
+    # # Assert equality by comparing columns
+    # assert (weight_df.select(pl.all()).collect() == expected_weight_df.select(pl.all()).collect()), \
+    #     f"Expected:\n{expected_weight_df}\nGot:\n{weight_df}"
+
+
+def test_get_zonal_faces_weight_at_constLat_on_pole_faces():
+    # There will be 4 faces touching the pole, so the weight should be 0.25 for each face
+    face_edges_cart = np.array([
+        [
+            [[5.22644277e-02, -5.22644277e-02, 9.97264689e-01], [5.23359562e-02, 0.00000000e+00, 9.98629535e-01]],
+            [[5.23359562e-02, 0.00000000e+00, 9.98629535e-01], [6.12323400e-17, 0.00000000e+00, 1.00000000e+00]],
+            [[6.12323400e-17, 0.00000000e+00, 1.00000000e+00], [3.20465306e-18, -5.23359562e-02, 9.98629535e-01]],
+            [[3.20465306e-18, -5.23359562e-02, 9.98629535e-01], [5.22644277e-02, -5.22644277e-02, 9.97264689e-01]]
+        ],
+        [
+            [[5.23359562e-02, 0.00000000e+00, 9.98629535e-01], [5.22644277e-02, 5.22644277e-02, 9.97264689e-01]],
+            [[5.22644277e-02, 5.22644277e-02, 9.97264689e-01], [3.20465306e-18, 5.23359562e-02, 9.98629535e-01]],
+            [[3.20465306e-18, 5.23359562e-02, 9.98629535e-01], [6.12323400e-17, 0.00000000e+00, 1.00000000e+00]],
+            [[6.12323400e-17, 0.00000000e+00, 1.00000000e+00], [5.23359562e-02, 0.00000000e+00, 9.98629535e-01]]
+        ],
+        [
+            [[3.20465306e-18, -5.23359562e-02, 9.98629535e-01], [6.12323400e-17, 0.00000000e+00, 1.00000000e+00]],
+            [[6.12323400e-17, 0.00000000e+00, 1.00000000e+00], [-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], [3.20465306e-18, -5.23359562e-02, 9.98629535e-01]]
+        ],
+        [
+            [[6.12323400e-17, 0.00000000e+00, 1.00000000e+00], [3.20465306e-18, 5.23359562e-02, 9.98629535e-01]],
+            [[3.20465306e-18, 5.23359562e-02, 9.98629535e-01], [-5.22644277e-02, 5.22644277e-02, 9.97264689e-01]],
+            [[-5.22644277e-02, 5.22644277e-02, 9.97264689e-01], [-5.23359562e-02, -6.40930613e-18, 9.98629535e-01]],
+            [[-5.23359562e-02, -6.40930613e-18, 9.98629535e-01], [6.12323400e-17, 0.00000000e+00, 1.00000000e+00]]
+        ]
+    ])
+
+    face_bounds = np.array([
+        [[1.4968158, 1.57079633], [4.71238898, 0.0]],
+        [[1.4968158, 1.57079633], [0.0, 1.57079633]],
+        [[1.4968158, 1.57079633], [3.14159265, 0.0]],
+        [[1.4968158, 1.57079633], [0.0, 3.14159265]]
+    ])
+
+    constLat_cart = 1.0
+
+    weights = _zonal_face_weights(face_edges_cart,
+                                  np.array([4, 4, 4, 4]),
+                                  z=constLat_cart)
+
+    expected_weights = np.array([0.25, 0.25, 0.25, 0.25])
+
+    nt.assert_array_equal(weights, expected_weights)
+
+    # weight_df = _get_zonal_faces_weight_at_constLat(face_edges_cart, constLat_cart, face_bounds)
+    #
+    # # Create expected Polars DataFrame
+    # expected_weight_df = pl.DataFrame(
+    #     {
+    #         'face_index': pl.Series([0, 1, 2, 3], dtype=pl.Int64),
+    #         'weight': pl.Series([0.25, 0.25, 0.25, 0.25], dtype=pl.Float64)
+    #     }
+    # )
+    #
+    # # Assert equality by comparing columns
+    # assert (weight_df.select(pl.all()).collect() == expected_weight_df.select(pl.all()).collect()), \
+    #     f"Expected:\n{expected_weight_df}\nGot:\n{weight_df}"
+
+
+def test_get_zonal_faces_weight_at_constLat_on_pole_one_face():
+    # The face is touching the pole, so the weight should be 1.0 since there's only 1 face
+    face_edges_cart = np.array([[
+        [[-5.22644277e-02, -5.22644277e-02, -9.97264689e-01],
+         [-5.23359562e-02, -6.40930613e-18, -9.98629535e-01]],
+        [[-5.23359562e-02, -6.40930613e-18, -9.98629535e-01],
+         [6.12323400e-17, 0.00000000e+00, -1.00000000e+00]],
+        [[6.12323400e-17, 0.00000000e+00, -1.00000000e+00],
+         [3.20465306e-18, -5.23359562e-02, -9.98629535e-01]],
         [[3.20465306e-18, -5.23359562e-02, -9.98629535e-01],
-            [-5.22644277e-02, -5.22644277e-02, -9.97264689e-01]]
+         [-5.22644277e-02, -5.22644277e-02, -9.97264689e-01]]
     ]])
 
     # Corrected face_bounds
@@ -561,16 +865,30 @@ def test_get_zonal_faces_weight_at_constLat_on_pole_one_face():
     ])
     constLat_cart = -1
 
-    weight_df = _get_zonal_faces_weight_at_constLat(face_edges_cart, constLat_cart, face_bounds)
-    # Define the expected DataFrame
-    expected_weight_df = pd.DataFrame({"face_index": [0], "weight": [1.0]})
+    # weight_df = _get_zonal_faces_weight_at_constLat(face_edges_cart, constLat_cart, face_bounds)
+    #
+    # # Create expected Polars DataFrame
+    # expected_weight_df = pl.DataFrame(
+    #     {
+    #         "face_index": pl.Series([0], dtype=pl.Int64),
+    #         "weight": pl.Series([1.0], dtype=pl.Float64)
+    #     }
+    # )
+    #
+    # assert_frame_equal(weight_df, expected_weight_df)
+
+    weights = _zonal_face_weights(face_edges_cart,
+                                  face_bounds,
+                                  np.array([4, 4, 4, 4]),
+                                  z=constLat_cart)
 
-    # Assert that the resulting should have weight is 1.0
-    pd.testing.assert_frame_equal(weight_df, expected_weight_df)
+    expected_weights = np.array([1.0, ])
+
+    nt.assert_array_equal(weights, expected_weights)
 
 
 def test_get_zonal_faces_weight_at_constLat_on_pole_faces():
-    #there will be 4 faces touching the pole, so the weight should be 0.25 for each face
+    # There will be 4 faces touching the pole, so the weight should be 0.25 for each face
     face_edges_cart = np.array([
         [
             [[5.22644277e-02, -5.22644277e-02, 9.97264689e-01], [5.23359562e-02, 0.00000000e+00, 9.98629535e-01]],
@@ -588,7 +906,7 @@ def test_get_zonal_faces_weight_at_constLat_on_pole_faces():
             [[3.20465306e-18, -5.23359562e-02, 9.98629535e-01], [6.12323400e-17, 0.00000000e+00, 1.00000000e+00]],
             [[6.12323400e-17, 0.00000000e+00, 1.00000000e+00], [-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]],
             [[-5.22644277e-02, -5.22644277e-02, 9.97264689e-01], [3.20465306e-18, -5.23359562e-02, 9.98629535e-01]]
         ],
         [
@@ -608,31 +926,29 @@ def test_get_zonal_faces_weight_at_constLat_on_pole_faces():
 
     constLat_cart = 1.0
 
-    weight_df = _get_zonal_faces_weight_at_constLat(face_edges_cart, constLat_cart, face_bounds)
-    # Define the expected DataFrame
-    expected_weight_df = pd.DataFrame({
-        'face_index': [0, 1, 2, 3],
-        'weight': [0.25, 0.25, 0.25, 0.25]
-    })
+    # weight_df = _get_zonal_faces_weight_at_constLat(face_edges_cart, constLat_cart, face_bounds)
 
-    # Assert that the DataFrame matches the expected DataFrame
-    pd.testing.assert_frame_equal(weight_df, expected_weight_df)
+    weights = _zonal_face_weights(face_edges_cart,
+                                  face_bounds,
+                                  np.array([4, 4, 4, 4]),
+                                  z=constLat_cart)
+
+    expected_weights = np.array([0.25, 0.25, 0.25, 0.25])
+
+    nt.assert_array_equal(weights, expected_weights)
 
 
 def test_get_zonal_face_interval_pole():
-    #The face is touching the pole
+    # The face is touching the pole
     face_edges_cart = np.array([
         [[-5.22644277e-02, -5.22644277e-02, -9.97264689e-01],
-            [-5.23359562e-02, -6.40930613e-18, -9.98629535e-01]],
-
+         [-5.23359562e-02, -6.40930613e-18, -9.98629535e-01]],
         [[-5.23359562e-02, -6.40930613e-18, -9.98629535e-01],
-            [6.12323400e-17, 0.00000000e+00, -1.00000000e+00]],
-
+         [6.12323400e-17, 0.00000000e+00, -1.00000000e+00]],
         [[6.12323400e-17, 0.00000000e+00, -1.00000000e+00],
-            [3.20465306e-18, -5.23359562e-02, -9.98629535e-01]],
-
+         [3.20465306e-18, -5.23359562e-02, -9.98629535e-01]],
         [[3.20465306e-18, -5.23359562e-02, -9.98629535e-01],
-            [-5.22644277e-02, -5.22644277e-02, -9.97264689e-01]]
+         [-5.22644277e-02, -5.22644277e-02, -9.97264689e-01]]
     ])
 
     # Corrected face_bounds
@@ -643,36 +959,38 @@ def test_get_zonal_face_interval_pole():
     constLat_cart = -0.9986295347545738
 
     weight_df = _get_zonal_face_interval(face_edges_cart, constLat_cart, face_bounds)
-    # No Nan values should be present in the weight_df
-    assert (not weight_df.isnull().values.any())
+    df_null_counts = weight_df.null_count()
+
+    total_nulls = df_null_counts.to_numpy().sum()
+
+    assert total_nulls == 0, f"Found {total_nulls} null values in the DataFrame"
 
 
 def test_get_zonal_faces_weight_at_constLat_latlonface():
     face_0 = [[np.deg2rad(350), np.deg2rad(40)], [np.deg2rad(350), np.deg2rad(20)],
-                [np.deg2rad(10), np.deg2rad(20)], [np.deg2rad(10), np.deg2rad(40)]]
+              [np.deg2rad(10), np.deg2rad(20)], [np.deg2rad(10), np.deg2rad(40)]]
     face_1 = [[np.deg2rad(5), np.deg2rad(20)], [np.deg2rad(5), np.deg2rad(10)],
-                [np.deg2rad(25), np.deg2rad(10)], [np.deg2rad(25), np.deg2rad(20)]]
+              [np.deg2rad(25), np.deg2rad(10)], [np.deg2rad(25), np.deg2rad(20)]]
     face_2 = [[np.deg2rad(30), np.deg2rad(40)], [np.deg2rad(30), np.deg2rad(20)],
-                [np.deg2rad(40), np.deg2rad(20)], [np.deg2rad(40), np.deg2rad(40)]]
+              [np.deg2rad(40), np.deg2rad(20)], [np.deg2rad(40), np.deg2rad(40)]]
 
     # Convert the face vertices to xyz coordinates
     face_0 = [_lonlat_rad_to_xyz(*v) for v in face_0]
     face_1 = [_lonlat_rad_to_xyz(*v) for v in face_1]
     face_2 = [_lonlat_rad_to_xyz(*v) for v in face_2]
 
-
     face_0_edge_nodes = np.array([[face_0[0], face_0[1]],
-                                    [face_0[1], face_0[2]],
-                                    [face_0[2], face_0[3]],
-                                    [face_0[3], face_0[0]]])
+                                  [face_0[1], face_0[2]],
+                                  [face_0[2], face_0[3]],
+                                  [face_0[3], face_0[0]]])
     face_1_edge_nodes = np.array([[face_1[0], face_1[1]],
-                                    [face_1[1], face_1[2]],
-                                    [face_1[2], face_1[3]],
-                                    [face_1[3], face_1[0]]])
+                                  [face_1[1], face_1[2]],
+                                  [face_1[2], face_1[3]],
+                                  [face_1[3], face_1[0]]])
     face_2_edge_nodes = np.array([[face_2[0], face_2[1]],
-                                    [face_2[1], face_2[2]],
-                                    [face_2[2], face_2[3]],
-                                    [face_2[3], face_2[0]]])
+                                  [face_2[1], face_2[2]],
+                                  [face_2[2], face_2[3]],
+                                  [face_2[3], face_2[0]]])
 
     face_0_latlon_bound = np.array([[np.deg2rad(20), np.deg2rad(40)],
                                     [np.deg2rad(350), np.deg2rad(10)]])
@@ -681,7 +999,6 @@ def test_get_zonal_faces_weight_at_constLat_latlonface():
     face_2_latlon_bound = np.array([[np.deg2rad(20), np.deg2rad(40)],
                                     [np.deg2rad(30), np.deg2rad(40)]])
 
-
     latlon_bounds = np.array([
         face_0_latlon_bound, face_1_latlon_bound, face_2_latlon_bound
     ])
@@ -689,22 +1006,65 @@ def test_get_zonal_faces_weight_at_constLat_latlonface():
     sum = 17.5 + 17.5 + 10
     expected_weight_df = pd.DataFrame({
         'face_index': [0, 1, 2],
-        'weight': [17.5 / sum, 17.5/sum, 10/sum]
+        'weight': [17.5 / sum, 17.5 / sum, 10 / sum]
     })
 
     # Assert the results is the same to the 3 decimal places
-    weight_df = _get_zonal_faces_weight_at_constLat(np.array([
+    weight_df = _zonal_face_weights_robust(np.array([
         face_0_edge_nodes, face_1_edge_nodes, face_2_edge_nodes
     ]), np.sin(np.deg2rad(20)), latlon_bounds, is_latlonface=True)
 
-
     nt.assert_array_almost_equal(weight_df, expected_weight_df, decimal=3)
 
-
-
     # A error will be raise if we don't set is_latlonface=True since the face_2 will be concave if
     # It's edges are all GCA
     with pytest.raises(ValueError):
-        _get_zonal_faces_weight_at_constLat(np.array([
+        _zonal_face_weights_robust(np.array([
             face_0_edge_nodes, face_1_edge_nodes, face_2_edge_nodes
         ]), np.deg2rad(20), latlon_bounds)
+
+
+def test_compare_zonal_weights():
+    """Compares the existing weights calculation (get_non_conservative_zonal_face_weights_at_const_lat_overlap) to
+    the faster implementation (get_non_conservative_zonal_face_weights_at_const_lat)"""
+    gridfiles = [current_path / "meshfiles" / "ugrid" / "outCSne30" / "outCSne30.ug",
+                 current_path / "meshfiles" / "scrip" / "outCSne8" / "outCSne8.nc",
+                 current_path / "meshfiles" / "ugrid" / "geoflow-small" / "grid.nc",]
+
+    lat = (-90, 90, 10)
+    latitudes = np.arange(lat[0], lat[1] + lat[2], lat[2])
+
+    for gridfile in gridfiles:
+        uxgrid = ux.open_grid(gridfile)
+        n_nodes_per_face = uxgrid.n_nodes_per_face.values
+        face_edge_nodes_xyz =  _get_cartesian_face_edge_nodes(
+                uxgrid.face_node_connectivity.values,
+                uxgrid.n_face,
+                uxgrid.n_max_face_edges,
+                uxgrid.node_x.values,
+                uxgrid.node_y.values,
+                uxgrid.node_z.values,
+            )
+        bounds = uxgrid.bounds.values
+
+        for i, lat in enumerate(latitudes):
+            face_indices = uxgrid.get_faces_at_constant_latitude(lat)
+            z = np.sin(np.deg2rad(lat))
+
+            face_edge_nodes_xyz_candidate = face_edge_nodes_xyz[face_indices, :, :, :]
+            n_nodes_per_face_candidate = n_nodes_per_face[face_indices]
+            bounds_candidate = bounds[face_indices]
+
+            new_weights = _zonal_face_weights(face_edge_nodes_xyz_candidate,
+                                              bounds_candidate,
+                                              n_nodes_per_face_candidate,
+                                              z)
+
+            existing_weights = _zonal_face_weights_robust(
+                face_edge_nodes_xyz_candidate, z, bounds_candidate
+            )["weight"].to_numpy()
+
+            abs_diff = np.abs(new_weights - existing_weights)
+
+            # For each latitude, make sure the aboslute difference is below our error tollerance
+            assert abs_diff.max() < ERROR_TOLERANCE
diff --git a/test/test_intersections.py b/test/test_intersections.py
index efe516dce..f2fbfe1ee 100644
--- a/test/test_intersections.py
+++ b/test/test_intersections.py
@@ -2,9 +2,9 @@
 import pytest
 import uxarray as ux
 from uxarray.constants import ERROR_TOLERANCE
-from uxarray.grid.arcs import _extreme_gca_latitude_cartesian
+from uxarray.grid.arcs import extreme_gca_z
 from uxarray.grid.coordinates import _lonlat_rad_to_xyz, _xyz_to_lonlat_rad,_xyz_to_lonlat_rad_scalar
-from uxarray.grid.intersections import gca_gca_intersection, gca_const_lat_intersection, _gca_gca_intersection_cartesian
+from uxarray.grid.intersections import gca_gca_intersection, gca_const_lat_intersection, _gca_gca_intersection_cartesian, get_number_of_intersections
 
 def test_get_GCA_GCA_intersections_antimeridian():
     GCA1 = _lonlat_rad_to_xyz(np.deg2rad(170.0), np.deg2rad(89.99))
@@ -233,7 +233,7 @@ def test_GCA_constLat_intersections_antimeridian():
         _lonlat_rad_to_xyz(np.deg2rad(170.0), np.deg2rad(10.0))
     ])
 
-    res = gca_const_lat_intersection(GCR1_cart, np.sin(np.deg2rad(60.0)), verbose=True)
+    res = gca_const_lat_intersection(GCR1_cart, np.sin(np.deg2rad(60.0)))
     res_lonlat_rad = _xyz_to_lonlat_rad(*(res[0].tolist()))
     assert np.allclose(res_lonlat_rad, np.array([np.deg2rad(170.0), np.deg2rad(60.0)]))
 
@@ -243,27 +243,17 @@ def test_GCA_constLat_intersections_empty():
         _lonlat_rad_to_xyz(np.deg2rad(170.0), np.deg2rad(10.0))
     ])
 
-    res = gca_const_lat_intersection(GCR1_cart, np.sin(np.deg2rad(-10.0)), verbose=False)
-    assert res.size == 0
+    res = gca_const_lat_intersection(GCR1_cart, np.sin(np.deg2rad(-10.0)))
+    assert get_number_of_intersections(res) == 0
 
 def test_GCA_constLat_intersections_two_pts():
     GCR1_cart = np.array([
         _lonlat_rad_to_xyz(np.deg2rad(10.0), np.deg2rad(10)),
         _lonlat_rad_to_xyz(np.deg2rad(170.0), np.deg2rad(10.0))
     ])
-    max_lat = _extreme_gca_latitude_cartesian(GCR1_cart, 'max')
+    max_lat = extreme_gca_z(GCR1_cart, 'max')
 
     query_lat = (np.deg2rad(10.0) + max_lat) / 2.0
 
-    res = gca_const_lat_intersection(GCR1_cart, np.sin(query_lat), verbose=False)
+    res = gca_const_lat_intersection(GCR1_cart, np.sin(query_lat))
     assert res.shape[0] == 2
-
-def test_GCA_constLat_intersections_no_converge():
-    GCR1_cart = np.array([[-0.59647278, 0.59647278, -0.53706651],
-                          [-0.61362973, 0.61362973, -0.49690755]])
-
-    constZ = -0.5150380749100542
-
-    with pytest.warns(UserWarning):
-        res = gca_const_lat_intersection(GCR1_cart, constZ, verbose=False)
-        assert res.shape[0] == 1