From a6b96ebb5839d8c72ac54d919805a657db8dba62 Mon Sep 17 00:00:00 2001
From: Sam Evans <47793072+Sevans711@users.noreply.github.com>
Date: Fri, 2 Oct 2026 09:46:05 -0400
Subject: [PATCH 01/10] avoid hit_buf tiny array

Checked an example raster is equivalent (via casting to xr.DataArray then checking .equals) before and after this commit.
---
 uxarray/grid/point_in_face.py | 30 +++++++++++++++++-------------
 1 file changed, 17 insertions(+), 13 deletions(-)

diff --git a/uxarray/grid/point_in_face.py b/uxarray/grid/point_in_face.py
index ed6faf7fd..d949ee4ac 100644
--- a/uxarray/grid/point_in_face.py
+++ b/uxarray/grid/point_in_face.py
@@ -77,7 +77,9 @@ def _face_contains_point(face_edges: np.ndarray, point: np.ndarray) -> bool:
 
 
 @njit(cache=True)
-def _get_faces_containing_point(
+def _set_faces_containing_point(
+    result: np.ndarray,
+    i: int,
     point: np.ndarray,
     candidate_indices: np.ndarray,
     face_node_connectivity: np.ndarray,
@@ -85,12 +87,17 @@ def _get_faces_containing_point(
     node_x: np.ndarray,
     node_y: np.ndarray,
     node_z: np.ndarray,
-) -> np.ndarray:
+) -> int:
     """
-    Test each candidate face to see if it contains the query point.
+    Test each candidate face to see if it contains the query point,
+    setting result[i, j] = candidate_indices[k] for all k where the
+    point is inside the face. j starts at 0 and increments by 1 for
+    each hit. Returns the total number of hits.
 
     Parameters
     ----------
+    result : np.ndarray, shape (n_points, max_candidates)
+        Preallocated array to store face indices for each point.
     point : np.ndarray, shape (3,)
         Cartesian unit-vector of the query point.
     candidate_indices : np.ndarray, shape (k,)
@@ -104,10 +111,9 @@ def _get_faces_containing_point(
 
     Returns
     -------
-    hits : np.ndarray, shape (h,)
-        Subset of `candidate_indices` for which the point is inside the face.
+    n_hits : int
+        Number of candidate faces that contain the point.
     """
-    hit_buf = np.empty(candidate_indices.shape[0], dtype=INT_DTYPE)
     count = 0
     for k in range(candidate_indices.shape[0]):
         fidx = candidate_indices[k]
@@ -115,9 +121,9 @@ def _get_faces_containing_point(
             fidx, face_node_connectivity, n_nodes_per_face, node_x, node_y, node_z
         )
         if _face_contains_point(face_edges, point):
-            hit_buf[count] = fidx
+            result[i, count] = fidx
             count += 1
-    return hit_buf[:count]
+    return count
 
 
 @njit(cache=True, parallel=True, nogil=True)
@@ -165,12 +171,10 @@ def _batch_point_in_face(
         p = points[i]
         cands = flat_candidate_indices[start:end]
 
-        hits = _get_faces_containing_point(
-            p, cands, face_node_connectivity, n_nodes_per_face, node_x, node_y, node_z
+        n_hits = _set_faces_containing_point(
+            results, i, p, cands, face_node_connectivity, n_nodes_per_face, node_x, node_y, node_z
         )
-        for j, fi in enumerate(hits):
-            results[i, j] = fi
-        counts[i] = hits.shape[0]
+        counts[i] = n_hits
 
     return results, counts
 

From 594739b74758cecb83c14f67e4a089efa5035b28 Mon Sep 17 00:00:00 2001
From: Sam Evans <47793072+Sevans711@users.noreply.github.com>
Date: Fri, 2 Oct 2026 10:42:19 -0400
Subject: [PATCH 02/10] add _point_in_face (from nodes, not edges)

renames old _face_contains_point to _face_contains_point_from_edges, for now. Next, delete it entirely, and rewrite call sites appropriately to use the new _point_in_face syntax instead.

Confirmed that a raster is exactly the same, before and after these changes. Only a small speedup so far (~4.70s to ~4.45s), but there are still more optimizations to do!
---
 test/grid/geometry/test_geometry.py      |   8 +-
 test/grid/geometry/test_point_in_face.py |  12 +--
 uxarray/grid/geometry.py                 |   4 +-
 uxarray/grid/point_in_face.py            | 110 +++++++++++++++++++++--
 uxarray/utils/numba_math.py              |  10 +++
 5 files changed, 126 insertions(+), 18 deletions(-)

diff --git a/test/grid/geometry/test_geometry.py b/test/grid/geometry/test_geometry.py
index e602e1e0d..ec9607919 100644
--- a/test/grid/geometry/test_geometry.py
+++ b/test/grid/geometry/test_geometry.py
@@ -6,7 +6,7 @@
 from uxarray.constants import ERROR_TOLERANCE, INT_FILL_VALUE
 from uxarray.grid.coordinates import _lonlat_rad_to_xyz, _normalize_xyz
 from uxarray.grid.geometry import haversine_distance, _pole_point_inside_polygon_cartesian
-from uxarray.grid.point_in_face import _face_contains_point
+from uxarray.grid.point_in_face import _face_contains_point_from_edges
 from uxarray.grid.utils import _get_cartesian_face_edge_nodes_array
 
 
@@ -132,7 +132,7 @@ def test_face_at_antimeridian():
         grid.node_z.values,
     )
 
-    assert _face_contains_point(faces_edges_cartesian[0], point)
+    assert _face_contains_point_from_edges(faces_edges_cartesian[0], point)
 
 
 def test_face_at_pole():
@@ -155,7 +155,7 @@ def test_face_at_pole():
         grid.node_z.values,
     )
 
-    assert _face_contains_point(faces_edges_cartesian[0], point)
+    assert _face_contains_point_from_edges(faces_edges_cartesian[0], point)
 
 
 def test_face_normal_face():
@@ -178,7 +178,7 @@ def test_face_normal_face():
         grid.node_z.values,
     )
 
-    assert _face_contains_point(faces_edges_cartesian[0], point)
+    assert _face_contains_point_from_edges(faces_edges_cartesian[0], point)
 
 
 def test_haversine_distance_creation():
diff --git a/test/grid/geometry/test_point_in_face.py b/test/grid/geometry/test_point_in_face.py
index 2a0cafc87..8bf7b0f89 100644
--- a/test/grid/geometry/test_point_in_face.py
+++ b/test/grid/geometry/test_point_in_face.py
@@ -6,7 +6,7 @@
 from uxarray.constants import INT_FILL_VALUE
 from uxarray.grid.coordinates import _lonlat_rad_to_xyz
 from uxarray.grid.utils import _get_cartesian_face_edge_nodes_array, _get_cartesian_face_edge_nodes
