diff --git a/benchmarks/mpas_ocean.py b/benchmarks/mpas_ocean.py
index f10fa6918..3343c969f 100644
--- a/benchmarks/mpas_ocean.py
+++ b/benchmarks/mpas_ocean.py
@@ -144,6 +144,22 @@ class HoleEdgeIndices(DatasetBenchmark):
     def time_construct_hole_edge_indices(self, resolution):
         ux.grid.geometry._construct_hole_edge_indices(self.uxds.uxgrid.edge_face_connectivity)
 
+class ConstructFaceLatLon(GridBenchmark):
+    def time_welzl(self, resolution):
+        from uxarray.grid.coordinates import _construct_face_centerpoint
+        _construct_face_centerpoints(self.uxgrid.node_lon.values,
+                                     self.uxgrid.node_lat.values,
+                                     self.uxgrid.face_node_connectivity.values,
+                                     self.uxgrid.n_nodes_per_face.values)
+
+    def time_cartesian_averaging(self, resolution):
+        from uxarray.grid.coordinates import _construct_face_centroids
+        _construct_face_centroids(self.uxgrid.node_x.values,
+                                  self.uxgrid.node_y.values,
+                                  self.uxgrid.node_z.values,
+                                  self.uxgrid.face_node_connectivity.values,
+                                  self.uxgrid.n_nodes_per_face.values)
+
 class CheckNorm:
     param_names = ['resolution']
     params = ['480km', '120km']
diff --git a/benchmarks/quad_hexagon.py b/benchmarks/quad_hexagon.py
index 4364f39b0..aeb7559bb 100644
--- a/benchmarks/quad_hexagon.py
+++ b/benchmarks/quad_hexagon.py
@@ -34,3 +34,39 @@ def time_open_dataset(self):
     def peakmem_open_dataset(self):
         """Peak memory usage of a `UxDataset`"""
         uxds = ux.open_dataset(grid_path, data_path)
+
+
+# from uxarray.grid.coordinates import _construct_face_centerpoints, _construct_face_centroids
+
+class ConstructFaceLatLon:
+
+
+    def time_welzl(self):
+        uxgrid = ux.open_grid("/Users/mbook/cp_uxarray/benchmarks/oQU480.grid.nc")
+
+        ux.grid.coordinates._construct_face_centerpoints(uxgrid.node_lon.values,
+                                     uxgrid.node_lat.values,
+                                     uxgrid.face_node_connectivity.values,
+                                     uxgrid.n_nodes_per_face.values)
+
+
+
+
+
+class QuadHexagon2:
+    def time_open_grid(self):
+        """Time to open a `Grid`"""
+        ux.open_grid(grid_path)
+
+    # def mem_open_grid(self):
+    #     """Memory Occupied by a `Grid`"""
+    #     return ux.open_grid(grid_path)
+
+    def peakmem_open_grid2(self):
+        """Peak memory usage of a `Grid`"""
+        uxgrid = ux.open_grid(grid_path)
+
+
+    def time_open_dataset(self):
+        """Time to open a `UxDataset`"""
+        ux.open_dataset(grid_path, data_path)
diff --git a/docs/internal_api/index.rst b/docs/internal_api/index.rst
index 122a1d96e..4c6e82dbd 100644
--- a/docs/internal_api/index.rst
+++ b/docs/internal_api/index.rst
@@ -148,10 +148,6 @@ Coordinates
 .. autosummary::
    :toctree: generated/
 
-   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
@@ -159,6 +155,13 @@ Coordinates
    grid.coordinates._construct_face_centroids
    grid.coordinates._construct_edge_centroids
    grid.coordinates._set_desired_longitude_range
+   grid.coordinates._populate_face_centerpoints
+   grid.coordinates._circle_from_two_points
+   grid.coordinates._circle_from_three_points
+   grid.coordinates._is_inside_circle
+   grid.coordinates._welzl_recursive
+   grid.coordinates._smallest_enclosing_circle
+   grid.coordinates._construct_face_centerpoints
 
 
 Arcs
