Skip to content

Feature/ux morton hashing #2173

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
fluidnumerics-joe wants to merge 3 commits into from
Closed

Feature/ux morton hashing #2173

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
a306bf2
Add support for using UXGrid in spatialhash.__init__()
fluidnumerics-joe Sep 2, 2025
4048685
Require _xc,_yc_,_zc private attributes are 2-d
fluidnumerics-joe Sep 2, 2025
95c8a31
Move get_spatial_hash() to basegrid and update uxgrid.search
fluidnumerics-joe Sep 2, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
26 changes: 26 additions & 0 deletions parcels/basegrid.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,8 @@

import numpy as np

from parcels.spatialhash import SpatialHash

if TYPE_CHECKING:
import numpy as np

Expand Down Expand Up @@ -140,6 +142,30 @@ def unravel_index(self, ei: int) -> dict[str, int]:
indices = _unravel(dims, ei)
return dict(zip(self.axes, indices, strict=True))

def get_spatial_hash(
self,
reconstruct=False,
):
"""Get the SpatialHash data structure of this Grid that allows for
fast face search queries. Face searches are used to find the faces that
a list of points, in spherical coordinates, are contained within.

Parameters
----------
reconstruct : bool, default=False
If true, reconstructs the spatial hash

Returns
-------
self._spatialhash : parcels.spatialhash.SpatialHash
SpatialHash instance

"""
if self._spatialhash is None or reconstruct:
self._spatialhash = SpatialHash(self, reconstruct)

return self._spatialhash

@property
@abstractmethod
def axes(self) -> list[str]:
Expand Down
206 changes: 87 additions & 119 deletions parcels/spatialhash.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,7 @@
import numpy as np

import parcels


