Add spatial hash describe function

docs/development/unstructured_grid_search.md

@@ -5,6 +5,8 @@ This page documents the algorithm used in Parcels to locate which grid cell a pa
 - `src/parcels/_core/spatialhash.py` — `SpatialHash` class and Morton encoding utilities
 - `src/parcels/_core/index_search.py` — point-in-cell tests and the high-level search dispatch
 
+For debugging information regarding the spatial hash grid that underlies a curvilinear or unstructured grid, call `SpatialHash.describe()`, which prints a summary of the hash table's statistics. This table includes the total number of occupied hash cells, the percentage of hash cells that are occupied, the minimum, mean, and maximum number of mesh faces in a hash cell, and more.
+
 ---
 
 ## Motivation

src/parcels/_core/spatialhash.py

@@ -1,4 +1,6 @@
+import sys
 import warnings
+from typing import IO
 
 import numpy as np
 
@@ -10,6 +12,7 @@
 )
 from parcels._core.warnings import FieldSetWarning
 from parcels._python import isinstance_noimport
+from parcels._reprs import spatialhash_describe
 
 # Budget on the total number of (face, hash cell) pairs in the hash table:
 # max(_HASH_ENTRIES_PER_FACE * nfaces, _HASH_ENTRY_BUDGET_MIN).
@@ -524,6 +527,22 @@ def query(self, y, x):
             coords_best.reshape((num_queries, coordinates.shape[1])),
         )
 
+    def describe(self, buf: IO | None = None) -> None:
+        """
+        Summary of the SpatialHash's hash-table statistics (resolution, occupancy,
+        entry counts).
+
+        Parameters
+        ----------
+        buf : file-like, default: sys.stdout
+            writable buffer
+        """
+        if buf is None:
+            buf = sys.stdout
+        assert buf is not None
+
+        buf.write(spatialhash_describe(self))
+
 
 def _dilate_bits(n):
     """

src/parcels/_reprs.py

@@ -18,6 +18,7 @@
     from parcels import Field, FieldSet, ParticleSet
     from parcels._core.field import VectorField
     from parcels._core.model import ModelData
+    from parcels._core.spatialhash import SpatialHash
     from parcels._core.utils.time import TimeInterval
 
 
@@ -280,6 +281,36 @@ def fieldset_describe(fieldset: FieldSet) -> str:
     )
 
 
+def spatialhash_describe(spatialhash: SpatialHash) -> str:
+    grid = spatialhash._source_grid
+    hash_table = spatialhash._hash_table
+    counts = hash_table["counts"]
+
+    n_faces = int(np.size(spatialhash._xlow))
+    n_entries = int(hash_table["faces"].size)
+    n_occupied_cells = int(hash_table["keys"].size)
+    n_total_cells = (spatialhash._bitwidth + 1) ** 3
+
+    rows = {
+        "Grid type": type(grid).__name__,
+        "Mesh": grid._mesh,
+        "Total mesh faces": f"{n_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}%",
+        "Total (hash cell --> gird 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)": (
+            f"{counts.min():,d} / {counts.mean():.2f} / {counts.max():,d}" if n_occupied_cells else "-"
+        ),
+    }
+    key_width = max(len(k) for k in rows)
+    table = "\n".join(f"{k.ljust(key_width)} : {v}" for k, v in rows.items())
+
+    return "Spatial Hash Grid Statistics" + "\n" + table + "\n"
+
+
 def _get_parent_model(field: Field | VectorField) -> ModelData:
     if isinstance_noimport(field, "Field"):
         return field.model  # type:ignore[union-attr]

tests/test_spatialhash.py

@@ -1,3 +1,5 @@
+from io import StringIO
+
 import numpy as np
 
 from parcels._core.fieldset import FieldSet
@@ -20,6 +22,30 @@ def test_spatialhash_init():
     assert spatialhash is not None
 
 
+def test_spatialhash_describe():
+    ds = datasets["2d_left_rotated"]
+    grid = FieldSet.from_sgrid_conventions(ds, mesh="flat").data_g.grid
+    spatialhash = grid.get_spatial_hash()
+
+    io = StringIO()
+    expected = """\
+Spatial Hash Grid Statistics
+Grid type                                       : XGrid
+Mesh                                            : FlatMesh()
+Total 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%
+Total (hash cell --> gird 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
+"""
+    spatialhash.describe(io)
+    actual = io.getvalue()
+    assert actual == expected
+
+
 def test_invalid_positions():
     ds = datasets["2d_left_rotated"]
     grid = FieldSet.from_sgrid_conventions(ds, mesh="flat").data_g.grid