Repository: UXARRAY/uxarray

Issue:
Implement `Grid.get_faces_containing_point()`
### Proposed new feature or change:

Add a function that finds, for a given point, the faces that contain that point.

Hints from issue comments:

Solution patch:
diff --git a/benchmarks/mpas_ocean.py b/benchmarks/mpas_ocean.py
index aadbd45a1..7f7c78a03 100644
--- a/benchmarks/mpas_ocean.py
+++ b/benchmarks/mpas_ocean.py
@@ -3,6 +3,7 @@
 from pathlib import Path
 
 import numpy as np
+from polars.testing.parametric import dtypes
 
 import uxarray as ux
 
@@ -194,6 +195,31 @@ def time_const_lat(self, resolution, lat_step):
             self.uxgrid.cross_section.constant_latitude(lat)
 
 
+class PointInPolygon:
+    param_names = ['resolution']
+    params = ['480km', '120km']
+
+    def setup(self, resolution):
+        self.uxgrid = ux.open_grid(file_path_dict[resolution][0])
+        self.uxgrid.normalize_cartesian_coordinates()
+
+        # Construct variables needed to ensure that the benchmark doesn't measure construction time
+        _ = self.uxgrid.face_edge_connectivity
+        _ = self.uxgrid.face_x.values
+        _ = self.uxgrid.face_lon.values
+
+        point = np.array([0.0, 0.0, 1.0])
+        res = self.uxgrid.get_faces_containing_point(point)
+
+    def teardown(self, resolution):
+        del self.uxgrid
+
+    def time_face_search(self, resolution):
+        point_xyz = np.array([self.uxgrid.face_x[0].values, self.uxgrid.face_y[0].values, self.uxgrid.face_z[0].values], dtype=np.float64)
+        point_lonlat = np.array([self.uxgrid.face_lon[0].values, self.uxgrid.face_lat.values[0]], dtype=np.float64)
+        self.uxgrid.get_faces_containing_point(point_xyz=point_xyz, point_lonlat=point_lonlat)
+
+
 class ZonalAverage(DatasetBenchmark):
     def setup(self, resolution, *args, **kwargs):
         self.uxds = ux.open_dataset(file_path_dict[resolution][0], file_path_dict[resolution][1])
diff --git a/docs/api.rst b/docs/api.rst
index e5ca434aa..2c9170447 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -139,6 +139,7 @@ Descriptors
    Grid.partial_sphere_coverage
    Grid.global_sphere_coverage
    Grid.triangular
+   Grid.max_face_radius
 
 Attributes
 ~~~~~~~~~~
@@ -159,6 +160,7 @@ Methods
    Grid.calculate_total_face_area
    Grid.normalize_cartesian_coordinates
    Grid.construct_face_centers
+   Grid.get_faces_containing_point
 
 Inheritance of Xarray Functionality
 -----------------------------------
diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py
index 2d78b978a..714de0fce 100644
--- a/uxarray/grid/coordinates.py
+++ b/uxarray/grid/coordinates.py
@@ -100,8 +100,8 @@ def _xyz_to_lonlat_rad(
     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.
+    """Converts Cartesian x, y, z coordinates in Spherical longitude and
+    latitude coordinates in radians.
 
     Parameters
     ----------
@@ -279,7 +279,7 @@ def _populate_face_centroids(grid, repopulate=False):
         # 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 = _lonlat_rad_to_xyz(
-            centroid_lon, centroid_lat
+            np.deg2rad(centroid_lon), np.deg2rad(centroid_lat)
         )
 
     # Populate the centroids
@@ -648,7 +648,7 @@ def _populate_edge_centroids(grid, repopulate=False):
         # 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 = _lonlat_rad_to_xyz(
-            centroid_lon, centroid_lat
+            np.deg2rad(centroid_lon), np.deg2rad(centroid_lat)
         )
 
     # Populate the centroids
diff --git a/uxarray/grid/geometry.py b/uxarray/grid/geometry.py
index f39120fed..1c5660da6 100644
--- a/uxarray/grid/geometry.py
+++ b/uxarray/grid/geometry.py
@@ -18,9 +18,16 @@
     INT_DTYPE,
     INT_FILL_VALUE,
 )
-from uxarray.grid.arcs import extreme_gca_latitude, point_within_gca
+from uxarray.grid.arcs import (
+    extreme_gca_latitude,
+    point_within_gca,
+)
+
 from uxarray.grid.coordinates import _xyz_to_lonlat_rad
-from uxarray.grid.intersections import gca_gca_intersection
+
+from uxarray.grid.intersections import (
+    gca_gca_intersection,
+)
 from uxarray.grid.utils import (
     _get_cartesian_face_edge_nodes,
     _get_lonlat_rad_face_edge_nodes,
@@ -31,6 +38,10 @@
     "North": np.array([0.0, 0.0, 1.0]),
     "South": np.array([0.0, 0.0, -1.0]),
 }
+
+REF_POINT_NORTH_XYZ = np.array([0.01745241, 0.0, 0.9998477], dtype=np.float64)
+REF_POINT_SOUTH_XYZ = np.array([0.01745241, 0.0, -0.9998477], dtype=np.float64)
+
 POLE_POINTS_LONLAT = {
     "North": np.array([0.0, np.pi / 2]),
     "South": np.array([0.0, -np.pi / 2]),
@@ -1561,3 +1572,234 @@ def inverse_stereographic_projection(x, y, central_lon, central_lat):
     )
 
     return lon, lat
+
+
+@njit(cache=True)
+def point_in_face(
+    edges_xyz,
+    point_xyz,
+    inclusive=True,
+):
+    """Determines if a point lies inside a face.
+
+    Parameters
+    ----------
+        edges_xyz : numpy.ndarray
+            Cartesian coordinates of each point in the face
+        point_xyz : numpy.ndarray
+            Cartesian coordinate of the point
+        inclusive : bool
+            Flag to determine whether to include points on the nodes and edges of the face
+
+    Returns
+    -------
+    bool
+        True if point is inside face, False otherwise
+    """
+
+    # Validate the inputs
+    if len(edges_xyz[0][0]) != 3:
+        raise ValueError("`edges_xyz` vertices must be in Cartesian coordinates.")
+
+    if len(point_xyz) != 3:
+        raise ValueError("`point_xyz` must be a single [3] Cartesian coordinate.")
+
+    # Initialize the intersection count
+    intersection_count = 0
+
+    # Set to hold unique intersections
+    unique_intersections = set()
+
+    location = _classify_polygon_location(edges_xyz)
+
+    if location == 1:
+        ref_point_xyz = REF_POINT_SOUTH_XYZ
+    elif location == -1:
+        ref_point_xyz = REF_POINT_NORTH_XYZ
+    else:
+        ref_point_xyz = REF_POINT_SOUTH_XYZ
+
+    # Initialize the points arc between the point and the reference point
+    gca_cart = np.empty((2, 3), dtype=np.float64)
+    gca_cart[0] = point_xyz
+    gca_cart[1] = ref_point_xyz
+
+    # Loop through the face's edges, checking each one for intersection
+    for ind in range(len(edges_xyz)):
+        # If the point lies on an edge, return True if inclusive
+        if point_within_gca(
+            point_xyz,
+            edges_xyz[ind][0],
+            edges_xyz[ind][1],
+        ):
+            if inclusive:
+                return True
+            else:
+                return False
+
+        # Get the number of intersections between the edge and the point arc
+        intersections = gca_gca_intersection(edges_xyz[ind], gca_cart)
+
+        # Add any unique intersections to the intersection_count
+        for intersection in intersections:
+            intersection_tuple = (
+                intersection[0],
+                intersection[1],
+                intersection[2],
+            )
+            if intersection_tuple not in unique_intersections:
+                unique_intersections.add(intersection_tuple)
+                intersection_count += 1
+
+    # Return True if the number of intersections is odd, False otherwise
+    return intersection_count % 2 == 1
+
+
+@njit(cache=True)
+def _find_faces(face_edge_cartesian, point_xyz, inverse_indices):
+    """Finds the faces that contain a given point, inside a subset `face_edge_cartesian`
+    Parameters
+    ----------
+        face_edge_cartesian : numpy.ndarray
+            Cartesian coordinates of all the faces according to their edges
+        point_xyz : numpy.ndarray
+            Cartesian coordinate of the point
+        inverse_indices : numpy.ndarray
+           The original indices of the subsetted grid
+
+    Returns
+    -------
+    index : array
+        The index of the face that contains the point
+    """
+
+    index = []
+
+    # Run for each face in the subset
+    for i, face in enumerate(inverse_indices):
+        # Check to see if the face contains the point
+        contains_point = point_in_face(
+            face_edge_cartesian[i],
+            point_xyz,
+            inclusive=True,
+        )
+
+        # If the point is found, add it to the index array
+        if contains_point:
+            index.append(face)
+
+    # Return the index array
+    return index
+
+
+def _populate_max_face_radius(self):
+    """Populates `max_face_radius`
+
+    Returns
+    -------
+    max_distance : np.float64
+        The max distance from a node to a face center
+    """
+
+    # Parse all variables needed for `njit` functions
+    face_node_connectivity = self.face_node_connectivity.values
+    node_lats_rad = np.deg2rad(self.node_lat.values)
+    node_lons_rad = np.deg2rad(self.node_lon.values)
+    face_lats_rad = np.deg2rad(self.face_lat.values)
+    face_lons_rad = np.deg2rad(self.face_lon.values)
+
+    # Get the max distance
+    max_distance = calculate_max_face_radius(
+        face_node_connectivity,
+        node_lats_rad,
+        node_lons_rad,
+        face_lats_rad,
+        face_lons_rad,
+    )
+
+    # Return the max distance, which is the `max_face_radius`
+    return np.rad2deg(max_distance)
+
+
+@njit(cache=True)
+def calculate_max_face_radius(
+    face_node_connectivity, node_lats_rad, node_lons_rad, face_lats_rad, face_lons_rad
+):
+    """Finds the max face radius in the mesh.
+    Parameters
+    ----------
+        face_node_connectivity : numpy.ndarray
+            Cartesian coordinates of all the faces according to their edges
+        node_lats_rad : numpy.ndarray
+            The `Grid.node_lat` array in radians
+        node_lons_rad : numpy.ndarray
+           The `Grid.node_lon` array in radians
+        face_lats_rad : numpy.ndarray
+           The `Grid.face_lat` array in radians
+        face_lons_rad : numpy.ndarray
+           The `Grid.face_lon` array in radians
+
+    Returns
+    -------
+    The max distance from a node to a face center
+    """
+
+    # Array to store all distances of each face to it's furthest node.
+    end_distances = np.zeros(len(face_node_connectivity))
+
+    # Loop over each face and its nodes
+    for ind, face in enumerate(face_node_connectivity):
+        # Filter out INT_FILL_VALUE
+        valid_nodes = face[face != INT_FILL_VALUE]
+
+        # Get the face lat/lon of this face
+        face_lat = face_lats_rad[ind]
+        face_lon = face_lons_rad[ind]
+
+        # Get the node lat/lon of this face
+        node_lat_rads = node_lats_rad[valid_nodes]
+        node_lon_rads = node_lons_rad[valid_nodes]
+
+        # Calculate Haversine distances for all nodes in this face
+        distances = haversine_distance(node_lon_rads, node_lat_rads, face_lon, face_lat)
+
+        # Store the max distance for this face
+        end_distances[ind] = np.max(distances)
+
+    # Return the maximum distance found across all faces
+    return np.max(end_distances)
+
+
+@njit(cache=True)
+def haversine_distance(lon_a, lat_a, lon_b, lat_b):
+    """Calculates the haversine distance between two points.
+
+    Parameters
+    ----------
+    lon_a : np.float64
+        The longitude of the first point
+    lat_a : np.float64
+        The latitude of the first point
+    lon_b : np.float64
+        The longitude of the second point
+    lat_b : np.float64
+        The latitude of the second point
+
+    Returns
+    -------
+    distance :  np.float64
+        The distance between the two points
+    """
+
+    # Differences in latitudes and longitudes
+    dlat = lat_b - lat_a
+    dlon = lon_b - lon_a
+
+    # Haversine formula
+    equation_in_sqrt = (np.sin(dlat / 2) ** 2) + np.cos(lat_a) * np.cos(lat_b) * (
+        np.sin(dlon / 2) ** 2
+    )
+    distance = 2 * np.arcsin(np.sqrt(equation_in_sqrt))
+
+    # Return the gotten distance
+    return distance
diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py
index f927dc4f9..44bdc37fa 100644
--- a/uxarray/grid/grid.py
+++ b/uxarray/grid/grid.py
@@ -14,6 +14,8 @@
     Tuple,
 )
 
+from uxarray.grid.utils import _get_cartesian_face_edge_nodes
+
 # reader and writer imports
 from uxarray.io._exodus import _read_exodus, _encode_exodus
 from uxarray.io._mpas import _read_mpas
@@ -51,6 +53,8 @@
     _populate_node_xyz,
     _normalize_xyz,
     prepare_points,
+    _lonlat_rad_to_xyz,
+    _xyz_to_lonlat_deg,
 )
 from uxarray.grid.connectivity import (
     _populate_edge_node_connectivity,
@@ -69,6 +73,8 @@
     _populate_bounds,
     _construct_boundary_edge_indices,
     compute_temp_latlon_array,
+    _find_faces,
+    _populate_max_face_radius,
 )
 
 from uxarray.grid.neighbors import (
@@ -118,7 +124,7 @@
 import copy
 
 
-from uxarray.constants import INT_FILL_VALUE
+from uxarray.constants import INT_FILL_VALUE, ERROR_TOLERANCE
 from uxarray.grid.dual import construct_dual
 
 
@@ -1582,6 +1588,13 @@ def is_subset(self):
         """Returns `True` if the Grid is a subset, 'False' otherwise."""
         return self._is_subset
 
+    @property
+    def max_face_radius(self):
+        """Maximum face radius of the grid in degrees"""
+        if "max_face_radius" not in self._ds:
+            self._ds["max_face_radius"] = _populate_max_face_radius(self)
+        return self._ds["max_face_radius"]
+
     def chunk(self, n_node="auto", n_edge="auto", n_face="auto"):
         """Converts all arrays to dask arrays with given chunks across grid
         dimensions in-place.
@@ -2485,3 +2498,118 @@ def get_faces_between_latitudes(self, lats: Tuple[float, float]):
 
         """
         return faces_within_lat_bounds(lats, self.face_bounds_lat.values)
+
+    def get_faces_containing_point(
+        self, point_xyz=None, point_lonlat=None, tolerance=ERROR_TOLERANCE
+    ):
+        """Identifies the indices of faces that contain a given point.
+
+        Parameters
+        ----------
+        point_xyz : numpy.ndarray
+            A point in cartesian coordinates. Best performance if
+        point_lonlat : numpy.ndarray
+            A point in spherical coordinates.
+        tolerance : numpy.ndarray
+            An optional error tolerance for points that lie on the nodes of a face.
+
+        Returns
+        -------
+        index : numpy.ndarray
+            Array of the face indices containing point. Empty if no face is found. This function will typically return
+            a single face, unless the point falls directly on a corner or edge, where there will be multiple values.
+
+        Examples
+        --------
+        Open a grid from a file path:
+
+        >>> import uxarray as ux
+        >>> uxgrid = ux.open_grid("grid_filename.nc")
+
+        Define a spherical point:
+
+        >>> import numpy as np
+        >>> point_lonlat = np.array([45.2, 32.6], dtype=np.float64)
+
+        Define a cartesian point:
+
+        >>> point_xyz = np.array([0.0, 0.0, 1.0], dtype=np.float64)
+
+        Find the indices of the faces that contain the given point:
+
+        >>> lonlat_point_face_indices = uxgrid.get_faces_containing_point(
+        ...     point_lonlat=point_lonlat
+        ... )
+        >>> xyz_point_face_indices = uxgrid.get_faces_containing_point(
+        ...     point_xyz=point_xyz
+        ... )
+
+        """
+        if point_xyz is None and point_lonlat is None:
+            raise ValueError("Either `point_xyz` or `point_lonlat` must be passed in.")
+
+        # Depending on the provided point coordinates, convert to get all needed coordinate systems
+        if point_xyz is None:
+            point_lonlat = np.asarray(point_lonlat, dtype=np.float64)
+            point_xyz = np.array(
+                _lonlat_rad_to_xyz(*np.deg2rad(point_lonlat)), dtype=np.float64
+            )
+        elif point_lonlat is None:
+            point_xyz = np.asarray(point_xyz, dtype=np.float64)
+            point_lonlat = np.array(_xyz_to_lonlat_deg(*point_xyz), dtype=np.float64)
+
+        # Get the maximum face radius of the grid, plus a small adjustment for if the point is this exact radius away
+        max_face_radius = self.max_face_radius.values + 0.0001
+
+        # Try to find a subset in which the point resides
+        try:
+            subset = self.subset.bounding_circle(
+                r=max_face_radius,
+                center_coord=point_lonlat,
+                element="face centers",
+                inverse_indices=True,
+            )
+        # If no subset is found, warn the user
+        except ValueError:
+            # If the grid is partial, let the user know the point likely lies outside the grid region
+            if self.partial_sphere_coverage:
+                warn(
+                    "No faces found. The grid has partial spherical coverage, and the point may be outside the defined region of the grid."
+                )
+            else:
+                warn("No faces found. Try adjusting the tolerance.")
+            return np.empty(0, dtype=np.int64)
+
+        # Get the faces in terms of their edges
+        face_edge_nodes_xyz = _get_cartesian_face_edge_nodes(
+            subset.face_node_connectivity.values,
+            subset.n_face,
+            subset.n_max_face_nodes,
+            subset.node_x.values,
+            subset.node_y.values,
+            subset.node_z.values,
+        )
+
+        # Get the original face indices from the subset
+        inverse_indices = subset.inverse_indices.face.values
+
+        # Check to see if the point is on the nodes of any face
+        lies_on_node = np.isclose(
+            face_edge_nodes_xyz,
+            point_xyz[None, None, :],  # Expands dimensions for broadcasting
+            rtol=tolerance,
+            atol=tolerance,
+        )
+
+        edge_matches = np.all(lies_on_node, axis=-1)
+        face_matches = np.any(edge_matches, axis=1)
+        face_indices = inverse_indices[np.any(face_matches, axis=1)]
+
+        # If a face is in face_indices, return that as the point was found to lie on a node
+        if len(face_indices) != 0:
+            return face_indices
+        else:
+            # Check if any of the faces in the subset contain the point
+            face_indices = _find_faces(face_edge_nodes_xyz, point_xyz, inverse_indices)
+
+            return face_indices


Test patch:
diff --git a/test/test_geometry.py b/test/test_geometry.py
index 7de043479..c066f612a 100644
--- a/test/test_geometry.py
+++ b/test/test_geometry.py
@@ -9,12 +9,15 @@
 import uxarray as ux
 from uxarray.constants import ERROR_TOLERANCE, INT_FILL_VALUE
 import uxarray.utils.computing as ac_utils
-from uxarray.grid.coordinates import _populate_node_latlon, _lonlat_rad_to_xyz, _normalize_xyz, _xyz_to_lonlat_rad
+from uxarray.grid.coordinates import _populate_node_latlon, _lonlat_rad_to_xyz, _normalize_xyz, _xyz_to_lonlat_rad, \
+    _xyz_to_lonlat_deg, _xyz_to_lonlat_rad_scalar
 from uxarray.grid.arcs import extreme_gca_latitude, extreme_gca_z
 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, _pole_point_inside_polygon_cartesian, stereographic_projection, inverse_stereographic_projection
 
+from uxarray.grid.geometry import _populate_face_latlon_bound, _populate_bounds, _pole_point_inside_polygon_cartesian, \
+    stereographic_projection, inverse_stereographic_projection, point_in_face, haversine_distance
 
+from sklearn.metrics.pairwise import haversine_distances
 
 current_path = Path(os.path.dirname(os.path.realpath(__file__)))
 
@@ -31,9 +34,12 @@
 grid_geoflow = current_path / "meshfiles" / "ugrid" / "geoflow-small" / "grid.nc"
 grid_mpas = current_path / "meshfiles" / "mpas" / "QU" / "oQU480.231010.nc"
 
+grid_mpas_2 = current_path / "meshfiles" / "mpas" / "QU" / "mesh.QU.1920km.151026.nc"
+
 # List of grid files to test
 grid_files_latlonBound = [grid_quad_hex, grid_geoflow, gridfile_CSne8, grid_mpas]
 
+
 def test_antimeridian_crossing():
     verts = [[[-170, 40], [180, 30], [165, 25], [-170, 20]]]
 
@@ -44,6 +50,7 @@ def test_antimeridian_crossing():
     assert len(uxgrid.antimeridian_face_indices) == 1
     assert len(gdf['geometry']) == 1
 
+
 def test_antimeridian_point_on():
     verts = [[[-170, 40], [180, 30], [-170, 20]]]
 
@@ -51,10 +58,12 @@ def test_antimeridian_point_on():
 
     assert len(uxgrid.antimeridian_face_indices) == 1
 
+
 def test_linecollection_execution():
     uxgrid = ux.open_grid(gridfile_CSne8)
     lines = uxgrid.to_linecollection()
 
+
 def test_pole_point_inside_polygon_from_vertice_north():
     vertices = [[0.5, 0.5, 0.5], [-0.5, 0.5, 0.5], [-0.5, -0.5, 0.5],
                 [0.5, -0.5, 0.5]]
@@ -74,6 +83,7 @@ def test_pole_point_inside_polygon_from_vertice_north():
     result = _pole_point_inside_polygon_cartesian('South', face_edge_cart)
     assert not result, "South pole should not be inside the polygon"
 
+
 def test_pole_point_inside_polygon_from_vertice_south():
     vertices = [[0.5, 0.5, -0.5], [-0.5, 0.5, -0.5], [0.0, 0.0, -1.0]]
 
@@ -91,6 +101,7 @@ def test_pole_point_inside_polygon_from_vertice_south():
     result = _pole_point_inside_polygon_cartesian('South', face_edge_cart)
     assert result, "South pole should be inside the polygon"
 
+
 def test_pole_point_inside_polygon_from_vertice_pole():
     vertices = [[0, 0, 1], [-0.5, 0.5, 0.5], [-0.5, -0.5, 0.5],
                 [0.5, -0.5, 0.5]]
@@ -110,6 +121,7 @@ def test_pole_point_inside_polygon_from_vertice_pole():
     result = _pole_point_inside_polygon_cartesian('South', face_edge_cart)
     assert not result, "South pole should not be inside the polygon"
 
+
 def test_pole_point_inside_polygon_from_vertice_cross():
     vertices = [[0.6, -0.3, 0.5], [0.2, 0.2, -0.2], [-0.5, 0.1, -0.2],
                 [-0.1, -0.2, 0.2]]
@@ -198,7 +210,7 @@ def _max_latitude_rad_iterative(gca_cart):
             # Update maximum latitude and section if needed
             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[
+                      w1_lonlat[1]) <= ERROR_TOLERANCE or w1_lonlat[
                 1] == max_lat == w2_lonlat[1]:
                 max_section = [w1_new, w2_new]
                 break
@@ -215,6 +227,7 @@ def _max_latitude_rad_iterative(gca_cart):
 
     return np.average([b_lonlat[1], c_lonlat[1]])
 
+
 def _min_latitude_rad_iterative(gca_cart):
     """Calculate the minimum latitude of a great circle arc defined by two
     points.
@@ -285,7 +298,7 @@ def _min_latitude_rad_iterative(gca_cart):
             # Update minimum latitude and section if needed
             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[
+                      w1_lonlat[1]) <= ERROR_TOLERANCE or w1_lonlat[
                 1] == min_lat == w2_lonlat[1]:
                 min_section = [w1_new, w2_new]
                 break
@@ -302,6 +315,7 @@ def _min_latitude_rad_iterative(gca_cart):
 
     return np.average([b_lonlat[1], c_lonlat[1]])
 
+
 def test_extreme_gca_latitude_max():
     gca_cart = np.array([
         _normalize_xyz(*[0.5, 0.5, 0.5]),
@@ -317,6 +331,7 @@ def test_extreme_gca_latitude_max():
     expected_max_latitude = 1.0
     assert np.isclose(max_latitude, expected_max_latitude, atol=ERROR_TOLERANCE)
 
+
 def test_extreme_gca_latitude_max_short():
     # Define a great circle arc in 3D space that has a small span
     gca_cart = np.array([[0.65465367, -0.37796447, -0.65465367], [0.6652466, -0.33896007, -0.6652466]])
@@ -346,6 +361,7 @@ def test_extreme_gca_latitude_min():
     expected_min_latitude = -np.pi / 2
     assert np.isclose(min_latitude, expected_min_latitude, atol=ERROR_TOLERANCE)
 
+
 def test_get_latlonbox_width():
     gca_latlon = np.array([[0.0, 0.0], [0.0, 3.0]])
     width = ux.grid.geometry._get_latlonbox_width(gca_latlon)
@@ -355,6 +371,7 @@ def test_get_latlonbox_width():
     width = ux.grid.geometry._get_latlonbox_width(gca_latlon)
     assert width == 2.0
 
+
 def test_insert_pt_in_latlonbox_non_periodic():
     old_box = np.array([[0.1, 0.2], [0.3, 0.4]])  # Radians
     new_pt = np.array([0.15, 0.35])
@@ -362,6 +379,7 @@ def test_insert_pt_in_latlonbox_non_periodic():
     result = ux.grid.geometry.insert_pt_in_latlonbox(old_box, new_pt, False)
     np.testing.assert_array_equal(result, expected)
 
+
 def test_insert_pt_in_latlonbox_periodic():
     old_box = np.array([[0.1, 0.2], [6.0, 0.1]])  # Radians, periodic
     new_pt = np.array([0.15, 6.2])
@@ -369,6 +387,7 @@ def test_insert_pt_in_latlonbox_periodic():
     result = ux.grid.geometry.insert_pt_in_latlonbox(old_box, new_pt, True)
     np.testing.assert_array_equal(result, expected)
 
+
 def test_insert_pt_in_latlonbox_pole():
     old_box = np.array([[0.1, 0.2], [0.3, 0.4]])
     new_pt = np.array([np.pi / 2, np.nan])  # Pole point
@@ -376,6 +395,7 @@ def test_insert_pt_in_latlonbox_pole():
     result = ux.grid.geometry.insert_pt_in_latlonbox(old_box, new_pt)
     np.testing.assert_array_equal(result, expected)
 
+
 def test_insert_pt_in_empty_state():
     old_box = np.array([[np.nan, np.nan],
                         [np.nan, np.nan]])  # Empty state
@@ -385,7 +405,8 @@ def test_insert_pt_in_empty_state():
     np.testing.assert_array_equal(result, expected)
 
 
-def _get_cartesian_face_edge_nodes_testcase_helper_latlon_bounds_gca(face_nodes_ind, face_edges_ind, edge_nodes_grid, node_x, node_y, node_z):
+def _get_cartesian_face_edge_nodes_testcase_helper_latlon_bounds_gca(face_nodes_ind, face_edges_ind, edge_nodes_grid,
+                                                                     node_x, node_y, node_z):
     """Construct an array to hold the edge Cartesian coordinates connectivity for a face in a grid."""
     mask = face_edges_ind != INT_FILL_VALUE
     valid_edges = face_edges_ind[mask]
@@ -406,7 +427,9 @@ def _get_cartesian_face_edge_nodes_testcase_helper_latlon_bounds_gca(face_nodes_
 
     return cartesian_coordinates
 
-def _get_lonlat_rad_face_edge_nodes_testcase_helper_latlon_bounds_gca(face_nodes_ind, face_edges_ind, edge_nodes_grid, node_lon, node_lat):
+
+def _get_lonlat_rad_face_edge_nodes_testcase_helper_latlon_bounds_gca(face_nodes_ind, face_edges_ind, edge_nodes_grid,
+                                                                      node_lon, node_lat):
     """Construct an array to hold the edge lat lon in radian connectivity for a face in a grid."""
     mask = face_edges_ind != INT_FILL_VALUE
     valid_edges = face_edges_ind[mask]
@@ -433,6 +456,7 @@ def _get_lonlat_rad_face_edge_nodes_testcase_helper_latlon_bounds_gca(face_nodes
 
     return lonlat_coordinates
 
+
 def test_populate_bounds_normal_latlon_bounds_gca():
     vertices_lonlat = [[10.0, 60.0], [10.0, 10.0], [50.0, 10.0], [50.0, 60.0]]
     vertices_lonlat = np.array(vertices_lonlat)
@@ -460,6 +484,7 @@ def test_populate_bounds_normal_latlon_bounds_gca():
     bounds = _populate_face_latlon_bound(face_edges_connectivity_cartesian, face_edges_connectivity_lonlat)
     np.testing.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE)
 
+
 def test_populate_bounds_antimeridian_latlon_bounds_gca():
     vertices_lonlat = [[350, 60.0], [350, 10.0], [50.0, 10.0], [50.0, 60.0]]
     vertices_lonlat = np.array(vertices_lonlat)
@@ -487,6 +512,7 @@ def test_populate_bounds_antimeridian_latlon_bounds_gca():
     bounds = _populate_face_latlon_bound(face_edges_connectivity_cartesian, face_edges_connectivity_lonlat)
     np.testing.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE)
 
+
 def test_populate_bounds_equator_latlon_bounds_gca():
     face_edges_cart = np.array([
         [[0.99726469, -0.05226443, -0.05226443], [0.99862953, 0.0, -0.05233596]],
@@ -501,6 +527,7 @@ def test_populate_bounds_equator_latlon_bounds_gca():
     expected_bounds = np.array([[-0.05235988, 0], [6.23082543, 0]])
     np.testing.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE)
 
+
 def test_populate_bounds_south_sphere_latlon_bounds_gca():
     face_edges_cart = np.array([
         [[-1.04386773e-01, -5.20500333e-02, -9.93173799e-01], [-1.04528463e-01, -1.28010448e-17, -9.94521895e-01]],
@@ -516,6 +543,7 @@ def test_populate_bounds_south_sphere_latlon_bounds_gca():
     expected_bounds = np.array([[-1.51843645, -1.45388627], [3.14159265, 3.92699082]])
     np.testing.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE)
 
+
 def test_populate_bounds_near_pole_latlon_bounds_gca():
     face_edges_cart = np.array([
         [[3.58367950e-01, 0.00000000e+00, -9.33580426e-01], [3.57939780e-01, 4.88684203e-02, -9.32465008e-01]],
@@ -531,6 +559,7 @@ def test_populate_bounds_near_pole_latlon_bounds_gca():
     expected_bounds = np.array([[-1.20427718, -1.14935491], [0, 0.13568803]])
     np.testing.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE)
 
+
 def test_populate_bounds_near_pole2_latlon_bounds_gca():
     face_edges_cart = np.array([
         [[3.57939780e-01, -4.88684203e-02, -9.32465008e-01], [3.58367950e-01, 0.00000000e+00, -9.33580426e-01]],
@@ -546,6 +575,7 @@ def test_populate_bounds_near_pole2_latlon_bounds_gca():
     expected_bounds = np.array([[-1.20427718, -1.14935491], [6.147497, 4.960524e-16]])
     np.testing.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE)
 
+
 def test_populate_bounds_long_face_latlon_bounds_gca():
     face_edges_cart = np.array([
         [[9.9999946355819702e-01, -6.7040475551038980e-04, 8.0396590055897832e-04],
@@ -570,6 +600,7 @@ def test_populate_bounds_long_face_latlon_bounds_gca():
     # The expected bounds should not contain the south pole [0,-0.5*np.pi]
     assert bounds[1][0] != 0.0
 
+
 def test_populate_bounds_node_on_pole_latlon_bounds_gca():
     vertices_lonlat = [[10.0, 90.0], [10.0, 10.0], [50.0, 10.0], [50.0, 60.0]]
     vertices_lonlat = np.array(vertices_lonlat)
@@ -596,6 +627,7 @@ def test_populate_bounds_node_on_pole_latlon_bounds_gca():
     bounds = _populate_face_latlon_bound(face_edges_connectivity_cartesian, face_edges_connectivity_lonlat)
     np.testing.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE)
 
+
 def test_populate_bounds_edge_over_pole_latlon_bounds_gca():
     vertices_lonlat = [[210.0, 80.0], [350.0, 60.0], [10.0, 60.0], [30.0, 80.0]]
     vertices_lonlat = np.array(vertices_lonlat)
@@ -622,6 +654,7 @@ def test_populate_bounds_edge_over_pole_latlon_bounds_gca():
     bounds = _populate_face_latlon_bound(face_edges_connectivity_cartesian, face_edges_connectivity_lonlat)
     np.testing.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE)
 
+
 def test_populate_bounds_pole_inside_latlon_bounds_gca():
     vertices_lonlat = [[200.0, 80.0], [350.0, 60.0], [10.0, 60.0], [40.0, 80.0]]
     vertices_lonlat = np.array(vertices_lonlat)
@@ -649,7 +682,8 @@ def test_populate_bounds_pole_inside_latlon_bounds_gca():
     np.testing.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE)
 
 
-def _get_cartesian_face_edge_nodes_testcase_helper_latlon_bounds_latlonface(face_nodes_ind, face_edges_ind, edge_nodes_grid, node_x, node_y, node_z):
+def _get_cartesian_face_edge_nodes_testcase_helper_latlon_bounds_latlonface(face_nodes_ind, face_edges_ind,
+                                                                            edge_nodes_grid, node_x, node_y, node_z):
     """Construct an array to hold the edge Cartesian coordinates connectivity for a face in a grid."""
     mask = face_edges_ind != INT_FILL_VALUE
     valid_edges = face_edges_ind[mask]
@@ -670,7 +704,9 @@ def _get_cartesian_face_edge_nodes_testcase_helper_latlon_bounds_latlonface(face
 
     return cartesian_coordinates
 
-def _get_lonlat_rad_face_edge_nodes_testcase_helper_latlon_bounds_latlonface(face_nodes_ind, face_edges_ind, edge_nodes_grid, node_lon, node_lat):
+
+def _get_lonlat_rad_face_edge_nodes_testcase_helper_latlon_bounds_latlonface(face_nodes_ind, face_edges_ind,
+                                                                             edge_nodes_grid, node_lon, node_lat):
     """Construct an array to hold the edge lat lon in radian connectivity for a face in a grid."""
     mask = face_edges_ind != INT_FILL_VALUE
     valid_edges = face_edges_ind[mask]
@@ -697,6 +733,7 @@ def _get_lonlat_rad_face_edge_nodes_testcase_helper_latlon_bounds_latlonface(fac
 
     return lonlat_coordinates
 
+
 def test_populate_bounds_normal_latlon_bounds_latlonface():
     vertices_lonlat = [[10.0, 60.0], [10.0, 10.0], [50.0, 10.0], [50.0, 60.0]]
     vertices_lonlat = np.array(vertices_lonlat)
@@ -718,9 +755,11 @@ def test_populate_bounds_normal_latlon_bounds_latlonface():
         grid.edge_node_connectivity.values, grid.node_lon.values,
         grid.node_lat.values)
     expected_bounds = np.array([[lat_min, lat_max], [lon_min, lon_max]])
-    bounds = _populate_face_latlon_bound(face_edges_connectivity_cartesian, face_edges_connectivity_lonlat, is_latlonface=True)
+    bounds = _populate_face_latlon_bound(face_edges_connectivity_cartesian, face_edges_connectivity_lonlat,
+                                         is_latlonface=True)
     np.testing.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE)
 
+
 def test_populate_bounds_antimeridian_latlon_bounds_latlonface():
     vertices_lonlat = [[350, 60.0], [350, 10.0], [50.0, 10.0], [50.0, 60.0]]
     vertices_lonlat = np.array(vertices_lonlat)
@@ -742,9 +781,11 @@ def test_populate_bounds_antimeridian_latlon_bounds_latlonface():
         grid.edge_node_connectivity.values, grid.node_lon.values,
         grid.node_lat.values)
     expected_bounds = np.array([[lat_min, lat_max], [lon_min, lon_max]])
-    bounds = _populate_face_latlon_bound(face_edges_connectivity_cartesian, face_edges_connectivity_lonlat, is_latlonface=True)
+    bounds = _populate_face_latlon_bound(face_edges_connectivity_cartesian, face_edges_connectivity_lonlat,
+                                         is_latlonface=True)
     np.testing.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE)
 
+
 def test_populate_bounds_node_on_pole_latlon_bounds_latlonface():
     vertices_lonlat = [[10.0, 90.0], [10.0, 10.0], [50.0, 10.0], [50.0, 60.0]]
     vertices_lonlat = np.array(vertices_lonlat)
@@ -766,9 +807,11 @@ def test_populate_bounds_node_on_pole_latlon_bounds_latlonface():
         grid.edge_node_connectivity.values, grid.node_lon.values,
         grid.node_lat.values)
     expected_bounds = np.array([[lat_min, lat_max], [lon_min, lon_max]])
-    bounds = _populate_face_latlon_bound(face_edges_connectivity_cartesian, face_edges_connectivity_lonlat, is_latlonface=True)
+    bounds = _populate_face_latlon_bound(face_edges_connectivity_cartesian, face_edges_connectivity_lonlat,
+                                         is_latlonface=True)
     np.testing.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE)
 
+
 def test_populate_bounds_edge_over_pole_latlon_bounds_latlonface():
     vertices_lonlat = [[210.0, 80.0], [350.0, 60.0], [10.0, 60.0], [30.0, 80.0]]
     vertices_lonlat = np.array(vertices_lonlat)
@@ -790,9 +833,11 @@ def test_populate_bounds_edge_over_pole_latlon_bounds_latlonface():
         grid.edge_node_connectivity.values, grid.node_lon.values,
         grid.node_lat.values)
     expected_bounds = np.array([[lat_min, lat_max], [lon_min, lon_max]])
-    bounds = _populate_face_latlon_bound(face_edges_connectivity_cartesian, face_edges_connectivity_lonlat, is_latlonface=True)
+    bounds = _populate_face_latlon_bound(face_edges_connectivity_cartesian, face_edges_connectivity_lonlat,
+                                         is_latlonface=True)
     np.testing.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE)
 
+
 def test_populate_bounds_pole_inside_latlon_bounds_latlonface():
     vertices_lonlat = [[200.0, 80.0], [350.0, 60.0], [10.0, 60.0], [40.0, 80.0]]
     vertices_lonlat = np.array(vertices_lonlat)
@@ -814,11 +859,13 @@ def test_populate_bounds_pole_inside_latlon_bounds_latlonface():
         grid.edge_node_connectivity.values, grid.node_lon.values,
         grid.node_lat.values)
     expected_bounds = np.array([[lat_min, lat_max], [lon_min, lon_max]])
