From 26371a41b9b5819f0ac0e2bc573b668c9ea53b06 Mon Sep 17 00:00:00 2001 From: Wyatt Sieminski Date: Tue, 4 Aug 2026 22:08:06 +0000 Subject: [PATCH 01/11] added a filter to remove any cells in contact with NaN nodes from the SpatialHash table --- src/parcels/_core/spatialhash.py | 25 +++++++++++++++++++++++-- 1 file changed, 23 insertions(+), 2 deletions(-) diff --git a/src/parcels/_core/spatialhash.py b/src/parcels/_core/spatialhash.py index 208cbfe95..12cdfcf9d 100644 --- a/src/parcels/_core/spatialhash.py +++ b/src/parcels/_core/spatialhash.py @@ -258,7 +258,17 @@ def _total_hash_entries(self, bitwidth): nx = xqhigh.astype(np.int64) - xqlow + 1 ny = yqhigh.astype(np.int64) - yqlow + 1 nz = zqhigh.astype(np.int64) - zqlow + 1 - return int((nx * ny * nz).sum()) + # NaN values are not allowed in the SpatialHash table, so faces with a NaN + # bounding box do not contribute to the entry count + invalid_face = ( + np.isnan(self._xlow) + | np.isnan(self._xhigh) + | np.isnan(self._ylow) + | np.isnan(self._yhigh) + | np.isnan(self._zlow) + | np.isnan(self._zhigh) + ) + return int(np.where(invalid_face, 0, nx * ny * nz).sum()) def _initialize_hash_table(self): """Create a mapping that relates unstructured grid faces to hash indices by determining @@ -299,7 +309,18 @@ def _initialize_hash_table(self): nx = (xqhigh - xqlow + 1).astype(np.int32, copy=False) ny = (yqhigh - yqlow + 1).astype(np.int32, copy=False) nz = (zqhigh - zqlow + 1).astype(np.int32, copy=False) - num_hash_per_face = (nx * ny * nz).astype( + + # prevent NaN values from entering the SpatialHash table by setting their + # num_hash_per_face equal to 0 + invalid_face = ( + np.isnan(self._xlow) + | np.isnan(self._xhigh) + | np.isnan(self._ylow) + | np.isnan(self._yhigh) + | np.isnan(self._zlow) + | np.isnan(self._zhigh) + ).ravel() + num_hash_per_face = np.where(invalid_face, 0, nx * ny * nz).astype( np.int32, copy=False ) # Since nx, ny, nz are in the 10-bit range, their product fits in int32 # Sums over faces can exceed int32, so accumulate in int64 From 958bb6c4e7cff78fcfbdd58042c0dd9abdbfe6c6 Mon Sep 17 00:00:00 2001 From: Wyatt Sieminski Date: Tue, 4 Aug 2026 22:29:10 +0000 Subject: [PATCH 02/11] Added testing for filtering NaN nodes from the SpatialHash table --- tests/test_spatialhash.py | 48 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) diff --git a/tests/test_spatialhash.py b/tests/test_spatialhash.py index 76cdf88d4..153fa3c28 100644 --- a/tests/test_spatialhash.py +++ b/tests/test_spatialhash.py @@ -93,3 +93,51 @@ def test_mixed_positions(): assert i[0] == 14 # Actual value for 2d_left_rotated center assert j[1] == -3 assert i[1] == -3 + + +def test_nan_node_invalidates_touching_faces(): + """ + Any mesh face that touches a NaN node should not be added to the HashTable. + """ + ds = datasets["2d_left_rotated"] + grid = FieldSet.from_sgrid_conventions(ds, mesh="flat").data_g.grid + clat, clon, jj, ii = _cell_centers(grid) + + # Set one interior node to NaN, and calculate the indexes of faces that touch it. + nj, ni = 10, 10 + touching = [ + (nj - 1, ni - 1), + (nj - 1, ni), + (nj, ni - 1), + (nj, ni), + ] + grid._ds["lon"].values[nj, ni] = np.nan + grid._ds["lat"].values[nj, ni] = np.nan + spatialhash = grid.get_spatial_hash(reconstruct=True) + + # From the indexes of the faces that touch the NaN node, calculate their + # face_ids. + n_faces_x = clon.shape[1] + invalid_ids = set() + for j, i in touching: + invalid_ids.add(j * n_faces_x + i) + + # Get a set of all of the face_ids that are in the table, and assert that the + # ones touching the NaN node are not among them. + faces_in_table = set(np.unique(spatialhash._hash_table["faces"]).tolist()) + assert invalid_ids.isdisjoint(faces_in_table) + + # Queries landing on those 4 faces should return GridSearchErrors (-3). + touching_lat = np.array([clat[j, i] for j, i in touching]) + touching_lon = np.array([clon[j, i] for j, i in touching]) + j_touch, i_touch, _ = spatialhash.query(touching_lat, touching_lon) + assert np.all(j_touch == -3) + assert np.all(i_touch == -3) + + # All mesh cells not contacting the NaN node should resolve queries. + mask = np.ones(clat.shape, dtype=bool) + for j, i in touching: + mask[j, i] = False + j_rest, i_rest, _ = spatialhash.query(clat[mask], clon[mask]) + assert np.array_equal(j_rest, jj[mask]) + assert np.array_equal(i_rest, ii[mask]) From 99897923fab92b12ab90692532cae5fefe388a6a Mon Sep 17 00:00:00 2001 From: Wyatt Sieminski Date: Tue, 4 Aug 2026 23:14:01 +0000 Subject: [PATCH 03/11] Changed test_nan_node_invalidates_touching_faces to deep copy the lat/lon grid before injecting nan values. Otherwise the nan values break other tests. --- tests/test_spatialhash.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/tests/test_spatialhash.py b/tests/test_spatialhash.py index 153fa3c28..5d2e44195 100644 --- a/tests/test_spatialhash.py +++ b/tests/test_spatialhash.py @@ -103,6 +103,11 @@ def test_nan_node_invalidates_touching_faces(): grid = FieldSet.from_sgrid_conventions(ds, mesh="flat").data_g.grid clat, clon, jj, ii = _cell_centers(grid) + # `grid._ds` shares its lon/lat arrays with the module-level `datasets` fixture + # (from_sgrid_conventions does not copy them), so deep-copy before mutating, + # otherwise the NaN injected below leaks into every other test in the session. + grid._ds = grid._ds.copy(deep=True) + # Set one interior node to NaN, and calculate the indexes of faces that touch it. nj, ni = 10, 10 touching = [ From 8202e71b925319adb09a9929d3ef2a12d121da5e Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Wed, 5 Aug 2026 08:36:03 +0200 Subject: [PATCH 04/11] Fix pre-commit issue --- tests/test_spatialhash.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/tests/test_spatialhash.py b/tests/test_spatialhash.py index e2c672fed..df0532f06 100644 --- a/tests/test_spatialhash.py +++ b/tests/test_spatialhash.py @@ -122,9 +122,7 @@ def test_mixed_positions(): def test_nan_node_invalidates_touching_faces(): - """ - Any mesh face that touches a NaN node should not be added to the HashTable. - """ + """Any mesh face that touches a NaN node should not be added to the HashTable.""" ds = datasets["2d_left_rotated"] grid = FieldSet.from_sgrid_conventions(ds, mesh="flat").data_g.grid clat, clon, jj, ii = _cell_centers(grid) From 3f35124f30f7b82944f3852214614f2119928658 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Wed, 5 Aug 2026 08:37:50 +0200 Subject: [PATCH 05/11] Fixed typo in hastable.describe --- src/parcels/_reprs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/parcels/_reprs.py b/src/parcels/_reprs.py index 03323b2a7..99ba23d5e 100644 --- a/src/parcels/_reprs.py +++ b/src/parcels/_reprs.py @@ -298,7 +298,7 @@ def spatialhash_describe(spatialhash: SpatialHash) -> str: "Bitwidth (current / max)": f"{spatialhash._bitwidth} / 1023 (higher = finer resolution hash grid)", "Total hash cells": f"{n_total_cells:,d}", "Occupied hash cells": f"{n_occupied_cells:,d}, {n_occupied_cells / n_total_cells * 100:.4f}%", - "Total (hash cell --> gird face) entries": f"{n_entries:,d}", + "Total (hash cell --> grid face) entries": f"{n_entries:,d}", "Entries per occupied hash cell (avg)": f"{n_entries / n_occupied_cells:.2f}" if n_occupied_cells else "-", "Entries per face (avg)": f"{n_entries / n_faces:.2f}" if n_faces else "-", "Faces per occupied hash cell (min / mean / max)": ( From 3c5c41414d1e8faf06b7a27d05bfd5d88d3ba780 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Wed, 5 Aug 2026 08:48:18 +0200 Subject: [PATCH 06/11] Also fixing typo in unit test --- tests/test_spatialhash.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_spatialhash.py b/tests/test_spatialhash.py index df0532f06..d6cecde29 100644 --- a/tests/test_spatialhash.py +++ b/tests/test_spatialhash.py @@ -36,7 +36,7 @@ def test_spatialhash_describe(): Bitwidth (current / max) : 1023 / 1023 (higher = finer resolution hash grid) Total hash cells : 1,073,741,824 Occupied hash cells : 796,054, 0.0741% -Total (hash cell --> gird face) entries : 1,080,194 +Total (hash cell --> grid face) entries : 1,080,194 Entries per occupied hash cell (avg) : 1.36 Entries per face (avg) : 631.32 Faces per occupied hash cell (min / mean / max) : 1 / 1.36 / 4 From 62e95723ac645a190d4ea8d1e2a9dd7e56fd677c Mon Sep 17 00:00:00 2001 From: Wyatt Sieminski Date: Wed, 5 Aug 2026 13:55:23 +0000 Subject: [PATCH 07/11] Refactor NaN mask computation to be done in a helper function --- src/parcels/_core/spatialhash.py | 54 ++++++++++++++++++++++---------- 1 file changed, 37 insertions(+), 17 deletions(-) diff --git a/src/parcels/_core/spatialhash.py b/src/parcels/_core/spatialhash.py index 12cdfcf9d..fdc09d5c5 100644 --- a/src/parcels/_core/spatialhash.py +++ b/src/parcels/_core/spatialhash.py @@ -260,15 +260,8 @@ def _total_hash_entries(self, bitwidth): nz = zqhigh.astype(np.int64) - zqlow + 1 # NaN values are not allowed in the SpatialHash table, so faces with a NaN # bounding box do not contribute to the entry count - invalid_face = ( - np.isnan(self._xlow) - | np.isnan(self._xhigh) - | np.isnan(self._ylow) - | np.isnan(self._yhigh) - | np.isnan(self._zlow) - | np.isnan(self._zhigh) - ) - return int(np.where(invalid_face, 0, nx * ny * nz).sum()) + valid_face = _generate_valid_mask(self._xlow, self._xhigh, self._ylow, self._yhigh, self._zlow, self._zhigh) + return int(np.where(valid_face, nx * ny * nz, 0).sum()) def _initialize_hash_table(self): """Create a mapping that relates unstructured grid faces to hash indices by determining @@ -312,15 +305,10 @@ def _initialize_hash_table(self): # prevent NaN values from entering the SpatialHash table by setting their # num_hash_per_face equal to 0 - invalid_face = ( - np.isnan(self._xlow) - | np.isnan(self._xhigh) - | np.isnan(self._ylow) - | np.isnan(self._yhigh) - | np.isnan(self._zlow) - | np.isnan(self._zhigh) + valid_face = _generate_valid_mask( + self._xlow, self._xhigh, self._ylow, self._yhigh, self._zlow, self._zhigh ).ravel() - num_hash_per_face = np.where(invalid_face, 0, nx * ny * nz).astype( + num_hash_per_face = np.where(valid_face, nx * ny * nz, 0).astype( np.int32, copy=False ) # Since nx, ny, nz are in the 10-bit range, their product fits in int32 # Sums over faces can exceed int32, so accumulate in int64 @@ -758,3 +746,35 @@ def _encode_morton3d(x, y, z, xmin, xmax, ymin, ymax, zmin, zmax, bitwidth=1023) # Since our compact type fits in 30 bits, uint32 is enough. return code.astype(np.uint32) + + +def _generate_valid_mask(xlow, xhigh, ylow, yhigh, zlow, zhigh): + """ + Flag faces whose bounding box is fully defined, i.e. none of their 6 bounds + is NaN (a NaN indicates a corner node with a missing/masked coordinate). + + Parameters + ---------- + xlow, xhigh : array_like + Per-face bounding box in x. + ylow, yhigh : array_like + Per-face bounding box in y. + zlow, zhigh : array_like + Per-face bounding box in z. + + Returns + ------- + valid_face : ndarray of bool + Same shape as the inputs; True where the face's bounding box is finite, + False where it contains a NaN. + """ + invalid_face = ( + np.isnan(xlow) + | np.isnan(xhigh) + | np.isnan(ylow) + | np.isnan(yhigh) + | np.isnan(zlow) + | np.isnan(zhigh) + ) + + return ~invalid_face From 1405b02f3d455129c77087c764056f0e64629490 Mon Sep 17 00:00:00 2001 From: Wyatt Sieminski Date: Wed, 5 Aug 2026 14:10:03 +0000 Subject: [PATCH 08/11] Added assertion stating total number of faces in the spatialhash table is less than the total number of mesh faces when NaNs are in the mesh. --- tests/test_spatialhash.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/tests/test_spatialhash.py b/tests/test_spatialhash.py index d6cecde29..77149402f 100644 --- a/tests/test_spatialhash.py +++ b/tests/test_spatialhash.py @@ -156,6 +156,11 @@ def test_nan_node_invalidates_touching_faces(): faces_in_table = set(np.unique(spatialhash._hash_table["faces"]).tolist()) assert invalid_ids.isdisjoint(faces_in_table) + # The total number of mesh faces should be greater than the number in the table, + # since the NaN faces are filtered from the table. + n_total_faces = jj.size + assert n_total_faces > len(faces_in_table) + # Queries landing on those 4 faces should return GridSearchErrors (-3). touching_lat = np.array([clat[j, i] for j, i in touching]) touching_lon = np.array([clon[j, i] for j, i in touching]) From 1a1485b83f0ef8e6e4d845ad58c90c0a1c644a8b Mon Sep 17 00:00:00 2001 From: Wyatt Sieminski Date: Wed, 5 Aug 2026 14:17:17 +0000 Subject: [PATCH 09/11] Modified SpatialHash.describe() to also show total non-NaN faces. --- src/parcels/_reprs.py | 2 ++ tests/test_spatialhash.py | 1 + 2 files changed, 3 insertions(+) diff --git a/src/parcels/_reprs.py b/src/parcels/_reprs.py index 99ba23d5e..9377cbc41 100644 --- a/src/parcels/_reprs.py +++ b/src/parcels/_reprs.py @@ -287,6 +287,7 @@ def spatialhash_describe(spatialhash: SpatialHash) -> str: counts = hash_table["counts"] n_faces = int(np.size(spatialhash._xlow)) + n_valid_faces = int(np.unique(hash_table["faces"]).size) n_entries = int(hash_table["faces"].size) n_occupied_cells = int(hash_table["keys"].size) n_total_cells = (spatialhash._bitwidth + 1) ** 3 @@ -295,6 +296,7 @@ def spatialhash_describe(spatialhash: SpatialHash) -> str: "Grid type": type(grid).__name__, "Mesh": grid._mesh, "Total mesh faces": f"{n_faces:,d}", + "Valid (non-NaN) mesh faces": f"{n_valid_faces:,d}", "Bitwidth (current / max)": f"{spatialhash._bitwidth} / 1023 (higher = finer resolution hash grid)", "Total hash cells": f"{n_total_cells:,d}", "Occupied hash cells": f"{n_occupied_cells:,d}, {n_occupied_cells / n_total_cells * 100:.4f}%", diff --git a/tests/test_spatialhash.py b/tests/test_spatialhash.py index 77149402f..2f6222f97 100644 --- a/tests/test_spatialhash.py +++ b/tests/test_spatialhash.py @@ -33,6 +33,7 @@ def test_spatialhash_describe(): Grid type : XGrid Mesh : FlatMesh() Total mesh faces : 1,711 +Valid (non-NaN) mesh faces : 1,711 Bitwidth (current / max) : 1023 / 1023 (higher = finer resolution hash grid) Total hash cells : 1,073,741,824 Occupied hash cells : 796,054, 0.0741% From 2cad923c7b31a342af316400b88e806a897ab663 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 5 Aug 2026 14:17:36 +0000 Subject: [PATCH 10/11] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/parcels/_core/spatialhash.py | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/src/parcels/_core/spatialhash.py b/src/parcels/_core/spatialhash.py index 6b76ca03b..6624be588 100644 --- a/src/parcels/_core/spatialhash.py +++ b/src/parcels/_core/spatialhash.py @@ -788,12 +788,7 @@ def _generate_valid_mask(xlow, xhigh, ylow, yhigh, zlow, zhigh): False where it contains a NaN. """ invalid_face = ( - np.isnan(xlow) - | np.isnan(xhigh) - | np.isnan(ylow) - | np.isnan(yhigh) - | np.isnan(zlow) - | np.isnan(zhigh) + np.isnan(xlow) | np.isnan(xhigh) | np.isnan(ylow) | np.isnan(yhigh) | np.isnan(zlow) | np.isnan(zhigh) ) return ~invalid_face From ef8c7e167addf5843caa73a5d82a11147825c4fb Mon Sep 17 00:00:00 2001 From: Wyatt Sieminski Date: Wed, 5 Aug 2026 10:27:09 -0400 Subject: [PATCH 11/11] Move .ravel() to np.where line for readability Co-authored-by: Erik van Sebille --- src/parcels/_core/spatialhash.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/parcels/_core/spatialhash.py b/src/parcels/_core/spatialhash.py index 6624be588..15e4e47d0 100644 --- a/src/parcels/_core/spatialhash.py +++ b/src/parcels/_core/spatialhash.py @@ -308,10 +308,8 @@ def _initialize_hash_table(self): # prevent NaN values from entering the SpatialHash table by setting their # num_hash_per_face equal to 0 - valid_face = _generate_valid_mask( - self._xlow, self._xhigh, self._ylow, self._yhigh, self._zlow, self._zhigh - ).ravel() - num_hash_per_face = np.where(valid_face, nx * ny * nz, 0).astype( + valid_face = _generate_valid_mask(self._xlow, self._xhigh, self._ylow, self._yhigh, self._zlow, self._zhigh) + num_hash_per_face = np.where(valid_face.ravel(), nx * ny * nz, 0).astype( np.int32, copy=False ) # Since nx, ny, nz are in the 10-bit range, their product fits in int32 # Sums over faces can exceed int32, so accumulate in int64