diff --git a/test/test_cross_sections.py b/test/test_cross_sections.py
index 55068b6c5..1847e69b3 100644
--- a/test/test_cross_sections.py
+++ b/test/test_cross_sections.py
@@ -1,5 +1,6 @@
 import uxarray as ux
 import pytest
+import numpy as np
 from pathlib import Path
 import os
 
@@ -10,7 +11,9 @@
 quad_hex_grid_path = current_path / 'meshfiles' / "ugrid" / "quad-hexagon" / 'grid.nc'
 quad_hex_data_path = current_path / 'meshfiles' / "ugrid" / "quad-hexagon" / 'data.nc'
 
-cube_sphere_grid = current_path / "meshfiles" / "geos-cs" / "c12" / "test-c12.native.nc4"
+cube_sphere_grid = current_path / "meshfiles" / "ugrid" / "outCSne30" / "outCSne30.ug"
+
+from uxarray.grid.intersections import constant_lat_intersections_face_bounds
 
 
 
@@ -32,33 +35,38 @@ class TestQuadHex:
     All four faces intersect a constant latitude of 0.0
     """
 
-    def test_constant_lat_cross_section_grid(self):
+    @pytest.mark.parametrize("use_spherical_bounding_box", [True, False])
+    def test_constant_lat_cross_section_grid(self, use_spherical_bounding_box):
+
+
+
         uxgrid = ux.open_grid(quad_hex_grid_path)
 
-        grid_top_two = uxgrid.cross_section.constant_latitude(lat=0.1)
+        grid_top_two = uxgrid.cross_section.constant_latitude(lat=0.1, use_spherical_bounding_box=use_spherical_bounding_box)
 
         assert grid_top_two.n_face == 2
 
-        grid_bottom_two = uxgrid.cross_section.constant_latitude(lat=-0.1)
+        grid_bottom_two = uxgrid.cross_section.constant_latitude(lat=-0.1, use_spherical_bounding_box=use_spherical_bounding_box)
 
         assert grid_bottom_two.n_face == 2
 
-        grid_all_four = uxgrid.cross_section.constant_latitude(lat=0.0)
+        grid_all_four = uxgrid.cross_section.constant_latitude(lat=0.0, use_spherical_bounding_box=use_spherical_bounding_box)
 
         assert grid_all_four.n_face == 4
 
         with pytest.raises(ValueError):
             # no intersections found at this line
-            uxgrid.cross_section.constant_latitude(lat=10.0)
+            uxgrid.cross_section.constant_latitude(lat=10.0, use_spherical_bounding_box=use_spherical_bounding_box)
 
-    def test_constant_lon_cross_section_grid(self):
+    @pytest.mark.parametrize("use_spherical_bounding_box", [False])
+    def test_constant_lon_cross_section_grid(self, use_spherical_bounding_box):
         uxgrid = ux.open_grid(quad_hex_grid_path)
 
-        grid_left_two = uxgrid.cross_section.constant_longitude(lon=-0.1)
+        grid_left_two = uxgrid.cross_section.constant_longitude(lon=-0.1, use_spherical_bounding_box=use_spherical_bounding_box)
 
         assert grid_left_two.n_face == 2
 
-        grid_right_two = uxgrid.cross_section.constant_longitude(lon=0.2)
+        grid_right_two = uxgrid.cross_section.constant_longitude(lon=0.2, use_spherical_bounding_box=use_spherical_bounding_box)
 
         assert grid_right_two.n_face == 2
 
@@ -66,59 +74,110 @@ def test_constant_lon_cross_section_grid(self):
             # no intersections found at this line
             uxgrid.cross_section.constant_longitude(lon=10.0)
 
-
-    def test_constant_lat_cross_section_uxds(self):
+    @pytest.mark.parametrize("use_spherical_bounding_box", [False])
+    def test_constant_lat_cross_section_uxds(self, use_spherical_bounding_box):
         uxds = ux.open_dataset(quad_hex_grid_path, quad_hex_data_path)
+        uxds.uxgrid.normalize_cartesian_coordinates()
 
-        da_top_two = uxds['t2m'].cross_section.constant_latitude(lat=0.1)
+        da_top_two = uxds['t2m'].cross_section.constant_latitude(lat=0.1, use_spherical_bounding_box=use_spherical_bounding_box)
 
         nt.assert_array_equal(da_top_two.data, uxds['t2m'].isel(n_face=[1, 2]).data)
 
-        da_bottom_two = uxds['t2m'].cross_section.constant_latitude(lat=-0.1)
+        da_bottom_two = uxds['t2m'].cross_section.constant_latitude(lat=-0.1, use_spherical_bounding_box=use_spherical_bounding_box)
 
         nt.assert_array_equal(da_bottom_two.data, uxds['t2m'].isel(n_face=[0, 3]).data)
 
-        da_all_four = uxds['t2m'].cross_section.constant_latitude(lat=0.0)
+        da_all_four = uxds['t2m'].cross_section.constant_latitude(lat=0.0, use_spherical_bounding_box=use_spherical_bounding_box)
 
         nt.assert_array_equal(da_all_four.data , uxds['t2m'].data)
 
         with pytest.raises(ValueError):
             # no intersections found at this line
-            uxds['t2m'].cross_section.constant_latitude(lat=10.0)
+            uxds['t2m'].cross_section.constant_latitude(lat=10.0, use_spherical_bounding_box=use_spherical_bounding_box)
 
-    def test_constant_lon_cross_section_uxds(self):
+    @pytest.mark.parametrize("use_spherical_bounding_box", [False])
+    def test_constant_lon_cross_section_uxds(self, use_spherical_bounding_box):
         uxds = ux.open_dataset(quad_hex_grid_path, quad_hex_data_path)
+        uxds.uxgrid.normalize_cartesian_coordinates()
 
-        da_left_two = uxds['t2m'].cross_section.constant_longitude(lon=-0.1)
+        da_left_two = uxds['t2m'].cross_section.constant_longitude(lon=-0.1, use_spherical_bounding_box=use_spherical_bounding_box)
 
         nt.assert_array_equal(da_left_two.data, uxds['t2m'].isel(n_face=[0, 2]).data)
 
-        da_right_two = uxds['t2m'].cross_section.constant_longitude(lon=0.2)
+        da_right_two = uxds['t2m'].cross_section.constant_longitude(lon=0.2, use_spherical_bounding_box=use_spherical_bounding_box)
 
         nt.assert_array_equal(da_right_two.data, uxds['t2m'].isel(n_face=[1, 3]).data)
 
         with pytest.raises(ValueError):
             # no intersections found at this line
-            uxds['t2m'].cross_section.constant_longitude(lon=10.0)
+            uxds['t2m'].cross_section.constant_longitude(lon=10.0, use_spherical_bounding_box=use_spherical_bounding_box)
 
 
-class TestGeosCubeSphere:
-    def test_north_pole(self):
+class TestCubeSphere:
+    @pytest.mark.parametrize("use_spherical_bounding_box", [True, False])
+    def test_north_pole(self, use_spherical_bounding_box):
         uxgrid = ux.open_grid(cube_sphere_grid)
 
         lats = [89.85, 89.9, 89.95, 89.99]
 
         for lat in lats:
-            cross_grid = uxgrid.cross_section.constant_latitude(lat=lat)
+            cross_grid = uxgrid.cross_section.constant_latitude(lat=lat, use_spherical_bounding_box=use_spherical_bounding_box)
             # Cube sphere grid should have 4 faces centered around the pole
             assert cross_grid.n_face == 4
-
-    def test_south_pole(self):
+    @pytest.mark.parametrize("use_spherical_bounding_box", [True, False])
+    def test_south_pole(self, use_spherical_bounding_box):
         uxgrid = ux.open_grid(cube_sphere_grid)
 
         lats = [-89.85, -89.9, -89.95, -89.99]
 
         for lat in lats:
-            cross_grid = uxgrid.cross_section.constant_latitude(lat=lat)
+            cross_grid = uxgrid.cross_section.constant_latitude(lat=lat, use_spherical_bounding_box=use_spherical_bounding_box)
             # Cube sphere grid should have 4 faces centered around the pole
             assert cross_grid.n_face == 4
+
+
+
+class TestCandidateFacesUsingBounds:
+
+    def test_constant_lat(self):
+        bounds = np.array([
+            [[-45, 45], [0, 360]],
+            [[-90, -45], [0, 360]],
+            [[45, 90], [0, 360]],
+        ])
+
+        bounds_rad = np.deg2rad(bounds)
+
+        const_lat = 0
+
+        candidate_faces = constant_lat_intersections_face_bounds(
+            lat=const_lat,
+            face_min_lat_rad=bounds_rad[:, 0, 0],
+            face_max_lat_rad=bounds_rad[:, 0, 1],
+        )
+
+        # Expected output
+        expected_faces = np.array([0])
+
+        # Test the function output
+        nt.assert_array_equal(candidate_faces, expected_faces)
+
+    def test_constant_lat_out_of_bounds(self):
+
+        bounds = np.array([
+            [[-45, 45], [0, 360]],
+            [[-90, -45], [0, 360]],
+            [[45, 90], [0, 360]],
+        ])
+
+        bounds_rad = np.deg2rad(bounds)
+
+        const_lat = 100
+
+        candidate_faces = constant_lat_intersections_face_bounds(
+            lat=const_lat,
+            face_min_lat_rad=bounds_rad[:, 0, 0],
+            face_max_lat_rad=bounds_rad[:, 0, 1],
+        )
+
+        assert len(candidate_faces) == 0
diff --git a/uxarray/__init__.py b/uxarray/__init__.py
index f14c9961d..857c65f6c 100644
--- a/uxarray/__init__.py
+++ b/uxarray/__init__.py
@@ -33,6 +33,9 @@ def disable_fma():
     uxarray.constants.ENABLE_FMA = False
 
 
+disable_fma()
+
+
 __all__ = (
     "open_grid",
     "open_dataset",
diff --git a/uxarray/cross_sections/dataarray_accessor.py b/uxarray/cross_sections/dataarray_accessor.py
index 299451d01..de35739ab 100644
--- a/uxarray/cross_sections/dataarray_accessor.py
+++ b/uxarray/cross_sections/dataarray_accessor.py
@@ -21,7 +21,7 @@ def __repr__(self):
 
         return prefix + methods_heading
 
-    def constant_latitude(self, lat: float, method="fast"):
+    def constant_latitude(self, lat: float, use_spherical_bounding_box=False):
         """Extracts a cross-section of the data array at a specified constant
         latitude.
 
