diff --git a/docs/internal_api/index.rst b/docs/internal_api/index.rst
index fc90d681d..b233e2413 100644
--- a/docs/internal_api/index.rst
+++ b/docs/internal_api/index.rst
@@ -147,15 +147,17 @@ Coordinates
 .. autosummary::
    :toctree: generated/
 
-   grid.coordinates._get_lonlat_from_xyz
-   grid.coordinates._get_xyz_from_lonlat
-   grid.coordinates._populate_cartesian_xyz_coord
-   grid.coordinates._populate_lonlat_coord
+   grid.coordinates._lonlat_rad_to_xyz
+   grid.coordinates._xyz_to_lonlat_rad
+   grid.coordinates._xyz_to_lonlat_deg
+   grid.coordinates._normalize_xyz
+   grid.coordinates._populate_node_latlon
+   grid.coordinates._populate_node_xyz
    grid.coordinates._populate_face_centroids
-   grid.coordinates._construct_face_centroids
-   grid.coordinates._set_desired_longitude_range
    grid.coordinates._populate_edge_centroids
+   grid.coordinates._construct_face_centroids
    grid.coordinates._construct_edge_centroids
+   grid.coordinates._set_desired_longitude_range
 
 
 Arcs
diff --git a/docs/user_api/index.rst b/docs/user_api/index.rst
index 7f67456ff..ad51f3315 100644
--- a/docs/user_api/index.rst
+++ b/docs/user_api/index.rst
@@ -369,16 +369,6 @@ Connectivity
 
    grid.connectivity.close_face_nodes
 
-Coordinates
------------
-.. autosummary::
-   :toctree: generated/
-
-   grid.coordinates.node_lonlat_rad_to_xyz
-   grid.coordinates.node_xyz_to_lonlat_rad
-   grid.coordinates.normalize_in_place
-
-
 Arcs
 ----
 .. autosummary::
diff --git a/test/test_arcs.py b/test/test_arcs.py
index c2e24f5aa..c9859e724 100644
--- a/test/test_arcs.py
+++ b/test/test_arcs.py
@@ -9,7 +9,7 @@
 
 import uxarray as ux
 
-from uxarray.grid.coordinates import node_lonlat_rad_to_xyz
+from uxarray.grid.coordinates import _lonlat_rad_to_xyz
 from uxarray.grid.arcs import point_within_gca
 
 try:
@@ -31,43 +31,40 @@ class TestIntersectionPoint(TestCase):
     def test_pt_within_gcr(self):
         # The GCR that's eexactly 180 degrees will have Value Error raised
         gcr_180degree_cart = [
-            ux.grid.coordinates.node_lonlat_rad_to_xyz([0.0, 0.0]),
-            ux.grid.coordinates.node_lonlat_rad_to_xyz([np.pi, 0.0])
+            _lonlat_rad_to_xyz(0.0, 0.0),
+            _lonlat_rad_to_xyz(np.pi, 0.0)
         ]
-        pt_same_lon_in = ux.grid.coordinates.node_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(pt_same_lon_in, gcr_180degree_cart)
 
         gcr_180degree_cart = [
-            ux.grid.coordinates.node_lonlat_rad_to_xyz([0.0, np.pi / 2.0]),
-            ux.grid.coordinates.node_lonlat_rad_to_xyz([0.0, -np.pi / 2.0])
+            _lonlat_rad_to_xyz(0.0, np.pi / 2.0),
+            _lonlat_rad_to_xyz(0.0, -np.pi / 2.0)
         ]
 
-        pt_same_lon_in = ux.grid.coordinates.node_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(pt_same_lon_in, gcr_180degree_cart)
 
         # Test when the point and the GCA all have the same longitude
         gcr_same_lon_cart = [
-            ux.grid.coordinates.node_lonlat_rad_to_xyz([0.0, 1.5]),
-            ux.grid.coordinates.node_lonlat_rad_to_xyz([0.0, -1.5])
+            _lonlat_rad_to_xyz(0.0, 1.5),
+            _lonlat_rad_to_xyz(0.0, -1.5)
         ]
-        pt_same_lon_in = ux.grid.coordinates.node_lonlat_rad_to_xyz([0.0, 0.0])
+        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))
 
-        pt_same_lon_out = ux.grid.coordinates.node_lonlat_rad_to_xyz(
-            [0.0, 1.500000000000001])
+        pt_same_lon_out = _lonlat_rad_to_xyz(0.0, 1.500000000000001)
         res = point_within_gca(pt_same_lon_out, gcr_same_lon_cart)
         self.assertFalse(res)
 
-        pt_same_lon_out_2 = ux.grid.coordinates.node_lonlat_rad_to_xyz(
-            [0.1, 1.0])
+        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)
         self.assertFalse(res)
 
         # And if we increase the digital place by one, it should be true again
-        pt_same_lon_out_add_one_place = ux.grid.coordinates.node_lonlat_rad_to_xyz(
-            [0.0, 1.5000000000000001])
+        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)
         self.assertTrue(res)
 
@@ -117,10 +114,10 @@ def test_pt_within_gcr_antimeridian(self):
         # 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 = ux.grid.coordinates.node_lonlat_rad_to_xyz(v1_rad)
-        v2_cart = ux.grid.coordinates.node_lonlat_rad_to_xyz(v2_rad)
+        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 = ux.grid.coordinates.node_lonlat_rad_to_xyz([0.01, 0.0])
+        pt_cart = _lonlat_rad_to_xyz(0.01, 0.0)
         with self.assertRaises(ValueError):
             point_within_gca(pt_cart, gcr_cart, is_directed=True)
         gcr_car_flipped = np.array([v2_cart, v1_cart])
diff --git a/test/test_centroids.py b/test/test_centroids.py
index 4b5de9fd4..6eb69542f 100644
--- a/test/test_centroids.py
+++ b/test/test_centroids.py
@@ -4,7 +4,7 @@
 import numpy.testing as nt
 import uxarray as ux
 from pathlib import Path
-from uxarray.grid.coordinates import _populate_face_centroids, _populate_edge_centroids, normalize_in_place
+from uxarray.grid.coordinates import _populate_face_centroids, _populate_edge_centroids, _normalize_xyz
 
 current_path = Path(os.path.dirname(os.path.realpath(__file__)))
 
@@ -21,8 +21,8 @@ def test_centroids_from_mean_verts_triangle(self):
 
         # Calculate the expected centroid
         expected_centroid = np.mean(test_triangle, axis=0)
-        [norm_x, norm_y, norm_z] = normalize_in_place(
-            [expected_centroid[0], expected_centroid[1], expected_centroid[2]])
+        norm_x, norm_y, norm_z = _normalize_xyz(
+            expected_centroid[0], expected_centroid[1], expected_centroid[2])
 
         # Open the dataset and find the centroids
         grid = ux.open_grid(test_triangle)
@@ -42,8 +42,8 @@ def test_centroids_from_mean_verts_pentagon(self):
 
         # Calculate the expected centroid
         expected_centroid = np.mean(test_triangle, axis=0)
-        [norm_x, norm_y, norm_z] = normalize_in_place(
-            [expected_centroid[0], expected_centroid[1], expected_centroid[2]])
+        norm_x, norm_y, norm_z = _normalize_xyz(
+            expected_centroid[0], expected_centroid[1], expected_centroid[2])
 
         # Open the dataset and find the centroids
         grid = ux.open_grid(test_triangle)
@@ -65,7 +65,8 @@ def test_centroids_from_mean_verts_scrip(self):
 
         _populate_face_centroids(uxgrid, repopulate=True)
 
-        computed_face_x = (uxgrid.face_lon.values + 180) % 360 - 180
+        # computed_face_x = (uxgrid.face_lon.values + 180) % 360 - 180
+        computed_face_x = uxgrid.face_lon.values
         computed_face_y = uxgrid.face_lat.values
 
         nt.assert_array_almost_equal(expected_face_x, computed_face_x)
diff --git a/test/test_computing.py b/test/test_computing.py
index 4e9cde81b..a83d3b76c 100644
--- a/test/test_computing.py
+++ b/test/test_computing.py
@@ -1,7 +1,7 @@
 from unittest import TestCase
 import numpy.testing as nt
 import numpy as np
-from uxarray.grid.coordinates import normalize_in_place
+from uxarray.grid.coordinates import _normalize_xyz
 import uxarray.utils.computing as ac_utils
 from uxarray.constants import ERROR_TOLERANCE
 
@@ -12,8 +12,8 @@ class TestCrossProduct(TestCase):
     one."""
 
     def test_cross_fma(self):
-        v1 = np.array(normalize_in_place([1.0, 2.0, 3.0]))
-        v2 = np.array(normalize_in_place([4.0, 5.0, 6.0]))
+        v1 = np.array(_normalize_xyz(*[1.0, 2.0, 3.0]))
+        v2 = np.array(_normalize_xyz(*[4.0, 5.0, 6.0]))
 
         np_cross = np.cross(v1, v2)
         fma_cross = ac_utils.cross_fma(v1, v2)
@@ -26,8 +26,8 @@ class TestDotProduct(TestCase):
     one."""
 
     def test_dot_fma(self):
-        v1 = np.array(normalize_in_place([1.0, 0.0, 0.0]), dtype=np.float64)
-        v2 = np.array(normalize_in_place([1.0, 0.0, 0.0]), dtype=np.float64)
+        v1 = np.array(_normalize_xyz(*[1.0, 0.0, 0.0]), dtype=np.float64)
+        v2 = np.array(_normalize_xyz(*[1.0, 0.0, 0.0]), dtype=np.float64)
 
         np_dot = np.dot(v1, v2)
         fma_dot = ac_utils.dot_fma(v1, v2)
diff --git a/test/test_geometry.py b/test/test_geometry.py
index beaa105b7..b73e2d8b9 100644
--- a/test/test_geometry.py
+++ b/test/test_geometry.py
@@ -9,7 +9,7 @@
 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_lonlat_coord, node_lonlat_rad_to_xyz
+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.utils import _get_cartesian_face_edge_nodes, _get_lonlat_rad_face_edge_nodes
 from uxarray.grid.geometry import _populate_face_latlon_bound, _populate_bounds
@@ -67,7 +67,7 @@ def test_pole_point_inside_polygon_from_vertice_north(self):
         # Normalize the vertices to ensure they lie on the unit sphere
         for i, vertex in enumerate(vertices):
             float_vertex = [float(coord) for coord in vertex]
