Add support for defining planetary radius

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",

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(

src/parcels/_core/mesh.py

@@ -0,0 +1,28 @@
+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):
+        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
+    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})"

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,11 @@ 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)

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)
@@ -73,6 +79,14 @@ def get_axis_dim(self, axis: _UXGRID_AXES) -> int:
         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):
         """
         Search for the grid cell (face) and vertical layer that contains the given points.

src/parcels/_core/xgrid.py

@@ -12,6 +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._sgrid.accessor import _get_dim_to_axis_mapping
 from parcels._sgrid.core import SGRID_PADDING_TO_XGCM_POSITION
 
@@ -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 = None
         self._spatialhash = None
         ds = model_data
 
@@ -249,6 +255,14 @@ def _datetimes(self):
     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:
         return self.get_axis_dim("X")

src/parcels/_typing.py

@@ -14,6 +14,8 @@
 import numpy as np
 from cftime import datetime as cftime_datetime
 
+from parcels._core.mesh import SphericalMesh
+
 if TYPE_CHECKING:
     import xgcm
 
@@ -73,4 +75,6 @@ def _validate_against_pure_literal(value, typing_literal):
 
 # Assertion functions to clean user input
 def assert_valid_mesh(value: Any):
+    if isinstance(value, SphericalMesh):
+        return
     _validate_against_pure_literal(value, Mesh)

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)

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
 

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
 

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)