From 0c4fa4549179a69a55b1ed423733478e6f71da56 Mon Sep 17 00:00:00 2001 From: Maureen Cohen Date: Tue, 7 Jul 2026 14:08:54 +0100 Subject: [PATCH 01/24] Added SphericalMesh object to new mesh.py module. Contains only planetary radius at this time. --- src/parcels/_core/mesh.py | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) create mode 100644 src/parcels/_core/mesh.py diff --git a/src/parcels/_core/mesh.py b/src/parcels/_core/mesh.py new file mode 100644 index 000000000..7a48bb3d8 --- /dev/null +++ b/src/parcels/_core/mesh.py @@ -0,0 +1,24 @@ +import numpy as np + +class SphericalMesh: + """ Spherical mesh object with configurable planetary radius. + + Pass to FieldSet object as ``mesh=SphericalMesh(radius=...)``. + radius is in meters; None reverts degree to meter conversion + to 1852 * 60 .""" + + def __init__(self, radius: float | None = None): + self.radius = radius + + @property + def deg2m(self) -> float: + """ Meters per degree of arc.""" + if self.radius is None: + return 1852 * 60.0 + else: + return self.radius * np.pi / 180.0 + + def __repr__(self) -> str: + return f"SphericalMesh(radius={self.radius})" + + \ No newline at end of file From 3ddff6051bf5cc2bc605c09a699cf8631cbc4788 Mon Sep 17 00:00:00 2001 From: Maureen Cohen Date: Tue, 7 Jul 2026 14:10:41 +0100 Subject: [PATCH 02/24] Add SphericalMesh object to imports in __init__.py. --- src/parcels/__init__.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/parcels/__init__.py b/src/parcels/__init__.py index 78227247b..81c901ca9 100644 --- a/src/parcels/__init__.py +++ b/src/parcels/__init__.py @@ -22,6 +22,7 @@ from parcels._core.basegrid import BaseGrid from parcels._core.uxgrid import UxGrid from parcels._core.xgrid import XGrid +from parcels._core.mesh import SphericalMesh from parcels._core.statuscodes import ( AllParcelsErrorCodes, @@ -54,6 +55,7 @@ "BaseGrid", "UxGrid", "XGrid", + "SphericalMesh", # Status codes and errors "AllParcelsErrorCodes", "FieldInterpolationError", From ffa6d858ccaebfd086c8fb5f253a9e2e64597326 Mon Sep 17 00:00:00 2001 From: Maureen Cohen Date: Tue, 7 Jul 2026 14:28:24 +0100 Subject: [PATCH 03/24] Added SphericalMesh object unpacking to Xgrid. --- src/parcels/_core/xgrid.py | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/src/parcels/_core/xgrid.py b/src/parcels/_core/xgrid.py index 91a48af5d..95bc97fd1 100644 --- a/src/parcels/_core/xgrid.py +++ b/src/parcels/_core/xgrid.py @@ -14,6 +14,7 @@ from parcels._core.index_search import _search_1d_array, _search_indices_curvilinear_2d from parcels._sgrid.accessor import _get_dim_to_axis_mapping from parcels._sgrid.core import SGRID_PADDING_TO_XGCM_POSITION +from parcels._core.mesh import SphericalMesh _FIELD_DATA_ORDERING: Sequence[ptyping.XgcmAxisDirection] = "TZYX" _XGRID_AXES_ORDERING: Sequence[ptyping.XgridAxis] = "ZYX" @@ -169,7 +170,12 @@ def __init__(self, model_data: xr.Dataset, mesh): self._ds = model_data grid = XgcmLikeGrid(self.sgrid_metadata, model_data) self.xgcm_grid = grid - self._mesh = mesh + if isinstance(mesh, SphericalMesh): + self._mesh = "spherical" + self._radius = mesh.radius + else: + self._mesh = mesh + self._radius = mesh.radius self._spatialhash = None ds = model_data @@ -248,6 +254,14 @@ def _datetimes(self): @property def time(self): return self._datetimes.astype(np.float64) / 1e9 + + @property + def deg2m(self) -> float: + """ Metres per degree of arc for this grid's mesh. """ + if self._radius is None: + return 1852 * 60.0 + else: + return self._radius * np.pi / 180.0 @cached_property def xdim(self) -> int: From f566fe10c3221f83ca9171ea34f2c03a266ecf22 Mon Sep 17 00:00:00 2001 From: Maureen Cohen Date: Tue, 7 Jul 2026 14:33:13 +0100 Subject: [PATCH 04/24] Switched self._radius = mesh.radius to None in second branch of if/else block. --- src/parcels/_core/xgrid.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/parcels/_core/xgrid.py b/src/parcels/_core/xgrid.py index 95bc97fd1..d1ce68666 100644 --- a/src/parcels/_core/xgrid.py +++ b/src/parcels/_core/xgrid.py @@ -175,7 +175,7 @@ def __init__(self, model_data: xr.Dataset, mesh): self._radius = mesh.radius else: self._mesh = mesh - self._radius = mesh.radius + self._radius = None self._spatialhash = None ds = model_data From 69fa293e16139b4f1dc62e2a6c8d5d040f247203 Mon Sep 17 00:00:00 2001 From: Maureen Cohen Date: Tue, 7 Jul 2026 16:11:30 +0100 Subject: [PATCH 05/24] Applied SphericalMesh/radius imports to Uxgrid as well. --- src/parcels/_core/uxgrid.py | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/src/parcels/_core/uxgrid.py b/src/parcels/_core/uxgrid.py index 39201c2e8..970738bd5 100644 --- a/src/parcels/_core/uxgrid.py +++ b/src/parcels/_core/uxgrid.py @@ -7,6 +7,7 @@ from parcels._core.basegrid import BaseGrid from parcels._core.index_search import GRID_SEARCH_ERROR, _search_1d_array, uxgrid_point_in_cell +from parcels._core.mesh import SphericalMesh from parcels._typing import assert_valid_mesh _UXGRID_AXES = Literal["Z", "FACE"] @@ -41,7 +42,12 @@ def __init__(self, grid: ux.grid.Grid, z: ux.UxDataArray, mesh) -> None: if z.ndim != 1: raise ValueError("z must be a 1D array of vertical coordinates") self.z = z - self._mesh = mesh + if isinstance(mesh, SphericalMesh): + self._mesh = "spherical" + self._radius = mesh.radius + else: + self._mesh = mesh + self._radius = None self._spatialhash = None assert_valid_mesh(mesh) @@ -72,6 +78,14 @@ def get_axis_dim(self, axis: _UXGRID_AXES) -> int: return len(self.z.values) elif axis == "FACE": return self.uxgrid.n_face + + @property + def deg2m(self) -> float: + """ Metres per degree of arc for this grid's mesh. """ + if self._radius is None: + return 1852 * 60.0 + else: + return self._radius * np.pi / 180.0 def search(self, z, y, x, ei=None, tol=1e-6): """ From 6ec7382c27ed99fffa4a3e95b94c2edc2ddc4156 Mon Sep 17 00:00:00 2001 From: Maureen Cohen Date: Tue, 7 Jul 2026 16:14:19 +0100 Subject: [PATCH 06/24] Adjust assert_valid_mesh in typing to allow SphericalMesh to pass. --- src/parcels/_typing.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/parcels/_typing.py b/src/parcels/_typing.py index e8993cb05..57f84bb60 100644 --- a/src/parcels/_typing.py +++ b/src/parcels/_typing.py @@ -73,4 +73,7 @@ def _validate_against_pure_literal(value, typing_literal): # Assertion functions to clean user input def assert_valid_mesh(value: Any): + from parcels._core.mesh import SphericalMesh + if isinstance(value, SphericalMesh) + return _validate_against_pure_literal(value, Mesh) From b10e32c32ebcaca7e6c9149a325ed191a1d1bf13 Mon Sep 17 00:00:00 2001 From: Maureen Cohen Date: Tue, 7 Jul 2026 16:16:11 +0100 Subject: [PATCH 07/24] Bugfix: missed colon in if isinstance(value, SphericalMesh) statement. --- src/parcels/_typing.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/parcels/_typing.py b/src/parcels/_typing.py index 57f84bb60..9c32aa111 100644 --- a/src/parcels/_typing.py +++ b/src/parcels/_typing.py @@ -74,6 +74,6 @@ def _validate_against_pure_literal(value, typing_literal): # Assertion functions to clean user input def assert_valid_mesh(value: Any): from parcels._core.mesh import SphericalMesh - if isinstance(value, SphericalMesh) + if isinstance(value, SphericalMesh): return _validate_against_pure_literal(value, Mesh) From 493ef8fd5974a46e8a537829e3b8420d34ee47f7 Mon Sep 17 00:00:00 2001 From: Maureen Cohen Date: Tue, 7 Jul 2026 16:19:54 +0100 Subject: [PATCH 08/24] Moved SphericalMesh import to top rather than within function. --- src/parcels/_typing.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/parcels/_typing.py b/src/parcels/_typing.py index 9c32aa111..779fcdffb 100644 --- a/src/parcels/_typing.py +++ b/src/parcels/_typing.py @@ -10,6 +10,7 @@ from collections.abc import Callable, Mapping from datetime import datetime from typing import TYPE_CHECKING, Any, Literal, get_args +from parcels._core.mesh import SphericalMesh import numpy as np from cftime import datetime as cftime_datetime @@ -73,7 +74,6 @@ def _validate_against_pure_literal(value, typing_literal): # Assertion functions to clean user input def assert_valid_mesh(value: Any): - from parcels._core.mesh import SphericalMesh if isinstance(value, SphericalMesh): return _validate_against_pure_literal(value, Mesh) From 08fec2c4d72f4b36d0b3d3ab960f5d1477c99ad7 Mon Sep 17 00:00:00 2001 From: Maureen Cohen Date: Tue, 7 Jul 2026 16:52:12 +0100 Subject: [PATCH 09/24] Replaced hard-coded instances of 1852 * 60 with calls to deg2m attribute. --- src/parcels/_core/kernel.py | 2 +- src/parcels/_core/utils/interpolation.py | 12 +++---- src/parcels/interpolators/_uxinterpolators.py | 4 +-- src/parcels/interpolators/_xinterpolators.py | 16 ++++----- src/parcels/kernels/_advection.py | 10 +++--- src/parcels/kernels/_advectiondiffusion.py | 36 +++++++++---------- 6 files changed, 40 insertions(+), 40 deletions(-) diff --git a/src/parcels/_core/kernel.py b/src/parcels/_core/kernel.py index f50c2876e..6a526461f 100644 --- a/src/parcels/_core/kernel.py +++ b/src/parcels/_core/kernel.py @@ -143,7 +143,7 @@ def check_fieldsets_in_kernels(self, kernel): # TODO v4: this can go into anoth self.fieldset.add_context("RK45_tol", 10) if self.fieldset.U.grid._mesh == "spherical": self.fieldset.RK45_tol /= ( - 1852 * 60 + self.fieldset.U.grid.deg2m ) # TODO does not account for zonal variation in meter -> degree conversion if not hasattr(self.fieldset, "RK45_min_dt"): warnings.warn( diff --git a/src/parcels/_core/utils/interpolation.py b/src/parcels/_core/utils/interpolation.py index a48c0ff21..f8c477500 100644 --- a/src/parcels/_core/utils/interpolation.py +++ b/src/parcels/_core/utils/interpolation.py @@ -61,12 +61,12 @@ def dphidxsi3D_lin(zeta: float, eta: float, xsi: float) -> tuple[list[float], li def dxdxsi3D_lin( - hexa_z: list[float], hexa_y: list[float], hexa_x: list[float], zeta: float, eta: float, xsi: float, mesh: Mesh + hexa_z: list[float], hexa_y: list[float], hexa_x: list[float], zeta: float, eta: float, xsi: float, mesh: Mesh, + deg2m: float = 1852 * 60.0 ) -> tuple[float, float, float, float, float, float, float, float, float]: dphidxsi, dphideta, dphidzet = dphidxsi3D_lin(zeta, eta, xsi) if mesh == 'spherical': - deg2m = 1852 * 60. rad = np.pi / 180. lat = (1-xsi) * (1-eta) * hexa_y[0] + \ xsi * (1-eta) * hexa_y[1] + \ @@ -92,9 +92,10 @@ def dxdxsi3D_lin( def jacobian3D_lin( - hexa_z: list[float], hexa_y: list[float], hexa_x: list[float], zeta: float, eta: float, xsi: float, mesh: Mesh + hexa_z: list[float], hexa_y: list[float], hexa_x: list[float], zeta: float, eta: float, xsi: float, mesh: Mesh, + deg2m: float = 1852 * 60.0 ) -> float: - dxdxsi, dxdeta, dxdzet, dydxsi, dydeta, dydzet, dzdxsi, dzdeta, dzdzet = dxdxsi3D_lin(hexa_z, hexa_y, hexa_x, zeta, eta, xsi, mesh) + dxdxsi, dxdeta, dxdzet, dydxsi, dydeta, dydzet, dzdxsi, dzdeta, dzdzet = dxdxsi3D_lin(hexa_z, hexa_y, hexa_x, zeta, eta, xsi, mesh, deg2m) jac = ( dxdxsi * (dydeta * dzdzet - dzdeta * dydzet) @@ -174,10 +175,9 @@ def interpolate(phi: Callable[[float], list[float]], f: list[float], xsi: float) return np.dot(phi(xsi), f) -def _geodetic_distance(lat1: float, lat2: float, lon1: float, lon2: float, mesh: Mesh, lat: float) -> float: +def _geodetic_distance(lat1: float, lat2: float, lon1: float, lon2: float, mesh: Mesh, lat: float, deg2m: float = 1852 * 60.0) -> float: if mesh == "spherical": rad = np.pi / 180.0 - deg2m = 1852 * 60.0 return np.sqrt(((lon2 - lon1) * deg2m * np.cos(rad * lat)) ** 2 + ((lat2 - lat1) * deg2m) ** 2) else: return np.sqrt((lon2 - lon1) ** 2 + (lat2 - lat1) ** 2) diff --git a/src/parcels/interpolators/_uxinterpolators.py b/src/parcels/interpolators/_uxinterpolators.py index 80e804475..0b04b1de6 100644 --- a/src/parcels/interpolators/_uxinterpolators.py +++ b/src/parcels/interpolators/_uxinterpolators.py @@ -171,8 +171,8 @@ def interp( u = vectorfield.U.interp_method.interp(particle_positions, grid_positions, vectorfield.U) v = vectorfield.V.interp_method.interp(particle_positions, grid_positions, vectorfield.V) if vectorfield.grid._mesh == "spherical": - u /= 1852 * 60 * np.cos(np.deg2rad(particle_positions["y"])) - v /= 1852 * 60 + u /= vectorfield.grid.deg2m * np.cos(np.deg2rad(particle_positions["y"])) + v /= vectorfield.grid.deg2m if "3D" in vectorfield.vector_type: w = vectorfield.W.interp_method.interp(particle_positions, grid_positions, vectorfield.W) diff --git a/src/parcels/interpolators/_xinterpolators.py b/src/parcels/interpolators/_xinterpolators.py index d20313ffb..d851e42c5 100644 --- a/src/parcels/interpolators/_xinterpolators.py +++ b/src/parcels/interpolators/_xinterpolators.py @@ -149,8 +149,8 @@ def interp( u = _xlinear.interp(particle_positions, grid_positions, vectorfield.U) v = _xlinear.interp(particle_positions, grid_positions, vectorfield.V) if vectorfield.grid._mesh == "spherical": - u /= 1852 * 60 * np.cos(np.deg2rad(particle_positions["y"])) - v /= 1852 * 60 + u /= vectorfield.grid.deg2m * np.cos(np.deg2rad(particle_positions["y"])) + v /= vectorfield.grid.deg2m if vectorfield.W: w = _xlinear.interp(particle_positions, grid_positions, vectorfield.W) @@ -201,16 +201,16 @@ def interp( px[1:] = np.where(px[1:] - px[0] > 180, px[1:] - 360, px[1:]) px[1:] = np.where(-px[1:] + px[0] > 180, px[1:] + 360, px[1:]) c1 = i_u._geodetic_distance( - py[0], py[1], px[0], px[1], grid._mesh, np.einsum("ij,ji->i", i_u.phi2D_lin(0.0, xsi), py) + py[0], py[1], px[0], px[1], grid._mesh, np.einsum("ij,ji->i", i_u.phi2D_lin(0.0, xsi), py), grid.deg2m ) c2 = i_u._geodetic_distance( - py[1], py[2], px[1], px[2], grid._mesh, np.einsum("ij,ji->i", i_u.phi2D_lin(eta, 1.0), py) + py[1], py[2], px[1], px[2], grid._mesh, np.einsum("ij,ji->i", i_u.phi2D_lin(eta, 1.0), py), grid.deg2m ) c3 = i_u._geodetic_distance( - py[2], py[3], px[2], px[3], grid._mesh, np.einsum("ij,ji->i", i_u.phi2D_lin(1.0, xsi), py) + py[2], py[3], px[2], px[3], grid._mesh, np.einsum("ij,ji->i", i_u.phi2D_lin(1.0, xsi), py), grid.deg2m ) c4 = i_u._geodetic_distance( - py[3], py[0], px[3], px[0], grid._mesh, np.einsum("ij,ji->i", i_u.phi2D_lin(eta, 0.0), py) + py[3], py[0], px[3], px[0], grid._mesh, np.einsum("ij,ji->i", i_u.phi2D_lin(eta, 0.0), py), grid.deg2m ) def _create_selection_dict(dims, zdir=False): @@ -283,7 +283,7 @@ def _compute_corner_data(data, selection_dict) -> np.ndarray: Vvel = (1 - eta) * V0 + eta * V1 if grid._mesh == "spherical": - jac = i_u._compute_jacobian_determinant(py, px, eta, xsi) * 1852 * 60.0 + jac = i_u._compute_jacobian_determinant(py, px, eta, xsi) * grid.deg2m else: jac = i_u._compute_jacobian_determinant(py, px, eta, xsi) @@ -304,7 +304,7 @@ def _compute_corner_data(data, selection_dict) -> np.ndarray: v = v.compute() if grid._mesh == "spherical": - conversion = 1852 * 60.0 * np.cos(np.deg2rad(particle_positions["y"])) + conversion = grid.deg2m * np.cos(np.deg2rad(particle_positions["y"])) u /= conversion v /= conversion diff --git a/src/parcels/kernels/_advection.py b/src/parcels/kernels/_advection.py index 10e7ee676..bcf133361 100644 --- a/src/parcels/kernels/_advection.py +++ b/src/parcels/kernels/_advection.py @@ -243,12 +243,12 @@ def AdvectionAnalytical(particles, fieldset): # pragma: no cover else: dz = 1.0 - c1 = i_u._geodetic_distance(py[0], py[1], px[0], px[1], grid.mesh, np.dot(i_u.phi2D_lin(0.0, xsi), py)) - c2 = i_u._geodetic_distance(py[1], py[2], px[1], px[2], grid.mesh, np.dot(i_u.phi2D_lin(eta, 1.0), py)) - c3 = i_u._geodetic_distance(py[2], py[3], px[2], px[3], grid.mesh, np.dot(i_u.phi2D_lin(1.0, xsi), py)) - c4 = i_u._geodetic_distance(py[3], py[0], px[3], px[0], grid.mesh, np.dot(i_u.phi2D_lin(eta, 0.0), py)) + c1 = i_u._geodetic_distance(py[0], py[1], px[0], px[1], grid.mesh, np.dot(i_u.phi2D_lin(0.0, xsi), py), grid.deg2m) + c2 = i_u._geodetic_distance(py[1], py[2], px[1], px[2], grid.mesh, np.dot(i_u.phi2D_lin(eta, 1.0), py), grid.deg2m) + c3 = i_u._geodetic_distance(py[2], py[3], px[2], px[3], grid.mesh, np.dot(i_u.phi2D_lin(1.0, xsi), py), grid.deg2m) + c4 = i_u._geodetic_distance(py[3], py[0], px[3], px[0], grid.mesh, np.dot(i_u.phi2D_lin(eta, 0.0), py), grid.deg2m) rad = np.pi / 180.0 - deg2m = 1852 * 60.0 + deg2m = grid.deg2m meshJac = (deg2m * deg2m * math.cos(rad * particles.y)) if grid.mesh == "spherical" else 1 dxdy = i_u._compute_jacobian_determinant(py, px, eta, xsi) * meshJac diff --git a/src/parcels/kernels/_advectiondiffusion.py b/src/parcels/kernels/_advectiondiffusion.py index 3593e2538..7b7a96276 100644 --- a/src/parcels/kernels/_advectiondiffusion.py +++ b/src/parcels/kernels/_advectiondiffusion.py @@ -8,14 +8,14 @@ __all__ = ["AdvectionDiffusionEM", "AdvectionDiffusionM1", "DiffusionUniformKh"] -def meters_to_degrees_zonal(deg, lat): # pragma: no cover +def meters_to_degrees_zonal(deg, lat, deg2m): # pragma: no cover """Convert square meters to square degrees longitude at a given latitude.""" - return deg / pow(1852 * 60.0 * np.cos(lat * np.pi / 180), 2) + return deg / pow(deg2m * np.cos(lat * np.pi / 180), 2) -def meters_to_degrees_meridional(deg): # pragma: no cover +def meters_to_degrees_meridional(deg, deg2m): # pragma: no cover """Convert square meters to square degrees latitude.""" - return deg / pow(1852 * 60.0, 2) + return deg / pow(deg2m, 2) def AdvectionDiffusionM1(particles, fieldset): # pragma: no cover @@ -40,26 +40,26 @@ def AdvectionDiffusionM1(particles, fieldset): # pragma: no cover Kxp1 = fieldset.Kh_zonal[particles.t, particles.z, particles.y, particles.x + fieldset.dres, particles] Kxm1 = fieldset.Kh_zonal[particles.t, particles.z, particles.y, particles.x - fieldset.dres, particles] if fieldset.Kh_zonal.grid._mesh == "spherical": - Kxp1 = meters_to_degrees_zonal(Kxp1, particles.y) - Kxm1 = meters_to_degrees_zonal(Kxm1, particles.y) + Kxp1 = meters_to_degrees_zonal(Kxp1, particles.y, fieldset.Kh_zonal.grid.deg2m) + Kxm1 = meters_to_degrees_zonal(Kxm1, particles.y, fieldset.Kh_zonal.grid.deg2m) dKdx = (Kxp1 - Kxm1) / (2 * fieldset.dres) u, v = fieldset.UV[particles.t, particles.z, particles.y, particles.x, particles] kh_zonal = fieldset.Kh_zonal[particles.t, particles.z, particles.y, particles.x, particles] if fieldset.Kh_zonal.grid._mesh == "spherical": - kh_zonal = meters_to_degrees_zonal(kh_zonal, particles.y) + kh_zonal = meters_to_degrees_zonal(kh_zonal, particles.y, fieldset.Kh_zonal.grid.deg2m) bx = np.sqrt(2 * kh_zonal) Kyp1 = fieldset.Kh_meridional[particles.t, particles.z, particles.y + fieldset.dres, particles.x, particles] Kym1 = fieldset.Kh_meridional[particles.t, particles.z, particles.y - fieldset.dres, particles.x, particles] if fieldset.Kh_meridional.grid._mesh == "spherical": - Kyp1 = meters_to_degrees_meridional(Kyp1) - Kym1 = meters_to_degrees_meridional(Kym1) + Kyp1 = meters_to_degrees_meridional(Kyp1, fieldset.Kh_meridional.grid.deg2m) + Kym1 = meters_to_degrees_meridional(Kym1, fieldset.Kh_meridional.grid.deg2m) dKdy = (Kyp1 - Kym1) / (2 * fieldset.dres) kh_meridional = fieldset.Kh_meridional[particles.t, particles.z, particles.y, particles.x, particles] if fieldset.Kh_meridional.grid._mesh == "spherical": - kh_meridional = meters_to_degrees_meridional(kh_meridional) + kh_meridional = meters_to_degrees_meridional(kh_meridional, fieldset.Kh_meridional.grid.deg2m) by = np.sqrt(2 * kh_meridional) # Particle positions are updated only after evaluating all terms. @@ -89,27 +89,27 @@ def AdvectionDiffusionEM(particles, fieldset): # pragma: no cover Kxp1 = fieldset.Kh_zonal[particles.t, particles.z, particles.y, particles.x + fieldset.dres, particles] Kxm1 = fieldset.Kh_zonal[particles.t, particles.z, particles.y, particles.x - fieldset.dres, particles] if fieldset.Kh_zonal.grid._mesh == "spherical": - Kxp1 = meters_to_degrees_zonal(Kxp1, particles.y) - Kxm1 = meters_to_degrees_zonal(Kxm1, particles.y) + Kxp1 = meters_to_degrees_zonal(Kxp1, particles.y, fieldset.Kh_zonal.grid.deg2m) + Kxm1 = meters_to_degrees_zonal(Kxm1, particles.y, fieldset.Kh_zonal.grid.deg2m) dKdx = (Kxp1 - Kxm1) / (2 * fieldset.dres) ax = u + dKdx kh_zonal = fieldset.Kh_zonal[particles.t, particles.z, particles.y, particles.x, particles] if fieldset.Kh_zonal.grid._mesh == "spherical": - kh_zonal = meters_to_degrees_zonal(kh_zonal, particles.y) + kh_zonal = meters_to_degrees_zonal(kh_zonal, particles.y, fieldset.Kh_zonal.grid.deg2m) bx = np.sqrt(2 * kh_zonal) Kyp1 = fieldset.Kh_meridional[particles.t, particles.z, particles.y + fieldset.dres, particles.x, particles] Kym1 = fieldset.Kh_meridional[particles.t, particles.z, particles.y - fieldset.dres, particles.x, particles] if fieldset.Kh_meridional.grid._mesh == "spherical": - Kyp1 = meters_to_degrees_meridional(Kyp1) - Kym1 = meters_to_degrees_meridional(Kym1) + Kyp1 = meters_to_degrees_meridional(Kyp1, fieldset.Kh_meridional.grid.deg2m) + Kym1 = meters_to_degrees_meridional(Kym1, fieldset.Kh_meridional.grid.deg2m) dKdy = (Kyp1 - Kym1) / (2 * fieldset.dres) ay = v + dKdy kh_meridional = fieldset.Kh_meridional[particles.t, particles.z, particles.y, particles.x, particles] if fieldset.Kh_meridional.grid._mesh == "spherical": - kh_meridional = meters_to_degrees_meridional(kh_meridional) + kh_meridional = meters_to_degrees_meridional(kh_meridional, fieldset.Kh_meridional.grid.deg2m) by = np.sqrt(2 * kh_meridional) # Particle positions are updated only after evaluating all terms. @@ -143,8 +143,8 @@ def DiffusionUniformKh(particles, fieldset): # pragma: no cover kh_meridional = fieldset.Kh_meridional[particles] if fieldset.Kh_zonal.grid._mesh == "spherical": - kh_zonal = meters_to_degrees_zonal(kh_zonal, particles.y) - kh_meridional = meters_to_degrees_meridional(kh_meridional) + kh_zonal = meters_to_degrees_zonal(kh_zonal, particles.y, fieldset.Kh_zonal.grid.deg2m) + kh_meridional = meters_to_degrees_meridional(kh_meridional, fieldset.Kh_meridional.grid.deg2m) bx = np.sqrt(2 * kh_zonal) by = np.sqrt(2 * kh_meridional) From a2779136791019b0a5a6e9221d43a1f980b7e458 Mon Sep 17 00:00:00 2001 From: Maureen Cohen Date: Tue, 7 Jul 2026 17:37:52 +0100 Subject: [PATCH 10/24] Added three tests: test that deg2m calculation works, test calculation for different input args (string vs SphericalMesh object, radius given or not), short advection test with radii None, Mars, Venus, Earth. --- tests/test_mesh.py | 44 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 44 insertions(+) create mode 100644 tests/test_mesh.py diff --git a/tests/test_mesh.py b/tests/test_mesh.py new file mode 100644 index 000000000..df4b24412 --- /dev/null +++ b/tests/test_mesh.py @@ -0,0 +1,44 @@ +import numpy as np +import pytest + +from parcels import FieldSet, ParticleSet, SphericalMesh +from parcels.kernels import AdvectionRK4 +from parcels._datasets.structured.generated import simple_UV_dataset + +EARTH_DEG2M = 1852 * 60.0 + +def test_spherical_mesh_deg2m(): + assert SphericalMesh().radius is None + assert SphericalMesh().deg2m == EARTH_DEG2M + r = 3389500.0 # Mars radius + assert SphericalMesh(radius=r).deg2m == pytest.approx(r * np.pi / 180) + +@pytest.mark.parametrize("mesh, exp_radius, exp_deg2m", + [ + ("spherical", None, EARTH_DEG2M), + (SphericalMesh(), None, EARTH_DEG2M), + (SphericalMesh(radius=3389500.0), 3389500.0, 3389500.0 * np.pi / 180), + ], + ) +def test_xgrid_radius_and_deg2m(mesh, exp_radius, exp_deg2m): + grid = FieldSet.from_sgrid_conventions(simple_UV_dataset(), mesh=mesh).U.grid + assert grid._mesh == "spherical" + assert grid._radius == exp_radius + assert grid.deg2m == pytest.approx(exp_deg2m) + +@pytest.mark.parametrize("radius", [None, 3389500.0, 6051800.0, 6371000.0]) # Mars, Venus, Earth +def test_advection_uses_custom_radius(radius, npart=10): + ds = simple_UV_dataset() + ds["U"].data[:] = 1.0 + fieldset = FieldSet.from_sgrid_conventions(ds, mesh=SphericalMesh(radius=radius)) + + runtime = 7200 + startlat = np.linspace(0, 80, npart) + startlon = 20.0 + np.zeros(npart) + pset = ParticleSet(fieldset, x=startlon, y=startlat) + pset.execute(AdvectionRK4, runtime=runtime, dt=np.timedelta64(15, "m")) + + deg2m = EARTH_DEG2M if radius is None else radius * np.pi / 180 + expected_dlon = runtime / (deg2m * np.cos(np.deg2rad(pset.y))) + np.testing.assert_allclose(pset.x - startlon, expected_dlon, atol=1e-5) + np.testing.assert_allclose(pset.y, startlat, atol=1e-5) From da32a50635169fe041c18ffc442150c8b15f82ab Mon Sep 17 00:00:00 2001 From: Maureen Cohen Date: Tue, 7 Jul 2026 17:53:21 +0100 Subject: [PATCH 11/24] Added a check to SphericalMesh init to catch bad values of radius: nonon-numeric, negative, or 0. Added matching tests to test_mesh. Co-authored by: Claude ). --- src/parcels/_core/mesh.py | 5 ++++- tests/test_mesh.py | 10 ++++++++++ 2 files changed, 14 insertions(+), 1 deletion(-) diff --git a/src/parcels/_core/mesh.py b/src/parcels/_core/mesh.py index 7a48bb3d8..334674613 100644 --- a/src/parcels/_core/mesh.py +++ b/src/parcels/_core/mesh.py @@ -8,6 +8,10 @@ class SphericalMesh: to 1852 * 60 .""" def __init__(self, radius: float | None = None): + if radius is not None and not isinstance(radius, (int, float, np.number)): + raise TypeError(f"radius must be a number of None, got {type(radius).__name__}") + if radius is not None and radius <= 0: + raise ValueError(f"radius must be positive, got {radius}") self.radius = radius @property @@ -21,4 +25,3 @@ def deg2m(self) -> float: def __repr__(self) -> str: return f"SphericalMesh(radius={self.radius})" - \ No newline at end of file diff --git a/tests/test_mesh.py b/tests/test_mesh.py index df4b24412..62baa490b 100644 --- a/tests/test_mesh.py +++ b/tests/test_mesh.py @@ -42,3 +42,13 @@ def test_advection_uses_custom_radius(radius, npart=10): expected_dlon = runtime / (deg2m * np.cos(np.deg2rad(pset.y))) np.testing.assert_allclose(pset.x - startlon, expected_dlon, atol=1e-5) np.testing.assert_allclose(pset.y, startlat, atol=1e-5) + +@pytest.mark.parametrize("bad_radius", ["6371000", [6371000], (1, 2), {}]) +def test_spherical_mesh_rejects_non_numeric_radius(bad_radius): + with pytest.raises(TypeError): + SphericalMesh(radius=bad_radius) + +@pytest.mark.parametrize("bad_radius", [0, -1.0, -6371000]) +def test_spherical_mesh_rejects_nonpos_radius(bad_radius): + with pytest.raises(ValueError): + SphericalMesh(radius=bad_radius) \ No newline at end of file From 72be313963f4e369d37272923db1c2048515411e Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 9 Jul 2026 16:11:41 +0000 Subject: [PATCH 12/24] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/parcels/_core/mesh.py | 21 +++++++++--------- src/parcels/_core/utils/interpolation.py | 4 +++- src/parcels/_core/uxgrid.py | 4 ++-- src/parcels/_core/xgrid.py | 6 ++--- src/parcels/_typing.py | 3 ++- tests/test_mesh.py | 28 ++++++++++++++---------- 6 files changed, 38 insertions(+), 28 deletions(-) diff --git a/src/parcels/_core/mesh.py b/src/parcels/_core/mesh.py index 334674613..b9ca044af 100644 --- a/src/parcels/_core/mesh.py +++ b/src/parcels/_core/mesh.py @@ -1,12 +1,14 @@ import numpy as np + class SphericalMesh: - """ Spherical mesh object with configurable planetary radius. - - Pass to FieldSet object as ``mesh=SphericalMesh(radius=...)``. - radius is in meters; None reverts degree to meter conversion - to 1852 * 60 .""" - + """Spherical mesh object with configurable planetary radius. + + Pass to FieldSet object as ``mesh=SphericalMesh(radius=...)``. + radius is in meters; None reverts degree to meter conversion + to 1852 * 60 . + """ + def __init__(self, radius: float | None = None): if radius is not None and not isinstance(radius, (int, float, np.number)): raise TypeError(f"radius must be a number of None, got {type(radius).__name__}") @@ -15,13 +17,12 @@ def __init__(self, radius: float | None = None): self.radius = radius @property - def deg2m(self) -> float: - """ Meters per degree of arc.""" + def deg2m(self) -> float: + """Meters per degree of arc.""" if self.radius is None: return 1852 * 60.0 else: return self.radius * np.pi / 180.0 - + def __repr__(self) -> str: return f"SphericalMesh(radius={self.radius})" - diff --git a/src/parcels/_core/utils/interpolation.py b/src/parcels/_core/utils/interpolation.py index f8c477500..214ba1933 100644 --- a/src/parcels/_core/utils/interpolation.py +++ b/src/parcels/_core/utils/interpolation.py @@ -175,7 +175,9 @@ def interpolate(phi: Callable[[float], list[float]], f: list[float], xsi: float) return np.dot(phi(xsi), f) -def _geodetic_distance(lat1: float, lat2: float, lon1: float, lon2: float, mesh: Mesh, lat: float, deg2m: float = 1852 * 60.0) -> float: +def _geodetic_distance( + lat1: float, lat2: float, lon1: float, lon2: float, mesh: Mesh, lat: float, deg2m: float = 1852 * 60.0 +) -> float: if mesh == "spherical": rad = np.pi / 180.0 return np.sqrt(((lon2 - lon1) * deg2m * np.cos(rad * lat)) ** 2 + ((lat2 - lat1) * deg2m) ** 2) diff --git a/src/parcels/_core/uxgrid.py b/src/parcels/_core/uxgrid.py index 970738bd5..8cfb41ec5 100644 --- a/src/parcels/_core/uxgrid.py +++ b/src/parcels/_core/uxgrid.py @@ -78,10 +78,10 @@ def get_axis_dim(self, axis: _UXGRID_AXES) -> int: return len(self.z.values) elif axis == "FACE": return self.uxgrid.n_face - + @property def deg2m(self) -> float: - """ Metres per degree of arc for this grid's mesh. """ + """Metres per degree of arc for this grid's mesh.""" if self._radius is None: return 1852 * 60.0 else: diff --git a/src/parcels/_core/xgrid.py b/src/parcels/_core/xgrid.py index d1ce68666..dbd0a3c0e 100644 --- a/src/parcels/_core/xgrid.py +++ b/src/parcels/_core/xgrid.py @@ -12,9 +12,9 @@ import parcels._typing as ptyping from parcels._core.basegrid import BaseGrid from parcels._core.index_search import _search_1d_array, _search_indices_curvilinear_2d +from parcels._core.mesh import SphericalMesh from parcels._sgrid.accessor import _get_dim_to_axis_mapping from parcels._sgrid.core import SGRID_PADDING_TO_XGCM_POSITION -from parcels._core.mesh import SphericalMesh _FIELD_DATA_ORDERING: Sequence[ptyping.XgcmAxisDirection] = "TZYX" _XGRID_AXES_ORDERING: Sequence[ptyping.XgridAxis] = "ZYX" @@ -254,10 +254,10 @@ def _datetimes(self): @property def time(self): return self._datetimes.astype(np.float64) / 1e9 - + @property def deg2m(self) -> float: - """ Metres per degree of arc for this grid's mesh. """ + """Metres per degree of arc for this grid's mesh.""" if self._radius is None: return 1852 * 60.0 else: diff --git a/src/parcels/_typing.py b/src/parcels/_typing.py index 779fcdffb..61fc28160 100644 --- a/src/parcels/_typing.py +++ b/src/parcels/_typing.py @@ -10,11 +10,12 @@ from collections.abc import Callable, Mapping from datetime import datetime from typing import TYPE_CHECKING, Any, Literal, get_args -from parcels._core.mesh import SphericalMesh import numpy as np from cftime import datetime as cftime_datetime +from parcels._core.mesh import SphericalMesh + if TYPE_CHECKING: import xgcm diff --git a/tests/test_mesh.py b/tests/test_mesh.py index 62baa490b..296b0291a 100644 --- a/tests/test_mesh.py +++ b/tests/test_mesh.py @@ -2,31 +2,35 @@ import pytest from parcels import FieldSet, ParticleSet, SphericalMesh -from parcels.kernels import AdvectionRK4 from parcels._datasets.structured.generated import simple_UV_dataset +from parcels.kernels import AdvectionRK4 EARTH_DEG2M = 1852 * 60.0 + def test_spherical_mesh_deg2m(): assert SphericalMesh().radius is None assert SphericalMesh().deg2m == EARTH_DEG2M - r = 3389500.0 # Mars radius + r = 3389500.0 # Mars radius assert SphericalMesh(radius=r).deg2m == pytest.approx(r * np.pi / 180) -@pytest.mark.parametrize("mesh, exp_radius, exp_deg2m", - [ - ("spherical", None, EARTH_DEG2M), - (SphericalMesh(), None, EARTH_DEG2M), - (SphericalMesh(radius=3389500.0), 3389500.0, 3389500.0 * np.pi / 180), - ], - ) + +@pytest.mark.parametrize( + "mesh, exp_radius, exp_deg2m", + [ + ("spherical", None, EARTH_DEG2M), + (SphericalMesh(), None, EARTH_DEG2M), + (SphericalMesh(radius=3389500.0), 3389500.0, 3389500.0 * np.pi / 180), + ], +) def test_xgrid_radius_and_deg2m(mesh, exp_radius, exp_deg2m): grid = FieldSet.from_sgrid_conventions(simple_UV_dataset(), mesh=mesh).U.grid assert grid._mesh == "spherical" assert grid._radius == exp_radius assert grid.deg2m == pytest.approx(exp_deg2m) -@pytest.mark.parametrize("radius", [None, 3389500.0, 6051800.0, 6371000.0]) # Mars, Venus, Earth + +@pytest.mark.parametrize("radius", [None, 3389500.0, 6051800.0, 6371000.0]) # Mars, Venus, Earth def test_advection_uses_custom_radius(radius, npart=10): ds = simple_UV_dataset() ds["U"].data[:] = 1.0 @@ -43,12 +47,14 @@ def test_advection_uses_custom_radius(radius, npart=10): np.testing.assert_allclose(pset.x - startlon, expected_dlon, atol=1e-5) np.testing.assert_allclose(pset.y, startlat, atol=1e-5) + @pytest.mark.parametrize("bad_radius", ["6371000", [6371000], (1, 2), {}]) def test_spherical_mesh_rejects_non_numeric_radius(bad_radius): with pytest.raises(TypeError): SphericalMesh(radius=bad_radius) + @pytest.mark.parametrize("bad_radius", [0, -1.0, -6371000]) def test_spherical_mesh_rejects_nonpos_radius(bad_radius): with pytest.raises(ValueError): - SphericalMesh(radius=bad_radius) \ No newline at end of file + SphericalMesh(radius=bad_radius) From 0c3c1f047713e099634f27bd836fc5fa86d4a561 Mon Sep 17 00:00:00 2001 From: Maureen Cohen <53348730+maureenjcohen@users.noreply.github.com> Date: Fri, 10 Jul 2026 10:30:15 +0100 Subject: [PATCH 13/24] Update src/parcels/_core/mesh.py Co-authored-by: Erik van Sebille --- src/parcels/_core/mesh.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/parcels/_core/mesh.py b/src/parcels/_core/mesh.py index b9ca044af..28f160b21 100644 --- a/src/parcels/_core/mesh.py +++ b/src/parcels/_core/mesh.py @@ -11,7 +11,7 @@ class SphericalMesh: def __init__(self, radius: float | None = None): if radius is not None and not isinstance(radius, (int, float, np.number)): - raise TypeError(f"radius must be a number of None, got {type(radius).__name__}") + raise TypeError(f"radius must be a number or None, got {type(radius).__name__}") if radius is not None and radius <= 0: raise ValueError(f"radius must be positive, got {radius}") self.radius = radius From ddac215d621fedf93257a2258bd4f446e4392fa2 Mon Sep 17 00:00:00 2001 From: Maureen Cohen <53348730+maureenjcohen@users.noreply.github.com> Date: Fri, 10 Jul 2026 10:30:40 +0100 Subject: [PATCH 14/24] Update src/parcels/_core/mesh.py Co-authored-by: Erik van Sebille --- src/parcels/_core/mesh.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/parcels/_core/mesh.py b/src/parcels/_core/mesh.py index 28f160b21..fc9f10147 100644 --- a/src/parcels/_core/mesh.py +++ b/src/parcels/_core/mesh.py @@ -5,8 +5,9 @@ class SphericalMesh: """Spherical mesh object with configurable planetary radius. Pass to FieldSet object as ``mesh=SphericalMesh(radius=...)``. - radius is in meters; None reverts degree to meter conversion - to 1852 * 60 . + radius is in meters; None reverts to default for Earth, where + arcdegree to meter conversion is defined as 1852 * 60 + (1852 meters per arcminute * 60 arcminutes per arcdegree). """ def __init__(self, radius: float | None = None): From d28fa4c814a36e67b2f8db14118935e0cd6e1784 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 10 Jul 2026 09:30:49 +0000 Subject: [PATCH 15/24] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/parcels/_core/mesh.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/parcels/_core/mesh.py b/src/parcels/_core/mesh.py index fc9f10147..4bbd96aad 100644 --- a/src/parcels/_core/mesh.py +++ b/src/parcels/_core/mesh.py @@ -5,7 +5,7 @@ class SphericalMesh: """Spherical mesh object with configurable planetary radius. Pass to FieldSet object as ``mesh=SphericalMesh(radius=...)``. - radius is in meters; None reverts to default for Earth, where + radius is in meters; None reverts to default for Earth, where arcdegree to meter conversion is defined as 1852 * 60 (1852 meters per arcminute * 60 arcminutes per arcdegree). """ From ce76ee3a3ad0358d5e1808b5e35aae7b84a252d5 Mon Sep 17 00:00:00 2001 From: Maureen Cohen <53348730+maureenjcohen@users.noreply.github.com> Date: Fri, 10 Jul 2026 10:31:11 +0100 Subject: [PATCH 16/24] Update src/parcels/_core/uxgrid.py Co-authored-by: Erik van Sebille --- src/parcels/_core/uxgrid.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/parcels/_core/uxgrid.py b/src/parcels/_core/uxgrid.py index 8cfb41ec5..8cd86508e 100644 --- a/src/parcels/_core/uxgrid.py +++ b/src/parcels/_core/uxgrid.py @@ -81,7 +81,7 @@ def get_axis_dim(self, axis: _UXGRID_AXES) -> int: @property def deg2m(self) -> float: - """Metres per degree of arc for this grid's mesh.""" + """Metres per arcdegree for this grid's mesh.""" if self._radius is None: return 1852 * 60.0 else: From f781eebe849447073c7ca97bf09072af481f00fa Mon Sep 17 00:00:00 2001 From: Maureen Cohen Date: Mon, 13 Jul 2026 15:30:16 +0100 Subject: [PATCH 17/24] Modified radius setting so default is EARTH_RADIUS (specified in mesh.py and imported wherever else needed). Radius = None leads to a crash in flat mesh cases, so still kept an if/else to set deg2m = 1.0 in xgrid/uxgrid if no radius is specified, i.e. if SphericalMesh or spherical are not set. --- src/parcels/_core/mesh.py | 19 ++++++++--------- src/parcels/_core/uxgrid.py | 11 +++++----- src/parcels/_core/xgrid.py | 11 +++++----- tests/test_mesh.py | 41 ++++++++++++++++++++++++++----------- 4 files changed, 47 insertions(+), 35 deletions(-) diff --git a/src/parcels/_core/mesh.py b/src/parcels/_core/mesh.py index 4bbd96aad..10311a3d0 100644 --- a/src/parcels/_core/mesh.py +++ b/src/parcels/_core/mesh.py @@ -1,29 +1,26 @@ import numpy as np +EARTH_RADIUS = 6366707.019493707 + class SphericalMesh: """Spherical mesh object with configurable planetary radius. Pass to FieldSet object as ``mesh=SphericalMesh(radius=...)``. - radius is in meters; None reverts to default for Earth, where - arcdegree to meter conversion is defined as 1852 * 60 - (1852 meters per arcminute * 60 arcminutes per arcdegree). + radius is in meters; defaults to Earth radius. """ - def __init__(self, radius: float | None = None): - if radius is not None and not isinstance(radius, (int, float, np.number)): - raise TypeError(f"radius must be a number or None, got {type(radius).__name__}") - if radius is not None and radius <= 0: + def __init__(self, radius: float = EARTH_RADIUS): + if not isinstance(radius, (int, float, np.number)): + raise TypeError(f"radius must be a number, got {type(radius).__name__}") + if radius <= 0: raise ValueError(f"radius must be positive, got {radius}") self.radius = radius @property def deg2m(self) -> float: """Meters per degree of arc.""" - if self.radius is None: - return 1852 * 60.0 - else: - return self.radius * np.pi / 180.0 + return self.radius * np.pi / 180.0 def __repr__(self) -> str: return f"SphericalMesh(radius={self.radius})" diff --git a/src/parcels/_core/uxgrid.py b/src/parcels/_core/uxgrid.py index 8cd86508e..7caf79993 100644 --- a/src/parcels/_core/uxgrid.py +++ b/src/parcels/_core/uxgrid.py @@ -7,7 +7,7 @@ from parcels._core.basegrid import BaseGrid from parcels._core.index_search import GRID_SEARCH_ERROR, _search_1d_array, uxgrid_point_in_cell -from parcels._core.mesh import SphericalMesh +from parcels._core.mesh import SphericalMesh, EARTH_RADIUS from parcels._typing import assert_valid_mesh _UXGRID_AXES = Literal["Z", "FACE"] @@ -47,7 +47,7 @@ def __init__(self, grid: ux.grid.Grid, z: ux.UxDataArray, mesh) -> None: self._radius = mesh.radius else: self._mesh = mesh - self._radius = None + self._radius = EARTH_RADIUS if mesh == "spherical" else None self._spatialhash = None assert_valid_mesh(mesh) @@ -82,10 +82,9 @@ def get_axis_dim(self, axis: _UXGRID_AXES) -> int: @property def deg2m(self) -> float: """Metres per arcdegree for this grid's mesh.""" - if self._radius is None: - return 1852 * 60.0 - else: - return self._radius * np.pi / 180.0 + if self._radius is None: # flat mesh; None causes crash in advection + return 1.0 + return self._radius * np.pi / 180.0 def search(self, z, y, x, ei=None, tol=1e-6): """ diff --git a/src/parcels/_core/xgrid.py b/src/parcels/_core/xgrid.py index dbd0a3c0e..841e839af 100644 --- a/src/parcels/_core/xgrid.py +++ b/src/parcels/_core/xgrid.py @@ -12,7 +12,7 @@ import parcels._typing as ptyping from parcels._core.basegrid import BaseGrid from parcels._core.index_search import _search_1d_array, _search_indices_curvilinear_2d -from parcels._core.mesh import SphericalMesh +from parcels._core.mesh import SphericalMesh, EARTH_RADIUS from parcels._sgrid.accessor import _get_dim_to_axis_mapping from parcels._sgrid.core import SGRID_PADDING_TO_XGCM_POSITION @@ -175,7 +175,7 @@ def __init__(self, model_data: xr.Dataset, mesh): self._radius = mesh.radius else: self._mesh = mesh - self._radius = None + self._radius = EARTH_RADIUS if mesh == "spherical" else None self._spatialhash = None ds = model_data @@ -258,10 +258,9 @@ def time(self): @property def deg2m(self) -> float: """Metres per degree of arc for this grid's mesh.""" - if self._radius is None: - return 1852 * 60.0 - else: - return self._radius * np.pi / 180.0 + if self._radius is None: # flat mesh; None causes crash in advection + return 1.0 + return self._radius * np.pi / 180.0 @cached_property def xdim(self) -> int: diff --git a/tests/test_mesh.py b/tests/test_mesh.py index 296b0291a..2f42ad33f 100644 --- a/tests/test_mesh.py +++ b/tests/test_mesh.py @@ -2,24 +2,23 @@ import pytest from parcels import FieldSet, ParticleSet, SphericalMesh +from parcels._core.mesh import EARTH_RADIUS from parcels._datasets.structured.generated import simple_UV_dataset from parcels.kernels import AdvectionRK4 -EARTH_DEG2M = 1852 * 60.0 - +EARTH_DEG2M = EARTH_RADIUS * np.pi / 180 def test_spherical_mesh_deg2m(): - assert SphericalMesh().radius is None + assert SphericalMesh().radius == EARTH_RADIUS assert SphericalMesh().deg2m == EARTH_DEG2M r = 3389500.0 # Mars radius assert SphericalMesh(radius=r).deg2m == pytest.approx(r * np.pi / 180) - @pytest.mark.parametrize( "mesh, exp_radius, exp_deg2m", [ - ("spherical", None, EARTH_DEG2M), - (SphericalMesh(), None, EARTH_DEG2M), + ("spherical", EARTH_RADIUS, EARTH_DEG2M), + (SphericalMesh(), EARTH_RADIUS, EARTH_DEG2M), (SphericalMesh(radius=3389500.0), 3389500.0, 3389500.0 * np.pi / 180), ], ) @@ -29,12 +28,18 @@ def test_xgrid_radius_and_deg2m(mesh, exp_radius, exp_deg2m): assert grid._radius == exp_radius assert grid.deg2m == pytest.approx(exp_deg2m) - -@pytest.mark.parametrize("radius", [None, 3389500.0, 6051800.0, 6371000.0]) # Mars, Venus, Earth -def test_advection_uses_custom_radius(radius, npart=10): +@pytest.mark.parametrize("mesh, deg2m", + [ + (SphericalMesh(), EARTH_DEG2M), # No radius entered + (SphericalMesh(radius=3389500.0), 3389500.0 * np.pi / 180), # Mars + (SphericalMesh(radius=6051800.0), 6051800.0 * np.pi / 180), # Venus + (SphericalMesh(radius=EARTH_RADIUS), EARTH_DEG2M), # Explicit Earth + ], + ) +def test_advection_uses_custom_radius(mesh, deg2m, npart=10): ds = simple_UV_dataset() ds["U"].data[:] = 1.0 - fieldset = FieldSet.from_sgrid_conventions(ds, mesh=SphericalMesh(radius=radius)) + fieldset = FieldSet.from_sgrid_conventions(ds, mesh=mesh) runtime = 7200 startlat = np.linspace(0, 80, npart) @@ -42,18 +47,30 @@ def test_advection_uses_custom_radius(radius, npart=10): pset = ParticleSet(fieldset, x=startlon, y=startlat) pset.execute(AdvectionRK4, runtime=runtime, dt=np.timedelta64(15, "m")) - deg2m = EARTH_DEG2M if radius is None else radius * np.pi / 180 expected_dlon = runtime / (deg2m * np.cos(np.deg2rad(pset.y))) np.testing.assert_allclose(pset.x - startlon, expected_dlon, atol=1e-5) np.testing.assert_allclose(pset.y, startlat, atol=1e-5) +def test_advection_flat_mesh(npart=10): + ds = simple_UV_dataset(mesh="flat") + ds["U"].data[:] = 1.0 + fieldset = FieldSet.from_sgrid_conventions(ds, mesh="flat") + + runtime = 7200 + startlat = np.linspace(0, 80, npart) + startlon = 20.0 + np.zeros(npart) + pset = ParticleSet(fieldset, x=startlon, y=startlat) + pset.execute(AdvectionRK4, runtime=runtime, dt=np.timedelta64(15, "m")) + + assert fieldset.U.grid.deg2m == 1.0 # flat mesh deg2m + np.testing.assert_allclose(pset.x - startlon, runtime, atol=1e-5) + np.testing.assert_allclose(pset.y, startlat, atol=1e-5) @pytest.mark.parametrize("bad_radius", ["6371000", [6371000], (1, 2), {}]) def test_spherical_mesh_rejects_non_numeric_radius(bad_radius): with pytest.raises(TypeError): SphericalMesh(radius=bad_radius) - @pytest.mark.parametrize("bad_radius", [0, -1.0, -6371000]) def test_spherical_mesh_rejects_nonpos_radius(bad_radius): with pytest.raises(ValueError): From 9ac08fc8b44525d7f150d262dbf908723ed1d429 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 13 Jul 2026 14:31:38 +0000 Subject: [PATCH 18/24] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/parcels/_core/uxgrid.py | 4 ++-- src/parcels/_core/xgrid.py | 4 ++-- tests/test_mesh.py | 27 +++++++++++++++++---------- 3 files changed, 21 insertions(+), 14 deletions(-) diff --git a/src/parcels/_core/uxgrid.py b/src/parcels/_core/uxgrid.py index 7caf79993..05adca6d7 100644 --- a/src/parcels/_core/uxgrid.py +++ b/src/parcels/_core/uxgrid.py @@ -7,7 +7,7 @@ from parcels._core.basegrid import BaseGrid from parcels._core.index_search import GRID_SEARCH_ERROR, _search_1d_array, uxgrid_point_in_cell -from parcels._core.mesh import SphericalMesh, EARTH_RADIUS +from parcels._core.mesh import EARTH_RADIUS, SphericalMesh from parcels._typing import assert_valid_mesh _UXGRID_AXES = Literal["Z", "FACE"] @@ -82,7 +82,7 @@ def get_axis_dim(self, axis: _UXGRID_AXES) -> int: @property def deg2m(self) -> float: """Metres per arcdegree for this grid's mesh.""" - if self._radius is None: # flat mesh; None causes crash in advection + if self._radius is None: # flat mesh; None causes crash in advection return 1.0 return self._radius * np.pi / 180.0 diff --git a/src/parcels/_core/xgrid.py b/src/parcels/_core/xgrid.py index 841e839af..ba2c4a8db 100644 --- a/src/parcels/_core/xgrid.py +++ b/src/parcels/_core/xgrid.py @@ -12,7 +12,7 @@ import parcels._typing as ptyping from parcels._core.basegrid import BaseGrid from parcels._core.index_search import _search_1d_array, _search_indices_curvilinear_2d -from parcels._core.mesh import SphericalMesh, EARTH_RADIUS +from parcels._core.mesh import EARTH_RADIUS, SphericalMesh from parcels._sgrid.accessor import _get_dim_to_axis_mapping from parcels._sgrid.core import SGRID_PADDING_TO_XGCM_POSITION @@ -258,7 +258,7 @@ def time(self): @property def deg2m(self) -> float: """Metres per degree of arc for this grid's mesh.""" - if self._radius is None: # flat mesh; None causes crash in advection + if self._radius is None: # flat mesh; None causes crash in advection return 1.0 return self._radius * np.pi / 180.0 diff --git a/tests/test_mesh.py b/tests/test_mesh.py index 2f42ad33f..61be7516d 100644 --- a/tests/test_mesh.py +++ b/tests/test_mesh.py @@ -8,12 +8,14 @@ EARTH_DEG2M = EARTH_RADIUS * np.pi / 180 + def test_spherical_mesh_deg2m(): assert SphericalMesh().radius == EARTH_RADIUS assert SphericalMesh().deg2m == EARTH_DEG2M r = 3389500.0 # Mars radius assert SphericalMesh(radius=r).deg2m == pytest.approx(r * np.pi / 180) + @pytest.mark.parametrize( "mesh, exp_radius, exp_deg2m", [ @@ -28,14 +30,16 @@ def test_xgrid_radius_and_deg2m(mesh, exp_radius, exp_deg2m): assert grid._radius == exp_radius assert grid.deg2m == pytest.approx(exp_deg2m) -@pytest.mark.parametrize("mesh, deg2m", - [ - (SphericalMesh(), EARTH_DEG2M), # No radius entered - (SphericalMesh(radius=3389500.0), 3389500.0 * np.pi / 180), # Mars - (SphericalMesh(radius=6051800.0), 6051800.0 * np.pi / 180), # Venus - (SphericalMesh(radius=EARTH_RADIUS), EARTH_DEG2M), # Explicit Earth - ], - ) + +@pytest.mark.parametrize( + "mesh, deg2m", + [ + (SphericalMesh(), EARTH_DEG2M), # No radius entered + (SphericalMesh(radius=3389500.0), 3389500.0 * np.pi / 180), # Mars + (SphericalMesh(radius=6051800.0), 6051800.0 * np.pi / 180), # Venus + (SphericalMesh(radius=EARTH_RADIUS), EARTH_DEG2M), # Explicit Earth + ], +) def test_advection_uses_custom_radius(mesh, deg2m, npart=10): ds = simple_UV_dataset() ds["U"].data[:] = 1.0 @@ -51,26 +55,29 @@ def test_advection_uses_custom_radius(mesh, deg2m, npart=10): np.testing.assert_allclose(pset.x - startlon, expected_dlon, atol=1e-5) np.testing.assert_allclose(pset.y, startlat, atol=1e-5) + def test_advection_flat_mesh(npart=10): ds = simple_UV_dataset(mesh="flat") ds["U"].data[:] = 1.0 fieldset = FieldSet.from_sgrid_conventions(ds, mesh="flat") - + runtime = 7200 startlat = np.linspace(0, 80, npart) startlon = 20.0 + np.zeros(npart) pset = ParticleSet(fieldset, x=startlon, y=startlat) pset.execute(AdvectionRK4, runtime=runtime, dt=np.timedelta64(15, "m")) - assert fieldset.U.grid.deg2m == 1.0 # flat mesh deg2m + assert fieldset.U.grid.deg2m == 1.0 # flat mesh deg2m np.testing.assert_allclose(pset.x - startlon, runtime, atol=1e-5) np.testing.assert_allclose(pset.y, startlat, atol=1e-5) + @pytest.mark.parametrize("bad_radius", ["6371000", [6371000], (1, 2), {}]) def test_spherical_mesh_rejects_non_numeric_radius(bad_radius): with pytest.raises(TypeError): SphericalMesh(radius=bad_radius) + @pytest.mark.parametrize("bad_radius", [0, -1.0, -6371000]) def test_spherical_mesh_rejects_nonpos_radius(bad_radius): with pytest.raises(ValueError): From 93a364b802de05eb199498ac97a2dbb2ab9b1977 Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Wed, 22 Jul 2026 16:22:12 +0200 Subject: [PATCH 19/24] Refactorings after review - Adds a FlatMesh object - Remove assert_valid_mesh --- src/parcels/_core/basegrid.py | 4 +- src/parcels/_core/fieldset.py | 4 +- src/parcels/_core/index_search.py | 4 +- src/parcels/_core/kernel.py | 2 +- src/parcels/_core/mesh.py | 48 ++++++++++++++++++- src/parcels/_core/model.py | 11 +++-- src/parcels/_core/particlefile.py | 5 +- src/parcels/_core/spatialhash.py | 6 +-- src/parcels/_core/utils/interpolation.py | 14 +++--- src/parcels/_core/uxgrid.py | 20 +++----- src/parcels/_core/xgrid.py | 16 ++----- src/parcels/_typing.py | 12 +---- src/parcels/interpolators/_uxinterpolators.py | 2 +- src/parcels/interpolators/_xinterpolators.py | 8 ++-- src/parcels/kernels/_advectiondiffusion.py | 18 +++---- tests/test_mesh.py | 4 +- tests/test_typing.py | 20 -------- tests/test_uxgrid.py | 3 +- tests/test_xgrid.py | 2 +- 19 files changed, 106 insertions(+), 97 deletions(-) diff --git a/src/parcels/_core/basegrid.py b/src/parcels/_core/basegrid.py index ce4f88e3f..7b7c7e9c4 100644 --- a/src/parcels/_core/basegrid.py +++ b/src/parcels/_core/basegrid.py @@ -8,7 +8,7 @@ import numpy as np -import parcels._typing as ptyping +from parcels._core.mesh import FlatMesh, SphericalMesh from parcels._core.spatialhash import SpatialHash if TYPE_CHECKING: @@ -26,7 +26,7 @@ class BaseGrid(ABC): """Base class for parcels.XGrid and parcels.UxGrid defining common methods and properties""" _spatialhash: SpatialHash | None - _mesh: ptyping.Mesh + _mesh: FlatMesh | SphericalMesh @abstractmethod def search(self, z: float, y: float, x: float, ei=None) -> dict[str, tuple[int, float | np.ndarray]]: diff --git a/src/parcels/_core/fieldset.py b/src/parcels/_core/fieldset.py index 9ed2f9460..e141743d4 100644 --- a/src/parcels/_core/fieldset.py +++ b/src/parcels/_core/fieldset.py @@ -152,7 +152,7 @@ def add_field(self, field: Field, name: str | None = None): self.fields[name] = field - def add_constant_field(self, name: str, value, mesh: ptyping.Mesh = "spherical"): + def add_constant_field(self, name: str, value, mesh: ptyping.TMesh = "spherical"): """Wrapper function to add a Field that is constant in space, useful e.g. when using constant horizontal diffusivity @@ -249,7 +249,7 @@ def from_ugrid_conventions( def from_sgrid_conventions( cls, ds: xr.Dataset, - mesh: ptyping.Mesh | None = None, + mesh: ptyping.TMesh | None = None, vector_fields: ptyping.VectorFields | NotSetType = NOTSET, ): # TODO: Update mesh to be discovered from the dataset metadata """Create a FieldSet from a dataset using SGRID convention metadata. diff --git a/src/parcels/_core/index_search.py b/src/parcels/_core/index_search.py index 5e84eee2e..cafd58b4c 100644 --- a/src/parcels/_core/index_search.py +++ b/src/parcels/_core/index_search.py @@ -109,7 +109,7 @@ def curvilinear_point_in_cell(grid, y: np.ndarray, x: np.ndarray, yi: np.ndarray dtype=float, ) - if grid._mesh == "spherical": + if grid._mesh.is_spherical(): xsi, eta = _bilinear_inverse_tangent_plane(clon, clat, x, y) is_in_cell = np.where((xsi >= 0) & (xsi <= 1) & (eta >= 0) & (eta <= 1), 1, 0) else: @@ -318,7 +318,7 @@ def uxgrid_point_in_cell(grid, y: np.ndarray, x: np.ndarray, yi: np.ndarray, xi: coords : np.ndarray Barycentric coordinates of the points within their respective cells. """ - if grid._mesh == "spherical": + if grid._mesh.is_spherical(): lon_rad = np.deg2rad(x) lat_rad = np.deg2rad(y) x_cart, y_cart, z_cart = _latlon_rad_to_xyz(lat_rad, lon_rad) diff --git a/src/parcels/_core/kernel.py b/src/parcels/_core/kernel.py index 6a526461f..f68c8bf09 100644 --- a/src/parcels/_core/kernel.py +++ b/src/parcels/_core/kernel.py @@ -141,7 +141,7 @@ def check_fieldsets_in_kernels(self, kernel): # TODO v4: this can go into anoth stacklevel=2, ) self.fieldset.add_context("RK45_tol", 10) - if self.fieldset.U.grid._mesh == "spherical": + if self.fieldset.U.grid._mesh.is_spherical(): self.fieldset.RK45_tol /= ( self.fieldset.U.grid.deg2m ) # TODO does not account for zonal variation in meter -> degree conversion diff --git a/src/parcels/_core/mesh.py b/src/parcels/_core/mesh.py index 10311a3d0..a5ba61cbf 100644 --- a/src/parcels/_core/mesh.py +++ b/src/parcels/_core/mesh.py @@ -1,9 +1,19 @@ +from abc import ABC, abstractmethod +from typing import Literal + import numpy as np EARTH_RADIUS = 6366707.019493707 -class SphericalMesh: +class BaseMesh(ABC): + radius: float | None + + @abstractmethod + def is_spherical(self) -> bool: ... + + +class SphericalMesh(BaseMesh): """Spherical mesh object with configurable planetary radius. Pass to FieldSet object as ``mesh=SphericalMesh(radius=...)``. @@ -22,5 +32,41 @@ def deg2m(self) -> float: """Meters per degree of arc.""" return self.radius * np.pi / 180.0 + def is_spherical(self): + return True + def __repr__(self) -> str: return f"SphericalMesh(radius={self.radius})" + + +class FlatMesh(BaseMesh): + """Flat mesh object.""" + + def __init__(self): + self.radius = None + return + + def __repr__(self) -> str: + return "FlatMesh()" + + def is_spherical(self): + return False + + +TMesh = SphericalMesh | Literal["spherical", "flat"] # corresponds with `mesh` + + +def get_mesh(mesh: TMesh): + if isinstance(mesh, SphericalMesh): + return mesh + if mesh == "flat": + return FlatMesh() + if mesh == "spherical": + return SphericalMesh(EARTH_RADIUS) + raise ValueError(f"mesh must be 'flat', 'spherical', or a SphericalMesh object. Got {mesh=!r}") + + +def is_spherical(mesh: FlatMesh | SphericalMesh): + if isinstance(mesh, SphericalMesh): + return True + return False diff --git a/src/parcels/_core/model.py b/src/parcels/_core/model.py index c44026370..514a2ba79 100644 --- a/src/parcels/_core/model.py +++ b/src/parcels/_core/model.py @@ -21,7 +21,6 @@ ) from parcels._logger import logger from parcels._python import NOTSET, NotSetType -from parcels._typing import Mesh from parcels.convert import _ds_rename_using_standard_names from parcels.interpolators import ( CGrid_Velocity, @@ -83,7 +82,7 @@ def preprocess_sgrid_model_data(ds: xr.Dataset) -> xr.Dataset: class StructuredModelData(ModelData): - def __init__(self, data: xr.Dataset, mesh: Mesh, vector_field_components: ptyping.VectorFields): + def __init__(self, data: xr.Dataset, mesh: ptyping.TMesh, vector_field_components: ptyping.VectorFields): if not isinstance(data, xr.Dataset): raise ValueError(f"Expected `data` to be an xarray.Dataset . Got {type(data)}") @@ -132,7 +131,7 @@ def construct_fields(self) -> list[Field | VectorField]: @classmethod def from_sgrid_conventions( - cls, ds: xr.Dataset, mesh: Mesh | None, vector_fields: ptyping.VectorFields | NotSetType + cls, ds: xr.Dataset, mesh: ptyping.TMesh | None, vector_fields: ptyping.VectorFields | NotSetType ) -> Self: ds = ds.copy() if mesh is None: @@ -279,7 +278,9 @@ def scalar_field_names(self) -> list[str]: return list(self.data.data_vars) @classmethod - def from_ugrid_conventions(cls, ds: ux.UxDataset, mesh: Mesh, vector_fields: ptyping.VectorFields | NotSetType): + def from_ugrid_conventions( + cls, ds: ux.UxDataset, mesh: ptyping.TMesh, vector_fields: ptyping.VectorFields | NotSetType + ): ds_dims = list(ds.dims) if not all(dim in ds_dims for dim in ["time", "zf", "zc"]): raise ValueError( @@ -303,7 +304,7 @@ def from_ugrid_conventions(cls, ds: ux.UxDataset, mesh: Mesh, vector_fields: pty # TODO: Refactor later into something like `parcels._metadata.discover(dataset)` helper that can be used to discover important metadata like this. I think this whole metadata handling should be refactored into its own module. -def _get_mesh_type_from_sgrid_dataset(ds_sgrid: xr.Dataset) -> Mesh: +def _get_mesh_type_from_sgrid_dataset(ds_sgrid: xr.Dataset) -> ptyping.TMesh: """Small helper to inspect SGRID metadata and dataset metadata to determine mesh type.""" sgrid_metadata = ds_sgrid.sgrid.metadata diff --git a/src/parcels/_core/particlefile.py b/src/parcels/_core/particlefile.py index 64eb45e66..50bdf309c 100644 --- a/src/parcels/_core/particlefile.py +++ b/src/parcels/_core/particlefile.py @@ -14,6 +14,7 @@ import xarray as xr import parcels +from parcels._core.mesh import BaseMesh from parcels._core.particle import ParticleClass from parcels._core.particlesetview import ParticleSetView from parcels._core.utils.time import timedelta_to_float @@ -119,14 +120,14 @@ def __init__( def __repr__(self) -> str: return particlefile_repr(self) - def set_metadata(self, parcels_grid_mesh: Literal["spherical", "flat"]): + def set_metadata(self, parcels_grid_mesh: BaseMesh): self.metadata.update( { "feature_type": "trajectory", "Conventions": "CF-1.6/CF-1.7", "ncei_template_version": "NCEI_NetCDF_Trajectory_Template_v2.0", "parcels_version": parcels.__version__, - "parcels_grid_mesh": parcels_grid_mesh, + "parcels_grid_mesh": repr(parcels_grid_mesh), } ) diff --git a/src/parcels/_core/spatialhash.py b/src/parcels/_core/spatialhash.py index c4f4d0a64..5dc886200 100644 --- a/src/parcels/_core/spatialhash.py +++ b/src/parcels/_core/spatialhash.py @@ -45,7 +45,7 @@ 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": + if self._source_grid._mesh.is_spherical(): # Boundaries of the hash grid are the unit cube self._xmin = -1.0 self._ymin = -1.0 @@ -148,7 +148,7 @@ 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": + if self._source_grid._mesh.is_spherical(): # Boundaries of the hash grid are the unit cube self._xmin = -1.0 self._ymin = -1.0 @@ -346,7 +346,7 @@ def query(self, y, x): y = np.asarray(y) x = np.asarray(x) - if self._source_grid._mesh == "spherical": + if self._source_grid._mesh.is_spherical(): # Convert coords to Cartesian coordinates (x, y, z) lat = np.deg2rad(y) lon = np.deg2rad(x) diff --git a/src/parcels/_core/utils/interpolation.py b/src/parcels/_core/utils/interpolation.py index 214ba1933..5900c2380 100644 --- a/src/parcels/_core/utils/interpolation.py +++ b/src/parcels/_core/utils/interpolation.py @@ -3,7 +3,7 @@ import numpy as np -from parcels._typing import Mesh +from parcels._core.mesh import BaseMesh __all__ = [] @@ -61,12 +61,12 @@ def dphidxsi3D_lin(zeta: float, eta: float, xsi: float) -> tuple[list[float], li def dxdxsi3D_lin( - hexa_z: list[float], hexa_y: list[float], hexa_x: list[float], zeta: float, eta: float, xsi: float, mesh: Mesh, + hexa_z: list[float], hexa_y: list[float], hexa_x: list[float], zeta: float, eta: float, xsi: float, mesh: BaseMesh, deg2m: float = 1852 * 60.0 ) -> tuple[float, float, float, float, float, float, float, float, float]: dphidxsi, dphideta, dphidzet = dphidxsi3D_lin(zeta, eta, xsi) - if mesh == 'spherical': + if mesh.is_spherical(): rad = np.pi / 180. lat = (1-xsi) * (1-eta) * hexa_y[0] + \ xsi * (1-eta) * hexa_y[1] + \ @@ -92,7 +92,7 @@ def dxdxsi3D_lin( def jacobian3D_lin( - hexa_z: list[float], hexa_y: list[float], hexa_x: list[float], zeta: float, eta: float, xsi: float, mesh: Mesh, + hexa_z: list[float], hexa_y: list[float], hexa_x: list[float], zeta: float, eta: float, xsi: float, mesh: BaseMesh, deg2m: float = 1852 * 60.0 ) -> float: dxdxsi, dxdeta, dxdzet, dydxsi, dydeta, dydzet, dzdxsi, dzdeta, dzdzet = dxdxsi3D_lin(hexa_z, hexa_y, hexa_x, zeta, eta, xsi, mesh, deg2m) @@ -113,7 +113,7 @@ def jacobian3D_lin_face( eta: float, xsi: float, orientation: Literal["zonal", "meridional", "vertical"], - mesh: Mesh, + mesh: BaseMesh, ) -> float: dxdxsi, dxdeta, dxdzet, dydxsi, dydeta, dydzet, dzdxsi, dzdeta, dzdzet = dxdxsi3D_lin(hexa_z, hexa_y, hexa_x, zeta, eta, xsi, mesh) @@ -176,9 +176,9 @@ def interpolate(phi: Callable[[float], list[float]], f: list[float], xsi: float) def _geodetic_distance( - lat1: float, lat2: float, lon1: float, lon2: float, mesh: Mesh, lat: float, deg2m: float = 1852 * 60.0 + lat1: float, lat2: float, lon1: float, lon2: float, mesh: BaseMesh, lat: float, deg2m: float = 1852 * 60.0 ) -> float: - if mesh == "spherical": + if mesh.is_spherical(): rad = np.pi / 180.0 return np.sqrt(((lon2 - lon1) * deg2m * np.cos(rad * lat)) ** 2 + ((lat2 - lat1) * deg2m) ** 2) else: diff --git a/src/parcels/_core/uxgrid.py b/src/parcels/_core/uxgrid.py index 05adca6d7..7f9360ab9 100644 --- a/src/parcels/_core/uxgrid.py +++ b/src/parcels/_core/uxgrid.py @@ -7,8 +7,7 @@ from parcels._core.basegrid import BaseGrid from parcels._core.index_search import GRID_SEARCH_ERROR, _search_1d_array, uxgrid_point_in_cell -from parcels._core.mesh import EARTH_RADIUS, SphericalMesh -from parcels._typing import assert_valid_mesh +from parcels._core.mesh import SphericalMesh, get_mesh _UXGRID_AXES = Literal["Z", "FACE"] @@ -19,7 +18,9 @@ class UxGrid(BaseGrid): for interpolation on unstructured grids. """ - def __init__(self, grid: ux.grid.Grid, z: ux.UxDataArray, mesh) -> None: + def __init__( + self, grid: ux.grid.Grid, z: ux.UxDataArray, mesh: Literal["flat", "spherical"] | SphericalMesh + ) -> None: """ Initializes the UxGrid with a uxarray grid and vertical coordinate array. @@ -42,16 +43,9 @@ def __init__(self, grid: ux.grid.Grid, z: ux.UxDataArray, mesh) -> None: if z.ndim != 1: raise ValueError("z must be a 1D array of vertical coordinates") self.z = z - if isinstance(mesh, SphericalMesh): - self._mesh = "spherical" - self._radius = mesh.radius - else: - self._mesh = mesh - self._radius = EARTH_RADIUS if mesh == "spherical" else None + self._mesh = get_mesh(mesh) self._spatialhash = None - assert_valid_mesh(mesh) - @property def depth(self): """ @@ -82,9 +76,9 @@ def get_axis_dim(self, axis: _UXGRID_AXES) -> int: @property def deg2m(self) -> float: """Metres per arcdegree for this grid's mesh.""" - if self._radius is None: # flat mesh; None causes crash in advection + if not self._mesh.is_spherical(): return 1.0 - return self._radius * np.pi / 180.0 + return self._mesh.deg2m def search(self, z, y, x, ei=None, tol=1e-6): """ diff --git a/src/parcels/_core/xgrid.py b/src/parcels/_core/xgrid.py index ba2c4a8db..20bd695f1 100644 --- a/src/parcels/_core/xgrid.py +++ b/src/parcels/_core/xgrid.py @@ -12,7 +12,7 @@ import parcels._typing as ptyping from parcels._core.basegrid import BaseGrid from parcels._core.index_search import _search_1d_array, _search_indices_curvilinear_2d -from parcels._core.mesh import EARTH_RADIUS, SphericalMesh +from parcels._core.mesh import SphericalMesh, get_mesh from parcels._sgrid.accessor import _get_dim_to_axis_mapping from parcels._sgrid.core import SGRID_PADDING_TO_XGCM_POSITION @@ -165,17 +165,12 @@ class XGrid(BaseGrid): """ - def __init__(self, model_data: xr.Dataset, mesh): + def __init__(self, model_data: xr.Dataset, mesh: Literal["flat", "spherical"] | SphericalMesh): self.sgrid_metadata = model_data.sgrid.metadata self._ds = model_data grid = XgcmLikeGrid(self.sgrid_metadata, model_data) self.xgcm_grid = grid - if isinstance(mesh, SphericalMesh): - self._mesh = "spherical" - self._radius = mesh.radius - else: - self._mesh = mesh - self._radius = EARTH_RADIUS if mesh == "spherical" else None + self._mesh = get_mesh(mesh) self._spatialhash = None ds = model_data @@ -191,7 +186,6 @@ def __init__(self, model_data: xr.Dataset, mesh): if "Z" in grid.axes: assert_valid_depth(ds["depth"]) - ptyping.assert_valid_mesh(mesh) self._ds = ds # def __repr__(self): @@ -258,9 +252,9 @@ def time(self): @property def deg2m(self) -> float: """Metres per degree of arc for this grid's mesh.""" - if self._radius is None: # flat mesh; None causes crash in advection + if not self._mesh.is_spherical(): return 1.0 - return self._radius * np.pi / 180.0 + return self._mesh.deg2m @cached_property def xdim(self) -> int: diff --git a/src/parcels/_typing.py b/src/parcels/_typing.py index 61fc28160..87e23a53c 100644 --- a/src/parcels/_typing.py +++ b/src/parcels/_typing.py @@ -9,12 +9,12 @@ import os from collections.abc import Callable, Mapping from datetime import datetime -from typing import TYPE_CHECKING, Any, Literal, get_args +from typing import TYPE_CHECKING, Literal, get_args import numpy as np from cftime import datetime as cftime_datetime -from parcels._core.mesh import SphericalMesh +from parcels._core.mesh import TMesh # noqa: F401 if TYPE_CHECKING: import xgcm @@ -35,7 +35,6 @@ InterpMethodOption | dict[str, InterpMethodOption] ) # corresponds with `interp_method` (which can also be dict mapping field names to method) PathLike = str | os.PathLike -Mesh = Literal["spherical", "flat"] # corresponds with `mesh` VectorType = Literal["3D", "3DSigma", "2D"] | None # corresponds with `vector_type` GridIndexingType = Literal["pop", "mom5", "mitgcm", "nemo", "croco"] # corresponds with `gridindexingtype` NetcdfEngine = Literal["netcdf4", "xarray", "scipy"] @@ -71,10 +70,3 @@ def _validate_against_pure_literal(value, typing_literal): if value not in get_args(typing_literal): msg = f"Invalid value {value!r}. Valid options are {get_args(typing_literal)!r}" raise ValueError(msg) - - -# Assertion functions to clean user input -def assert_valid_mesh(value: Any): - if isinstance(value, SphericalMesh): - return - _validate_against_pure_literal(value, Mesh) diff --git a/src/parcels/interpolators/_uxinterpolators.py b/src/parcels/interpolators/_uxinterpolators.py index 0b04b1de6..17249bfdc 100644 --- a/src/parcels/interpolators/_uxinterpolators.py +++ b/src/parcels/interpolators/_uxinterpolators.py @@ -170,7 +170,7 @@ def interp( ): u = vectorfield.U.interp_method.interp(particle_positions, grid_positions, vectorfield.U) v = vectorfield.V.interp_method.interp(particle_positions, grid_positions, vectorfield.V) - if vectorfield.grid._mesh == "spherical": + if vectorfield.grid._mesh.is_spherical(): u /= vectorfield.grid.deg2m * np.cos(np.deg2rad(particle_positions["y"])) v /= vectorfield.grid.deg2m diff --git a/src/parcels/interpolators/_xinterpolators.py b/src/parcels/interpolators/_xinterpolators.py index d851e42c5..38b10d2ea 100644 --- a/src/parcels/interpolators/_xinterpolators.py +++ b/src/parcels/interpolators/_xinterpolators.py @@ -148,7 +148,7 @@ def interp( _xlinear = XLinear() u = _xlinear.interp(particle_positions, grid_positions, vectorfield.U) v = _xlinear.interp(particle_positions, grid_positions, vectorfield.V) - if vectorfield.grid._mesh == "spherical": + if vectorfield.grid._mesh.is_spherical(): u /= vectorfield.grid.deg2m * np.cos(np.deg2rad(particle_positions["y"])) v /= vectorfield.grid.deg2m @@ -196,7 +196,7 @@ def interp( px = np.array([grid.lon[yi, xi], grid.lon[yi, xi + 1], grid.lon[yi + 1, xi + 1], grid.lon[yi + 1, xi]]) py = np.array([grid.lat[yi, xi], grid.lat[yi, xi + 1], grid.lat[yi + 1, xi + 1], grid.lat[yi + 1, xi]]) - if grid._mesh == "spherical": + if grid._mesh.is_spherical(): px = ((px + 180.0) % 360.0) - 180.0 px[1:] = np.where(px[1:] - px[0] > 180, px[1:] - 360, px[1:]) px[1:] = np.where(-px[1:] + px[0] > 180, px[1:] + 360, px[1:]) @@ -282,7 +282,7 @@ def _compute_corner_data(data, selection_dict) -> np.ndarray: V1 = corner_data[1, :] * c3 Vvel = (1 - eta) * V0 + eta * V1 - if grid._mesh == "spherical": + if grid._mesh.is_spherical(): jac = i_u._compute_jacobian_determinant(py, px, eta, xsi) * grid.deg2m else: jac = i_u._compute_jacobian_determinant(py, px, eta, xsi) @@ -303,7 +303,7 @@ def _compute_corner_data(data, selection_dict) -> np.ndarray: u = u.compute() v = v.compute() - if grid._mesh == "spherical": + if grid._mesh.is_spherical(): conversion = grid.deg2m * np.cos(np.deg2rad(particle_positions["y"])) u /= conversion v /= conversion diff --git a/src/parcels/kernels/_advectiondiffusion.py b/src/parcels/kernels/_advectiondiffusion.py index 7b7a96276..b5d327de9 100644 --- a/src/parcels/kernels/_advectiondiffusion.py +++ b/src/parcels/kernels/_advectiondiffusion.py @@ -39,26 +39,26 @@ def AdvectionDiffusionM1(particles, fieldset): # pragma: no cover Kxp1 = fieldset.Kh_zonal[particles.t, particles.z, particles.y, particles.x + fieldset.dres, particles] Kxm1 = fieldset.Kh_zonal[particles.t, particles.z, particles.y, particles.x - fieldset.dres, particles] - if fieldset.Kh_zonal.grid._mesh == "spherical": + if fieldset.Kh_zonal.grid._mesh.is_spherical(): Kxp1 = meters_to_degrees_zonal(Kxp1, particles.y, fieldset.Kh_zonal.grid.deg2m) Kxm1 = meters_to_degrees_zonal(Kxm1, particles.y, fieldset.Kh_zonal.grid.deg2m) dKdx = (Kxp1 - Kxm1) / (2 * fieldset.dres) u, v = fieldset.UV[particles.t, particles.z, particles.y, particles.x, particles] kh_zonal = fieldset.Kh_zonal[particles.t, particles.z, particles.y, particles.x, particles] - if fieldset.Kh_zonal.grid._mesh == "spherical": + if fieldset.Kh_zonal.grid._mesh.is_spherical(): kh_zonal = meters_to_degrees_zonal(kh_zonal, particles.y, fieldset.Kh_zonal.grid.deg2m) bx = np.sqrt(2 * kh_zonal) Kyp1 = fieldset.Kh_meridional[particles.t, particles.z, particles.y + fieldset.dres, particles.x, particles] Kym1 = fieldset.Kh_meridional[particles.t, particles.z, particles.y - fieldset.dres, particles.x, particles] - if fieldset.Kh_meridional.grid._mesh == "spherical": + if fieldset.Kh_meridional.grid._mesh.is_spherical(): Kyp1 = meters_to_degrees_meridional(Kyp1, fieldset.Kh_meridional.grid.deg2m) Kym1 = meters_to_degrees_meridional(Kym1, fieldset.Kh_meridional.grid.deg2m) dKdy = (Kyp1 - Kym1) / (2 * fieldset.dres) kh_meridional = fieldset.Kh_meridional[particles.t, particles.z, particles.y, particles.x, particles] - if fieldset.Kh_meridional.grid._mesh == "spherical": + if fieldset.Kh_meridional.grid._mesh.is_spherical(): kh_meridional = meters_to_degrees_meridional(kh_meridional, fieldset.Kh_meridional.grid.deg2m) by = np.sqrt(2 * kh_meridional) @@ -88,27 +88,27 @@ def AdvectionDiffusionEM(particles, fieldset): # pragma: no cover Kxp1 = fieldset.Kh_zonal[particles.t, particles.z, particles.y, particles.x + fieldset.dres, particles] Kxm1 = fieldset.Kh_zonal[particles.t, particles.z, particles.y, particles.x - fieldset.dres, particles] - if fieldset.Kh_zonal.grid._mesh == "spherical": + if fieldset.Kh_zonal.grid._mesh.is_spherical(): Kxp1 = meters_to_degrees_zonal(Kxp1, particles.y, fieldset.Kh_zonal.grid.deg2m) Kxm1 = meters_to_degrees_zonal(Kxm1, particles.y, fieldset.Kh_zonal.grid.deg2m) dKdx = (Kxp1 - Kxm1) / (2 * fieldset.dres) ax = u + dKdx kh_zonal = fieldset.Kh_zonal[particles.t, particles.z, particles.y, particles.x, particles] - if fieldset.Kh_zonal.grid._mesh == "spherical": + if fieldset.Kh_zonal.grid._mesh.is_spherical(): kh_zonal = meters_to_degrees_zonal(kh_zonal, particles.y, fieldset.Kh_zonal.grid.deg2m) bx = np.sqrt(2 * kh_zonal) Kyp1 = fieldset.Kh_meridional[particles.t, particles.z, particles.y + fieldset.dres, particles.x, particles] Kym1 = fieldset.Kh_meridional[particles.t, particles.z, particles.y - fieldset.dres, particles.x, particles] - if fieldset.Kh_meridional.grid._mesh == "spherical": + if fieldset.Kh_meridional.grid._mesh.is_spherical(): Kyp1 = meters_to_degrees_meridional(Kyp1, fieldset.Kh_meridional.grid.deg2m) Kym1 = meters_to_degrees_meridional(Kym1, fieldset.Kh_meridional.grid.deg2m) dKdy = (Kyp1 - Kym1) / (2 * fieldset.dres) ay = v + dKdy kh_meridional = fieldset.Kh_meridional[particles.t, particles.z, particles.y, particles.x, particles] - if fieldset.Kh_meridional.grid._mesh == "spherical": + if fieldset.Kh_meridional.grid._mesh.is_spherical(): kh_meridional = meters_to_degrees_meridional(kh_meridional, fieldset.Kh_meridional.grid.deg2m) by = np.sqrt(2 * kh_meridional) @@ -142,7 +142,7 @@ def DiffusionUniformKh(particles, fieldset): # pragma: no cover kh_zonal = fieldset.Kh_zonal[particles] kh_meridional = fieldset.Kh_meridional[particles] - if fieldset.Kh_zonal.grid._mesh == "spherical": + if fieldset.Kh_zonal.grid._mesh.is_spherical(): kh_zonal = meters_to_degrees_zonal(kh_zonal, particles.y, fieldset.Kh_zonal.grid.deg2m) kh_meridional = meters_to_degrees_meridional(kh_meridional, fieldset.Kh_meridional.grid.deg2m) diff --git a/tests/test_mesh.py b/tests/test_mesh.py index 61be7516d..9a0a4b709 100644 --- a/tests/test_mesh.py +++ b/tests/test_mesh.py @@ -26,8 +26,8 @@ def test_spherical_mesh_deg2m(): ) def test_xgrid_radius_and_deg2m(mesh, exp_radius, exp_deg2m): grid = FieldSet.from_sgrid_conventions(simple_UV_dataset(), mesh=mesh).U.grid - assert grid._mesh == "spherical" - assert grid._radius == exp_radius + assert isinstance(grid._mesh, SphericalMesh) + assert grid._mesh.radius == exp_radius assert grid.deg2m == pytest.approx(exp_deg2m) diff --git a/tests/test_typing.py b/tests/test_typing.py index 4f1d3e0c5..e69de29bb 100644 --- a/tests/test_typing.py +++ b/tests/test_typing.py @@ -1,20 +0,0 @@ -import numpy as np -import pytest -import xarray as xr - -from parcels._typing import ( - assert_valid_mesh, -) - - -def test_invalid_assert_valid_mesh(): - with pytest.raises(ValueError, match="Invalid value"): - assert_valid_mesh("invalid option") - - ds = xr.Dataset({"A": (("a", "b"), np.arange(20).reshape(4, 5))}) - with pytest.raises(ValueError, match="Invalid input type"): - assert_valid_mesh(ds) - - -def test_assert_valid_mesh(): - assert_valid_mesh("spherical") diff --git a/tests/test_uxgrid.py b/tests/test_uxgrid.py index 5e3038cf3..63df6a3ca 100644 --- a/tests/test_uxgrid.py +++ b/tests/test_uxgrid.py @@ -22,7 +22,8 @@ def test_uxgrid_axes(uxds): @pytest.mark.parametrize("mesh", ["flat", "spherical"]) def test_uxgrid_mesh(uxds, mesh): grid = UxGrid(uxds.uxgrid, z=uxds.coords["zf"], mesh=mesh) - assert grid._mesh == mesh + + assert mesh in grid._mesh.__class__.__name__.lower() @pytest.mark.parametrize("uxds", [uxdatasets["stommel_gyre_delaunay"]]) diff --git a/tests/test_xgrid.py b/tests/test_xgrid.py index 2ba04e115..d957ab0b9 100644 --- a/tests/test_xgrid.py +++ b/tests/test_xgrid.py @@ -72,7 +72,7 @@ def test_xgrid_axes(fieldset): @pytest.mark.parametrize("mesh", ["flat", "spherical"]) def test_uxgrid_mesh(ds, mesh): grid = FieldSet.from_sgrid_conventions(ds, mesh=mesh).data_g.grid - assert grid._mesh == mesh + assert mesh in grid._mesh.__class__.__name__.lower() @pytest.mark.skip( From d466ccd81c986d1e089ff4bddf3bfd14aa4f98d8 Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Wed, 22 Jul 2026 17:42:04 +0200 Subject: [PATCH 20/24] Fix typing --- src/parcels/_core/mesh.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/parcels/_core/mesh.py b/src/parcels/_core/mesh.py index a5ba61cbf..b60b0d6c8 100644 --- a/src/parcels/_core/mesh.py +++ b/src/parcels/_core/mesh.py @@ -30,6 +30,7 @@ def __init__(self, radius: float = EARTH_RADIUS): @property def deg2m(self) -> float: """Meters per degree of arc.""" + assert self.radius is not None return self.radius * np.pi / 180.0 def is_spherical(self): From 27e6456d78db64b74d686f3233f47ffd0112d238 Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Wed, 22 Jul 2026 17:49:54 +0200 Subject: [PATCH 21/24] Fix tests --- tests/test_fieldset.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/tests/test_fieldset.py b/tests/test_fieldset.py index abdd26062..995021c18 100644 --- a/tests/test_fieldset.py +++ b/tests/test_fieldset.py @@ -454,8 +454,8 @@ def test_fieldset_describe_backends(tmp_path): | UV | VectorField | 0 | CGrid_Velocity(...) | - | | UVW | VectorField | 0 | CGrid_Velocity(...) | - | -mesh: spherical -time interval: (np.datetime64('2000-01-02T12:00:00.000000000'), np.datetime64('2000-01-12T12:00:00.000000000')) +mesh: SphericalMesh(radius=6366707.019493707) +time interval: (np.datetime64('2000-01-02T12:00:00.000000000'), np.datetime64('2000-01-27T12:00:00.000000000')) """ fieldset.describe(io) actual = io.getvalue() @@ -474,8 +474,8 @@ def test_fieldset_describe_backends(tmp_path): | UV | VectorField | 0 | CGrid_Velocity(...) | - | | UVW | VectorField | 0 | CGrid_Velocity(...) | - | -mesh: spherical -time interval: (np.datetime64('2000-01-02T12:00:00.000000000'), np.datetime64('2000-01-12T12:00:00.000000000')) +mesh: SphericalMesh(radius=6366707.019493707) +time interval: (np.datetime64('2000-01-02T12:00:00.000000000'), np.datetime64('2000-01-27T12:00:00.000000000')) """ fieldset.describe(io) actual = io.getvalue() @@ -496,8 +496,8 @@ def test_fieldset_describe_backends(tmp_path): | UV | VectorField | 0 | CGrid_Velocity(...) | - | | UVW | VectorField | 0 | CGrid_Velocity(...) | - | -mesh: spherical -time interval: (np.datetime64('2000-01-02T12:00:00.000000000'), np.datetime64('2000-01-12T12:00:00.000000000')) +mesh: SphericalMesh(radius=6366707.019493707) +time interval: (np.datetime64('2000-01-02T12:00:00.000000000'), np.datetime64('2000-01-27T12:00:00.000000000')) """ fieldset.describe(io) actual = io.getvalue() From a50c0ff0ce5479b3fd703c39be333aee51333f9b Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Thu, 23 Jul 2026 13:26:05 +0200 Subject: [PATCH 22/24] Revert "Fix tests" This reverts commit 27e6456d78db64b74d686f3233f47ffd0112d238. --- tests/test_fieldset.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/tests/test_fieldset.py b/tests/test_fieldset.py index 995021c18..abdd26062 100644 --- a/tests/test_fieldset.py +++ b/tests/test_fieldset.py @@ -454,8 +454,8 @@ def test_fieldset_describe_backends(tmp_path): | UV | VectorField | 0 | CGrid_Velocity(...) | - | | UVW | VectorField | 0 | CGrid_Velocity(...) | - | -mesh: SphericalMesh(radius=6366707.019493707) -time interval: (np.datetime64('2000-01-02T12:00:00.000000000'), np.datetime64('2000-01-27T12:00:00.000000000')) +mesh: spherical +time interval: (np.datetime64('2000-01-02T12:00:00.000000000'), np.datetime64('2000-01-12T12:00:00.000000000')) """ fieldset.describe(io) actual = io.getvalue() @@ -474,8 +474,8 @@ def test_fieldset_describe_backends(tmp_path): | UV | VectorField | 0 | CGrid_Velocity(...) | - | | UVW | VectorField | 0 | CGrid_Velocity(...) | - | -mesh: SphericalMesh(radius=6366707.019493707) -time interval: (np.datetime64('2000-01-02T12:00:00.000000000'), np.datetime64('2000-01-27T12:00:00.000000000')) +mesh: spherical +time interval: (np.datetime64('2000-01-02T12:00:00.000000000'), np.datetime64('2000-01-12T12:00:00.000000000')) """ fieldset.describe(io) actual = io.getvalue() @@ -496,8 +496,8 @@ def test_fieldset_describe_backends(tmp_path): | UV | VectorField | 0 | CGrid_Velocity(...) | - | | UVW | VectorField | 0 | CGrid_Velocity(...) | - | -mesh: SphericalMesh(radius=6366707.019493707) -time interval: (np.datetime64('2000-01-02T12:00:00.000000000'), np.datetime64('2000-01-27T12:00:00.000000000')) +mesh: spherical +time interval: (np.datetime64('2000-01-02T12:00:00.000000000'), np.datetime64('2000-01-12T12:00:00.000000000')) """ fieldset.describe(io) actual = io.getvalue() From a884bf81d643d276cf85a435aa38e337d860795d Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Thu, 23 Jul 2026 13:28:18 +0200 Subject: [PATCH 23/24] Fix tests --- tests/test_fieldset.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/test_fieldset.py b/tests/test_fieldset.py index abdd26062..015d4110e 100644 --- a/tests/test_fieldset.py +++ b/tests/test_fieldset.py @@ -454,7 +454,7 @@ def test_fieldset_describe_backends(tmp_path): | UV | VectorField | 0 | CGrid_Velocity(...) | - | | UVW | VectorField | 0 | CGrid_Velocity(...) | - | -mesh: spherical +mesh: SphericalMesh(radius=6366707.019493707) time interval: (np.datetime64('2000-01-02T12:00:00.000000000'), np.datetime64('2000-01-12T12:00:00.000000000')) """ fieldset.describe(io) @@ -474,7 +474,7 @@ def test_fieldset_describe_backends(tmp_path): | UV | VectorField | 0 | CGrid_Velocity(...) | - | | UVW | VectorField | 0 | CGrid_Velocity(...) | - | -mesh: spherical +mesh: SphericalMesh(radius=6366707.019493707) time interval: (np.datetime64('2000-01-02T12:00:00.000000000'), np.datetime64('2000-01-12T12:00:00.000000000')) """ fieldset.describe(io) @@ -496,7 +496,7 @@ def test_fieldset_describe_backends(tmp_path): | UV | VectorField | 0 | CGrid_Velocity(...) | - | | UVW | VectorField | 0 | CGrid_Velocity(...) | - | -mesh: spherical +mesh: SphericalMesh(radius=6366707.019493707) time interval: (np.datetime64('2000-01-02T12:00:00.000000000'), np.datetime64('2000-01-12T12:00:00.000000000')) """ fieldset.describe(io) From c82a69ea66467c63d56bfec00ecfb1e3c7ecc5c1 Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Thu, 23 Jul 2026 13:29:17 +0200 Subject: [PATCH 24/24] Review feedback --- src/parcels/_core/uxgrid.py | 6 +++--- src/parcels/_core/xgrid.py | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/parcels/_core/uxgrid.py b/src/parcels/_core/uxgrid.py index 7f9360ab9..1b1e778ab 100644 --- a/src/parcels/_core/uxgrid.py +++ b/src/parcels/_core/uxgrid.py @@ -76,9 +76,9 @@ def get_axis_dim(self, axis: _UXGRID_AXES) -> int: @property def deg2m(self) -> float: """Metres per arcdegree for this grid's mesh.""" - if not self._mesh.is_spherical(): - return 1.0 - return self._mesh.deg2m + if self._mesh.is_spherical(): + return self._mesh.deg2m + return 1.0 def search(self, z, y, x, ei=None, tol=1e-6): """ diff --git a/src/parcels/_core/xgrid.py b/src/parcels/_core/xgrid.py index 94724e207..de213f823 100644 --- a/src/parcels/_core/xgrid.py +++ b/src/parcels/_core/xgrid.py @@ -259,9 +259,9 @@ def time(self): @property def deg2m(self) -> float: """Metres per degree of arc for this grid's mesh.""" - if not self._mesh.is_spherical(): - return 1.0 - return self._mesh.deg2m + if self._mesh.is_spherical(): + return self._mesh.deg2m + return 1.0 @cached_property def xdim(self) -> int: