From 0c4fa4549179a69a55b1ed423733478e6f71da56 Mon Sep 17 00:00:00 2001 From: Maureen Cohen Date: Tue, 7 Jul 2026 14:08:54 +0100 Subject: [PATCH 01/16] 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/16] 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/16] 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/16] 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/16] 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/16] 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/16] 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/16] 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/16] 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/16] 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/16] 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/16] [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/16] 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/16] 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/16] [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/16] 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: