diff --git a/test/test_geometry.py b/test/test_geometry.py
index c066f612a..0c729374e 100644
--- a/test/test_geometry.py
+++ b/test/test_geometry.py
@@ -12,7 +12,7 @@
 from uxarray.grid.coordinates import _populate_node_latlon, _lonlat_rad_to_xyz, _normalize_xyz, _xyz_to_lonlat_rad, \
     _xyz_to_lonlat_deg, _xyz_to_lonlat_rad_scalar
 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.utils import _get_cartesian_faces_edge_nodes, _get_lonlat_rad_faces_edge_nodes
 
 from uxarray.grid.geometry import _populate_face_latlon_bound, _populate_bounds, _pole_point_inside_polygon_cartesian, \
     stereographic_projection, inverse_stereographic_projection, point_in_face, haversine_distance
@@ -1184,14 +1184,9 @@ def test_point_inside():
     grid = ux.open_grid(grid_mpas_2)
 
     # Get the face edges of all faces in the grid
-    faces_edges_cartesian = _get_cartesian_face_edge_nodes(
-        grid.face_node_connectivity.values,
-        grid.n_face,
-        grid.n_max_face_edges,
-        grid.node_x.values,
-        grid.node_y.values,
-        grid.node_z.values,
-    )
+    faces_edges_cartesian = _get_cartesian_faces_edge_nodes(grid.face_node_connectivity.values, grid.n_face,
+                                                            grid.n_max_face_edges, grid.node_x.values,
+                                                            grid.node_y.values, grid.node_z.values)
 
     # Loop through each face
     for i in range(grid.n_face):
@@ -1209,14 +1204,9 @@ def test_point_outside():
     grid = ux.open_grid(grid_mpas_2)
 
     # Get the face edges of all faces in the grid
-    faces_edges_cartesian = _get_cartesian_face_edge_nodes(
-        grid.face_node_connectivity.values,
-        grid.n_face,
-        grid.n_max_face_edges,
-        grid.node_x.values,
-        grid.node_y.values,
-        grid.node_z.values,
-    )
+    faces_edges_cartesian = _get_cartesian_faces_edge_nodes(grid.face_node_connectivity.values, grid.n_face,
+                                                            grid.n_max_face_edges, grid.node_x.values,
+                                                            grid.node_y.values, grid.node_z.values)
 
     # Set the point as the face center of a different face than the face tested
     point_xyz = np.array([grid.face_x[1].values, grid.face_y[1].values, grid.face_z[1].values])
@@ -1232,14 +1222,9 @@ def test_point_on_node():
     grid = ux.open_grid(grid_mpas_2)
 
     # Get the face edges of all faces in the grid
-    faces_edges_cartesian = _get_cartesian_face_edge_nodes(
-        grid.face_node_connectivity.values,
-        grid.n_face,
-        grid.n_max_face_edges,
-        grid.node_x.values,
-        grid.node_y.values,
-        grid.node_z.values,
-    )
+    faces_edges_cartesian = _get_cartesian_faces_edge_nodes(grid.face_node_connectivity.values, grid.n_face,
+                                                            grid.n_max_face_edges, grid.node_x.values,
+                                                            grid.node_y.values, grid.node_z.values)
 
     # Set the point as a node
     point_xyz = np.array([*faces_edges_cartesian[0][0][0]])
@@ -1263,14 +1248,9 @@ def test_point_inside_close():
 
     # Create the grid and face edges
     grid = ux.Grid.from_face_vertices(vertices_lonlat, latlon=True)
-    faces_edges_cartesian = _get_cartesian_face_edge_nodes(
-        grid.face_node_connectivity.values,
-        grid.n_face,
-        grid.n_max_face_edges,
-        grid.node_x.values,
-        grid.node_y.values,
-        grid.node_z.values,
-    )
+    faces_edges_cartesian = _get_cartesian_faces_edge_nodes(grid.face_node_connectivity.values, grid.n_face,
+                                                            grid.n_max_face_edges, grid.node_x.values,
+                                                            grid.node_y.values, grid.node_z.values)
 
     # Use point in face to determine if the point is inside or out of the face
     assert point_in_face(faces_edges_cartesian[0], point_xyz=point, inclusive=False)
@@ -1288,14 +1268,9 @@ def test_point_outside_close():
 
     # Create the grid and face edges
     grid = ux.Grid.from_face_vertices(vertices_lonlat, latlon=True)
-    faces_edges_cartesian = _get_cartesian_face_edge_nodes(
-        grid.face_node_connectivity.values,
-        grid.n_face,
-        grid.n_max_face_edges,
-        grid.node_x.values,
-        grid.node_y.values,
-        grid.node_z.values,
-    )
+    faces_edges_cartesian = _get_cartesian_faces_edge_nodes(grid.face_node_connectivity.values, grid.n_face,
+                                                            grid.n_max_face_edges, grid.node_x.values,
+                                                            grid.node_y.values, grid.node_z.values)
 
     # Use point in face to determine if the point is inside or out of the face
     assert not point_in_face(faces_edges_cartesian[0], point_xyz=point, inclusive=False)
