Merge branch 'v4-dev' into removing_explicit_field_chunking

Author: erikvansebille

docs/examples/tutorial_diffusion.ipynb

@@ -118,7 +118,9 @@
     "import trajan as ta\n",
     "import xarray as xr\n",
     "\n",
-    "import parcels"
+    "import parcels\n",
+    "from parcels.grid import GridType\n",
+    "from parcels.tools.converters import Geographic, GeographicPolar, UnitConverter"
    ]
   },
   {
@@ -471,8 +473,40 @@
     "x = fieldset.U.grid.lon\n",
     "y = fieldset.U.grid.lat\n",
     "\n",
+    "\n",
+    "def calc_cell_edge_sizes(grid):\n",
+    "    \"\"\"Calculate cell sizes based on numpy.gradient method.\n",
+    "\n",
+    "    Currently only works for Rectilinear Grids. Operates in place adding a `cell_edge_sizes`\n",
+    "    attribute to the grid.\n",
+    "    \"\"\"\n",
+    "    assert grid._gtype in (GridType.RectilinearZGrid, GridType.RectilinearSGrid), (\n",
+    "        f\"_cell_edge_sizes() not implemented for {grid._gtype} grids. \"\n",
+    "        \"You can provide cell_edge_sizes yourself by in, e.g., \"\n",
+    "        \"NEMO using the e1u fields etc from the mesh_mask.nc file.\"\n",
+    "    )\n",
+    "\n",
+    "    cell_edge_sizes_x = np.zeros((grid.ydim, grid.xdim), dtype=np.float32)\n",
+    "    cell_edge_sizes_y = np.zeros((grid.ydim, grid.xdim), dtype=np.float32)\n",
+    "\n",
+    "    x_conv = GeographicPolar() if grid.mesh == \"spherical\" else UnitConverter()\n",
+    "    y_conv = Geographic() if grid.mesh == \"spherical\" else UnitConverter()\n",
+    "    for y, (lat, dy) in enumerate(zip(grid.lat, np.gradient(grid.lat), strict=False)):\n",
+    "        for x, (lon, dx) in enumerate(\n",
+    "            zip(grid.lon, np.gradient(grid.lon), strict=False)\n",
+    "        ):\n",
+    "            cell_edge_sizes_x[y, x] = x_conv.to_source(dx, grid.depth[0], lat, lon)\n",
+    "            cell_edge_sizes_y[y, x] = y_conv.to_source(dy, grid.depth[0], lat, lon)\n",
+    "    return cell_edge_sizes_x, cell_edge_sizes_y\n",
+    "\n",
+    "\n",
+    "def calc_cell_areas(grid):\n",
+    "    cell_edge_sizes_x, cell_edge_sizes_y = calc_cell_edge_sizes(grid)\n",
+    "    return cell_edge_sizes_x * cell_edge_sizes_y\n",
+    "\n",
+    "\n",
     "cell_areas = parcels.Field(\n",
-    "    name=\"cell_areas\", data=fieldset.U.cell_areas(), lon=x, lat=y\n",
+    "    name=\"cell_areas\", data=calc_cell_areas(fieldset.U.grid), lon=x, lat=y\n",
     ")\n",
     "fieldset.add_field(cell_areas)\n",
     "\n",
@@ -580,7 +614,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "parcels",
+   "display_name": "parcels-dev",
    "language": "python",
    "name": "python3"
   },
@@ -594,7 +628,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.12.3"
+   "version": "3.10.15"
   }
  },
  "nbformat": 4,

parcels/_index_search.py

