diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 026408959..c8c1c5f00 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -14,7 +14,7 @@ repos:
 
 - repo: https://github.com/astral-sh/ruff-pre-commit
   # Ruff version.
-  rev: v0.8.2
+  rev: v0.8.3
   hooks:
     # Run the linter.
     - id: ruff
diff --git a/ci/asv.yml b/ci/asv.yml
index cbe442b33..2c4089646 100644
--- a/ci/asv.yml
+++ b/ci/asv.yml
@@ -17,6 +17,7 @@ dependencies:
   - pandas
   - pathlib
   - pre_commit
+  - polars
   - pyarrow
   - pytest
   - pytest-cov
diff --git a/ci/docs.yml b/ci/docs.yml
index a7411af73..6fcac9726 100644
--- a/ci/docs.yml
+++ b/ci/docs.yml
@@ -34,6 +34,7 @@ dependencies:
     - pandas
     - geocat-datafiles
     - spatialpandas
+    - polars
     - geopandas
     - pip:
       - antimeridian
diff --git a/ci/environment.yml b/ci/environment.yml
index a696e397a..32cc0c8e9 100644
--- a/ci/environment.yml
+++ b/ci/environment.yml
@@ -21,6 +21,7 @@ dependencies:
   - pandas
   - pathlib
   - pre_commit
+  - polars
   - pyarrow
   - pytest
   - pip
diff --git a/pyproject.toml b/pyproject.toml
index 169689fcd..6a0ed0b37 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -41,6 +41,7 @@ dependencies = [
   "geopandas",
   "xarray",
   "hvplot",
+  "polars",
 ]
 # minimal dependencies end
 
diff --git a/test/test_arcs.py b/test/test_arcs.py
index c139e3907..4e9262982 100644
--- a/test/test_arcs.py
+++ b/test/test_arcs.py
@@ -10,7 +10,7 @@
 import uxarray as ux
 
 from uxarray.grid.coordinates import _lonlat_rad_to_xyz
-from uxarray.grid.arcs import point_within_gca, _point_within_gca_cartesian
+from uxarray.grid.arcs import point_within_gca
 
 try:
     import constants
@@ -31,48 +31,31 @@ class TestIntersectionPoint(TestCase):
     def test_pt_within_gcr(self):
         # The GCR that's eexactly 180 degrees will have Value Error raised
 
-        gcr_180degree_cart = [
+        gcr_180degree_cart = np.asarray([
             _lonlat_rad_to_xyz(0.0, np.pi / 2.0),
             _lonlat_rad_to_xyz(0.0, -np.pi / 2.0)
-        ]
+        ])
+        pt_same_lon_in = np.asarray(_lonlat_rad_to_xyz(0.0, 0.0))
 
-        pt_same_lon_in = _lonlat_rad_to_xyz(0.0, 0.0)
         with self.assertRaises(ValueError):
-            _point_within_gca_cartesian(pt_same_lon_in, gcr_180degree_cart)
-
+            point_within_gca(pt_same_lon_in, gcr_180degree_cart[0],gcr_180degree_cart[1] )
+        #
         # Test when the point and the GCA all have the same longitude
-        gcr_same_lon_cart = [
+        gcr_same_lon_cart = np.asarray([
             _lonlat_rad_to_xyz(0.0, 1.5),
             _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_cartesian(pt_same_lon_in, gcr_same_lon_cart))
+        ])
+        pt_same_lon_in = np.asarray(_lonlat_rad_to_xyz(0.0, 0.0))
+        self.assertTrue(point_within_gca(pt_same_lon_in, gcr_same_lon_cart[0], gcr_same_lon_cart[1]))
 
-        pt_same_lon_out = _lonlat_rad_to_xyz(0.0, 1.500000000000001)
-        res = _point_within_gca_cartesian(pt_same_lon_out, gcr_same_lon_cart)
+        pt_same_lon_out = np.asarray(_lonlat_rad_to_xyz(0.0, 1.5000001))
+        res = point_within_gca(pt_same_lon_out, gcr_same_lon_cart[0], gcr_same_lon_cart[1])
         self.assertFalse(res)
 
-        pt_same_lon_out_2 = _lonlat_rad_to_xyz(0.1, 1.0)
-        res = _point_within_gca_cartesian(pt_same_lon_out_2, gcr_same_lon_cart)
+        pt_same_lon_out_2 = np.asarray(_lonlat_rad_to_xyz(0.1, 1.0))
+        res = point_within_gca(pt_same_lon_out_2, gcr_same_lon_cart[0], gcr_same_lon_cart[1])
         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_cartesian(pt_same_lon_out_add_one_place, gcr_same_lon_cart)
-        self.assertTrue(res)
-
-        # Normal case
-        # GCR vertex0 in radian : [1.3003315590159483, -0.007004587172323237],
-        # GCR vertex1 in radian : [3.5997458123873827, -1.4893379576608758]
-        # Point in radian : [1.3005410084914981, -0.010444274637648326]
-        gcr_cart_2 = np.array([[0.267, 0.963, -0.007], [-0.073, -0.036,
-                                                        -0.997]])
-        pt_cart_within = np.array(
-            [0.25616109352676675, 0.9246590335292105, -0.010021496695000144])
-        self.assertTrue(_point_within_gca_cartesian(pt_cart_within, gcr_cart_2, True))
-
-        # Test other more complicate cases : The anti-meridian case
-
     def test_pt_within_gcr_antimeridian(self):
         # GCR vertex0 in radian : [5.163808182822441, 0.6351384888657234],
         # GCR vertex1 in radian : [0.8280410325693055, 0.42237025187091526]
@@ -80,16 +63,14 @@ def test_pt_within_gcr_antimeridian(self):
         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_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
+        self.assertTrue(
+            point_within_gca(pt_cart, gcr_cart[0], gcr_cart[1]))
+
         gcr_cart_flip = np.array([[0.617, 0.672, 0.410], [0.351, -0.724,
                                                           0.593]])
-        with self.assertRaises(ValueError):
-            _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_cartesian(pt_cart, gcr_cart_flip, is_directed=False))
+            point_within_gca(pt_cart, gcr_cart_flip[0], gcr_cart_flip[1]))
 
         # 2nd anti-meridian case
         # GCR vertex0 in radian : [4.104711496596806, 0.5352983676533828],
@@ -100,22 +81,9 @@ def test_pt_within_gcr_antimeridian(self):
         pt_cart_within = np.array(
             [0.6136726305712109, 0.28442243941920053, -0.365605190899831])
         self.assertFalse(
-            _point_within_gca_cartesian(pt_cart_within, gcr_cart_1, is_directed=True))
-        self.assertFalse(
-            _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]
-        v2_rad = [2 * np.pi - 0.1, 0.0]
-        v1_cart = _lonlat_rad_to_xyz(v1_rad[0], v1_rad[1])
-        v2_cart = _lonlat_rad_to_xyz(v2_rad[0], v1_rad[1])
-        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_cartesian(pt_cart, gcr_cart, is_directed=True)
-        gcr_car_flipped = np.array([v2_cart, v1_cart])
-        self.assertTrue(
-            _point_within_gca_cartesian(pt_cart, gcr_car_flipped, is_directed=True))
+            point_within_gca(pt_cart_within, gcr_cart_1[0], gcr_cart_1[1]))
+
+
 
     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,15 +93,4 @@ def test_pt_within_gcr_cross_pole(self):
         # 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_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(