-    bounds = _populate_face_latlon_bound(face_edges_connectivity_cartesian, face_edges_connectivity_lonlat, is_latlonface=True)
+    bounds = _populate_face_latlon_bound(face_edges_connectivity_cartesian, face_edges_connectivity_lonlat,
+                                         is_latlonface=True)
     np.testing.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE)
 
 
-def _get_cartesian_face_edge_nodes_testcase_helper_latlon_bounds_gca_list(face_nodes_ind, face_edges_ind, edge_nodes_grid, node_x, node_y, node_z):
+def _get_cartesian_face_edge_nodes_testcase_helper_latlon_bounds_gca_list(face_nodes_ind, face_edges_ind,
+                                                                          edge_nodes_grid, node_x, node_y, node_z):
     """Construct an array to hold the edge Cartesian coordinates connectivity for a face in a grid."""
     mask = face_edges_ind != INT_FILL_VALUE
     valid_edges = face_edges_ind[mask]
@@ -839,7 +886,9 @@ def _get_cartesian_face_edge_nodes_testcase_helper_latlon_bounds_gca_list(face_n
 
     return cartesian_coordinates
 
-def _get_lonlat_rad_face_edge_nodes_testcase_helper_latlon_bounds_gca_list(face_nodes_ind, face_edges_ind, edge_nodes_grid, node_lon, node_lat):
+
+def _get_lonlat_rad_face_edge_nodes_testcase_helper_latlon_bounds_gca_list(face_nodes_ind, face_edges_ind,
+                                                                           edge_nodes_grid, node_lon, node_lat):
     """Construct an array to hold the edge lat lon in radian connectivity for a face in a grid."""
     mask = face_edges_ind != INT_FILL_VALUE
     valid_edges = face_edges_ind[mask]
@@ -866,6 +915,7 @@ def _get_lonlat_rad_face_edge_nodes_testcase_helper_latlon_bounds_gca_list(face_
 
     return lonlat_coordinates
 
+
 def test_populate_bounds_normal_latlon_bounds_gca_list():
     vertices_lonlat = [[10.0, 60.0], [10.0, 10.0], [50.0, 10.0], [50.0, 60.0]]
     vertices_lonlat = np.array(vertices_lonlat)
@@ -891,6 +941,7 @@ def test_populate_bounds_normal_latlon_bounds_gca_list():
                                          is_GCA_list=[True, False, True, False])
     np.testing.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE)
 