-from uxarray.grid.point_in_face import _face_contains_point
+from uxarray.grid.point_in_face import _face_contains_point_from_edges
 
 
 def test_point_inside(gridpath):
@@ -25,7 +25,7 @@ def test_point_inside(gridpath):
         # 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 _face_contains_point(face_edges, point_xyz)
+        assert _face_contains_point_from_edges(face_edges, point_xyz)
 
 
 def test_point_outside(gridpath):
@@ -48,7 +48,7 @@ def test_point_outside(gridpath):
     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 _face_contains_point(faces_edges_cartesian[0], point_xyz)
+    assert not _face_contains_point_from_edges(faces_edges_cartesian[0], point_xyz)
 
 
 def test_point_on_node(gridpath):
@@ -71,7 +71,7 @@ def test_point_on_node(gridpath):
     point_xyz = np.array([*faces_edges_cartesian[0][0][0]])
 
     # Assert that the point is in the face when inclusive is true
-    assert _face_contains_point(faces_edges_cartesian[0], point_xyz)
+    assert _face_contains_point_from_edges(faces_edges_cartesian[0], point_xyz)
 
 
 def test_point_inside_close():
@@ -96,7 +96,7 @@ def test_point_inside_close():
     )
 
     # Use point in face to determine if the point is inside or out of the face
-    assert _face_contains_point(faces_edges_cartesian[0], point)
+    assert _face_contains_point_from_edges(faces_edges_cartesian[0], point)
 
 
 def test_point_outside_close():
@@ -121,4 +121,4 @@ def test_point_outside_close():
     )
 
     # Use point in face to determine if the point is inside or out of the face
-    assert not _face_contains_point(faces_edges_cartesian[0], point)
+    assert not _face_contains_point_from_edges(faces_edges_cartesian[0], point)
diff --git a/uxarray/grid/geometry.py b/uxarray/grid/geometry.py
index 31481956d..0d850b846 100644
--- a/uxarray/grid/geometry.py
+++ b/uxarray/grid/geometry.py
@@ -13,7 +13,7 @@
 from uxarray.grid.intersections import (
     gca_gca_intersection,
 )
-from uxarray.grid.point_in_face import _face_contains_point
+from uxarray.grid.point_in_face import _face_contains_point_from_edges
 from uxarray.grid.utils import _get_cartesian_face_edge_nodes
 from uxarray.utils.imports import _raise_hint_if_optional_deps_missing
 
@@ -1271,7 +1271,7 @@ def barycentric_coordinates_cartesian(polygon_xyz, point_xyz):
             )
 
             # Check to see if the point lies within the current triangle