-            [0.10, 0.0, 0.8])
-
-        # 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_cartesian(pt_cart, gcr_cart, is_directed=False))
-        with self.assertRaises(ValueError):
-            _point_within_gca_cartesian(pt_cart, gcr_cart, is_directed=True)
+        self.assertTrue(point_within_gca(pt_cart, gcr_cart[0], gcr_cart[1]))
diff --git a/test/test_helpers.py b/test/test_helpers.py
index da9b8c3eb..491351e14 100644
--- a/test/test_helpers.py
+++ b/test/test_helpers.py
@@ -238,6 +238,31 @@ def test_angle_of_2_vectors(self):
         v2 = np.array([1.0, 0.0, 0.0])
         self.assertAlmostEqual(_angle_of_2_vectors(v1, v2), 0.0)
 
+    def test_angle_of_2_vectors_180_degree(self):
+        GCR1_cart = np.array([
+            _lonlat_rad_to_xyz(np.deg2rad(0.0),
+                                    np.deg2rad(0.0)),
+            _lonlat_rad_to_xyz(np.deg2rad(181.0),
+                                    np.deg2rad(0.0))
+        ])
+
+        res = _angle_of_2_vectors( GCR1_cart[0],  GCR1_cart[1])
+
+        # The angle between the two vectors should be 181 degree
+        self.assertAlmostEqual(res, np.deg2rad(181.0), places=8)
+
+        GCR1_cart = np.array([
+            _lonlat_rad_to_xyz(np.deg2rad(170.0),
+                                    np.deg2rad(89.0)),
+            _lonlat_rad_to_xyz(np.deg2rad(170.0),
+                                    np.deg2rad(-10.0))
+        ])
+
+        res = _angle_of_2_vectors( GCR1_cart[0],  GCR1_cart[1])
+
+        # The angle between the two vectors should be 181 degree
+        self.assertAlmostEqual(res, np.deg2rad(89.0+10.0), places=8)
+
 
 class TestFaceEdgeConnectivityHelper(TestCase):
 
diff --git a/test/test_integrate.py b/test/test_integrate.py
index 8fa5b2eb8..db327959e 100644
--- a/test/test_integrate.py
+++ b/test/test_integrate.py
@@ -81,10 +81,9 @@ def test_get_faces_constLat_intersection_info_one_intersection(self):
         ])
 
         latitude_cart = -0.8660254037844386
-        is_directed=False
         is_latlonface=False
         is_GCA_list=None
-        unique_intersections, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, is_GCA_list, is_latlonface, is_directed)
+        unique_intersections, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, is_GCA_list, is_latlonface)
         # The expected unique_intersections length is 1
         self.assertEqual(len(unique_intersections), 1)
 
@@ -111,10 +110,10 @@ def test_get_faces_constLat_intersection_info_encompass_pole(self):
         latitude_rad = np.arcsin(latitude_cart)
         latitude_deg = np.rad2deg(latitude_rad)
         print(latitude_deg)
-        is_directed=False
+
         is_latlonface=False
         is_GCA_list=None
-        unique_intersections, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, is_GCA_list, is_latlonface, is_directed)
+        unique_intersections, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, is_GCA_list, is_latlonface)
         # The expected unique_intersections length should be no greater than 2* n_edges
         self.assertLessEqual(len(unique_intersections), 2*len(face_edges_cart))
 
@@ -133,10 +132,9 @@ def test_get_faces_constLat_intersection_info_on_pole(self):
              [-5.2264427688714095e-02, -5.2264427688714102e-02, -9.9726468863423734e-01]]
         ])
         latitude_cart = -0.9998476951563913
-        is_directed=False
         is_latlonface=False
         is_GCA_list=None
-        unique_intersections, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, is_GCA_list, is_latlonface, is_directed)
+        unique_intersections, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, is_GCA_list, is_latlonface)
         # The expected unique_intersections length is 2
         self.assertEqual(len(unique_intersections), 2)
 
@@ -156,10 +154,9 @@ def test_get_faces_constLat_intersection_info_near_pole(self):
         latitude_cart = -0.9876883405951378
         latitude_rad = np.arcsin(latitude_cart)
         latitude_deg = np.rad2deg(latitude_rad)
-        is_directed=False
         is_latlonface=False
         is_GCA_list=None
-        unique_intersections, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, is_GCA_list, is_latlonface, is_directed)
+        unique_intersections, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, is_GCA_list, is_latlonface)
         # The expected unique_intersections length is 2
         self.assertEqual(len(unique_intersections), 1)
 
@@ -182,10 +179,9 @@ def test_get_faces_constLat_intersection_info_2(self):
                                      [0.6546536707079771, -0.37796447300922714, -0.6546536707079772]]])
 
         latitude_cart = -0.6560590289905073
-        is_directed=False
         is_latlonface=False
         is_GCA_list=None
-        unique_intersections, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, is_GCA_list, is_latlonface, is_directed)
+        unique_intersections, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, is_GCA_list, is_latlonface)
         # The expected unique_intersections length is 2
         self.assertEqual(len(unique_intersections), 2)
 
@@ -208,10 +204,9 @@ def test_get_faces_constLat_intersection_info_2(self):
                                      [0.6546536707079771, -0.37796447300922714, -0.6546536707079772]]])
 
         latitude_cart = -0.6560590289905073
-        is_directed=False
         is_latlonface=False
         is_GCA_list=None
-        unique_intersections, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, is_GCA_list, is_latlonface, is_directed)
+        unique_intersections, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, is_GCA_list, is_latlonface)
         # The expected unique_intersections length is 2
         self.assertEqual(len(unique_intersections), 2)
 
@@ -240,8 +235,7 @@ def test_get_zonal_face_interval(self):
         # The latlon bounds for the latitude is not necessarily correct below since we don't use the latitudes bound anyway
         interval_df = _get_zonal_face_interval(face_edge_nodes, constZ,
                                                np.array([[-0.25 * np.pi, 0.25 * np.pi], [1.6 * np.pi,
-                                                                                         0.4 * np.pi]]),
-                                               is_directed=False)
+                                                                                         0.4 * np.pi]]))
         expected_interval_df = pd.DataFrame({
             'start': [1.6 * np.pi, 0.0],
             'end': [2.0 * np.pi, 00.4 * np.pi]
@@ -285,7 +279,7 @@ def test_get_zonal_face_interval_empty_interval(self):
             [3.14159265, 3.2321175]
         ])
 
-        res = _get_zonal_face_interval(face_edges_cart, latitude_cart, face_latlon_bounds, is_directed=False)
+        res = _get_zonal_face_interval(face_edges_cart, latitude_cart, face_latlon_bounds)
         expected_res = pd.DataFrame({"start": [0.0], "end": [0.0]})
         pd.testing.assert_frame_equal(res, expected_res)
 
@@ -322,7 +316,7 @@ def test_get_zonal_face_interval_encompass_pole(self):
         })
 
         # Call the function to get the result
-        res = _get_zonal_face_interval(face_edges_cart, latitude_cart, face_latlon_bounds, is_directed=False)
+        res = _get_zonal_face_interval(face_edges_cart, latitude_cart, face_latlon_bounds)
 
         # Assert the result matches the expected DataFrame
         pd.testing.assert_frame_equal(res,  expected_df)
@@ -349,8 +343,7 @@ def test_get_zonal_face_interval_FILL_VALUE(self):
         # The latlon bounds for the latitude is not necessarily correct below since we don't use the latitudes bound anyway
         interval_df = _get_zonal_face_interval(face_edge_nodes, constZ,
                                                np.array([[-0.25 * np.pi, 0.25 * np.pi], [1.6 * np.pi,
-                                                                                         0.4 * np.pi]]),
-                                               is_directed=False)
+                                                                                         0.4 * np.pi]]))
         expected_interval_df = pd.DataFrame({
             'start': [1.6 * np.pi, 0.0],
             'end': [2.0 * np.pi, 00.4 * np.pi]
@@ -383,7 +376,7 @@ def test_get_zonal_face_interval_GCA_constLat(self):
         interval_df = _get_zonal_face_interval(face_edge_nodes, constZ,
                                                np.array([[-0.25 * np.pi, 0.25 * np.pi], [1.6 * np.pi,
                                                                                          0.4 * np.pi]]),
-                                               is_directed=False, is_GCA_list=np.array([True, False, True, False]))
+                                               is_GCA_list=np.array([True, False, True, False]))
         expected_interval_df = pd.DataFrame({
             'start': [1.6 * np.pi, 0.0],
             'end': [2.0 * np.pi, 00.4 * np.pi]
@@ -415,7 +408,7 @@ def test_get_zonal_face_interval_equator(self):
         interval_df = _get_zonal_face_interval(face_edge_nodes, 0.0,
                                                np.array([[-0.25 * np.pi, 0.25 * np.pi], [1.6 * np.pi,
                                                                                          0.4 * np.pi]]),
-                                               is_directed=False, is_GCA_list=np.array([True, True, True, True]))
+                                               is_GCA_list=np.array([True, True, True, True]))
         expected_interval_df = pd.DataFrame({
             'start': [1.6 * np.pi, 0.0],
             'end': [2.0 * np.pi, 00.4 * np.pi]
@@ -434,7 +427,7 @@ def test_get_zonal_face_interval_equator(self):
         interval_df = _get_zonal_face_interval(face_edge_nodes, 0.0,
                                                np.array([[-0.25 * np.pi, 0.25 * np.pi], [1.6 * np.pi,
                                                                                          0.4 * np.pi]]),
-                                               is_directed=False, is_GCA_list=np.array([True, False, True, False]))
+                                               is_GCA_list=np.array([True, False, True, False]))
         expected_interval_df = pd.DataFrame({
             'start': [1.6 * np.pi, 0.0],
             'end': [2.0 * np.pi, 00.4 * np.pi]
@@ -605,7 +598,7 @@ def test_get_zonal_faces_weight_at_constLat_equator(self):
         weight_df = _get_zonal_faces_weight_at_constLat(np.array([
             face_0_edge_nodes, face_1_edge_nodes, face_2_edge_nodes,
             face_3_edge_nodes
-        ]), 0.0, latlon_bounds, is_directed=False)
+        ]), 0.0, latlon_bounds)
 
         nt.assert_array_almost_equal(weight_df, expected_weight_df, decimal=3)
 
@@ -666,7 +659,7 @@ def test_get_zonal_faces_weight_at_constLat_regular(self):
         weight_df = _get_zonal_faces_weight_at_constLat(np.array([
             face_0_edge_nodes, face_1_edge_nodes, face_2_edge_nodes,
             face_3_edge_nodes
-        ]), np.sin(0.1 * np.pi), latlon_bounds, is_directed=False)
+        ]), np.sin(0.1 * np.pi), latlon_bounds)
 
         nt.assert_array_almost_equal(weight_df, expected_weight_df, decimal=3)
 
@@ -693,7 +686,7 @@ def test_get_zonal_faces_weight_at_constLat_on_pole_one_face(self):
         ])
         constLat_cart = -1
 
-        weight_df = _get_zonal_faces_weight_at_constLat(face_edges_cart, constLat_cart, face_bounds, is_directed=False)
+        weight_df = _get_zonal_faces_weight_at_constLat(face_edges_cart, constLat_cart, face_bounds)
         # Define the expected DataFrame
         expected_weight_df = pd.DataFrame({"face_index": [0], "weight": [1.0]})
 
@@ -740,7 +733,7 @@ def test_get_zonal_faces_weight_at_constLat_on_pole_faces(self):
 
         constLat_cart = 1.0
 
-        weight_df = _get_zonal_faces_weight_at_constLat(face_edges_cart, constLat_cart, face_bounds, is_directed=False)
+        weight_df = _get_zonal_faces_weight_at_constLat(face_edges_cart, constLat_cart, face_bounds)
         # Define the expected DataFrame
         expected_weight_df = pd.DataFrame({
             'face_index': [0, 1, 2, 3],
@@ -774,7 +767,7 @@ def test_get_zonal_face_interval_pole(self):
         ])
         constLat_cart = -0.9986295347545738
 
-        weight_df = _get_zonal_face_interval(face_edges_cart, constLat_cart, face_bounds, is_directed=False)
+        weight_df = _get_zonal_face_interval(face_edges_cart, constLat_cart, face_bounds)
         # No Nan values should be present in the weight_df
         self.assertFalse(weight_df.isnull().values.any())
 
@@ -827,7 +820,7 @@ def test_get_zonal_faces_weight_at_constLat_latlonface(self):
         # Assert the results is the same to the 3 decimal places
         weight_df = _get_zonal_faces_weight_at_constLat(np.array([
             face_0_edge_nodes, face_1_edge_nodes, face_2_edge_nodes
-        ]), np.sin(np.deg2rad(20)), latlon_bounds, is_directed=False, is_latlonface=True)
+        ]), np.sin(np.deg2rad(20)), latlon_bounds, is_latlonface=True)
 
 
         nt.assert_array_almost_equal(weight_df, expected_weight_df, decimal=3)
@@ -839,4 +832,4 @@ def test_get_zonal_faces_weight_at_constLat_latlonface(self):
         with self.assertRaises(ValueError):
             _get_zonal_faces_weight_at_constLat(np.array([
                 face_0_edge_nodes, face_1_edge_nodes, face_2_edge_nodes
-            ]), np.deg2rad(20), latlon_bounds, is_directed=False)
+            ]), np.deg2rad(20), latlon_bounds)
diff --git a/test/test_intersections.py b/test/test_intersections.py
index c28c5e0bc..a9051d466 100644
--- a/test/test_intersections.py
+++ b/test/test_intersections.py
@@ -6,7 +6,7 @@
 # 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.coordinates import _lonlat_rad_to_xyz, _xyz_to_lonlat_rad,_xyz_to_lonlat_rad_scalar
 from uxarray.grid.intersections import gca_gca_intersection, gca_const_lat_intersection, _gca_gca_intersection_cartesian
 
 
@@ -14,7 +14,6 @@ class TestGCAGCAIntersection(TestCase):
 
     def test_get_GCA_GCA_intersections_antimeridian(self):
         # Test the case where the two GCAs are on the antimeridian
-        GCA1 = _lonlat_rad_to_xyz(np.deg2rad(170.0), np.deg2rad(89.99))
         GCR1_cart = np.array([
             _lonlat_rad_to_xyz(np.deg2rad(170.0),
                                     np.deg2rad(89.99)),
@@ -36,6 +35,8 @@ def test_get_GCA_GCA_intersections_antimeridian(self):
             _lonlat_rad_to_xyz(np.deg2rad(170.0),
                                     np.deg2rad(-10.0))
         ])