+
 def test_populate_bounds_antimeridian_latlon_bounds_gca_list():
     vertices_lonlat = [[350, 60.0], [350, 10.0], [50.0, 10.0], [50.0, 60.0]]
     vertices_lonlat = np.array(vertices_lonlat)
@@ -916,6 +967,7 @@ def test_populate_bounds_antimeridian_latlon_bounds_gca_list():
                                          is_GCA_list=[True, False, True, False])
     np.testing.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE)
 
+
 def test_populate_bounds_node_on_pole_latlon_bounds_gca_list():
     vertices_lonlat = [[10.0, 90.0], [10.0, 10.0], [50.0, 10.0], [50.0, 60.0]]
     vertices_lonlat = np.array(vertices_lonlat)
@@ -941,6 +993,7 @@ def test_populate_bounds_node_on_pole_latlon_bounds_gca_list():
                                          is_GCA_list=[True, False, True, False])
     np.testing.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE)
 
+
 def test_populate_bounds_edge_over_pole_latlon_bounds_gca_list():
     vertices_lonlat = [[210.0, 80.0], [350.0, 60.0], [10.0, 60.0], [30.0, 80.0]]
     vertices_lonlat = np.array(vertices_lonlat)
@@ -966,6 +1019,7 @@ def test_populate_bounds_edge_over_pole_latlon_bounds_gca_list():
                                          is_GCA_list=[True, False, True, False])
     np.testing.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE)
 