@@ -29,12 +29,8 @@ def constant_latitude(self, lat: float, method="fast"):
         ----------
         lat : float
             The latitude at which to extract the cross-section, in degrees.
-        method : str, optional
-            The internal method to use when identifying faces at the constant latitude.
-            Options are:
-            - 'fast': Uses a faster but potentially less accurate method for face identification.
-            - 'accurate': Uses a slower but more accurate method.
-            Default is 'fast'.
+        use_spherical_bounding_box : bool, optional
+            If True, uses a spherical bounding box for intersection calculations.
 
         Raises
         ------
@@ -44,17 +40,14 @@ def constant_latitude(self, lat: float, method="fast"):
         Examples
         --------
         >>> uxda.constant_latitude_cross_section(lat=-15.5)
-
-        Notes
-        -----
-        The accuracy and performance of the function can be controlled using the `method` parameter.
-        For higher precision requreiments, consider using method='acurate'.
         """
-        faces = self.uxda.uxgrid.get_faces_at_constant_latitude(lat, method)
+        faces = self.uxda.uxgrid.get_faces_at_constant_latitude(
+            lat, use_spherical_bounding_box
+        )
 
         return self.uxda.isel(n_face=faces)
 
-    def constant_longitude(self, lon: float, method="fast"):
+    def constant_longitude(self, lon: float, use_spherical_bounding_box=False):
         """Extracts a cross-section of the data array at a specified constant
         longitude.
 
