Fix `tutorial_nested_grids.ipynb` - grid searching

[VeckoTheGecko]

OK - running into a separate (unrelated) problem surfaced by this change which I think runs deeper...

docs/user_guide/examples/tutorial_nestedgrids.ipynb combines different datasets together - but it was specifying "flat" for the gridIDs and "spherical" for the data itself (which although somehow worked - doesn't make sense to me, and doesn't align with the work we've been doing).

This is in lon/lat space, so I updated to use "spherical" and got grid searching errors.

---------------------------------------------------------------------------
GridSearchingError                        Traceback (most recent call last)
Cell In[10], [line 18](vscode-notebook-cell:?execution_count=10&line=18)
     14 def SampleGridID(particles, fieldset):
     15     particles.gridID = fieldset.GridID[particles]
     16 
     17 
---> [18](vscode-notebook-cell:?execution_count=10&line=18) pset.execute(
     19     SampleGridID,
     20     runtime=np.timedelta64(1, "s"),
     21     dt=np.timedelta64(1, "s"),

File ~/coding/repos/parcels/src/parcels/_core/particleset.py:453, in ParticleSet.execute(self, kernels, dt, endtime, runtime, output_file, verbose_progress)
    450 else:
    451     next_time = end_time
--> [453](https://file+.vscode-resource.vscode-cdn.net/Users/Hodgs004/coding/repos/parcels/docs/user_guide/examples/~/coding/repos/parcels/src/parcels/_core/particleset.py:453) self._kernel.execute(self, endtime=next_time, dt=dt)
    455 if next_output is not None:
    456     if np.abs(next_time - next_output) < 0.001:

File ~/coding/repos/parcels/src/parcels/_core/kernel.py:244, in Kernel.execute(self, pset, endtime, dt)
    242                 error_func(pset[inds].t)
    243             else:
--> [244](https://file+.vscode-resource.vscode-cdn.net/Users/Hodgs004/coding/repos/parcels/docs/user_guide/examples/~/coding/repos/parcels/src/parcels/_core/kernel.py:244)                 error_func(pset[inds].z, pset[inds].y, pset[inds].x)
    246 return pset

File ~/coding/repos/parcels/src/parcels/_core/statuscodes.py:85, in _raise_grid_searching_error(z, y, x)
     84 def _raise_grid_searching_error(z, y, x):
---> [85](https://file+.vscode-resource.vscode-cdn.net/Users/Hodgs004/coding/repos/parcels/docs/user_guide/examples/~/coding/repos/parcels/src/parcels/_core/statuscodes.py:85)     raise GridSearchingError(f"Grid searching failed at (z={z}, y={y}, x={x})")

GridSearchingError: Grid searching failed at (z=array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0.], dtype=float32), y=array([-18.        , -18.        , -18.        , -18.        ,
       -18.        , -15.052631  , -15.052631  , -15.052631  ,
       -15.052631  , -12.105263  , -12.105263  ,  -9.157895  ,
        -9.157895  ,  -9.157895  ,  -6.2105265 ,  -3.2631578 ,
        -3.2631578 ,  -3.2631578 ,  -3.2631578 ,  -3.2631578 ,
        -3.2631578 ,  -0.31578946,  -0.31578946,   2.631579  ,
         5.5789475 ,   8.526316  ,   8.526316  ,  11.473684  ,
        11.473684  ,  11.473684  ,  14.421053  ,  14.421053  ,
        14.421053  ,  17.368422  ,  17.368422  ,  17.368422  ,
        17.368422  ,  17.368422  ,  20.31579   ,  20.31579   ,
        23.263159  ,  26.210526  ,  29.157894  ,  29.157894  ,
        29.157894  ,  29.157894  ,  32.105263  ,  35.05263   ,
        35.05263   ,  38.        ,  38.        ,  38.        ],
      dtype=float32), x=array([ 0.25,  3.  ,  5.75,  8.5 , 11.25,  8.5 , 11.25, 14.  , 52.5 ,
       16.75, 49.75, -2.5 , 44.25, 47.  , 41.5 , -8.  ,  3.  , 16.75,
       19.5 , 30.5 , 38.75, -8.  , 33.25, -8.  , -8.  , -8.  , -5.25,
       -8.  , -5.25, 25.  , -8.  , -5.25, 25.  , -8.  , -5.25, -2.5 ,
       25.  , 44.25, 25.  , 47.  , 47.  , 49.75, 16.75, 22.25, 49.75,
       52.5 , 52.5 , 49.75, 55.25, 22.25, 55.25, 58.  ], dtype=float32))

Adapting the plot I was able to find the grid search failure locations

Using an LLM it was talking about there being a bug in the point-in-cell check, but to be honest I don't know how accurate that is.

LLM output

Spherical UxGrid search fails for large triangles

When using mesh="spherical" with UxGrid, grid searching fails for points inside large triangular faces. Two independent bugs:

Bug 1: SpatialHash Cartesian bounding box too tight

SpatialHash.__init__ (spatialhash.py:182-188) computes each face's bounding box from the Cartesian (x,y,z) coordinates of its vertices only. On a sphere, great-circle edges curve outward relative to the chord, so for large faces a query point on the sphere surface can lie inside the face but outside the vertex-only Cartesian bounding box. The hash never considers the face as a candidate.

Suggested fix: Include great-circle edge midpoints (average vertex Cartesian coords, renormalize to unit sphere) in the bounding box computation.

Bug 2: uxgrid_point_in_cell barycentric check too strict

uxgrid_point_in_cell (index_search.py:298-372) projects the query onto the face plane and computes barycentric coordinates using unsigned 3D triangle areas. For large spherical triangles, projection distortion causes sum(bcoords) to exceed 1 (e.g. 1.023 for a ~70°-wide face), failing the np.isclose(..., 1.0, rtol=1e-3) check at line 370.

Minimal reproducer

import numpy as np
import uxarray as ux
import xarray as xr
from shapely.geometry import Point, Polygon
from triangle import triangulate

from parcels._core.index_search import uxgrid_point_in_cell
from parcels._core.uxgrid import UxGrid

# Build a multi-face constrained Delaunay triangulation from three nested polygons
polygons = [
    np.array([(10, 15), (25, 10), (25, 25), (17, 36), (10, 32)]),
    np.array([(0, -5), (35, 0), (35, 25), (0, 20)]),
    np.array([(-10, -20), (60, -20), (60, 40), (-10, 40)]),
]
verts, segments = [], []
offset = 0
for poly in polygons:
    Ni = len(poly)
    verts.extend(poly.tolist())
    segments.extend([[offset + j, offset + ((j + 1) % Ni)] for j in range(Ni)])
    offset += Ni
B = triangulate(
    {"vertices": np.asarray(verts, dtype=float), "segments": np.asarray(segments, dtype=int)},
    "p",
)
points, face_tris = B["vertices"], B["triangles"].astype(int)

ds = xr.Dataset(
    {
        "node_lon": ("n_node", points[:, 0]),
        "node_lat": ("n_node", points[:, 1]),
        "face_node_connectivity": (("n_face", "n_max_face_nodes"), face_tris),
    },
    coords={
        "n_node": np.arange(points.shape[0]),
        "n_face": np.arange(face_tris.shape[0]),
    },
    attrs={"Conventions": "UGRID-1.0"},
)
uxgrid = ux.Grid(ds)
z = ux.UxDataArray(np.array([0.0, 1.0]), dims="z", uxgrid=uxgrid)
grid = UxGrid(uxgrid, z, mesh="spherical")

# --- Bug 1: SpatialHash misses points inside a large face ---
sh = grid.get_spatial_hash()
for lon, lat in [(0.25, -18.0), (3.0, -18.0), (5.75, -18.0)]:
    _, face_ids, _ = sh.query(
        np.array([lat], dtype=np.float32),
        np.array([lon], dtype=np.float32),
    )
    print(f"Bug 1 — ({lon}, {lat}): face_id={face_ids[0]} (expected >=0, got GRID_SEARCH_ERROR)")

# --- Bug 2: point-in-cell rejects a point near the edge of a large face ---
# Find which face contains (52.5, -15) by brute force — none pass the check
for fi in range(face_tris.shape[0]):
    is_in, coords = uxgrid_point_in_cell(
        grid,
        np.array([-15.0], dtype=np.float32),
        np.array([52.5], dtype=np.float32),
        np.array([fi]),
        np.array([fi]),
    )
    if abs(coords[0].sum() - 1.0) < 0.05:  # close but fails the rtol=1e-3 check
        print(
            f"Bug 2 — (52.5, -15) face {fi}: is_in={is_in[0]}, "
            f"bcoords sum={coords[0].sum():.4f} (fails isclose rtol=1e-3)"
        )

Output:

Bug 1 — (0.25, -18.0): face_id=-3 (expected >=0, got GRID_SEARCH_ERROR)
Bug 1 — (3.0, -18.0): face_id=-3 (expected >=0, got GRID_SEARCH_ERROR)
Bug 1 — (5.75, -18.0): face_id=-3 (expected >=0, got GRID_SEARCH_ERROR)
Bug 2 — (52.5, -15) face 12: is_in=0, bcoords sum=1.0234 (fails isclose rtol=1e-3)
Bug 2 — (52.5, -15) face 13: is_in=0, bcoords sum=1.0492 (fails isclose rtol=1e-3)

Bug 1 is the primary blocker — it prevents tutorial_nestedgrids.ipynb from running with mesh="spherical" on the triangulation FieldSet.

@erikvansebille did you run into something similar when authoring the notebook? cc @fluidnumericsJoe

Originally posted by @VeckoTheGecko in #2848 (comment)

---

cc @erikvansebille, @wyatt-fluidnumerics, @fluidnumericsJoe - nicer to continue here than in the original PR :)