+
 def test_populate_bounds_pole_inside_latlon_bounds_gca_list():
     vertices_lonlat = [[200.0, 80.0], [350.0, 60.0], [10.0, 60.0], [40.0, 80.0]]
     vertices_lonlat = np.array(vertices_lonlat)
@@ -1010,6 +1064,7 @@ def test_populate_bounds_GCA_mix_latlon_bounds_mix():
     for i in range(len(faces)):
         np.testing.assert_allclose(face_bounds[i], expected_bounds[i], atol=ERROR_TOLERANCE)
 
+
 def test_populate_bounds_LatlonFace_mix_latlon_bounds_mix():
     face_1 = [[10.0, 60.0], [10.0, 10.0], [50.0, 10.0], [50.0, 60.0]]
     face_2 = [[350, 60.0], [350, 10.0], [50.0, 10.0], [50.0, 60.0]]
@@ -1029,6 +1084,7 @@ def test_populate_bounds_LatlonFace_mix_latlon_bounds_mix():
     for i in range(len(faces)):
         np.testing.assert_allclose(face_bounds[i], expected_bounds[i], atol=ERROR_TOLERANCE)
 
+
 def test_populate_bounds_GCAList_mix_latlon_bounds_mix():
     face_1 = [[10.0, 60.0], [10.0, 10.0], [50.0, 10.0], [50.0, 60.0]]
     face_2 = [[350, 60.0], [350, 10.0], [50.0, 10.0], [50.0, 60.0]]
@@ -1051,6 +1107,7 @@ def test_populate_bounds_GCAList_mix_latlon_bounds_mix():
     for i in range(len(faces)):
         np.testing.assert_allclose(face_bounds[i], expected_bounds[i], atol=ERROR_TOLERANCE)
 
+
 def test_face_bounds_latlon_bounds_files():
     """Test to ensure ``Grid.face_bounds`` works correctly for all grid files."""
     for grid_path in grid_files_latlonBound:
@@ -1071,16 +1128,19 @@ def test_face_bounds_latlon_bounds_files():
             # Clean up the grid object
             del uxgrid
 
+
 def test_engine_geodataframe():
     uxgrid = ux.open_grid(gridfile_geoflow)
     for engine in ['geopandas', 'spatialpandas']:
         gdf = uxgrid.to_geodataframe(engine=engine)
 
+
 def test_periodic_elements_geodataframe():
     uxgrid = ux.open_grid(gridfile_geoflow)
     for periodic_elements in ['ignore', 'exclude', 'split']:
         gdf = uxgrid.to_geodataframe(periodic_elements=periodic_elements)
 
+
 def test_to_gdf_geodataframe():
     uxgrid = ux.open_grid(gridfile_geoflow)
 
@@ -1088,6 +1148,7 @@ def test_to_gdf_geodataframe():
 
     gdf_without_am = uxgrid.to_geodataframe(exclude_antimeridian=True)
 
+
 def test_cache_and_override_geodataframe():
     """Tests the cache and override functionality for GeoDataFrame conversion."""
     uxgrid = ux.open_grid(gridfile_geoflow)
@@ -1114,6 +1175,200 @@ def test_cache_and_override_geodataframe():
 
     assert gdf_f is not gdf_e
 