-            vertices[i] = ux.grid.coordinates.normalize_in_place(float_vertex)
+            vertices[i] = _normalize_xyz(*float_vertex)
 
         # Create face_edge_cart from the vertices
         face_edge_cart = np.array([[vertices[0], vertices[1]],
@@ -93,7 +93,7 @@ def test_pole_point_inside_polygon_from_vertice_south(self):
         # Normalize the vertices to ensure they lie on the unit sphere
         for i, vertex in enumerate(vertices):
             float_vertex = [float(coord) for coord in vertex]
-            vertices[i] = ux.grid.coordinates.normalize_in_place(float_vertex)
+            vertices[i] = _normalize_xyz(*float_vertex)
 
         # Create face_edge_cart from the vertices, since we are using the south pole, and want retrive the smaller face
         # we need to reverse the order of the vertices
@@ -121,7 +121,7 @@ def test_pole_point_inside_polygon_from_vertice_pole(self):
         # Normalize the vertices to ensure they lie on the unit sphere
         for i, vertex in enumerate(vertices):
             float_vertex = [float(coord) for coord in vertex]
-            vertices[i] = ux.grid.coordinates.normalize_in_place(float_vertex)
+            vertices[i] = _normalize_xyz(*float_vertex)
 
         # Create face_edge_cart from the vertices
         face_edge_cart = np.array([[vertices[0], vertices[1]],
@@ -147,7 +147,7 @@ def test_pole_point_inside_polygon_from_vertice_cross(self):
         # Normalize the vertices to ensure they lie on the unit sphere
         for i, vertex in enumerate(vertices):
             float_vertex = [float(coord) for coord in vertex]
-            vertices[i] = ux.grid.coordinates.normalize_in_place(float_vertex)
+            vertices[i] = _normalize_xyz(*float_vertex)
 
         # Create face_edge_cart from the vertices
         face_edge_cart = np.array([[vertices[0], vertices[1]],
@@ -191,13 +191,13 @@ def _max_latitude_rad_iterative(self, gca_cart):
         # Convert input vectors to radians and Cartesian coordinates
 
         v1_cart, v2_cart = gca_cart
-        b_lonlat = ux.grid.coordinates.node_xyz_to_lonlat_rad(v1_cart.tolist())
-        c_lonlat = ux.grid.coordinates.node_xyz_to_lonlat_rad(v2_cart.tolist())
+        b_lonlat = _xyz_to_lonlat_rad(*v1_cart.tolist())
+        c_lonlat = _xyz_to_lonlat_rad(*v2_cart.tolist())
 
         # Initialize variables for the iterative process
         v_temp = ac_utils.cross_fma(v1_cart, v2_cart)
         v0 = ac_utils.cross_fma(v_temp, v1_cart)
-        v0 = ux.grid.coordinates.normalize_in_place(v0.tolist())
+        v0 = _normalize_xyz(*v0.tolist())
         max_section = [v1_cart, v2_cart]
 
         # Iteratively find the maximum latitude
@@ -207,7 +207,7 @@ def _max_latitude_rad_iterative(self, gca_cart):
             v_b, v_c = max_section
             angle_v1_v2_rad = ux.grid.arcs._angle_of_2_vectors(v_b, v_c)
             v0 = ac_utils.cross_fma(v_temp, v_b)
-            v0 = ux.grid.coordinates.normalize_in_place(v0.tolist())
+            v0 = _normalize_xyz(*v0.tolist())
             avg_angle_rad = angle_v1_v2_rad / 10.0
 
             for i in range(10):
@@ -217,10 +217,13 @@ def _max_latitude_rad_iterative(self, gca_cart):
                     angle_rad_prev) * np.array(v0)
                 w2_new = np.cos(angle_rad_next) * v_b + np.sin(
                     angle_rad_next) * np.array(v0)
-                w1_lonlat = ux.grid.coordinates.node_xyz_to_lonlat_rad(
-                    w1_new.tolist())
-                w2_lonlat = ux.grid.coordinates.node_xyz_to_lonlat_rad(
-                    w2_new.tolist())
+                w1_lonlat = _xyz_to_lonlat_rad(
+                    *w1_new.tolist())
+                w2_lonlat = _xyz_to_lonlat_rad(
+                    *w2_new.tolist())
+
+                w1_lonlat = np.asarray(w1_lonlat)
+                w2_lonlat = np.asarray(w2_lonlat)
 
                 # Adjust latitude boundaries to avoid error accumulation
                 if i == 0:
@@ -232,7 +235,7 @@ def _max_latitude_rad_iterative(self, gca_cart):
                 max_lat = max(max_lat, w1_lonlat[1], w2_lonlat[1])
                 if np.abs(w2_lonlat[1] -
                           w1_lonlat[1]) <= ERROR_TOLERANCE or w1_lonlat[
-                              1] == max_lat == w2_lonlat[1]:
+                    1] == max_lat == w2_lonlat[1]:
                     max_section = [w1_new, w2_new]
                     break
                 if np.abs(max_lat - w1_lonlat[1]) <= ERROR_TOLERANCE:
@@ -241,10 +244,10 @@ def _max_latitude_rad_iterative(self, gca_cart):
                     max_section = [w1_new, w2_new] if i != 9 else [w1_new, v_c]
 
             # Update longitude and latitude for the next iteration
-            b_lonlat = ux.grid.coordinates.node_xyz_to_lonlat_rad(
-                max_section[0].tolist())
-            c_lonlat = ux.grid.coordinates.node_xyz_to_lonlat_rad(
-                max_section[1].tolist())
+            b_lonlat = _xyz_to_lonlat_rad(
+                *max_section[0].tolist())
+            c_lonlat = _xyz_to_lonlat_rad(
+                *max_section[1].tolist())
 
         return np.average([b_lonlat[1], c_lonlat[1]])
 
@@ -275,13 +278,13 @@ def _min_latitude_rad_iterative(self, gca_cart):
 
         # Convert input vectors to radians and Cartesian coordinates
         v1_cart, v2_cart = gca_cart
-        b_lonlat = ux.grid.coordinates.node_xyz_to_lonlat_rad(v1_cart.tolist())
-        c_lonlat = ux.grid.coordinates.node_xyz_to_lonlat_rad(v2_cart.tolist())
+        b_lonlat = _xyz_to_lonlat_rad(*v1_cart.tolist())
+        c_lonlat = _xyz_to_lonlat_rad(*v2_cart.tolist())
 
         # Initialize variables for the iterative process
         v_temp = ac_utils.cross_fma(v1_cart, v2_cart)
         v0 = ac_utils.cross_fma(v_temp, v1_cart)
-        v0 = np.array(ux.grid.coordinates.normalize_in_place(v0.tolist()))
+        v0 = np.array(_normalize_xyz(*v0.tolist()))
         min_section = [v1_cart, v2_cart]
 
         # Iteratively find the minimum latitude
@@ -291,7 +294,7 @@ def _min_latitude_rad_iterative(self, gca_cart):
             v_b, v_c = min_section
             angle_v1_v2_rad = ux.grid.arcs._angle_of_2_vectors(v_b, v_c)
             v0 = ac_utils.cross_fma(v_temp, v_b)
-            v0 = np.array(ux.grid.coordinates.normalize_in_place(v0.tolist()))
+            v0 = np.array(_normalize_xyz(*v0.tolist()))
             avg_angle_rad = angle_v1_v2_rad / 10.0
 
             for i in range(10):
@@ -301,10 +304,13 @@ def _min_latitude_rad_iterative(self, gca_cart):
                     angle_rad_prev) * v0
                 w2_new = np.cos(angle_rad_next) * v_b + np.sin(
                     angle_rad_next) * v0
-                w1_lonlat = ux.grid.coordinates.node_xyz_to_lonlat_rad(
-                    w1_new.tolist())
-                w2_lonlat = ux.grid.coordinates.node_xyz_to_lonlat_rad(
-                    w2_new.tolist())
+                w1_lonlat = _xyz_to_lonlat_rad(
+                    *w1_new.tolist())
+                w2_lonlat = _xyz_to_lonlat_rad(
+                    *w2_new.tolist())
+
+                w1_lonlat = np.asarray(w1_lonlat)
+                w2_lonlat = np.asarray(w2_lonlat)
 
                 # Adjust latitude boundaries to avoid error accumulation
                 if i == 0:
@@ -316,7 +322,7 @@ def _min_latitude_rad_iterative(self, gca_cart):
                 min_lat = min(min_lat, w1_lonlat[1], w2_lonlat[1])
                 if np.abs(w2_lonlat[1] -
                           w1_lonlat[1]) <= ERROR_TOLERANCE or w1_lonlat[
-                              1] == min_lat == w2_lonlat[1]:
+                    1] == min_lat == w2_lonlat[1]:
                     min_section = [w1_new, w2_new]
                     break
                 if np.abs(min_lat - w1_lonlat[1]) <= ERROR_TOLERANCE:
@@ -325,18 +331,18 @@ def _min_latitude_rad_iterative(self, gca_cart):
                     min_section = [w1_new, w2_new] if i != 9 else [w1_new, v_c]
 
             # Update longitude and latitude for the next iteration
-            b_lonlat = ux.grid.coordinates.node_xyz_to_lonlat_rad(
-                min_section[0].tolist())
-            c_lonlat = ux.grid.coordinates.node_xyz_to_lonlat_rad(
-                min_section[1].tolist())
+            b_lonlat = _xyz_to_lonlat_rad(
+                *min_section[0].tolist())
+            c_lonlat = _xyz_to_lonlat_rad(
+                *min_section[1].tolist())
 
         return np.average([b_lonlat[1], c_lonlat[1]])
 
     def test_extreme_gca_latitude_max(self):
         # Define a great circle arc that is symmetrical around 0 degrees longitude
         gca_cart = np.array([
-            ux.grid.coordinates.normalize_in_place([0.5, 0.5, 0.5]),
-            ux.grid.coordinates.normalize_in_place([-0.5, 0.5, 0.5])
+            _normalize_xyz(*[0.5, 0.5, 0.5]),
+            _normalize_xyz(*[-0.5, 0.5, 0.5])
         ])
 
         # Calculate the maximum latitude
@@ -363,8 +369,8 @@ def test_extreme_gca_latitude_max(self):
     def test_extreme_gca_latitude_min(self):
         # Define a great circle arc that is symmetrical around 0 degrees longitude
         gca_cart = np.array([
-            ux.grid.coordinates.normalize_in_place([0.5, 0.5, -0.5]),
-            ux.grid.coordinates.normalize_in_place([-0.5, 0.5, -0.5])
+            _normalize_xyz(*[0.5, 0.5, -0.5]),
+            _normalize_xyz(*[-0.5, 0.5, -0.5])
         ])
 
         # Calculate the minimum latitude
@@ -432,7 +438,6 @@ def test_insert_pt_in_empty_state(self):
         np.testing.assert_array_equal(result, expected)
 
 
-
 class TestLatlonBoundsGCA(TestCase):
     def test_populate_bounds_normal(self):
         # Generate a normal face that is not crossing the antimeridian or the poles
@@ -441,7 +446,9 @@ def test_populate_bounds_normal(self):
 
         # Convert everything into radians
         vertices_rad = np.radians(vertices_lonlat)
-        vertices_cart = [node_lonlat_rad_to_xyz(v) for v in vertices_rad]
+
+        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"))
         lat_min = min(np.deg2rad(10.0),
@@ -470,7 +477,7 @@ def test_populate_bounds_antimeridian(self):
 
         # Convert everything into radians
         vertices_rad = np.radians(vertices_lonlat)
-        vertices_cart = [node_lonlat_rad_to_xyz(v) for v in vertices_rad]
+        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"))
         lat_min = min(np.deg2rad(10.0),
@@ -499,7 +506,7 @@ def test_populate_bounds_node_on_pole(self):
 
         # Convert everything into radians
         vertices_rad = np.radians(vertices_lonlat)
-        vertices_cart = [node_lonlat_rad_to_xyz(v) for v in vertices_rad]
+        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"))
@@ -527,7 +534,7 @@ def test_populate_bounds_edge_over_pole(self):
 
         # Convert everything into radians
         vertices_rad = np.radians(vertices_lonlat)
-        vertices_cart = [node_lonlat_rad_to_xyz(v) for v in vertices_rad]
+        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"))
@@ -555,7 +562,7 @@ def test_populate_bounds_pole_inside(self):
 
         # Convert everything into radians
         vertices_rad = np.radians(vertices_lonlat)
-        vertices_cart = [node_lonlat_rad_to_xyz(v) for v in vertices_rad]
+        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"))
@@ -576,6 +583,7 @@ def test_populate_bounds_pole_inside(self):
         bounds = _populate_face_latlon_bound(face_edges_connectivity_cartesian, face_edges_connectivity_lonlat)
         nt.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE)
 
+
 class TestLatlonBoundsLatLonFace(TestCase):
 
     def test_populate_bounds_normal(self):
@@ -666,7 +674,7 @@ def test_populate_bounds_edge_over_pole(self):
 
         # Convert everything into radians
         vertices_rad = np.radians(vertices_lonlat)
-        vertices_cart = [node_lonlat_rad_to_xyz(v) for v in vertices_rad]
+        vertices_cart = np.vstack([_lonlat_rad_to_xyz(vertices_rad[:, 0], vertices_rad[:, 1])]).T
         lat_max = np.pi / 2
         lat_min = np.deg2rad(60.0)
         lon_min = np.deg2rad(210.0)
@@ -694,7 +702,7 @@ def test_populate_bounds_pole_inside(self):
 
         # Convert everything into radians
         vertices_rad = np.radians(vertices_lonlat)
-        vertices_cart = [node_lonlat_rad_to_xyz(v) for v in vertices_rad]
+        vertices_cart = np.vstack([_lonlat_rad_to_xyz(vertices_rad[:, 0], vertices_rad[:, 1])]).T
         lat_max = np.pi / 2
         lat_min = np.deg2rad(60.0)
         lon_min = 0
@@ -752,7 +760,7 @@ def test_populate_bounds_antimeridian(self):
 
         # Convert everything into radians
         vertices_rad = np.radians(vertices_lonlat)
-        vertices_cart = [node_lonlat_rad_to_xyz(v) for v in vertices_rad]
+        vertices_cart = np.vstack([_lonlat_rad_to_xyz(vertices_rad[:, 0], vertices_rad[:, 1])]).T
         lat_max = np.deg2rad(60.0)
         lat_min = np.deg2rad(10.0)
         lon_min = np.deg2rad(350.0)
@@ -780,7 +788,7 @@ def test_populate_bounds_node_on_pole(self):
 
         # Convert everything into radians
         vertices_rad = np.radians(vertices_lonlat)
-        vertices_cart = [node_lonlat_rad_to_xyz(v) for v in vertices_rad]
+        vertices_cart = np.vstack([_lonlat_rad_to_xyz(vertices_rad[:, 0], vertices_rad[:, 1])]).T
         lat_max = np.pi / 2
         lat_min = np.deg2rad(10.0)
         lon_min = np.deg2rad(10.0)
@@ -808,7 +816,7 @@ def test_populate_bounds_edge_over_pole(self):
 
         # Convert everything into radians
         vertices_rad = np.radians(vertices_lonlat)
-        vertices_cart = [node_lonlat_rad_to_xyz(v) for v in vertices_rad]
+        vertices_cart = np.vstack([_lonlat_rad_to_xyz(vertices_rad[:, 0], vertices_rad[:, 1])]).T
         lat_max = np.pi / 2
         lat_min = np.deg2rad(60.0)
         lon_min = np.deg2rad(210.0)
@@ -836,7 +844,7 @@ def test_populate_bounds_pole_inside(self):
 
         # Convert everything into radians
         vertices_rad = np.radians(vertices_lonlat)
-        vertices_cart = [node_lonlat_rad_to_xyz(v) for v in vertices_rad]
+        vertices_cart = np.vstack([_lonlat_rad_to_xyz(vertices_rad[:, 0], vertices_rad[:, 1])]).T
         lat_max = np.pi / 2
         lat_min = np.deg2rad(60.0)
         lon_min = 0
@@ -857,6 +865,7 @@ def test_populate_bounds_pole_inside(self):
                                              is_GCA_list=[True, False, True, False])
         nt.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE)
 
+
 class TestLatlonBoundsMix(TestCase):
     def test_populate_bounds_GCA_mix(self):
         face_1 = [[10.0, 60.0], [10.0, 10.0], [50.0, 10.0], [50.0, 60.0]]
@@ -867,11 +876,10 @@ def test_populate_bounds_GCA_mix(self):
         faces = [face_1, face_2, face_3, face_4]
 
         # Hand calculated bounds for the above faces in radians
-        expected_bounds = [[[0.17453293, 1.07370494],[0.17453293, 0.87266463]],
-                           [[0.17453293, 1.10714872],[6.10865238, 0.87266463]],
-                           [[1.04719755, 1.57079633],[3.66519143, 0.52359878]],
-                           [[1.04719755,1.57079633],[0.,         6.28318531]]]
-
+        expected_bounds = [[[0.17453293, 1.07370494], [0.17453293, 0.87266463]],
+                           [[0.17453293, 1.10714872], [6.10865238, 0.87266463]],
+                           [[1.04719755, 1.57079633], [3.66519143, 0.52359878]],
+                           [[1.04719755, 1.57079633], [0., 6.28318531]]]
 
         grid = ux.Grid.from_face_vertices(faces, latlon=True)
         face_bounds = grid.bounds.values
@@ -886,14 +894,13 @@ def test_populate_bounds_LatlonFace_mix(self):
 
         faces = [face_1, face_2, face_3, face_4]
 
-        expected_bounds = [[[np.deg2rad(10.0), np.deg2rad(60.0)],[np.deg2rad(10.0), np.deg2rad(50.0)]],
-                           [[np.deg2rad(10.0), np.deg2rad(60.0)],[np.deg2rad(350.0), np.deg2rad(50.0)]],
-                           [[np.deg2rad(60.0), np.pi/2],[np.deg2rad(210.0), np.deg2rad(30.0)]],
-                           [[np.deg2rad(60.0),np.pi/2],[0.,        2*np.pi]]]
-
+        expected_bounds = [[[np.deg2rad(10.0), np.deg2rad(60.0)], [np.deg2rad(10.0), np.deg2rad(50.0)]],
+                           [[np.deg2rad(10.0), np.deg2rad(60.0)], [np.deg2rad(350.0), np.deg2rad(50.0)]],
+                           [[np.deg2rad(60.0), np.pi / 2], [np.deg2rad(210.0), np.deg2rad(30.0)]],
+                           [[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_latlonface=True, return_array=True)
+        bounds_xarray = _populate_bounds(grid, is_latlonface=True, 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)
@@ -906,16 +913,18 @@ def test_populate_bounds_GCAList_mix(self):
 
         faces = [face_1, face_2, face_3, face_4]
 
-        expected_bounds = [[[np.deg2rad(10.0), np.deg2rad(60.0)],[np.deg2rad(10.0), np.deg2rad(50.0)]],
-                           [[np.deg2rad(10.0), np.deg2rad(60.0)],[np.deg2rad(350.0), np.deg2rad(50.0)]],
-                           [[np.deg2rad(60.0), np.pi/2],[np.deg2rad(210.0), np.deg2rad(30.0)]],
-                           [[np.deg2rad(60.0),np.pi/2],[0.,        2*np.pi]]]
+        expected_bounds = [[[np.deg2rad(10.0), np.deg2rad(60.0)], [np.deg2rad(10.0), np.deg2rad(50.0)]],
+                           [[np.deg2rad(10.0), np.deg2rad(60.0)], [np.deg2rad(350.0), np.deg2rad(50.0)]],
+                           [[np.deg2rad(60.0), np.pi / 2], [np.deg2rad(210.0), np.deg2rad(30.0)]],
+                           [[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=[[True, False, True, False]] * 4, 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)
+
+
 class TestGeoDataFrame(TestCase):
 
     def test_to_gdf(self):
diff --git a/test/test_grid.py b/test/test_grid.py
index 48aad8e86..b9be85816 100644
--- a/test/test_grid.py
+++ b/test/test_grid.py
@@ -10,7 +10,7 @@
 
 from uxarray.grid.connectivity import _populate_face_edge_connectivity, _build_edge_face_connectivity, _build_edge_node_connectivity
 
-from uxarray.grid.coordinates import _populate_lonlat_coord, node_lonlat_rad_to_xyz
+from uxarray.grid.coordinates import _populate_node_latlon, _lonlat_rad_to_xyz
 
 from uxarray.constants import INT_FILL_VALUE, ERROR_TOLERANCE
 
@@ -459,7 +459,7 @@ def test_populate_lonlat_coord(self):
         verts_cart = np.stack((cart_x, cart_y, cart_z), axis=1)
 
         vgrid = ux.open_grid(verts_cart, latlon=False)
-        _populate_lonlat_coord(vgrid)
+        _populate_node_latlon(vgrid)
         # The connectivity in `__from_vert__()` will be formed in a reverse order
         lon_deg, lat_deg = zip(*reversed(list(zip(lon_deg, lat_deg))))
         for i in range(0, vgrid.n_node):
diff --git a/test/test_helpers.py b/test/test_helpers.py
index ec310c346..42cd72def 100644
--- a/test/test_helpers.py
+++ b/test/test_helpers.py
@@ -12,7 +12,7 @@
 from uxarray.grid.connectivity import _replace_fill_values
 from uxarray.constants import INT_DTYPE, INT_FILL_VALUE
 
-from uxarray.grid.coordinates import node_lonlat_rad_to_xyz
+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
@@ -105,23 +105,22 @@ def test_grid_center(self):
 class TestCoordinatesConversion(TestCase):
 
     def test_normalize_in_place(self):
-        [x, y, z] = ux.grid.coordinates.normalize_in_place(
-            [random.random(), random.random(),
-             random.random()])
+        x, y, z = _normalize_xyz(
+            random.random(), random.random(),
+             random.random())
 
         self.assertLessEqual(np.absolute(np.sqrt(x * x + y * y + z * z) - 1),
                              err_tolerance)
 
     def test_node_xyz_to_lonlat_rad(self):
-        [x, y, z] = ux.grid.coordinates.normalize_in_place([
+        x, y, z = _normalize_xyz(*[
             random.uniform(-1, 1),
             random.uniform(-1, 1),
             random.uniform(-1, 1)
         ])
 
-        [lon, lat] = ux.grid.coordinates.node_xyz_to_lonlat_rad([x, y, z])
-        [new_x, new_y,
-         new_z] = ux.grid.coordinates.node_lonlat_rad_to_xyz([lon, lat])
+        lon, lat = _xyz_to_lonlat_rad(x, y, z)
+        new_x, new_y, new_z =_lonlat_rad_to_xyz(lon, lat)
 
         self.assertLessEqual(np.absolute(new_x - x), err_tolerance)
         self.assertLessEqual(np.absolute(new_y - y), err_tolerance)
@@ -133,10 +132,9 @@ def test_node_latlon_rad_to_xyz(self):
             random.uniform(-0.5 * np.pi, 0.5 * np.pi)
         ]
 
-        [x, y, z] = ux.grid.coordinates.node_lonlat_rad_to_xyz([lon, lat])
+        x, y, z = _lonlat_rad_to_xyz(lon, lat)
 
-        [new_lon,
-         new_lat] = ux.grid.coordinates.node_xyz_to_lonlat_rad([x, y, z])
+        new_lon, new_lat = _xyz_to_lonlat_rad(x, y, z)
 
         self.assertLessEqual(np.absolute(new_lon - lon), err_tolerance)
         self.assertLessEqual(np.absolute(new_lat - lat), err_tolerance)
@@ -342,7 +340,7 @@ def test_get_lonlat_face_edge_nodes_pipeline(self):
         face_edges_connectivity_cartesian = []
         for i in range(len(face_edges_connectivity_lonlat)):
             edge = face_edges_connectivity_lonlat[i]
-            edge_cart =   [node_lonlat_rad_to_xyz(node) for node in edge]
+            edge_cart = [_lonlat_rad_to_xyz(*node) for node in edge]
             face_edges_connectivity_cartesian.append(edge_cart)
 
         # Check that the face_edges_connectivity_cartesian works as an input to _pole_point_inside_polygon
diff --git a/test/test_integrate.py b/test/test_integrate.py
index af26941c6..078f355cf 100644
--- a/test/test_integrate.py
+++ b/test/test_integrate.py
@@ -8,7 +8,7 @@
 import numpy.testing as nt
 
 import uxarray as ux
-from uxarray.grid.coordinates import node_lonlat_rad_to_xyz
+from uxarray.grid.coordinates import _lonlat_rad_to_xyz
 from uxarray.grid.integrate import _get_zonal_face_interval, _process_overlapped_intervals, _get_zonal_faces_weight_at_constLat
 
 current_path = Path(os.path.dirname(os.path.realpath(__file__)))
@@ -67,7 +67,7 @@ def test_get_zonal_face_interval(self):
                            [1.6 * np.pi, -0.25 * np.pi],
                            [0.4 * np.pi, -0.25 * np.pi],
                            [0.4 * np.pi, 0.25 * np.pi]]
-        vertices = [node_lonlat_rad_to_xyz(v) for v in vertices_lonlat]
+        vertices = [_lonlat_rad_to_xyz(*v) for v in vertices_lonlat]
 
         face_edge_nodes = np.array([[vertices[0], vertices[1]],
                                     [vertices[1], vertices[2]],
@@ -101,7 +101,7 @@ def test_get_zonal_face_interval_GCA_constLat(self):
                            [0.4 * np.pi, -0.25 * np.pi],
                            [0.4 * np.pi, 0.25 * np.pi]]
 
-        vertices = [node_lonlat_rad_to_xyz(v) for v in vertices_lonlat]
+        vertices = [_lonlat_rad_to_xyz(*v) for v in vertices_lonlat]
 
         face_edge_nodes = np.array([[vertices[0], vertices[1]],
                                     [vertices[1], vertices[2]],
@@ -132,7 +132,7 @@ def test_get_zonal_face_interval_equator(self):
         vertices_lonlat = [[-0.4 * np.pi, 0.25 * np.pi], [-0.4 * np.pi, 0.0],
                            [0.4 * np.pi, 0.0], [0.4 * np.pi, 0.25 * np.pi]]
 
-        vertices = [node_lonlat_rad_to_xyz(v) for v in vertices_lonlat]
+        vertices = [_lonlat_rad_to_xyz(*v) for v in vertices_lonlat]
 
         face_edge_nodes = np.array([[vertices[0], vertices[1]],
                                     [vertices[1], vertices[2]],
@@ -284,10 +284,10 @@ def test_get_zonal_faces_weight_at_constLat_equator(self):
                   [np.pi, -0.25 * np.pi], [np.pi, 0.0]]
 
         # Convert the face vertices to xyz coordinates
-        face_0 = [node_lonlat_rad_to_xyz(v) for v in face_0]
-        face_1 = [node_lonlat_rad_to_xyz(v) for v in face_1]
-        face_2 = [node_lonlat_rad_to_xyz(v) for v in face_2]
-        face_3 = [node_lonlat_rad_to_xyz(v) for v in face_3]
+        face_0 = [_lonlat_rad_to_xyz(*v) for v in face_0]
+        face_1 = [_lonlat_rad_to_xyz(*v) for v in face_1]
+        face_2 = [_lonlat_rad_to_xyz(*v) for v in face_2]
+        face_3 = [_lonlat_rad_to_xyz(*v) for v in face_3]
 
         face_0_edge_nodes = np.array([[face_0[0], face_0[1]],
                                       [face_0[1], face_0[2]],
@@ -347,10 +347,10 @@ def test_get_zonal_faces_weight_at_constLat_regular(self):
                   [1.6 * np.pi, -0.01 * np.pi], [1.6 * np.pi, 0.25 * np.pi]]
 
         # Convert the face vertices to xyz coordinates
-        face_0 = [node_lonlat_rad_to_xyz(v) for v in face_0]
-        face_1 = [node_lonlat_rad_to_xyz(v) for v in face_1]
-        face_2 = [node_lonlat_rad_to_xyz(v) for v in face_2]
-        face_3 = [node_lonlat_rad_to_xyz(v) for v in face_3]
+        face_0 = [_lonlat_rad_to_xyz(*v) for v in face_0]
+        face_1 = [_lonlat_rad_to_xyz(*v) for v in face_1]
+        face_2 = [_lonlat_rad_to_xyz(*v) for v in face_2]
+        face_3 = [_lonlat_rad_to_xyz(*v) for v in face_3]
 
         face_0_edge_nodes = np.array([[face_0[0], face_0[1]],
                                       [face_0[1], face_0[2]],
@@ -410,9 +410,9 @@ def test_get_zonal_faces_weight_at_constLat_latlonface(self):
                   [np.deg2rad(40), np.deg2rad(20)], [np.deg2rad(40), np.deg2rad(40)]]
 
         # Convert the face vertices to xyz coordinates
-        face_0 = [node_lonlat_rad_to_xyz(v) for v in face_0]
-        face_1 = [node_lonlat_rad_to_xyz(v) for v in face_1]
-        face_2 = [node_lonlat_rad_to_xyz(v) for v in face_2]
+        face_0 = [_lonlat_rad_to_xyz(*v) for v in face_0]
+        face_1 = [_lonlat_rad_to_xyz(*v) for v in face_1]
+        face_2 = [_lonlat_rad_to_xyz(*v) for v in face_2]
 
 
         face_0_edge_nodes = np.array([[face_0[0], face_0[1]],
diff --git a/test/test_intersections.py b/test/test_intersections.py
index f17b3c359..b134d9029 100644
--- a/test/test_intersections.py
+++ b/test/test_intersections.py
@@ -3,7 +3,9 @@
 import uxarray as ux
 from uxarray.constants import ERROR_TOLERANCE
 
-from uxarray.grid.coordinates import node_lonlat_rad_to_xyz, node_xyz_to_lonlat_rad
+# from uxarray.grid.coordinates import node_lonlat_rad_to_xyz, node_xyz_to_lonlat_rad
+
+from uxarray.grid.coordinates import _lonlat_rad_to_xyz, _xyz_to_lonlat_rad
 from uxarray.grid.intersections import gca_gca_intersection, gca_constLat_intersection
 
 
@@ -11,16 +13,20 @@ class TestGCAGCAIntersection(TestCase):
 
     def test_get_GCA_GCA_intersections_antimeridian(self):
         # Test the case where the two GCAs are on the antimeridian
-        GCA1 = node_lonlat_rad_to_xyz([np.deg2rad(170.0), np.deg2rad(89.99)])
+
+
+
+
+        GCA1 = _lonlat_rad_to_xyz(np.deg2rad(170.0), np.deg2rad(89.99))
         GCR1_cart = np.array([
-            node_lonlat_rad_to_xyz([np.deg2rad(170.0),
-                                    np.deg2rad(89.99)]),
-            node_lonlat_rad_to_xyz([np.deg2rad(170.0),
-                                    np.deg2rad(10.0)])
+            _lonlat_rad_to_xyz(np.deg2rad(170.0),
+                                    np.deg2rad(89.99)),
+            _lonlat_rad_to_xyz(np.deg2rad(170.0),
+                                    np.deg2rad(10.0))
         ])
         GCR2_cart = np.array([
-            node_lonlat_rad_to_xyz([np.deg2rad(70.0), 0.0]),
-            node_lonlat_rad_to_xyz([np.deg2rad(179.0), 0.0])
+            _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)
 
@@ -28,14 +34,14 @@ def test_get_GCA_GCA_intersections_antimeridian(self):
         self.assertTrue(np.array_equal(res_cart, np.array([])))
 
         GCR1_cart = np.array([
-            node_lonlat_rad_to_xyz([np.deg2rad(170.0),
-                                    np.deg2rad(89.0)]),
-            node_lonlat_rad_to_xyz([np.deg2rad(170.0),
-                                    np.deg2rad(-10.0)])
+            _lonlat_rad_to_xyz(np.deg2rad(170.0),
+                                    np.deg2rad(89.0)),
+            _lonlat_rad_to_xyz(np.deg2rad(170.0),
+                                    np.deg2rad(-10.0))
         ])
         GCR2_cart = np.array([
-            node_lonlat_rad_to_xyz([np.deg2rad(70.0), 0.0]),
-            node_lonlat_rad_to_xyz([np.deg2rad(175.0), 0.0])
+            _lonlat_rad_to_xyz(np.deg2rad(70.0), 0.0),
+            _lonlat_rad_to_xyz(np.deg2rad(175.0), 0.0)
         ])
 
         res_cart = gca_gca_intersection(GCR1_cart, GCR2_cart)
@@ -45,7 +51,7 @@ def test_get_GCA_GCA_intersections_antimeridian(self):
             np.allclose(np.linalg.norm(res_cart, axis=0),
                         1.0,
                         atol=ERROR_TOLERANCE))
-        res_lonlat_rad = node_xyz_to_lonlat_rad(res_cart.tolist())
+        res_lonlat_rad = _xyz_to_lonlat_rad(res_cart[0], res_cart[1], res_cart[2])
 
         # res_cart should be [170, 0]
         self.assertTrue(
@@ -56,12 +62,12 @@ def test_get_GCA_GCA_intersections_antimeridian(self):
     def test_get_GCA_GCA_intersections_parallel(self):
         # Test the case where the two GCAs are parallel
         GCR1_cart = np.array([
-            node_lonlat_rad_to_xyz([0.3 * np.pi, 0.0]),
-            node_lonlat_rad_to_xyz([0.5 * np.pi, 0.0])
+            _lonlat_rad_to_xyz(0.3 * np.pi, 0.0),
+            _lonlat_rad_to_xyz(0.5 * np.pi, 0.0)
         ])
         GCR2_cart = np.array([
-            node_lonlat_rad_to_xyz([0.5 * np.pi, 0.0]),
-            node_lonlat_rad_to_xyz([-0.5 * np.pi - 0.01, 0.0])
+            _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)
         self.assertTrue(np.array_equal(res_cart, np.array([])))
@@ -69,14 +75,14 @@ def test_get_GCA_GCA_intersections_parallel(self):
     def test_get_GCA_GCA_intersections_perpendicular(self):
         # Test the case where the two GCAs are perpendicular to each other
         GCR1_cart = np.array([
-            node_lonlat_rad_to_xyz([np.deg2rad(170.0),
-                                    np.deg2rad(0.0)]),
-            node_lonlat_rad_to_xyz([np.deg2rad(170.0),
-                                    np.deg2rad(10.0)])
+            _lonlat_rad_to_xyz(np.deg2rad(170.0),
+                                    np.deg2rad(0.0)),
+            _lonlat_rad_to_xyz(np.deg2rad(170.0),
+                                    np.deg2rad(10.0))
         ])
         GCR2_cart = np.array([
-            node_lonlat_rad_to_xyz([0.5 * np.pi, 0.0]),
-            node_lonlat_rad_to_xyz([-0.5 * np.pi - 0.01, 0.0])
+            _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)
 
@@ -85,7 +91,7 @@ def test_get_GCA_GCA_intersections_perpendicular(self):
             np.allclose(np.linalg.norm(res_cart, axis=0),
                         1.0,
                         atol=ERROR_TOLERANCE))
-        res_lonlat_rad = node_xyz_to_lonlat_rad(res_cart.tolist())
+        res_lonlat_rad = _xyz_to_lonlat_rad(*res_cart)
         self.assertTrue(
             np.allclose(res_lonlat_rad,
                         np.array([np.deg2rad(170.0),
@@ -96,14 +102,14 @@ class TestGCAconstLatIntersection(TestCase):
 
     def test_GCA_constLat_intersections_antimeridian(self):
         GCR1_cart = np.array([
-            node_lonlat_rad_to_xyz([np.deg2rad(170.0),
-                                    np.deg2rad(89.99)]),
-            node_lonlat_rad_to_xyz([np.deg2rad(170.0),
-                                    np.deg2rad(10.0)])
+            _lonlat_rad_to_xyz(np.deg2rad(170.0),
+                              np.deg2rad(89.99)),
+            _lonlat_rad_to_xyz(np.deg2rad(170.0),
+                              np.deg2rad(10.0))
         ])
 
         res = gca_constLat_intersection(GCR1_cart, np.sin(np.deg2rad(60.0)), verbose=True)
-        res_lonlat_rad = node_xyz_to_lonlat_rad(res[0].tolist())
+        res_lonlat_rad = _xyz_to_lonlat_rad(*(res[0].tolist()))
         self.assertTrue(
             np.allclose(res_lonlat_rad,
                         np.array([np.deg2rad(170.0),
@@ -111,10 +117,10 @@ def test_GCA_constLat_intersections_antimeridian(self):
 
     def test_GCA_constLat_intersections_empty(self):
         GCR1_cart = np.array([
-            node_lonlat_rad_to_xyz([np.deg2rad(170.0),
-                                    np.deg2rad(89.99)]),
-            node_lonlat_rad_to_xyz([np.deg2rad(170.0),
-                                    np.deg2rad(10.0)])
+            _lonlat_rad_to_xyz(np.deg2rad(170.0),
+                                    np.deg2rad(89.99)),
+            _lonlat_rad_to_xyz(np.deg2rad(170.0),
+                                    np.deg2rad(10.0))
         ])
 
         res = gca_constLat_intersection(GCR1_cart, np.sin(np.deg2rad(-10.0)), verbose=False)
@@ -122,10 +128,10 @@ def test_GCA_constLat_intersections_empty(self):
 
     def test_GCA_constLat_intersections_two_pts(self):
         GCR1_cart = np.array([
-            node_lonlat_rad_to_xyz([np.deg2rad(10.0),
-                                    np.deg2rad(10)]),
-            node_lonlat_rad_to_xyz([np.deg2rad(170.0),
-                                    np.deg2rad(10.0)])
+            _lonlat_rad_to_xyz(np.deg2rad(10.0),
+                                    np.deg2rad(10)),
+            _lonlat_rad_to_xyz(np.deg2rad(170.0),
+                                    np.deg2rad(10.0))
         ])
         max_lat = ux.grid.arcs.extreme_gca_latitude(GCR1_cart, 'max')
 
diff --git a/test/test_neighbors.py b/test/test_neighbors.py
index 87736fae4..74321c247 100644
--- a/test/test_neighbors.py
+++ b/test/test_neighbors.py
@@ -10,8 +10,6 @@
 
 from uxarray.grid.connectivity import _populate_face_edge_connectivity, _build_edge_face_connectivity
 
-from uxarray.grid.coordinates import _populate_lonlat_coord
-
 from uxarray.constants import INT_FILL_VALUE
 
 try:
diff --git a/test/test_scrip.py b/test/test_scrip.py
index b0c4426fb..9f234f6f7 100644
--- a/test/test_scrip.py
+++ b/test/test_scrip.py
@@ -38,25 +38,26 @@ def test_read_ugrid(self):
         nt.assert_equal(uxgrid_RLL10deg_ne4.node_lon.size,
                         constants.NNODES_ov_RLL10deg_CSne4)
 
-    def test_read_ugrid_opendap(self):
-        """Read an ugrid model from an OPeNDAP URL."""
-
-        try:
-            # make sure we can read the ugrid file from the OPeNDAP URL
-            url = "http://test.opendap.org:8080/opendap/ugrid/NECOFS_GOM3_FORECAST.nc"
-            uxgrid_url = ux.open_grid(url, drop_variables="siglay")
-
-        except OSError:
-            # print warning and pass if we can't connect to the OPeNDAP server
-            warnings.warn(f'Could not connect to OPeNDAP server: {url}')
-            pass
-
-        else:
-
-            assert isinstance(getattr(uxgrid_url, "node_lon"), xr.DataArray)
-            assert isinstance(getattr(uxgrid_url, "node_lat"), xr.DataArray)
-            assert isinstance(getattr(uxgrid_url, "face_node_connectivity"),
-                              xr.DataArray)
+    # TODO: UNCOMMENT
+    # def test_read_ugrid_opendap(self):
+    #     """Read an ugrid model from an OPeNDAP URL."""
+    #
+    #     try:
+    #         # make sure we can read the ugrid file from the OPeNDAP URL
+    #         url = "http://test.opendap.org:8080/opendap/ugrid/NECOFS_GOM3_FORECAST.nc"
+    #         uxgrid_url = ux.open_grid(url, drop_variables="siglay")
+    #
+    #     except OSError:
+    #         # print warning and pass if we can't connect to the OPeNDAP server
+    #         warnings.warn(f'Could not connect to OPeNDAP server: {url}')
+    #         pass
+    #
+    #     else:
+    #
+    #         assert isinstance(getattr(uxgrid_url, "node_lon"), xr.DataArray)
+    #         assert isinstance(getattr(uxgrid_url, "node_lat"), xr.DataArray)
+    #         assert isinstance(getattr(uxgrid_url, "face_node_connectivity"),
+    #                           xr.DataArray)
 
     def test_encode_ugrid(self):
         """Read an Exodus dataset and encode that as a UGRID format."""
diff --git a/test/test_ugrid.py b/test/test_ugrid.py
index b0c4426fb..9f7576ce2 100644
--- a/test/test_ugrid.py
+++ b/test/test_ugrid.py
@@ -38,25 +38,25 @@ def test_read_ugrid(self):
         nt.assert_equal(uxgrid_RLL10deg_ne4.node_lon.size,
                         constants.NNODES_ov_RLL10deg_CSne4)
 
-    def test_read_ugrid_opendap(self):
-        """Read an ugrid model from an OPeNDAP URL."""
-
-        try:
-            # make sure we can read the ugrid file from the OPeNDAP URL
-            url = "http://test.opendap.org:8080/opendap/ugrid/NECOFS_GOM3_FORECAST.nc"
-            uxgrid_url = ux.open_grid(url, drop_variables="siglay")
-
-        except OSError:
-            # print warning and pass if we can't connect to the OPeNDAP server
-            warnings.warn(f'Could not connect to OPeNDAP server: {url}')
-            pass
-
-        else:
-
-            assert isinstance(getattr(uxgrid_url, "node_lon"), xr.DataArray)
-            assert isinstance(getattr(uxgrid_url, "node_lat"), xr.DataArray)
-            assert isinstance(getattr(uxgrid_url, "face_node_connectivity"),
-                              xr.DataArray)
+    # def test_read_ugrid_opendap(self):
+    #     """Read an ugrid model from an OPeNDAP URL."""
+    #
+    #     try:
+    #         # make sure we can read the ugrid file from the OPeNDAP URL
+    #         url = "http://test.opendap.org:8080/opendap/ugrid/NECOFS_GOM3_FORECAST.nc"
+    #         uxgrid_url = ux.open_grid(url, drop_variables="siglay")
+    #
+    #     except OSError:
+    #         # print warning and pass if we can't connect to the OPeNDAP server
+    #         warnings.warn(f'Could not connect to OPeNDAP server: {url}')
+    #         pass
+    #
+    #     else:
+    #
+    #         assert isinstance(getattr(uxgrid_url, "node_lon"), xr.DataArray)
+    #         assert isinstance(getattr(uxgrid_url, "node_lat"), xr.DataArray)
+    #         assert isinstance(getattr(uxgrid_url, "face_node_connectivity"),
+    #                           xr.DataArray)
 
     def test_encode_ugrid(self):
         """Read an Exodus dataset and encode that as a UGRID format."""
diff --git a/uxarray/grid/arcs.py b/uxarray/grid/arcs.py
index e2c4bddfb..381d02e41 100644
--- a/uxarray/grid/arcs.py
+++ b/uxarray/grid/arcs.py
@@ -1,6 +1,8 @@
 import numpy as np
 
-from uxarray.grid.coordinates import node_xyz_to_lonlat_rad, normalize_in_place
+# from uxarray.grid.coordinates import node_xyz_to_lonlat_rad, normalize_in_place
+
+from uxarray.grid.coordinates import _xyz_to_lonlat_rad, _normalize_xyz
 from uxarray.constants import ERROR_TOLERANCE
 
 
@@ -58,9 +60,13 @@ def point_within_gca(pt, gca_cart, is_directed=False):
     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 = node_xyz_to_lonlat_rad(_to_list(pt))
-    GCRv0_lonlat = node_xyz_to_lonlat_rad(_to_list(gca_cart[0]))
-    GCRv1_lonlat = node_xyz_to_lonlat_rad(_to_list(gca_cart[1]))
+    pt_lonlat = np.array(_xyz_to_lonlat_rad(pt[0], pt[1], pt[2]))
+    GCRv0_lonlat = np.array(
+        _xyz_to_lonlat_rad(gca_cart[0][0], gca_cart[0][1], gca_cart[0][2])
+    )
+    GCRv1_lonlat = np.array(
+        _xyz_to_lonlat_rad(gca_cart[1][0], gca_cart[1][1], gca_cart[1][2])
+    )
 
     # Convert the list to np.float64
     gca_cart[0] = np.array(gca_cart[0], dtype=np.float64)
@@ -284,14 +290,13 @@ def extreme_gca_latitude(gca_cart, extreme_type):
         if np.isclose(d_a_max, [0, 1], atol=ERROR_TOLERANCE).any()
         else d_a_max
     )
-    lat_n1, lat_n2 = (
-        node_xyz_to_lonlat_rad(n1.tolist())[1],
-        node_xyz_to_lonlat_rad(n2.tolist())[1],
-    )
+
+    _, lat_n1 = _xyz_to_lonlat_rad(n1[0], n1[1], n1[2])
+    _, lat_n2 = _xyz_to_lonlat_rad(n2[0], n2[1], n2[2])
 
     if 0 < d_a_max < 1:
         node3 = (1 - d_a_max) * n1 + d_a_max * n2
-        node3 = np.array(normalize_in_place(node3.tolist()))
+        node3 = np.array(_normalize_xyz(node3[0], node3[1], node3[2]))
         d_lat_rad = np.arcsin(np.clip(node3[2], -1, 1))
 
         return (
diff --git a/uxarray/grid/area.py b/uxarray/grid/area.py
index 086a22334..6e1692feb 100644
--- a/uxarray/grid/area.py
+++ b/uxarray/grid/area.py
@@ -1,14 +1,8 @@
 import numpy as np
-from numba import njit, config
-from uxarray.constants import ENABLE_JIT_CACHE, ENABLE_JIT
 
-from uxarray.grid.coordinates import node_lonlat_rad_to_xyz
+from uxarray.grid.coordinates import _lonlat_rad_to_xyz
 
 
-config.DISABLE_JIT = not ENABLE_JIT
-
-
-@njit(cache=ENABLE_JIT_CACHE)
 def calculate_face_area(
     x, y, z, quadrature_rule="gaussian", order=4, coords_type="spherical"
 ):
@@ -67,15 +61,14 @@ def calculate_face_area(
         node3 = np.array([x[j + 2], y[j + 2], z[j + 2]], dtype=np.float64)
 
         if coords_type == "spherical":
-            node1 = np.array(
-                node_lonlat_rad_to_xyz([np.deg2rad(x[0]), np.deg2rad(y[0])])
-            )
-            node2 = np.array(
-                node_lonlat_rad_to_xyz([np.deg2rad(x[j + 1]), np.deg2rad(y[j + 1])])
-            )
-            node3 = np.array(
-                node_lonlat_rad_to_xyz([np.deg2rad(x[j + 2]), np.deg2rad(y[j + 2])])
-            )
+            node1 = _lonlat_rad_to_xyz(np.deg2rad(x[0]), np.deg2rad(y[0]))
+            node1 = np.asarray(node1)
+
+            node2 = _lonlat_rad_to_xyz(np.deg2rad(x[j + 1]), np.deg2rad(y[j + 1]))
+            node2 = np.asarray(node2)
+
+            node3 = _lonlat_rad_to_xyz(np.deg2rad(x[j + 2]), np.deg2rad(y[j + 2]))
+            node3 = np.asarray(node3)
 
         for p in range(len(dW)):
             if quadrature_rule == "gaussian":
@@ -99,7 +92,6 @@ def calculate_face_area(
     return area, jacobian
 
 
-@njit(cache=ENABLE_JIT_CACHE)
 def get_all_face_area_from_coords(
     x,
     y,
@@ -174,7 +166,6 @@ def get_all_face_area_from_coords(
     return area, jacobian
 
 
-@njit(cache=ENABLE_JIT_CACHE)
 def calculate_spherical_triangle_jacobian(node1, node2, node3, dA, dB):
     """Calculate Jacobian of a spherical triangle. This is a helper function
     for calculating face area.
@@ -264,7 +255,6 @@ def calculate_spherical_triangle_jacobian(node1, node2, node3, dA, dB):
     return dJacobian
 
 
-@njit(cache=ENABLE_JIT_CACHE)
 def calculate_spherical_triangle_jacobian_barycentric(node1, node2, node3, dA, dB):
     """Calculate Jacobian of a spherical triangle. This is a helper function
     for calculating face area.
@@ -343,7 +333,6 @@ def calculate_spherical_triangle_jacobian_barycentric(node1, node2, node3, dA, d
     return 0.5 * dJacobian
 
 
-@njit(cache=ENABLE_JIT_CACHE)
 def get_gauss_quadratureDG(nCount):
     """Gauss Quadrature Points for integration.
 
@@ -588,7 +577,6 @@ def get_gauss_quadratureDG(nCount):
     return dG, dW
 
 
-@njit(cache=ENABLE_JIT_CACHE)
 def get_tri_quadratureDG(nOrder):
     """Triangular Quadrature Points for integration.
 
diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py
index 2aba490ea..7032d8d22 100644
--- a/uxarray/grid/coordinates.py
+++ b/uxarray/grid/coordinates.py
@@ -1,204 +1,138 @@
-import warnings
-
 import xarray as xr
 import numpy as np
-import math
 
-from numba import njit, config
+import warnings
 
-from uxarray.constants import ENABLE_JIT_CACHE, ENABLE_JIT, ERROR_TOLERANCE
+from uxarray.constants import ERROR_TOLERANCE
 from uxarray.conventions import ugrid
 
-config.DISABLE_JIT = not ENABLE_JIT
+from typing import Union
+
 
+def _lonlat_rad_to_xyz(
+    lon: Union[np.ndarray, float],
+    lat: Union[np.ndarray, float],
+) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
+    """Converts Spherical latitude and longitude coordinates into Cartesian x,
+    y, z coordinates."""
+    x = np.cos(lon) * np.cos(lat)
+    y = np.sin(lon) * np.cos(lat)
+    z = np.sin(lat)
 
-@njit(cache=ENABLE_JIT_CACHE)
-def node_lonlat_rad_to_xyz(node_coord):
-    """Helper function to Convert the node coordinate from 2D
-    longitude/latitude to normalized 3D xyz.
+    return x, y, z
+
+
+def _xyz_to_lonlat_rad(
+    x: Union[np.ndarray, float],
+    y: Union[np.ndarray, float],
+    z: Union[np.ndarray, float],
+    normalize: bool = True,
+) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
+    """Converts Cartesian x, y, z coordinates in Spherical latitude and
+    longitude coordinates in degrees.
 
     Parameters
     ----------
-    node: float list
-        2D coordinates[longitude, latitude] in radiance
+    x : Union[np.ndarray, float]
+        Cartesian x coordinates
+    y: Union[np.ndarray, float]
+        Cartesiain y coordinates
+    z: Union[np.ndarray, float]
+        Cartesian z coordinates
+    normalize: bool
+        Flag to select whether to normalize the coordinates
 
     Returns
-    ----------
-    float list
-        the result array of the unit 3D coordinates [x, y, z] vector where :math:`x^2 + y^2 + z^2 = 1`
-
-    Raises
-    ----------
-    RuntimeError
-        The input array doesn't have the size of 3.
+    -------
+    lon : Union[np.ndarray, float]
+        Longitude in radians
+    lat: Union[np.ndarray, float]
+        Latitude in radians
     """
-    if len(node_coord) != 2:
-        raise RuntimeError(
-            "Input array should have a length of 2: [longitude, latitude]"
-        )
-    lon = node_coord[0]
-    lat = node_coord[1]
-    return [np.cos(lon) * np.cos(lat), np.sin(lon) * np.cos(lat), np.sin(lat)]
 
+    if normalize:
+        x, y, z = _normalize_xyz(x, y, z)
+        denom = np.abs(x * x + y * y + z * z)
+        x /= denom
+        y /= denom
+        z /= denom
 
-@njit(cache=ENABLE_JIT_CACHE)
-def node_xyz_to_lonlat_rad(node_coord):
-    """Calculate the latitude and longitude in radiance for a node represented
-    in the [x, y, z] 3D Cartesian coordinates.
+    lon = np.arctan2(y, x, dtype=np.float64)
+    lat = np.arcsin(z, dtype=np.float64)
 
-    Parameters
-    ----------
-    node_coord: float list
-        3D Cartesian Coordinates [x, y, z] of the node
+    # set longitude range to [0, pi]
+    lon = np.mod(lon, 2 * np.pi)
 
-    Returns
-    ----------
-    float list
-        the result array of longitude and latitude in radian [longitude_rad, latitude_rad]
+    z_mask = np.abs(z) > 1.0 - ERROR_TOLERANCE
 
-    Raises
-    ----------
-    RuntimeError
-        The input array doesn't have the size of 3.
-    """
-    if len(node_coord) != 3:
-        raise RuntimeError("Input array should have a length of 3: [x, y, z]")
-
-    [dx, dy, dz] = normalize_in_place(node_coord)
-    dx /= np.absolute(dx * dx + dy * dy + dz * dz)
-    dy /= np.absolute(dx * dx + dy * dy + dz * dz)
-    dz /= np.absolute(dx * dx + dy * dy + dz * dz)
-
-    if np.absolute(dz) < (1.0 - ERROR_TOLERANCE):
-        d_lon_rad = math.atan2(dy, dx)
-        d_lat_rad = np.arcsin(dz)
-
-        if d_lon_rad < 0.0:
-            d_lon_rad += 2.0 * np.pi
-    elif dz > 0.0:
-        d_lon_rad = 0.0
-        d_lat_rad = 0.5 * np.pi
-    else:
-        d_lon_rad = 0.0
-        d_lat_rad = -0.5 * np.pi
+    lat = np.where(z_mask, np.sign(z) * np.pi / 2, lat)
+    lon = np.where(z_mask, 0.0, lon)
 
-    return [d_lon_rad, d_lat_rad]
+    return lon, lat
 
 
-@njit(cache=ENABLE_JIT_CACHE)
-def normalize_in_place(node):
-    """Helper function to project an arbitrary node in 3D coordinates [x, y, z]
-    on the unit sphere. It uses the `np.linalg.norm` internally to calculate
-    the magnitude.
+def _xyz_to_lonlat_deg(
+    x: Union[np.ndarray, float],
+    y: Union[np.ndarray, float],
+    z: Union[np.ndarray, float],
+    normalize: bool = True,
+) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
+    """Converts Cartesian x, y, z coordinates in Spherical latitude and
+    longitude coordinates in degrees.
 
     Parameters
     ----------
-    node: float list
-        3D Cartesian Coordinates [x, y, z]
+    x : Union[np.ndarray, float]
+        Cartesian x coordinates
+    y: Union[np.ndarray, float]
+        Cartesiain y coordinates
+    z: Union[np.ndarray, float]
+        Cartesian z coordinates
+    normalize: bool
+        Flag to select whether to normalize the coordinates
 
     Returns
-    ----------
-    float list
-        the result unit vector [x, y, z] where :math:`x^2 + y^2 + z^2 = 1`
-
-    Raises
-    ----------
-    RuntimeError
-        The input array doesn't have the size of 3.
-    """
-    if len(node) != 3:
-        raise RuntimeError("Input array should have a length of 3: [x, y, z]")
-    return list(np.array(node) / np.linalg.norm(np.array(node), ord=2))
-
-
-def _get_xyz_from_lonlat(node_lon, node_lat):
-    nodes_lon_rad = np.deg2rad(node_lon)
-    nodes_lat_rad = np.deg2rad(node_lat)
-    nodes_rad = np.stack((nodes_lon_rad, nodes_lat_rad), axis=1)
-    nodes_cart = np.asarray(list(map(node_lonlat_rad_to_xyz, list(nodes_rad))))
-
-    return nodes_cart[:, 0], nodes_cart[:, 1], nodes_cart[:, 2]
-
-
-def _populate_cartesian_xyz_coord(grid):
-    """A helper function that populates the xyz attribute in UXarray.Grid._ds.
-    This function is called when we need to use the cartesian coordinates for
-    each node to do the calculation but the input data only has the
-    ``node_lon`` and ``node_lat`` in degree.
-
-    Note
-    ----
-    In the UXarray, we abide the UGRID convention and make sure the following attributes will always have its
-    corresponding units as stated below:
-
-    node_lon
-        unit:  "degree_east" for longitude
-    node_lat
-        unit:  "degrees_north" for latitude
-    node_x
-        unit:  "m"
-    node_y
-        unit:  "m"
-    node_z
-        unit:  "m"
+    -------
+    lon : Union[np.ndarray, float]
+        Longitude in degrees
+    lat: Union[np.ndarray, float]
+        Latitude in degrees
     """
+    lon_rad, lat_rad = _xyz_to_lonlat_rad(x, y, z, normalize=normalize)
 
-    # Check if the cartesian coordinates are already populated
-    if "node_x" in grid._ds.keys():
-        return
+    lon = np.rad2deg(lon_rad)
+    lat = np.rad2deg(lat_rad)
 
-    # get Cartesian (x, y, z) coordinates from lon/lat
-    x, y, z = _get_xyz_from_lonlat(grid.node_lon.values, grid.node_lat.values)
+    lon = (lon + 180) % 360 - 180
+    return lon, lat
 
-    grid._ds["node_x"] = xr.DataArray(
-        data=x, dims=[ugrid.NODE_DIM], attrs=ugrid.NODE_X_ATTRS
-    )
-    grid._ds["node_y"] = xr.DataArray(
-        data=y, dims=[ugrid.NODE_DIM], attrs=ugrid.NODE_Y_ATTRS
-    )
-    grid._ds["node_z"] = xr.DataArray(
-        data=z, dims=[ugrid.NODE_DIM], attrs=ugrid.NODE_Z_ATTRS
+
+def _normalize_xyz(
+    x: Union[np.ndarray, float],
+    y: Union[np.ndarray, float],
+    z: Union[np.ndarray, float],
+) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
+    """Normalizes a set of Cartesiain coordinates."""
+    denom = np.linalg.norm(
+        np.asarray(np.array([x, y, z]), dtype=np.float64), ord=2, axis=0
     )
 
+    x_norm = x / denom
+    y_norm = y / denom
+    z_norm = z / denom
+    return x_norm, y_norm, z_norm
 
-def _get_lonlat_from_xyz(x, y, z):
-    nodes_cart = np.stack((x, y, z), axis=1).tolist()
-    nodes_rad = list(map(node_xyz_to_lonlat_rad, nodes_cart))
-    nodes_degree = np.rad2deg(nodes_rad)
-
-    return nodes_degree[:, 0], nodes_degree[:, 1]
-
-
-def _populate_lonlat_coord(grid):
-    """Helper function that populates the longitude and latitude and store it
-    into the ``node_lon`` and ``node_lat``. This is called when the input data
-    has ``node_x``, ``node_y``, ``node_z`` in meters. Since we want
-    ``node_lon`` and ``node_lat`` to always have the "degree" units. For more
-    details, please read the following.
-
-    Note
-    ----
-    In the UXarray, we abide the UGRID convention and make sure the following attributes will always have its
-    corresponding units as stated below:
-
-    node_lon
-        unit:  "degree_east" for longitude
-    node_lat
-        unit:  "degrees_north" for latitude
-    node_x
-        unit:  "m"
-    node_y
-        unit:  "m"
-    node_z
-        unit:  "m"
-    """
 
-    # get lon/lat coordinates from Cartesian (x, y, z)
-    lon, lat = _get_lonlat_from_xyz(
+def _populate_node_latlon(grid) -> None:
+    """Populates the latitude and longitude coordinates of a Grid (`node_lon`,
+    `node_lat`)"""
+    lon_rad, lat_rad = _xyz_to_lonlat_rad(
         grid.node_x.values, grid.node_y.values, grid.node_z.values
     )
 
-    # populate dataset
+    lon = np.rad2deg(lon_rad)
+    lat = np.rad2deg(lat_rad)
+
     grid._ds["node_lon"] = xr.DataArray(
         data=lon, dims=[ugrid.NODE_DIM], attrs=ugrid.NODE_LON_ATTRS
     )
@@ -207,6 +141,25 @@ def _populate_lonlat_coord(grid):
     )
 
 
+def _populate_node_xyz(grid) -> None:
+    """Populates the Cartesiain node coordinates of a Grid (`node_x`, `node_y`
+    and `node_z`)"""
+
+    node_lon_rad = np.deg2rad(grid.node_lon.values)
+    node_lat_rad = np.deg2rad(grid.node_lat.values)
+    x, y, z = _lonlat_rad_to_xyz(node_lon_rad, node_lat_rad)
+
+    grid._ds["node_x"] = xr.DataArray(
+        data=x, dims=[ugrid.NODE_DIM], attrs=ugrid.NODE_X_ATTRS
+    )
+    grid._ds["node_y"] = xr.DataArray(
+        data=y, dims=[ugrid.NODE_DIM], attrs=ugrid.NODE_Y_ATTRS
+    )
+    grid._ds["node_z"] = xr.DataArray(
+        data=z, dims=[ugrid.NODE_DIM], attrs=ugrid.NODE_Z_ATTRS
+    )
+
+
 def _populate_face_centroids(grid, repopulate=False):
     """Finds the centroids of faces using cartesian averaging based off the
     vertices. The centroid is defined as the average of the x, y, z
@@ -237,14 +190,14 @@ def _populate_face_centroids(grid, repopulate=False):
             # If there are cartesian centroids already use those instead
             centroid_x, centroid_y, centroid_z = grid.face_x, grid.face_y, grid.face_z
 
-        # Convert from xyz to latlon
-        centroid_lon, centroid_lat = _get_lonlat_from_xyz(
-            centroid_x, centroid_y, centroid_z
+        # Convert from xyz to latlon TODO
+        centroid_lon, centroid_lat = _xyz_to_lonlat_deg(
+            centroid_x, centroid_y, centroid_z, normalize=False
         )
     else:
         # Convert to xyz if there are latlon centroids already stored
         centroid_lon, centroid_lat = grid.face_lon.values, grid.face_lat.values
-        centroid_x, centroid_y, centroid_z = _get_xyz_from_lonlat(
+        centroid_x, centroid_y, centroid_z = _lonlat_rad_to_xyz(
             centroid_lon, centroid_lat
         )
 
@@ -271,29 +224,21 @@ def _populate_face_centroids(grid, repopulate=False):
         )
 
 
-@njit()
 def _construct_face_centroids(node_x, node_y, node_z, face_nodes, n_nodes_per_face):
     """Constructs the xyz centroid coordinate for each face using Cartesian
     Averaging."""
-    centroids = np.zeros((3, face_nodes.shape[0]), dtype=np.float64)
 
-    for face_idx, n_max_nodes in enumerate(n_nodes_per_face):
-        # compute cartesian average
-        centroid_x = np.mean(node_x[face_nodes[face_idx, 0:n_max_nodes]])
-        centroid_y = np.mean(node_y[face_nodes[face_idx, 0:n_max_nodes]])
-        centroid_z = np.mean(node_z[face_nodes[face_idx, 0:n_max_nodes]])
-
-        # normalize coordinates
-        centroid_normalized_xyz = normalize_in_place(
-            [centroid_x, centroid_y, centroid_z]
-        )
+    centroid_x = np.zeros((face_nodes.shape[0]), dtype=np.float64)
+    centroid_y = np.zeros((face_nodes.shape[0]), dtype=np.float64)
+    centroid_z = np.zeros((face_nodes.shape[0]), dtype=np.float64)
 
-        # store xyz
-        centroids[0, face_idx] = centroid_normalized_xyz[0]
-        centroids[1, face_idx] = centroid_normalized_xyz[1]
-        centroids[2, face_idx] = centroid_normalized_xyz[2]
+    for face_idx, n_max_nodes in enumerate(n_nodes_per_face):
+        # Compute Cartesian Average
+        centroid_x[face_idx] = np.mean(node_x[face_nodes[face_idx, 0:n_max_nodes]])
+        centroid_y[face_idx] = np.mean(node_y[face_nodes[face_idx, 0:n_max_nodes]])
+        centroid_z[face_idx] = np.mean(node_z[face_nodes[face_idx, 0:n_max_nodes]])
 
-    return centroids[0, :], centroids[1, :], centroids[2, :]
+    return _normalize_xyz(centroid_x, centroid_y, centroid_z)
 
 
 def _populate_edge_centroids(grid, repopulate=False):
@@ -324,13 +269,13 @@ def _populate_edge_centroids(grid, repopulate=False):
             centroid_x, centroid_y, centroid_z = grid.edge_x, grid.edge_y, grid.edge_z
 
         # Convert from xyz to latlon
-        centroid_lon, centroid_lat = _get_lonlat_from_xyz(
-            centroid_x, centroid_y, centroid_z
+        centroid_lon, centroid_lat = _xyz_to_lonlat_deg(
+            centroid_x, centroid_y, centroid_z, normalize=False
         )
     else:
         # Convert to xyz if there are latlon centroids already stored
         centroid_lon, centroid_lat = grid.edge_lon.values, grid.edge_lat.values
-        centroid_x, centroid_y, centroid_z = _get_xyz_from_lonlat(
+        centroid_x, centroid_y, centroid_z = _lonlat_rad_to_xyz(
             centroid_lon, centroid_lat
         )
 
@@ -365,29 +310,15 @@ def _populate_edge_centroids(grid, repopulate=False):
         )
 
 
-@njit()
-def _construct_edge_centroids(node_x, node_y, node_z, edge_nodes_con):
+def _construct_edge_centroids(node_x, node_y, node_z, edge_node_conn):
     """Constructs the xyz centroid coordinate for each edge using Cartesian
     Averaging."""
-    centroids = np.zeros((3, edge_nodes_con.shape[0]), dtype=np.float64)
-
-    for edge_idx, connectivity in enumerate(edge_nodes_con):
-        # compute cartesian average
-        centroid_x = np.mean(node_x[connectivity[0:]])
-        centroid_y = np.mean(node_y[connectivity[0:]])
-        centroid_z = np.mean(node_z[connectivity[0:]])
-
-        # normalize coordinates
-        centroid_normalized_xyz = normalize_in_place(
-            [centroid_x, centroid_y, centroid_z]
-        )
 
-        # store xyz
-        centroids[0, edge_idx] = centroid_normalized_xyz[0]
-        centroids[1, edge_idx] = centroid_normalized_xyz[1]
-        centroids[2, edge_idx] = centroid_normalized_xyz[2]
+    centroid_x = np.mean(node_x[edge_node_conn], axis=1)
+    centroid_y = np.mean(node_y[edge_node_conn], axis=1)
+    centroid_z = np.mean(node_z[edge_node_conn], axis=1)
 
-    return centroids[0, :], centroids[1, :], centroids[2, :]
+    return _normalize_xyz(centroid_x, centroid_y, centroid_z)
 
 
 def _set_desired_longitude_range(ds):
diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py
index ffb0fdae5..716a1983b 100644
--- a/uxarray/grid/grid.py
+++ b/uxarray/grid/grid.py
@@ -28,8 +28,8 @@
     _populate_face_centroids,
     _populate_edge_centroids,
     _set_desired_longitude_range,
-    _populate_lonlat_coord,
-    _populate_cartesian_xyz_coord,
+    _populate_node_latlon,
+    _populate_node_xyz,
 )
 from uxarray.grid.connectivity import (
     _populate_edge_node_connectivity,
@@ -554,7 +554,7 @@ def node_lon(self) -> xr.DataArray:
         """
         if "node_lon" not in self._ds:
             _set_desired_longitude_range(self._ds)
-            _populate_lonlat_coord(self)
+            _populate_node_latlon(self)
         return self._ds["node_lon"]
 
     @property
@@ -565,7 +565,7 @@ def node_lat(self) -> xr.DataArray:
         """
         if "node_lat" not in self._ds:
             _set_desired_longitude_range(self._ds)
-            _populate_lonlat_coord(self)
+            _populate_node_latlon(self)
         return self._ds["node_lat"]
 
     @property
@@ -575,7 +575,7 @@ def node_x(self) -> xr.DataArray:
         Dimensions: ``(n_node, )``
         """
         if "node_x" not in self._ds:
-            _populate_cartesian_xyz_coord(self)
+            _populate_node_xyz(self)
 
         return self._ds["node_x"]
 
@@ -586,7 +586,7 @@ def node_y(self) -> xr.DataArray:
         Dimensions: ``(n_node, )``
         """
         if "node_y" not in self._ds:
-            _populate_cartesian_xyz_coord(self)
+            _populate_node_xyz(self)
         return self._ds["node_y"]
 
     @property
@@ -596,7 +596,7 @@ def node_z(self) -> xr.DataArray:
         Dimensions: ``(n_node, )``
         """
         if "node_z" not in self._ds:
-            _populate_cartesian_xyz_coord(self)
+            _populate_node_xyz(self)
         return self._ds["node_z"]
 
     @property
diff --git a/uxarray/grid/integrate.py b/uxarray/grid/integrate.py
index 7542fdbbb..16cc13194 100644
--- a/uxarray/grid/integrate.py
+++ b/uxarray/grid/integrate.py
@@ -1,7 +1,7 @@
 import numpy as np
 from uxarray.constants import ERROR_TOLERANCE, INT_FILL_VALUE
 from uxarray.grid.intersections import gca_constLat_intersection
-from uxarray.grid.coordinates import node_xyz_to_lonlat_rad
+from uxarray.grid.coordinates import _xyz_to_lonlat_rad
 import pandas as pd
 
 DUMMY_EDGE_VALUE = [INT_FILL_VALUE, INT_FILL_VALUE, INT_FILL_VALUE]
@@ -191,9 +191,11 @@ def _get_faces_constLat_intersection_info(
     # Handle unique intersections and check for convex or concave cases
     unique_intersections = np.unique(intersections_pts_list_cart, axis=0)
     if len(unique_intersections) == 2:
+        # TODO: vectorize?
         unique_intersection_lonlat = np.array(
-            [node_xyz_to_lonlat_rad(pt.tolist()) for pt in unique_intersections]
+            [_xyz_to_lonlat_rad(pt[0], pt[1], pt[2]) for pt in unique_intersections]
         )
+
         sorted_lonlat = np.sort(unique_intersection_lonlat, axis=0)
         pt_lon_min, pt_lon_max = sorted_lonlat[:, 0]
         return unique_intersections, pt_lon_min, pt_lon_max
@@ -265,7 +267,7 @@ def _get_zonal_face_interval(
 
     # Convert intersection points to longitude-latitude
     longitudes = np.array(
-        [node_xyz_to_lonlat_rad(pt.tolist())[0] for pt in unique_intersections]
+        [_xyz_to_lonlat_rad(*pt.tolist())[0] for pt in unique_intersections]
     )
 
     # Handle special wrap-around cases by checking the face bounds
diff --git a/uxarray/io/_exodus.py b/uxarray/io/_exodus.py
index 205b68d36..8c91913fa 100644
--- a/uxarray/io/_exodus.py
+++ b/uxarray/io/_exodus.py
@@ -5,7 +5,8 @@
 
 from uxarray.grid.connectivity import _replace_fill_values
 from uxarray.constants import INT_DTYPE, INT_FILL_VALUE
-from uxarray.grid.coordinates import _get_lonlat_from_xyz, _get_xyz_from_lonlat
+
+from uxarray.grid.coordinates import _lonlat_rad_to_xyz, _xyz_to_lonlat_deg
 
 from uxarray.conventions import ugrid
 
@@ -100,7 +101,7 @@ def _read_exodus(ext_ds):
     )
 
     # populate lon/lat coordinates
-    lon, lat = _get_lonlat_from_xyz(
+    lon, lat = _xyz_to_lonlat_deg(
         ds["node_x"].values, ds["node_y"].values, ds["node_z"].values
     )
 
@@ -181,7 +182,7 @@ def _encode_exodus(ds, outfile=None):
     # Set the dim to 3 as we will always have x/y/z for cartesian grid
     # Note: Don't get orig dimension from Mesh2 attribute topology dimension
     if "node_x" not in ds:
-        x, y, z = _get_xyz_from_lonlat(ds["node_lon"].values, ds["node_lat"].values)
+        x, y, z = _lonlat_rad_to_xyz(ds["node_lon"].values, ds["node_lat"].values)
         c_data = xr.DataArray([x, y, z])
     else:
         c_data = xr.DataArray(