Remove add_periodic_halo() methods on Grid, Field, FieldSet

- Remove Grid.meriodonal_halo and Grid.zonal_halo
- Update tests (removal and marking) accordingly

Author: VeckoTheGecko

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

parcels/field.py

@@ -629,12 +629,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), ())))
@@ -760,28 +760,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):
@@ -935,54 +929,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

@@ -1361,26 +1361,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.
 
@@ -1448,9 +1428,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:
@@ -1469,9 +1447,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

@@ -429,17 +429,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
 

tests/test_fieldset_sampling.py

@@ -112,6 +112,8 @@ def test_fieldset_sample_eval(fieldset):
     assert np.allclose(u_s, lat, rtol=1e-5)
 
 
+@pytest.mark.v4remove
+@pytest.mark.xfail(reason="Test is directly testing adding the halo. This test should either be adapted or removed.")
 def test_fieldset_polar_with_halo(fieldset_geometric_polar):
     fieldset_geometric_polar.add_periodic_halo(zonal=5)
     pset = ParticleSet(fieldset_geometric_polar, pclass=pclass(), lon=0, lat=0)

tests/test_grids.py

@@ -149,6 +149,8 @@ def test_avoid_repeated_grids():
     assert fieldset.V.grid is not fieldset.U.grid
 
 
+@pytest.mark.v4alpha
+@pytest.mark.xfail(reason="Calls fieldset.add_periodic_halo(). Should adapt this test case.")
 def test_multigrids_pointer():
     lon_g0 = np.linspace(0, 1e4, 21, dtype=np.float32)
     lat_g0 = np.linspace(0, 1000, 2, dtype=np.float32)

tests/tools/test_warnings.py

@@ -17,17 +17,6 @@
 from tests.utils import TEST_DATA
 
 
-def test_fieldset_warning_halo():
-    # halo with inconsistent boundaries
-    lat = [0, 1, 5, 10]
-    lon = [0, 1, 5, 10]
-    u = [[1, 1, 1, 1] for _ in range(4)]
-    v = [[1, 1, 1, 1] for _ in range(4)]
-    fieldset = FieldSet.from_data(data={"U": u, "V": v}, dimensions={"lon": lon, "lat": lat}, transpose=True)
-    with pytest.warns(FieldSetWarning, match="The zonal halo is located at the east and west of current grid.*"):
-        fieldset.add_periodic_halo(meridional=True, zonal=True)
-
-
 def test_fieldset_warning_latflipped():
     # flipping lats warning
     lat = [0, 1, 5, -5]