+
+# Test point_in_face function
+def test_point_inside():
+    """Test the function `point_in_face`, where the points are all inside the face"""
+
+    # Open grid
+    grid = ux.open_grid(grid_mpas_2)
+
+    # Get the face edges of all faces in the grid
+    faces_edges_cartesian = _get_cartesian_face_edge_nodes(
+        grid.face_node_connectivity.values,
+        grid.n_face,
+        grid.n_max_face_edges,
+        grid.node_x.values,
+        grid.node_y.values,
+        grid.node_z.values,
+    )
+
+    # Loop through each face
+    for i in range(grid.n_face):
+        # Set the point as the face center of the polygon
+        point_xyz = np.array([grid.face_x[i].values, grid.face_y[i].values, grid.face_z[i].values])
+
+        # Assert that the point is in the polygon
+        assert point_in_face(faces_edges_cartesian[i], point_xyz, inclusive=True)
+
+
+def test_point_outside():
+    """Test the function `point_in_face`, where the point is outside the face"""
+
+    # Open grid
+    grid = ux.open_grid(grid_mpas_2)
+
+    # Get the face edges of all faces in the grid
+    faces_edges_cartesian = _get_cartesian_face_edge_nodes(
+        grid.face_node_connectivity.values,
+        grid.n_face,
+        grid.n_max_face_edges,
+        grid.node_x.values,
+        grid.node_y.values,
+        grid.node_z.values,
+    )
+
+    # Set the point as the face center of a different face than the face tested
+    point_xyz = np.array([grid.face_x[1].values, grid.face_y[1].values, grid.face_z[1].values])
+
+    # Assert that the point is not in the face tested
+    assert not point_in_face(faces_edges_cartesian[0], point_xyz, inclusive=True)
+
+
+def test_point_on_node():
+    """Test the function `point_in_face`, when the point is on the node of the polygon"""
+
+    # Open grid
+    grid = ux.open_grid(grid_mpas_2)
+
+    # Get the face edges of all faces in the grid
+    faces_edges_cartesian = _get_cartesian_face_edge_nodes(
+        grid.face_node_connectivity.values,
+        grid.n_face,
+        grid.n_max_face_edges,
+        grid.node_x.values,
+        grid.node_y.values,
+        grid.node_z.values,
+    )
+
+    # Set the point as a node
+    point_xyz = np.array([*faces_edges_cartesian[0][0][0]])
+
+    # Assert that the point is in the face when inclusive is true
+    assert point_in_face(faces_edges_cartesian[0], point_xyz, inclusive=True)
+
+    # Assert that the point is not in the face when inclusive is false
+    assert not point_in_face(faces_edges_cartesian[0], point_xyz, inclusive=False)
+
+
+def test_point_inside_close():
+    """Test the function `point_in_face`, where the point is inside the face, but very close to the edge"""
+
+    # Create a square
+    vertices_lonlat = [[-10.0, 10.0], [-10.0, -10.0], [10.0, -10.0], [10.0, 10.0]]
+    vertices_lonlat = np.array(vertices_lonlat)
+
+    # Choose a point just inside the square
+    point = np.array(_lonlat_rad_to_xyz(np.deg2rad(0.0), np.deg2rad(-9.8)))
+
+    # Create the grid and face edges
+    grid = ux.Grid.from_face_vertices(vertices_lonlat, latlon=True)
+    faces_edges_cartesian = _get_cartesian_face_edge_nodes(
+        grid.face_node_connectivity.values,
+        grid.n_face,
+        grid.n_max_face_edges,
+        grid.node_x.values,
+        grid.node_y.values,
+        grid.node_z.values,
+    )
+
+    # Use point in face to determine if the point is inside or out of the face
+    assert point_in_face(faces_edges_cartesian[0], point_xyz=point, inclusive=False)
+
+
+def test_point_outside_close():
+    """Test the function `point_in_face`, where the point is outside the face, but very close to the edge"""
+
+    # Create a square
+    vertices_lonlat = [[-10.0, 10.0], [-10.0, -10.0], [10.0, -10.0], [10.0, 10.0]]
+    vertices_lonlat = np.array(vertices_lonlat)
+
+    # Choose a point just outside the square
+    point = np.array(_lonlat_rad_to_xyz(np.deg2rad(0.0), np.deg2rad(-10.2)))
+
+    # Create the grid and face edges
+    grid = ux.Grid.from_face_vertices(vertices_lonlat, latlon=True)
+    faces_edges_cartesian = _get_cartesian_face_edge_nodes(
+        grid.face_node_connectivity.values,
+        grid.n_face,
+        grid.n_max_face_edges,
+        grid.node_x.values,
+        grid.node_y.values,
+        grid.node_z.values,
+    )
+
+    # Use point in face to determine if the point is inside or out of the face
+    assert not point_in_face(faces_edges_cartesian[0], point_xyz=point, inclusive=False)
+
+
+def test_face_at_pole():
+    """Test the function `point_in_face`, when the face is at the North Pole"""
+
+    # Generate a face that is at a pole
+    vertices_lonlat = [[10.0, 90.0], [10.0, 10.0], [50.0, 10.0], [50.0, 60.0]]
+    vertices_lonlat = np.array(vertices_lonlat)
+
+    point = np.array(_lonlat_rad_to_xyz(np.deg2rad(25), np.deg2rad(30)))
+
+    # Create the grid and face edges
+    grid = ux.Grid.from_face_vertices(vertices_lonlat, latlon=True)
+    faces_edges_cartesian = _get_cartesian_face_edge_nodes(
+        grid.face_node_connectivity.values,
+        grid.n_face,
+        grid.n_max_face_edges,
+        grid.node_x.values,
+        grid.node_y.values,
+        grid.node_z.values,
+    )
+
+    assert point_in_face(faces_edges_cartesian[0], point_xyz=point, inclusive=True)
+
+
+def test_face_at_antimeridian():
+    """Test the function `point_in_face`, where the face crosses the antimeridian"""
+
+    # Generate a face crossing the antimeridian
+    vertices_lonlat = [[350, 60.0], [350, 10.0], [50.0, 10.0], [50.0, 60.0]]
+    vertices_lonlat = np.array(vertices_lonlat)
+    point = np.array(_lonlat_rad_to_xyz(np.deg2rad(25), np.deg2rad(30)))
+
+    # Create the grid and face edges
+    grid = ux.Grid.from_face_vertices(vertices_lonlat, latlon=True)
+    faces_edges_cartesian = _get_cartesian_face_edge_nodes(
+        grid.face_node_connectivity.values,
+        grid.n_face,
+        grid.n_max_face_edges,
+        grid.node_x.values,
+        grid.node_y.values,
+        grid.node_z.values,
+    )
+
+    assert point_in_face(faces_edges_cartesian[0], point_xyz=point, inclusive=True)
+
+
+def test_face_normal_face():
+    """Test the function `point_in_face`, where the face is a normal face, not crossing the antimeridian or the
+    poles"""
+
+    # Generate a normal face that is not crossing the antimeridian or the poles
+    vertices_lonlat = [[10.0, 60.0], [10.0, 10.0], [50.0, 10.0], [50.0, 60.0]]
+    vertices_lonlat = np.array(vertices_lonlat)
+    point = np.array(_lonlat_rad_to_xyz(np.deg2rad(25), np.deg2rad(30)))
+
+    # Create the grid and face edges
+    grid = ux.Grid.from_face_vertices(vertices_lonlat, latlon=True)
+    faces_edges_cartesian = _get_cartesian_face_edge_nodes(
+        grid.face_node_connectivity.values,
+        grid.n_face,
+        grid.n_max_face_edges,
+        grid.node_x.values,
+        grid.node_y.values,
+        grid.node_z.values,
+    )
+
+    assert point_in_face(faces_edges_cartesian[0], point_xyz=point, inclusive=True)
+
+
 def test_stereographic_projection_stereographic_projection():
     lon = np.array(0)
     lat = np.array(0)
@@ -1128,3 +1383,17 @@ def test_stereographic_projection_stereographic_projection():
     assert np.array_equal(lon, new_lon)
     assert np.array_equal(lat, new_lat)
     assert np.array_equal(x, y) and x == 0
+
+
+def test_haversine_distance_creation():
+    """Tests the use of `haversine_distance`"""
+
+    # Create two points
+    point_a = [np.deg2rad(-34.8), np.deg2rad(-58.5)]
+    point_b = [np.deg2rad(49.0), np.deg2rad(2.6)]
+
+    result = haversine_distances([point_a, point_b])
+
+    distance = haversine_distance(point_a[1], point_a[0], point_b[1], point_b[0])
+
+    assert np.isclose(result[0][1], distance, atol=1e-6)
diff --git a/test/test_grid.py b/test/test_grid.py
index 4c8d98c76..df9b9dddd 100644
--- a/test/test_grid.py
+++ b/test/test_grid.py
@@ -13,13 +13,14 @@
 from uxarray.grid.connectivity import _populate_face_edge_connectivity, _build_edge_face_connectivity, \
     _build_edge_node_connectivity, _build_face_face_connectivity, _populate_face_face_connectivity
 
-from uxarray.grid.coordinates import _populate_node_latlon, _lonlat_rad_to_xyz
+from uxarray.grid.coordinates import _populate_node_latlon, _lonlat_rad_to_xyz, _xyz_to_lonlat_rad_scalar
 
 from uxarray.constants import INT_FILL_VALUE, ERROR_TOLERANCE
 
 from uxarray.grid.arcs import extreme_gca_latitude
 
 from uxarray.grid.validation import _find_duplicate_nodes
+from .test_gradient import quad_hex_grid_path
 
 try:
     import constants
@@ -34,8 +35,6 @@
 gridfile_CSne30 = current_path / "meshfiles" / "ugrid" / "outCSne30" / "outCSne30.ug"
 gridfile_fesom = current_path / "meshfiles" / "ugrid" / "fesom" / "fesom.mesh.diag.nc"
 gridfile_geoflow = current_path / "meshfiles" / "ugrid" / "geoflow-small" / "grid.nc"
-gridfile_mpas = current_path / 'meshfiles' / "mpas" / "QU" / 'mesh.QU.1920km.151026.nc'
-gridfile_mpas_two = current_path / 'meshfiles' / "mpas" / "QU" / 'oQU480.231010.nc'
 gridfile_geos = current_path / 'meshfiles' / "geos-cs" / "c12" / 'test-c12.native.nc4'
 gridfile_mpas_holes = current_path / 'meshfiles' / "mpas" / "QU" / 'oQU480.231010.nc'
 
@@ -44,13 +43,10 @@
 
 shp_filename = current_path / "meshfiles" / "shp" / "grid_fire.shp"
 
-
-
 grid_CSne30 = ux.open_grid(gridfile_CSne30)
 grid_RLL1deg = ux.open_grid(gridfile_RLL1deg)
 grid_RLL10deg_CSne4 = ux.open_grid(gridfile_RLL10deg_CSne4)
 
-
 mpas_filepath = current_path / "meshfiles" / "mpas" / "QU" / "mesh.QU.1920km.151026.nc"
 exodus_filepath = current_path / "meshfiles" / "exodus" / "outCSne8" / "outCSne8.g"
 ugrid_filepath_01 = current_path / "meshfiles" / "ugrid" / "outCSne30" / "outCSne30.ug"
@@ -80,6 +76,7 @@
 
 gridfile_ugrid = current_path / "meshfiles" / "ugrid" / "geoflow-small" / "grid.nc"
 gridfile_mpas = current_path / "meshfiles" / "mpas" / "QU" / "mesh.QU.1920km.151026.nc"
+gridfile_mpas_two = current_path / 'meshfiles' / "mpas" / "QU" / 'oQU480.231010.nc'
 gridfile_exodus = current_path / "meshfiles" / "exodus" / "outCSne8" / "outCSne8.g"
 gridfile_scrip = current_path / "meshfiles" / "scrip" / "outCSne8" / "outCSne8.nc"
 
@@ -89,6 +86,7 @@ def test_grid_validate():
     grid_mpas = ux.open_grid(gridfile_mpas)
     assert grid_mpas.validate()
 
+
 def test_grid_with_holes():
     """Test _holes_in_mesh function."""
     grid_without_holes = ux.open_grid(gridfile_mpas)