@@ -62,12 +55,8 @@ def constant_longitude(self, lon: float, method="fast"):
         ----------
         lon : float
             The longitude at which to extract the cross-section, in degrees.
-        method : str, optional
-            The internal method to use when identifying faces at the constant longitude.
-            Options are:
-            - 'fast': Uses a faster but potentially less accurate method for face identification.
-            - 'accurate': Uses a slower but more accurate method.
-            Default is 'fast'.
+        use_spherical_bounding_box : bool, optional
+            If True, uses a spherical bounding box for intersection calculations.
 
         Raises
         ------
@@ -83,7 +72,9 @@ def constant_longitude(self, lon: float, method="fast"):
         The accuracy and performance of the function can be controlled using the `method` parameter.
         For higher precision requreiments, consider using method='acurate'.
         """
-        faces = self.uxda.uxgrid.get_faces_at_constant_longitude(lon, method)
+        faces = self.uxda.uxgrid.get_faces_at_constant_longitude(
+            lon, use_spherical_bounding_box
+        )
 
         return self.uxda.isel(n_face=faces)
 
diff --git a/uxarray/cross_sections/grid_accessor.py b/uxarray/cross_sections/grid_accessor.py
index a193ec8e0..5b9055166 100644
--- a/uxarray/cross_sections/grid_accessor.py
+++ b/uxarray/cross_sections/grid_accessor.py
@@ -20,15 +20,12 @@ def __repr__(self):
         methods_heading += "  * constant_latitude(lat, )\n"
         return prefix + methods_heading
 
-    def constant_latitude(self, lat: float, return_face_indices=False, method="fast"):
+    def constant_latitude(
+        self, lat: float, return_face_indices=False, use_spherical_bounding_box=False
+    ):
         """Extracts a cross-section of the grid at a specified constant
         latitude.
 