@@ -180,7 +183,12 @@ Utils
    grid.utils._swap_first_fill_value_with_last
    grid.utils._get_cartesiain_face_edge_nodes
    grid.utils._get_lonlat_rad_face_edge_nodes
-
+   grid.utils._lonlat_rad_to_xyz
+   grid.utils._xyz_to_lonlat_rad
+   grid.utils._xyz_to_lonlat_rad_no_norm
+   grid.utils._xyz_to_lonlat_deg
+   grid.utils._normalize_xyz
+   grid.utils._normalize_xyz_scalar
 
 
 Validation
diff --git a/docs/user_api/index.rst b/docs/user_api/index.rst
index 57e40d832..3e5e8b568 100644
--- a/docs/user_api/index.rst
+++ b/docs/user_api/index.rst
@@ -238,6 +238,7 @@ Methods
    Grid.get_kd_tree
    Grid.copy
    Grid.isel
+   Grid.construct_face_centers
    Grid.chunk
 
 
diff --git a/test/test_centroids.py b/test/test_centroids.py
index 6eb69542f..70c5cb208 100644
--- a/test/test_centroids.py
+++ b/test/test_centroids.py
@@ -4,7 +4,8 @@
 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_xyz
+from uxarray.grid.coordinates import _populate_face_centroids, _populate_edge_centroids, _populate_face_centerpoints, _is_inside_circle, _circle_from_three_points, _circle_from_two_points
+from uxarray.grid.coordinates import _normalize_xyz
 
 current_path = Path(os.path.dirname(os.path.realpath(__file__)))
 
@@ -63,7 +64,8 @@ def test_centroids_from_mean_verts_scrip(self):
         expected_face_x = uxgrid.face_lon.values
         expected_face_y = uxgrid.face_lat.values
 
-        _populate_face_centroids(uxgrid, repopulate=True)
+        # _populate_face_centroids(uxgrid, repopulate=True)
+        uxgrid.construct_face_centers(method="cartesian average")
 
         # computed_face_x = (uxgrid.face_lon.values + 180) % 360 - 180
         computed_face_x = uxgrid.face_lon.values
@@ -107,3 +109,60 @@ def test_edge_centroids_from_mpas(self):
 
         nt.assert_array_almost_equal(expected_edge_lon, computed_edge_lon)
         nt.assert_array_almost_equal(expected_edge_lat, computed_edge_lat)
+
+class TestCenterPoints(TestCase):
+
+    def test_circle_from_two_points(self):
+        """Test creation of circle from 2 points."""
+        p1 = (0, 0)
+        p2 = (0, 90)
+        center, radius = _circle_from_two_points(p1, p2)
+
+        # The expected radius in radians should be half the angle between the two vectors
+        expected_center = (0.0, 45.0)
+        expected_radius = np.deg2rad(45.0)
+
+        assert np.allclose(center, expected_center), f"Expected center {expected_center}, but got {center}"
+        assert np.allclose(radius, expected_radius), f"Expected radius {expected_radius}, but got {radius}"
+
+    def test_circle_from_three_points(self):
+        """Test creation of circle from 3 points."""
+        p1 = (0, 0)
+        p2 = (0, 90)
+        p3 = (90, 0)
+        center, radius = _circle_from_three_points(p1, p2, p3)
+        expected_radius = np.deg2rad(45.0)
+        expected_center = (30.0, 30.0)
+
+        assert np.allclose(center, expected_center), f"Expected center {expected_center}, but got {center}"
+        assert np.allclose(radius, expected_radius), f"Expected radius {expected_radius}, but got {radius}"
+
+    def test_is_inside_circle(self):
+        """Test if a points is inside the circle."""
+        # Define the circle
+        circle = ((0.0, 0.0), 1)  # Center at lon/lat with a radius in radians (angular measure of the radius)
+
+        # Define test points
+        point_inside = (30.0, 30.0)  # Should be inside the circle
+        point_outside = (90.0, 0.0)  # Should be outside the circle
+
+        # Test _is_inside_circle function
+        assert _is_inside_circle(circle, point_inside), f"Point {point_inside} should be inside the circle."
+        assert not _is_inside_circle(circle, point_outside), f"Point {point_outside} should be outside the circle."
+
+    def test_face_centerpoint(self):
+        """Use points from an actual spherical face and get the centerpoint."""
+
+        points = np.array([(-35.26438968, -45.0), (-36.61769496, -42.0), (-33.78769181, -42.0), (-32.48416571, -45.0)])
+        uxgrid = ux.open_grid(points, latlon=True)
+
+        # Uses the @property from get face_lon/lat - default is average or centroid
+        ctr_lon = uxgrid.face_lon.values[0]
+        ctr_lat = uxgrid.face_lat.values[0]
+
+        # now explicitly get the centerpoints stored to face_lon/lat using welzl's centerpoint algorithm
+        uxgrid.construct_face_centers(method = "welzl")
+
+        # Test the values of the calculated centerpoint, giving high tolerance of two decimal place
+        nt.assert_array_almost_equal(ctr_lon, uxgrid.face_lon.values[0], decimal=2)
+        nt.assert_array_almost_equal(ctr_lat, uxgrid.face_lat.values[0], decimal=2)
diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py
index f25553321..cc8b84318 100644
--- a/uxarray/grid/coordinates.py
+++ b/uxarray/grid/coordinates.py
@@ -1,16 +1,14 @@
 import xarray as xr
 import numpy as np
-
+import math
 import warnings
 
-from uxarray.constants import ERROR_TOLERANCE
 from uxarray.conventions import ugrid
-
-from typing import Union
+from uxarray.grid.arcs import _angle_of_2_vectors
 
 from numba import njit
-
-import math
+from uxarray.constants import ERROR_TOLERANCE
+from typing import Union
 
 
 @njit(cache=True)
@@ -228,7 +226,7 @@ def _normalize_xyz_scalar(x: float, y: float, z: float):
 
 
 def _populate_node_latlon(grid) -> None:
-    """Populates the latitude and longitude coordinates of a Grid (`node_lon`,
+    """Populates the lon and lat 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
@@ -330,8 +328,26 @@ def _populate_face_centroids(grid, repopulate=False):
 
 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."""
+    Averaging.
 
+    Parameters
+    ----------
+    node_x : numpy.ndarray
+        X coordinates of the nodes.
+    node_y : numpy.ndarray
+        Y coordinates of the nodes.
+    node_z : numpy.ndarray
+        Z coordinates of the nodes.
+    face_nodes : numpy.ndarray
+        Indices of nodes per face.
+    n_nodes_per_face : numpy.ndarray
+        Number of nodes per face.
+
+    Returns
+    -------
+    tuple
+        The x, y, and z coordinates of the centroids.
+    """
     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)
