Repository: UXARRAY/uxarray

Issue:
Performance Issue: Unnecessary Re-conversion of Coordinates in `point_within_gca`
### Proposed new feature or change:

The current implementation of point_within_gca unnecessarily re-converts lat/lon coordinates from Cartesian coordinates, which negatively impacts performance.

The following lines re-compute lat/lon from Cartesian coordinates:

python
```
GCRv0_lonlat = np.array(
        _xyz_to_lonlat_rad_no_norm(gca_cart[0][0], gca_cart[0][1], gca_cart[0][2])
    )
GCRv1_lonlat = np.array(
    _xyz_to_lonlat_rad_no_norm(gca_cart[1][0], gca_cart[1][1], gca_cart[1][2])
)
```

However, Uxarray already stores lat/lon coordinates, which can be passed directly to avoid redundant computation. This approach reduces both time and space complexity while aligning better with the overall architecture.

Proposed Solution:

- Pass the existing lat/lon coordinates directly to eliminate the unnecessary conversion.
- Update upstream function calls to `point_within_gca`

Discussion:
philipc2: @hongyuchen1030 

I still need to finish adjusting all the tests and ensure that everything is passing, but I wanted to quickly updated you on some initial timings. 
- 30km (~40,000 mostly hexagonal faces): 4 seconds

This was without any parallelization. Since the bounding box for each face can be computed interpedently, my next step will be to parallelize the loop that calls `_populate_face_latlon_bound`


github-actions[bot]: 
## ASV Benchmarking

<details>
  <summary>Benchmark Comparison Results</summary>

  
Benchmarks that have improved:

| Change   | Before [01efecde]    | After [9fe8df08]    |   Ratio | Benchmark (Parameter)                           |
|----------|----------------------|---------------------|---------|-------------------------------------------------|
| -        | 465M                 | 372M                |     0.8 | mpas_ocean.Integrate.peakmem_integrate('480km') |

Benchmarks that have stayed the same:

| Change   | Before [01efecde]    | After [9fe8df08]    | Ratio   | Benchmark (Parameter)                                                                |
|----------|----------------------|---------------------|---------|--------------------------------------------------------------------------------------|
|          | 2.81±0.01s           | 2.81±0.03s          | 1.00    | import.Imports.timeraw_import_uxarray                                                |
|          | 660±9μs              | 641±9μs             | 0.97    | mpas_ocean.CheckNorm.time_check_norm('120km')                                        |
|          | 428±6μs              | 425±3μs             | 0.99    | mpas_ocean.CheckNorm.time_check_norm('480km')                                        |
|          | 615±10ms             | 625±7ms             | 1.02    | mpas_ocean.ConnectivityConstruction.time_face_face_connectivity('120km')             |
|          | 39.5±0.8ms           | 40.5±0.3ms          | 1.03    | mpas_ocean.ConnectivityConstruction.time_face_face_connectivity('480km')             |
|          | 1.71±0.01ms          | 1.70±0.03ms         | 0.99    | mpas_ocean.ConnectivityConstruction.time_n_nodes_per_face('120km')                   |
|          | 488±5μs              | 517±40μs            | 1.06    | mpas_ocean.ConnectivityConstruction.time_n_nodes_per_face('480km')                   |
|          | 5.79±0.04ms          | 5.73±0.06ms         | 0.99    | mpas_ocean.ConstructFaceLatLon.time_cartesian_averaging('120km')                     |
|          | 3.40±0.03ms          | 3.40±0.09ms         | 1.00    | mpas_ocean.ConstructFaceLatLon.time_cartesian_averaging('480km')                     |
|          | 3.53±0s              | 3.53±0.01s          | 1.00    | mpas_ocean.ConstructFaceLatLon.time_welzl('120km')                                   |
|          | 227±1ms              | 225±0.6ms           | 0.99    | mpas_ocean.ConstructFaceLatLon.time_welzl('480km')                                   |
|          | 1.14±0μs             | 1.09±0μs            | 0.95    | mpas_ocean.ConstructTreeStructures.time_ball_tree('120km')                           |
|          | 262±0.9ns            | 278±3ns             | 1.06    | mpas_ocean.ConstructTreeStructures.time_ball_tree('480km')                           |
|          | 718±2ns              | 749±2ns             | 1.04    | mpas_ocean.ConstructTreeStructures.time_kd_tree('120km')                             |
|          | 261±3ns              | 271±8ns             | 1.04    | mpas_ocean.ConstructTreeStructures.time_kd_tree('480km')                             |
|          | failed               | failed              | n/a     | mpas_ocean.CrossSections.time_constant_lat_fast('120km', 1)                          |
|          | failed               | failed              | n/a     | mpas_ocean.CrossSections.time_constant_lat_fast('120km', 2)                          |
|          | failed               | failed              | n/a     | mpas_ocean.CrossSections.time_constant_lat_fast('120km', 4)                          |
|          | failed               | failed              | n/a     | mpas_ocean.CrossSections.time_constant_lat_fast('120km', 8)                          |
|          | failed               | failed              | n/a     | mpas_ocean.CrossSections.time_constant_lat_fast('480km', 1)                          |
|          | failed               | failed              | n/a     | mpas_ocean.CrossSections.time_constant_lat_fast('480km', 2)                          |
|          | failed               | failed              | n/a     | mpas_ocean.CrossSections.time_constant_lat_fast('480km', 4)                          |
|          | failed               | failed              | n/a     | mpas_ocean.CrossSections.time_constant_lat_fast('480km', 8)                          |
|          | 122±0.3ms            | 123±0.6ms           | 1.00    | mpas_ocean.DualMesh.time_dual_mesh_construction('120km')                             |
|          | 8.48±0.2ms           | 8.25±0.04ms         | 0.97    | mpas_ocean.DualMesh.time_dual_mesh_construction('480km')                             |
|          | 1.01±0s              | 1.02±0s             | 1.01    | mpas_ocean.GeoDataFrame.time_to_geodataframe('120km', False)                         |
|          | 52.3±0.6ms           | 52.4±0.6ms          | 1.00    | mpas_ocean.GeoDataFrame.time_to_geodataframe('120km', True)                          |
|          | 76.3±0.9ms           | 76.1±0.5ms          | 1.00    | mpas_ocean.GeoDataFrame.time_to_geodataframe('480km', False)                         |
|          | 5.01±0.2ms           | 4.92±0.1ms          | 0.98    | mpas_ocean.GeoDataFrame.time_to_geodataframe('480km', True)                          |
|          | 318M                 | 318M                | 1.00    | mpas_ocean.Gradient.peakmem_gradient('120km')                                        |
|          | 295M                 | 296M                | 1.00    | mpas_ocean.Gradient.peakmem_gradient('480km')                                        |
|          | 2.71±0.02ms          | 2.74±0.03ms         | 1.01    | mpas_ocean.Gradient.time_gradient('120km')                                           |
|          | 311±1μs              | 318±6μs             | 1.02    | mpas_ocean.Gradient.time_gradient('480km')                                           |
|          | failed               | failed              | n/a     | mpas_ocean.HoleEdgeIndices.time_construct_hole_edge_indices('120km')                 |
|          | failed               | failed              | n/a     | mpas_ocean.HoleEdgeIndices.time_construct_hole_edge_indices('480km')                 |
|          | 389M                 | 389M                | 1.00    | mpas_ocean.Integrate.peakmem_integrate('120km')                                      |
|          | 174±0.9ms            | 172±0.5ms           | 0.99    | mpas_ocean.Integrate.time_integrate('120km')                                         |
|          | 12.0±0.2ms           | 11.9±0.08ms         | 0.99    | mpas_ocean.Integrate.time_integrate('480km')                                         |
|          | 333±3ms              | 332±2ms             | 1.00    | mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('120km', 'exclude') |
|          | 336±3ms              | 330±2ms             | 0.98    | mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('120km', 'include') |
|          | 331±0.9ms            | 335±2ms             | 1.01    | mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('120km', 'split')   |
|          | 21.2±0.3ms           | 21.4±0.1ms          | 1.01    | mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('480km', 'exclude') |
|          | 21.1±0.3ms           | 21.5±0.2ms          | 1.02    | mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('480km', 'include') |
|          | 21.4±0.1ms           | 21.3±0.1ms          | 1.00    | mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('480km', 'split')   |
|          | 55.2±0.09ms          | 55.1±0.06ms         | 1.00    | mpas_ocean.RemapDownsample.time_inverse_distance_weighted_remapping                  |
|          | 45.0±0.08ms          | 44.9±0.1ms          | 1.00    | mpas_ocean.RemapDownsample.time_nearest_neighbor_remapping                           |
|          | 359±0.8ms            | 359±0.7ms           | 1.00    | mpas_ocean.RemapUpsample.time_inverse_distance_weighted_remapping                    |
|          | 264±0.2ms            | 264±0.7ms           | 1.00    | mpas_ocean.RemapUpsample.time_nearest_neighbor_remapping                             |
|          | 291M                 | 291M                | 1.00    | quad_hexagon.QuadHexagon.peakmem_open_dataset                                        |
|          | 291M                 | 291M                | 1.00    | quad_hexagon.QuadHexagon.peakmem_open_grid                                           |
|          | 6.15±0.06ms          | 6.20±0.08ms         | 1.01    | quad_hexagon.QuadHexagon.time_open_dataset                                           |
|          | 5.29±0.01ms          | 5.29±0.03ms         | 1.00    | quad_hexagon.QuadHexagon.time_open_grid                                              |

Benchmarks that have got worse:

| Change   | Before [01efecde]    | After [9fe8df08]    | Ratio   | Benchmark (Parameter)                                                                                                                 |
|----------|----------------------|---------------------|---------|---------------------------------------------------------------------------------------------------------------------------------------|
| !        | 375M                 | failed              | n/a     | face_bounds.FaceBounds.peakmem_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/mpas/QU/oQU480.231010.nc'))    |
| !        | 375M                 | failed              | n/a     | face_bounds.FaceBounds.peakmem_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/scrip/outCSne8/outCSne8.nc'))  |
| !        | 400M                 | failed              | n/a     | face_bounds.FaceBounds.peakmem_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/geoflow-small/grid.nc')) |
| !        | 445M                 | failed              | n/a     | face_bounds.FaceBounds.peakmem_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/quad-hexagon/grid.nc'))  |
| !        | 1.57±0.02s           | failed              | n/a     | face_bounds.FaceBounds.time_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/mpas/QU/oQU480.231010.nc'))       |
| !        | 220±0.3ms            | failed              | n/a     | face_bounds.FaceBounds.time_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/scrip/outCSne8/outCSne8.nc'))     |
| !        | 1.99±0.01s           | failed              | n/a     | face_bounds.FaceBounds.time_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/geoflow-small/grid.nc'))    |
| !        | 7.61±0.1ms           | failed              | n/a     | face_bounds.FaceBounds.time_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/quad-hexagon/grid.nc'))     |

</details>

philipc2: @hongyuchen1030 

I might scrap the previous methods for computing the intersections and exclusively support the methods that use the bounding box now that the performance has significantly improved. 


hongyuchen1030: > @hongyuchen1030
> 
> I might scrap the previous methods for computing the intersections and exclusively support the methods that use the bounding box now that the performance has significantly improved.

sounds great!
philipc2: > Did you remove all face_latlon_bounds usages in this PR?
> 
> Please prioritize merging the face_latlon_bounds here. You can set the `use_spherical_bounds` by default as `false`, but please try to at least keep the implementation in this PR

I did the opposite. I removed the usages of the naive implementation. All constant longitude and constant latitude queries now use the `face_latlon_bounds`.

There is no longer a `use_spherical_bounds` parameter. I hope this clarifies. 
hongyuchen1030: > > Did you remove all face_latlon_bounds usages in this PR?
> > Please prioritize merging the face_latlon_bounds here. You can set the `use_spherical_bounds` by default as `false`, but please try to at least keep the implementation in this PR
> 
> I did the opposite. I removed the usages of the naive implementation. All constant longitude and constant latitude queries now use the `face_latlon_bounds`.
> 
> There is no longer a `use_spherical_bounds` parameter. I hope this clarifies.

Ohh I see. I guess I was misled by those “No implementation warning ” for `use_spherical_bounds` . Is there any reason we kept those?
philipc2: > > > Did you remove all face_latlon_bounds usages in this PR?
> > > Please prioritize merging the face_latlon_bounds here. You can set the `use_spherical_bounds` by default as `false`, but please try to at least keep the implementation in this PR
> > 
> > 
> > I did the opposite. I removed the usages of the naive implementation. All constant longitude and constant latitude queries now use the `face_latlon_bounds`.
> > There is no longer a `use_spherical_bounds` parameter. I hope this clarifies.
> 
> Ohh I see. I guess I was misled by those “No implementation warning ” for `use_spherical_bounds` . Is there any reason we kept those?

Those are the `edge` queries. Right now I am only using the `face_bounds` to determine if a face intersects with a line of longitude or latitude. 

Patch:
diff --git a/docs/api.rst b/docs/api.rst
index 1b03b5d3..f9146e69 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -126,6 +126,8 @@ Descriptors
    Grid.descriptors
    Grid.face_areas
    Grid.bounds
+   Grid.face_bounds_lon
+   Grid.face_bounds_lat
    Grid.edge_node_distances
    Grid.edge_face_distances
    Grid.antimeridian_face_indices
diff --git a/uxarray/__init__.py b/uxarray/__init__.py
index 857c65f6..c9cdff9e 100644
--- a/uxarray/__init__.py
+++ b/uxarray/__init__.py
@@ -1,4 +1,9 @@
 import uxarray.constants
+import sys
+#
+# # TODO: numba recursion limit ?
+
+
 from .core.api import open_grid, open_dataset, open_mfdataset
 
 from .core.dataset import UxDataset
@@ -34,6 +39,7 @@ def disable_fma():
 
 
 disable_fma()