-        This method identifies and returns all faces (or grid elements) that intersect
-        with a given latitude. The returned cross-section can include either just the grid
-        or both the grid elements and the corresponding face indices, depending
-        on the `return_face_indices` parameter.
-
         Parameters
         ----------
         lat : float
@@ -37,12 +34,9 @@ def constant_latitude(self, lat: float, return_face_indices=False, method="fast"
             If True, returns both the grid at the specified latitude and the indices
             of the intersecting faces. If False, only the grid is returned.
             Default is False.
-        method : str, optional
-            The internal method to use when identifying faces at the constant latitude.
-            Options are:
-            - 'fast': Uses a faster but potentially less accurate method for face identification.
-            - 'accurate': Uses a slower but more accurate method.
-            Default is 'fast'.
+        use_spherical_bounding_box : bool, optional
+            If True, uses a spherical bounding box for intersection calculations.
+
 
         Returns
         -------
@@ -69,7 +63,9 @@ def constant_latitude(self, lat: float, return_face_indices=False, method="fast"
         The accuracy and performance of the function can be controlled using the `method` parameter.
         For higher precision requreiments, consider using method='acurate'.
         """
-        faces = self.uxgrid.get_faces_at_constant_latitude(lat, method)
+        faces = self.uxgrid.get_faces_at_constant_latitude(
+            lat, use_spherical_bounding_box
+        )
 
         if len(faces) == 0:
             raise ValueError(f"No intersections found at lat={lat}.")
@@ -81,29 +77,22 @@ def constant_latitude(self, lat: float, return_face_indices=False, method="fast"
         else:
             return grid_at_constant_lat
 
-    def constant_longitude(self, lon: float, return_face_indices=False, method="fast"):
+    def constant_longitude(
+        self, lon: float, use_spherical_bounding_box=False, return_face_indices=False
+    ):
         """Extracts a cross-section of the grid at a specified constant
         longitude.
 
-        This method identifies and returns all faces (or grid elements) that intersect
-        with a given longitude. The returned cross-section can include either just the grid
-        or both the grid elements and the corresponding face indices, depending
-        on the `return_face_indices` parameter.
-
         Parameters
         ----------
         lon : float
             The longitude at which to extract the cross-section, in degrees.
+        use_spherical_bounding_box : bool, optional
+            If True, uses a spherical bounding box for intersection calculations.
         return_face_indices : bool, optional
             If True, returns both the grid at the specified longitude and the indices
             of the intersecting faces. If False, only the grid is returned.
             Default is False.
-        method : str, optional
-            The internal method to use when identifying faces at the constant longitude.
-            Options are:
-            - 'fast': Uses a faster but potentially less accurate method for face identification.
-            - 'accurate': Uses a slower but more accurate method.
-            Default is 'fast'.
 
         Returns
         -------
@@ -130,10 +119,12 @@ def constant_longitude(self, lon: float, return_face_indices=False, method="fast
         The accuracy and performance of the function can be controlled using the `method` parameter.
         For higher precision requreiments, consider using method='acurate'.
         """
-        faces = self.uxgrid.get_faces_at_constant_longitude(lon, method)
+        faces = self.uxgrid.get_faces_at_constant_longitude(
+            lon, use_spherical_bounding_box
+        )
 
         if len(faces) == 0:
-            raise ValueError(f"No intersections found at lon={lon}.")
+            raise ValueError(f"No intersections found at lon={lon}")
 
         grid_at_constant_lon = self.uxgrid.isel(n_face=faces)
 
diff --git a/uxarray/grid/geometry.py b/uxarray/grid/geometry.py
index 9a6566f67..1ae95d99e 100644
--- a/uxarray/grid/geometry.py
+++ b/uxarray/grid/geometry.py
@@ -34,8 +34,7 @@
 
 
 def _unique_points(points, tolerance=ERROR_TOLERANCE):
-    """Identify unique intersection points from a list of points, considering
-    floating point precision errors.
+    """Identify unique intersection points from a list of points, considering floating point precision errors.
 
     Parameters
     ----------
diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py
index d81fa5f51..e7007af26 100644
--- a/uxarray/grid/grid.py
+++ b/uxarray/grid/grid.py
@@ -69,8 +69,9 @@
 )
 
 from uxarray.grid.intersections import (
-    fast_constant_lat_intersections,
-    fast_constant_lon_intersections,
+    constant_lat_intersections_no_extreme,
+    constant_lon_intersections_no_extreme,
+    constant_lat_intersections_face_bounds,
 )
 
 from spatialpandas import GeoDataFrame
@@ -1359,7 +1360,8 @@ def bounds(self):
         """
         if "bounds" not in self._ds:
             warn(
-                "Constructing of `Grid.bounds` has not been optimized, which may lead to a long execution time."
+                "Computing 'Grid.bounds' for the first time. This may take some time...",
+                UserWarning,
             )
             _populate_bounds(self)
 
@@ -2181,132 +2183,153 @@ def isel(self, **dim_kwargs):
                 "Indexing must be along a grid dimension: ('n_node', 'n_edge', 'n_face')"
             )
 
-    def get_edges_at_constant_latitude(self, lat, method="fast"):
-        """Identifies the edges of the grid that intersect with a specified
-        constant latitude.
-
-        This method computes the intersection of grid edges with a given latitude and
-        returns a collection of edges that cross or are aligned with that latitude.
-        The method used for identifying these edges can be controlled by the `method`
-        parameter.
+    def get_edges_at_constant_latitude(self, lat, use_spherical_bounding_box=False):
+        """Identifies the indices of edges that intersect with a line of constant latitude.
 
         Parameters
         ----------
-        lat : float
+        lon : float
             The latitude at which to identify intersecting edges, in degrees.
-        method : str, optional
-            The computational method used to determine edge intersections. Options are:
-            - 'fast': Uses a faster but potentially less accurate method for determining intersections.
-            - 'accurate': Uses a slower but more precise method.
-            Default is 'fast'.
+                use_spherical_bounding_box : bool, optional
+            If `True`,
+            computes the bounding box for each face using great circle arcs for edges
+            and considers extreme minimums or maximums to increase accuracy.
+            Defaults to `False`.
 
         Returns
         -------
-        edges : array
-            A squeezed array of edges that intersect the specified constant latitude.
+        faces : numpy.ndarray
+            An array of edge indices that intersect with the specified latitude.
         """
-        if method == "fast":
-            edges = fast_constant_lat_intersections(
-                lat, self.edge_node_z.values, self.n_edge
+
+        if lat > 90.0 or lat < -90.0:
+            raise ValueError(
+                f"Latitude must be between -90 and 90 degrees. Received {lat}"
+            )
+
+        if use_spherical_bounding_box:
+            raise NotImplementedError(
+                "Computing the intersection using the spherical bounding box"
+                "is not yet supported."
             )
-        elif method == "accurate":
-            raise NotImplementedError("Accurate method not yet implemented.")
         else:
-            raise ValueError(f"Invalid method: {method}.")
+            edges = constant_lat_intersections_no_extreme(
+                lat, self.edge_node_z.values, self.n_edge
+            )
+
         return edges.squeeze()
 
-    def get_faces_at_constant_latitude(self, lat, method="fast"):
-        """Identifies the faces of the grid that intersect with a specified
-        constant latitude.
+    def get_faces_at_constant_latitude(self, lat, use_spherical_bounding_box=False):
+        """
+        Identifies the indices of faces that intersect with a line of constant latitude.
 
-        This method finds the faces (or cells) of the grid that intersect a given latitude
-        by first identifying the intersecting edges and then determining the faces connected
-        to these edges. The method used for identifying edges can be adjusted with the `method`
-        parameter.
+        When `use_spherical_bounding_box` is set to `True`,
+        the bounding box for each face is computed by representing each edge as a great circle arc.
+        This approach takes into account the extreme minimums or maximums along the arcs.
 
         Parameters
         ----------
         lat : float
             The latitude at which to identify intersecting faces, in degrees.
-        method : str, optional
-            The computational method used to determine intersecting edges. Options are:
-            - 'fast': Uses a faster but potentially less accurate method for determining intersections.
-            - 'accurate': Uses a slower but more precise method.
-            Default is 'fast'.
+        use_spherical_bounding_box : bool, optional
+            If `True`,
+            computes the bounding box for each face using great circle arcs for edges
+            and considers extreme minimums or maximums to increase accuracy.
+            Defaults to `False`.
 
         Returns
         -------
-        faces : array
-            An array of unique face indices that intersect the specified latitude.
-            Faces that are invalid or missing (e.g., with a fill value) are excluded
-            from the result.
+        faces : numpy.ndarray
+            An array of face indices that intersect with the specified latitude.
         """
-        edges = self.get_edges_at_constant_latitude(lat, method)
-        faces = np.unique(self.edge_face_connectivity[edges].data.ravel())
 
-        return faces[faces != INT_FILL_VALUE]
+        if lat > 90.0 or lat < -90.0:
+            raise ValueError(
+                f"Latitude must be between -90 and 90 degrees. Received {lat}"
+            )
+
+        if use_spherical_bounding_box:
+            faces = constant_lat_intersections_face_bounds(
+                lat=lat,
+                face_min_lat_rad=self.bounds.values[:, 0, 0],
+                face_max_lat_rad=self.bounds.values[:, 0, 1],
+            )
+            return faces
+        else:
+            edges = self.get_edges_at_constant_latitude(lat, use_spherical_bounding_box)
+            faces = np.unique(self.edge_face_connectivity[edges].data.ravel())
 
-    def get_edges_at_constant_longitude(self, lon, method="fast"):
-        """Identifies the edges of the grid that intersect with a specified
-        constant longitude.
+            return faces[faces != INT_FILL_VALUE]
 
-        This method computes the intersection of grid edges with a given longitude and
-        returns a collection of edges that cross or are aligned with that longitude.
-        The method used for identifying these edges can be controlled by the `method`
-        parameter.
+    def get_edges_at_constant_longitude(self, lon, use_spherical_bounding_box=False):
+        """
+        Identifies the indices of edges that intersect with a line of constant longitude.
 
         Parameters
         ----------
         lon : float
             The longitude at which to identify intersecting edges, in degrees.
-        method : str, optional
-            The computational method used to determine edge intersections. Options are:
-            - 'fast': Uses a faster but potentially less accurate method for determining intersections.
-            - 'accurate': Uses a slower but more precise method.
-            Default is 'fast'.
+        use_spherical_bounding_box : bool, optional
+            If `True`,
+            computes the bounding box for each face using great circle arcs for edges
+            and considers extreme minimums or maximums to increase accuracy.
+            Defaults to `False`.
 
         Returns
         -------
-        edges : array
-            A squeezed array of edges that intersect the specified constant longitude.
+        faces : numpy.ndarray
+            An array of edge indices that intersect with the specified longitude.
         """
-        if method == "fast":
-            edges = fast_constant_lon_intersections(
-                lon, self.edge_node_x.values, self.edge_node_y.values, self.n_edge
+
+        if lon > 180.0 or lon < -180.0:
+            raise ValueError(
+                f"Longitude must be between -180 and 180 degrees. Received {lon}"
+            )
+
+        if use_spherical_bounding_box:
+            raise NotImplementedError(
+                "Computing the intersection using the spherical bounding box"
+                "is not yet supported."
             )
-        elif method == "accurate":
-            raise NotImplementedError("Accurate method not yet implemented.")
         else:
-            raise ValueError(f"Invalid method: {method}.")
-        return edges.squeeze()
+            edges = constant_lon_intersections_no_extreme(
+                lon, self.edge_node_x.values, self.edge_node_y.values, self.n_edge
+            )
+            return edges.squeeze()
 
-    def get_faces_at_constant_longitude(self, lon, method="fast"):
-        """Identifies the faces of the grid that intersect with a specified
-        constant longitude.
+    def get_faces_at_constant_longitude(self, lon, use_spherical_bounding_box=False):
+        """
+        Identifies the indices of faces that intersect with a line of constant longitude.
 
-        This method finds the faces (or cells) of the grid that intersect a given longitude
-        by first identifying the intersecting edges and then determining the faces connected
-        to these edges. The method used for identifying edges can be adjusted with the `method`
-        parameter.
+        When `use_spherical_bounding_box` is set to `True`,
+        the bounding box for each face is computed by representing each edge as a great circle arc.
+        This approach takes into account the extreme minimums or maximums along the arcs.
 
         Parameters
         ----------
         lon : float
             The longitude at which to identify intersecting faces, in degrees.
-        method : str, optional
-            The computational method used to determine intersecting edges. Options are:
-            - 'fast': Uses a faster but potentially less accurate method for determining intersections.
-            - 'accurate': Uses a slower but more precise method.
-            Default is 'fast'.
+        use_spherical_bounding_box : bool, optional
+            If `True`,
+            computes the bounding box for each face using great circle arcs for edges
+            and considers extreme minimums or maximums to increase accuracy.
+            Defaults to `False`.
 
         Returns
         -------
-        faces : array
-            An array of unique face indices that intersect the specified longitude.
-            Faces that are invalid or missing (e.g., with a fill value) are excluded
-            from the result.
+        faces : numpy.ndarray
+            An array of face indices that intersect with the specified longitude.
         """
-        edges = self.get_edges_at_constant_longitude(lon, method)
-        faces = np.unique(self.edge_face_connectivity[edges].data.ravel())
 
-        return faces[faces != INT_FILL_VALUE]
+        if use_spherical_bounding_box:
+            raise NotImplementedError(
+                "Computing the intersection using the spherical bounding box is not"
+                "yet supported."
+            )
+        else:
+            edges = self.get_edges_at_constant_longitude(
+                lon, use_spherical_bounding_box
+            )
+            faces = np.unique(self.edge_face_connectivity[edges].data.ravel())
+
+            return faces[faces != INT_FILL_VALUE]
diff --git a/uxarray/grid/intersections.py b/uxarray/grid/intersections.py
index 3d3e9d96e..147a4b156 100644
--- a/uxarray/grid/intersections.py
+++ b/uxarray/grid/intersections.py
@@ -11,16 +11,19 @@
 
 
 @njit(parallel=True, nogil=True, cache=True)
-def fast_constant_lat_intersections(lat, edge_node_z, n_edge):
-    """Determine which edges intersect a constant line of latitude on a sphere,
-    including edges that lie exactly along the latitude.
+def constant_lat_intersections_no_extreme(lat, edge_node_z, n_edge):
+    """Determine which edges intersect a constant line of latitude on a
+    sphere, without wrapping to the opposite longitude, with extremes
+    along each great circle arc not considered.
 
     Parameters
     ----------
     lat:
         Constant latitude value in degrees.
-    edge_node_z:
-        Array of shape (n_edge, 2) containing z-coordinates of the edge nodes.
+    edge_node_x:
+        Array of shape (n_edge, 2) containing x-coordinates of the edge nodes.
+    edge_node_y:
+        Array of shape (n_edge, 2) containing y-coordinates of the edge nodes.
     n_edge:
         Total number of edges to check.
 
@@ -38,15 +41,8 @@ def fast_constant_lat_intersections(lat, edge_node_z, n_edge):
 
     # Iterate through each edge and check for intersections
     for i in prange(n_edge):
-        # Get the z-coordinates of the edge's nodes
-        z0 = edge_node_z[i, 0]
-        z1 = edge_node_z[i, 1]
-
         # Check if the edge crosses the constant latitude or lies exactly on it
-        if (z0 - z_constant) * (z1 - z_constant) < 0.0 or (
-            abs(z0 - z_constant) < ERROR_TOLERANCE
-            and abs(z1 - z_constant) < ERROR_TOLERANCE
-        ):
+        if edge_intersects_constant_lat_no_extreme(edge_node_z[i], z_constant):
             intersecting_edges_mask[i] = 1
 
     intersecting_edges = np.argwhere(intersecting_edges_mask)
@@ -54,10 +50,28 @@ def fast_constant_lat_intersections(lat, edge_node_z, n_edge):
     return np.unique(intersecting_edges)
 
 
+@njit(cache=True, nogil=True)
+def edge_intersects_constant_lat_no_extreme(edge_node_z, z_constant):
+    """Helper to compute whether an edge intersects a line of constant latitude."""
+
+    # z coordinate of edge nodes
+    z0 = edge_node_z[0]
+    z1 = edge_node_z[1]
+
+    if (z0 - z_constant) * (z1 - z_constant) < 0.0 or (
+        abs(z0 - z_constant) < ERROR_TOLERANCE
+        and abs(z1 - z_constant) < ERROR_TOLERANCE
+    ):
+        return True
+    else:
+        return False
+
+
 @njit(parallel=True, nogil=True, cache=True)
-def fast_constant_lon_intersections(lon, edge_node_x, edge_node_y, n_edge):
+def constant_lon_intersections_no_extreme(lon, edge_node_x, edge_node_y, n_edge):
     """Determine which edges intersect a constant line of longitude on a
-    sphere, without wrapping to the opposite longitude.
+    sphere, without wrapping to the opposite longitude, with extremes
+    along each great circle arc not considered.
 
     Parameters
     ----------
@@ -108,6 +122,70 @@ def fast_constant_lon_intersections(lon, edge_node_x, edge_node_y, n_edge):
     return np.unique(intersecting_edges)
 
 
+@njit
+def constant_lat_intersections_face_bounds(lat, face_min_lat_rad, face_max_lat_rad):
+    """Identifies the candidate faces on a grid that intersect with a given
+    constant latitude.
+
+    This function checks whether the specified latitude, `lat`, in degrees lies within
+    the latitude bounds of grid faces, defined by `face_min_lat_rad` and `face_max_lat_rad`,
+    which are given in radians. The function returns the indices of the faces where the
+    latitude is within these bounds.
+
+    Parameters
+    ----------
+    lat : float
+        The latitude in degrees for which to find intersecting faces.
+    face_min_lat_rad : numpy.ndarray
+        A 1D array containing the minimum latitude bounds (in radians) of each face.
+    face_max_lat_rad : numpy.ndarray
+        A 1D array containing the maximum latitude bounds (in radians) of each face.
+
+    Returns
+    -------
+    candidate_faces : numpy.ndarray
+        A 1D array containing the indices of the faces that intersect with the given latitude.
+    """
+    lat = np.deg2rad(lat)
+    within_bounds = (face_min_lat_rad <= lat) & (face_max_lat_rad >= lat)
+    candidate_faces = np.where(within_bounds)[0]
+    return candidate_faces
+
+
+@njit(cache=True)
+def constant_lon_intersections_face_bounds(lon, face_min_lon_rad, face_max_lon_rad):
+    """Identifies the candidate faces on a grid that intersect with a given
+    constant longitude.
+
+    This function checks whether the specified longitude, `lon`, in degrees lies within
+    the longitude bounds of grid faces, defined by `face_min_lon_rad` and `face_max_lon_rad`,
+    which are given in radians. The function returns the indices of the faces where the
+    longitude is within these bounds.
+
+    Parameters
+    ----------
+    lon : float
+        The longitude in degrees for which to find intersecting faces.
+    face_min_lon_rad : numpy.ndarray
+        A 1D array containing the minimum longitude bounds (in radians) of each face.
+    face_max_lon_rad : numpy.ndarray
+        A 1D array containing the maximum longitude bounds (in radians) of each face.
+
+    Returns
+    -------
+    candidate_faces : numpy.ndarray
+        A 1D array containing the indices of the faces that intersect with the given longitude.
+    """
+
+    # lon = np.deg2rad(lon)
+    # lon = (lon + 2 * np.pi) % (2 * np.pi)
+    # within_bounds = (face_min_lon_rad <= lon) & (face_max_lon_rad >= lon)
+    # candidate_faces = np.where(within_bounds)[0]
+    # return candidate_faces
+
+    raise NotImplementedError
+
+
 def gca_gca_intersection(gca1_cart, gca2_cart, fma_disabled=True):
     """Calculate the intersection point(s) of two Great Circle Arcs (GCAs) in a
     Cartesian coordinate system.