@@ -97,6 +95,7 @@ def test_grid_with_holes():
     assert grid_with_holes.partial_sphere_coverage
     assert grid_without_holes.global_sphere_coverage
 
+
 def test_grid_encode_as():
     """Reads a ugrid file and encodes it as `xarray.Dataset` in various types."""
     grid_CSne30.encode_as("UGRID")
@@ -107,6 +106,7 @@ def test_grid_encode_as():
     grid_RLL1deg.encode_as("Exodus")
     grid_RLL10deg_CSne4.encode_as("Exodus")
 
+
 def test_grid_init_verts():
     """Create a uxarray grid from multiple face vertices with duplicate nodes and saves a ugrid file."""
     cart_x = [
@@ -131,7 +131,7 @@ def test_grid_init_verts():
         [5, 4, 7, 6],  # back face
         [4, 0, 3, 7],  # left face
         [3, 2, 6, 7],  # top face
-        [4, 5, 1, 0]   # bottom face
+        [4, 5, 1, 0]  # bottom face
     ]
 
     faces_coords = []
@@ -163,6 +163,7 @@ def test_grid_init_verts():
     assert vgrid.n_node == 6
     vgrid.encode_as("UGRID")
 
+
 def test_grid_init_verts_different_input_datatype():
     """Create a uxarray grid from multiple face vertices with different datatypes (ndarray, list, tuple) and saves a ugrid file."""
     faces_verts_ndarray = np.array([
@@ -176,8 +177,8 @@ def test_grid_init_verts_different_input_datatype():
     vgrid.encode_as("UGRID")
 
     faces_verts_list = [[[150, 10], [160, 20], [150, 30], [135, 30], [125, 20], [135, 10]],
-                         [[125, 20], [135, 30], [125, 60], [110, 60], [100, 30], [105, 20]],
-                         [[95, 10], [105, 20], [100, 30], [85, 30], [75, 20], [85, 10]]]
+                        [[125, 20], [135, 30], [125, 60], [110, 60], [100, 30], [105, 20]],
+                        [[95, 10], [105, 20], [100, 30], [85, 30], [75, 20], [85, 10]]]
     vgrid = ux.open_grid(faces_verts_list, latlon=True)
     assert vgrid.n_face == 3
     assert vgrid.n_node == 14
@@ -195,6 +196,7 @@ def test_grid_init_verts_different_input_datatype():
     assert vgrid.validate()
     vgrid.encode_as("UGRID")
 
+
 def test_grid_init_verts_fill_values():
     faces_verts_filled_values = [[[150, 10], [160, 20], [150, 30],
                                   [135, 30], [125, 20], [135, 10]],
@@ -208,6 +210,7 @@ def test_grid_init_verts_fill_values():
     assert vgrid.n_face == 3
     assert vgrid.n_node == 12
 
+
 def test_grid_properties():
     """Tests to see if accessing variables through set properties is equal to using the dict."""
     xr.testing.assert_equal(grid_CSne30.node_lon, grid_CSne30._ds["node_lon"])
@@ -234,27 +237,32 @@ def test_grid_properties():
     assert n_faces == grid_geoflow.n_face
     assert n_face_nodes == grid_geoflow.n_max_face_nodes
 
+
 def test_read_shpfile():
     """Reads a shape file and write ugrid file."""
     with pytest.raises(ValueError):
         grid_shp = ux.open_grid(shp_filename)
 
+
 def test_read_scrip():
     """Reads a scrip file."""
     grid_CSne8 = ux.open_grid(gridfile_CSne8)  # tests from scrip
 
+
 def test_operators_eq():
     """Test Equals ('==') operator."""
     grid_CSne30_01 = ux.open_grid(gridfile_CSne30)
     grid_CSne30_02 = ux.open_grid(gridfile_CSne30)
     assert grid_CSne30_01 == grid_CSne30_02
 
+
 def test_operators_ne():
     """Test Not Equals ('!=') operator."""
     grid_CSne30_01 = ux.open_grid(gridfile_CSne30)
     grid_RLL1deg = ux.open_grid(gridfile_RLL1deg)
     assert grid_CSne30_01 != grid_RLL1deg
 
+
 def test_face_areas_calculate_total_face_area_triangle():
     """Create a uxarray grid from vertices and saves an exodus file."""
     verts = [[[0.57735027, -5.77350269e-01, -0.57735027],
@@ -275,11 +283,13 @@ def test_face_areas_calculate_total_face_area_triangle():
         quadrature_rule="triangular", order=4)
     nt.assert_almost_equal(area_triangular, constants.TRI_AREA, decimal=1)
 
+
 def test_face_areas_calculate_total_face_area_file():
     """Create a uxarray grid from vertices and saves an exodus file."""
     area = ux.open_grid(gridfile_CSne30).calculate_total_face_area()
     nt.assert_almost_equal(area, constants.MESH30_AREA, decimal=3)
 
+
 def test_face_areas_calculate_total_face_area_sphere():
     """Computes the total face area of an MPAS mesh that lies on a unit sphere, with an expected total face area of 4pi."""
     mpas_grid_path = current_path / 'meshfiles' / "mpas" / "QU" / 'mesh.QU.1920km.151026.nc'
@@ -293,11 +303,13 @@ def test_face_areas_calculate_total_face_area_sphere():
     nt.assert_almost_equal(primal_face_area, constants.UNIT_SPHERE_AREA, decimal=3)
     nt.assert_almost_equal(dual_face_area, constants.UNIT_SPHERE_AREA, decimal=3)
 
+
 def test_face_areas_compute_face_areas_geoflow_small():
     """Checks if the GeoFlow Small can generate a face areas output."""
     grid_geoflow = ux.open_grid(gridfile_geoflow)
     grid_geoflow.compute_face_areas()
 
+
 def test_face_areas_verts_calc_area():
     faces_verts_ndarray = np.array([
         np.array([[150, 10, 0], [160, 20, 0], [150, 30, 0], [135, 30, 0],
@@ -311,11 +323,12 @@ def test_face_areas_verts_calc_area():
     face_verts_areas = verts_grid.face_areas
     nt.assert_almost_equal(face_verts_areas.sum(), constants.FACE_VERTS_AREA, decimal=3)
 
+
 def test_populate_coordinates_populate_cartesian_xyz_coord():
     # The following testcases are generated through the matlab cart2sph/sph2cart functions
     lon_deg = [
         45.0001052295749, 45.0001052295749, 360 - 45.0001052295749,
-        360 - 45.0001052295749
+                                            360 - 45.0001052295749
     ]
     lat_deg = [
         35.2655522903022, -35.2655522903022, 35.2655522903022,
@@ -342,10 +355,11 @@ def test_populate_coordinates_populate_cartesian_xyz_coord():
         nt.assert_almost_equal(vgrid.node_y.values[i], cart_y[i], decimal=12)
         nt.assert_almost_equal(vgrid.node_z.values[i], cart_z[i], decimal=12)
 
+
 def test_populate_coordinates_populate_lonlat_coord():
     lon_deg = [
         45.0001052295749, 45.0001052295749, 360 - 45.0001052295749,
-        360 - 45.0001052295749
+                                            360 - 45.0001052295749
     ]
     lat_deg = [
         35.2655522903022, -35.2655522903022, 35.2655522903022,
@@ -374,8 +388,8 @@ def test_populate_coordinates_populate_lonlat_coord():
 
 
 def _revert_edges_conn_to_face_nodes_conn(edge_nodes_connectivity: np.ndarray,
-                                           face_edges_connectivity: np.ndarray,
-                                           original_face_nodes_connectivity: np.ndarray):
+                                          face_edges_connectivity: np.ndarray,
+                                          original_face_nodes_connectivity: np.ndarray):
     """Utilize the edge_nodes_connectivity and face_edges_connectivity to
     generate the res_face_nodes_connectivity in the counter-clockwise
     order. The counter-clockwise order will be enforced by the passed in
@@ -440,6 +454,7 @@ def _revert_edges_conn_to_face_nodes_conn(edge_nodes_connectivity: np.ndarray,
 
     return np.array(res_face_nodes_connectivity)
 
+
 def test_connectivity_build_n_nodes_per_face():
     """Tests the construction of the ``n_nodes_per_face`` variable."""
     grids = [grid_mpas, grid_exodus, grid_ugrid]
@@ -457,6 +472,7 @@ def test_connectivity_build_n_nodes_per_face():
     expected_nodes_per_face = np.array([6, 3, 4, 6, 6, 4, 4], dtype=int)
     nt.assert_equal(grid_from_verts.n_nodes_per_face.values, expected_nodes_per_face)
 
+
 def test_connectivity_edge_nodes_euler():
     """Verifies that (``n_edge``) follows euler's formula."""
     grid_paths = [exodus_filepath, ugrid_filepath_01, ugrid_filepath_02, ugrid_filepath_03]
@@ -470,6 +486,7 @@ def test_connectivity_edge_nodes_euler():
 
         assert (n_face == n_edge - n_node + 2)
 
+
 def test_connectivity_build_face_edges_connectivity_mpas():
     """Tests the construction of (``Mesh2_edge_nodes``) on an MPAS grid with known edge nodes."""
     from uxarray.grid.connectivity import _build_edge_node_connectivity
@@ -492,6 +509,7 @@ def test_connectivity_build_face_edges_connectivity_mpas():
 
     assert (n_face == n_edge - n_node + 2)
 
+
 def test_connectivity_build_face_edges_connectivity():
     """Generates Grid.Mesh2_edge_nodes from Grid.face_node_connectivity."""
     ug_filename_list = [ugrid_filepath_01, ugrid_filepath_02, ugrid_filepath_03]
@@ -522,6 +540,7 @@ def test_connectivity_build_face_edges_connectivity():
         for i in range(len(reverted_mesh2_edge_nodes)):
             assert np.array_equal(reverted_mesh2_edge_nodes[i], original_face_nodes_connectivity[i])
 
+
 def test_connectivity_build_face_edges_connectivity_fillvalues():
     verts = [f0_deg, f1_deg, f2_deg, f3_deg, f4_deg, f5_deg, f6_deg]
     uds = ux.open_grid(verts)
@@ -543,6 +562,7 @@ def test_connectivity_build_face_edges_connectivity_fillvalues():
 
     assert np.array_equal(res_face_nodes_connectivity, uds._ds["face_node_connectivity"].values)
 
+
 def test_connectivity_node_face_connectivity_from_verts():
     """Test generating Grid.Mesh2_node_faces from array input."""
     face_nodes_conn_lonlat_degree = [[162., 30], [216., 30], [70., 30],
@@ -573,6 +593,7 @@ def test_connectivity_node_face_connectivity_from_verts():
 
     assert np.array_equal(vgrid.node_face_connectivity.values, expected)
 
+
 def test_connectivity_node_face_connectivity_from_files():
     """Test generating Grid.Mesh2_node_faces from file input."""
     grid_paths = [exodus_filepath, ugrid_filepath_01, ugrid_filepath_02, ugrid_filepath_03]
@@ -593,12 +614,14 @@ def test_connectivity_node_face_connectivity_from_files():
 
         for i in range(grid_ux.n_node):
             face_index_from_sparse_matrix = grid_ux.node_face_connectivity.values[i]
-            valid_face_index_from_sparse_matrix = face_index_from_sparse_matrix[face_index_from_sparse_matrix != grid_ux.node_face_connectivity.attrs["_FillValue"]]
+            valid_face_index_from_sparse_matrix = face_index_from_sparse_matrix[
+                face_index_from_sparse_matrix != grid_ux.node_face_connectivity.attrs["_FillValue"]]
             valid_face_index_from_sparse_matrix.sort()
             face_index_from_dict = node_face_connectivity[i]
             face_index_from_dict.sort()
             assert np.array_equal(valid_face_index_from_sparse_matrix, face_index_from_dict)
 
+
 def test_connectivity_edge_face_connectivity_mpas():
     """Tests the construction of ``Mesh2_face_edges`` to the expected results of an MPAS grid."""
     uxgrid = ux.open_grid(mpas_filepath)
@@ -611,6 +634,7 @@ def test_connectivity_edge_face_connectivity_mpas():
 
     nt.assert_array_equal(edge_faces_output, edge_faces_gold)
 
+
 def test_connectivity_edge_face_connectivity_sample():
     """Tests the construction of ``Mesh2_face_edges`` on an example with one shared edge, and the remaining edges only being part of one face."""
     verts = [[(0.0, -90.0), (180, 0.0), (0.0, 90)],
@@ -633,6 +657,7 @@ def test_connectivity_edge_face_connectivity_sample():
     assert n_solo == uxgrid.n_edge - n_shared
     assert n_invalid == 0
 
+
 def test_connectivity_face_face_connectivity_construction():
     """Tests the construction of face-face connectivity."""
     grid = ux.open_grid(mpas_filepath)
@@ -645,6 +670,7 @@ def test_connectivity_face_face_connectivity_construction():
 
     nt.assert_array_equal(face_face_conn_new_sorted, face_face_conn_old_sorted)
 
+
 def test_class_methods_from_dataset():
     # UGRID
     xrds = xr.open_dataset(gridfile_ugrid)
@@ -663,6 +689,7 @@ def test_class_methods_from_dataset():
     xrds = xr.open_dataset(gridfile_scrip)
     uxgrid = ux.Grid.from_dataset(xrds)
 
+
 def test_class_methods_from_face_vertices():
     single_face_latlon = [(0.0, 90.0), (-180, 0.0), (0.0, -90)]
     uxgrid = ux.Grid.from_face_vertices(single_face_latlon, latlon=True)
@@ -673,6 +700,7 @@ def test_class_methods_from_face_vertices():
 
     single_face_cart = [(0.0,)]
 
+
 def test_latlon_bounds_populate_bounds_GCA_mix():
     gridfile_mpas = current_path / "meshfiles" / "mpas" / "QU" / "oQU480.231010.nc"
     face_1 = [[10.0, 60.0], [10.0, 10.0], [50.0, 10.0], [50.0, 60.0]]
@@ -691,11 +719,13 @@ def test_latlon_bounds_populate_bounds_GCA_mix():
     bounds_xarray = grid.bounds
     nt.assert_allclose(bounds_xarray.values, expected_bounds, atol=ERROR_TOLERANCE)
 
+
 def test_latlon_bounds_populate_bounds_MPAS():
     gridfile_mpas = current_path / "meshfiles" / "mpas" / "QU" / "oQU480.231010.nc"
     uxgrid = ux.open_grid(gridfile_mpas)
     bounds_xarray = uxgrid.bounds
 
+
 def test_dual_mesh_mpas():
     grid = ux.open_grid(gridfile_mpas, use_dual=False)
     mpas_dual = ux.open_grid(gridfile_mpas, use_dual=True)
@@ -708,11 +738,13 @@ def test_dual_mesh_mpas():
 
     nt.assert_equal(dual.face_node_connectivity.values, mpas_dual.face_node_connectivity.values)
 
+
 def test_dual_duplicate():
     dataset = ux.open_dataset(gridfile_geoflow, gridfile_geoflow)
     with pytest.raises(RuntimeError):
         dataset.get_dual()
 
+
 def test_normalize_existing_coordinates_non_norm_initial():
     gridfile_mpas = current_path / "meshfiles" / "mpas" / "QU" / "mesh.QU.1920km.151026.nc"
     from uxarray.grid.validation import _check_normalization
@@ -726,9 +758,82 @@ def test_normalize_existing_coordinates_non_norm_initial():
     uxgrid.normalize_cartesian_coordinates()
     assert _check_normalization(uxgrid)
 
+
 def test_normalize_existing_coordinates_norm_initial():
     gridfile_CSne30 = current_path / "meshfiles" / "ugrid" / "outCSne30" / "outCSne30.ug"
     from uxarray.grid.validation import _check_normalization
     uxgrid = ux.open_grid(gridfile_CSne30)
 
     assert _check_normalization(uxgrid)
+
+
+def test_number_of_faces_found():
+    """Test function for `self.get_face_containing_point`,
+    to ensure the correct number of faces is found, depending on where the point is."""
+    grid = ux.open_grid(gridfile_mpas)
+    partial_grid = ux.open_grid(quad_hex_grid_path)
+
+    # For a face center only one face should be found
+    point_xyz = np.array([grid.face_x[100].values, grid.face_y[100].values, grid.face_z[100].values], dtype=np.float64)
+
+    assert len(grid.get_faces_containing_point(point_xyz=point_xyz)) == 1
+
+    # For an edge two faces should be found
+    point_xyz = np.array([grid.edge_x[100].values, grid.edge_y[100].values, grid.edge_z[100].values], dtype=np.float64)
+
+    assert len(grid.get_faces_containing_point(point_xyz=point_xyz)) == 2
+
+    # For a node three faces should be found
+    point_xyz = np.array([grid.node_x[100].values, grid.node_y[100].values, grid.node_z[100].values], dtype=np.float64)
+
+    assert len(grid.get_faces_containing_point(point_xyz=point_xyz)) == 3
+
+    partial_grid.normalize_cartesian_coordinates()
+
+    # Test for a node on the edge where only 2 faces should be found
+    point_xyz = np.array([partial_grid.node_x[1].values, partial_grid.node_y[1].values, partial_grid.node_z[1].values], dtype=np.float64)
+
+    assert len(partial_grid.get_faces_containing_point(point_xyz=point_xyz)) == 2
+
+
+def test_whole_grid():
+    """Tests `self.get_faces_containing_point`on an entire grid,
+    checking that for each face center, one face is found to contain it"""
+
+    grid = ux.open_grid(gridfile_mpas_two)
+    grid.normalize_cartesian_coordinates()
+
+    # Ensure a face is found on the grid for every face center
+    for i in range(len(grid.face_x.values)):
+        point_xyz = np.array([grid.face_x[i].values, grid.face_y[i].values, grid.face_z[i].values], dtype=np.float64)
+
+        assert len(grid.get_faces_containing_point(point_xyz=point_xyz)) == 1
+
+def test_point_types():
+    """Tests that `self.get_faces_containing_point` works with cartesian and lonlat"""
+
+    # Open the grid
+    grid = ux.open_grid(gridfile_mpas)
+
+    # Assign a cartesian point and a lon/lat point
+    point_xyz = np.array([grid.node_x[100].values, grid.node_y[100].values, grid.node_z[100].values], dtype=np.float64)
+    point_lonlat = np.array([grid.node_lon[100].values, grid.node_lat[100].values])
+
+    # Test both points find faces
+    assert len(grid.get_faces_containing_point(point_xyz=point_xyz)) != 0
+    assert len(grid.get_faces_containing_point(point_lonlat=point_lonlat)) !=0
+
+def test_point_along_arc():
+    node_lon = np.array([-40, -40, 40, 40])
+    node_lat = np.array([-20, 20, 20, -20])
+    face_node_connectivity = np.array([[0, 1, 2, 3]], dtype=np.int64)
+
+    uxgrid = ux.Grid.from_topology(node_lon, node_lat, face_node_connectivity)
+
+    # point at exactly 20 degrees latitude
+    out1 = uxgrid.get_faces_containing_point(point_lonlat=np.array([0, 20], dtype=np.float64))
+
+    # point at 25.41 degrees latitude (max along the great circle arc)
+    out2 = uxgrid.get_faces_containing_point(point_lonlat=np.array([0, 25.41], dtype=np.float64))
+
+    nt.assert_array_equal(out1, out2)