merge conflict fixed

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/field.py

@@ -49,7 +49,7 @@
     DeferredNetcdfFileBuffer,
     NetcdfFileBuffer,
 )
-from .grid import Grid, GridType, _calc_cell_areas
+from .grid import Grid, GridType
 
 if TYPE_CHECKING:
     from parcels.fieldset import FieldSet
@@ -337,10 +337,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
@@ -870,13 +866,6 @@ 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, particle=None, search2D=False):
         tau, ti = self._time_index(time)
 

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
@@ -136,10 +135,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,
@@ -633,34 +628,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"]

tests/test_fieldset.py

@@ -423,22 +423,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}
@@ -472,21 +456,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
 

tests/test_grids.py

@@ -19,8 +19,6 @@
     UnitConverter,
     Variable,
 )
-from parcels.grid import Grid, _calc_cell_edge_sizes
-from parcels.tools.converters import TimeConverter
 from tests.utils import TEST_DATA
 
 
@@ -975,27 +973,3 @@ def VelocityInterpolator(particle, fieldset, time):  # pragma: no cover
                 assert np.allclose(pset.Wvel[0], 0, atol=1e-9)
             else:
                 assert np.allclose(pset.Wvel[0], -w * convfactor)
-
-
-@pytest.mark.parametrize(
-    "lon, lat",
-    [
-        (np.arange(0.0, 20.0, 1.0), np.arange(0.0, 10.0, 1.0)),
-    ],
-)
-@pytest.mark.parametrize("mesh", ["flat", "spherical"])
-def test_grid_celledgesizes(lon, lat, mesh):
-    grid = Grid.create_grid(
-        lon=lon, lat=lat, depth=np.array([0]), time=np.array([0]), time_origin=TimeConverter(0), mesh=mesh
-    )
-
-    _calc_cell_edge_sizes(grid)
-    D_meridional = grid.cell_edge_sizes["y"]
-    D_zonal = grid.cell_edge_sizes["x"]
-    assert np.allclose(
-        D_meridional.flatten(), D_meridional[0, 0]
-    )  # all meridional distances should be the same in either mesh
-    if mesh == "flat":
-        assert np.allclose(D_zonal.flatten(), D_zonal[0, 0])  # all zonal distances should be the same in flat mesh
-    else:
-        assert all((np.gradient(D_zonal, axis=0) < 0).flatten())  # zonal distances should decrease in spherical mesh