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.

Author: maureenjcohen

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})"

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):
         """

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:

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,31 +28,49 @@ 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)
     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)
 
+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):