class SpatialHash:
"""Custom data structure that is used for performing grid searches using Spatial Hashing. This class constructs an overlying
Expand All @@ -24,89 +26,100 @@ def __init__(
grid,
reconstruct=False,
):
# TODO : Enforce grid to be an instance of parcels.xgrid.XGrid
# Currently, this is not done due to circular import with parcels.xgrid

self._source_grid = grid
self.reconstruct = reconstruct

if self._source_grid._mesh == "spherical":
# Boundaries of the hash grid are the unit cube
if isinstance(grid, parcels.xgrid.XGrid):
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)
_xbound = np.stack(
(
x[:-1, :-1],
x[:-1, 1:],
x[1:, 1:],
x[1:, :-1],
),
axis=-1,
)
_ybound = np.stack(
(
y[:-1, :-1],
y[:-1, 1:],
y[1:, 1:],
y[1:, :-1],
),
axis=-1,
)
_zbound = np.stack(
(
z[:-1, :-1],
z[:-1, 1:],
z[1:, 1:],
z[1:, :-1],
),
axis=-1,
)
# Compute centroid locations of each cells
self._xc = np.mean(_xbound, axis=-1)
self._yc = np.mean(_ybound, axis=-1)
self._zc = np.mean(_zbound, axis=-1)

else:
# Boundaries of the hash grid are the bounding box of the source grid
self._xmin = self._source_grid.lon.min()
self._xmax = self._source_grid.lon.max()
self._ymin = self._source_grid.lat.min()
self._ymax = self._source_grid.lat.max()
# setting min and max below is needed for mesh="flat"
self._zmin = 0.0
self._zmax = 0.0
x = self._source_grid.lon
y = self._source_grid.lat

_xbound = np.stack(
(
x[:-1, :-1],
x[:-1, 1:],
x[1:, 1:],
x[1:, :-1],
),
axis=-1,
)
_ybound = np.stack(
(
y[:-1, :-1],
y[:-1, 1:],
y[1:, 1:],
y[1:, :-1],
),
axis=-1,
)
# Compute centroid locations of each cells
self._xc = np.mean(_xbound, axis=-1)
self._yc = np.mean(_ybound, axis=-1)
self._zc = np.zeros_like(self._xc)
else:
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)
_xbound = np.stack(
(
x[:-1, :-1],
x[:-1, 1:],
x[1:, 1:],
x[1:, :-1],
),
axis=-1,
)
_ybound = np.stack(
(
y[:-1, :-1],
y[:-1, 1:],
y[1:, 1:],
y[1:, :-1],
),
axis=-1,
)
_zbound = np.stack(
(
z[:-1, :-1],
z[:-1, 1:],
z[1:, 1:],
z[1:, :-1],
),
axis=-1,
)
# Compute centroid locations of each cells
self._xc = np.mean(_xbound, axis=-1)
self._yc = np.mean(_ybound, axis=-1)
self._zc = np.mean(_zbound, axis=-1)
self._zmax = 1.0

else:
# Boundaries of the hash grid are the bounding box of the source grid
self._xmin = self._source_grid.lon.min()
self._xmax = self._source_grid.lon.max()
self._ymin = self._source_grid.lat.min()
self._ymax = self._source_grid.lat.max()
# setting min and max below is needed for mesh="flat"
self._zmin = 0.0
self._zmax = 0.0
x = self._source_grid.lon
y = self._source_grid.lat

_xbound = np.stack(
(
x[:-1, :-1],
x[:-1, 1:],
x[1:, 1:],
x[1:, :-1],
),
axis=-1,
)
_ybound = np.stack(
(
y[:-1, :-1],
y[:-1, 1:],
y[1:, 1:],
y[1:, :-1],
),
axis=-1,
)
# Compute centroid locations of each cells
self._xc = np.mean(_xbound, axis=-1)
self._yc = np.mean(_ybound, axis=-1)
self._zc = np.zeros_like(self._xc)
# Here, we force _xc, _yc, _zc to be 2D arrays to
# mininimizes code change requirements between curvilinear and unstructured grids
self._xc = np.atleast_2d(self._source_grid.uxgrid.face_x.values)
self._yc = np.atleast_2d(self._source_grid.uxgrid.face_y.values)
self._zc = np.atleast_2d(self._source_grid.uxgrid.face_z.values)

# Generate the mapping from the hash indices to unstructured grid elements
self._hash_table = None
Expand Down Expand Up @@ -268,51 +281,6 @@ def query(
return (j_best.reshape(query_codes.shape), i_best.reshape(query_codes.shape))


def _triangle_area(A, B, C):
"""Compute the area of a triangle given by three points."""
d1 = B - A
d2 = C - A
d3 = np.cross(d1, d2)
return 0.5 * np.linalg.norm(d3)


def _barycentric_coordinates(nodes, point, min_area=1e-8):
"""
Compute the barycentric coordinates of a point P inside a convex polygon using area-based weights.
So that this method generalizes to n-sided polygons, we use the Waschpress points as the generalized
barycentric coordinates, which is only valid for convex polygons.

Parameters
----------
nodes : numpy.ndarray
Spherical coordinates (lat,lon) of each corner node of a face
point : numpy.ndarray
Spherical coordinates (lat,lon) of the point

Returns
-------
numpy.ndarray
Barycentric coordinates corresponding to each vertex.

"""
n = len(nodes)
sum_wi = 0
w = []

for i in range(0, n):
vim1 = nodes[i - 1]
vi = nodes[i]
vi1 = nodes[(i + 1) % n]
a0 = _triangle_area(vim1, vi, vi1)
a1 = max(_triangle_area(point, vim1, vi), min_area)
a2 = max(_triangle_area(point, vi, vi1), min_area)
sum_wi += a0 / (a1 * a2)
w.append(a0 / (a1 * a2))
barycentric_coords = [w_i / sum_wi for w_i in w]

return barycentric_coords


def _latlon_rad_to_xyz(lat, lon):
"""Converts Spherical latitude and longitude coordinates into Cartesian x,
y, z coordinates.
Expand Down
Loading
Loading