@@ -1312,14 +1287,9 @@ def test_face_at_pole():
 
     # Create the grid and face edges
     grid = ux.Grid.from_face_vertices(vertices_lonlat, latlon=True)
-    faces_edges_cartesian = _get_cartesian_face_edge_nodes(
-        grid.face_node_connectivity.values,
-        grid.n_face,
-        grid.n_max_face_edges,
-        grid.node_x.values,
-        grid.node_y.values,
-        grid.node_z.values,
-    )
+    faces_edges_cartesian = _get_cartesian_faces_edge_nodes(grid.face_node_connectivity.values, grid.n_face,
+                                                            grid.n_max_face_edges, grid.node_x.values,
+                                                            grid.node_y.values, grid.node_z.values)
 
     assert point_in_face(faces_edges_cartesian[0], point_xyz=point, inclusive=True)
 
@@ -1334,14 +1304,9 @@ def test_face_at_antimeridian():
 
     # Create the grid and face edges
     grid = ux.Grid.from_face_vertices(vertices_lonlat, latlon=True)
-    faces_edges_cartesian = _get_cartesian_face_edge_nodes(
-        grid.face_node_connectivity.values,
-        grid.n_face,
-        grid.n_max_face_edges,
-        grid.node_x.values,
-        grid.node_y.values,
-        grid.node_z.values,
-    )
+    faces_edges_cartesian = _get_cartesian_faces_edge_nodes(grid.face_node_connectivity.values, grid.n_face,
+                                                            grid.n_max_face_edges, grid.node_x.values,
+                                                            grid.node_y.values, grid.node_z.values)
 
     assert point_in_face(faces_edges_cartesian[0], point_xyz=point, inclusive=True)
 
@@ -1357,14 +1322,9 @@ def test_face_normal_face():
 
     # Create the grid and face edges
     grid = ux.Grid.from_face_vertices(vertices_lonlat, latlon=True)
-    faces_edges_cartesian = _get_cartesian_face_edge_nodes(
-        grid.face_node_connectivity.values,
-        grid.n_face,
-        grid.n_max_face_edges,
-        grid.node_x.values,
-        grid.node_y.values,
-        grid.node_z.values,
-    )
+    faces_edges_cartesian = _get_cartesian_faces_edge_nodes(grid.face_node_connectivity.values, grid.n_face,
+                                                            grid.n_max_face_edges, grid.node_x.values,
+                                                            grid.node_y.values, grid.node_z.values)
 
     assert point_in_face(faces_edges_cartesian[0], point_xyz=point, inclusive=True)
 
diff --git a/test/test_helpers.py b/test/test_helpers.py
index 76e352487..e8a4adbcd 100644
--- a/test/test_helpers.py
+++ b/test/test_helpers.py
@@ -11,7 +11,7 @@
 from uxarray.constants import INT_DTYPE, INT_FILL_VALUE
 from uxarray.grid.coordinates import _lonlat_rad_to_xyz, _normalize_xyz, _xyz_to_lonlat_rad
 from uxarray.grid.arcs import point_within_gca, _angle_of_2_vectors, in_between
-from uxarray.grid.utils import _get_cartesian_face_edge_nodes, _get_lonlat_rad_face_edge_nodes
+from uxarray.grid.utils import _get_cartesian_faces_edge_nodes, _get_lonlat_rad_faces_edge_nodes
 from uxarray.grid.geometry import pole_point_inside_polygon, _pole_point_inside_polygon_cartesian
 
 try:
@@ -273,9 +273,8 @@ def test_get_cartesian_face_edge_nodes_pipeline():
     node_y = grid.node_y.values
     node_z = grid.node_z.values
 
