From f144abb27a393ffaf5429c123427bb434baacc78 Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Fri, 10 Jul 2026 16:10:47 -0400 Subject: [PATCH 1/3] Tighten bounding cube for regional spherical domains When the mesh type is spherical, we now use the actual x,y,z bounding box rather than the unit cube. For regional runs that leverage the spherical barycentric coordinates and built in unit conversion in the velocity field interpolators, this reduces the hit count per hash cell and can improve the search time in the spatialhash.query() Note that quantize_coordinates now clips in float space before casting to uint32. This matters now that spherical queries can fall outside the hash-grid bounds (previously impossible with the unit cube). A negative normalized value cast to uint32 wraps to a huge number, which would have mapped below-range queries to the top bin (with morton code 2^32) erroneously. For spherical grids, points outside the regional box now quantize to edge bins rather than interior bins. In this case, they either fail the exact Morton-key match or are rejected by the point-in-cell check; either case results in the same outcome as before (GRID_SEARCH_ERROR) --- src/parcels/_core/spatialhash.py | 46 +++++++++++++++++++------------- 1 file changed, 27 insertions(+), 19 deletions(-) diff --git a/src/parcels/_core/spatialhash.py b/src/parcels/_core/spatialhash.py index c4f4d0a64..5a6564cad 100644 --- a/src/parcels/_core/spatialhash.py +++ b/src/parcels/_core/spatialhash.py @@ -46,16 +46,18 @@ def __init__( if isinstance_noimport(grid, "XGrid"): self._coord_dim = 2 # Number of computational coordinates is 2 (bilinear interpolation) if self._source_grid._mesh == "spherical": - # Boundaries of the hash grid are the unit cube - self._xmin = -1.0 - self._ymin = -1.0 - self._zmin = -1.0 - self._xmax = 1.0 - self._ymax = 1.0 - self._zmax = 1.0 # Compute the cell centers of the source grid (for now, assuming Xgrid) lon = np.deg2rad(self._source_grid.lon) lat = np.deg2rad(self._source_grid.lat) x, y, z = _latlon_rad_to_xyz(lat, lon) + # Boundaries of the hash grid are the Cartesian bounding box of the + # transformed grid, so that regional domains retain full quantization + # resolution instead of spreading it over the whole unit cube + self._xmin = x.min() + self._xmax = x.max() + self._ymin = y.min() + self._ymax = y.max() + self._zmin = z.min() + self._zmax = z.max() _xbound = np.stack( ( x[:-1, :-1], @@ -149,19 +151,20 @@ def __init__( elif isinstance_noimport(grid, "UxGrid"): self._coord_dim = grid.uxgrid.n_max_face_nodes # Number of barycentric coordinates if self._source_grid._mesh == "spherical": - # Boundaries of the hash grid are the unit cube - self._xmin = -1.0 - self._ymin = -1.0 - self._zmin = -1.0 - self._xmax = 1.0 - self._ymax = 1.0 - self._zmax = 1.0 # Compute the cell centers of the source grid (for now, assuming Xgrid) # Reshape node coordinates to (nfaces, nnodes_per_face) nids = self._source_grid.uxgrid.face_node_connectivity.values lon = self._source_grid.uxgrid.node_lon.values[nids] lat = self._source_grid.uxgrid.node_lat.values[nids] - x, y, z = _latlon_rad_to_xyz(np.deg2rad(lat), np.deg2rad(lon)) _xbound, _ybound, _zbound = _latlon_rad_to_xyz(np.deg2rad(lat), np.deg2rad(lon)) + # Boundaries of the hash grid are the Cartesian bounding box of the + # transformed grid, so that regional domains retain full quantization + # resolution instead of spreading it over the whole unit cube + self._xmin = _xbound.min() + self._xmax = _xbound.max() + self._ymin = _ybound.min() + self._ymax = _ybound.max() + self._zmin = _zbound.min() + self._zmax = _zbound.max() # Compute bounding box of each face self._xlow = np.atleast_2d(np.min(_xbound, axis=-1)) @@ -588,10 +591,15 @@ def quantize_coordinates(x, y, z, xmin, xmax, ymin, ymax, zmin, zmax, bitwidth=1 zn = np.where(dz != 0, (z - zmin) / dz, 0.0) # --- 2) Quantize to (0..bitwidth). --- - # Multiply by bitwidth, round down, and clip to be safe against slight overshoot. - xq = np.clip((xn * bitwidth).astype(np.uint32), 0, bitwidth) - yq = np.clip((yn * bitwidth).astype(np.uint32), 0, bitwidth) - zq = np.clip((zn * bitwidth).astype(np.uint32), 0, bitwidth) + # Multiply by bitwidth, round down, and clip to be safe against overshoot. + # Clip in float space before casting: out-of-range queries (e.g., points outside + # a regional domain) would otherwise wrap around when a negative float is cast to uint32. + # NaN queries produce arbitrary codes here; they are discarded downstream by the + # finite-coordinate mask in SpatialHash.query. + with np.errstate(invalid="ignore"): + xq = np.clip(xn * bitwidth, 0, bitwidth).astype(np.uint32) + yq = np.clip(yn * bitwidth, 0, bitwidth).astype(np.uint32) + zq = np.clip(zn * bitwidth, 0, bitwidth).astype(np.uint32) return xq, yq, zq From 312901d2676b3073639dd8cf1da3fc73d765e447 Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Fri, 10 Jul 2026 16:17:49 -0400 Subject: [PATCH 2/3] Add structured grid spherical spatialhash query tests Tests for both points inside the domain and explicitly outside the domain to confirm correct behavior for out of bounds PIC checks --- tests/test_spatialhash.py | 34 ++++++++++++++++++++++++++++++++++ 1 file changed, 34 insertions(+) diff --git a/tests/test_spatialhash.py b/tests/test_spatialhash.py index 6233feb17..a0634fc9e 100644 --- a/tests/test_spatialhash.py +++ b/tests/test_spatialhash.py @@ -20,6 +20,40 @@ def test_invalid_positions(): assert np.all(i == -3) +def test_spherical_regional_bounds(): + """Hash-grid bounds for spherical meshes are the Cartesian bounding box of the + (regional) grid, not the whole unit cube, so quantization resolution is not + wasted on parts of the sphere the grid does not cover. + """ + ds = datasets["2d_left_rotated"] + grid = FieldSet.from_sgrid_conventions(ds, mesh="spherical").data_g.grid + spatialhash = grid.get_spatial_hash() + + extents = np.array( + [ + spatialhash._xmax - spatialhash._xmin, + spatialhash._ymax - spatialhash._ymin, + spatialhash._zmax - spatialhash._zmin, + ] + ) + assert np.all(extents > 0.0) + assert np.all(extents < 2.0) # strictly tighter than the unit cube + + # Queries at cell centers must still resolve to the correct cell + lon, lat = grid.lon, grid.lat + clon = 0.25 * (lon[:-1, :-1] + lon[:-1, 1:] + lon[1:, :-1] + lon[1:, 1:]) + clat = 0.25 * (lat[:-1, :-1] + lat[:-1, 1:] + lat[1:, :-1] + lat[1:, 1:]) + jj, ii = np.meshgrid(np.arange(clat.shape[0]), np.arange(clat.shape[1]), indexing="ij") + j, i, _ = spatialhash.query(clat.ravel(), clon.ravel()) + assert np.array_equal(j, jj.ravel()) + assert np.array_equal(i, ii.ravel()) + + # Points far outside the regional domain must not match any cell + j, i, _ = spatialhash.query([-60.0, 80.0], [120.0, -150.0]) + assert np.all(j == -3) + assert np.all(i == -3) + + def test_mixed_positions(): ds = datasets["2d_left_rotated"] grid = FieldSet.from_sgrid_conventions(ds, mesh="flat").data_g.grid From ac22b349e78f6bd92c92970f1bdf9c02cbc1db98 Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Fri, 10 Jul 2026 16:18:59 -0400 Subject: [PATCH 3/3] Update documentation on the unstructured grid search to note change to spherical bounding box --- docs/development/unstructured_grid_search.md | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/docs/development/unstructured_grid_search.md b/docs/development/unstructured_grid_search.md index 170221fac..bedc1aa5e 100644 --- a/docs/development/unstructured_grid_search.md +++ b/docs/development/unstructured_grid_search.md @@ -42,12 +42,12 @@ The algorithm has two phases: **initialisation** (once, when the grid is first u The hash grid is three-dimensional regardless of the source grid: -| Source grid | Mesh type | Hash-grid space | -| ----------- | --------- | ------------------------------------- | -| `XGrid` | spherical | Cartesian unit cube (lon/lat → x,y,z) | -| `XGrid` | flat | 2-D lon/lat bounding box (z set to 0) | -| `UxGrid` | spherical | Cartesian unit cube | -| `UxGrid` | flat | 2-D lon/lat bounding box | +| Source grid | Mesh type | Hash-grid space | +| ----------- | --------- | ---------------------------------------- | +| `XGrid` | spherical | Cartesian bounding box (lon/lat → x,y,z) | +| `XGrid` | flat | 2-D lon/lat bounding box (z set to 0) | +| `UxGrid` | spherical | Cartesian bounding box (lon/lat → x,y,z) | +| `UxGrid` | flat | 2-D lon/lat bounding box | For spherical grids, node coordinates are converted from degrees to radians and then to Cartesian (x, y, z) on the unit sphere using `parcels._core.index_search._latlon_rad_to_xyz`: @@ -57,7 +57,7 @@ y = sin(lon) * cos(lat) z = sin(lat) ``` -The hash grid then spans the unit cube `[-1, 1]³`. Working in Cartesian space avoids the longitude wrap-around discontinuity that would otherwise cause the bounding boxes of cells crossing the antimeridian to erroneously span the entire domain. +The hash grid then spans the axis-aligned Cartesian bounding box of the transformed grid nodes. For a global grid this is (nearly) the unit cube `[-1, 1]³`; for a regional grid it is much tighter, so the full quantization resolution (1024 bins per axis) is spent on the region actually covered by the grid rather than the whole sphere. Working in Cartesian space avoids the longitude wrap-around discontinuity that would otherwise cause the bounding boxes of cells crossing the antimeridian to erroneously span the entire domain. For flat meshes the hash grid simply spans `[lon_min, lon_max] × [lat_min, lat_max]`, with z fixed at 0.