Merge pull request #1901 from OceanParcels/1898

Remove add_periodic_halo() methods on Grid, Field, FieldSet

Author: VeckoTheGecko

.github/workflows/ci.yml

@@ -76,8 +76,8 @@ jobs:
         with:
           environment-file: environment.yml
       - name: Integration test
-        run: |
-          coverage run -m pytest -v -s --nbval-lax -k "not documentation" --html="${{ matrix.os }}_${{ matrix.python-version }}_integration_test_report.html" --self-contained-html docs/examples
+        run: | # TODO v4: Re-enable `tutorial_periodic_boundaries` notebook
+          coverage run -m pytest -v -s --nbval-lax -k "not documentation and not tutorial_periodic_boundaries" --html="${{ matrix.os }}_${{ matrix.python-version }}_integration_test_report.html" --self-contained-html docs/examples
           coverage xml
       - name: Codecov
         uses: codecov/[email protected]

docs/examples/example_mitgcm.py

@@ -1,6 +1,8 @@
 from datetime import timedelta
 from pathlib import Path
 
+import pytest
+
 import parcels
 
 
@@ -49,6 +51,8 @@ def periodicBC(particle, fieldset, time):  # pragma: no cover
     )
 
 
+@pytest.mark.v4alpha
+@pytest.mark.xfail(reason="Test uses add_periodic_halo(). Update to not use it.")
 def test_mitgcm_output(tmpdir):
     def get_path() -> Path:
         return tmpdir / "MIT_particles.zarr"

docs/examples/example_moving_eddies.py

@@ -250,6 +250,10 @@ def test_moving_eddies_file(mesh, tmpdir):
         assert pset[1].lon < 2.0 and 48.8 < pset[1].lat < 48.85
 
 
+@pytest.mark.v4alpha
+@pytest.mark.xfail(
+    reason="Calls fieldset.add_periodic_halo(). In v4, interpolation should work without adding halo."
+)
 def test_periodic_and_computeTimeChunk_eddies():
     data_folder = parcels.download_example_dataset("MovingEddies_data")
     filename = str(data_folder / "moving_eddies")

docs/examples/tutorial_analyticaladvection.ipynb

@@ -467,9 +467,12 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "fieldsetBJ.add_constant(\"halo_west\", fieldsetBJ.U.grid.lon[0])\n",
-    "fieldsetBJ.add_constant(\"halo_east\", fieldsetBJ.U.grid.lon[-1])\n",
-    "fieldsetBJ.add_periodic_halo(zonal=True)\n",
+    "fieldsetBJ.add_constant(\n",
+    "    \"halo_west\", fieldsetBJ.U.grid.lon[1]\n",
+    ")  # TODO v4: Change index back to 0 when periodic interpolation is done\n",
+    "fieldsetBJ.add_constant(\n",
+    "    \"halo_east\", fieldsetBJ.U.grid.lon[-2]\n",
+    ")  # TODO v4: Change index back to -1 when periodic interpolation is done\n",
     "\n",
     "\n",
     "def ZonalBC(particle, fieldset, time):\n",
@@ -608,7 +611,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "parcels",
+   "display_name": "parcels-dev",
    "language": "python",
    "name": "python3"
   },
@@ -622,7 +625,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.12.3"
+   "version": "3.10.15"
   }
  },
  "nbformat": 4,

docs/v4/api.md

@@ -36,6 +36,8 @@ Here, important things to note are:
 
 - Interpolators (which would implement the `Interpolator` protocol) are responsible for the actual interpolation of the data, and performance considerations. There will be interpolation and indexing utilities that can be made available to the interpolators, allowing for code re-use.
 