@@ -137,7 +137,7 @@ def search_indices_vertical_s(
 
 
 def _search_indices_rectilinear(
-    field: Field, time: float, z: float, y: float, x: float, ti=-1, particle=None, search2D=False
+    field: Field, time: float, z: float, y: float, x: float, ti: int, particle=None, search2D=False
 ):
     grid = field.grid
 
@@ -223,7 +223,7 @@ def _search_indices_rectilinear(
     return (zeta, eta, xsi, zi, yi, xi)
 
 
-def _search_indices_curvilinear(field: Field, time, z, y, x, ti=-1, particle=None, search2D=False):
+def _search_indices_curvilinear(field: Field, time, z, y, x, ti, particle=None, search2D=False):
     if particle:
         zi, yi, xi = field.unravel_index(particle.ei)
     else:

parcels/application_kernels/advection.py

@@ -185,7 +185,7 @@ def AdvectionAnalytical(particle, fieldset, time):  # pragma: no cover
         ds_t = min(ds_t, time_i[np.where(time - fieldset.U.grid.time[ti] < time_i)[0][0]])
 
     zeta, eta, xsi, zi, yi, xi = fieldset.U._search_indices(
-        -1, particle.depth, particle.lat, particle.lon, particle=particle
+        time, particle.depth, particle.lat, particle.lon, ti, particle=particle
     )
     if withW:
         if abs(xsi - 1) < tol:

parcels/field.py

@@ -46,7 +46,7 @@
     DeferredNetcdfFileBuffer,
     NetcdfFileBuffer,
 )
-from .grid import Grid, GridType, _calc_cell_areas
+from .grid import Grid, GridType
 
 if TYPE_CHECKING:
     from parcels.fieldset import FieldSet
@@ -325,10 +325,6 @@ def depth(self):
         """Depth defined on the Grid object"""
         return self.grid.depth
 
-    @property
-    def cell_edge_sizes(self):
-        return self.grid.cell_edge_sizes
-
     @property
     def interp_method(self):
         return self._interp_method
@@ -837,22 +833,15 @@ def set_depth_from_field(self, field):
         if self.grid != field.grid:
             field.grid.depth_field = field
 
-    def cell_areas(self):
-        """Method to calculate cell sizes based on cell_edge_sizes.
-
-        Only works for Rectilinear Grids
-        """
-        return _calc_cell_areas(self.grid)
-
-    def _search_indices(self, time, z, y, x, ti=-1, particle=None, search2D=False):
+    def _search_indices(self, time, z, y, x, ti, particle=None, search2D=False):
         if self.grid._gtype in [GridType.RectilinearSGrid, GridType.RectilinearZGrid]:
             return _search_indices_rectilinear(self, time, z, y, x, ti, particle=particle, search2D=search2D)
         else:
             return _search_indices_curvilinear(self, time, z, y, x, ti, particle=particle, search2D=search2D)
 
-    def _interpolator2D(self, ti, z, y, x, particle=None):
+    def _interpolator2D(self, time, z, y, x, ti, particle=None):
         """Impelement 2D interpolation with coordinate transformations as seen in Delandmeter and Van Sebille (2019), 10.5194/gmd-12-3571-2019.."""
-        (_, eta, xsi, _, yi, xi) = self._search_indices(-1, z, y, x, particle=particle)
+        (_, eta, xsi, _, yi, xi) = self._search_indices(time, z, y, x, ti, particle=particle)
         ctx = InterpolationContext2D(self.data, eta, xsi, ti, yi, xi)
 
         try:
@@ -866,7 +855,7 @@ def _interpolator2D(self, ti, z, y, x, particle=None):
                 raise RuntimeError(self.interp_method + " is not implemented for 2D grids")
         return f(ctx)
 
-    def _interpolator3D(self, ti, z, y, x, time, particle=None):
+    def _interpolator3D(self, time, z, y, x, ti, particle=None):
         """Impelement 3D interpolation with coordinate transformations as seen in Delandmeter and Van Sebille (2019), 10.5194/gmd-12-3571-2019.."""
         (zeta, eta, xsi, zi, yi, xi) = self._search_indices(time, z, y, x, ti, particle=particle)
         ctx = InterpolationContext3D(self.data, zeta, eta, xsi, ti, zi, yi, xi, self.gridindexingtype)
@@ -877,34 +866,13 @@ def _interpolator3D(self, ti, z, y, x, time, particle=None):
             raise RuntimeError(self.interp_method + " is not implemented for 3D grids")
         return f(ctx)
 
-    def temporal_interpolate_fullfield(self, ti, time):
-        """Calculate the data of a field between two snapshots using linear interpolation.
-
-        Parameters
-        ----------
-        ti :
-            Index in time array associated with time (via :func:`time_index`)
-        time :
-            Time to interpolate to
-        """
-        t0 = self.grid.time[ti]
-        if time == t0:
-            return self.data[ti, :]
-        elif ti + 1 >= len(self.grid.time):
-            raise TimeExtrapolationError(time, field=self)
-        else:
-            t1 = self.grid.time[ti + 1]
-            f0 = self.data[ti, :]
-            f1 = self.data[ti + 1, :]
-            return f0 + (f1 - f0) * ((time - t0) / (t1 - t0))
-
-    def _spatial_interpolation(self, ti, z, y, x, time, particle=None):
-        """Interpolate horizontal field values using a SciPy interpolator."""
+    def _spatial_interpolation(self, time, z, y, x, ti, particle=None):
+        """Interpolate spatial field values."""
         try:
             if self.grid.zdim == 1:
-                val = self._interpolator2D(ti, z, y, x, particle=particle)
+                val = self._interpolator2D(time, z, y, x, ti, particle=particle)
             else:
-                val = self._interpolator3D(ti, z, y, x, time, particle=particle)
+                val = self._interpolator3D(time, z, y, x, ti, particle=particle)
 
             if np.isnan(val):
                 # Detect Out-of-bounds sampling and raise exception
@@ -966,16 +934,16 @@ def eval(self, time, z, y, x, particle=None, applyConversion=True):
         if self.gridindexingtype == "croco" and self not in [self.fieldset.H, self.fieldset.Zeta]:
             z = _croco_from_z_to_sigma_scipy(self.fieldset, time, z, y, x, particle=particle)
         if ti < self.grid.tdim - 1 and time > self.grid.time[ti]:
-            f0 = self._spatial_interpolation(ti, z, y, x, time, particle=particle)
-            f1 = self._spatial_interpolation(ti + 1, z, y, x, time, particle=particle)
+            f0 = self._spatial_interpolation(time, z, y, x, ti, particle=particle)
+            f1 = self._spatial_interpolation(time, z, y, x, ti + 1, particle=particle)
             t0 = self.grid.time[ti]
             t1 = self.grid.time[ti + 1]
             value = f0 + (f1 - f0) * ((time - t0) / (t1 - t0))
         else:
             # Skip temporal interpolation if time is outside
             # of the defined time range or if we have hit an
             # exact value in the time array.
-            value = self._spatial_interpolation(ti, z, y, x, self.grid.time[ti], particle=particle)
+            value = self._spatial_interpolation(self.grid.time[ti], z, y, x, ti, particle=particle)
 
         if applyConversion:
             return self.units.to_target(value, z, y, x)

parcels/grid.py

@@ -6,7 +6,7 @@
 import numpy.typing as npt
 
 from parcels._typing import Mesh, UpdateStatus, assert_valid_mesh
-from parcels.tools.converters import Geographic, GeographicPolar, TimeConverter, UnitConverter
+from parcels.tools.converters import TimeConverter
 from parcels.tools.warnings import FieldSetWarning
 
 __all__ = [
@@ -67,7 +67,6 @@ def __init__(
         assert isinstance(self.time_origin, TimeConverter), "time_origin needs to be a TimeConverter object"
         assert_valid_mesh(mesh)
         self._mesh = mesh
-        self._cell_edge_sizes: dict[str, npt.NDArray] = {}
         self._zonal_periodic = False
         self._zonal_halo = 0
         self._meridional_halo = 0
@@ -133,10 +132,6 @@ def zonal_halo(self):
     def defer_load(self):
         return self._defer_load
 
-    @property
-    def cell_edge_sizes(self):
-        return self._cell_edge_sizes
-
     @staticmethod
     def create_grid(
         lon: npt.ArrayLike,
@@ -610,34 +605,3 @@ def __init__(
     @property
     def zdim(self):
         return self.depth.shape[-3]
-
-
-def _calc_cell_edge_sizes(grid: RectilinearGrid) -> None:
-    """Method to calculate cell sizes based on numpy.gradient method.
-
-    Currently only works for Rectilinear Grids. Operates in place adding a `cell_edge_sizes`
-    attribute to the grid.
-    """
-    if not grid.cell_edge_sizes:
-        if grid._gtype in (GridType.RectilinearZGrid, GridType.RectilinearSGrid):  # type: ignore[attr-defined]
-            grid.cell_edge_sizes["x"] = np.zeros((grid.ydim, grid.xdim), dtype=np.float32)
-            grid.cell_edge_sizes["y"] = np.zeros((grid.ydim, grid.xdim), dtype=np.float32)
-
-            x_conv = GeographicPolar() if grid.mesh == "spherical" else UnitConverter()
-            y_conv = Geographic() if grid.mesh == "spherical" else UnitConverter()
-            for y, (lat, dy) in enumerate(zip(grid.lat, np.gradient(grid.lat), strict=False)):
-                for x, (lon, dx) in enumerate(zip(grid.lon, np.gradient(grid.lon), strict=False)):
-                    grid.cell_edge_sizes["x"][y, x] = x_conv.to_source(dx, grid.depth[0], lat, lon)
-                    grid.cell_edge_sizes["y"][y, x] = y_conv.to_source(dy, grid.depth[0], lat, lon)
-        else:
-            raise ValueError(
-                f"_cell_edge_sizes() not implemented for {grid._gtype} grids. "  # type: ignore[attr-defined]
-                "You can provide Field.grid.cell_edge_sizes yourself by in, e.g., "
-                "NEMO using the e1u fields etc from the mesh_mask.nc file."
-            )
-
-
-def _calc_cell_areas(grid: RectilinearGrid) -> np.ndarray:
-    if not grid.cell_edge_sizes:
-        _calc_cell_edge_sizes(grid)
-    return grid.cell_edge_sizes["x"] * grid.cell_edge_sizes["y"]

parcels/particle.py

@@ -149,9 +149,6 @@ def __init__(self, lon, lat, pid, fieldset=None, ngrids=None, depth=0.0, time=0.
                 initial = v.initial
             setattr(self, v.name, v.dtype(initial))
 
-    def __del__(self):
-        pass  # superclass is 'object', and object itself has no destructor, hence 'pass'
-
     def __repr__(self):
         time_string = "not_yet_set" if self.time is None or np.isnan(self.time) else f"{self.time:f}"
         p_string = f"P[{self.id}](lon={self.lon:f}, lat={self.lat:f}, depth={self.depth:f}, "

parcels/particledata.py

@@ -122,8 +122,8 @@ def __init__(self, pclass, lon, lat, depth, time, lonlatdepth_dtype, pid_orig, n
         self._ncount = len(lon)
 
         for v in self.ptype.variables:
-            if v.name in ["ei", "ti"]:
-                self._data[v.name] = np.empty((len(lon), ngrid), dtype=v.dtype)
+            if v.name == "ei":
+                self._data[v.name] = np.empty((len(lon), ngrid), dtype=v.dtype)  # TODO len(lon) can be self._ncount?
             else:
                 self._data[v.name] = np.empty(self._ncount, dtype=v.dtype)
 
@@ -182,9 +182,6 @@ def __init__(self, pclass, lon, lat, depth, time, lonlatdepth_dtype, pid_orig, n
         else:
             raise ValueError("Latitude and longitude required for generating ParticleSet")
 
-    def __del__(self):
-        pass
-
     @property
     def pu_indicators(self):
         """

parcels/particleset.py

@@ -128,13 +128,11 @@ def ArrayClass_init(self, *args, **kwargs):
                 self.ngrids = type(self).ngrids.initial
                 if self.ngrids >= 0:
                     self.ei = np.zeros(self.ngrids, dtype=np.int32)
-                    self.ti = -1 * np.ones(self.ngrids, dtype=np.int32)
                 super(type(self), self).__init__(*args, **kwargs)
 
             array_class_vdict = {
                 "ngrids": Variable("ngrids", dtype=np.int32, to_write=False, initial=-1),
                 "ei": Variable("ei", dtype=np.int32, to_write=False),
-                "ti": Variable("ti", dtype=np.int32, to_write=False, initial=-1),
                 "__init__": ArrayClass_init,
             }
             array_class = type(class_name, (pclass,), array_class_vdict)
@@ -719,7 +717,6 @@ def from_particlefile(
                 v.name
                 not in [
                     "ei",
-                    "ti",
                     "dt",
                     "depth",
                     "id",

tests/test_fieldset.py

@@ -407,22 +407,6 @@ def test_fieldset_samegrids_from_data():
     assert fieldset1.U.grid == fieldset1.B.grid
 
 
-@pytest.mark.parametrize("dx, dy", [("e1u", "e2u"), ("e1v", "e2v")])
-def test_fieldset_celledgesizes_curvilinear(dx, dy):
-    fname = str(TEST_DATA / "mask_nemo_cross_180lon.nc")
-    filenames = {"dx": fname, "dy": fname, "mesh_mask": fname}
-    variables = {"dx": dx, "dy": dy}
-    dimensions = {"dx": {"lon": "glamu", "lat": "gphiu"}, "dy": {"lon": "glamu", "lat": "gphiu"}}
-    fieldset = FieldSet.from_nemo(filenames, variables, dimensions)
-
-    # explicitly setting cell_edge_sizes from e1u and e2u etc
-    fieldset.dx.grid.cell_edge_sizes["x"] = fieldset.dx.data
-    fieldset.dx.grid.cell_edge_sizes["y"] = fieldset.dy.data
-
-    A = fieldset.dx.cell_areas()
-    assert np.allclose(A, fieldset.dx.data * fieldset.dy.data)
-
-
 def test_fieldset_write_curvilinear(tmpdir):
     fname = str(TEST_DATA / "mask_nemo_cross_180lon.nc")
     filenames = {"dx": fname, "mesh_mask": fname}
@@ -456,21 +440,6 @@ def test_curv_fieldset_add_periodic_halo():
         fieldset.add_periodic_halo(zonal=3, meridional=2)
 
 
-@pytest.mark.parametrize("mesh", ["flat", "spherical"])
-def test_fieldset_cellareas(mesh):
-    data, dimensions = generate_fieldset_data(10, 7)
-    fieldset = FieldSet.from_data(data, dimensions, mesh=mesh)
-    cell_areas = fieldset.V.cell_areas()
-    if mesh == "flat":
-        assert np.allclose(cell_areas.flatten(), cell_areas[0, 0], rtol=1e-3)
-    else:
-        assert all(
-            (np.gradient(cell_areas, axis=0) < 0).flatten()
-        )  # areas should decrease with latitude in spherical mesh
-        for y in range(cell_areas.shape[0]):
-            assert np.allclose(cell_areas[y, :], cell_areas[y, 0], rtol=1e-3)
-
-
 def addConst(particle, fieldset, time):  # pragma: no cover
     particle.lon = particle.lon + fieldset.movewest + fieldset.moveeast
 
@@ -555,7 +524,8 @@ def test_fieldset_write(tmp_zarrfile):
     def UpdateU(particle, fieldset, time):  # pragma: no cover
         tmp1, tmp2 = fieldset.UV[particle]
         _, yi, xi = fieldset.U.unravel_index(particle.ei)
-        fieldset.U.data[particle.ti, yi, xi] += 1
+        ti = fieldset.U._time_index(time)
+        fieldset.U.data[ti, yi, xi] += 1
         fieldset.U.grid.time[0] = time
 
     pset = ParticleSet(fieldset, pclass=Particle, lon=5, lat=5)