-            contains_point = _face_contains_point(
+            contains_point = _face_contains_point_from_edges(
                 face_edge,
                 point_xyz,
             )
diff --git a/uxarray/grid/point_in_face.py b/uxarray/grid/point_in_face.py
index d949ee4ac..d84b109d5 100644
--- a/uxarray/grid/point_in_face.py
+++ b/uxarray/grid/point_in_face.py
@@ -8,6 +8,13 @@
 from uxarray.constants import ERROR_TOLERANCE, INT_DTYPE, INT_FILL_VALUE
 from uxarray.grid.arcs import point_within_gca
 from uxarray.grid.utils import _get_cartesian_face_edge_nodes, _small_angle_of_2_vectors
+from uxarray.utils.numba_math import (
+    _numba_allclose3,
+    _numba_cross3,
+    _numba_dot3,
+    _numba_norm3,
+    _numba_sub3,
+)
 
 if TYPE_CHECKING:
     from numpy.typing import ArrayLike
@@ -16,7 +23,7 @@
 
 
 @njit(cache=True)
-def _face_contains_point(face_edges: np.ndarray, point: np.ndarray) -> bool:
+def _face_contains_point_from_edges(face_edges: np.ndarray, point: np.ndarray) -> bool:
     """
     Determine whether a point lies within a face using the spherical winding-number method.
 
@@ -76,6 +83,95 @@ def _face_contains_point(face_edges: np.ndarray, point: np.ndarray) -> bool:
     return np.abs(total) > np.pi
 
 
+@njit(cache=True)
+def _point_in_face(
+    point: np.ndarray,
+    nodes_idx: np.ndarray,
+    node_x: np.ndarray,
+    node_y: np.ndarray,
+    node_z: np.ndarray,
+) -> bool:
+    """Returns whether this point lies within the face formed by these nodes.
+
+    Uses the spherical winding-number method, which
+    sums the signed central angles between successive vertices of the face
+    as seen from `point`.  If the total absolute winding exceeds π, the point is inside.
+    Points exactly on a node or edge also count as inside.
+
+    Parameters
+    ----------
+    point : np.ndarray, shape (3,)
+        3D unit-vector of the query point on the unit sphere.
+    nodes_idx : np.ndarray, shape (n_nodes,)
+        Node indices (within node_x, node_y, node_z) for precisely all nodes in this face.
+        Likely from `face_node_connectivity[fidx][:n_nodes_per_face[fidx]]`.
+    node_x, node_y, node_z : np.ndarray, shape (n_nodes,)
+        Cartesian coordinates of all nodes.
+        (This method uses the values at indices indicated by nodes_idx.)
+
+    Returns
+    -------
+    inside : bool
+        True if the point is inside the face or lies exactly on a node/edge; False otherwise.
+    """
+    # Rewritten from _face_contains_point_from_edges to avoids creating tiny numpy arrays.
+    # Creating tiny numpy arrays from scratch inside numba is very inefficient.
+    # (Creating tiny numpy arrays from indexing larger arrays is fine, though.)
+    # This is the main reason to provide inputs as nodes_idx and node_x, ..., instead of
+    # simply asking to provide a single array of x, y, z coordinates for all nodes;
+    # the former avoids any need to create a tiny numpy array to store the coordinates.
+
+    n_nodes = len(nodes_idx)
+    max_i_node = n_nodes - 1
+
+    # Check for an exact hit with any of the corner nodes
+    for i in range(n_nodes):
+        node_idx = nodes_idx[i]
+        node_xyz = (node_x[node_idx], node_y[node_idx], node_z[node_idx])
+        if _numba_allclose3(node_xyz, point, rtol=ERROR_TOLERANCE, atol=ERROR_TOLERANCE):
+            return True
+
+    # Check whether point lies on any edge of the face
+    # (edges are great-circle arcs between successive nodes)
+    for i in range(n_nodes):
+        # edge is formed by nodes (a, b)
+        ai = nodes_idx[i]
+        bi = nodes_idx[i + 1] if i < max_i_node else nodes_idx[0]
+
+        # TODO: avoid tiny numpy arrays, after rewriting point_within_gca to accept tuples
+        a = np.array([node_x[ai], node_y[ai], node_z[ai]])
+        b = np.array([node_x[bi], node_y[bi], node_z[bi]])
+        if point_within_gca(point, a, b):
+            return True
+
+    # Apply spherical winding-number method:
+    total = 0.0
+    for i in range(n_nodes):
+        # edge is formed by nodes (a, b).
+        ai = nodes_idx[i]
+        bi = nodes_idx[i + 1] if i < max_i_node else nodes_idx[0]
+
+        a = (node_x[ai], node_y[ai], node_z[ai])
+        b = (node_x[bi], node_y[bi], node_z[bi])
+
+        vi = _numba_sub3(a, point)
+        vj = _numba_sub3(b, point)
+
+        # check if you’re right on a vertex
+        if _numba_norm3(vi) < ERROR_TOLERANCE or _numba_norm3(vj) < ERROR_TOLERANCE:
+            return True
+
+        ang = _small_angle_of_2_vectors(vi, vj)
+
+        # determine sign from cross
+        c = _numba_cross3(vi, vj)
+        sign = 1.0 if _numba_dot3(c, point) >= 0.0 else -1.0
+
+        total += sign * ang
+
+    return np.abs(total) > np.pi
+
+
 @njit(cache=True)
 def _set_faces_containing_point(
     result: np.ndarray,
@@ -96,8 +192,12 @@ def _set_faces_containing_point(
 
     Parameters
     ----------
-    result : np.ndarray, shape (n_points, max_candidates)
+    result : np.ndarray, shape (n_points, max_possible_hits)
         Preallocated array to store face indices for each point.
+        max_possible_hits can be much less than max_candidates;
+        the maximum number of hits is n_max_face_nodes, because
+        the "worst case" of point being a node would lead to hits
+        of all faces it is a part of, but nothing else.
     point : np.ndarray, shape (3,)
         Cartesian unit-vector of the query point.
     candidate_indices : np.ndarray, shape (k,)
@@ -117,10 +217,8 @@ def _set_faces_containing_point(
     count = 0
     for k in range(candidate_indices.shape[0]):
         fidx = candidate_indices[k]
-        face_edges = _get_cartesian_face_edge_nodes(
-            fidx, face_node_connectivity, n_nodes_per_face, node_x, node_y, node_z
-        )
-        if _face_contains_point(face_edges, point):
+        nodes_idx = face_node_connectivity[fidx][:n_nodes_per_face[fidx]]
+        if _point_in_face(point, nodes_idx, node_x, node_y, node_z):
             result[i, count] = fidx
             count += 1
     return count
diff --git a/uxarray/utils/numba_math.py b/uxarray/utils/numba_math.py
index c2f0ad673..7f75f2da2 100644
--- a/uxarray/utils/numba_math.py
+++ b/uxarray/utils/numba_math.py
@@ -124,3 +124,13 @@ def _numba_allfinite3(u):
         int(math.isfinite(u[0])) * int(math.isfinite(u[1])) * int(math.isfinite(u[2]))
     )
     # (use `*` instead of `and` to avoid branching logic)
+
+
+@njit(cache=True)
+def _numba_allclose3(u, v, rtol=1e-05, atol=1e-08):
+    """Return (as 1 or 0) whether all components of two 3-vectors are close to each other."""
+    return (
+        int(np.isclose(u[0], v[0], rtol=rtol, atol=atol))
+        * int(np.isclose(u[1], v[1], rtol=rtol, atol=atol))
+        * int(np.isclose(u[2], v[2], rtol=rtol, atol=atol))
+    )

From 5eaaa4154e98d699b3d0a0fd695f60b6b0dfb508 Mon Sep 17 00:00:00 2001
From: Sam Evans <47793072+Sevans711@users.noreply.github.com>
Date: Fri, 2 Oct 2026 11:12:14 -0400
Subject: [PATCH 03/10] add _point_in_face_from_grid; rmv ..._from_edges

Something in test_cross_sections.py is breaking now (numba python types compiler errors) but all other tests pass.
---
 test/grid/geometry/test_geometry.py      | 44 ++---------------
 test/grid/geometry/test_point_in_face.py | 63 ++++--------------------
 uxarray/grid/geometry.py                 | 20 +++-----
 uxarray/grid/point_in_face.py            | 13 +++++
 4 files changed, 33 insertions(+), 107 deletions(-)

diff --git a/test/grid/geometry/test_geometry.py b/test/grid/geometry/test_geometry.py
index ec9607919..9fa83ca7e 100644
--- a/test/grid/geometry/test_geometry.py
+++ b/test/grid/geometry/test_geometry.py
@@ -6,7 +6,7 @@
 from uxarray.constants import ERROR_TOLERANCE, INT_FILL_VALUE
 from uxarray.grid.coordinates import _lonlat_rad_to_xyz, _normalize_xyz
 from uxarray.grid.geometry import haversine_distance, _pole_point_inside_polygon_cartesian
-from uxarray.grid.point_in_face import _face_contains_point_from_edges
+from uxarray.grid.point_in_face import _point_in_face_from_grid
 from uxarray.grid.utils import _get_cartesian_face_edge_nodes_array
 
 
@@ -114,25 +114,14 @@ def test_cache_and_override_geodataframe(gridpath):
 
 
 def test_face_at_antimeridian():
-    """Test the function `point_in_face`, where the face crosses the 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_array(
-        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 _face_contains_point_from_edges(faces_edges_cartesian[0], point)
+    assert _point_in_face_from_grid(point, grid, fidx=0)
 
 
 def test_face_at_pole():
@@ -141,21 +130,9 @@ def test_face_at_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_array(
-        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 _face_contains_point_from_edges(faces_edges_cartesian[0], point)
+    assert _point_in_face_from_grid(point, grid, fidx=0)
 
 
 def test_face_normal_face():
@@ -166,19 +143,8 @@ def test_face_normal_face():
     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_array(
-        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 _face_contains_point_from_edges(faces_edges_cartesian[0], point)
+    assert _point_in_face_from_grid(point, grid, fidx=0)
 
 
 def test_haversine_distance_creation():
diff --git a/test/grid/geometry/test_point_in_face.py b/test/grid/geometry/test_point_in_face.py
index 8bf7b0f89..1829c7eeb 100644
--- a/test/grid/geometry/test_point_in_face.py
+++ b/test/grid/geometry/test_point_in_face.py
@@ -6,72 +6,45 @@
 from uxarray.constants import INT_FILL_VALUE
 from uxarray.grid.coordinates import _lonlat_rad_to_xyz
 from uxarray.grid.utils import _get_cartesian_face_edge_nodes_array, _get_cartesian_face_edge_nodes
-from uxarray.grid.point_in_face import _face_contains_point_from_edges
+from uxarray.grid.point_in_face import _point_in_face_from_grid
 
 
 def test_point_inside(gridpath):
     """Test the function `point_in_face`, where the points are all inside the face"""
 
-    # Open grid
     grid = ux.open_grid(gridpath("mpas", "QU", "mesh.QU.1920km.151026.nc"))
     grid.normalize_cartesian_coordinates()
 
     # Loop through each face
     for i in range(grid.n_face):
-        face_edges = _get_cartesian_face_edge_nodes(
-            i, grid.face_node_connectivity.values, grid.n_nodes_per_face.values, grid.node_x.values, grid.node_y.values, grid.node_z.values
-        )
-
         # 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 _face_contains_point_from_edges(face_edges, point_xyz)
+        assert _point_in_face_from_grid(point_xyz, grid, i)
 
 
 def test_point_outside(gridpath):
     """Test the function `point_in_face`, where the point is outside the face"""
 
-    # Open grid
     grid = ux.open_grid(gridpath("mpas", "QU", "mesh.QU.1920km.151026.nc"))
 
-    # Get the face edges of all faces in the grid
-    faces_edges_cartesian = _get_cartesian_face_edge_nodes_array(
-        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 _face_contains_point_from_edges(faces_edges_cartesian[0], point_xyz)
+    assert not _point_in_face_from_grid(point_xyz, grid, fidx=0)
 
 
 def test_point_on_node(gridpath):
     """Test the function `point_in_face`, when the point is on the node of the polygon"""
 
-    # Open grid
     grid = ux.open_grid(gridpath("mpas", "QU", "mesh.QU.1920km.151026.nc"))
 
-    # Get the face edges of all faces in the grid
-    faces_edges_cartesian = _get_cartesian_face_edge_nodes_array(
-        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]])
+    point_xyz = np.array([grid.node_x[0].values, grid.node_y[0].values, grid.node_z[0].values])
 
-    # Assert that the point is in the face when inclusive is true
-    assert _face_contains_point_from_edges(faces_edges_cartesian[0], point_xyz)
+    # Assert that the point is in the face
+    assert _point_in_face_from_grid(point_xyz, grid, fidx=0)
 
 
 def test_point_inside_close():
@@ -86,17 +59,8 @@ def test_point_inside_close():
 
     # Create the grid and face edges
     grid = ux.Grid.from_face_vertices(vertices_lonlat, latlon=True)
-    faces_edges_cartesian = _get_cartesian_face_edge_nodes_array(
-        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 _face_contains_point_from_edges(faces_edges_cartesian[0], point)
+    assert _point_in_face_from_grid(point, grid, fidx=0)
 
 
 def test_point_outside_close():
@@ -111,14 +75,5 @@ def test_point_outside_close():
 
     # Create the grid and face edges
     grid = ux.Grid.from_face_vertices(vertices_lonlat, latlon=True)
-    faces_edges_cartesian = _get_cartesian_face_edge_nodes_array(
-        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 _face_contains_point_from_edges(faces_edges_cartesian[0], point)
+
+    assert not _point_in_face_from_grid(point, grid, fidx=0)
diff --git a/uxarray/grid/geometry.py b/uxarray/grid/geometry.py
index 0d850b846..a8aa54013 100644
--- a/uxarray/grid/geometry.py
+++ b/uxarray/grid/geometry.py
@@ -13,7 +13,7 @@
 from uxarray.grid.intersections import (
     gca_gca_intersection,
 )
-from uxarray.grid.point_in_face import _face_contains_point_from_edges
+from uxarray.grid.point_in_face import _point_in_face
 from uxarray.grid.utils import _get_cartesian_face_edge_nodes
 from uxarray.utils.imports import _raise_hint_if_optional_deps_missing
 
@@ -1260,21 +1260,13 @@ def barycentric_coordinates_cartesian(polygon_xyz, point_xyz):
                 polygon_xyz[i + 2],
             )
 
-            # Get the triangle in terms of its edges for the `point_in_face` check
-            face_edge = _get_cartesian_face_edge_nodes(
-                face_idx=0,
-                face_node_connectivity=np.array([[0, 1, 2]]),
-                n_edges_per_face=np.array([3]),
-                node_x=np.array([node_0[0], node_1[0], node_2[0]], dtype=np.float64),
-                node_y=np.array([node_0[1], node_1[1], node_2[1]], dtype=np.float64),
-                node_z=np.array([node_0[2], node_1[2], node_2[2]], dtype=np.float64),
-            )
+            node_x = np.array([node_0[0], node_1[0], node_2[0]], dtype=np.float64)
+            node_y = np.array([node_0[1], node_1[1], node_2[1]], dtype=np.float64)
+            node_z = np.array([node_0[2], node_1[2], node_2[2]], dtype=np.float64)
+            nodes_idx = np.array([0, 1, 2])
 
             # Check to see if the point lies within the current triangle
-            contains_point = _face_contains_point_from_edges(
-                face_edge,
-                point_xyz,
-            )
+            contains_point = _point_in_face(point_xyz, nodes_idx, node_x, node_y, node_z)
 
             # If the point is in the current triangle, get the weights for that triangle
             if contains_point:
diff --git a/uxarray/grid/point_in_face.py b/uxarray/grid/point_in_face.py
index d84b109d5..94233461b 100644
--- a/uxarray/grid/point_in_face.py
+++ b/uxarray/grid/point_in_face.py
@@ -83,6 +83,19 @@ def _face_contains_point_from_edges(face_edges: np.ndarray, point: np.ndarray) -
     return np.abs(total) > np.pi
 
 
+def _point_in_face_from_grid(point: np.ndarray, grid: Grid, fidx: int):
+    """Returns whether this point lies within the indicated face of this grid.
+    Helper function providing convenient entry point into `_point_in_face`;
+    see `_point_in_face` for full docstring.
+    """
+    n_nodes_in_face = grid.n_nodes_per_face[fidx].item()
+    nodes_idx = grid.face_node_connectivity[fidx][:n_nodes_in_face].values
+    nodes_x = grid.node_x.values
+    nodes_y = grid.node_y.values
+    nodes_z = grid.node_z.values
+    return _point_in_face(point, nodes_idx, nodes_x, nodes_y, nodes_z)
+
+
 @njit(cache=True)
 def _point_in_face(
     point: np.ndarray,

From 738bd181315cfcd118d7c7fe844c423d73181305 Mon Sep 17 00:00:00 2001
From: Sam Evans <47793072+Sevans711@users.noreply.github.com>
Date: Mon, 5 Oct 2026 16:16:14 -0400
Subject: [PATCH 04/10] delete _face_contains_point_from_edges

(there are no more references to it; forgot to remove it during previous commit.)
---
 uxarray/grid/point_in_face.py | 61 -----------------------------------
 1 file changed, 61 deletions(-)

diff --git a/uxarray/grid/point_in_face.py b/uxarray/grid/point_in_face.py
index 94233461b..97df0d32e 100644
--- a/uxarray/grid/point_in_face.py
+++ b/uxarray/grid/point_in_face.py
@@ -22,67 +22,6 @@
     from uxarray.grid.grid import Grid
 
 
-@njit(cache=True)
-def _face_contains_point_from_edges(face_edges: np.ndarray, point: np.ndarray) -> bool:
-    """
-    Determine whether a point lies within a face using the spherical winding-number method.
-
-    This function sums the signed central angles between successive vertices of the face
-    as seen from `point`.  If the total absolute winding exceeds π, the point is inside.
-    Points exactly on a node or edge also count as inside.
-
-    Parameters
-    ----------
-    face_edges : np.ndarray, shape (n_edges, 2, 3)
-        Cartesian coordinates (unit-vectors) of each great-circle edge of the face.
-        Each row is [start_xyz, end_xyz].
-    point : np.ndarray, shape (3,)
-        3D unit-vector of the query point on the unit sphere.
-
-    Returns
-    -------
-    inside : bool
-        True if the point is inside the face or lies exactly on a node/edge; False otherwise.
-    """
-    # Check for an exact hit with any of the corner nodes
-    for e in range(face_edges.shape[0]):
-        if np.allclose(
-            face_edges[e, 0], point, rtol=ERROR_TOLERANCE, atol=ERROR_TOLERANCE
-        ):
-            return True
-        if np.allclose(
-            face_edges[e, 1], point, rtol=ERROR_TOLERANCE, atol=ERROR_TOLERANCE
-        ):
-            return True
-        if point_within_gca(point, face_edges[e, 0], face_edges[e, 1]):
-            return True
-
-    n = face_edges.shape[0]
-
-    total = 0.0
-    p = point
-    for i in range(n):
-        a = face_edges[i, 0]
-        b = face_edges[i + 1, 0] if i + 1 < n else face_edges[0, 0]
-
-        vi = a - p
-        vj = b - p
-
-        # check if you’re right on a vertex
-        if np.linalg.norm(vi) < ERROR_TOLERANCE or np.linalg.norm(vj) < ERROR_TOLERANCE:
-            return True
-
-        ang = _small_angle_of_2_vectors(vi, vj)
-
-        # determine sign from cross
-        c = np.cross(vi, vj)
-        sign = 1.0 if (c[0] * p[0] + c[1] * p[1] + c[2] * p[2]) >= 0.0 else -1.0
-
-        total += sign * ang
-
-    return np.abs(total) > np.pi
-
-
 def _point_in_face_from_grid(point: np.ndarray, grid: Grid, fidx: int):
     """Returns whether this point lies within the indicated face of this grid.
     Helper function providing convenient entry point into `_point_in_face`;

From a645a29ee68bae5152660df75e0aa1df6a3ba608 Mon Sep 17 00:00:00 2001
From: Sam Evans <47793072+Sevans711@users.noreply.github.com>
Date: Mon, 5 Oct 2026 16:28:51 -0400
Subject: [PATCH 05/10] =?UTF-8?q?optimize:=20point=5Fwithin=5Fgca,=20=5Fan?=
 =?UTF-8?q?gle=5Fof=5F2=5Fvectors,=20=E2=80=A6?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

and _point_in_face, by avoiding tiny numpy arrays.

Still only a small speedup for to_raster though... ~4.30s compared to ~4.45s from 594739b7 (and ~4.70s on main). Confirmed that the pytest suite passes locally, and that a produced raster is exactly the same as on main, via .equals().
---
 uxarray/grid/arcs.py          | 25 +++++++++++++++----------
 uxarray/grid/point_in_face.py |  5 ++---
 uxarray/grid/utils.py         | 13 ++++++++-----
 3 files changed, 25 insertions(+), 18 deletions(-)

diff --git a/uxarray/grid/arcs.py b/uxarray/grid/arcs.py
index 9ec3df40d..a20d11e55 100644
--- a/uxarray/grid/arcs.py
+++ b/uxarray/grid/arcs.py
@@ -9,6 +9,11 @@
 )
 from uxarray.grid.utils import _angle_of_2_vectors
 from uxarray.utils.computing import _cdp8, accucross
+from uxarray.utils.numba_math import (
+    _numba_cross3,
+    _numba_dot3,
+    _numba_sub3,
+)
 
 # Magnitude below which orient3d_on_sphere classifies a result as zero. For
 # double-precision unit-vector inputs this covers rounding error in the
@@ -40,11 +45,11 @@ def point_within_gca(pt_xyz, gca_a_xyz, gca_b_xyz):
 
     Parameters
     ----------
-    pt_xyz : numpy.ndarray
+    pt_xyz : iterable of length 3
         Cartesian coordinates of the point.
-    gca_a_xyz : numpy.ndarray
+    gca_a_xyz : iterable of length 3
         Cartesian coordinates of the first endpoint of the Great Circle Arc.
-    gca_b_xyz : numpy.ndarray
+    gca_b_xyz : iterable of length 3
         Cartesian coordinates of the second endpoint of the Great Circle Arc.
 
     Returns
@@ -66,25 +71,25 @@ def point_within_gca(pt_xyz, gca_a_xyz, gca_b_xyz):
     """
     # 1. Check if the input GCA spans exactly 180 degrees
     angle_ab = _angle_of_2_vectors(gca_a_xyz, gca_b_xyz)
-    if np.allclose(angle_ab, np.pi, rtol=0.0, atol=MACHINE_EPSILON):
+    if np.isclose(angle_ab, np.pi, rtol=0.0, atol=MACHINE_EPSILON):
         raise ValueError(
             "The input Great Circle Arc spans exactly 180 degrees, which can correspond to multiple planes. "
             "Consider breaking the Great Circle Arc into two smaller arcs."
         )  # (numba complains about f-strings, so don't put actual values in message.)
 
     # 2. Verify if the point lies on the plane of the GCA
-    cross_product = np.cross(gca_a_xyz, gca_b_xyz)
-    if not np.allclose(
-        np.dot(cross_product, pt_xyz), 0, rtol=MACHINE_EPSILON, atol=MACHINE_EPSILON
+    cross_product = _numba_cross3(gca_a_xyz, gca_b_xyz)
+    if not np.isclose(
+        _numba_dot3(cross_product, pt_xyz), 0, rtol=MACHINE_EPSILON, atol=MACHINE_EPSILON
     ):
         return False
 
     # 3. Check if the point lies within the Great Circle Arc interval
-    pt_a = gca_a_xyz - pt_xyz
-    pt_b = gca_b_xyz - pt_xyz
+    pt_a = _numba_sub3(gca_a_xyz, pt_xyz)
+    pt_b = _numba_sub3(gca_b_xyz, pt_xyz)
 
     # Use the dot product to determine the sign of the angle between pt_a and pt_b
-    cos_theta = np.dot(pt_a, pt_b)
+    cos_theta = _numba_dot3(pt_a, pt_b)
 
     # Return True if the point lies within the interval (smaller arc)
     if cos_theta < 0:
diff --git a/uxarray/grid/point_in_face.py b/uxarray/grid/point_in_face.py
index 97df0d32e..1799aee94 100644
--- a/uxarray/grid/point_in_face.py
+++ b/uxarray/grid/point_in_face.py
@@ -90,9 +90,8 @@ def _point_in_face(
         ai = nodes_idx[i]
         bi = nodes_idx[i + 1] if i < max_i_node else nodes_idx[0]
 
-        # TODO: avoid tiny numpy arrays, after rewriting point_within_gca to accept tuples
-        a = np.array([node_x[ai], node_y[ai], node_z[ai]])
-        b = np.array([node_x[bi], node_y[bi], node_z[bi]])
+        a = (node_x[ai], node_y[ai], node_z[ai])
+        b = (node_x[bi], node_y[bi], node_z[bi])
         if point_within_gca(point, a, b):
             return True
 
diff --git a/uxarray/grid/utils.py b/uxarray/grid/utils.py
index 4d01e5049..bfa389371 100644
--- a/uxarray/grid/utils.py
+++ b/uxarray/grid/utils.py
@@ -5,6 +5,8 @@
 from uxarray.constants import EDGE_NODE_SORT_THRESHOLD, INT_DTYPE, INT_FILL_VALUE
 from uxarray.utils.numba_math import (
     _numba_add3,
+    _numba_cross3,
+    _numba_dot3,
     _numba_mul3_scalar,
     _numba_norm3,
     _numba_sub3,
@@ -46,9 +48,9 @@ def _angle_of_2_vectors(u, v):
 
     Parameters
     ----------
-    u : numpy.ndarray
+    u : iterable of length 3
         The first 3D vector (float), originating from the center of the unit sphere.
-    v : numpy.ndarray
+    v : iterable of length 3
         The second 3D vector (float), originating from the center of the unit sphere.
 
     Returns
@@ -62,13 +64,14 @@ def _angle_of_2_vectors(u, v):
     - Special cases such as vectors aligned along the same longitude are handled explicitly.
     """
     # Compute the cross product to determine the direction of the normal
-    normal = np.cross(u, v)
+    normal = _numba_cross3(u, v)
 
     # Calculate the angle using arctangent of cross and dot products
-    angle_u_v_rad = np.arctan2(np.linalg.norm(normal), np.dot(u, v))
+    angle_u_v_rad = np.arctan2(_numba_norm3(normal), _numba_dot3(u, v))
 
     # Determine the direction of the angle
-    normal_z = np.dot(normal, np.array([0.0, 0.0, 1.0]))
+    normal_z = normal[2]
+
     if normal_z > 0:
         # Counterclockwise direction
         return angle_u_v_rad

From c878947c5f716f367eda37ff1f6fe0e0a94f8257 Mon Sep 17 00:00:00 2001
From: Sam Evans <47793072+Sevans711@users.noreply.github.com>
Date: Mon, 5 Oct 2026 16:33:26 -0400
Subject: [PATCH 06/10] _point_in_face support tuple point

Although, this actually provides no noticeable speed improvement to to_raster()...
---
 uxarray/grid/point_in_face.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/uxarray/grid/point_in_face.py b/uxarray/grid/point_in_face.py
index 1799aee94..6400876bd 100644
--- a/uxarray/grid/point_in_face.py
+++ b/uxarray/grid/point_in_face.py
@@ -37,7 +37,7 @@ def _point_in_face_from_grid(point: np.ndarray, grid: Grid, fidx: int):
 
 @njit(cache=True)
 def _point_in_face(
-    point: np.ndarray,
+    point: np.ndarray | tuple[float, float, float],
     nodes_idx: np.ndarray,
     node_x: np.ndarray,
     node_y: np.ndarray,
@@ -52,7 +52,7 @@ def _point_in_face(
 
     Parameters
     ----------
-    point : np.ndarray, shape (3,)
+    point : iterable of length 3
         3D unit-vector of the query point on the unit sphere.
     nodes_idx : np.ndarray, shape (n_nodes,)
         Node indices (within node_x, node_y, node_z) for precisely all nodes in this face.
@@ -127,7 +127,7 @@ def _point_in_face(
 def _set_faces_containing_point(
     result: np.ndarray,
     i: int,
-    point: np.ndarray,
+    point: np.ndarray | tuple[float, float, float],
     candidate_indices: np.ndarray,
     face_node_connectivity: np.ndarray,
     n_nodes_per_face: np.ndarray,
@@ -149,7 +149,7 @@ def _set_faces_containing_point(
         the maximum number of hits is n_max_face_nodes, because
         the "worst case" of point being a node would lead to hits
         of all faces it is a part of, but nothing else.
-    point : np.ndarray, shape (3,)
+    point : iterable of length 3
         Cartesian unit-vector of the query point.
     candidate_indices : np.ndarray, shape (k,)
         Array of face indices to test (e.g., from a k-d tree cull).
@@ -217,7 +217,7 @@ def _batch_point_in_face(
     for i in prange(n_points):
         start = offsets[i]
         end = offsets[i + 1]
-        p = points[i]
+        p = (points[i][0], points[i][1], points[i][2])
         cands = flat_candidate_indices[start:end]
 
         n_hits = _set_faces_containing_point(

From 1d4537df723dc310562fe23adf0eb96e70b5c961 Mon Sep 17 00:00:00 2001
From: Sam Evans <47793072+Sevans711@users.noreply.github.com>
Date: Mon, 5 Oct 2026 20:39:07 -0400
Subject: [PATCH 07/10] benchmarks for to_raster

---
 benchmarks/to_raster.py | 162 ++++++++++++++++++++++++++++++++++++++++
 1 file changed, 162 insertions(+)
 create mode 100644 benchmarks/to_raster.py

diff --git a/benchmarks/to_raster.py b/benchmarks/to_raster.py
new file mode 100644
index 000000000..a23c728fd
--- /dev/null
+++ b/benchmarks/to_raster.py
@@ -0,0 +1,162 @@
+import cartopy.crs as ccrs
+import matplotlib
+matplotlib.use("Agg")  # default backend causes ASV to crash on MacOS.
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy.spatial import ConvexHull
+import uxarray as ux
+
+
+
+from .helpers._memsize import dataset_nbytes, grid_nbytes
+from .helpers._peakmem import numba_threads, peak_allocated
+
+
+def make_sphere_grid(n_faces, area_ratio, *, seed=0):
+    """
+    Returns a variable-resolution triangular grid on the unit sphere.
+
+    n_faces    : target number of triangular faces (result is approximate)
+    area_ratio : (face area at poles) / (face area at equator)
+    seed       : RNG seed for the per-ring longitude offsets
+
+    Returns
+    -------
+    node_lon, node_lat : (n_nodes,) degrees (lon in [-180, 180))
+    face_node_connectivity : (n_faces, 3) int array, counter-clockwise
+                             when viewed from outside the sphere.
+    """
+    # note: Claude made this function, aside from some small edits by hand.
+    # The goal was to generate ANY "reasonable-looking" grid with varying face areas,
+    # because to_raster() performance was bad for grids with significant area variations.
+
+    # You can check it produces reasonable-looking results by doing something like:
+    # areas = make_sphere_grid(1e4, 100).compute_face_areas(as_uxarray=True)
+    # areas.plot()    # to view the actual areas
+    # (areas*0 + np.arange(len(areas.n_face))).plot()   # to view face indices.
+
+    edge_ratio = np.sqrt(area_ratio)          # edge length ~ sqrt(area)
+
+    def h_of(lat, h_pole):                    # target edge length vs latitude
+        return (h_pole / edge_ratio) * edge_ratio ** (np.abs(lat) / (np.pi / 2))
+
+    def layout(h_pole):
+        """Ring latitudes (all, south to north) and nodes per ring."""
+        lats = [0.0]
+        while True:
+            nxt = lats[-1] + 0.866 * h_of(lats[-1], h_pole)
+            if nxt > np.pi / 2 - 0.6 * h_pole:   # leave room for the pole node
+                break
+            lats.append(nxt)
+        lats = np.array(lats)
+        lats = np.concatenate([-lats[:0:-1], lats])
+        counts = np.maximum(
+            3, np.round(2 * np.pi * np.cos(lats) / h_of(lats, h_pole))
+        ).astype(int)
+        return lats, counts
+
+    # --- choose h_pole so that the face count is close to n_faces
+    # analytic estimate: faces ~ C / h_pole^2 (equilateral triangles)
+    phi = np.linspace(-np.pi / 2, np.pi / 2, 20001)
+    C = np.sum(2 * np.pi * np.cos(phi) / ((np.sqrt(3) / 4) * h_of(phi, 1.0) ** 2)) \
+        * (phi[1] - phi[0])
+    h0 = np.sqrt(C / n_faces)
+
+    def faces_for(h_pole):                    # triangles = 2*nodes - 4 on a sphere
+        return 2 * (layout(h_pole)[1].sum() + 2) - 4
+
+    lo, hi = np.log(h0 / 2), np.log(h0 * 2)   # refine with the true ring counts
+    for _ in range(30):
+        mid = 0.5 * (lo + hi)
+        if faces_for(np.exp(mid)) > n_faces:
+            lo = mid                          # too many faces -> larger h
+        else:
+            hi = mid
+    h_pole = np.exp(0.5 * (lo + hi))
+
+    # --- nodes on each ring
+    rng = np.random.default_rng(seed)
+    lats, counts = layout(h_pole)
+    lon_list, lat_list = [], []
+    for lat, n in zip(lats, counts):
+        lon_list.append(rng.uniform(0, 2 * np.pi) + 2 * np.pi * np.arange(n) / n)
+        lat_list.append(np.full(n, lat))
+    lon_list += [np.array([0.0]), np.array([0.0])]            # poles
+    lat_list += [np.array([np.pi / 2]), np.array([-np.pi / 2])]
+
+    lon = np.concatenate(lon_list)
+    lat = np.concatenate(lat_list)
+
+    # --- spherical Delaunay via 3D convex hull
+    xyz = np.column_stack([np.cos(lat) * np.cos(lon),
+                           np.cos(lat) * np.sin(lon),
+                           np.sin(lat)])
+    faces = ConvexHull(xyz).simplices.copy()
+
+    # orient counter-clockwise as seen from outside
+    a, b, c = xyz[faces[:, 0]], xyz[faces[:, 1]], xyz[faces[:, 2]]
+    flip = np.einsum("ij,ij->i", np.cross(b - a, c - a), a + b + c) < 0
+    faces[flip] = faces[flip][:, [0, 2, 1]]
+
+    node_lon = (np.degrees(lon) + 180.0) % 360.0 - 180.0
+    node_lat = np.degrees(lat)
+    connectivity = faces.astype(np.int64)
+    return ux.Grid.from_topology(node_lon, node_lat, connectivity)
+
+
+
+class ToRaster:
+    """Benchmark the time it takes to call to_raster(),
+    on a variable-resolution grid.
+    """
+    param_names = ["n_faces_and_area_ratio",]
+    params = [[
+        (1e3, 10.0),
+        (1e3, 100.0),
+        # (1e3, 1000.0) isn't properly resolved;
+        # it would just make a bunch of faces stretching from pole to equator.
+        (1e4, 10.0),
+        #(1e4, 100.0),
+        #(1e4, 1000.0),
+        #(1e5, 10.0),
+        #(1e5, 100.0),
+        #(1e5, 1000.0),
+    ]]
+
+    # if any combination takes longer than 3 mins, give up.
+    timeout = 180
+
+    def _warmup(self):
+        # warm up numba, using a small example
+        grid = make_sphere_grid(1e3, 10.0)
+        data = grid.compute_face_areas(as_uxarray=True)
+        # data could be any values; let's use areas for fun :)
+        fig, ax = plt.subplots(subplot_kw={"projection": ccrs.Robinson()})
+        ax.set_global()
+        data.to_raster(ax=ax)
+        plt.close()
+
+    def setup(self, n_faces_and_area_ratio):
+        n_faces, area_ratio = n_faces_and_area_ratio
+        grid = make_sphere_grid(n_faces, area_ratio)
+        self.data = grid.compute_face_areas(as_uxarray=True)
+        self._warmup()
+        self.fig, self.ax = plt.subplots(subplot_kw={"projection": ccrs.Robinson()})
+        self.ax.set_global()
+
+    def teardown(self, n_faces_and_area_ratio):
+        del self.data
+        del self.fig
+        del self.ax
+        plt.close()
+
+    def time_to_raster(self, n_faces_and_area_ratio):
+        """Time to call to_raster() on a variable-resolution grid."""
+        self.data.to_raster(ax=self.ax)
+
+    def track_peakmem(self, n_faces_and_area_ratio):
+        """High-water allocation of to_raster()"""
+        with numba_threads(1):
+            return peak_allocated(lambda: self.data.to_raster(ax=self.ax))
+
+    track_peakmem.unit = "bytes"

From b99bb4f9d3d6aef579f67ac993cfda28cc65da36 Mon Sep 17 00:00:00 2001
From: Sam Evans <47793072+Sevans711@users.noreply.github.com>
Date: Mon, 5 Oct 2026 20:40:16 -0400
Subject: [PATCH 08/10] ruff formatting

---
 uxarray/grid/arcs.py          |  5 ++++-
 uxarray/grid/geometry.py      |  5 +++--
 uxarray/grid/point_in_face.py | 18 ++++++++++++++----
 3 files changed, 21 insertions(+), 7 deletions(-)

diff --git a/uxarray/grid/arcs.py b/uxarray/grid/arcs.py
index a20d11e55..dfe410988 100644
--- a/uxarray/grid/arcs.py
+++ b/uxarray/grid/arcs.py
@@ -80,7 +80,10 @@ def point_within_gca(pt_xyz, gca_a_xyz, gca_b_xyz):
     # 2. Verify if the point lies on the plane of the GCA
     cross_product = _numba_cross3(gca_a_xyz, gca_b_xyz)
     if not np.isclose(
-        _numba_dot3(cross_product, pt_xyz), 0, rtol=MACHINE_EPSILON, atol=MACHINE_EPSILON
+        _numba_dot3(cross_product, pt_xyz),
+        0,
+        rtol=MACHINE_EPSILON,
+        atol=MACHINE_EPSILON,
     ):
         return False
 
diff --git a/uxarray/grid/geometry.py b/uxarray/grid/geometry.py
index a8aa54013..26f8defc9 100644
--- a/uxarray/grid/geometry.py
+++ b/uxarray/grid/geometry.py
@@ -14,7 +14,6 @@
     gca_gca_intersection,
 )
 from uxarray.grid.point_in_face import _point_in_face
-from uxarray.grid.utils import _get_cartesian_face_edge_nodes
 from uxarray.utils.imports import _raise_hint_if_optional_deps_missing
 
 POLE_POINTS_XYZ = {
@@ -1266,7 +1265,9 @@ def barycentric_coordinates_cartesian(polygon_xyz, point_xyz):
             nodes_idx = np.array([0, 1, 2])
 
             # Check to see if the point lies within the current triangle
-            contains_point = _point_in_face(point_xyz, nodes_idx, node_x, node_y, node_z)
+            contains_point = _point_in_face(
+                point_xyz, nodes_idx, node_x, node_y, node_z
+            )
 
             # If the point is in the current triangle, get the weights for that triangle
             if contains_point:
diff --git a/uxarray/grid/point_in_face.py b/uxarray/grid/point_in_face.py
index 6400876bd..d2f361aec 100644
--- a/uxarray/grid/point_in_face.py
+++ b/uxarray/grid/point_in_face.py
@@ -7,7 +7,7 @@
 
 from uxarray.constants import ERROR_TOLERANCE, INT_DTYPE, INT_FILL_VALUE
 from uxarray.grid.arcs import point_within_gca
-from uxarray.grid.utils import _get_cartesian_face_edge_nodes, _small_angle_of_2_vectors
+from uxarray.grid.utils import _small_angle_of_2_vectors
 from uxarray.utils.numba_math import (
     _numba_allclose3,
     _numba_cross3,
@@ -80,7 +80,9 @@ def _point_in_face(
     for i in range(n_nodes):
         node_idx = nodes_idx[i]
         node_xyz = (node_x[node_idx], node_y[node_idx], node_z[node_idx])
-        if _numba_allclose3(node_xyz, point, rtol=ERROR_TOLERANCE, atol=ERROR_TOLERANCE):
+        if _numba_allclose3(
+            node_xyz, point, rtol=ERROR_TOLERANCE, atol=ERROR_TOLERANCE
+        ):
             return True
 
     # Check whether point lies on any edge of the face
@@ -168,7 +170,7 @@ def _set_faces_containing_point(
     count = 0
     for k in range(candidate_indices.shape[0]):
         fidx = candidate_indices[k]
-        nodes_idx = face_node_connectivity[fidx][:n_nodes_per_face[fidx]]
+        nodes_idx = face_node_connectivity[fidx][: n_nodes_per_face[fidx]]
         if _point_in_face(point, nodes_idx, node_x, node_y, node_z):
             result[i, count] = fidx
             count += 1
@@ -221,7 +223,15 @@ def _batch_point_in_face(
         cands = flat_candidate_indices[start:end]
 
         n_hits = _set_faces_containing_point(
-            results, i, p, cands, face_node_connectivity, n_nodes_per_face, node_x, node_y, node_z
+            results,
+            i,
+            p,
+            cands,
+            face_node_connectivity,
+            n_nodes_per_face,
+            node_x,
+            node_y,
+            node_z,
         )
         counts[i] = n_hits
 

From a38a1f4f3a5f0429386e840e612c598c11be3665 Mon Sep 17 00:00:00 2001
From: Sam Evans <47793072+Sevans711@users.noreply.github.com>
Date: Tue, 6 Oct 2026 09:33:57 -0400
Subject: [PATCH 09/10] uncomment the other to_raster benchmark cases

---
 benchmarks/to_raster.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/benchmarks/to_raster.py b/benchmarks/to_raster.py
index a23c728fd..ef2e31e60 100644
--- a/benchmarks/to_raster.py
+++ b/benchmarks/to_raster.py
@@ -116,11 +116,11 @@ class ToRaster:
         # (1e3, 1000.0) isn't properly resolved;
         # it would just make a bunch of faces stretching from pole to equator.
         (1e4, 10.0),
-        #(1e4, 100.0),
-        #(1e4, 1000.0),
-        #(1e5, 10.0),
-        #(1e5, 100.0),
-        #(1e5, 1000.0),
+        (1e4, 100.0),
+        (1e4, 1000.0),
+        (1e5, 10.0),
+        (1e5, 100.0),
+        (1e5, 1000.0),
     ]]
 
     # if any combination takes longer than 3 mins, give up.

From 30ff14eb2e56e442dfcda15d12239a8faf405387 Mon Sep 17 00:00:00 2001
From: Sam Evans <47793072+Sevans711@users.noreply.github.com>
Date: Tue, 6 Oct 2026 12:00:44 -0400
Subject: [PATCH 10/10] fewer to_raster benchmark param combinations

---
 benchmarks/to_raster.py | 15 +++++++++------
 1 file changed, 9 insertions(+), 6 deletions(-)

diff --git a/benchmarks/to_raster.py b/benchmarks/to_raster.py
index ef2e31e60..fb0cb45b6 100644
--- a/benchmarks/to_raster.py
+++ b/benchmarks/to_raster.py
@@ -111,16 +111,19 @@ class ToRaster:
     """
     param_names = ["n_faces_and_area_ratio",]
     params = [[
-        (1e3, 10.0),
-        (1e3, 100.0),
-        # (1e3, 1000.0) isn't properly resolved;
+        ## commented-out below are some additional reasonable parameter combinations,
+        ## excluded from default benchmark suite for brevity, but kept because they
+        ## could be uncommented for local runs or more detailed profiling.
+        # (1e3, 10.0),
+        # (1e3, 100.0),
+        ## never do (1e3, 1000.0), it isn't properly resolved;
         # it would just make a bunch of faces stretching from pole to equator.
         (1e4, 10.0),
         (1e4, 100.0),
-        (1e4, 1000.0),
-        (1e5, 10.0),
+        # (1e4, 1000.0),
+        # (1e5, 10.0),
         (1e5, 100.0),
-        (1e5, 1000.0),
+        # (1e5, 1000.0),
     ]]
 
     # if any combination takes longer than 3 mins, give up.