+  - Interpolators of the data should handle spatial periodicity and, for the case of rectilinear structured grids, without pre-computing a halo for the FieldSet and Grid ([issue](https://github.com/OceanParcels/Parcels/issues/1898)).
+
 - In the `Field` class, not all combinations of `data`, `grid`, and `interpolator` will logically make sense (e.g., a `xr.DataArray` on a `ux.Grid`, or `ux.DataArray` on a `parcels.Grid`). It's up to the `Interpolator.assert_is_compatible(Field)` to define what is and is not compatible, and raise `ValueError` / `TypeError` on incompatible data types. The `.assert_is_compatible()` method also acts as developer documentation, defining clearly for the `.interpolate()` method what assumptions it is working on. The `.assert_is_compatible()` method should be lightweight as it will be called on `Field` initialisation.
 
 - The `grid` object, in the case of unstructured grids, will be the `Grid` class from UXarray. For structured `Grid`s, it will be an object similar to that of `xgcm.Grid` (note that it will be very different from the v3 `Grid` object hierarchy).

parcels/field.py

@@ -596,12 +596,12 @@ def from_netcdf(
                     if len(buffer_data.shape) == 4:
                         errormessage = (
                             f"Field {filebuffer.name} expecting a data shape of [tdim={grid.tdim}, zdim={grid.zdim}, "
-                            f"ydim={grid.ydim - 2 * grid.meridional_halo}, xdim={grid.xdim - 2 * grid.zonal_halo}] "
+                            f"ydim={grid.ydim}, xdim={grid.xdim }] "
                             f"but got shape {buffer_data.shape}. Flag transpose=True could help to reorder the data."
                         )
                         assert buffer_data.shape[0] == grid.tdim, errormessage
-                        assert buffer_data.shape[2] == grid.ydim - 2 * grid.meridional_halo, errormessage
-                        assert buffer_data.shape[3] == grid.xdim - 2 * grid.zonal_halo, errormessage
+                        assert buffer_data.shape[2] == grid.ydim, errormessage
+                        assert buffer_data.shape[3] == grid.xdim, errormessage
 
                     if len(buffer_data.shape) == 2:
                         data_list.append(buffer_data.reshape(sum(((len(tslice), 1), buffer_data.shape), ())))
@@ -723,28 +723,22 @@ def _reshape(self, data, transpose=False):
                 "Flag transpose=True could help to reorder the data."
             )
             assert data.shape[0] == self.grid.tdim, errormessage
-            assert data.shape[2] == self.grid.ydim - 2 * self.grid.meridional_halo, errormessage
-            assert data.shape[3] == self.grid.xdim - 2 * self.grid.zonal_halo, errormessage
+            assert data.shape[2] == self.grid.ydim, errormessage
+            assert data.shape[3] == self.grid.xdim, errormessage
             if self.gridindexingtype == "pop":
                 assert data.shape[1] == self.grid.zdim or data.shape[1] == self.grid.zdim - 1, errormessage
             else:
                 assert data.shape[1] == self.grid.zdim, errormessage
         else:
             assert data.shape == (
                 self.grid.tdim,
-                self.grid.ydim - 2 * self.grid.meridional_halo,
-                self.grid.xdim - 2 * self.grid.zonal_halo,
+                self.grid.ydim,
+                self.grid.xdim,
             ), (
                 f"Field {self.name} expecting a data shape of [tdim, ydim, xdim]. "
                 "Flag transpose=True could help to reorder the data."
             )
-        if self.grid.meridional_halo > 0 or self.grid.zonal_halo > 0:
-            data = self.add_periodic_halo(
-                zonal=self.grid.zonal_halo > 0,
-                meridional=self.grid.meridional_halo > 0,
-                halosize=max(self.grid.meridional_halo, self.grid.zonal_halo),
-                data=data,
-            )
+
         return data
 
     def set_scaling_factor(self, factor):
@@ -898,54 +892,6 @@ def eval(self, time, z, y, x, particle=None, applyConversion=True):
         else:
             return value
 
-    def add_periodic_halo(self, zonal, meridional, halosize=5, data=None):
-        """Add a 'halo' to all Fields in a FieldSet.
-
-        Add a 'halo' to all Fields in a FieldSet, through extending the Field (and lon/lat)
-        by copying a small portion of the field on one side of the domain to the other.
-        Before adding a periodic halo to the Field, it has to be added to the Grid on which the Field depends
-
-        See `this tutorial <../examples/tutorial_periodic_boundaries.ipynb>`__
-        for a detailed explanation on how to set up periodic boundaries
-
-        Parameters
-        ----------
-        zonal : bool
-            Create a halo in zonal direction.
-        meridional : bool
-            Create a halo in meridional direction.
-        halosize : int
-            Size of the halo (in grid points). Default is 5 grid points
-        data :
-            if data is not None, the periodic halo will be achieved on data instead of self.data and data will be returned (Default value = None)
-        """
-        dataNone = not isinstance(data, np.ndarray)
-        if self.grid.defer_load and dataNone:
-            return
-        data = self.data if dataNone else data
-        if zonal:
-            if len(data.shape) == 3:
-                data = np.concatenate((data[:, :, -halosize:], data, data[:, :, 0:halosize]), axis=len(data.shape) - 1)
-                assert data.shape[2] == self.grid.xdim, "Third dim must be x."
-            else:
-                data = np.concatenate(
-                    (data[:, :, :, -halosize:], data, data[:, :, :, 0:halosize]), axis=len(data.shape) - 1
-                )
-                assert data.shape[3] == self.grid.xdim, "Fourth dim must be x."
-        if meridional:
-            if len(data.shape) == 3:
-                data = np.concatenate((data[:, -halosize:, :], data, data[:, 0:halosize, :]), axis=len(data.shape) - 2)
-                assert data.shape[1] == self.grid.ydim, "Second dim must be y."
-            else:
-                data = np.concatenate(
-                    (data[:, :, -halosize:, :], data, data[:, :, 0:halosize, :]), axis=len(data.shape) - 2
-                )
-                assert data.shape[2] == self.grid.ydim, "Third dim must be y."
-        if dataNone:
-            self.data = data
-        else:
-            return data
-
     def write(self, filename, varname=None):
         """Write a :class:`Field` to a netcdf file.
 

parcels/fieldset.py

@@ -1358,26 +1358,6 @@ def add_constant(self, name, value):
         """
         setattr(self, name, value)
 
-    def add_periodic_halo(self, zonal=False, meridional=False, halosize=5):
-        """Add a 'halo' to all :class:`parcels.field.Field` objects in a FieldSet,
-        through extending the Field (and lon/lat) by copying a small portion
-        of the field on one side of the domain to the other.
-
-        Parameters
-        ----------
-        zonal : bool
-            Create a halo in zonal direction (Default value = False)
-        meridional : bool
-            Create a halo in meridional direction (Default value = False)
-        halosize : int
-            size of the halo (in grid points). Default is 5 grid points
-        """
-        for grid in self.gridset.grids:
-            grid.add_periodic_halo(zonal, meridional, halosize)
-        for value in self.__dict__.values():
-            if isinstance(value, Field):
-                value.add_periodic_halo(zonal, meridional, halosize)
-
     def write(self, filename):
         """Write FieldSet to NetCDF file using NEMO convention.
 
@@ -1445,9 +1425,7 @@ def computeTimeChunk(self, time=0.0, dt=1):
                     zd = g.zdim - 1
                 else:
                     zd = g.zdim
-                data = np.empty(
-                    (g.tdim, zd, g.ydim - 2 * g.meridional_halo, g.xdim - 2 * g.zonal_halo), dtype=np.float32
-                )
+                data = np.empty((g.tdim, zd, g.ydim, g.xdim), dtype=np.float32)
                 f._loaded_time_indices = range(2)
                 for tind in f._loaded_time_indices:
                     for fb in f.filebuffers:
@@ -1466,9 +1444,7 @@ def computeTimeChunk(self, time=0.0, dt=1):
                     zd = g.zdim - 1
                 else:
                     zd = g.zdim
-                data = np.empty(
-                    (g.tdim, zd, g.ydim - 2 * g.meridional_halo, g.xdim - 2 * g.zonal_halo), dtype=np.float32
-                )
+                data = np.empty((g.tdim, zd, g.ydim, g.xdim), dtype=np.float32)
                 if signdt >= 0:
                     f._loaded_time_indices = [1]
                     if f.filebuffers[0] is not None:

parcels/grid.py

@@ -68,8 +68,6 @@ def __init__(
         assert_valid_mesh(mesh)
         self._mesh = mesh
         self._zonal_periodic = False
-        self._zonal_halo = 0
-        self._meridional_halo = 0
         self._lat_flipped = False
         self._defer_load = False
         self._lonlat_minmax = np.array(
@@ -108,10 +106,6 @@ def negate_depth(self):
     def mesh(self):
         return self._mesh
 
-    @property
-    def meridional_halo(self):
-        return self._meridional_halo
-
     @property
     def lonlat_minmax(self):
         return self._lonlat_minmax
@@ -124,10 +118,6 @@ def time_origin(self):
     def zonal_periodic(self):
         return self._zonal_periodic
 
-    @property
-    def zonal_halo(self):
-        return self._zonal_halo
-
     @property
     def defer_load(self):
         return self._defer_load
@@ -279,50 +269,6 @@ def xdim(self):
     def ydim(self):
         return self.lat.size
 
-    def add_periodic_halo(self, zonal: bool, meridional: bool, halosize: int = 5):
-        """Add a 'halo' to the Grid, through extending the Grid (and lon/lat)
-        similarly to the halo created for the Fields
-
-        Parameters
-        ----------
-        zonal : bool
-            Create a halo in zonal direction
-        meridional : bool
-            Create a halo in meridional direction
-        halosize : int
-            size of the halo (in grid points). Default is 5 grid points
-        """
-        if zonal:
-            lonshift = self.lon[-1] - 2 * self.lon[0] + self.lon[1]
-            if not np.allclose(self.lon[1] - self.lon[0], self.lon[-1] - self.lon[-2]):
-                warnings.warn(
-                    "The zonal halo is located at the east and west of current grid, "
-                    "with a dx = lon[1]-lon[0] between the last nodes of the original grid and the first ones of the halo. "
-                    "In your grid, lon[1]-lon[0] != lon[-1]-lon[-2]. Is the halo computed as you expect?",
-                    FieldSetWarning,
-                    stacklevel=2,
-                )
-            self._lon = np.concatenate((self.lon[-halosize:] - lonshift, self.lon, self.lon[0:halosize] + lonshift))
-            self._zonal_periodic = True
-            self._zonal_halo = halosize
-        if meridional:
-            if not np.allclose(self.lat[1] - self.lat[0], self.lat[-1] - self.lat[-2]):
-                warnings.warn(
-                    "The meridional halo is located at the north and south of current grid, "
-                    "with a dy = lat[1]-lat[0] between the last nodes of the original grid and the first ones of the halo. "
-                    "In your grid, lat[1]-lat[0] != lat[-1]-lat[-2]. Is the halo computed as you expect?",
-                    FieldSetWarning,
-                    stacklevel=2,
-                )
-            latshift = self.lat[-1] - 2 * self.lat[0] + self.lat[1]
-            self._lat = np.concatenate((self.lat[-halosize:] - latshift, self.lat, self.lat[0:halosize] + latshift))
-            self._meridional_halo = halosize
-        self._lonlat_minmax = np.array(
-            [np.nanmin(self.lon), np.nanmax(self.lon), np.nanmin(self.lat), np.nanmax(self.lat)], dtype=np.float32
-        )
-        if isinstance(self, RectilinearSGrid):
-            self._add_Sdepth_periodic_halo(zonal, meridional, halosize)
-
 
 class RectilinearZGrid(RectilinearGrid):
     """Rectilinear Z Grid.
@@ -469,23 +415,6 @@ def xdim(self):
     def ydim(self):
         return self.lon.shape[0]
 
-    def add_periodic_halo(self, zonal, meridional, halosize=5):
-        """Add a 'halo' to the Grid, through extending the Grid (and lon/lat)
-        similarly to the halo created for the Fields
-
-        Parameters
-        ----------
-        zonal : bool
-            Create a halo in zonal direction
-        meridional : bool
-            Create a halo in meridional direction
-        halosize : int
-            size of the halo (in grid points). Default is 5 grid points
-        """
-        raise NotImplementedError(
-            "CurvilinearGrid does not support add_periodic_halo. See https://github.com/OceanParcels/Parcels/pull/1811"
-        )
-
 
 class CurvilinearZGrid(CurvilinearGrid):
     """Curvilinear Z Grid.

tests/test_advection.py

@@ -275,6 +275,8 @@ def periodicBC(particle, fieldset, time):  # pragma: no cover
     particle.lat = math.fmod(particle.lat, 1)
 
 
+@pytest.mark.v4alpha
+@pytest.mark.xfail(reason="Calls fieldset.add_periodic_halo(). In v4, interpolation should work without adding halo.")
 def test_advection_periodic_zonal():
     xdim, ydim, halosize = 100, 100, 3
     fieldset = create_periodic_fieldset(xdim, ydim, uvel=1.0, vvel=0.0)
@@ -286,6 +288,8 @@ def test_advection_periodic_zonal():
     assert abs(pset.lon[0] - 0.15) < 0.1
 
 
+@pytest.mark.v4alpha
+@pytest.mark.xfail(reason="Calls fieldset.add_periodic_halo(). In v4, interpolation should work without adding halo.")
 def test_advection_periodic_meridional():
     xdim, ydim = 100, 100
     fieldset = create_periodic_fieldset(xdim, ydim, uvel=0.0, vvel=1.0)
@@ -297,6 +301,8 @@ def test_advection_periodic_meridional():
     assert abs(pset.lat[0] - 0.15) < 0.1
 
 
+@pytest.mark.v4alpha
+@pytest.mark.xfail(reason="Calls fieldset.add_periodic_halo(). In v4, interpolation should work without adding halo.")
 def test_advection_periodic_zonal_meridional():
     xdim, ydim = 100, 100
     fieldset = create_periodic_fieldset(xdim, ydim, uvel=1.0, vvel=1.0)

tests/test_fieldset.py

@@ -397,17 +397,6 @@ def test_fieldset_write_curvilinear(tmpdir):
         assert np.allclose(getattr(fieldset2.dx, var), getattr(fieldset.dx, var))
 
 
-def test_curv_fieldset_add_periodic_halo():
-    fname = str(TEST_DATA / "mask_nemo_cross_180lon.nc")
-    filenames = {"dx": fname, "dy": fname, "mesh_mask": fname}
-    variables = {"dx": "e1u", "dy": "e1v"}
-    dimensions = {"dx": {"lon": "glamu", "lat": "gphiu"}, "dy": {"lon": "glamu", "lat": "gphiu"}}
-    fieldset = FieldSet.from_nemo(filenames, variables, dimensions)
-
-    with pytest.raises(NotImplementedError):
-        fieldset.add_periodic_halo(zonal=3, meridional=2)
-
-
 def addConst(particle, fieldset, time):  # pragma: no cover
     particle.lon = particle.lon + fieldset.movewest + fieldset.moveeast