@@ -345,6 +361,270 @@ def _construct_face_centroids(node_x, node_y, node_z, face_nodes, n_nodes_per_fa
     return _normalize_xyz(centroid_x, centroid_y, centroid_z)
 
 
+def _populate_face_centerpoints(grid, repopulate=False):
+    """Calculates the face centerpoints using Welzl's algorithm. It is a
+    randomized algorithm for finding the center and radius of the smallest
+    circle that encloses a set of points. It is here adapted to work on a unit
+    sphere. Also, this algorithm cannot be guaranteed to work on concave
+    polygons.
+
+    Parameters
+    ----------
+    grid : Grid
+        The grid containing the nodes and faces.
+    repopulate : bool, optional
+        Bool used to turn on/off repopulating the face coordinates of the centerpoints, default is False.
+
+    Returns
+    -------
+    None, populates the grid with the face centerpoints: face_lon, face_lat
+    """
+    # warnings.warn("This cannot be guaranteed to work correctly on concave polygons")
+
+    node_lon = grid.node_lon.values
+    node_lat = grid.node_lat.values
+
+    centerpoint_lat = []
+    centerpoint_lon = []
+
+    face_nodes = grid.face_node_connectivity.values
+    n_nodes_per_face = grid.n_nodes_per_face.values
+
+    # Check if the centerpoints are already populated
+    if "face_lon" not in grid._ds or repopulate:
+        centerpoint_lon, centerpoint_lat = _construct_face_centerpoints(
+            node_lon, node_lat, face_nodes, n_nodes_per_face
+        )
+    # get the cartesian coordinates of the centerpoints
+    ctrpt_x, ctrpt_y, ctrpt_z = _lonlat_rad_to_xyz(centerpoint_lon, centerpoint_lat)
+
+    # set the grid variables for centerpoints
+    if "face_lon" not in grid._ds or repopulate:
+        grid._ds["face_lon"] = xr.DataArray(
+            centerpoint_lon, dims=[ugrid.FACE_DIM], attrs=ugrid.FACE_LON_ATTRS
+        )
+        grid._ds["face_lat"] = xr.DataArray(
+            centerpoint_lat, dims=[ugrid.FACE_DIM], attrs=ugrid.FACE_LAT_ATTRS
+        )
+
+    if "face_x" not in grid._ds or repopulate:
+        grid._ds["face_x"] = xr.DataArray(
+            ctrpt_x, dims=[ugrid.FACE_DIM], attrs=ugrid.FACE_X_ATTRS
+        )
+
+        grid._ds["face_y"] = xr.DataArray(
+            ctrpt_y, dims=[ugrid.FACE_DIM], attrs=ugrid.FACE_Y_ATTRS
+        )
+
+        grid._ds["face_z"] = xr.DataArray(
+            ctrpt_z, dims=[ugrid.FACE_DIM], attrs=ugrid.FACE_Z_ATTRS
+        )
+
+
+@njit
+def _circle_from_two_points(p1, p2):
+    """Calculate the smallest circle that encloses two points on a unit sphere.
+
+    Parameters
+    ----------
+    p1 : tuple
+        The first point as a tuple of (latitude, longitude).
+    p2 : tuple
+        The second point as a tuple of (latitude, longitude).
+
+    Returns
+    -------
+    tuple
+        A tuple containing the center (as a tuple of lon and lat) and the radius of the circle.
+    """
+    center_lon = (p1[0] + p2[0]) / 2
+    center_lat = (p1[1] + p2[1]) / 2
+    center = (center_lon, center_lat)
+
+    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)
+    radius = distance / 2
+
+    return center, radius
+
+
+@njit
+def _circle_from_three_points(p1, p2, p3):
+    """Calculate the smallest circle that encloses three points on a unit
+    sphere. This is a placeholder implementation.
+
+    Parameters
+    ----------
+    p1 : tuple
+        The first point.
+    p2 : tuple
+        The second point.
+    p3 : tuple
+        The third point.
+
+    Returns
+    -------
+    tuple
+        A tuple containing the center (as a tuple of lon and lat) and the radius of the circle.
+    """
+    # Placeholder implementation for three-point circle calculation
+    center_lon = (p1[0] + p2[0] + p3[0]) / 3
+    center_lat = (p1[1] + p2[1] + p3[1]) / 3
+    center = (center_lon, center_lat)
+
+    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])))
+    v3 = np.array(_lonlat_rad_to_xyz(np.radians(p3[0]), np.radians(p3[1])))
+
+    radius = (
+        max(
+            _angle_of_2_vectors(v1, v2),
+            _angle_of_2_vectors(v1, v3),
+            _angle_of_2_vectors(v2, v3),
+        )
+        / 2
+    )
+
+    return center, radius
+
+
+@njit
+def _is_inside_circle(circle, point):
+    """Check if a point is inside a given circle on a unit sphere.
+
+    Parameters
+    ----------
+    circle : tuple
+        A tuple containing the center (as a tuple of lon and lat) and the radius of the circle.
+    point : tuple
+        The point to check, as a tuple of (lon, lat).
+
+    Returns
+    -------
+    bool
+        True if the point is inside the circle, False otherwise.
+    """
+    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)
+    return distance <= radius
+
+
+@njit
+def _welzl_recursive(points, boundary, R):
+    """Recursive helper function for Welzl's algorithm to find the smallest
+    enclosing circle.
+
+    Parameters
+    ----------
+    points : numpy.ndarray
+        The set of points to consider.
+    boundary : numpy.ndarray
+        The current boundary points of the minimal enclosing circle.
+    R : tuple
+        The current minimal enclosing circle.
+
+    Returns
+    -------
+    tuple
+        The smallest enclosing circle as a tuple of center and radius.
+    """
+    # Base case: no points or boundary has 3 points
+    if len(points) == 0 or len(boundary) == 3:
+        # Construct the minimal circle based on the number of boundary points
+        if len(boundary) == 0:
+            # Return a default circle if no boundary points are available
+            return ((0.0, 0.0), 0.0)
+        elif len(boundary) == 1:
+            return _circle_from_two_points(boundary[0], boundary[0])
+        elif len(boundary) == 2:
+            return _circle_from_two_points(boundary[0], boundary[1])
+        elif len(boundary) == 3:
+            return _circle_from_three_points(boundary[0], boundary[1], boundary[2])
+
+    p = points[-1]
+    temp_points = points[:-1]
+    circle = _welzl_recursive(temp_points, boundary, R)
+
+    # Check if the chosen point is inside the current circle
+    if circle and _is_inside_circle(circle, p):
+        return circle
+    # If not, the point must be on the boundary of the minimal enclosing circle
+    else:
+        new_boundary = np.empty(
+            (boundary.shape[0] + 1, boundary.shape[1]), dtype=boundary.dtype
+        )
+        new_boundary[:-1] = boundary
+        new_boundary[-1] = p
+        return _welzl_recursive(temp_points, new_boundary, R)
+
+
+def _smallest_enclosing_circle(points):
+    """Find the smallest circle that encloses all given points on a unit sphere
+    using Welzl's algorithm.
+
+    Parameters
+    ----------
+    points : numpy.ndarray
+        An array of points as tuples of (lon, lat).
+
+    Returns
+    -------
+    tuple
+        The smallest enclosing circle as a tuple of center and radius.
+    """
+    np.random.shuffle(
+        points
+    )  # Randomize the input to increase the chance of an optimal solution
+    return _welzl_recursive(points, np.empty((0, 2)), None)
+
+
+def _construct_face_centerpoints(node_lon, node_lat, face_nodes, n_nodes_per_face):
+    """Constructs the face centerpoint using Welzl's algorithm.
+
+    Parameters
+    ----------
+    node_lon : array_like
+        Longitudes of the nodes.
+    node_lat : array_like
+        Latitudes of the nodes.
+    face_nodes : array_like
+        Indices of nodes per face.
+    n_nodes_per_face : array_like
+        Number of nodes per face.
+
+    Returns
+    -------
+    tuple of numpy.ndarray
+        Two arrays containing the longitudes and latitudes of the centerpoints.
+    """
+    num_faces = face_nodes.shape[0]
+    ctrpt_lon = np.zeros(num_faces, dtype=np.float64)
+    ctrpt_lat = np.zeros(num_faces, dtype=np.float64)
+
+    # Pre-compute all points arrays
+    points_arrays = [
+        np.column_stack(
+            (
+                node_lon[face_nodes[face_idx, :n_max_nodes]],
+                node_lat[face_nodes[face_idx, :n_max_nodes]],
+            )
+        )
+        for face_idx, n_max_nodes in enumerate(n_nodes_per_face)
+    ]
+
+    # Compute circles for all faces
+    circles = [_smallest_enclosing_circle(points) for points in points_arrays]
+
+    # Extract centerpoints
+    ctrpt_lon, ctrpt_lat = zip(*[circle[0] for circle in circles])
+
+    return np.array(ctrpt_lon), np.array(ctrpt_lat)
+
+
 def _populate_edge_centroids(grid, repopulate=False):
     """Finds the centroids using cartesian averaging of the edges based off the
     vertices. The centroid is defined as the average of the x, y, z
@@ -432,3 +712,168 @@ def _set_desired_longitude_range(ds):
         if lon_name in ds:
             if ds[lon_name].max() > 180:
                 ds[lon_name].data = (ds[lon_name].data + 180) % 360 - 180
+
+
+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
+    ----------
+    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
+    -------
+    lon : Union[np.ndarray, float]
+        Longitude in radians
+    lat: Union[np.ndarray, float]
+        Latitude in radians
+    """
+
+    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
+
+    lon = np.arctan2(y, x, dtype=np.float64)
+    lat = np.arcsin(z, dtype=np.float64)
+
+    # set longitude range to [0, pi]
+    lon = np.mod(lon, 2 * np.pi)
+
+    z_mask = np.abs(z) > 1.0 - ERROR_TOLERANCE
+
+    lat = np.where(z_mask, np.sign(z) * np.pi / 2, lat)
+    lon = np.where(z_mask, 0.0, lon)
+
+    return lon, lat
+
+
+@njit
+def _xyz_to_lonlat_rad_no_norm(
+    x: Union[np.ndarray, float],
+    y: Union[np.ndarray, float],
+    z: Union[np.ndarray, float],
+):
+    """Converts a Cartesian x,y,z coordinates into Spherical latitude and
+    longitude without normalization, decorated with Numba.
+
+    Parameters
+    ----------
+    x : float
+        Cartesian x coordinate
+    y: float
+        Cartesiain y coordinate
+    z: float
+        Cartesian z coordinate
+
+
+    Returns
+    -------
+    lon : float
+        Longitude in radians
+    lat: float
+        Latitude in radians
+    """
+
+    lon = math.atan2(y, x)
+    lat = math.asin(z)
+
+    # set longitude range to [0, pi]
+    lon = np.mod(lon, 2 * np.pi)
+
+    z_mask = np.abs(z) > 1.0 - ERROR_TOLERANCE
+
+    lat = np.where(z_mask, np.sign(z) * np.pi / 2, lat)
+    lon = np.where(z_mask, 0.0, lon)
+
+    return lon, lat
+
+
+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
+
+
+@njit(cache=True)
+def _lonlat_rad_to_xyz(
+    lon: Union[np.ndarray, float],
+    lat: Union[np.ndarray, float],
+) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
+    """Converts Spherical lon and lat 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)
+
+    return x, y, z
+
+
+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
+    ----------
+    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
+    -------
+    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)
+
+    lon = np.rad2deg(lon_rad)
+    lat = np.rad2deg(lat_rad)
+
+    lon = (lon + 180) % 360 - 180
+    return lon, lat
+
+
+@njit
+def _normalize_xyz_scalar(x: float, y: float, z: float):
+    denom = np.linalg.norm(np.asarray(np.array([x, y, z]), dtype=np.float64), ord=2)
+    x_norm = x / denom
+    y_norm = y / denom
+    z_norm = z / denom
+    return x_norm, y_norm, z_norm
diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py
index 64d6ee653..8fa0ea554 100644
--- a/uxarray/grid/grid.py
+++ b/uxarray/grid/grid.py
@@ -35,6 +35,7 @@
 from uxarray.grid.coordinates import (
     _populate_face_centroids,
     _populate_edge_centroids,
+    _populate_face_centerpoints,
     _set_desired_longitude_range,
     _populate_node_latlon,
     _populate_node_xyz,
@@ -403,6 +404,39 @@ def validate(self):
         else:
             raise RuntimeError("Mesh validation failed.")
 
+    def construct_face_centers(self, method="cartesian average"):
+        """Constructs face centers, this method provides users direct control
+        of the method for constructing the face centers, the default method is
+        "cartesian average", but a more accurate method is "welzl" that is
+        based on the recursive Welzl algorithm. It must be noted that this
+        method can override the parsed/recompute the original parsed face
+        centers.
+
+        Parameters
+        ----------
+        method : str, default="cartesian average"
+            Supported methods are "cartesian average" and "welzl"
+
+        Returns
+        -------
+        None
+            This method constructs the face_lon and face_lat attributes for the grid object.
+
+        Usage
+        -----
+        >>> import uxarray as ux
+        >>> uxgrid = ux.open_grid("GRID_FILE_NAME")
+        >>> face_lat = uxgrid.construct_face_center(method="welzl")
+        """
+        if method == "cartesian average":
+            _populate_face_centroids(self, repopulate=True)
+        elif method == "welzl":
+            _populate_face_centerpoints(self, repopulate=True)
+        else:
+            raise ValueError(
+                f"Unknown method for face center calculation. Expected one of ['cartesian average', 'welzl'] but received {method}"
+            )
+
     def __repr__(self):
         """Constructs a string representation of the contents of a ``Grid``."""
 
diff --git a/uxarray/grid/utils.py b/uxarray/grid/utils.py
index 33f9f75da..9701cef0b 100644
--- a/uxarray/grid/utils.py
+++ b/uxarray/grid/utils.py
@@ -4,70 +4,6 @@
 import uxarray.utils.computing as ac_utils
 
 
-def _replace_fill_values(grid_var, original_fill, new_fill, new_dtype=None):
-    """Replaces all instances of the current fill value (``original_fill``) in
-    (``grid_var``) with (``new_fill``) and converts to the dtype defined by
-    (``new_dtype``)
-
-    Parameters
-    ----------
-    grid_var : np.ndarray
-        grid variable to be modified
-    original_fill : constant
-        original fill value used in (``grid_var``)
-    new_fill : constant
-        new fill value to be used in (``grid_var``)
-    new_dtype : np.dtype, optional
-        new data type to convert (``grid_var``) to
-
-    Returns
-    ----------
-    grid_var : xarray.Dataset
-        Input Dataset with correct fill value and dtype
-    """
-
-    # locations of fill values
-    if original_fill is not None and np.isnan(original_fill):
-        fill_val_idx = np.isnan(grid_var)
-    else:
-        fill_val_idx = grid_var == original_fill
-
-    # convert to new data type
-    if new_dtype != grid_var.dtype and new_dtype is not None:
-        grid_var = grid_var.astype(new_dtype)
-
-    # ensure fill value can be represented with current integer data type
-    if np.issubdtype(new_dtype, np.integer):
-        int_min = np.iinfo(grid_var.dtype).min
-        int_max = np.iinfo(grid_var.dtype).max
-        # ensure new_fill is in range [int_min, int_max]
-        if new_fill < int_min or new_fill > int_max:
-            raise ValueError(
-                f"New fill value: {new_fill} not representable by"
-                f" integer dtype: {grid_var.dtype}"
-            )
-
-    # ensure non-nan fill value can be represented with current float data type
-    elif np.issubdtype(new_dtype, np.floating) and not np.isnan(new_fill):
-        float_min = np.finfo(grid_var.dtype).min
-        float_max = np.finfo(grid_var.dtype).max
-        # ensure new_fill is in range [float_min, float_max]
-        if new_fill < float_min or new_fill > float_max:
-            raise ValueError(
-                f"New fill value: {new_fill} not representable by"
-                f" float dtype: {grid_var.dtype}"
-            )
-    else:
-        raise ValueError(
-            f"Data type {grid_var.dtype} not supported" f"for grid variables"
-        )
-
-    # replace all zeros with a fill value
-    grid_var[fill_val_idx] = new_fill
-
-    return grid_var
-
-
 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.