+
+
         GCR2_cart = np.array([
             _lonlat_rad_to_xyz(np.deg2rad(70.0), 0.0),
             _lonlat_rad_to_xyz(np.deg2rad(175.0), 0.0)
@@ -82,21 +83,172 @@ def test_get_GCA_GCA_intersections_perpendicular(self):
                                     np.deg2rad(10.0))
         ])
         GCR2_cart = np.array([
-            _lonlat_rad_to_xyz(*[0.5 * np.pi, 0.0]),
-            _lonlat_rad_to_xyz(*[-0.5 * np.pi - 0.01, 0.0])
+            _lonlat_rad_to_xyz(*[0.5 * np.pi - 0.01, 0.0]),
+            _lonlat_rad_to_xyz(*[-0.5 * np.pi + 0.01, 0.0])
         ])
         res_cart = _gca_gca_intersection_cartesian(GCR1_cart, GCR2_cart)
-        res_cart = res_cart[0]
-        # Test if the result is normalized
-        self.assertTrue(
-            np.allclose(np.linalg.norm(res_cart, axis=0),
-                        1.0,
-                        atol=ERROR_TOLERANCE))
-        res_lonlat_rad = _xyz_to_lonlat_rad(*res_cart)
-        self.assertTrue(
-            np.allclose(res_lonlat_rad,
-                        np.array([np.deg2rad(170.0),
-                                  np.deg2rad(0.0)])))
+
+        # rest_cart should be empty since these two GCAs are not intersecting
+        self.assertTrue(len(res_cart) == 0)
+
+
+    # def test_GCA_GCA_single_edge_to_pole(self):
+    #     # GCA_a - Face Center connected to South Pole
+    #     # Point A - South Pole
+    #     ref_point_lonlat = np.deg2rad(np.array([0.0, -90.0]))
+    #     ref_point_xyz = np.array(_lonlat_rad_to_xyz(*ref_point_lonlat))
+    #     # Point B - Face Center
+    #     face_lonlat = np.deg2rad(np.array([-175, 26.5]))
+    #     face_xyz = np.array(_lonlat_rad_to_xyz(*face_lonlat))
+    #     gca_a_xyz = np.array([face_xyz, ref_point_xyz])
+    #
+    #     # GCA_b - Single Face Edge
+    #     # Point A - First Edge Point
+    #     edge_a_lonlat = np.deg2rad(np.array((-175, -24.5)))
+    #     edge_b_lonlat = np.deg2rad(np.array((-173, 28.7)))
+    #
+    #     # Point B - Second Edge Point
+    #     edge_a_xyz = np.array(_lonlat_rad_to_xyz(*edge_a_lonlat))
+    #     edge_b_xyz = np.array(_lonlat_rad_to_xyz(*edge_b_lonlat))
+    #     gca_b_xyz = np.array([edge_a_xyz, edge_b_xyz])
+    #
+    #     # The edge should intersect
+    #     self.assertTrue(len(gca_gca_intersection(gca_a_xyz, gca_b_xyz)))
+
+    def test_GCA_GCA_south_pole(self):
+
+        # GCA_a - Face Center connected to South Pole
+        # Point A - South Pole
+        ref_point_lonlat = np.deg2rad(np.array([0.0, -90.0]))
+        ref_point_xyz = np.array(_lonlat_rad_to_xyz(*ref_point_lonlat))
+        # Point B - Face Center
+        face_lonlat = np.deg2rad(np.array([0.0, 0.0]))
+        face_xyz = np.array(_lonlat_rad_to_xyz(*face_lonlat))
+        gca_a_xyz = np.array([face_xyz, ref_point_xyz])
+
+        # GCA_b - Single Face Edge
+        # Point A - First Edge Point
+        edge_a_lonlat = np.deg2rad(np.array((-45, -1.0)))
+        edge_b_lonlat = np.deg2rad(np.array((45, -1.0)))
+
+        # Point B - Second Edge Point
+        edge_a_xyz = np.array(_lonlat_rad_to_xyz(*edge_a_lonlat))
+        edge_b_xyz = np.array(_lonlat_rad_to_xyz(*edge_b_lonlat))
+        gca_b_xyz = np.array([edge_a_xyz, edge_b_xyz])
+
+        # The edge should intersect
+        self.assertTrue(len(gca_gca_intersection(gca_a_xyz, gca_b_xyz)))
+
+    def test_GCA_GCA_north_pole(self):
+        # GCA_a - Face Center connected to South Pole
+        ref_point_lonlat = np.deg2rad(np.array([0.0, 90.0]))
+        ref_point_xyz = np.array(_lonlat_rad_to_xyz(*ref_point_lonlat))
+        face_lonlat = np.deg2rad(np.array([0.0, 0.0]))
+        face_xyz = np.array(_lonlat_rad_to_xyz(*face_lonlat))
+        gca_a_xyz = np.array([face_xyz, ref_point_xyz])
+
+        # GCA_b - Single Face Edge
+        edge_a_lonlat = np.deg2rad(np.array((-45, 1.0)))
+        edge_b_lonlat = np.deg2rad(np.array((45, 1.0)))
+
+        edge_a_xyz = np.array(_lonlat_rad_to_xyz(*edge_a_lonlat))
+        edge_b_xyz = np.array(_lonlat_rad_to_xyz(*edge_b_lonlat))
+        gca_b_xyz = np.array([edge_a_xyz, edge_b_xyz])
+
+        # The edge should intersect
+        self.assertTrue(len(gca_gca_intersection(gca_a_xyz, gca_b_xyz)))
+
+    def test_GCA_GCA_north_pole_angled(self):
+        # GCA_a
+        ref_point_lonlat = np.deg2rad(np.array([0.0, 90.0]))
+        ref_point_xyz = np.array(_lonlat_rad_to_xyz(*ref_point_lonlat))
+        face_lonlat = np.deg2rad(np.array([-45.0, 45.0]))
+        face_xyz = np.array(_lonlat_rad_to_xyz(*face_lonlat))
+        gca_a_xyz = np.array([face_xyz, ref_point_xyz])
+
+        # GCA_b
+        edge_a_lonlat = np.deg2rad(np.array((-45.0, 50.0)))
+        edge_b_lonlat = np.deg2rad(np.array((-40.0, 45.0)))
+
+        # Point B - Second Edge Point
+        edge_a_xyz = np.array(_lonlat_rad_to_xyz(*edge_a_lonlat))
+        edge_b_xyz = np.array(_lonlat_rad_to_xyz(*edge_b_lonlat))
+        gca_b_xyz = np.array([edge_a_xyz, edge_b_xyz])
+
+        # The edge should intersect
+        self.assertTrue(len(gca_gca_intersection(gca_a_xyz, gca_b_xyz)))
+
+
+    def test_GCA_edge_intersection_count(self):
+
+        from uxarray.grid.utils import _get_cartesian_face_edge_nodes
+
+        # Generate a normal face that is not crossing the antimeridian or the poles
+        vertices_lonlat = [[29.5, 11.0], [29.5, 10.0], [30.5, 10.0], [30.5, 11.0]]
+        vertices_lonlat = np.array(vertices_lonlat)
+
+        grid = ux.Grid.from_face_vertices(vertices_lonlat, latlon=True)
+        face_edge_nodes_cartesian = _get_cartesian_face_edge_nodes(
+            grid.face_node_connectivity.values,
+            grid.n_face,
+            grid.n_max_face_edges,
+            grid.node_x.values,
+            grid.node_y.values,
+            grid.node_z.values)
+
+        face_center_xyz = np.array([grid.face_x.values[0], grid.face_y.values[0], grid.face_z.values[0]], dtype=np.float64)
+        north_pole_xyz = np.array([0.0, 0.0, 1.0], dtype=np.float64)
+        south_pole_xyz = np.array([0.0, 0.0, -1.0], dtype=np.float64)
+
+        gca_face_center_north_pole = np.array([face_center_xyz, north_pole_xyz], dtype=np.float64)
+        gca_face_center_south_pole = np.array([face_center_xyz, south_pole_xyz], dtype=np.float64)
+
+        intersect_north_pole_count = 0
+        intersect_south_pole_count = 0
+
+        for edge in face_edge_nodes_cartesian[0]:
+            res1 = gca_gca_intersection(edge, gca_face_center_north_pole)
+            res2 = gca_gca_intersection(edge, gca_face_center_south_pole)
+
+            if len(res1):
+                intersect_north_pole_count += 1
+            if len(res2):
+                intersect_south_pole_count += 1
+
+        print(intersect_north_pole_count, intersect_south_pole_count)
+        self.assertTrue(intersect_north_pole_count == 1)
+        self.assertTrue(intersect_south_pole_count == 1)
+
+    def test_GCA_GCA_single_edge_to_pole(self):
+        # GCA_a - Face Center connected to South Pole
+        # Point A - South Pole
+        ref_point_lonlat_exact = np.deg2rad(np.array([0.0, -90.0]))
+        ref_point_lonlat_close = np.deg2rad(np.array([0.0, -89.99999]))
+        ref_point_xyz_exact = np.array(_lonlat_rad_to_xyz(*ref_point_lonlat_exact))
+        ref_point_xyz_close = np.array(_lonlat_rad_to_xyz(*ref_point_lonlat_close))
+
+        # Point B - Face Center
+        face_lonlat = np.deg2rad(np.array([-175.0, 26.5]))
+        face_xyz = np.array(_lonlat_rad_to_xyz(*face_lonlat))
+        gca_a_xyz_close = np.array([face_xyz, ref_point_xyz_close])
+        gca_a_xyz_exact = np.array([face_xyz, ref_point_xyz_exact])
+
+        # GCA_b - Single Face Edge
+        # Point A - First Edge Point
+        edge_a_lonlat = np.deg2rad(np.array((-175.0, -24.5)))
+        edge_b_lonlat = np.deg2rad(np.array((-173.0, 28.7)))
+
+        # Point B - Second Edge Point
+        edge_a_xyz = np.array(_lonlat_rad_to_xyz(*edge_a_lonlat))
+        edge_b_xyz = np.array(_lonlat_rad_to_xyz(*edge_b_lonlat))
+        gca_b_xyz = np.array([edge_a_xyz, edge_b_xyz])
+
+        # The edge should intersect
+        self.assertTrue(len(gca_gca_intersection(gca_a_xyz_close, gca_b_xyz)))
+        self.assertTrue(len(gca_gca_intersection(gca_a_xyz_exact, gca_b_xyz)))
+
+
+
 
 
 class TestGCAconstLatIntersection(TestCase):