+sys.setrecursionlimit(10000)
 
 
 __all__ = (
diff --git a/uxarray/cross_sections/dataarray_accessor.py b/uxarray/cross_sections/dataarray_accessor.py
index de35739a..d8c357e5 100644
--- a/uxarray/cross_sections/dataarray_accessor.py
+++ b/uxarray/cross_sections/dataarray_accessor.py
@@ -17,63 +17,89 @@ class UxDataArrayCrossSectionAccessor:
         prefix = "<uxarray.UxDataArray.cross_section>\n"
         methods_heading = "Supported Methods:\n"
 
-        methods_heading += "  * constant_latitude(center_coord, k, element, **kwargs)\n"
+        methods_heading += "  * constant_latitude(lat)\n"
+        methods_heading += "  * constant_longitude(lon)\n"
 
         return prefix + methods_heading
 
-    def constant_latitude(self, lat: float, use_spherical_bounding_box=False):
-        """Extracts a cross-section of the data array at a specified constant
-        latitude.
+    def constant_latitude(self, lat: float):
+        """Extracts a cross-section of the data array by selecting all faces that
+        intersect with a specified line of constant latitude.
 
         Parameters
         ----------
         lat : float
             The latitude 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.
+            Must be between -90.0 and 90.0
+
+        Returns
+        -------
+        uxarray.UxDataArray
+            A subset of the original data array containing only the faces that intersect
+            with the specified latitude.
 
         Raises
         ------
         ValueError
-            If no intersections are found at the specified latitude, a ValueError is raised.
+            If no intersections are found at the specified longitude or the data variable is not face-centered.
 
         Examples
         --------
-        >>> uxda.constant_latitude_cross_section(lat=-15.5)
+        >>> # Extract data at 15.5°S latitude
+        >>> cross_section = uxda.constant_latitude(lat=-15.5)
+
+        Notes
+        -----
+        The initial execution time may be significantly longer than subsequent runs
+        due to Numba's just-in-time compilation. Subsequent calls will be faster due to caching.
         """
-        faces = self.uxda.uxgrid.get_faces_at_constant_latitude(
-            lat, use_spherical_bounding_box
-        )
+        if not self.uxda._face_centered():
+            raise ValueError(
+                "Cross sections are only supported for face-centered data variables."
+            )
+
+        faces = self.uxda.uxgrid.get_faces_at_constant_latitude(lat)
 
         return self.uxda.isel(n_face=faces)
 
-    def constant_longitude(self, lon: float, use_spherical_bounding_box=False):
-        """Extracts a cross-section of the data array at a specified constant
-        longitude.
+    def constant_longitude(self, lon: float):
+        """Extracts a cross-section of the data array by selecting all faces that
+        intersect with a specified line of constant longitude.
 
         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.
+            The latitude at which to extract the cross-section, in degrees.
+            Must be between -180.0 and 180.0
+
+        Returns
+        -------
+        uxarray.UxDataArray
+            A subset of the original data array containing only the faces that intersect
+            with the specified longitude.
 
         Raises
         ------
         ValueError
-            If no intersections are found at the specified longitude, a ValueError is raised.
+            If no intersections are found at the specified longitude or the data variable is not face-centered.
 
         Examples
         --------
-        >>> uxda.constant_longitude_cross_section(lon=-15.5)
+        >>> # Extract data at 0° longitude
+        >>> cross_section = uxda.constant_latitude(lon=0.0)
 
         Notes
         -----
-        The accuracy and performance of the function can be controlled using the `method` parameter.
-        For higher precision requreiments, consider using method='acurate'.
+        The initial execution time may be significantly longer than subsequent runs
+        due to Numba's just-in-time compilation. Subsequent calls will be faster due to caching.
         """
+        if not self.uxda._face_centered():
+            raise ValueError(
+                "Cross sections are only supported for face-centered data variables."
+            )
+
         faces = self.uxda.uxgrid.get_faces_at_constant_longitude(
-            lon, use_spherical_bounding_box
+            lon,
         )
 
         return self.uxda.isel(n_face=faces)
diff --git a/uxarray/cross_sections/grid_accessor.py b/uxarray/cross_sections/grid_accessor.py
index 5b905516..146da10a 100644
--- a/uxarray/cross_sections/grid_accessor.py
+++ b/uxarray/cross_sections/grid_accessor.py
@@ -17,54 +17,50 @@ class GridCrossSectionAccessor:
         prefix = "<uxarray.Grid.cross_section>\n"
         methods_heading = "Supported Methods:\n"
 
-        methods_heading += "  * constant_latitude(lat, )\n"
+        methods_heading += "  * constant_latitude(lat, return_face_indices)\n"
+        methods_heading += "  * constant_longitude(lon, return_face_indices)\n"
         return prefix + methods_heading
 
     def constant_latitude(
-        self, lat: float, return_face_indices=False, use_spherical_bounding_box=False
+        self,
+        lat: float,
+        return_face_indices: bool = False,
     ):
-        """Extracts a cross-section of the grid at a specified constant
-        latitude.
+        """Extracts a cross-section of the grid by selecting all faces that
+        intersect with a specified line of constant latitude.
 
         Parameters
         ----------
         lat : float
             The latitude at which to extract the cross-section, in degrees.
+            Must be between -90.0 and 90.0
         return_face_indices : bool, optional
-            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.
-        use_spherical_bounding_box : bool, optional
-            If True, uses a spherical bounding box for intersection calculations.
-
+            If True, also returns the indices of the faces that intersect with the line of constant latitude.
 
         Returns
         -------
-        grid_at_constant_lat : Grid
-            The grid with faces that interesected at a given lattitude
-        faces : array, optional
-            The indices of the faces that intersect with the specified latitude. This is only
-            returned if `return_face_indices` is set to True.
+        uxarray.Grid
+            A subset of the original grid containing only the faces that intersect
+            with the specified latitude.
 
         Raises
         ------
         ValueError
-            If no intersections are found at the specified latitude, a ValueError is raised.
+            If no intersections are found at the specified latitude.
 
         Examples
         --------
-        >>> grid, indices = grid.cross_section.constant_latitude(
-        ...     lat=30.0, return_face_indices=True
-        ... )
-        >>> grid = grid.cross_section.constant_latitude(lat=-15.5)
+        >>> # Extract data at 15.5°S latitude
+        >>> cross_section = grid.constant_latitude(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'.
+        The initial execution time may be significantly longer than subsequent runs
+        due to Numba's just-in-time compilation. Subsequent calls will be faster due to caching.
         """
+
         faces = self.uxgrid.get_faces_at_constant_latitude(
-            lat, use_spherical_bounding_box
+            lat,
         )
 
         if len(faces) == 0:
@@ -78,49 +74,44 @@ class GridCrossSectionAccessor:
             return grid_at_constant_lat
 
     def constant_longitude(
-        self, lon: float, use_spherical_bounding_box=False, return_face_indices=False
+        self,
+        lon: float,
+        return_face_indices: bool = False,
     ):
-        """Extracts a cross-section of the grid at a specified constant
-        longitude.
+        """Extracts a cross-section of the grid by selecting all faces that
+        intersect with a specified line of constant longitude.
 
         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.
+            Must be between -90.0 and 90.0
         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.
+            If True, also returns the indices of the faces that intersect with the line of constant longitude.
 
         Returns
         -------
-        grid_at_constant_lon : Grid
-            The grid with faces that interesected at a given longitudes
-        faces : array, optional
-            The indices of the faces that intersect with the specified longitude. This is only
-            returned if `return_face_indices` is set to True.
+        uxarray.Grid
+            A subset of the original grid containing only the faces that intersect
+            with the specified longitude.
 
         Raises
         ------
         ValueError
-            If no intersections are found at the specified longitude, a ValueError is raised.
+            If no intersections are found at the specified longitude.
 
         Examples
         --------
-        >>> grid, indices = grid.cross_section.constant_longitude(
-        ...     lat=0.0, return_face_indices=True
-        ... )
-        >>> grid = grid.cross_section.constant_longitude(lat=20.0)
+        >>> # Extract data at 0° longitude
+        >>> cross_section = grid.constant_latitude(lon=0.0)
 
         Notes
         -----
-        The accuracy and performance of the function can be controlled using the `method` parameter.
-        For higher precision requreiments, consider using method='acurate'.
+        The initial execution time may be significantly longer than subsequent runs
+        due to Numba's just-in-time compilation. Subsequent calls will be faster due to caching.
         """
         faces = self.uxgrid.get_faces_at_constant_longitude(
-            lon, use_spherical_bounding_box
+            lon,
         )
 
         if len(faces) == 0:
diff --git a/uxarray/grid/arcs.py b/uxarray/grid/arcs.py
index 6be15800..3ccb4501 100644
--- a/uxarray/grid/arcs.py
+++ b/uxarray/grid/arcs.py
@@ -1,9 +1,10 @@
 import numpy as np
+import math
 
-# from uxarray.grid.coordinates import node_xyz_to_lonlat_rad, normalize_in_place
+
+from uxarray.grid.coordinates import _xyz_to_lonlat_rad_scalar
 
 from uxarray.grid.coordinates import (
-    _xyz_to_lonlat_rad_scalar,
     _normalize_xyz_scalar,
 )
 
@@ -11,7 +12,7 @@ from uxarray.grid.utils import _angle_of_2_vectors
 
 from uxarray.constants import ERROR_TOLERANCE, MACHINE_EPSILON
 
-from uxarray.utils.computing import isclose, cross, dot, allclose
+from uxarray.utils.computing import isclose, dot
 
 from numba import njit
 
@@ -27,111 +28,206 @@ def _to_list(obj):
     return obj
 
 
+def _point_within_gca_cartesian(pt_xyz, gca_xyz, is_directed=False):
+    pt_xyz = np.asarray(pt_xyz)
+    gca_xyz = np.asarray(gca_xyz)
+
+    pt_lonlat = np.array(
+        _xyz_to_lonlat_rad_scalar(pt_xyz[0], pt_xyz[1], pt_xyz[2], normalize=False)
+    )
+    gca_a_xyz = gca_xyz[0]
+
+    gca_a_lonlat = np.array(
+        _xyz_to_lonlat_rad_scalar(
+            gca_xyz[0][0], gca_xyz[0][1], gca_xyz[0][2], normalize=False
+        )
+    )
+    gca_b_xyz = gca_xyz[1]
+
+    gca_b_lonlat = np.array(
+        _xyz_to_lonlat_rad_scalar(
+            gca_xyz[1][0], gca_xyz[1][1], gca_xyz[1][2], normalize=False
+        )
+    )
+
+    return point_within_gca(
+        pt_xyz,
+        pt_lonlat,
+        gca_a_xyz,
+        gca_a_lonlat,
+        gca_b_xyz,
+        gca_b_lonlat,
+        is_directed=is_directed,
+    )
+
+
 @njit(cache=True)
-def _point_within_gca_body(
-    angle, gca_cart, pt, GCRv0_lonlat, GCRv1_lonlat, pt_lonlat, is_directed
+def point_within_gca(
+    pt_xyz,
+    pt_lonlat,
+    gca_a_xyz,
+    gca_a_lonlat,
+    gca_b_xyz,
+    gca_b_lonlat,
+    is_directed=False,
 ):
-    angle = _angle_of_2_vectors(gca_cart[0], gca_cart[1])
-    if allclose(angle, np.pi, rtol=0.0, atol=MACHINE_EPSILON):
+    """Check if a point lies on a given Great Circle Arc (GCA). The anti-
+    meridian case is also considered.
+
+    Parameters
+    ----------
+    pt : numpy.ndarray (float)
+        Cartesian coordinates of the point.
+    gca_cart : numpy.ndarray of shape (2, 3), (np.float or gmpy2.mpfr)
+        Cartesian coordinates of the Great Circle Arc (GCR).
+    is_directed : bool, optional, default = False
+        If True, the GCA is considered to be directed, which means it can only from v0-->v1. If False, the GCA is undirected,
+        and we will always assume the small circle (The one less than 180 degree) side is the GCA. The default is False.
+        For the case of the anti-podal case, the direction is v_0--> the pole point that on the same hemisphere as v_0-->v_1
+
+    Returns
+    -------
+    bool
+        True if the point lies between the two endpoints of the GCR, False otherwise.
+
+    Raises
+    ------
+    ValueError
+        If the input GCR spans exactly 180 degrees (π radians), as this GCR can have multiple planes.
+        In such cases, consider breaking the GCR into two separate GCRs.
+
+    ValueError
+        If the input GCR spans more than 180 degrees (π radians).
+        In such cases, consider breaking the GCR into two separate GCRs.
+
+    Notes
+    -----
+    The function checks if the given point is on the Great Circle Arc by considering its cartesian coordinates and
+    accounting for the anti-meridian case.
+
+    The anti-meridian case occurs when the GCR crosses the anti-meridian (0 longitude).
+    In this case, the function handles scenarios where the GCA spans across more than 180 degrees, requiring specific operation.
+
+    The function relies on the `_angle_of_2_vectors` and `is_between` functions to perform the necessary calculations.
+
+    Please ensure that the input coordinates are in radians and adhere to the ERROR_TOLERANCE value for floating-point comparisons.
+    """
+
+    # ==================================================================================================================
+    # 1. If the input GCR is exactly 180 degree, we throw an exception, since this GCR can have multiple planes
+    angle = _angle_of_2_vectors(gca_a_xyz, gca_b_xyz)
+
+    if np.allclose(angle, np.pi, rtol=0.0, atol=MACHINE_EPSILON):
         raise ValueError(
             "The input Great Circle Arc is exactly 180 degree, this Great Circle Arc can have multiple planes. "
             "Consider breaking the Great Circle Arc"
             "into two Great Circle Arcs"
         )
-
+    # ==================================================================================================================
     # See if the point is on the plane of the GCA, because we are dealing with floating point numbers with np.dot now
     # just using the rtol=MACHINE_EPSILON, atol=MACHINE_EPSILON, but consider using the more proper error tolerance
     # in the future
-    cross_product = cross(np.asarray(gca_cart[0]), np.asarray(gca_cart[1]))
 
-    if not allclose(
-        dot(np.asarray(cross_product), np.asarray(pt)),  # Custom dot function
+    cross_product = np.cross(gca_a_xyz, gca_b_xyz)
+
+    if not np.allclose(
+        np.dot(cross_product, pt_xyz),
         0,
         rtol=MACHINE_EPSILON,
         atol=MACHINE_EPSILON,
     ):
         return False
-
-    if isclose(
-        GCRv0_lonlat[0], GCRv1_lonlat[0], rtol=MACHINE_EPSILON, atol=MACHINE_EPSILON
+    # ==================================================================================================================
+    if np.isclose(
+        gca_a_lonlat[0], gca_b_lonlat[0], rtol=MACHINE_EPSILON, atol=MACHINE_EPSILON
     ):
         # If the pt and the GCA are on the same longitude (the y coordinates are the same)
-        if isclose(
-            GCRv0_lonlat[0], pt_lonlat[0], rtol=MACHINE_EPSILON, atol=MACHINE_EPSILON
+        if np.isclose(
+            gca_a_lonlat[0], pt_lonlat[0], rtol=MACHINE_EPSILON, atol=MACHINE_EPSILON
         ):
             # Now use the latitude to determine if the pt falls between the interval
-            return in_between(GCRv0_lonlat[1], pt_lonlat[1], GCRv1_lonlat[1])
+            return in_between(gca_a_lonlat[1], pt_lonlat[1], gca_b_lonlat[1])
         else:
-            # If the pt and the GCA are not on the same longitude when the GCA is a longnitude arc, then the pt is not on the GCA
+            # If the pt and the GCA are not on the same longitude when the GCA is a longitude arc, then the pt is not on the GCA
             return False
-
-    # If the longnitude span is exactly 180 degree, then the GCA goes through the pole point
+    # ==================================================================================================================
+    # If the longitude span is exactly 180 degree, then the GCA goes through the pole point
     # Or if one of the endpoints is on the pole point, then the GCA goes through the pole point
     if (
-        isclose(
-            abs(GCRv1_lonlat[0] - GCRv0_lonlat[0]),
+        np.isclose(
+            np.abs(gca_b_lonlat[0] - gca_a_lonlat[0]),
             np.pi,
             rtol=0.0,
             atol=MACHINE_EPSILON,
         )
-        or isclose(
-            abs(GCRv0_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=ERROR_TOLERANCE
+        or np.isclose(
+            np.abs(gca_a_lonlat[1]),
+            np.pi / 2,
+            rtol=ERROR_TOLERANCE,
+            atol=ERROR_TOLERANCE,
         )
-        or isclose(
-            abs(GCRv1_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=ERROR_TOLERANCE
+        or np.isclose(
+            np.abs(gca_b_lonlat[1]),
+            np.pi / 2,
+            rtol=ERROR_TOLERANCE,
+            atol=ERROR_TOLERANCE,
         )
     ):
+        # ==============================================================================================================
         # Special case, if the pt is on the pole point, then set its longitude to the GCRv0_lonlat[0]
         # Since the point is our calculated properly, we use the atol=ERROR_TOLERANCE and rtol=ERROR_TOLERANCE
-        if isclose(
-            abs(pt_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=ERROR_TOLERANCE
+        if np.isclose(
+            np.abs(pt_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=ERROR_TOLERANCE
         ):
-            pt_lonlat[0] = GCRv0_lonlat[0]
+            pt_lonlat[0] = gca_a_lonlat[0]
 
+        # ==============================================================================================================
         # Special case, if one of the GCA endpoints is on the pole point, and another endpoint is not
         # then we need to check if the pt is on the GCA
-        if isclose(
-            abs(GCRv0_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=0.0
-        ) or isclose(abs(GCRv1_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=0.0):
+        if np.isclose(
+            abs(gca_a_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=0.0
+        ) or np.isclose(
+            abs(gca_b_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=0.0
+        ):
             # Identify the non-pole endpoint
             non_pole_endpoint = None
-            if not isclose(
-                abs(GCRv0_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=0.0
+            if not np.isclose(
+                abs(gca_a_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=0.0
             ):
-                non_pole_endpoint = GCRv0_lonlat
-            elif not isclose(
-                abs(GCRv1_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=0.0
+                non_pole_endpoint = gca_a_lonlat
+            elif not np.isclose(
+                abs(gca_b_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=0.0
             ):
-                non_pole_endpoint = GCRv1_lonlat
+                non_pole_endpoint = gca_b_lonlat
 
-            if non_pole_endpoint is not None and not isclose(
+            if non_pole_endpoint is not None and not np.isclose(
                 non_pole_endpoint[0], pt_lonlat[0], rtol=ERROR_TOLERANCE, atol=0.0
             ):
                 return False
-
-        if not isclose(
-            GCRv0_lonlat[0], pt_lonlat[0], rtol=ERROR_TOLERANCE, atol=0.0
-        ) and not isclose(
-            GCRv1_lonlat[0], pt_lonlat[0], rtol=ERROR_TOLERANCE, atol=0.0
+        # ==============================================================================================================
+        if not np.isclose(
+            gca_a_lonlat[0], pt_lonlat[0], rtol=ERROR_TOLERANCE, atol=0.0
+        ) and not np.isclose(
+            gca_b_lonlat[0], pt_lonlat[0], rtol=ERROR_TOLERANCE, atol=0.0
         ):
             return False
         else:
             # Determine the pole latitude and latitude extension
-            if (GCRv0_lonlat[1] > 0.0 and GCRv1_lonlat[1] > 0.0) or (
-                GCRv0_lonlat[1] < 0.0 and GCRv1_lonlat[1] < 0.0
+            if (gca_a_lonlat[1] > 0.0 and gca_b_lonlat[1] > 0.0) or (
+                gca_a_lonlat[1] < 0.0 and gca_b_lonlat[1] < 0.0
             ):
-                pole_lat = np.pi / 2 if GCRv0_lonlat[1] > 0.0 else -np.pi / 2
+                pole_lat = np.pi / 2 if gca_a_lonlat[1] > 0.0 else -np.pi / 2
                 lat_extend = (
-                    abs(np.pi / 2 - abs(GCRv0_lonlat[1]))
+                    abs(np.pi / 2 - abs(gca_a_lonlat[1]))
                     + np.pi / 2
-                    + abs(GCRv1_lonlat[1])
+                    + abs(gca_b_lonlat[1])
                 )
             else:
-                pole_lat = _decide_pole_latitude(GCRv0_lonlat[1], GCRv1_lonlat[1])
+                pole_lat = _decide_pole_latitude(gca_a_lonlat[1], gca_b_lonlat[1])
                 lat_extend = (
-                    abs(np.pi / 2 - abs(GCRv0_lonlat[1]))
+                    abs(np.pi / 2 - abs(gca_a_lonlat[1]))
                     + np.pi / 2
-                    + abs(GCRv1_lonlat[1])
+                    + abs(gca_b_lonlat[1])
                 )
 
             if is_directed and lat_extend >= np.pi:
@@ -139,109 +235,44 @@ def _point_within_gca_body(
                     "The input GCA spans more than 180 degrees. Please divide it into two."
                 )
 
-            return in_between(GCRv0_lonlat[1], pt_lonlat[1], pole_lat) or in_between(
-                pole_lat, pt_lonlat[1], GCRv1_lonlat[1]
+            return in_between(gca_a_lonlat[1], pt_lonlat[1], pole_lat) or in_between(
+                pole_lat, pt_lonlat[1], gca_b_lonlat[1]
             )
-
+        # endIf
+        # ==============================================================================================================
     if is_directed:
         # The anti-meridian case Sufficient condition: absolute difference between the longitudes of the two
         # vertices is greater than 180 degrees (π radians): abs(GCRv1_lon - GCRv0_lon) > π
-        if abs(GCRv1_lonlat[0] - GCRv0_lonlat[0]) > np.pi:
+        if abs(gca_b_lonlat[0] - gca_a_lonlat[0]) > np.pi:
             # The necessary condition: the pt longitude is on the opposite side of the anti-meridian
             # Case 1: where 0 --> x0--> 180 -->x1 -->0 case is lager than the 180degrees (pi radians)
-            if GCRv0_lonlat[0] <= np.pi <= GCRv1_lonlat[0]:
+            if gca_a_lonlat[0] <= np.pi <= gca_b_lonlat[0]:
                 raise ValueError(
                     "The input Great Circle Arc span is larger than 180 degree, please break it into two."
                 )
 
             # The necessary condition: the pt longitude is on the opposite side of the anti-meridian
             # Case 2: The anti-meridian case where 180 -->x0 --> 0 lon --> x1 --> 180 lon
-            elif 2 * np.pi > GCRv0_lonlat[0] > np.pi > GCRv1_lonlat[0] > 0.0:
+            elif 2 * np.pi > gca_a_lonlat[0] > np.pi > gca_b_lonlat[0] > 0.0:
                 return in_between(
-                    GCRv0_lonlat[0], pt_lonlat[0], 2 * np.pi
-                ) or in_between(0, pt_lonlat[0], GCRv1_lonlat[0])
+                    gca_a_lonlat[0], pt_lonlat[0], 2 * np.pi
+                ) or in_between(0, pt_lonlat[0], gca_b_lonlat[0])
 
         # The non-anti-meridian case.
         else:
-            return in_between(GCRv0_lonlat[0], pt_lonlat[0], GCRv1_lonlat[0])
+            return in_between(gca_a_lonlat[0], pt_lonlat[0], gca_b_lonlat[0])
     else:
         # The undirected case
         # sort the longitude
-        GCRv0_lonlat_min, GCRv1_lonlat_max = sorted([GCRv0_lonlat[0], GCRv1_lonlat[0]])
+        GCRv0_lonlat_min, GCRv1_lonlat_max = sorted([gca_a_lonlat[0], gca_b_lonlat[0]])
         if np.pi > GCRv1_lonlat_max - GCRv0_lonlat_min >= 0.0:
-            return in_between(GCRv0_lonlat[0], pt_lonlat[0], GCRv1_lonlat[0])
+            return in_between(gca_a_lonlat[0], pt_lonlat[0], gca_b_lonlat[0])
         else:
             return in_between(GCRv1_lonlat_max, pt_lonlat[0], 2 * np.pi) or in_between(
                 0.0, pt_lonlat[0], GCRv0_lonlat_min
             )
 
-
-def point_within_gca(pt, gca_cart, is_directed=False):
-    """Check if a point lies on a given Great Circle Arc (GCA). The anti-
-    meridian case is also considered.
-
-    Parameters
-    ----------
-    pt : numpy.ndarray (float)
-        Cartesian coordinates of the point.
-    gca_cart : numpy.ndarray of shape (2, 3), (np.float or gmpy2.mpfr)
-        Cartesian coordinates of the Great Circle Arc (GCR).
-    is_directed : bool, optional, default = False
-        If True, the GCA is considered to be directed, which means it can only from v0-->v1. If False, the GCA is undirected,
-        and we will always assume the small circle (The one less than 180 degree) side is the GCA. The default is False.
-        For the case of the anti-podal case, the direction is v_0--> the pole point that on the same hemisphere as v_0-->v_1
-
-    Returns
-    -------
-    bool
-        True if the point lies between the two endpoints of the GCR, False otherwise.
-
-    Raises
-    ------
-    ValueError
-        If the input GCR spans exactly 180 degrees (π radians), as this GCR can have multiple planes.
-        In such cases, consider breaking the GCR into two separate GCRs.
-
-    ValueError
-        If the input GCR spans more than 180 degrees (π radians).
-        In such cases, consider breaking the GCR into two separate GCRs.
-
-    Notes
-    -----
-    The function checks if the given point is on the Great Circle Arc by considering its cartesian coordinates and
-    accounting for the anti-meridian case.
-
-    The anti-meridian case occurs when the GCR crosses the anti-meridian (0 longitude).
-    In this case, the function handles scenarios where the GCA spans across more than 180 degrees, requiring specific operation.
-
-    The function relies on the `_angle_of_2_vectors` and `is_between` functions to perform the necessary calculations.
-
-    Please ensure that the input coordinates are in radians and adhere to the ERROR_TOLERANCE value for floating-point comparisons.
-    """
-    # Convert the cartesian coordinates to lonlat coordinates
-    pt_lonlat = np.array(
-        _xyz_to_lonlat_rad_scalar(pt[0], pt[1], pt[2], normalize=False)
-    )
-    GCRv0_lonlat = np.array(
-        _xyz_to_lonlat_rad_scalar(
-            gca_cart[0][0], gca_cart[0][1], gca_cart[0][2], normalize=False
-        )
-    )
-    GCRv1_lonlat = np.array(
-        _xyz_to_lonlat_rad_scalar(
-            gca_cart[1][0], gca_cart[1][1], gca_cart[1][2], normalize=False
-        )
-    )
-    gca_cart = np.asarray(gca_cart)
-
-    # First if the input GCR is exactly 180 degree, we throw an exception, since this GCR can have multiple planes
-    angle = _angle_of_2_vectors(gca_cart[0], gca_cart[1])
-
-    out = _point_within_gca_body(
-        angle, gca_cart, pt, GCRv0_lonlat, GCRv1_lonlat, pt_lonlat, is_directed
-    )
-
-    return out
+    return None
 
 
 @njit(cache=True)
@@ -306,8 +337,52 @@ def _decide_pole_latitude(lat1, lat2):
     return closest_pole
 
 
-def extreme_gca_latitude(gca_cart, extreme_type):
-    """Calculate the maximum or minimum latitude of a great circle arc defined
+@njit(cache=True)
+def max3(a, b, c):
+    if a >= b and a >= c:
+        return a
+    elif b >= c:
+        return b
+    else:
+        return c
+
+
+@njit(cache=True)
+def min3(a, b, c):
+    if a <= b and a <= c:
+        return a
+    elif b <= c:
+        return b
+    else:
+        return c
+
+
+@njit(cache=True)
+def clip_scalar(a, a_min, a_max):
+    if a < a_min:
+        return a_min
+    elif a > a_max:
+        return a_max
+    else:
+        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):
+    """
+    Calculate the maximum or minimum latitude of a great circle arc defined
     by two 3D points.
 
     Parameters
@@ -315,6 +390,9 @@ def extreme_gca_latitude(gca_cart, extreme_type):
     gca_cart : numpy.ndarray
         An array containing two 3D vectors that define a great circle arc.
 
+    gca_lonlat : numpy.ndarray
+        An array containing the longitude and latitude of the two points.
+
     extreme_type : str
         The type of extreme latitude to calculate. Must be either 'max' or 'min'.
 
@@ -328,36 +406,57 @@ def extreme_gca_latitude(gca_cart, extreme_type):
     ValueError
         If `extreme_type` is not 'max' or 'min'.
     """
-    extreme_type = extreme_type.lower()
-
-    if extreme_type not in ("max", "min"):
+    # Validate extreme_type
+    if (extreme_type != "max") and (extreme_type != "min"):
         raise ValueError("extreme_type must be either 'max' or 'min'")
 
-    n1, n2 = gca_cart
+    # Extract the two points
+    n1 = gca_cart[0]
+    n2 = gca_cart[1]
+
+    # Compute dot product
+    dot_n1_n2 = dot(n1, n2)
 
-    dot_n1_n2 = dot(np.asarray(n1), np.asarray(n2))
+    # Compute denominator
     denom = (n1[2] + n2[2]) * (dot_n1_n2 - 1.0)
-    d_a_max = (n1[2] * dot_n1_n2 - n2[2]) / denom
 
-    d_a_max = (
-        np.clip(d_a_max, 0, 1)
-        if isclose(d_a_max, 0, atol=ERROR_TOLERANCE)
-        or isclose(d_a_max, 1, atol=ERROR_TOLERANCE)
-        else d_a_max
-    )
-    # Before we make sure the grid coordinates are normalized, do not try to skip the normalization steps!
-    _, lat_n1 = _xyz_to_lonlat_rad_scalar(n1[0], n1[1], n1[2], normalize=True)
-    _, lat_n2 = _xyz_to_lonlat_rad_scalar(n2[0], n2[1], n2[2], normalize=True)
-
-    if 0 < d_a_max < 1:
-        node3 = (1 - d_a_max) * n1 + d_a_max * n2
-        node3 = np.array(_normalize_xyz_scalar(node3[0], node3[1], node3[2]))
-        d_lat_rad = np.arcsin(np.clip(node3[2], -1, 1))
-
-        return (
-            max(d_lat_rad, lat_n1, lat_n2)
-            if extreme_type == "max"
-            else min(d_lat_rad, lat_n1, lat_n2)
-        )
+    # Initialize latitudes
+    lon_n1, lat_n1 = gca_lonlat[0]
+    lon_n2, lat_n2 = gca_lonlat[1]
+
+    # Check if 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
+
+            # Compute latitude of the intermediate point
+            d_lat_rad = math.asin(clip_scalar(node3_normalized[2], -1.0, 1.0))
+
+            # Return the extreme latitude
+            if extreme_type == "max":
+                return max3(d_lat_rad, lat_n1, lat_n2)
+            else:
+                return min3(d_lat_rad, lat_n1, lat_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(lat_n1, lat_n2)
     else:
-        return max(lat_n1, lat_n2) if extreme_type == "max" else min(lat_n1, lat_n2)
+        return min(lat_n1, lat_n2)
diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py
index cb3b2cad..25e1d010 100644
--- a/uxarray/grid/coordinates.py
+++ b/uxarray/grid/coordinates.py
@@ -1,6 +1,7 @@
 import xarray as xr
 import numpy as np
 import warnings
+import math
 
 from uxarray.conventions import ugrid
 
@@ -67,36 +68,10 @@ def _xyz_to_lonlat_rad_no_norm(
 
 
 @njit(cache=True)
-def _xyz_to_lonlat_rad_scalar(
-    x: Union[np.ndarray, float],
-    y: Union[np.ndarray, float],
-    z: Union[np.ndarray, float],
-    normalize: bool = True,
-):
-    """Converts a Cartesian x,y,z coordinates into Spherical latitude and
-    longitude without normalization, decorated with Numba.
-
-    Parameters
-    ----------
-    x : float
-        Cartesian x coordinate
-    y: float
-        Cartesiain y coordinate
-    z: float
-        Cartesian z coordinate
-
-
-    Returns
-    -------
-    lon : float
-        Longitude in radians
-    lat: float
-        Latitude in radians
-    """
-
+def _xyz_to_lonlat_rad_scalar(x, y, z, normalize=True):
     if normalize:
         x, y, z = _normalize_xyz_scalar(x, y, z)
-        denom = np.abs(x * x + y * y + z * z)
+        denom = abs(x * x + y * y + z * z)
         x /= denom
         y /= denom
         z /= denom
@@ -104,15 +79,19 @@ def _xyz_to_lonlat_rad_scalar(
     lon = np.atan2(y, x)
     lat = np.asin(z)
 
-    # set longitude range to [0, pi]
-    lon = np.mod(lon, 2 * np.pi)
+    # Set longitude range to [0, 2*pi]
+    lon = lon % (2 * math.pi)
 
-    z_mask = np.abs(z) > 1.0 - ERROR_TOLERANCE
+    z_abs = abs(z)
+    if z_abs > (1.0 - ERROR_TOLERANCE):
+        lat = math.copysign(math.pi / 2, z)
+        lon = 0.0
 
-    lat = np.where(z_mask, np.sign(z) * np.pi / 2, lat)
-    lon = np.where(z_mask, 0.0, lon)
+    lonlat = np.empty(2)
+    lonlat[0] = lon
+    lonlat[1] = lat
 
-    return lon, lat
+    return lonlat
 
 
 def _xyz_to_lonlat_rad(
diff --git a/uxarray/grid/geometry.py b/uxarray/grid/geometry.py
index ed564aaf..8f05ee66 100644
--- a/uxarray/grid/geometry.py
+++ b/uxarray/grid/geometry.py
@@ -2,7 +2,7 @@ import antimeridian
 import cartopy.crs as ccrs
 import geopandas
 from matplotlib.collections import LineCollection, PolyCollection
-from numba import njit
+from numba import njit, prange
 import numpy as np
 import pandas as pd
 import shapely
@@ -11,6 +11,7 @@ from shapely import polygons as Polygons
 import spatialpandas
 from spatialpandas.geometry import MultiPolygonArray, PolygonArray
 import xarray as xr
+import math
 
 from uxarray.constants import (
     ERROR_TOLERANCE,
@@ -18,6 +19,7 @@ from uxarray.constants import (
     INT_FILL_VALUE,
 )
 from uxarray.grid.arcs import extreme_gca_latitude, point_within_gca
+from uxarray.grid.coordinates import _xyz_to_lonlat_rad
 from uxarray.grid.intersections import gca_gca_intersection
 from uxarray.grid.utils import (
     _get_cartesian_face_edge_nodes,
@@ -25,62 +27,67 @@ from uxarray.grid.utils import (
 )
 from uxarray.utils.computing import allclose, isclose
 
-POLE_POINTS = {"North": np.array([0.0, 0.0, 1.0]), "South": np.array([0.0, 0.0, -1.0])}
+POLE_POINTS_XYZ = {
+    "North": np.array([0.0, 0.0, 1.0]),
+    "South": np.array([0.0, 0.0, -1.0]),
+}
+POLE_POINTS_LONLAT = {
+    "North": np.array([0.0, np.pi / 2]),
+    "South": np.array([0.0, -np.pi / 2]),
+}
 
 # number of faces/polygons before raising a warning for performance
 GDF_POLYGON_THRESHOLD = 100000
 
-REFERENCE_POINT_EQUATOR = np.array([1.0, 0.0, 0.0])
+REFERENCE_POINT_EQUATOR_XYZ = np.array([1.0, 0.0, 0.0])
+REFERENCE_POINT_EQUATOR_LONLAT = np.array([0.0, 0.0])  # TODO ?
 
+POLE_NAME_TO_INT = {"North": 1, "Equator": 0, "South": -1}
 
+
+@njit(cache=True)
+def error_radius(p1, p2):
+    """Calculate the error radius between two points in 3D space."""
+    numerator = math.sqrt(
+        (p1[0] - p2[0]) ** 2 + (p1[1] - p2[1]) ** 2 + (p1[2] - p2[2]) ** 2
+    )
+    denominator = math.sqrt(p2[0] ** 2 + p2[1] ** 2 + p2[2] ** 2)
+    return numerator / denominator
+
+
+@njit(cache=True)
 def _unique_points(points, tolerance=ERROR_TOLERANCE):
     """Identify unique intersection points from a list of points, considering floating point precision errors.
 
     Parameters
     ----------
-    points : list of array-like
-        A list containing the intersection points, where each point is an array-like structure (e.g., list, tuple, or numpy array) containing x, y, and z coordinates.
-    tolerance : float, optional
-        The distance threshold within which two points are considered identical. Default is 1e-6.
+    points : np.ndarray
+        An array of shape (n_points, 3) containing the intersection points.
+    tolerance : float
+        The distance threshold within which two points are considered identical.
 
     Returns
     -------
-    list of np.ndarray
-        A list of unique snapped points in Cartesian coordinates.
-
-    Notes
-    -----
-    Given the nature of the mathematical equations and the spherical surface, it is more reasonable to calculate the "error radius" of the results using the following formula.
-    In the equation below, \(\tilde{P}_x\) and \(\tilde{P}_y\) are the calculated results, and \(P_x\) and \(P_y\) are the actual intersection points for the \(x\) and \(y\) coordinates, respectively.
-    The \(z\) coordinate is always the input \(z_0\), the constant latitude, so it is not included in this error calculation.
-
-    .. math::
-        \begin{aligned}
-            &\frac{\sqrt{(\tilde{v}_x - v_x)^2 + (\tilde{v}_y - v_y)^2 + (\tilde{v}_z - v_z)^2}}{\sqrt{v_x^2 + v_y^2 + v_z^2}}\\
-            &= \sqrt{(\tilde{v}_x - v_x)^2 + (\tilde{v}_y - v_y)^2 + (\tilde{v}_z - v_z)^2}
-        \end{aligned}
-
-    This method ensures that small numerical inaccuracies do not lead to multiple close points being considered different.
+    np.ndarray
+        An array of unique points in Cartesian coordinates.
     """
-    unique_points = []
-    points = [np.array(point) for point in points]  # Ensure all points are numpy arrays
-
-    def error_radius(p1, p2):
-        """Calculate the error radius between two points in 3D space."""
-        numerator = np.sqrt(
-            (p1[0] - p2[0]) ** 2 + (p1[1] - p2[1]) ** 2 + (p1[2] - p2[2]) ** 2
-        )
-        denominator = np.sqrt(p2[0] ** 2 + p2[1] ** 2 + p2[2] ** 2)
-        return numerator / denominator
+    n_points = points.shape[0]
+    unique_points = np.empty((n_points, 3), dtype=points.dtype)
+    unique_count = 0
 
-    for point in points:
-        if not any(
-            error_radius(point, unique_point) < tolerance
-            for unique_point in unique_points
-        ):
-            unique_points.append(point)
+    for i in range(n_points):
+        point = points[i]
+        is_unique = True
+        for j in range(unique_count):
+            unique_point = unique_points[j]
+            if error_radius(point, unique_point) < tolerance:
+                is_unique = False
+                break
+        if is_unique:
+            unique_points[unique_count] = point
+            unique_count += 1
 
-    return unique_points
+    return unique_points[:unique_count]
 
 
 @njit(cache=True)
@@ -670,145 +677,250 @@ def _convert_shells_to_polygons(shells):
     return polygons
 
 
-def _pole_point_inside_polygon(pole, face_edge_cart):
-    """Determines if a pole point is inside a polygon.
+def _pole_point_inside_polygon_cartesian(pole, face_edges_xyz):
+    if isinstance(pole, str):
+        pole = POLE_NAME_TO_INT[pole]
 
-    .. note::
-        - If the pole point is on the edge of the polygon, it will be considered "inside the polygon".
+    x = face_edges_xyz[:, :, 0]
+    y = face_edges_xyz[:, :, 1]
+    z = face_edges_xyz[:, :, 2]
 
-    Parameters
-    ----------
-    pole : str
-        Either 'North' or 'South'.
-    face_edge_cart : np.ndarray
-        A face polygon represented by edges in Cartesian coordinates. Shape: (n_edges, 2, 3)
+    lon, lat = _xyz_to_lonlat_rad(x, y, z)
 
-    Returns
-    -------
-    bool
-        True if pole point is inside polygon, False otherwise.
+    face_edges_lonlat = np.stack((lon, lat), axis=2)
 
-    Raises
-    ------
-    ValueError
-        If the provided pole is neither 'North' nor 'South'.
+    return pole_point_inside_polygon(pole, face_edges_xyz, face_edges_lonlat)
 
-    Warning
-    -------
-    UserWarning
-        Raised if the face contains both pole points.
-    """
-    if pole not in POLE_POINTS:
-        raise ValueError('Pole point must be either "North" or "South"')
+    pass
 
-    # Classify the polygon's location
-    location = _classify_polygon_location(face_edge_cart)
-    pole_point = POLE_POINTS[pole]
 
-    if location == pole:
-        ref_edge = np.array([pole_point, REFERENCE_POINT_EQUATOR])
-        return _check_intersection(ref_edge, face_edge_cart) % 2 != 0
-    elif location == "Equator":
-        # smallest offset I can obtain when using the float64 type
+@njit(cache=True)
+def pole_point_inside_polygon(pole, face_edges_xyz, face_edges_lonlat):
+    """Determines if a pole point is inside a polygon."""
 
-        ref_edge_north = np.array([pole_point, REFERENCE_POINT_EQUATOR])
-        ref_edge_south = np.array([-pole_point, REFERENCE_POINT_EQUATOR])
+    if pole != 1 and pole != -1:
+        raise ValueError("Pole must be 1 (North) or -1 (South)")
 
-        north_edges = face_edge_cart[np.any(face_edge_cart[:, :, 2] > 0, axis=1)]
-        south_edges = face_edge_cart[
-            ~np.isin(face_edge_cart, north_edges).all(axis=(1, 2))
-        ]
-        # The equator one is assigned to the south edges
-        return (
-            _check_intersection(ref_edge_north, north_edges)
-            + _check_intersection(ref_edge_south, south_edges)
-        ) % 2 != 0
-    elif (
-        location == "North"
-        and pole == "South"
-        or location == "South"
-        and pole == "North"
-    ):
+    # Define constants within the function
+    pole_point_xyz = np.empty(3, dtype=np.float64)
+    pole_point_xyz[0] = 0.0
+    pole_point_xyz[1] = 0.0
+    pole_point_xyz[2] = 1.0 * pole
+
+    pole_point_lonlat = np.empty(2, dtype=np.float64)
+    pole_point_lonlat[0] = 0.0
+    pole_point_lonlat[1] = (math.pi / 2) * pole
+
+    REFERENCE_POINT_EQUATOR_XYZ = np.empty(3, dtype=np.float64)
+    REFERENCE_POINT_EQUATOR_XYZ[0] = 1.0
+    REFERENCE_POINT_EQUATOR_XYZ[1] = 0.0
+    REFERENCE_POINT_EQUATOR_XYZ[2] = 0.0
+
+    REFERENCE_POINT_EQUATOR_LONLAT = np.empty(2, dtype=np.float64)
+    REFERENCE_POINT_EQUATOR_LONLAT[0] = 0.0
+    REFERENCE_POINT_EQUATOR_LONLAT[1] = 0.0
+
+    # Classify the polygon's location
+    location = _classify_polygon_location(face_edges_xyz)
+
+    if (location == 1 and pole == -1) or (location == -1 and pole == 1):
         return False
+
+    elif location == -1 or location == 1:
+        # Initialize ref_edge_xyz
+        ref_edge_xyz = np.empty((2, 3), dtype=np.float64)
+        ref_edge_xyz[0, 0] = pole_point_xyz[0]
+        ref_edge_xyz[0, 1] = pole_point_xyz[1]
+        ref_edge_xyz[0, 2] = pole_point_xyz[2]
+        ref_edge_xyz[1, :] = REFERENCE_POINT_EQUATOR_XYZ
+
+        # Initialize ref_edge_lonlat
+        ref_edge_lonlat = np.empty((2, 2), dtype=np.float64)
+        ref_edge_lonlat[0, 0] = pole_point_lonlat[0]
+        ref_edge_lonlat[0, 1] = pole_point_lonlat[1]
+        ref_edge_lonlat[1, :] = REFERENCE_POINT_EQUATOR_LONLAT
+
+        intersection_count = _check_intersection(
+            ref_edge_xyz, ref_edge_lonlat, face_edges_xyz, face_edges_lonlat
+        )
+        return (intersection_count % 2) != 0
+
+    elif location == 0:  # Equator
+        # Initialize ref_edge_north_xyz and ref_edge_north_lonlat
+        ref_edge_north_xyz = np.empty((2, 3), dtype=np.float64)
+        ref_edge_north_xyz[0, 0] = 0.0
+        ref_edge_north_xyz[0, 1] = 0.0
+        ref_edge_north_xyz[0, 2] = 1.0
+        ref_edge_north_xyz[1, :] = REFERENCE_POINT_EQUATOR_XYZ
+
+        ref_edge_north_lonlat = np.empty((2, 2), dtype=np.float64)
+        ref_edge_north_lonlat[0, 0] = 0.0
+        ref_edge_north_lonlat[0, 1] = math.pi / 2
+        ref_edge_north_lonlat[1, :] = REFERENCE_POINT_EQUATOR_LONLAT
+
+        # Initialize ref_edge_south_xyz and ref_edge_south_lonlat
+        ref_edge_south_xyz = np.empty((2, 3), dtype=np.float64)
+        ref_edge_south_xyz[0, 0] = 0.0
+        ref_edge_south_xyz[0, 1] = 0.0
+        ref_edge_south_xyz[0, 2] = -1.0
+        ref_edge_south_xyz[1, :] = REFERENCE_POINT_EQUATOR_XYZ
+
+        ref_edge_south_lonlat = np.empty((2, 2), dtype=np.float64)
+        ref_edge_south_lonlat[0, 0] = 0.0
+        ref_edge_south_lonlat[0, 1] = -math.pi / 2
+        ref_edge_south_lonlat[1, :] = REFERENCE_POINT_EQUATOR_LONLAT
+
+        # Classify edges based on z-coordinate
+        n_edges = face_edges_xyz.shape[0]
+        north_edges_xyz = np.empty((n_edges, 2, 3), dtype=np.float64)
+        north_edges_lonlat = np.empty((n_edges, 2, 2), dtype=np.float64)
+        south_edges_xyz = np.empty((n_edges, 2, 3), dtype=np.float64)
+        south_edges_lonlat = np.empty((n_edges, 2, 2), dtype=np.float64)
+        north_count = 0
+        south_count = 0
+
+        for i in range(n_edges):
+            edge_xyz = face_edges_xyz[i]
+            edge_lonlat = face_edges_lonlat[i]
+            if edge_xyz[0, 2] > 0 or edge_xyz[1, 2] > 0:
+                north_edges_xyz[north_count] = edge_xyz
+                north_edges_lonlat[north_count] = edge_lonlat
+                north_count += 1
+            elif edge_xyz[0, 2] < 0 or edge_xyz[1, 2] < 0:
+                south_edges_xyz[south_count] = edge_xyz
+                south_edges_lonlat[south_count] = edge_lonlat
+                south_count += 1
+            else:
+                # skip edges exactly on the equator
+                continue
+
+        if north_count > 0:
+            north_edges_xyz = north_edges_xyz[:north_count]
+            north_edges_lonlat = north_edges_lonlat[:north_count]
+        else:
+            north_edges_xyz = np.empty((0, 2, 3), dtype=np.float64)
+            north_edges_lonlat = np.empty((0, 2, 2), dtype=np.float64)
+
+        if south_count > 0:
+            south_edges_xyz = south_edges_xyz[:south_count]
+            south_edges_lonlat = south_edges_lonlat[:south_count]
+        else:
+            south_edges_xyz = np.empty((0, 2, 3), dtype=np.float64)
+            south_edges_lonlat = np.empty((0, 2, 2), dtype=np.float64)
+
+        # Count south intersections
+        north_intersections = _check_intersection(
+            ref_edge_north_xyz,
+            ref_edge_north_lonlat,
+            north_edges_xyz,
+            north_edges_lonlat,
+        )
+
+        # Count south intersections
+        south_intersections = _check_intersection(
+            ref_edge_south_xyz,
+            ref_edge_south_lonlat,
+            south_edges_xyz,
+            south_edges_lonlat,
+        )
+
+        return ((north_intersections + south_intersections) % 2) != 0
+
     else:
         raise ValueError(
-            "Invalid pole point query. Current location: {}, query pole point: {}".format(
-                location, pole
-            )
+            f"Invalid pole point query. Current location: {location}, query pole point: {pole}"
         )
 
 
-def _check_intersection(ref_edge, edges):
-    """Check the number of intersections of the reference edge with the given
-    edges.
+@njit(cache=True)
+def _classify_polygon_location(face_edge_cart):
+    """Classify the location of the polygon relative to the hemisphere."""
+    z_coords = face_edge_cart[:, :, 2]
+    if np.all(z_coords > 0):  # Use strict inequality
+        return 1  # North
+    elif np.all(z_coords < 0):  # Use strict inequality
+        return -1  # South
+    else:
+        return 0  # Equator
+
+
+@njit(cache=True)
+def _check_intersection(ref_edge_xyz, ref_edge_lonlat, edges_xyz, edges_lonlat):
+    """Check the number of intersections of the reference edge with the given edges.
 
     Parameters
     ----------
-    ref_edge : np.ndarray
-        Reference edge to check intersections against.
-    edges : np.ndarray
+    ref_edge_xyz : np.ndarray
+        Reference edge to check intersections against. Shape: (2, 3)
+    ref_edge_lonlat : np.ndarray
+        Reference edge longitude and latitude. Shape: (2, 2)
+    edges_xyz : np.ndarray
         Edges to check for intersections. Shape: (n_edges, 2, 3)
+    edges_lonlat : np.ndarray
+        Longitude and latitude of the edges. Shape: (n_edges, 2, 2)
 
     Returns
     -------
     int
         Count of intersections.
     """
-    pole_point, ref_point = ref_edge
-    intersection_points = []
+    pole_point_xyz = ref_edge_xyz[0]
+    n_edges = edges_xyz.shape[0]
 
-    for edge in edges:
-        intersection_point = gca_gca_intersection(ref_edge, edge)
+    # Assuming at most 2 intersections per edge
+    max_intersections = n_edges * 2
+    intersection_points = np.empty((max_intersections, 3), dtype=np.float64)
+    intersection_count = 0
 
-        if intersection_point.size != 0:
-            # Handle both single point and multiple points case
-            if (
-                intersection_point.ndim == 1
-            ):  # If there's only one point, make it a 2D array
-                intersection_point = [intersection_point]  # Convert to list of points
+    for i in range(n_edges):
+        edge_xyz = edges_xyz[i]
+        edge_lonlat = edges_lonlat[i]
 
-            # for each intersection point, check if it is a pole point
-            for point in intersection_point:
-                if allclose(point, pole_point, atol=ERROR_TOLERANCE):
-                    return True
-                intersection_points.append(point)
+        # compute intersection
+        intersection_point = gca_gca_intersection(
+            ref_edge_xyz, ref_edge_lonlat, edge_xyz, edge_lonlat
+        )
 
-    # Only return the unique intersection points, the unique tolerance is set to ERROR_TOLERANCE
-    intersection_points = _unique_points(intersection_points, tolerance=ERROR_TOLERANCE)
+        if intersection_point.size != 0:
+            if intersection_point.ndim == 1:
+                # Only one point
+                point = intersection_point
+                if allclose(point, pole_point_xyz, atol=ERROR_TOLERANCE):
+                    return True
+                intersection_points[intersection_count] = point
+                intersection_count += 1
+            else:
+                # Multiple points
+                num_points = intersection_point.shape[0]
+                for j in range(num_points):
+                    point = intersection_point[j]
+                    if allclose(point, pole_point_xyz, atol=ERROR_TOLERANCE):
+                        return True
+                    intersection_points[intersection_count] = point
+                    intersection_count += 1
+
+    if intersection_count == 0:
+        return 0
+
+    unique_intersection_points = _unique_points(
+        intersection_points[:intersection_count], tolerance=ERROR_TOLERANCE
+    )
+    unique_count = unique_intersection_points.shape[0]
 
-    # If the unique intersection point is one and it is exactly one of the nodes of the face, return 0
-    if len(intersection_points) == 1:
-        for edge in edges:
+    # If there's only one unique intersection point, check if it matches any edge nodes
+    if unique_count == 1:
+        intersection_point = unique_intersection_points[0]
+        for i in range(n_edges):
+            edge_xyz = edges_xyz[i]
             if allclose(
-                intersection_points[0], edge[0], atol=ERROR_TOLERANCE
-            ) or allclose(intersection_points[0], edge[1], atol=ERROR_TOLERANCE):
+                intersection_point, edge_xyz[0], atol=ERROR_TOLERANCE
+            ) or allclose(intersection_point, edge_xyz[1], atol=ERROR_TOLERANCE):
                 return 0
 
-    return len(intersection_points)
-
-
-def _classify_polygon_location(face_edge_cart):
-    """Classify the location of the polygon relative to the hemisphere.
-
-    Parameters
-    ----------
-    face_edge_cart : np.ndarray
-        A face polygon represented by edges in Cartesian coordinates. Shape: (n_edges, 2, 3)
-
-    Returns
-    -------
-    str
-        Returns either 'North', 'South' or 'Equator' based on the polygon's location.
-    """
-    z_coords = face_edge_cart[:, :, 2]
-    if np.all(z_coords > 0):
-        return "North"
-    elif np.all(z_coords < 0):
-        return "South"
-    else:
-        return "Equator"
+    return unique_count
 
 
+@njit(cache=True)
 def _get_latlonbox_width(latlonbox_rad):
     """Calculate the width of a latitude-longitude box in radians. The box
     should be represented by a 2x2 array in radians and lon0 represent the
@@ -847,18 +959,16 @@ def _get_latlonbox_width(latlonbox_rad):
         If the input longitude range is flagged as periodic but in the form [lon0, lon1] where lon0 < lon1.
         The function will automatically use the is_lon_periodic=False instead.
     """
-
+    # TODO: ensure this function works properly
     lon0, lon1 = latlonbox_rad[1]
 
-    # Check longitude range validity
-    # Normalize the longitude so that it is within the range [0, 2π]
-    # Normalize the longitude
     if lon0 != INT_FILL_VALUE:
         lon0 = np.mod(lon0, 2 * np.pi)
     if lon1 != INT_FILL_VALUE:
         lon1 = np.mod(lon1, 2 * np.pi)
     if (lon0 < 0.0 or lon0 > 2.0 * np.pi) and lon0 != INT_FILL_VALUE:
-        raise Exception("lon0 out of range ({} not in [0, 2π])".format(lon0))
+        # -1 used for exception
+        return -1
 
     if lon0 <= lon1:
         return lon1 - lon0
@@ -867,18 +977,33 @@ def _get_latlonbox_width(latlonbox_rad):
         return 2 * np.pi - lon0 + lon1
 
 
-def _insert_pt_in_latlonbox(old_box, new_pt, is_lon_periodic=True):
-    """Update the latitude-longitude box to include a new point in radians.
+@njit(cache=True)
+def all_elements_nan(arr):
+    """Check if all elements in an array are np.nan."""
+    for i in range(arr.shape[0]):
+        if not np.isnan(arr[i]):
+            return False
+    return True
 
-    This function compares the new point's latitude and longitude with the
-    existing latitude-longitude box and updates the box if necessary to include the new point.
+
+@njit(cache=True)
+def any_close_lat(lat_pt, atol):
+    """Check if the latitude point is close to either the North or South Pole."""
+    return isclose(lat_pt, 0.5 * math.pi, atol) or isclose(lat_pt, -0.5 * math.pi, atol)
+
+
+@njit(cache=True)
+def insert_pt_in_latlonbox(old_box, new_pt, is_lon_periodic=True):
+    """Update the latitude-longitude box to include a new point in radians.
 
     Parameters
     ----------
     old_box : np.ndarray
-        The original latitude-longitude box in radian, a 2x2 array: [min_lat, max_lat],[left_lon, right_lon]].
+        The original latitude-longitude box in radians, a 2x2 array:
+        [[min_lat, max_lat],
+         [left_lon, right_lon]].
     new_pt : np.ndarray
-        The new latitude-longitude point in radian, an array: [lat, lon].
+        The new latitude-longitude point in radians, an array: [lat, lon].
     is_lon_periodic : bool, optional
         Flag indicating if the latitude-longitude box is periodic in longitude (default is True).
 
@@ -891,301 +1016,340 @@ def _insert_pt_in_latlonbox(old_box, new_pt, is_lon_periodic=True):
     ------
     Exception
         If logic errors occur in the calculation process.
-
-    Examples
-    --------
-    >>> _insert_pt_in_latlonbox(
-    ...     np.array([[1.0, 2.0], [3.0, 4.0]]), np.array([1.5, 3.5])
-    ... )
-    array([[1.0, 2.0], [3.0, 4.0]])
     """
-    if np.all(new_pt == INT_FILL_VALUE):
+    # Check if the new point is a fill value
+    all_fill = all_elements_nan(new_pt)
+    if all_fill:
         return old_box
 
-    latlon_box = np.copy(old_box)  # Create a copy of the old box
-    latlon_box = np.array(
-        latlon_box, dtype=np.float64
-    )  # Cast to float64, otherwise the following update might fail
+    # Create a copy of the old box
+    latlon_box = np.copy(old_box)
 
-    lat_pt, lon_pt = new_pt
+    # Extract latitude and longitude from the new point
+    lat_pt = new_pt[0]
+    lon_pt = new_pt[1]
 
-    # Normalize the longitude
-    if lon_pt != INT_FILL_VALUE:
-        lon_pt = np.mod(lon_pt, 2 * np.pi)
+    # Normalize the longitude if it's not a fill value
+    if not np.isnan(lon_pt):
+        lon_pt = lon_pt % (2.0 * math.pi)
 
     # Check if the latitude range is uninitialized and update it
-    if old_box[0][0] == old_box[0][1] == INT_FILL_VALUE:
-        latlon_box[0] = np.array([lat_pt, lat_pt])
-
-    # Check if the longitude range is uninitialized and update it
-    if old_box[1][0] == old_box[1][1] == INT_FILL_VALUE:
-        latlon_box[1] = np.array([lon_pt, lon_pt])
-
-    if lon_pt != INT_FILL_VALUE and (lon_pt < 0.0 or lon_pt > 2.0 * np.pi):
-        raise Exception(f"lon_pt out of range ({lon_pt} not in [0, 2π])")
-
-    # Check for pole points and update latitudes
-    is_pole_point = (
-        lon_pt == INT_FILL_VALUE
-        and isclose(
-            new_pt[0], np.asarray([0.5 * np.pi, -0.5 * np.pi]), atol=ERROR_TOLERANCE
-        ).any()
-    )
-
-    if is_pole_point:
-        # Check if the new point is close to the North Pole
-        if isclose(new_pt[0], 0.5 * np.pi, atol=ERROR_TOLERANCE):
-            latlon_box[0][1] = 0.5 * np.pi
-
-        # Check if the new point is close to the South Pole
-        elif isclose(new_pt[0], -0.5 * np.pi, atol=ERROR_TOLERANCE):
-            latlon_box[0][0] = -0.5 * np.pi
-
-        return latlon_box
+    if np.isnan(old_box[0, 0]) and np.isnan(old_box[0, 1]):
+        latlon_box[0, 0] = lat_pt
+        latlon_box[0, 1] = lat_pt
     else:
-        latlon_box[0] = [min(latlon_box[0][0], lat_pt), max(latlon_box[0][1], lat_pt)]
+        # Update latitude range
+        if lat_pt < latlon_box[0, 0]:
+            latlon_box[0, 0] = lat_pt
+        if lat_pt > latlon_box[0, 1]:
+            latlon_box[0, 1] = lat_pt
 
-    # Update longitude range for non-periodic or periodic cases
-    if not is_lon_periodic:
-        latlon_box[1] = [min(latlon_box[1][0], lon_pt), max(latlon_box[1][1], lon_pt)]
+    # Check if the longitude range is uninitialized and update it
+    if np.isnan(old_box[1, 0]) and np.isnan(old_box[1, 1]):
+        latlon_box[1, 0] = lon_pt
+        latlon_box[1, 1] = lon_pt
     else:
-        if (
-            latlon_box[1][0] > latlon_box[1][1]
-            and (lon_pt < latlon_box[1][0] and lon_pt > latlon_box[1][1])
-        ) or (
-            latlon_box[1][0] <= latlon_box[1][1]
-            and not (latlon_box[1][0] <= lon_pt <= latlon_box[1][1])
-        ):
-            # Calculate and compare new box widths
-            box_a, box_b = np.copy(latlon_box), np.copy(latlon_box)
-            box_a[1][0], box_b[1][1] = lon_pt, lon_pt
-            d_width_a, d_width_b = (
-                _get_latlonbox_width(box_a),
-                _get_latlonbox_width(box_b),
-            )
-
-            # The width should not be negative, if so, raise an exception
-            if d_width_a < 0 or d_width_b < 0:
-                raise Exception("logic error")
-
-            # Return the arc with the smaller width
-            latlon_box = box_a if d_width_a < d_width_b else box_b
+        # Validate longitude point
+        if not np.isnan(lon_pt) and (lon_pt < 0.0 or lon_pt > 2.0 * math.pi):
+            raise Exception("Longitude point out of range")
+
+        # Check for pole points
+        is_pole_point = False
+        if np.isnan(lon_pt):
+            if any_close_lat(lat_pt, ERROR_TOLERANCE):
+                is_pole_point = True
+
+        if is_pole_point:
+            # Update latitude for pole points
+            if isclose(lat_pt, 0.5 * math.pi, ERROR_TOLERANCE):
+                latlon_box[0, 1] = 0.5 * math.pi  # Update max_lat for North Pole
+            elif isclose(lat_pt, -0.5 * math.pi, ERROR_TOLERANCE):
+                latlon_box[0, 0] = -0.5 * math.pi  # Update min_lat for South Pole
+        else:
+            # Update longitude range based on periodicity
+            if not is_lon_periodic:
+                # Non-periodic: straightforward min and max updates
+                if lon_pt < latlon_box[1, 0]:
+                    latlon_box[1, 0] = lon_pt
+                if lon_pt > latlon_box[1, 1]:
+                    latlon_box[1, 1] = lon_pt
+            else:
+                # Periodic longitude handling
+                # Determine if the new point extends the current longitude range
+                # considering the periodic boundary at 2*pi
+                left_lon = latlon_box[1, 0]
+                right_lon = latlon_box[1, 1]
+
+                # Check if the current box wraps around
+                wraps_around = left_lon > right_lon
+
+                if wraps_around:
+                    # If the box wraps around, check if the new point is outside the current range
+                    if not (left_lon <= lon_pt or lon_pt <= right_lon):
+                        # Decide to extend either the left or the right
+                        # Calculate the new width for both possibilities
+                        # Option 1: Extend the left boundary to lon_pt
+                        box_a = np.copy(latlon_box)
+                        box_a[1, 0] = lon_pt
+                        d_width_a = _get_latlonbox_width(box_a)
+
+                        # Option 2: Extend the right boundary to lon_pt
+                        box_b = np.copy(latlon_box)
+                        box_b[1, 1] = lon_pt
+                        d_width_b = _get_latlonbox_width(box_b)
+
+                        # Ensure widths are non-negative
+                        if (d_width_a < 0.0) or (d_width_b < 0.0):
+                            raise Exception(
+                                "Logic error in longitude box width calculation"
+                            )
+
+                        # Choose the box with the smaller width
+                        if d_width_a < d_width_b:
+                            latlon_box = box_a
+                        else:
+                            latlon_box = box_b
+                else:
+                    # If the box does not wrap around, simply update min or max longitude
+                    if lon_pt < left_lon or lon_pt > right_lon:
+                        # Calculate the new width for both possibilities
+                        # Option 1: Extend the left boundary to lon_pt
+                        box_a = np.copy(latlon_box)
+                        box_a[1, 0] = lon_pt
+                        d_width_a = _get_latlonbox_width(box_a)
+
+                        # Option 2: Extend the right boundary to lon_pt
+                        box_b = np.copy(latlon_box)
+                        box_b[1, 1] = lon_pt
+                        d_width_b = _get_latlonbox_width(box_b)
+
+                        # Ensure widths are non-negative
+                        if (d_width_a < 0.0) or (d_width_b < 0.0):
+                            raise Exception(
+                                "Logic error in longitude box width calculation"
+                            )
+
+                        # Choose the box with the smaller width
+                        if d_width_a < d_width_b:
+                            latlon_box = box_a
+                        else:
+                            latlon_box = box_b
 
     return latlon_box
 
 
+@njit(cache=True)
 def _populate_face_latlon_bound(
-    face_edges_cartesian,
-    face_edges_lonlat_rad,
-    is_latlonface: bool = False,
+    face_edges_xyz,
+    face_edges_lonlat,
+    is_latlonface=False,
     is_GCA_list=None,
 ):
-    """Populates the bounding box for each face in the grid by evaluating the
-    geographical bounds based on the Cartesian and latitudinal/longitudinal
-    edge connectivity. This function also considers the presence of pole points
-    within the face's bounds and adjusts the bounding box accordingly.
-
-    Parameters
-    ----------
-    face_edges_cartesian : np.ndarray, shape (n_edges, 2, 3)
-        An array holding the Cartesian coordinates for the edges of a face, where `n_edges`
-        is the number of edges for the specific face. Each edge is represented by two points
-        (start and end), and each point is a 3D vector (x, y, z) in Cartesian coordinates.
-
-    face_edges_lonlat_connectivity_rad : np.ndarray, shape (n_edges, 2, 2)
-        An array holding the longitude and latitude in radians for the edges of a face,
-        formatted similarly to `face_edges_cartesian`. Each edge's start and
-        end points are represented by their longitude and latitude values in radians.
-
-    is_latlonface : bool, optional
-        A flag indicating if the current face is a latitudinal/longitudinal (latlon) face,
-        meaning its edges align with lines of constant latitude or longitude. If `True`,
-        edges are treated as following constant latitudinal or longitudinal lines. If `False`,
-        edges are considered as great circle arcs (GCA). Default is `False`.
-
-    is_GCA_list : list or np.ndarray, optional
-        A list or an array of boolean values corresponding to each edge within the face,
-        indicating whether the edge is a GCA. If `None`, the determination of whether an
-        edge is a GCA or a constant latitudinal/longitudinal line is based on `is_latlonface`.
-        Default is `None`.
-
-    Returns
-    -------
-    face_latlon_array : np.ndarray, shape (2, 2)
-        An array representing the bounding box for the face in latitude and longitude
-        coordinates (in radians). The first row contains the minimum and maximum latitude
-        values, while the second row contains the minimum and maximum longitude values.
-
-    Notes
-    -----
-    This function evaluates the presence of North or South pole points within the face's
-    bounds by inspecting the Cartesian coordinates of the face's edges. It then constructs
-    the face's bounding box by considering the extreme latitude and longitude values found
-    among the face's edges, adjusting for the presence of pole points as necessary.
-
-    The bounding box is used to determine the face's geographical extent and is crucial
-    for spatial analyses involving the grid.
-
-    Example
-    -------
-    Assuming the existence of a grid face with edges defined in both Cartesian and
-    latitudinal/longitudinal coordinates:
-
-        face_edges_cartesian = np.array([...])  # Cartesian coords
-        face_edges_connectivity_rad = np.array([...])  # Lon/Lat coords in radians
-
-    Populate the bounding box for the face, treating it as a latlon face:
-
-        face_latlon_bound = _populate_face_latlon_bound(face_edges_cartesian,
-                                                        face_edges_lonlat_connectivity_rad,
-                                                        is_latlonface=True)
-    """
-
     # Check if face_edges contains pole points
-    has_north_pole = _pole_point_inside_polygon("North", face_edges_cartesian)
-    has_south_pole = _pole_point_inside_polygon("South", face_edges_cartesian)
+    has_north_pole = pole_point_inside_polygon(1, face_edges_xyz, face_edges_lonlat)
+    has_south_pole = pole_point_inside_polygon(-1, face_edges_xyz, face_edges_lonlat)
 
-    face_latlon_array = np.full((2, 2), INT_FILL_VALUE, dtype=np.float64)
+    # Initialize face_latlon_array with INT_FILL_VALUE
+    face_latlon_array = np.full((2, 2), np.nan, dtype=np.float64)
 
     if has_north_pole or has_south_pole:
         # Initial assumption that the pole point is inside the face
         is_center_pole = True
 
-        # Define the pole point based on the hemisphere
-        pole_point = (
-            np.array([0.0, 0.0, 1.0]) if has_north_pole else np.array([0.0, 0.0, -1.0])
-        )
-        # Pre-defined new point latitude based on the pole
-        new_pt_latlon = np.array(
-            [np.pi / 2 if has_north_pole else -np.pi / 2, INT_FILL_VALUE],
-            dtype=np.float64,
-        )
+        pole_point_xyz = np.zeros(3, dtype=np.float64)
+        pole_point_lonlat = np.zeros(2, dtype=np.float64)
+        new_pt_latlon = np.zeros(2, dtype=np.float64)
 
-        for i in range(face_edges_cartesian.shape[0]):
-            edge_cart = face_edges_cartesian[i]
-            edge_lonlat = face_edges_lonlat_rad[i]
+        if has_north_pole:
+            pole_point_xyz[2] = 1.0  # [0.0, 0.0, 1.0]
+            pole_point_lonlat[1] = math.pi / 2  # [0.0, pi/2]
+            new_pt_latlon[0] = math.pi / 2  # [pi/2, INT_FILL_VALUE]
+            new_pt_latlon[1] = np.nan
+        else:
+            pole_point_xyz[2] = -1.0  # [0.0, 0.0, -1.0]
+            pole_point_lonlat[1] = -math.pi / 2  # [0.0, -pi/2]
+            new_pt_latlon[0] = -math.pi / 2  # [-pi/2, INT_FILL_VALUE]
+            new_pt_latlon[1] = np.nan
 
-            # Skip processing if the edge_cart is marked as a dummy with a fill value
-            if np.any(edge_cart == INT_FILL_VALUE):
-                continue
+        for i in range(face_edges_xyz.shape[0]):
+            edge_xyz = face_edges_xyz[i]
+            edge_lonlat = face_edges_lonlat[i]
 
-            # Extract cartesian coordinates of the edge_cart's endpoints
-            n1_cart, n2_cart = edge_cart
-            n1_lonlat, n2_lonlat = edge_lonlat
+            # Skip processing if the edge is marked as a dummy with a fill value
+            if np.any(edge_xyz == INT_FILL_VALUE):
+                continue
 
-            # Convert latitudes and longitudes of the nodes to radians
-            node1_lon_rad, node1_lat_rad = n1_lonlat
+            # Extract cartesian coordinates of the edge's endpoints
+            n1_cart = edge_xyz[0]
+            n2_cart = edge_xyz[1]
+            n1_lonlat = edge_lonlat[0]
+            n2_lonlat = edge_lonlat[1]
 
-            # Determine if the edge_cart's extreme latitudes need to be considered using the corrected logic
-            is_GCA = (
-                is_GCA_list[i]
-                if is_GCA_list is not None
-                else not is_latlonface or n1_cart[2] != n2_cart[2]
-            )
+            # Extract latitudes and longitudes of the nodes
+            node1_lon_rad = n1_lonlat[0]
+            node1_lat_rad = n1_lonlat[1]
 
-            # Check if the node matches the pole point or if the pole point is within the edge_cart
-            if allclose(n1_cart, pole_point, atol=ERROR_TOLERANCE) or point_within_gca(
-                pole_point, np.array([n1_cart, n2_cart]), is_directed=False
+            # Determine if the edge's extreme latitudes need to be considered
+            if is_GCA_list is not None:
+                is_GCA = is_GCA_list[i]
+            else:
+                is_GCA = not is_latlonface or n1_cart[2] != n2_cart[2]
+
+            # Check if the node matches the pole point or if the pole point is within the edge
+            max_abs_diff = np.max(np.abs(n1_cart - pole_point_xyz))
+            if max_abs_diff <= ERROR_TOLERANCE or point_within_gca(
+                pole_point_xyz,
+                pole_point_lonlat,
+                n1_cart,
+                n1_lonlat,
+                n2_cart,
+                n2_lonlat,
+                is_directed=False,
             ):
                 is_center_pole = False
-                face_latlon_array = _insert_pt_in_latlonbox(
+                face_latlon_array = insert_pt_in_latlonbox(
                     face_latlon_array, new_pt_latlon
                 )
 
             # Insert the current node's lat/lon into the latlonbox
-            face_latlon_array = _insert_pt_in_latlonbox(
+            face_latlon_array = insert_pt_in_latlonbox(
                 face_latlon_array, np.array([node1_lat_rad, node1_lon_rad])
             )
 
+            # Create n1n2_cart and n1n2_lonlat arrays using preallocation
+            n1n2_cart = np.empty((2, 3), dtype=np.float64)
+            n1n2_cart[0, :] = n1_cart
+            n1n2_cart[1, :] = n2_cart
+
+            n1n2_lonlat = np.empty((2, 2), dtype=np.float64)
+            n1n2_lonlat[0, :] = n1_lonlat
+            n1n2_lonlat[1, :] = n2_lonlat
+
             # Determine extreme latitudes for GCA edges
-            lat_max, lat_min = (
-                (
-                    extreme_gca_latitude(np.array([n1_cart, n2_cart]), "max"),
-                    extreme_gca_latitude(np.array([n1_cart, n2_cart]), "min"),
-                )
-                if is_GCA
-                else (node1_lat_rad, node1_lat_rad)
-            )
+            if is_GCA:
+                lat_max = extreme_gca_latitude(n1n2_cart, n1n2_lonlat, "max")
+                lat_min = extreme_gca_latitude(n1n2_cart, n1n2_lonlat, "min")
+            else:
+                lat_max = node1_lat_rad
+                lat_min = node1_lat_rad
 
             # Insert latitudinal extremes based on pole presence
             if has_north_pole:
-                face_latlon_array = _insert_pt_in_latlonbox(
+                face_latlon_array = insert_pt_in_latlonbox(
                     face_latlon_array, np.array([lat_min, node1_lon_rad])
                 )
-                face_latlon_array[0][1] = (
-                    np.pi / 2
-                )  # Ensure north pole is the upper latitude bound
+                face_latlon_array[0, 1] = math.pi / 2  # Upper latitude bound
             else:
-                face_latlon_array = _insert_pt_in_latlonbox(
+                face_latlon_array = insert_pt_in_latlonbox(
                     face_latlon_array, np.array([lat_max, node1_lon_rad])
                 )
-                face_latlon_array[0][0] = (
-                    -np.pi / 2
-                )  # Ensure south pole is the lower latitude bound
+                face_latlon_array[0, 0] = -math.pi / 2  # Lower latitude bound
 
         # Adjust longitude bounds globally if the pole is centrally inside the polygon
         if is_center_pole:
-            face_latlon_array[1] = [0.0, 2 * np.pi]
+            face_latlon_array[1, 0] = 0.0
+            face_latlon_array[1, 1] = 2 * math.pi
 
     else:
         # Normal Face
-        # Iterate through each edge_cart of a face to update the bounding box (latlonbox) with extreme latitudes and longitudes
-        for i in range(face_edges_cartesian.shape[0]):
-            edge_cart = face_edges_cartesian[i]
-            edge_lonlat = face_edges_lonlat_rad[i]
+        for i in range(face_edges_xyz.shape[0]):
+            edge_xyz = face_edges_xyz[i]
+            edge_lonlat = face_edges_lonlat[i]
 
-            # Skip processing if the edge_cart is marked as a dummy with a fill value
-            if np.any(edge_cart == INT_FILL_VALUE):
+            # Skip processing if the edge is marked as a dummy with a fill value
+            if np.any(edge_xyz == INT_FILL_VALUE):
                 continue
 
-            # Extract cartesian coordinates of the edge_cart's endpoints
-            n1_cart, n2_cart = edge_cart
-            n1_lonlat, n2_lonlat = edge_lonlat
+            # Extract cartesian coordinates of the edge's endpoints
+            n1_cart = edge_xyz[0]
+            n2_cart = edge_xyz[1]
+            n1_lonlat = edge_lonlat[0]
+            n2_lonlat = edge_lonlat[1]
+
+            # Extract latitudes and longitudes of the nodes
+            node1_lon_rad = n1_lonlat[0]
+            node1_lat_rad = n1_lonlat[1]
+            # node2_lon_rad = n2_lonlat[0]
+            node2_lat_rad = n2_lonlat[1]
+
+            # Determine if the edge's extreme latitudes need to be considered
+            if is_GCA_list is not None:
+                is_GCA = is_GCA_list[i]
+            else:
+                is_GCA = not is_latlonface or n1_cart[2] != n2_cart[2]
 
-            # Convert latitudes and longitudes of the nodes to radians
-            node1_lon_rad, node1_lat_rad = n1_lonlat
-            node2_lon_rad, node2_lat_rad = n2_lonlat
+            # Create n1n2_cart and n1n2_lonlat arrays using preallocation
+            n1n2_cart = np.empty((2, 3), dtype=np.float64)
+            n1n2_cart[0, :] = n1_cart
+            n1n2_cart[1, :] = n2_cart
 
-            # Determine if the edge_cart's extreme latitudes need to be considered using the corrected logic
-            is_GCA = (
-                is_GCA_list[i]
-                if is_GCA_list is not None
-                else not is_latlonface or n1_cart[2] != n2_cart[2]
-            )
+            n1n2_lonlat = np.empty((2, 2), dtype=np.float64)
+            n1n2_lonlat[0, :] = n1_lonlat
+            n1n2_lonlat[1, :] = n2_lonlat
 
-            lat_max, lat_min = (
-                (
-                    extreme_gca_latitude(np.array([n1_cart, n2_cart]), "max"),
-                    extreme_gca_latitude(np.array([n1_cart, n2_cart]), "min"),
-                )
-                if is_GCA
-                else (node1_lat_rad, node1_lat_rad)
-            )
+            if is_GCA:
+                lat_max = extreme_gca_latitude(n1n2_cart, n1n2_lonlat, "max")
+                lat_min = extreme_gca_latitude(n1n2_cart, n1n2_lonlat, "min")
+            else:
+                lat_max = node1_lat_rad
+                lat_min = node1_lat_rad
 
-            # Insert extreme latitude points into the latlonbox if they differ from the node latitudes
-            if not isclose(
-                node1_lat_rad, lat_max, atol=ERROR_TOLERANCE
-            ) and not isclose(node2_lat_rad, lat_max, atol=ERROR_TOLERANCE):
-                # Insert the maximum latitude
-                face_latlon_array = _insert_pt_in_latlonbox(
+            # Insert extreme latitude points into the latlonbox
+            if (
+                abs(node1_lat_rad - lat_max) > ERROR_TOLERANCE
+                and abs(node2_lat_rad - lat_max) > ERROR_TOLERANCE
+            ):
+                face_latlon_array = insert_pt_in_latlonbox(
                     face_latlon_array, np.array([lat_max, node1_lon_rad])
                 )
-            elif not isclose(
-                node1_lat_rad, lat_min, atol=ERROR_TOLERANCE
-            ) and not isclose(node2_lat_rad, lat_min, atol=ERROR_TOLERANCE):
-                # Insert the minimum latitude
-                face_latlon_array = _insert_pt_in_latlonbox(
+            elif (
+                abs(node1_lat_rad - lat_min) > ERROR_TOLERANCE
+                and abs(node2_lat_rad - lat_min) > ERROR_TOLERANCE
+            ):
+                face_latlon_array = insert_pt_in_latlonbox(
                     face_latlon_array, np.array([lat_min, node1_lon_rad])
                 )
             else:
-                # Insert the node's latitude and longitude as it matches the extreme latitudes
-                face_latlon_array = _insert_pt_in_latlonbox(
+                face_latlon_array = insert_pt_in_latlonbox(
                     face_latlon_array, np.array([node1_lat_rad, node1_lon_rad])
                 )
 
     return face_latlon_array
 
 
+@njit(cache=True, parallel=True)
+def compute_temp_latlon_array(
+    face_node_connectivity,
+    faces_edges_cartesian,
+    faces_edges_lonlat_rad,
+    n_nodes_per_face,
+    is_latlonface,
+    is_face_GCA_list,
+    INT_FILL_VALUE,
+):
+    n_face = face_node_connectivity.shape[0]
+    temp_latlon_array = np.full((n_face, 2, 2), INT_FILL_VALUE, dtype=np.float64)
+    for face_idx in prange(n_face):
+        cur_face_edges_cartesian = faces_edges_cartesian[
+            face_idx, 0 : n_nodes_per_face[face_idx]
+        ]
+        cur_faces_edges_lonlat_rad = faces_edges_lonlat_rad[
+            face_idx, 0 : n_nodes_per_face[face_idx]
+        ]
+        if is_face_GCA_list is not None:
+            is_GCA_list = is_face_GCA_list[face_idx]
+        else:
+            is_GCA_list = None
+
+        temp_latlon_array[face_idx] = _populate_face_latlon_bound(
+            cur_face_edges_cartesian,
+            cur_faces_edges_lonlat_rad,
+            is_latlonface=is_latlonface,
+            is_GCA_list=is_GCA_list,
+        )
+    return temp_latlon_array
+
+
 def _populate_bounds(
     grid, is_latlonface: bool = False, is_face_GCA_list=None, return_array=False
 ):
@@ -1245,13 +1409,11 @@ def _populate_bounds(
 
     This will calculate and store the bounds for each face within the grid, adjusting for any special conditions such as crossing the antimeridian, and return them as a DataArray.
     """
-    temp_latlon_array = np.full((grid.n_face, 2, 2), INT_FILL_VALUE, dtype=np.float64)
 
-    # Because Pandas.IntervalIndex does not support naming for each interval, we need to create a mapping
-    # between the intervals and the face indices
-    intervals_tuple_list = []
-    intervals_name_list = []
+    # Ensure grid's cartesian coordinates are normalized
+    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,
@@ -1269,42 +1431,29 @@ def _populate_bounds(
         grid.node_lat.values,
     )
 
-    for face_idx, face_nodes in enumerate(grid.face_node_connectivity):
-        face_edges_cartesian = faces_edges_cartesian[face_idx]
-
-        # Remove the edge in the face that contains the fill value
-        face_edges_cartesian = face_edges_cartesian[
-            np.all(face_edges_cartesian != INT_FILL_VALUE, axis=(1, 2))
-        ]
-
-        face_edges_lonlat_rad = faces_edges_lonlat_rad[face_idx]
-
-        # Remove the edge in the face that contains the fill value
-        face_edges_lonlat_rad = face_edges_lonlat_rad[
-            np.all(face_edges_lonlat_rad != INT_FILL_VALUE, axis=(1, 2))
-        ]
-
-        is_GCA_list = (
-            is_face_GCA_list[face_idx] if is_face_GCA_list is not None else None
-        )
+    n_nodes_per_face = grid.n_nodes_per_face.values
 
-        temp_latlon_array[face_idx] = _populate_face_latlon_bound(
-            face_edges_cartesian,
-            face_edges_lonlat_rad,
-            is_latlonface=is_latlonface,
-            is_GCA_list=is_GCA_list,
-        )
+    # Call the Numba-compiled function
+    temp_latlon_array = compute_temp_latlon_array(
+        grid.face_node_connectivity.values,
+        faces_edges_cartesian,
+        faces_edges_lonlat_rad,
+        n_nodes_per_face,
+        is_latlonface,
+        is_face_GCA_list,
+        INT_FILL_VALUE,
+    )
 
+    # Process results outside Numba
+    intervals_tuple_list = []
+    intervals_name_list = []
+    for face_idx in range(grid.n_face):
         assert temp_latlon_array[face_idx][0][0] != temp_latlon_array[face_idx][0][1]
         assert temp_latlon_array[face_idx][1][0] != temp_latlon_array[face_idx][1][1]
         lat_array = temp_latlon_array[face_idx][0]
-
-        # Now store the latitude intervals in the tuples
         intervals_tuple_list.append((lat_array[0], lat_array[1]))
         intervals_name_list.append(face_idx)
 
-    # Because Pandas.IntervalIndex does not support naming for each interval, we need to create a mapping
-    # between the intervals and the face indices
     intervalsIndex = pd.IntervalIndex.from_tuples(intervals_tuple_list, closed="both")
     df_intervals_map = pd.DataFrame(
         index=intervalsIndex, data=intervals_name_list, columns=["face_id"]
@@ -1312,7 +1461,7 @@ def _populate_bounds(
 
     bounds = xr.DataArray(
         temp_latlon_array,
-        dims=["n_face", "Two", "Two"],
+        dims=["n_face", "lon_lat", "min_max"],
         attrs={
             "cf_role": "face_latlon_bounds",
             "_FillValue": INT_FILL_VALUE,
diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py
index 88a755b0..8e0a9527 100644
--- a/uxarray/grid/grid.py
+++ b/uxarray/grid/grid.py
@@ -62,6 +62,7 @@ from uxarray.grid.geometry import (
     _grid_to_matplotlib_linecollection,
     _populate_bounds,
     _construct_boundary_edge_indices,
+    compute_temp_latlon_array,
 )
 
 from uxarray.grid.neighbors import (
@@ -75,6 +76,7 @@ from uxarray.grid.intersections import (
     constant_lat_intersections_no_extreme,
     constant_lon_intersections_no_extreme,
     constant_lat_intersections_face_bounds,
+    constant_lon_intersections_face_bounds,
 )
 
 from spatialpandas import GeoDataFrame
@@ -1360,15 +1362,19 @@ class Grid:
 
     @property
     def bounds(self):
-        """Latitude Longitude Bounds for each Face in degrees.
+        """Latitude Longitude Bounds for each Face in radians.
 
         Dimensions ``(n_face", two, two)``
         """
         if "bounds" not in self._ds:
-            warn(
-                "Computing 'Grid.bounds' for the first time. This may take some time...",
-                UserWarning,
-            )
+            if hasattr(compute_temp_latlon_array, "inspect_llvm"):
+                if len(compute_temp_latlon_array.inspect_llvm()) == 0:
+                    warn(
+                        "Necessary functions for computing face bounds are not translated yet with Numba. This initial"
+                        "translation may take some time.",
+                        RuntimeWarning,
+                    )
+
             _populate_bounds(self)
 
         return self._ds["bounds"]
@@ -1379,6 +1385,38 @@ class Grid:
         assert isinstance(value, xr.DataArray)
         self._ds["bounds"] = value
 
+    @property
+    def face_bounds_lon(self):
+        """Longitude bounds for each face in degrees."""
+
+        if "face_bounds_lon" not in self._ds:
+            bounds = self.bounds.values
+
+            bounds_deg = np.rad2deg(bounds[:, 1, :])
+            bounds_normalized = (bounds_deg + 180.0) % 360.0 - 180.0
+            bounds_lon = bounds_normalized
+            mask_zero = (bounds_lon[:, 0] == 0) & (bounds_lon[:, 1] == 0)
+            # for faces that span all longitudes (i.e. pole faces)
+            bounds_lon[mask_zero] = [-180.0, 180.0]
+            self._ds["face_bounds_lon"] = xr.DataArray(
+                data=bounds_lon,
+                dims=["n_face", "min_max"],
+            )
+
+        return self._ds["face_bounds_lon"]
+
+    @property
+    def face_bounds_lat(self):
+        """Latitude bounds for each face in degrees."""
+        if "face_bounds_lat" not in self._ds:
+            bounds = self.bounds.values
+            bounds_lat = np.sort(np.rad2deg(bounds[:, 0, :]), axis=-1)
+            self._ds["face_bounds_lat"] = xr.DataArray(
+                data=bounds_lat,
+                dims=["n_face", "min_max"],
+            )
+        return self._ds["face_bounds_lat"]
+
     @property
     def face_jacobian(self):
         """Declare face_jacobian as a property."""
@@ -2189,18 +2227,16 @@ class Grid:
                 "Indexing must be along a grid dimension: ('n_node', 'n_edge', 'n_face')"
             )
 
-    def get_edges_at_constant_latitude(self, lat, use_spherical_bounding_box=False):
+    def get_edges_at_constant_latitude(self, lat: float, use_face_bounds: bool = False):
         """Identifies the indices of edges that intersect with a line of constant latitude.
 
         Parameters
         ----------
-        lon : float
-            The latitude at which to identify intersecting edges, in degrees.
-                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`.
+        lat : float
+            The latitude at which to extract the cross-section, in degrees.
+            Must be between -90.0 and 90.0
+        use_face_bounds : bool, optional
+            If True, uses the bounds of each face for computing intersections.
 
         Returns
         -------
@@ -2213,7 +2249,7 @@ class Grid:
                 f"Latitude must be between -90 and 90 degrees. Received {lat}"
             )
 
-        if use_spherical_bounding_box:
+        if use_face_bounds:
             raise NotImplementedError(
                 "Computing the intersection using the spherical bounding box"
                 "is not yet supported."
@@ -2225,23 +2261,18 @@ class Grid:
 
         return edges.squeeze()
 
-    def get_faces_at_constant_latitude(self, lat, use_spherical_bounding_box=False):
+    def get_faces_at_constant_latitude(
+        self,
+        lat: float,
+    ):
         """
         Identifies the indices of faces that intersect with a line of constant latitude.
 
-        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.
-        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`.
+            The latitude at which to extract the cross-section, in degrees.
+            Must be between -90.0 and 90.0
 
         Returns
         -------
@@ -2254,32 +2285,25 @@ class Grid:
                 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())
-
-            return faces[faces != INT_FILL_VALUE]
+        faces = constant_lat_intersections_face_bounds(
+            lat=lat,
+            face_bounds_lat=self.face_bounds_lat.values,
+        )
+        return faces
 
-    def get_edges_at_constant_longitude(self, lon, use_spherical_bounding_box=False):
+    def get_edges_at_constant_longitude(
+        self, lon: float, use_face_bounds: bool = 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.
-        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`.
+            The longitude at which to extract the cross-section, in degrees.
+            Must be between -90.0 and 90.0
+        use_face_bounds : bool, optional
+            If True, uses the bounds of each face for computing intersections.
 
         Returns
         -------
@@ -2292,7 +2316,7 @@ class Grid:
                 f"Longitude must be between -180 and 180 degrees. Received {lon}"
             )
 
-        if use_spherical_bounding_box:
+        if use_face_bounds:
             raise NotImplementedError(
                 "Computing the intersection using the spherical bounding box"
                 "is not yet supported."
@@ -2303,23 +2327,15 @@ class Grid:
             )
             return edges.squeeze()
 
-    def get_faces_at_constant_longitude(self, lon, use_spherical_bounding_box=False):
+    def get_faces_at_constant_longitude(self, lon: float):
         """
         Identifies the indices of faces that intersect with a line of constant longitude.
 
-        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.
-        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`.
+            The longitude at which to extract the cross-section, in degrees.
+            Must be between -90.0 and 90.0
 
         Returns
         -------
@@ -2327,15 +2343,10 @@ class Grid:
             An array of face indices that intersect with the specified longitude.
         """
 
-        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
+        if lon > 180.0 or lon < -180.0:
+            raise ValueError(
+                f"Longitude must be between -180 and 180 degrees. Received {lon}"
             )
-            faces = np.unique(self.edge_face_connectivity[edges].data.ravel())
 
-            return faces[faces != INT_FILL_VALUE]
+        faces = constant_lon_intersections_face_bounds(lon, self.face_bounds_lon.values)
+        return faces
diff --git a/uxarray/grid/intersections.py b/uxarray/grid/intersections.py
index 147a4b15..528304d3 100644
--- a/uxarray/grid/intersections.py
+++ b/uxarray/grid/intersections.py
@@ -1,10 +1,16 @@
 import numpy as np
-from uxarray.constants import MACHINE_EPSILON, ERROR_TOLERANCE
+from uxarray.constants import MACHINE_EPSILON, ERROR_TOLERANCE, INT_DTYPE
 from uxarray.grid.utils import _newton_raphson_solver_for_gca_constLat
-from uxarray.grid.arcs import point_within_gca, extreme_gca_latitude, in_between
+from uxarray.grid.arcs import (
+    point_within_gca,
+    in_between,
+    _extreme_gca_latitude_cartesian,
+    _point_within_gca_cartesian,
+)
+from uxarray.grid.coordinates import _xyz_to_lonlat_rad_scalar
 import platform
 import warnings
-from uxarray.utils.computing import cross_fma, allclose, dot, cross, norm
+from uxarray.utils.computing import cross_fma, allclose, cross, norm
 
 
 from numba import njit, prange
@@ -122,8 +128,8 @@ def constant_lon_intersections_no_extreme(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):
+@njit(cache=True)
+def constant_lat_intersections_face_bounds(lat, face_bounds_lat):
     """Identifies the candidate faces on a grid that intersect with a given
     constant latitude.
 
@@ -136,24 +142,24 @@ def constant_lat_intersections_face_bounds(lat, face_min_lat_rad, face_max_lat_r
     ----------
     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.
+    TODO:
 
     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)
+
+    face_bounds_lat_min = face_bounds_lat[:, 0]
+    face_bounds_lat_max = face_bounds_lat[:, 1]
+
+    within_bounds = (face_bounds_lat_min <= lat) & (face_bounds_lat_max >= 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):
+def constant_lon_intersections_face_bounds(lon, face_bounds_lon):
     """Identifies the candidate faces on a grid that intersect with a given
     constant longitude.
 
@@ -166,10 +172,7 @@ def constant_lon_intersections_face_bounds(lon, face_min_lon_rad, face_max_lon_r
     ----------
     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.
+    TODO:
 
     Returns
     -------
@@ -177,121 +180,113 @@ def constant_lon_intersections_face_bounds(lon, face_min_lon_rad, face_max_lon_r
         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.
-
-    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
-    ----------
-    gca1_cart : [n, 3] np.ndarray where n is the number of intersection points
-        Cartesian coordinates of the first GCA.
-    gca2_cart : [n, 3] np.ndarray where n is the number of intersection points
-        Cartesian coordinates of the second GCA.
-    fma_disabled : bool, optional (default=True)
-        If True, the FMA operation is disabled. And a naive `np.cross` is used instead.
-
-    Returns
-    -------
-    np.ndarray
-        Cartesian coordinates of the intersection point(s). Returns an empty array if no intersections,
-        a 2D array with one row if one intersection, and a 2D array with two rows if two intersections.
+    face_bounds_lon_min = face_bounds_lon[:, 0]
+    face_bounds_lon_max = face_bounds_lon[:, 1]
+    n_face = face_bounds_lon.shape[0]
+
+    candidate_faces = []
+    for i in range(n_face):
+        cur_face_bounds_lon_min = face_bounds_lon_min[i]
+        cur_face_bounds_lon_max = face_bounds_lon_max[i]
+
+        if cur_face_bounds_lon_min < cur_face_bounds_lon_max:
+            if (lon >= cur_face_bounds_lon_min) & (lon <= cur_face_bounds_lon_max):
+                candidate_faces.append(i)
+        else:
+            # antimeridian case
+            if (lon >= cur_face_bounds_lon_min) | (lon <= cur_face_bounds_lon_max):
+                candidate_faces.append(i)
+
+    return np.array(candidate_faces, dtype=INT_DTYPE)
+
+
+def _gca_gca_intersection_cartesian(gca_a_xyz, gca_b_xyz):
+    gca_a_xyz = np.asarray(gca_a_xyz)
+    gca_b_xyz = np.asarray(gca_b_xyz)
+
+    gca_a_lonlat = np.array(
+        [
+            _xyz_to_lonlat_rad_scalar(
+                gca_a_xyz[0, 0], gca_a_xyz[0, 1], gca_a_xyz[0, 2]
+            ),
+            _xyz_to_lonlat_rad_scalar(
+                gca_a_xyz[1, 0], gca_a_xyz[1, 1], gca_a_xyz[1, 2]
+            ),
+        ]
+    )
+    gca_b_lonlat = np.array(
+        [
+            _xyz_to_lonlat_rad_scalar(
+                gca_b_xyz[0, 0], gca_b_xyz[0, 1], gca_b_xyz[0, 2]
+            ),
+            _xyz_to_lonlat_rad_scalar(
+                gca_b_xyz[1, 0], gca_b_xyz[1, 1], gca_b_xyz[1, 2]
+            ),
+        ]
+    )
+    return gca_gca_intersection(gca_a_xyz, gca_a_lonlat, gca_b_xyz, gca_b_lonlat)
 
-    Raises
-    ------
-    ValueError
-        If the input GCAs are not in the cartesian [x, y, z] format.
-    """
 
-    # Support lists as an input
-    gca1_cart = np.asarray(gca1_cart)
-    gca2_cart = np.asarray(gca2_cart)
-    # Check if the two GCAs are in the cartesian format (size of three)
-    if gca1_cart.shape[1] != 3 or gca2_cart.shape[1] != 3:
+@njit(cache=True)
+def gca_gca_intersection(gca_a_xyz, gca_a_lonlat, gca_b_xyz, gca_b_lonlat):
+    if gca_a_xyz.shape[1] != 3 or gca_b_xyz.shape[1] != 3:
         raise ValueError("The two GCAs must be in the cartesian [x, y, z] format")
 
-    w0, w1 = gca1_cart
-    v0, v1 = gca2_cart
+    # Extract points
+    w0_xyz = gca_a_xyz[0]
+    w1_xyz = gca_a_xyz[1]
+    v0_xyz = gca_b_xyz[0]
+    v1_xyz = gca_b_xyz[1]
 
-    # Compute normals and orthogonal bases using FMA
-    if fma_disabled:
-        w0w1_norm = cross(w0, w1)
-        v0v1_norm = cross(v0, v1)
-        cross_norms = cross(w0w1_norm, v0v1_norm)
-    else:
-        w0w1_norm = cross_fma(w0, w1)
-        v0v1_norm = cross_fma(v0, v1)
-        cross_norms = cross_fma(w0w1_norm, v0v1_norm)
+    w0_lonlat = gca_a_lonlat[0]
+    w1_lonlat = gca_a_lonlat[1]
+    v0_lonlat = gca_b_lonlat[0]
+    v1_lonlat = gca_b_lonlat[1]
 
-        # 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."
-            )
+    # Compute normals and orthogonal bases
+    w0w1_norm = cross(w0_xyz, w1_xyz)
+    v0v1_norm = cross(v0_xyz, v1_xyz)
+    cross_norms = cross(w0w1_norm, v0v1_norm)
 
-    # Check perpendicularity conditions and floating-point arithmetic limitations
-    if not allclose(dot(w0w1_norm, w0), 0.0, atol=MACHINE_EPSILON) or not allclose(
-        dot(w0w1_norm, w1), 0.0, atol=MACHINE_EPSILON
-    ):
-        warnings.warn(
-            "The current input data cannot be computed accurately using floating-point arithmetic. Use with caution."
-        )
+    # Initialize result array and counter
+    res = np.empty((2, 3))
+    count = 0
 
-    if not allclose(dot(v0v1_norm, v0), 0.0, 0.0, atol=MACHINE_EPSILON) or not allclose(
-        dot(v0v1_norm, v1), 0.0, atol=MACHINE_EPSILON
-    ):
-        warnings.warn(
-            "The current input data cannot be computed accurately using floating-point arithmetic.  Use with caution. "
-        )
-
-    if not allclose(
-        dot(cross_norms, v0v1_norm), 0.0, atol=MACHINE_EPSILON
-    ) or not allclose(dot(cross_norms, w0w1_norm), 0.0, atol=MACHINE_EPSILON):
-        warnings.warn(
-            "The current input data cannot be computed accurately using floating-point arithmetic. Use with caution. "
-        )
-
-    # If the cross_norms is zero, the two GCAs are parallel
+    # Check if the two GCAs are parallel
     if allclose(cross_norms, 0.0, atol=MACHINE_EPSILON):
-        res = []
-        # Check if the two GCAs are overlapping
-        if point_within_gca(v0, [w0, w1]):
-            # The two GCAs are overlapping, return both end points
-            res.append(v0)
+        if point_within_gca(v0_xyz, v0_lonlat, w0_xyz, w0_lonlat, w1_xyz, w1_lonlat):
+            res[count, :] = v0_xyz
+            count += 1
+
+        if point_within_gca(v1_xyz, v1_lonlat, w0_xyz, w0_lonlat, w1_xyz, w1_lonlat):
+            res[count, :] = v1_xyz
+            count += 1
 
-        if point_within_gca(v1, [w0, w1]):
-            res.append(v1)
-        return np.array(res)
+        return res[:count, :]
 
     # Normalize the cross_norms
     cross_norms = cross_norms / norm(cross_norms)
-    x1 = cross_norms
-    x2 = -x1
+    x1_xyz = cross_norms
+    x2_xyz = -x1_xyz
 
-    res = []
+    # Get lon/lat arrays
+    x1_lonlat = _xyz_to_lonlat_rad_scalar(x1_xyz[0], x1_xyz[1], x1_xyz[2])
+    x2_lonlat = _xyz_to_lonlat_rad_scalar(x2_xyz[0], x2_xyz[1], x2_xyz[2])
 
-    # Determine which intersection point is within the GCAs range
-    if point_within_gca(x1, [w0, w1]) and point_within_gca(x1, [v0, v1]):
-        res.append(x1)
+    # Check intersection points
+    if point_within_gca(
+        x1_xyz, x1_lonlat, w0_xyz, w0_lonlat, w1_xyz, w1_lonlat
+    ) and point_within_gca(x1_xyz, x1_lonlat, v0_xyz, v0_lonlat, v1_xyz, v1_lonlat):
+        res[count, :] = x1_xyz
+        count += 1
 
-    if point_within_gca(x2, [w0, w1]) and point_within_gca(x2, [v0, v1]):
-        res.append(x2)
+    if point_within_gca(
+        x2_xyz, x2_lonlat, w0_xyz, w0_lonlat, w1_xyz, w1_lonlat
+    ) and point_within_gca(x2_xyz, x2_lonlat, v0_xyz, v0_lonlat, v1_xyz, v1_lonlat):
+        res[count, :] = x2_xyz
+        count += 1
 
-    return np.array(res)
+    return res[:count, :]
 
 
 def gca_const_lat_intersection(
@@ -351,8 +346,9 @@ def gca_const_lat_intersection(
         return res
 
     # If the constant latitude is not the same as the GCA endpoints, calculate the intersection point
-    lat_min = extreme_gca_latitude(gca_cart, extreme_type="min")
-    lat_max = extreme_gca_latitude(gca_cart, extreme_type="max")
+    # TODO:
+    lat_min = _extreme_gca_latitude_cartesian(gca_cart, extreme_type="min")
+    lat_max = _extreme_gca_latitude_cartesian(gca_cart, extreme_type="max")
 
     constLat_rad = np.arcsin(constZ)
 
@@ -389,7 +385,7 @@ def gca_const_lat_intersection(
     res = None
 
     # Now test which intersection point is within the GCA range
-    if point_within_gca(p1, gca_cart, is_directed=is_directed):
+    if _point_within_gca_cartesian(p1, gca_cart, is_directed=is_directed):
         try:
             converged_pt = _newton_raphson_solver_for_gca_constLat(
                 p1, gca_cart, verbose=verbose
@@ -411,7 +407,7 @@ def gca_const_lat_intersection(
         except RuntimeError:
             raise RuntimeError(f"Error encountered with initial guess: {p1}")
 
-    if point_within_gca(p2, gca_cart, is_directed=is_directed):
+    if _point_within_gca_cartesian(p2, gca_cart, is_directed=is_directed):
         try:
             converged_pt = _newton_raphson_solver_for_gca_constLat(
                 p2, gca_cart, verbose=verbose
diff --git a/uxarray/plot/accessor.py b/uxarray/plot/accessor.py
index 47cc5c9d..9c0adc9b 100644
--- a/uxarray/plot/accessor.py
+++ b/uxarray/plot/accessor.py
@@ -399,7 +399,6 @@ class UxDataArrayPlotAccessor:
 
         if projection is not None:
             kwargs["projection"] = projection
-            print("Hello")
             kwargs["geo"] = True
             if "crs" not in kwargs:
                 central_longitude = projection.proj4_params["lon_0"]
diff --git a/uxarray/utils/computing.py b/uxarray/utils/computing.py
index 588d921d..b3db2756 100644
--- a/uxarray/utils/computing.py
+++ b/uxarray/utils/computing.py
@@ -4,6 +4,16 @@ import sys
 from numba import njit
 
 
+@njit(cache=True)
+def clip(a, a_min, a_max):
+    return np.clip(a, a_min, a_max)
+
+
+@njit(cache=True)
+def arcsin(x):
+    return np.arcsin(x)
+
+
 @njit(cache=True)
 def all(a):
     """Numba decorated implementation of ``np.all()``


Test patch:
diff --git a/test/test_arcs.py b/test/test_arcs.py
index fa6bc849..c139e390 100644
--- a/test/test_arcs.py
+++ b/test/test_arcs.py
@@ -10,7 +10,7 @@ from pathlib import Path
 import uxarray as ux
 
 from uxarray.grid.coordinates import _lonlat_rad_to_xyz
-from uxarray.grid.arcs import point_within_gca
+from uxarray.grid.arcs import point_within_gca, _point_within_gca_cartesian
 
 try:
     import constants
@@ -38,7 +38,7 @@ class TestIntersectionPoint(TestCase):
 
         pt_same_lon_in = _lonlat_rad_to_xyz(0.0, 0.0)
         with self.assertRaises(ValueError):
-            point_within_gca(pt_same_lon_in, gcr_180degree_cart)
+            _point_within_gca_cartesian(pt_same_lon_in, gcr_180degree_cart)
 
         # Test when the point and the GCA all have the same longitude
         gcr_same_lon_cart = [
@@ -46,19 +46,19 @@ class TestIntersectionPoint(TestCase):
             _lonlat_rad_to_xyz(0.0, -1.5)
         ]
         pt_same_lon_in = _lonlat_rad_to_xyz(0.0, 0.0)
-        self.assertTrue(point_within_gca(pt_same_lon_in, gcr_same_lon_cart))
+        self.assertTrue(_point_within_gca_cartesian(pt_same_lon_in, gcr_same_lon_cart))
 
         pt_same_lon_out = _lonlat_rad_to_xyz(0.0, 1.500000000000001)
-        res = point_within_gca(pt_same_lon_out, gcr_same_lon_cart)
+        res = _point_within_gca_cartesian(pt_same_lon_out, gcr_same_lon_cart)
         self.assertFalse(res)
 
         pt_same_lon_out_2 = _lonlat_rad_to_xyz(0.1, 1.0)
-        res = point_within_gca(pt_same_lon_out_2, gcr_same_lon_cart)
+        res = _point_within_gca_cartesian(pt_same_lon_out_2, gcr_same_lon_cart)
         self.assertFalse(res)
 
         # And if we increase the digital place by one, it should be true again
         pt_same_lon_out_add_one_place = _lonlat_rad_to_xyz(0.0, 1.5000000000000001)
-        res = point_within_gca(pt_same_lon_out_add_one_place, gcr_same_lon_cart)
+        res = _point_within_gca_cartesian(pt_same_lon_out_add_one_place, gcr_same_lon_cart)
         self.assertTrue(res)
 
         # Normal case
@@ -69,7 +69,7 @@ class TestIntersectionPoint(TestCase):
                                                         -0.997]])
         pt_cart_within = np.array(
             [0.25616109352676675, 0.9246590335292105, -0.010021496695000144])
-        self.assertTrue(point_within_gca(pt_cart_within, gcr_cart_2, True))
+        self.assertTrue(_point_within_gca_cartesian(pt_cart_within, gcr_cart_2, True))
 
         # Test other more complicate cases : The anti-meridian case
 
@@ -80,16 +80,16 @@ class TestIntersectionPoint(TestCase):
         gcr_cart = np.array([[0.351, -0.724, 0.593], [0.617, 0.672, 0.410]])
         pt_cart = np.array(
             [0.9438777657502077, 0.1193199333436068, 0.922714737029319])
-        self.assertTrue(point_within_gca(pt_cart, gcr_cart, is_directed=True))
+        self.assertTrue(_point_within_gca_cartesian(pt_cart, gcr_cart, is_directed=True))
         # If we swap the gcr, it should throw a value error since it's larger than 180 degree
         gcr_cart_flip = np.array([[0.617, 0.672, 0.410], [0.351, -0.724,
                                                           0.593]])
         with self.assertRaises(ValueError):
-            point_within_gca(pt_cart, gcr_cart_flip, is_directed=True)
+            _point_within_gca_cartesian(pt_cart, gcr_cart_flip, is_directed=True)
 
         # If we flip the gcr in the undirected mode, it should still work
         self.assertTrue(
-            point_within_gca(pt_cart, gcr_cart_flip, is_directed=False))
+            _point_within_gca_cartesian(pt_cart, gcr_cart_flip, is_directed=False))
 
         # 2nd anti-meridian case
         # GCR vertex0 in radian : [4.104711496596806, 0.5352983676533828],
@@ -100,9 +100,9 @@ class TestIntersectionPoint(TestCase):
         pt_cart_within = np.array(
             [0.6136726305712109, 0.28442243941920053, -0.365605190899831])
         self.assertFalse(
-            point_within_gca(pt_cart_within, gcr_cart_1, is_directed=True))
+            _point_within_gca_cartesian(pt_cart_within, gcr_cart_1, is_directed=True))
         self.assertFalse(
-            point_within_gca(pt_cart_within, gcr_cart_1, is_directed=False))
+            _point_within_gca_cartesian(pt_cart_within, gcr_cart_1, is_directed=False))
 
         # The first case should not work and the second should work
         v1_rad = [0.1, 0.0]
@@ -112,10 +112,10 @@ class TestIntersectionPoint(TestCase):
         gcr_cart = np.array([v1_cart, v2_cart])
         pt_cart = _lonlat_rad_to_xyz(0.01, 0.0)
         with self.assertRaises(ValueError):
-            point_within_gca(pt_cart, gcr_cart, is_directed=True)
+            _point_within_gca_cartesian(pt_cart, gcr_cart, is_directed=True)
         gcr_car_flipped = np.array([v2_cart, v1_cart])
         self.assertTrue(
-            point_within_gca(pt_cart, gcr_car_flipped, is_directed=True))
+            _point_within_gca_cartesian(pt_cart, gcr_car_flipped, is_directed=True))
 
     def test_pt_within_gcr_cross_pole(self):
         gcr_cart = np.array([[0.351, 0.0, 0.3], [-0.351, 0.0, 0.3]])
@@ -125,7 +125,7 @@ class TestIntersectionPoint(TestCase):
         # Normalize the point abd the GCA
         pt_cart = pt_cart / np.linalg.norm(pt_cart)
         gcr_cart = np.array([x / np.linalg.norm(x) for x in gcr_cart])
-        self.assertTrue(point_within_gca(pt_cart, gcr_cart, is_directed=False))
+        self.assertTrue(_point_within_gca_cartesian(pt_cart, gcr_cart, is_directed=False))
 
         gcr_cart = np.array([[0.351, 0.0, 0.3], [-0.351, 0.0, -0.6]])
         pt_cart = np.array(
@@ -134,6 +134,6 @@ class TestIntersectionPoint(TestCase):
         # When the point is not within the GCA
         pt_cart = pt_cart / np.linalg.norm(pt_cart)
         gcr_cart = np.array([x / np.linalg.norm(x) for x in gcr_cart])
-        self.assertFalse(point_within_gca(pt_cart, gcr_cart, is_directed=False))
+        self.assertFalse(_point_within_gca_cartesian(pt_cart, gcr_cart, is_directed=False))
         with self.assertRaises(ValueError):
-            point_within_gca(pt_cart, gcr_cart, is_directed=True)
+            _point_within_gca_cartesian(pt_cart, gcr_cart, is_directed=True)
diff --git a/test/test_cross_sections.py b/test/test_cross_sections.py
index 1847e69b..6fbe8abc 100644
--- a/test/test_cross_sections.py
+++ b/test/test_cross_sections.py
@@ -35,38 +35,34 @@ class TestQuadHex:
     All four faces intersect a constant latitude of 0.0
     """
 
-    @pytest.mark.parametrize("use_spherical_bounding_box", [True, False])
-    def test_constant_lat_cross_section_grid(self, use_spherical_bounding_box):
-
-
+    def test_constant_lat_cross_section_grid(self):
 
         uxgrid = ux.open_grid(quad_hex_grid_path)
 
-        grid_top_two = uxgrid.cross_section.constant_latitude(lat=0.1, use_spherical_bounding_box=use_spherical_bounding_box)
+        grid_top_two = uxgrid.cross_section.constant_latitude(lat=0.1, )
 
         assert grid_top_two.n_face == 2
 
-        grid_bottom_two = uxgrid.cross_section.constant_latitude(lat=-0.1, use_spherical_bounding_box=use_spherical_bounding_box)
+        grid_bottom_two = uxgrid.cross_section.constant_latitude(lat=-0.1, )
 
         assert grid_bottom_two.n_face == 2
 
-        grid_all_four = uxgrid.cross_section.constant_latitude(lat=0.0, use_spherical_bounding_box=use_spherical_bounding_box)
+        grid_all_four = uxgrid.cross_section.constant_latitude(lat=0.0, )
 
         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, use_spherical_bounding_box=use_spherical_bounding_box)
+            uxgrid.cross_section.constant_latitude(lat=10.0, )
 
-    @pytest.mark.parametrize("use_spherical_bounding_box", [False])
-    def test_constant_lon_cross_section_grid(self, use_spherical_bounding_box):
+    def test_constant_lon_cross_section_grid(self):
         uxgrid = ux.open_grid(quad_hex_grid_path)
 
-        grid_left_two = uxgrid.cross_section.constant_longitude(lon=-0.1, use_spherical_bounding_box=use_spherical_bounding_box)
+        grid_left_two = uxgrid.cross_section.constant_longitude(lon=-0.1, )
 
         assert grid_left_two.n_face == 2
 
-        grid_right_two = uxgrid.cross_section.constant_longitude(lon=0.2, use_spherical_bounding_box=use_spherical_bounding_box)
+        grid_right_two = uxgrid.cross_section.constant_longitude(lon=0.2, )
 
         assert grid_right_two.n_face == 2
 
@@ -74,64 +70,63 @@ class TestQuadHex:
             # no intersections found at this line
             uxgrid.cross_section.constant_longitude(lon=10.0)
 
-    @pytest.mark.parametrize("use_spherical_bounding_box", [False])
-    def test_constant_lat_cross_section_uxds(self, use_spherical_bounding_box):
+    def test_constant_lat_cross_section_uxds(self):
         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, use_spherical_bounding_box=use_spherical_bounding_box)
+        da_top_two = uxds['t2m'].cross_section.constant_latitude(lat=0.1, )
 
         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, use_spherical_bounding_box=use_spherical_bounding_box)
+        da_bottom_two = uxds['t2m'].cross_section.constant_latitude(lat=-0.1, )
 
         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, use_spherical_bounding_box=use_spherical_bounding_box)
+        da_all_four = uxds['t2m'].cross_section.constant_latitude(lat=0.0, )
 
         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, use_spherical_bounding_box=use_spherical_bounding_box)
+            uxds['t2m'].cross_section.constant_latitude(lat=10.0, )
+
 
-    @pytest.mark.parametrize("use_spherical_bounding_box", [False])
-    def test_constant_lon_cross_section_uxds(self, use_spherical_bounding_box):
+    def test_constant_lon_cross_section_uxds(self):
         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, use_spherical_bounding_box=use_spherical_bounding_box)
+        da_left_two = uxds['t2m'].cross_section.constant_longitude(lon=-0.1, )
 
         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, use_spherical_bounding_box=use_spherical_bounding_box)
+        da_right_two = uxds['t2m'].cross_section.constant_longitude(lon=0.2, )
 
         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, use_spherical_bounding_box=use_spherical_bounding_box)
+            uxds['t2m'].cross_section.constant_longitude(lon=10.0, )
 
 
 class TestCubeSphere:
-    @pytest.mark.parametrize("use_spherical_bounding_box", [True, False])
-    def test_north_pole(self, use_spherical_bounding_box):
+
+    def test_north_pole(self):
         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, use_spherical_bounding_box=use_spherical_bounding_box)
+            cross_grid = uxgrid.cross_section.constant_latitude(lat=lat, )
             # Cube sphere grid should have 4 faces centered around the pole
             assert cross_grid.n_face == 4
-    @pytest.mark.parametrize("use_spherical_bounding_box", [True, False])
-    def test_south_pole(self, use_spherical_bounding_box):
+
+    def test_south_pole(self):
         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, use_spherical_bounding_box=use_spherical_bounding_box)
+            cross_grid = uxgrid.cross_section.constant_latitude(lat=lat, )
             # Cube sphere grid should have 4 faces centered around the pole
             assert cross_grid.n_face == 4
 
@@ -152,8 +147,7 @@ class TestCandidateFacesUsingBounds:
 
         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],
+            face_bounds_lat=bounds_rad[:, 0],
         )
 
         # Expected output
@@ -176,8 +170,7 @@ class TestCandidateFacesUsingBounds:
 
         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],
+            face_bounds_lat=bounds_rad[:, 0],
         )
 
         assert len(candidate_faces) == 0
diff --git a/test/test_geometry.py b/test/test_geometry.py
index 67e9369b..05e0930f 100644
--- a/test/test_geometry.py
+++ b/test/test_geometry.py
@@ -10,12 +10,11 @@ import uxarray as ux
 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
+from uxarray.grid.arcs import extreme_gca_latitude, _extreme_gca_latitude_cartesian
 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, stereographic_projection, \
-    inverse_stereographic_projection
+from uxarray.grid.geometry import _populate_face_latlon_bound, _populate_bounds, _pole_point_inside_polygon_cartesian, stereographic_projection, inverse_stereographic_projection
+
 
-from spatialpandas.geometry import MultiPolygon
 
 current_path = Path(os.path.dirname(os.path.realpath(__file__)))
 
@@ -84,12 +83,12 @@ class TestPredicate(TestCase):
                                    [vertices[3], vertices[0]]])
 
         # Check if the North pole is inside the polygon
-        result = ux.grid.geometry._pole_point_inside_polygon(
+        result = _pole_point_inside_polygon_cartesian(
             'North', face_edge_cart)
         self.assertTrue(result, "North pole should be inside the polygon")
 
         # Check if the South pole is inside the polygon
-        result = ux.grid.geometry._pole_point_inside_polygon(
+        result = _pole_point_inside_polygon_cartesian(
             'South', face_edge_cart)
         self.assertFalse(result, "South pole should not be inside the polygon")
 
@@ -111,12 +110,12 @@ class TestPredicate(TestCase):
                                    [vertices[2], vertices[0]]])
 
         # Check if the North pole is inside the polygon
-        result = ux.grid.geometry._pole_point_inside_polygon(
+        result = _pole_point_inside_polygon_cartesian(
             'North', face_edge_cart)
         self.assertFalse(result, "North pole should not be inside the polygon")
 
         # Check if the South pole is inside the polygon
-        result = ux.grid.geometry._pole_point_inside_polygon(
+        result = _pole_point_inside_polygon_cartesian(
             'South', face_edge_cart)
         self.assertTrue(result, "South pole should be inside the polygon")
 
@@ -138,12 +137,12 @@ class TestPredicate(TestCase):
                                    [vertices[3], vertices[0]]])
 
         # Check if the North pole is inside the polygon
-        result = ux.grid.geometry._pole_point_inside_polygon(
+        result = _pole_point_inside_polygon_cartesian(
             'North', face_edge_cart)
         self.assertTrue(result, "North pole should be inside the polygon")
 
         # Check if the South pole is inside the polygon
-        result = ux.grid.geometry._pole_point_inside_polygon(
+        result = _pole_point_inside_polygon_cartesian(
             'South', face_edge_cart)
         self.assertFalse(result, "South pole should not be inside the polygon")
 
@@ -164,7 +163,7 @@ class TestPredicate(TestCase):
                                    [vertices[3], vertices[0]]])
 
         # Check if the North pole is inside the polygon
-        result = ux.grid.geometry._pole_point_inside_polygon(
+        result = _pole_point_inside_polygon_cartesian(
             'North', face_edge_cart)
         self.assertTrue(result, "North pole should be inside the polygon")
 
@@ -354,7 +353,7 @@ class TestLatlonBoundUtils(TestCase):
         ])
 
         # Calculate the maximum latitude
-        max_latitude = ux.grid.arcs.extreme_gca_latitude(gca_cart, 'max')
+        max_latitude = _extreme_gca_latitude_cartesian(gca_cart, 'max')
 
         # Check if the maximum latitude is correct
         expected_max_latitude = self._max_latitude_rad_iterative(gca_cart)
@@ -366,7 +365,7 @@ class TestLatlonBoundUtils(TestCase):
         gca_cart = np.array([[0.0, 0.0, 1.0], [1.0, 0.0, 0.0]])
 
         # Calculate the maximum latitude
-        max_latitude = ux.grid.arcs.extreme_gca_latitude(gca_cart, 'max')
+        max_latitude = _extreme_gca_latitude_cartesian(gca_cart, 'max')
 
         # Check if the maximum latitude is correct
         expected_max_latitude = np.pi / 2  # 90 degrees in radians
@@ -379,7 +378,7 @@ class TestLatlonBoundUtils(TestCase):
         gca_cart = np.array([[0.65465367, -0.37796447, -0.65465367], [0.6652466, -0.33896007, -0.6652466]])
 
         # Calculate the maximum latitude
-        max_latitude = ux.grid.arcs.extreme_gca_latitude(gca_cart, 'max')
+        max_latitude = _extreme_gca_latitude_cartesian(gca_cart, 'max')
 
         # Check if the maximum latitude is correct
         expected_max_latitude = self._max_latitude_rad_iterative(gca_cart)
@@ -395,7 +394,7 @@ class TestLatlonBoundUtils(TestCase):
         ])
 
         # Calculate the minimum latitude
-        min_latitude = ux.grid.arcs.extreme_gca_latitude(gca_cart, 'min')
+        min_latitude = _extreme_gca_latitude_cartesian(gca_cart, 'min')
 
         # Check if the minimum latitude is correct
         expected_min_latitude = self._min_latitude_rad_iterative(gca_cart)
@@ -407,7 +406,7 @@ class TestLatlonBoundUtils(TestCase):
         gca_cart = np.array([[0.0, 0.0, -1.0], [1.0, 0.0, 0.0]])
 
         # Calculate the minimum latitude
-        min_latitude = ux.grid.arcs.extreme_gca_latitude(gca_cart, 'min')
+        min_latitude = _extreme_gca_latitude_cartesian(gca_cart, 'min')
 
         # Check if the minimum latitude is correct
         expected_min_latitude = -np.pi / 2  # 90 degrees in radians
@@ -432,7 +431,7 @@ class TestLatlonBoundUtils(TestCase):
         old_box = np.array([[0.1, 0.2], [0.3, 0.4]])  # Radians
         new_pt = np.array([0.15, 0.35])
         expected = np.array([[0.1, 0.2], [0.3, 0.4]])
-        result = ux.grid.geometry._insert_pt_in_latlonbox(
+        result = ux.grid.geometry.insert_pt_in_latlonbox(
             old_box, new_pt, False)
         np.testing.assert_array_equal(result, expected)
 
@@ -440,22 +439,22 @@ class TestLatlonBoundUtils(TestCase):
         old_box = np.array([[0.1, 0.2], [6.0, 0.1]])  # Radians, periodic
         new_pt = np.array([0.15, 6.2])
         expected = np.array([[0.1, 0.2], [6.0, 0.1]])
-        result = ux.grid.geometry._insert_pt_in_latlonbox(old_box, new_pt, True)
+        result = ux.grid.geometry.insert_pt_in_latlonbox(old_box, new_pt, True)
         np.testing.assert_array_equal(result, expected)
 
     def test_insert_pt_in_latlonbox_pole(self):
         old_box = np.array([[0.1, 0.2], [0.3, 0.4]])
-        new_pt = np.array([np.pi / 2, INT_FILL_VALUE])  # Pole point
+        new_pt = np.array([np.pi / 2, np.nan])  # Pole point
         expected = np.array([[0.1, np.pi / 2], [0.3, 0.4]])
-        result = ux.grid.geometry._insert_pt_in_latlonbox(old_box, new_pt)
+        result = ux.grid.geometry.insert_pt_in_latlonbox(old_box, new_pt)
         np.testing.assert_array_equal(result, expected)
 
     def test_insert_pt_in_empty_state(self):
-        old_box = np.array([[INT_FILL_VALUE, INT_FILL_VALUE],
-                            [INT_FILL_VALUE, INT_FILL_VALUE]])  # Empty state
+        old_box = np.array([[np.nan, np.nan],
+                            [np.nan, np.nan]])  # Empty state
         new_pt = np.array([0.15, 0.35])
         expected = np.array([[0.15, 0.15], [0.35, 0.35]])
-        result = ux.grid.geometry._insert_pt_in_latlonbox(old_box, new_pt)
+        result = ux.grid.geometry.insert_pt_in_latlonbox(old_box, new_pt)
         np.testing.assert_array_equal(result, expected)
 
 
@@ -593,9 +592,9 @@ class TestLatlonBoundsGCA(TestCase):
         vertices_cart = np.vstack([_lonlat_rad_to_xyz(vertices_rad[:, 0], vertices_rad[:, 1])]).T
         # vertices_cart = [_lonlat_rad_to_xyz(vertices_rad[0], vertices_rad[1])]
         lat_max = max(np.deg2rad(60.0),
-                      extreme_gca_latitude(np.array([vertices_cart[0], vertices_cart[3]]), extreme_type="max"))
+                      _extreme_gca_latitude_cartesian(np.array([vertices_cart[0], vertices_cart[3]]), extreme_type="max"))
         lat_min = min(np.deg2rad(10.0),
-                      extreme_gca_latitude(np.array([vertices_cart[1], vertices_cart[2]]), extreme_type="min"))
+                      _extreme_gca_latitude_cartesian(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)
@@ -622,9 +621,9 @@ class TestLatlonBoundsGCA(TestCase):
         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(np.array([vertices_cart[0], vertices_cart[3]]), extreme_type="max"))
+                      _extreme_gca_latitude_cartesian(np.array([vertices_cart[0], vertices_cart[3]]), extreme_type="max"))
         lat_min = min(np.deg2rad(10.0),
-                      extreme_gca_latitude(np.array([vertices_cart[1], vertices_cart[2]]), extreme_type="min"))
+                      _extreme_gca_latitude_cartesian(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)
@@ -745,7 +744,7 @@ class TestLatlonBoundsGCA(TestCase):
         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(np.array([vertices_cart[1], vertices_cart[2]]), extreme_type="min"))
+                      _extreme_gca_latitude_cartesian(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)
@@ -773,7 +772,7 @@ class TestLatlonBoundsGCA(TestCase):
         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(np.array([vertices_cart[1], vertices_cart[2]]), extreme_type="min"))
+                      _extreme_gca_latitude_cartesian(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)
@@ -801,7 +800,7 @@ class TestLatlonBoundsGCA(TestCase):
         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(np.array([vertices_cart[1], vertices_cart[2]]), extreme_type="min"))
+                      _extreme_gca_latitude_cartesian(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)
@@ -1396,7 +1395,10 @@ class TestLatlonBoundsMix(TestCase):
                            [[np.deg2rad(60.0), np.pi / 2], [0., 2 * np.pi]]]
 
         grid = ux.Grid.from_face_vertices(faces, latlon=True)
-        bounds_xarray = _populate_bounds(grid, is_face_GCA_list=[[True, False, True, False]] * 4, return_array=True)
+        bounds_xarray = _populate_bounds(grid, is_face_GCA_list=np.array([[True, False, True, False],
+                                                                          [True, False, True, False],
+                                                                          [True, False, True, False],
+                                                                          [True, False, True, False]]) , return_array=True)
         face_bounds = bounds_xarray.values
         for i in range(len(faces)):
             nt.assert_allclose(face_bounds[i], expected_bounds[i], atol=ERROR_TOLERANCE)
diff --git a/test/test_helpers.py b/test/test_helpers.py
index c5b923a2..da9b8c3e 100644
--- a/test/test_helpers.py
+++ b/test/test_helpers.py
@@ -15,7 +15,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.geometry import _pole_point_inside_polygon
+from uxarray.grid.geometry import pole_point_inside_polygon, _pole_point_inside_polygon_cartesian
 
 try:
     import constants
@@ -267,7 +267,7 @@ class TestFaceEdgeConnectivityHelper(TestCase):
         )
 
         # Check that the face_edges_connectivity_cartesian works as an input to _pole_point_inside_polygon
-        result = ux.grid.geometry._pole_point_inside_polygon(
+        result = _pole_point_inside_polygon_cartesian(
             'North', face_edges_connectivity_cartesian[0]
         )
 
@@ -300,7 +300,7 @@ class TestFaceEdgeConnectivityHelper(TestCase):
         )
 
         # Check that the face_edges_connectivity_cartesian works as an input to _pole_point_inside_polygon
-        result = ux.grid.geometry._pole_point_inside_polygon(
+        result = _pole_point_inside_polygon_cartesian(
             'North', face_edges_connectivity_cartesian[0]
         )
 
@@ -422,7 +422,7 @@ class TestFaceEdgeConnectivityHelper(TestCase):
             face_edges_connectivity_cartesian.append(edge_cart)
 
         # Check that the face_edges_connectivity_cartesian works as an input to _pole_point_inside_polygon
-        result = ux.grid.geometry._pole_point_inside_polygon(
+        result = _pole_point_inside_polygon_cartesian(
             'North', np.array(face_edges_connectivity_cartesian)
         )
 
@@ -461,7 +461,7 @@ class TestFaceEdgeConnectivityHelper(TestCase):
             face_edges_connectivity_cartesian.append(edge_cart)
 
         # Check that the face_edges_connectivity_cartesian works as an input to _pole_point_inside_polygon
-        result = ux.grid.geometry._pole_point_inside_polygon(
+        result = _pole_point_inside_polygon_cartesian(
             'North', np.array(face_edges_connectivity_cartesian)
         )
 
diff --git a/test/test_intersections.py b/test/test_intersections.py
index 0f3c6260..c28c5e0b 100644
--- a/test/test_intersections.py
+++ b/test/test_intersections.py
@@ -5,8 +5,9 @@ from uxarray.constants import ERROR_TOLERANCE
 
 # from uxarray.grid.coordinates import node_lonlat_rad_to_xyz, node_xyz_to_lonlat_rad
 
+from uxarray.grid.arcs import _extreme_gca_latitude_cartesian
 from uxarray.grid.coordinates import _lonlat_rad_to_xyz, _xyz_to_lonlat_rad
-from uxarray.grid.intersections import gca_gca_intersection, gca_const_lat_intersection
+from uxarray.grid.intersections import gca_gca_intersection, gca_const_lat_intersection, _gca_gca_intersection_cartesian
 
 
 class TestGCAGCAIntersection(TestCase):
@@ -24,10 +25,10 @@ class TestGCAGCAIntersection(TestCase):
             _lonlat_rad_to_xyz(np.deg2rad(70.0), 0.0),
             _lonlat_rad_to_xyz(np.deg2rad(179.0), 0.0)
         ])
-        res_cart = gca_gca_intersection(GCR1_cart, GCR2_cart)
+        res_cart = _gca_gca_intersection_cartesian(GCR1_cart, GCR2_cart)
 
         # res_cart should be empty since these two GCRs are not intersecting
-        self.assertTrue(np.array_equal(res_cart, np.array([])))
+        self.assertTrue(len(res_cart) == 0)
 
         GCR1_cart = np.array([
             _lonlat_rad_to_xyz(np.deg2rad(170.0),
@@ -40,7 +41,7 @@ class TestGCAGCAIntersection(TestCase):
             _lonlat_rad_to_xyz(np.deg2rad(175.0), 0.0)
         ])
 
-        res_cart = gca_gca_intersection(GCR1_cart, GCR2_cart)
+        res_cart = _gca_gca_intersection_cartesian(GCR1_cart, GCR2_cart)
         res_cart = res_cart[0]
 
         # Test if the result is normalized
@@ -66,7 +67,7 @@ class TestGCAGCAIntersection(TestCase):
             _lonlat_rad_to_xyz(0.5 * np.pi, 0.0),
             _lonlat_rad_to_xyz(-0.5 * np.pi - 0.01, 0.0)
         ])
-        res_cart = gca_gca_intersection(GCR1_cart, GCR2_cart)
+        res_cart = _gca_gca_intersection_cartesian(GCR1_cart, GCR2_cart)
         res_cart = res_cart[0]
         expected_res = np.array(_lonlat_rad_to_xyz(0.5 * np.pi, 0.0))
         # Test if two results are equal within the error tolerance
@@ -84,7 +85,7 @@ class TestGCAGCAIntersection(TestCase):
             _lonlat_rad_to_xyz(*[0.5 * np.pi, 0.0]),
             _lonlat_rad_to_xyz(*[-0.5 * np.pi - 0.01, 0.0])
         ])
-        res_cart = gca_gca_intersection(GCR1_cart, GCR2_cart)
+        res_cart = _gca_gca_intersection_cartesian(GCR1_cart, GCR2_cart)
         res_cart = res_cart[0]
         # Test if the result is normalized
         self.assertTrue(
@@ -133,7 +134,7 @@ class TestGCAconstLatIntersection(TestCase):
             _lonlat_rad_to_xyz(np.deg2rad(170.0),
                                     np.deg2rad(10.0))
         ])
-        max_lat = ux.grid.arcs.extreme_gca_latitude(GCR1_cart, 'max')
+        max_lat = _extreme_gca_latitude_cartesian(GCR1_cart, 'max')
 
         query_lat = (np.deg2rad(10.0) + max_lat) / 2.0