-    face_edges_connectivity_cartesian = _get_cartesian_face_edge_nodes(
-        face_node_conn, n_face, n_max_face_edges, node_x, node_y, node_z
-    )
+    face_edges_connectivity_cartesian = _get_cartesian_faces_edge_nodes(face_node_conn, n_face, n_max_face_edges,
+                                                                        node_x, node_y, node_z)
 
     result = _pole_point_inside_polygon_cartesian(
         'North', face_edges_connectivity_cartesian[0]
@@ -298,9 +297,8 @@ def test_get_cartesian_face_edge_nodes_filled_value():
     node_y = grid.node_y.values
     node_z = grid.node_z.values
 
-    face_edges_connectivity_cartesian = _get_cartesian_face_edge_nodes(
-        face_node_conn, n_face, n_max_face_edges, node_x, node_y, node_z
-    )
+    face_edges_connectivity_cartesian = _get_cartesian_faces_edge_nodes(face_node_conn, n_face, n_max_face_edges,
+                                                                        node_x, node_y, node_z)
 
     result = _pole_point_inside_polygon_cartesian(
         'North', face_edges_connectivity_cartesian[0]
@@ -335,9 +333,8 @@ def test_get_cartesian_face_edge_nodes_filled_value2():
     node_y = np.array([v0_cart[1],v1_cart[1],v2_cart[1],v3_cart[1],v4_cart[1]])
     node_z = np.array([v0_cart[2],v1_cart[2],v2_cart[2],v3_cart[2],v4_cart[2]])
 
-    face_edges_connectivity_cartesian = _get_cartesian_face_edge_nodes(
-        face_node_conn, n_face, n_max_face_edges, node_x, node_y, node_z
-    )
+    face_edges_connectivity_cartesian = _get_cartesian_faces_edge_nodes(face_node_conn, n_face, n_max_face_edges,
+                                                                        node_x, node_y, node_z)
 
     correct_result = np.array([
         [
@@ -368,9 +365,8 @@ def test_get_lonlat_face_edge_nodes_pipeline():
     node_lon = grid.node_lon.values
     node_lat = grid.node_lat.values
 
-    face_edges_connectivity_lonlat = _get_lonlat_rad_face_edge_nodes(
-        face_node_conn, n_face, n_max_face_edges, node_lon, node_lat
-    )
+    face_edges_connectivity_lonlat = _get_lonlat_rad_faces_edge_nodes(face_node_conn, n_face, n_max_face_edges,
+                                                                      node_lon, node_lat)
 
     face_edges_connectivity_lonlat = face_edges_connectivity_lonlat[0]
     face_edges_connectivity_cartesian = []
@@ -398,9 +394,8 @@ def test_get_lonlat_face_edge_nodes_filled_value():
     node_lon = grid.node_lon.values
     node_lat = grid.node_lat.values
 
-    face_edges_connectivity_lonlat = _get_lonlat_rad_face_edge_nodes(
-        face_node_conn, n_face, n_max_face_edges, node_lon, node_lat
-    )
+    face_edges_connectivity_lonlat = _get_lonlat_rad_faces_edge_nodes(face_node_conn, n_face, n_max_face_edges,
+                                                                      node_lon, node_lat)
 
     face_edges_connectivity_lonlat = face_edges_connectivity_lonlat[0]
     face_edges_connectivity_cartesian = []
@@ -434,9 +429,8 @@ def test_get_lonlat_face_edge_nodes_filled_value2():
     node_lon = np.array([v0_rad[0],v1_rad[0],v2_rad[0],v3_rad[0],v4_rad[0]])
     node_lat = np.array([v0_rad[1],v1_rad[1],v2_rad[1],v3_rad[1],v4_rad[1]])
 
-    face_edges_connectivity_lonlat = _get_lonlat_rad_face_edge_nodes(
-        face_node_conn, n_face, n_max_face_edges, node_lon, node_lat
-    )
+    face_edges_connectivity_lonlat = _get_lonlat_rad_faces_edge_nodes(face_node_conn, n_face, n_max_face_edges,
+                                                                      node_lon, node_lat)
 
     correct_result = np.array([
         [
diff --git a/test/test_integrate.py b/test/test_integrate.py
index ad9881c04..2f17c39ae 100644
--- a/test/test_integrate.py
+++ b/test/test_integrate.py
@@ -19,7 +19,7 @@
     _get_faces_constLat_intersection_info, _zonal_face_weights, \
     _zonal_face_weights_robust
 
-from uxarray.grid.utils import _get_cartesian_face_edge_nodes
+from uxarray.grid.utils import _get_cartesian_faces_edge_nodes
 
 current_path = Path(os.path.dirname(os.path.realpath(__file__)))
 
@@ -67,9 +67,9 @@ def test_get_faces_constLat_intersection_info_one_intersection():
     latitude_cart = -0.8660254037844386
     is_latlonface = False
     is_GCA_list = None
-    unique_intersections, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart,
+    unique_intersections_cart,unique_intersections_lonlat, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart,
                                                                                          is_GCA_list, is_latlonface)
-    assert len(unique_intersections) == 1
+    assert len(unique_intersections_cart) == 1
 
 
 def test_get_faces_constLat_intersection_info_encompass_pole():
@@ -97,9 +97,9 @@ 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)
-    assert len(unique_intersections) <= 2 * len(face_edges_cart)
+    unique_intersections_cart,unique_intersections_lonlat, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart,
+                                                                                              is_GCA_list, is_latlonface)
+    assert len(unique_intersections_cart) <= 2 * len(face_edges_cart)
 
 
 def test_get_faces_constLat_intersection_info_on_pole():
@@ -119,9 +119,9 @@ def test_get_faces_constLat_intersection_info_on_pole():
     latitude_cart = -0.9998476951563913
     is_latlonface = False
     is_GCA_list = None
-    unique_intersections, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart,
+    unique_intersections_cart,unique_intersections_lonlat, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart,
                                                                                          is_GCA_list, is_latlonface)
-    assert len(unique_intersections) == 2
+    assert len(unique_intersections_cart) == 2
 
 
 def test_get_faces_constLat_intersection_info_near_pole():
@@ -141,9 +141,9 @@ def test_get_faces_constLat_intersection_info_near_pole():
     latitude_deg = np.rad2deg(latitude_rad)
     is_latlonface = False
     is_GCA_list = None
-    unique_intersections, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart,
+    unique_intersections_cart,unique_intersections_lonlat, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart,
                                                                                          is_GCA_list, is_latlonface)
-    assert len(unique_intersections) == 1
+    assert len(unique_intersections_cart) == 1
 
 
 def test_get_zonal_face_interval():
@@ -1037,7 +1037,7 @@ def test_compare_zonal_weights():
     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(
+        face_edge_nodes_xyz =  _get_cartesian_faces_edge_nodes(
                 uxgrid.face_node_connectivity.values,
                 uxgrid.n_face,
                 uxgrid.n_max_face_edges,
diff --git a/test/test_intersections.py b/test/test_intersections.py
index f2fbfe1ee..6d52e4d46 100644
--- a/test/test_intersections.py
+++ b/test/test_intersections.py
@@ -161,20 +161,16 @@ def test_GCA_GCA_north_pole_angled():
 
 def test_GCA_edge_intersection_count():
 
-    from uxarray.grid.utils import _get_cartesian_face_edge_nodes
+    from uxarray.grid.utils import _get_cartesian_faces_edge_nodes
 
     # Generate a normal face that is not crossing the antimeridian or the poles
     vertices_lonlat = [[29.5, 11.0], [29.5, 10.0], [30.5, 10.0], [30.5, 11.0]]
     vertices_lonlat = np.array(vertices_lonlat)
 
     grid = ux.Grid.from_face_vertices(vertices_lonlat, latlon=True)
-    face_edge_nodes_cartesian = _get_cartesian_face_edge_nodes(
-        grid.face_node_connectivity.values,
-        grid.n_face,
-        grid.n_max_face_edges,
-        grid.node_x.values,
-        grid.node_y.values,
-        grid.node_z.values)
+    face_edge_nodes_cartesian = _get_cartesian_faces_edge_nodes(grid.face_node_connectivity.values, grid.n_face,
+                                                                grid.n_max_face_edges, grid.node_x.values,
+                                                                grid.node_y.values, grid.node_z.values)
 
     face_center_xyz = np.array([grid.face_x.values[0], grid.face_y.values[0], grid.face_z.values[0]], dtype=np.float64)
     north_pole_xyz = np.array([0.0, 0.0, 1.0], dtype=np.float64)
diff --git a/uxarray/core/zonal.py b/uxarray/core/zonal.py
index e57975b87..271d3d8a6 100644
--- a/uxarray/core/zonal.py
+++ b/uxarray/core/zonal.py
@@ -3,7 +3,6 @@
 
 
 from uxarray.grid.integrate import _zonal_face_weights, _zonal_face_weights_robust
-from uxarray.grid.utils import _get_cartesian_face_edge_nodes
 
 
 def _compute_non_conservative_zonal_mean(uxda, latitudes, use_robust_weights=False):
@@ -18,23 +17,15 @@ def _compute_non_conservative_zonal_mean(uxda, latitudes, use_robust_weights=Fal
         # Create a NumPy array for storing results
         result = np.zeros(shape, dtype=uxda.dtype)
 
-    faces_edge_nodes_xyz = _get_cartesian_face_edge_nodes(
-        uxgrid.face_node_connectivity.values,
-        uxgrid.n_face,
-        uxgrid.n_max_face_nodes,
-        uxgrid.node_x.values,
-        uxgrid.node_y.values,
-        uxgrid.node_z.values,
-    )
-
     bounds = uxgrid.bounds.values
+    faces_edges_cartesian = uxgrid.faces_edges_cartesian.values
 
     for i, lat in enumerate(latitudes):
         face_indices = uxda.uxgrid.get_faces_at_constant_latitude(lat)
 
         z = np.sin(np.deg2rad(lat))
 
-        faces_edge_nodes_xyz_candidate = faces_edge_nodes_xyz[face_indices, :, :, :]
+        faces_edge_nodes_xyz_candidate = faces_edges_cartesian[face_indices, :, :, :]
 
         n_nodes_per_face_candidate = n_nodes_per_face[face_indices]
 
diff --git a/uxarray/grid/arcs.py b/uxarray/grid/arcs.py
index dc53a2bd7..8bb12a04a 100644
--- a/uxarray/grid/arcs.py
+++ b/uxarray/grid/arcs.py
@@ -90,6 +90,55 @@ def point_within_gca(pt_xyz, gca_a_xyz, gca_b_xyz):
         return False
 
 
+@njit(cache=True)
+def _intersection_within_interval(pt_xyz, gca_a_xyz, gca_b_xyz):
+    """
+    Check if a an intersection point, which is certain to on the same plane as the CGA, lies on a given Great Circle Arc (GCA) interval, considering the smaller arc of the circle.
+    Handles the anti-meridian case as well.
+
+    Parameters
+    ----------
+    pt_xyz : numpy.ndarray
+        Cartesian coordinates of the point.
+    gca_a_xyz : numpy.ndarray
+        Cartesian coordinates of the first endpoint of the Great Circle Arc.
+    gca_b_xyz : numpy.ndarray
+        Cartesian coordinates of the second endpoint of the Great Circle Arc.
+
+    Returns
+    -------
+    bool
+        True if the point lies within the specified GCA interval, False otherwise.
+
+    Raises
+    ------
+    ValueError
+        In such cases, consider breaking the GCA into two separate arcs.
+
+    Notes
+    -----
+    - The function assumes that the point lies on the same plane as the GCA and it also assumes that the input will not be exactly 180 degree.
+    - It assumes the input represents the smaller arc of the Great Circle.
+    - The `_angle_of_2_vectors` and `_xyz_to_lonlat_rad_scalar` functions are used for calculations.
+    """
+
+    # Check if the point lies within the Great Circle Arc interval
+    pt_a = gca_a_xyz - pt_xyz
+    pt_b = gca_b_xyz - pt_xyz
+
+    # Use the dot product to determine the sign of the angle between pt_a and pt_b
+    cos_theta = np.dot(pt_a, pt_b)
+
+    # Return True if the point lies within the interval (smaller arc)
+    if cos_theta < 0:
+        return True
+    elif isclose(cos_theta, 0.0, atol=MACHINE_EPSILON):
+        # set error tolerance to 0.0
+        return True
+    else:
+        return False
+
+
 @njit(cache=True)
 def in_between(p, q, r) -> bool:
     """Determines whether the number q is between p and r.
diff --git a/uxarray/grid/geometry.py b/uxarray/grid/geometry.py
index 1c5660da6..e7ae36403 100644
--- a/uxarray/grid/geometry.py
+++ b/uxarray/grid/geometry.py
@@ -29,8 +29,8 @@
     gca_gca_intersection,
 )
 from uxarray.grid.utils import (
-    _get_cartesian_face_edge_nodes,
-    _get_lonlat_rad_face_edge_nodes,
+    _get_cartesian_faces_edge_nodes,
+    _get_lonlat_rad_faces_edge_nodes,
 )
 from uxarray.utils.computing import allclose, isclose
 
@@ -688,6 +688,7 @@ def _convert_shells_to_polygons(shells):
     return polygons
 
 
+# TODO: If we don't actually use this one, remove it
 def _pole_point_inside_polygon_cartesian(pole, face_edges_xyz):
     if isinstance(pole, str):
         pole = POLE_NAME_TO_INT[pole]
@@ -1346,6 +1347,59 @@ def compute_temp_latlon_array(
     return temp_latlon_array
 
 
+def _populate_faces_edges_cartesian(grid, return_array=False):
+    faces_edges_cartesian = _get_cartesian_faces_edge_nodes(
+        grid.face_node_connectivity.values,
+        grid.n_face,
+        grid.n_max_face_edges,
+        grid.node_x.values,
+        grid.node_y.values,
+        grid.node_z.values,
+    )
+
+    faces_edges_cartesian_xarray = xr.DataArray(
+        faces_edges_cartesian,
+        dims=["n_face", "n_max_face_edges", "two", "three"],
+        attrs={
+            "cf_role": "faces_edges_cartesian",
+            "_FillValue": INT_FILL_VALUE,
+            "long_name": "Provide the Cartesian coordinates of the edge for each face",
+            "start_index": INT_DTYPE(0),
+        },
+    )
+
+    if return_array:
+        return faces_edges_cartesian_xarray
+    else:
+        grid._ds["faces_edges_cartesian"] = faces_edges_cartesian_xarray
+
+
+def _populate_faces_edges_spherical(grid, return_array=False):
+    faces_edges_lonlat_rad = _get_lonlat_rad_faces_edge_nodes(
+        grid.face_node_connectivity.values,
+        grid.n_face,
+        grid.n_max_face_edges,
+        grid.node_lon.values,
+        grid.node_lat.values,
+    )
+
+    faces_edges_lonlat_rad_xarray = xr.DataArray(
+        faces_edges_lonlat_rad,
+        dims=["n_face", "n_max_face_edges", "two", "two"],
+        attrs={
+            "cf_role": "faces_edges_spherical",
+            "_FillValue": INT_FILL_VALUE,
+            "long_name": "Provide the Spherical coordinates (lon lat) in radians of the edge for each face",
+            "start_index": INT_DTYPE(0),
+        },
+    )
+
+    if return_array:
+        return faces_edges_lonlat_rad_xarray
+    else:
+        grid._ds["faces_edges_spherical"] = faces_edges_lonlat_rad_xarray
+
+
 def _populate_bounds(
     grid, is_latlonface: bool = False, is_face_GCA_list=None, return_array=False
 ):
@@ -1357,6 +1411,10 @@ def _populate_bounds(
 
     Parameters
     ----------
+    grid : object
+        The grid object whose internal representation will be updated with the computed
+        face bounds.
+
     is_latlonface : bool, optional
         A global flag that indicates if faces are latlon faces. If True, all faces
         are treated as latlon faces, meaning that all edges are either longitude or
@@ -1407,25 +1465,14 @@ def _populate_bounds(
     """
 
     # Ensure grid's cartesian coordinates are normalized
+    # TODO: Is it possible to have a more flexible normalization? (Avoid duplicated normalizations)
     grid.normalize_cartesian_coordinates()
 
-    # Prepare data for Numba functions
-    faces_edges_cartesian = _get_cartesian_face_edge_nodes(
-        grid.face_node_connectivity.values,
-        grid.n_face,
-        grid.n_max_face_edges,
-        grid.node_x.values,
-        grid.node_y.values,
-        grid.node_z.values,
-    )
+    # If the face_edges_cartesian and face_edges_lonlat_rad didn't exist,
+    # The following call will automatically populate them
+    faces_edges_cartesian = grid.faces_edges_cartesian.values
 
-    faces_edges_lonlat_rad = _get_lonlat_rad_face_edge_nodes(
-        grid.face_node_connectivity.values,
-        grid.n_face,
-        grid.n_max_face_edges,
-        grid.node_lon.values,
-        grid.node_lat.values,
-    )
+    faces_edges_lonlat_rad = grid.faces_edges_spherical.values
 
     n_nodes_per_face = grid.n_nodes_per_face.values
 
diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py
index 9cae36bac..3d0acaa0a 100644
--- a/uxarray/grid/grid.py
+++ b/uxarray/grid/grid.py
@@ -14,7 +14,7 @@
     Tuple,
 )
 
-from uxarray.grid.utils import _get_cartesian_face_edge_nodes
+from uxarray.grid.utils import _get_cartesian_faces_edge_nodes
 
 # reader and writer imports
 from uxarray.io._exodus import _read_exodus, _encode_exodus
@@ -71,6 +71,8 @@
     _grid_to_matplotlib_polycollection,
     _grid_to_matplotlib_linecollection,
     _populate_bounds,
+    _populate_faces_edges_cartesian,
+    _populate_faces_edges_spherical,
     _construct_boundary_edge_indices,
     compute_temp_latlon_array,
     _find_faces,
@@ -1443,6 +1445,26 @@ def face_areas(self, value):
         assert isinstance(value, xr.DataArray)
         self._ds["face_areas"] = value
 
+    @property
+    def faces_edges_cartesian(self):
+        """Cartesian Coordinates for each Face Edge.
+
+        Dimensions ``(n_face, n_max_face_edges, 2, 3)``
+        """
+        if "faces_edges_cartesian" not in self._ds:
+            _populate_faces_edges_cartesian(self)
+        return self._ds["faces_edges_cartesian"]
+
+    @property
+    def faces_edges_spherical(self):
+        """Latitude Longitude Coordinates for each Face in radians.
+
+        Dimensions ``(n_face, n_max_face_edges, 2, 2)``
+        """
+        if "faces_edges_spherical" not in self._ds:
+            _populate_faces_edges_spherical(self)
+        return self._ds["faces_edges_spherical"]
+
     @property
     def bounds(self):
         """Latitude Longitude Bounds for each Face in radians.
@@ -1456,6 +1478,10 @@ def bounds(self):
                     "This initial execution will be significantly longer.",
                     RuntimeWarning,
                 )
+            if "faces_edges_cartesian" not in self._ds:
+                _populate_faces_edges_cartesian(self)
+            if "faces_edges_spherical" not in self._ds:
+                _populate_faces_edges_spherical(self)
             _populate_bounds(self)
         return self._ds["bounds"]
 
@@ -2638,7 +2664,8 @@ def get_faces_containing_point(
             return np.empty(0, dtype=np.int64)
 
         # Get the faces in terms of their edges
-        face_edge_nodes_xyz = _get_cartesian_face_edge_nodes(
+        # Since this is a new subset, the cartesian_faces_edge_nodes need to be recalculated anyway
+        face_edge_nodes_xyz = _get_cartesian_faces_edge_nodes(
             subset.face_node_connectivity.values,
             subset.n_face,
             subset.n_max_face_nodes,
diff --git a/uxarray/grid/integrate.py b/uxarray/grid/integrate.py
index d0beade57..4ac4dba11 100644
--- a/uxarray/grid/integrate.py
+++ b/uxarray/grid/integrate.py
@@ -215,20 +215,21 @@ def _get_zonal_face_interval(
     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
-            )
+        (
+            unique_intersections_cart,
+            unique_intersection_lonlat,
+            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:
+        if len(unique_intersections_cart) == 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]
-        )
+        # Extract unique longitudes from intersection points
+        longitudes = np.array([pt[0] for pt in unique_intersection_lonlat])
 
         # Handle special wrap-around cases (crossing anti-meridian, etc.)
         if face_lon_bound_left >= face_lon_bound_right or (
@@ -378,7 +379,8 @@ def _get_faces_constLat_intersection_info(
     -------
     tuple
         A tuple containing:
-        - intersections_pts_list_cart (list): A list of intersection points where each point is where an edge intersects with the latitude.
+        - intersections_pts_list_cart (list): A list of intersection points in cartesian where each point is where an edge intersects with the latitude.
+        - intersections_pts_list_lon (list): A list of intersection points in lonlat where each point is where an edge intersects with the latitude.
         - pt_lon_min (float): The min longnitude of the interseted intercal in radian if any; otherwise, None..
         - pt_lon_max (float): The max longnitude of the interseted intercal in radian, if any; otherwise, None.
     """
@@ -418,28 +420,50 @@ def _get_faces_constLat_intersection_info(
                     intersections_pts_list_cart.append(intersections[0])
 
     # Find the unique intersection points
-    unique_intersections = np.unique(intersections_pts_list_cart, axis=0)
+    unique_intersections_cart = np.unique(intersections_pts_list_cart, axis=0)
 
-    if len(unique_intersections) == 2:
+    if len(unique_intersections_cart) == 2:
         unique_intersection_lonlat = np.array(
-            [_xyz_to_lonlat_rad(pt[0], pt[1], pt[2]) for pt in unique_intersections]
+            [
+                _xyz_to_lonlat_rad(pt[0], pt[1], pt[2])
+                for pt in unique_intersections_cart
+            ]
         )
 
         sorted_lonlat = np.sort(unique_intersection_lonlat, axis=0)
         pt_lon_min, pt_lon_max = sorted_lonlat[:, 0]
-        return unique_intersections, pt_lon_min, pt_lon_max
-    elif len(unique_intersections) == 1:
-        return unique_intersections, None, None
-    elif len(unique_intersections) != 0 and len(unique_intersections) != 1:
+        return (
+            unique_intersections_cart,
+            unique_intersection_lonlat,
+            pt_lon_min,
+            pt_lon_max,
+        )
+    elif len(unique_intersections_cart) == 1:
+        return (
+            unique_intersections_cart,
+            np.array(
+                _xyz_to_lonlat_rad(
+                    unique_intersections_cart[0][0],
+                    unique_intersections_cart[0][1],
+                    unique_intersections_cart[0][2],
+                )
+            ),
+            None,
+            None,
+        )
+    elif len(unique_intersections_cart) != 0 and len(unique_intersections_cart) != 1:
         # If the unique intersections numbers is larger than n_edges * 2, then it means the face is concave
-        if len(unique_intersections) > len(valid_edges) * 2:
+        if len(unique_intersections_cart) > len(valid_edges) * 2:
             raise ValueError(
                 "UXarray doesn't support concave face with intersections points as currently, please modify your grids accordingly"
             )
         else:
             # Now return all the intersections points and the pt_lon_min, pt_lon_max
             unique_intersection_lonlat = np.array(
-                [_xyz_to_lonlat_rad(pt[0], pt[1], pt[2]) for pt in unique_intersections]
+                [
+                    _xyz_to_lonlat_rad(pt[0], pt[1], pt[2])
+                    for pt in unique_intersections_cart
+                ]
             )
 
             sorted_lonlat = np.sort(unique_intersection_lonlat, axis=0)
@@ -449,8 +473,13 @@ def _get_faces_constLat_intersection_info(
                 np.max(sorted_lonlat[:, 0]),
             )
 
-            return unique_intersections, pt_lon_min, pt_lon_max
-    elif len(unique_intersections) == 0:
+            return (
+                unique_intersections_cart,
+                unique_intersection_lonlat,
+                pt_lon_min,
+                pt_lon_max,
+            )
+    elif len(unique_intersections_cart) == 0:
         raise ValueError(
             "No intersections are found for the face, please make sure the build_latlon_box generates the correct results"
         )
diff --git a/uxarray/grid/intersections.py b/uxarray/grid/intersections.py
index 0eabb7b26..dccf0ef06 100644
--- a/uxarray/grid/intersections.py
+++ b/uxarray/grid/intersections.py
@@ -7,6 +7,7 @@
     in_between,
     extreme_gca_z,
     point_within_gca,
+    _intersection_within_interval,
 )
 from uxarray.utils.computing import allclose, cross, norm
 
@@ -347,15 +348,15 @@ def gca_gca_intersection(gca_a_xyz, gca_b_xyz):
     x2_xyz = -x1_xyz
 
     # Check intersection points
-    if point_within_gca(x1_xyz, w0_xyz, w1_xyz) and point_within_gca(
-        x1_xyz, v0_xyz, v1_xyz
-    ):
+    if _intersection_within_interval(
+        x1_xyz, w0_xyz, w1_xyz
+    ) and _intersection_within_interval(x1_xyz, v0_xyz, v1_xyz):
         res[count, :] = x1_xyz
         count += 1
 
-    if point_within_gca(x2_xyz, w0_xyz, w1_xyz) and point_within_gca(
-        x2_xyz, v0_xyz, v1_xyz
-    ):
+    if _intersection_within_interval(
+        x2_xyz, w0_xyz, w1_xyz
+    ) and _intersection_within_interval(x2_xyz, v0_xyz, v1_xyz):
         res[count, :] = x2_xyz
         count += 1
 
@@ -417,17 +418,20 @@ def gca_const_lat_intersection(gca_cart, const_z):
 
     nx, ny, nz = n
 
-    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)
+    nx_square_ny_square = nx**2 + ny**2
+    n_square = nx_square_ny_square + nz**2
+
+    s_tilde = np.sqrt(nx_square_ny_square - n_square * const_z**2)
+    p1_x = -(1.0 / (nx_square_ny_square)) * (const_z * nx * nz + s_tilde * ny)
+    p2_x = -(1.0 / (nx_square_ny_square)) * (const_z * nx * nz - s_tilde * ny)
+    p1_y = -(1.0 / (nx_square_ny_square)) * (const_z * ny * nz - s_tilde * nx)
+    p2_y = -(1.0 / (nx_square_ny_square)) * (const_z * ny * nz + s_tilde * nx)
 
     p1 = np.array([p1_x, p1_y, const_z])
     p2 = np.array([p2_x, p2_y, const_z])
 
-    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])
+    p1_intersects_gca = _intersection_within_interval(p1, gca_cart[0], gca_cart[1])
+    p2_intersects_gca = _intersection_within_interval(p2, gca_cart[0], gca_cart[1])
 
     if p1_intersects_gca and p2_intersects_gca:
         res[0] = p1
diff --git a/uxarray/grid/utils.py b/uxarray/grid/utils.py
index d9a6621c9..13a053081 100644
--- a/uxarray/grid/utils.py
+++ b/uxarray/grid/utils.py
@@ -138,7 +138,7 @@ def _swap_first_fill_value_with_last(arr):
     return arr
 
 
-def _get_cartesian_face_edge_nodes(
+def _get_cartesian_faces_edge_nodes(
     face_node_conn, n_face, n_max_face_edges, node_x, node_y, node_z
 ):
     """Construct an array to hold the edge Cartesian coordinates connectivity
@@ -179,7 +179,7 @@ def _get_cartesian_face_edge_nodes(
     >>> node_x = np.array([0, 1, 1, 0, 1, 0])
     >>> node_y = np.array([0, 0, 1, 1, 2, 2])
     >>> node_z = np.array([0, 0, 0, 0, 1, 1])
-    >>> _get_cartesian_face_edge_nodes(
+    >>> _get_cartesian_faces_edge_nodes(
     ...     face_node_conn, n_face, n_max_face_edges, node_x, node_y, node_z
     ... )
     array([[[[    0,     0,     0],
@@ -259,7 +259,7 @@ def _get_cartesian_face_edge_nodes(
     return face_edges_cartesian.reshape(n_face, n_max_face_edges, 2, 3)
 
 
-def _get_lonlat_rad_face_edge_nodes(
+def _get_lonlat_rad_faces_edge_nodes(
     face_node_conn, n_face, n_max_face_edges, node_lon, node_lat
 ):
     """Construct an array to hold the edge latitude and longitude in radians