diff --git a/uxarray/core/utils.py b/uxarray/core/utils.py
index 7b2292308..fdd38c5fa 100644
--- a/uxarray/core/utils.py
+++ b/uxarray/core/utils.py
@@ -39,16 +39,21 @@ def _map_dims_to_ugrid(
             # drop dimensions not present in the original dataset
             _source_dims_dict.pop(key)
 
+        # only check edge dimension if it is present (to avoid overhead of computing connectivity)
+        if "n_edge" in grid._ds.dims:
+            n_edge = grid._ds.sizes["n_edge"]
+        else:
+            n_edge = None
+
         for dim in set(ds.dims) ^ _source_dims_dict.keys():
             # obtain dimensions that were not parsed source_dims_dict and attempt to match to a grid element
             if ds.sizes[dim] == grid.n_face:
                 _source_dims_dict[dim] = "n_face"
             elif ds.sizes[dim] == grid.n_node:
                 _source_dims_dict[dim] = "n_node"
-            elif ds.sizes[dim] == grid.n_edge:
-                _source_dims_dict[dim] = "n_edge"
-
-        # Possible Issue: https://github.com/UXARRAY/uxarray/issues/610
+            elif n_edge is not None:
+                if ds.sizes[dim] == n_edge:
+                    _source_dims_dict[dim] = "n_edge"
 
         # rename dimensions to follow the UGRID conventions
         ds = ds.swap_dims(_source_dims_dict)
diff --git a/uxarray/grid/arcs.py b/uxarray/grid/arcs.py
index 3ccb4501b..f01418986 100644
--- a/uxarray/grid/arcs.py
+++ b/uxarray/grid/arcs.py
@@ -28,251 +28,68 @@ 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(
-    pt_xyz,
-    pt_lonlat,
-    gca_a_xyz,
-    gca_a_lonlat,
-    gca_b_xyz,
-    gca_b_lonlat,
-    is_directed=False,
-):
-    """Check if a point lies on a given Great Circle Arc (GCA). The anti-
-    meridian case is also considered.
+def point_within_gca(pt_xyz, gca_a_xyz, gca_b_xyz):
+    """
+    Check if a point lies on a given Great Circle Arc (GCA) interval, considering the smaller arc of the circle.
+    Handles the anti-meridian case as well.
 
     Parameters
     ----------
-    pt : numpy.ndarray (float)
+    pt_xyz : numpy.ndarray
         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
+    gca_a_xyz : numpy.ndarray
+        Cartesian coordinates of the first endpoint of the Great Circle Arc.
+    gca_b_xyz : numpy.ndarray
+        Cartesian coordinates of the second endpoint of the Great Circle Arc.
 
     Returns
     -------
     bool
-        True if the point lies between the two endpoints of the GCR, False otherwise.
+        True if the point lies within the specified GCA interval, 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.
+        If the input GCA spans exactly 180 degrees (π radians), as this GCA can have multiple planes.
+        In such cases, consider breaking the GCA into two separate arcs.
 
     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.
+    - The function ensures that the point lies on the same plane as the GCA before performing interval checks.
+    - It assumes the input represents the smaller arc of the Great Circle.
+    - The `_angle_of_2_vectors` and `_xyz_to_lonlat_rad_scalar` functions are used for calculations.
     """
-
-    # ==================================================================================================================
-    # 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):
+    # 1. Check if the input GCA spans exactly 180 degrees
+    angle_ab = _angle_of_2_vectors(gca_a_xyz, gca_b_xyz)
+    if np.allclose(angle_ab, 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"
+            "The input Great Circle Arc spans exactly 180 degrees, which can correspond to multiple planes. "
+            "Consider breaking the Great Circle Arc into two smaller 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
 
+    # 2. Verify if the point lies on the plane of the GCA
     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,
+        np.dot(cross_product, pt_xyz), 0, rtol=MACHINE_EPSILON, atol=MACHINE_EPSILON
     ):
         return False
-    # ==================================================================================================================
-    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 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(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 longitude arc, then the pt is not on the GCA
-            return False
-    # ==================================================================================================================
-    # 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 (
-        np.isclose(
-            np.abs(gca_b_lonlat[0] - gca_a_lonlat[0]),
-            np.pi,
-            rtol=0.0,
-            atol=MACHINE_EPSILON,
-        )
-        or np.isclose(
-            np.abs(gca_a_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 np.isclose(
-            np.abs(pt_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=ERROR_TOLERANCE
-        ):
-            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 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 np.isclose(
-                abs(gca_a_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 = gca_b_lonlat
-
-            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 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 (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 gca_a_lonlat[1] > 0.0 else -np.pi / 2
-                lat_extend = (
-                    abs(np.pi / 2 - abs(gca_a_lonlat[1]))
-                    + np.pi / 2
-                    + abs(gca_b_lonlat[1])
-                )
-            else:
-                pole_lat = _decide_pole_latitude(gca_a_lonlat[1], gca_b_lonlat[1])
-                lat_extend = (
-                    abs(np.pi / 2 - abs(gca_a_lonlat[1]))
-                    + np.pi / 2
-                    + abs(gca_b_lonlat[1])
-                )
-
-            if is_directed and lat_extend >= np.pi:
-                raise ValueError(
-                    "The input GCA spans more than 180 degrees. Please divide it into two."
-                )
-
-            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(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 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 > gca_a_lonlat[0] > np.pi > gca_b_lonlat[0] > 0.0:
-                return in_between(
-                    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(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([gca_a_lonlat[0], gca_b_lonlat[0]])
-        if np.pi > GCRv1_lonlat_max - GCRv0_lonlat_min >= 0.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
-            )
 
-    return None
+    # 3. Check if the point lies within the Great Circle Arc interval
+    pt_a = gca_a_xyz - pt_xyz
+    pt_b = gca_b_xyz - pt_xyz
+
+    # Use the dot product to determine the sign of the angle between pt_a and pt_b
+    cos_theta = np.dot(pt_a, pt_b)
+
+    # Return True if the point lies within the interval (smaller arc)
+    if cos_theta < 0:
+        return True
+    elif isclose(cos_theta, 0.0, atol=MACHINE_EPSILON):
+        # set error tolerance to 0.0
+        return True
+    else:
+        return False
 
 
 @njit(cache=True)
diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py
index 25e1d0101..5c11c8281 100644
--- a/uxarray/grid/coordinates.py
+++ b/uxarray/grid/coordinates.py
@@ -9,7 +9,7 @@
 from uxarray.constants import ERROR_TOLERANCE
 from typing import Union
 
-from uxarray.grid.utils import _angle_of_2_vectors
+from uxarray.grid.utils import _small_angle_of_2_vectors
 
 
 @njit(cache=True)
@@ -431,7 +431,7 @@ def _circle_from_two_points(p1, p2):
     v1 = np.array(_lonlat_rad_to_xyz(np.radians(p1[0]), np.radians(p1[1])))
     v2 = np.array(_lonlat_rad_to_xyz(np.radians(p2[0]), np.radians(p2[1])))
 
-    distance = _angle_of_2_vectors(v1, v2)
+    distance = _small_angle_of_2_vectors(v1, v2)
     radius = distance / 2
 
     return center, radius
@@ -467,9 +467,9 @@ def _circle_from_three_points(p1, p2, p3):
 
     radius = (
         max(
-            _angle_of_2_vectors(v1, v2),
-            _angle_of_2_vectors(v1, v3),
-            _angle_of_2_vectors(v2, v3),
+            _small_angle_of_2_vectors(v1, v2),
+            _small_angle_of_2_vectors(v1, v3),
+            _small_angle_of_2_vectors(v2, v3),
         )
         / 2
     )
@@ -496,7 +496,7 @@ def _is_inside_circle(circle, point):
     center, radius = circle
     v1 = np.array(_lonlat_rad_to_xyz(np.radians(center[0]), np.radians(center[1])))
     v2 = np.array(_lonlat_rad_to_xyz(np.radians(point[0]), np.radians(point[1])))
-    distance = _angle_of_2_vectors(v1, v2)
+    distance = _small_angle_of_2_vectors(v1, v2)
     return distance <= radius
 
 
diff --git a/uxarray/grid/geometry.py b/uxarray/grid/geometry.py
index 92a018011..f39120fed 100644
--- a/uxarray/grid/geometry.py
+++ b/uxarray/grid/geometry.py
@@ -740,9 +740,7 @@ def pole_point_inside_polygon(pole, face_edges_xyz, face_edges_lonlat):
         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
-        )
+        intersection_count = _check_intersection(ref_edge_xyz, face_edges_xyz)
         return (intersection_count % 2) != 0
 
     elif location == 0:  # Equator
@@ -811,17 +809,13 @@ def pole_point_inside_polygon(pole, face_edges_xyz, face_edges_lonlat):
         # 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
@@ -843,7 +837,7 @@ def _classify_polygon_location(face_edge_cart):
 
 
 @njit(cache=True)
-def _check_intersection(ref_edge_xyz, ref_edge_lonlat, edges_xyz, edges_lonlat):
+def _check_intersection(ref_edge_xyz, edges_xyz):
     """Check the number of intersections of the reference edge with the given edges.
 
     Parameters
@@ -872,12 +866,9 @@ def _check_intersection(ref_edge_xyz, ref_edge_lonlat, edges_xyz, edges_lonlat):
 
     for i in range(n_edges):
         edge_xyz = edges_xyz[i]
-        edge_lonlat = edges_lonlat[i]
 
         # compute intersection
-        intersection_point = gca_gca_intersection(
-            ref_edge_xyz, ref_edge_lonlat, edge_xyz, edge_lonlat
-        )
+        intersection_point = gca_gca_intersection(ref_edge_xyz, edge_xyz)
 
         if intersection_point.size != 0:
             if intersection_point.ndim == 1:
@@ -1197,12 +1188,8 @@ def _populate_face_latlon_bound(
             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(
diff --git a/uxarray/grid/integrate.py b/uxarray/grid/integrate.py
index 657026760..cb2d9d1d2 100644
--- a/uxarray/grid/integrate.py
+++ b/uxarray/grid/integrate.py
@@ -13,7 +13,6 @@ def _get_zonal_faces_weight_at_constLat(
     faces_edges_cart_candidate,
     latitude_cart,
     face_latlon_bound_candidate,
-    is_directed=False,
     is_latlonface=False,
     is_face_GCA_list=None,
 ):
@@ -37,10 +36,6 @@ def _get_zonal_faces_weight_at_constLat(
         Shape: (n_faces(candidate),,2, 2),
         [...,[lat_min, lat_max], [lon_min, lon_max],...]
 
-    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.
-
     is_latlonface : bool, optional, default=False
         A global flag to indicate if faces are latlon face. If True, then treat all faces as latlon faces. Latlon face means
         That all edge is either a longitude or constant latitude line. If False, then all edges are GCA.
@@ -104,7 +99,6 @@ def _get_zonal_faces_weight_at_constLat(
             face_edges,
             latitude_cart,
             face_latlon_bound_candidate[face_index],
-            is_directed=is_directed,
             is_latlonface=is_latlonface,
             is_GCA_list=is_GCA_list,
         )
@@ -188,7 +182,7 @@ def _is_edge_gca(is_GCA_list, is_latlonface, edges_z):
 
 
 def _get_faces_constLat_intersection_info(
-    face_edges_cart, latitude_cart, is_GCA_list, is_latlonface, is_directed
+    face_edges_cart, latitude_cart, is_GCA_list, is_latlonface
 ):
     """Processes each edge of a face polygon in a vectorized manner to
     determine overlaps and calculate the intersections for a given latitude and
@@ -206,9 +200,7 @@ def _get_faces_constLat_intersection_info(
     is_latlonface : bool
         Flag indicating if all faces are considered as lat-lon faces, meaning all edges are either
         constant latitude or longitude lines. This parameter overwrites the `is_GCA_list` if set to True.
-    is_directed : bool
-        Flag indicating if the GCA should be considered as directed (from v0 to v1). If False,
-        the smaller circle (less than 180 degrees) side of the GCA is used.
+
 
     Returns:
     -------
@@ -244,9 +236,7 @@ def _get_faces_constLat_intersection_info(
         # Calculate intersections (assuming a batch-capable intersection function)
         for idx, edge in enumerate(valid_edges):
             if is_GCA[idx]:
-                intersections = gca_const_lat_intersection(
-                    edge, latitude_cart, is_directed=is_directed
-                )
+                intersections = gca_const_lat_intersection(edge, latitude_cart)
 
                 if intersections.size == 0:
                     continue
@@ -299,7 +289,6 @@ def _get_zonal_face_interval(
     face_edges_cart,
     latitude_cart,
     face_latlon_bound,
-    is_directed=False,
     is_latlonface=False,
     is_GCA_list=None,
 ):
@@ -323,9 +312,6 @@ def _get_zonal_face_interval(
         The latitude in cartesian, the normalized Z coordinates.
     face_latlon_bound : np.ndarray
         The latitude and longitude bounds of the face. Shape: (2, 2), [[lat_min, lat_max], [lon_min, lon_max]]
-    is_directed : bool, optional
-        If True, the GCA is considered to be directed (from v0 to v1). If False, the GCA is undirected,
-        and the smaller circle (less than 180 degrees) side of the GCA is used. Default is False.
     is_latlonface : bool, optional, default=False
         A global flag to indicate if faces are latlon face. If True, then treat all faces as latlon faces. Latlon face means
         That all edge is either a longitude or constant latitude line. If False, then all edges are GCA.
@@ -347,7 +333,7 @@ def _get_zonal_face_interval(
         # Call the vectorized function to process all edges
         unique_intersections, pt_lon_min, pt_lon_max = (
             _get_faces_constLat_intersection_info(
-                face_edges_cart, latitude_cart, is_GCA_list, is_latlonface, is_directed
+                face_edges_cart, latitude_cart, is_GCA_list, is_latlonface
             )
         )
 
diff --git a/uxarray/grid/intersections.py b/uxarray/grid/intersections.py
index 528304d35..4550598bb 100644
--- a/uxarray/grid/intersections.py
+++ b/uxarray/grid/intersections.py
@@ -1,13 +1,14 @@
 import numpy as np
 from uxarray.constants import MACHINE_EPSILON, ERROR_TOLERANCE, INT_DTYPE
-from uxarray.grid.utils import _newton_raphson_solver_for_gca_constLat
+from uxarray.grid.utils import (
+    _newton_raphson_solver_for_gca_constLat,
+    _angle_of_2_vectors,
+)
 from uxarray.grid.arcs import (
-    point_within_gca,
     in_between,
     _extreme_gca_latitude_cartesian,
-    _point_within_gca_cartesian,
+    point_within_gca,
 )
-from uxarray.grid.coordinates import _xyz_to_lonlat_rad_scalar
 import platform
 import warnings
 from uxarray.utils.computing import cross_fma, allclose, cross, norm
@@ -128,7 +129,7 @@ def constant_lon_intersections_no_extreme(lon, edge_node_x, edge_node_y, n_edge)
     return np.unique(intersecting_edges)
 
 
-@njit(cache=True)
+# @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.
@@ -204,31 +205,11 @@ 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)
+    return gca_gca_intersection(gca_a_xyz, gca_b_xyz)
 
 
 @njit(cache=True)
-def gca_gca_intersection(gca_a_xyz, gca_a_lonlat, gca_b_xyz, gca_b_lonlat):
+def gca_gca_intersection(gca_a_xyz, gca_b_xyz):
     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")
 
@@ -238,12 +219,15 @@ def gca_gca_intersection(gca_a_xyz, gca_a_lonlat, gca_b_xyz, gca_b_lonlat):
     v0_xyz = gca_b_xyz[0]
     v1_xyz = gca_b_xyz[1]
 
-    w0_lonlat = gca_a_lonlat[0]
-    w1_lonlat = gca_a_lonlat[1]
-    v0_lonlat = gca_b_lonlat[0]
-    v1_lonlat = gca_b_lonlat[1]
+    angle_w0w1 = _angle_of_2_vectors(w0_xyz, w1_xyz)
+    angle_v0v1 = _angle_of_2_vectors(v0_xyz, v1_xyz)
+
+    if angle_w0w1 > np.pi:
+        w0_xyz, w1_xyz = w1_xyz, w0_xyz
+
+    if angle_v0v1 > np.pi:
+        v0_xyz, v1_xyz = v1_xyz, v0_xyz
 
-    # 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)
@@ -254,11 +238,11 @@ def gca_gca_intersection(gca_a_xyz, gca_a_lonlat, gca_b_xyz, gca_b_lonlat):
 
     # Check if the two GCAs are parallel
     if allclose(cross_norms, 0.0, atol=MACHINE_EPSILON):
-        if point_within_gca(v0_xyz, v0_lonlat, w0_xyz, w0_lonlat, w1_xyz, w1_lonlat):
+        if point_within_gca(v0_xyz, w0_xyz, w1_xyz):
             res[count, :] = v0_xyz
             count += 1
 
-        if point_within_gca(v1_xyz, v1_lonlat, w0_xyz, w0_lonlat, w1_xyz, w1_lonlat):
+        if point_within_gca(v1_xyz, w0_xyz, w1_xyz):
             res[count, :] = v1_xyz
             count += 1
 
@@ -269,29 +253,23 @@ def gca_gca_intersection(gca_a_xyz, gca_a_lonlat, gca_b_xyz, gca_b_lonlat):
     x1_xyz = cross_norms
     x2_xyz = -x1_xyz
 
-    # 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])
-
     # 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):
+    if point_within_gca(x1_xyz, w0_xyz, w1_xyz) and point_within_gca(
+        x1_xyz, v0_xyz, v1_xyz
+    ):
         res[count, :] = x1_xyz
         count += 1
 
-    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):
+    if point_within_gca(x2_xyz, w0_xyz, w1_xyz) and point_within_gca(
+        x2_xyz, v0_xyz, v1_xyz
+    ):
         res[count, :] = x2_xyz
         count += 1
 
     return res[:count, :]
 
 
-def gca_const_lat_intersection(
-    gca_cart, constZ, fma_disabled=True, verbose=False, is_directed=False
-):
+def gca_const_lat_intersection(gca_cart, constZ, fma_disabled=True, verbose=False):
     """Calculate the intersection point(s) of a Great Circle Arc (GCA) and a
     constant latitude line in a Cartesian coordinate system.
 
@@ -308,9 +286,6 @@ def gca_const_lat_intersection(
         If True, the FMA operation is disabled. And a naive `np.cross` is used instead.
     verbose : bool, optional (default=False)
         If True, the function prints out the intermediate results.
-    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.
 
     Returns
     -------
@@ -385,7 +360,7 @@ def gca_const_lat_intersection(
     res = None
 
     # Now test which intersection point is within the GCA range
-    if _point_within_gca_cartesian(p1, gca_cart, is_directed=is_directed):
+    if point_within_gca(p1, gca_cart[0], gca_cart[1]):
         try:
             converged_pt = _newton_raphson_solver_for_gca_constLat(
                 p1, gca_cart, verbose=verbose
@@ -407,7 +382,7 @@ def gca_const_lat_intersection(
         except RuntimeError:
             raise RuntimeError(f"Error encountered with initial guess: {p1}")
 
-    if _point_within_gca_cartesian(p2, gca_cart, is_directed=is_directed):
+    if point_within_gca(p2, gca_cart[0], gca_cart[1]):
         try:
             converged_pt = _newton_raphson_solver_for_gca_constLat(
                 p2, gca_cart, verbose=verbose
diff --git a/uxarray/grid/utils.py b/uxarray/grid/utils.py
index c88a9e51c..228fd96fb 100644
--- a/uxarray/grid/utils.py
+++ b/uxarray/grid/utils.py
@@ -7,21 +7,21 @@
 
 
 @njit(cache=True)
-def _angle_of_2_vectors(u, v):
-    """Calculate the angle between two 3D vectors u and v in radians. Can be
-    used to calcualte the span of a GCR.
+def _small_angle_of_2_vectors(u, v):
+    """
+    Compute the smallest angle between two vectors using the new _angle_of_2_vectors.
 
     Parameters
     ----------
-    u : numpy.ndarray (float)
+    u : numpy.ndarray
         The first 3D vector.
-    v : numpy.ndarray (float)
+    v : numpy.ndarray
         The second 3D vector.
 
     Returns
     -------
     float
-        The angle between u and v in radians.
+        The smallest angle between `u` and `v` in radians.
     """
     v_norm_times_u = np.linalg.norm(v) * u
     u_norm_times_v = np.linalg.norm(u) * v
@@ -31,6 +31,55 @@ def _angle_of_2_vectors(u, v):
     return angle_u_v_rad
 
 
+@njit(cache=True)
+def _angle_of_2_vectors(u, v):
+    """
+    Calculate the angle between two 3D vectors `u` and `v` on the unit sphere in radians.
+
+    This function computes the angle between two vectors originating from the center of a unit sphere.
+    The result is returned in the range [0, 2π]. It can be used to calculate the span of a great circle arc (GCA).
+
+    Parameters
+    ----------
+    u : numpy.ndarray
+        The first 3D vector (float), originating from the center of the unit sphere.
+    v : numpy.ndarray
+        The second 3D vector (float), originating from the center of the unit sphere.
+
+    Returns
+    -------
+    float
+        The angle between `u` and `v` in radians, in the range [0, 2π].
+
+    Notes
+    -----
+    - The direction of the angle (clockwise or counter-clockwise) is determined using the cross product of `u` and `v`.
+    - Special cases such as vectors aligned along the same longitude are handled explicitly.
+    """
+    # Compute the cross product to determine the direction of the normal
+    normal = np.cross(u, v)
+
+    # Calculate the angle using arctangent of cross and dot products
+    angle_u_v_rad = np.arctan2(np.linalg.norm(normal), np.dot(u, v))
+
+    # Determine the direction of the angle
+    normal_z = np.dot(normal, np.array([0.0, 0.0, 1.0]))
+    if normal_z > 0:
+        # Counterclockwise direction
+        return angle_u_v_rad
+    elif normal_z == 0:
+        # Handle collinear vectors (same longitude)
+        if u[2] > v[2]:
+            return angle_u_v_rad
+        elif u[2] < v[2]:
+            return 2 * np.pi - angle_u_v_rad
+        else:
+            return 0.0  # u == v
+    else:
+        # Clockwise direction
+        return 2 * np.pi - angle_u_v_rad
+
+
 def _inv_jacobian(x0, x1, y0, y1, z0, z1, x_i_old, y_i_old):
     """Calculate the inverse Jacobian matrix for a given set of parameters.
 
diff --git a/uxarray/io/_scrip.py b/uxarray/io/_scrip.py
index c07824147..f6b874e93 100644
--- a/uxarray/io/_scrip.py
+++ b/uxarray/io/_scrip.py
@@ -1,6 +1,8 @@
 import xarray as xr
 import numpy as np
 
+import polars as pl
+
 from uxarray.grid.connectivity import _replace_fill_values
 from uxarray.constants import INT_DTYPE, INT_FILL_VALUE
 
@@ -11,43 +13,45 @@ def _to_ugrid(in_ds, out_ds):
     """If input dataset (``in_ds``) file is an unstructured SCRIP file,
     function will reassign SCRIP variables to UGRID conventions in output file
     (``out_ds``).
-
-    Parameters
-    ----------
-    in_ds : xarray.Dataset
-        Original scrip dataset of interest being used
-
-    out_ds : xarray.Variable
-        file to be returned by ``_populate_scrip_data``, used as an empty placeholder file
-        to store reassigned SCRIP variables in UGRID conventions
     """
 
     source_dims_dict = {}
 
     if in_ds["grid_area"].all():
         # Create node_lon & node_lat variables from grid_corner_lat/lon
-        # Turn latitude scrip array into 1D instead of 2D
+        # Turn latitude and longitude scrip arrays into 1D
         corner_lat = in_ds["grid_corner_lat"].values.ravel()
-
-        # Repeat above steps with longitude data instead
         corner_lon = in_ds["grid_corner_lon"].values.ravel()
 
-        # Combine flat lat and lon arrays
-        corner_lon_lat = np.vstack((corner_lon, corner_lat)).T
+        # Use Polars to find unique coordinate pairs
+        df = pl.DataFrame({"lon": corner_lon, "lat": corner_lat}).with_row_count(
+            "original_index"
+        )
+
+        # Get unique rows (first occurrence). This preserves the order in which they appear.
+        unique_df = df.unique(subset=["lon", "lat"], keep="first")
+
+        # unq_ind: The indices of the unique rows in the original array
+        unq_ind = unique_df["original_index"].to_numpy().astype(INT_DTYPE)
+
+        # To get the inverse index (unq_inv): map each original row back to its unique row index.
+        # Add a unique_id to the unique_df which will serve as the "inverse" mapping.
+        unique_df = unique_df.with_row_count("unique_id")
 
-        # Run numpy unique to determine which rows/values are actually unique
-        _, unq_ind, unq_inv = np.unique(
-            corner_lon_lat, return_index=True, return_inverse=True, axis=0
+        # Join original df with unique_df to find out which unique_id corresponds to each row
+        df_joined = df.join(
+            unique_df.drop("original_index"), on=["lon", "lat"], how="left"
         )
+        unq_inv = df_joined["unique_id"].to_numpy().astype(INT_DTYPE)
 
-        # Now, calculate unique lon and lat values to account for 'node_lon' and 'node_lat'
-        unq_lon = corner_lon_lat[unq_ind, :][:, 0]
-        unq_lat = corner_lon_lat[unq_ind, :][:, 1]
+        # Extract unique lon and lat values using unq_ind
+        unq_lon = corner_lon[unq_ind]
+        unq_lat = corner_lat[unq_ind]
 
         # Reshape face nodes array into original shape for use in 'face_node_connectivity'
         unq_inv = np.reshape(unq_inv, (len(in_ds.grid_size), len(in_ds.grid_corners)))
 
-        # Create node_lon & node_lat from unsorted, unique grid_corner_lat/lon
+        # Create node_lon & node_lat
         out_ds[ugrid.NODE_COORDINATES[0]] = xr.DataArray(
             unq_lon, dims=[ugrid.NODE_DIM], attrs=ugrid.NODE_LON_ATTRS
         )