Fix mask at boundaries (#165)

* Test that mask and topography along grid boundaries are identical to adjacent inland cells

* Handle mask and topography boundaries correctly

* Increase resolution of test grid due to modified mask

* Overwrite test data due to modified mask

roms_tools/setup/topography.py

@@ -61,6 +61,10 @@ def _add_topography_and_mask(
 
     # fill enclosed basins with land
     mask = _fill_enclosed_basins(mask.values)
+
+    # adjust mask boundaries by copying values from adjacent cells
+    mask = _handle_boundaries(mask)
+
     ds["mask_rho"] = xr.DataArray(mask.astype(np.int32), dims=("eta_rho", "xi_rho"))
     ds["mask_rho"].attrs = {
         "long_name": "Mask at rho-points",
@@ -233,10 +237,7 @@ def _smooth_topography_locally(h, hmin=5, rmax=0.2):
         )
 
         # No gradient at the domain boundaries
-        h_log[0, :] = h_log[1, :]
-        h_log[-1, :] = h_log[-2, :]
-        h_log[:, 0] = h_log[:, 1]
-        h_log[:, -1] = h_log[:, -2]
+        h_log = _handle_boundaries(h_log)
 
         # Update h
         h = hmin * np.exp(h_log)
@@ -252,6 +253,31 @@ def _smooth_topography_locally(h, hmin=5, rmax=0.2):
     return h
 
 
+def _handle_boundaries(field):
+    """
+    Adjust the boundaries of a 2D field by copying values from adjacent cells.
+
+    Parameters
+    ----------
+    field : numpy.ndarray or xarray.DataArray
+        A 2D array representing a field (e.g., topography or mask) whose boundary values
+        need to be adjusted.
+
+    Returns
+    -------
+    field : numpy.ndarray or xarray.DataArray
+        The input field with adjusted boundary values.
+
+    """
+
+    field[0, :] = field[1, :]
+    field[-1, :] = field[-2, :]
+    field[:, 0] = field[:, 1]
+    field[:, -1] = field[:, -2]
+
+    return field
+
+
 def _compute_rfactor(h):
     """
     Computes slope parameter (or r-factor) r = |Delta h| / 2h in both horizontal grid directions.

roms_tools/tests/test_setup/test_surface_forcing.py

@@ -13,8 +13,8 @@ def grid_that_straddles_dateline():
     Fixture for creating a domain that straddles the dateline and lies within the bounds of the regional ERA5 data.
     """
     grid = Grid(
-        nx=5,
-        ny=5,
+        nx=20,
+        ny=20,
         size_x=1800,
         size_y=2400,
         center_lon=-10,

roms_tools/tests/test_setup/test_topography.py

@@ -56,3 +56,43 @@ def test_hmin_criterion():
 
     assert grid.hmin == 10.0
     assert np.less_equal(grid.hmin, grid.ds.h.min())
+
+
+def test_mask_topography_boundary():
+    """
+    Test that the mask and topography along the grid boundaries (north, south, east, west)
+    are identical to the adjacent inland cells.
+    """
+
+    # Create a grid with some land along the northern boundary
+    grid = Grid(
+        nx=10, ny=10, size_x=1000, size_y=1000, center_lon=-20, center_lat=60, rot=0
+    )
+
+    # Toopography
+    np.testing.assert_array_equal(
+        grid.ds.h.isel(eta_rho=0).data, grid.ds.h.isel(eta_rho=1).data
+    )
+    np.testing.assert_array_equal(
+        grid.ds.h.isel(eta_rho=-1).data, grid.ds.h.isel(eta_rho=-2).data
+    )
+    np.testing.assert_array_equal(
+        grid.ds.h.isel(xi_rho=0).data, grid.ds.h.isel(xi_rho=1).data
+    )
+    np.testing.assert_array_equal(
+        grid.ds.h.isel(xi_rho=-1).data, grid.ds.h.isel(xi_rho=-2).data
+    )
+
+    # Mask
+    np.testing.assert_array_equal(
+        grid.ds.mask_rho.isel(eta_rho=0).data, grid.ds.mask_rho.isel(eta_rho=1).data
+    )
+    np.testing.assert_array_equal(
+        grid.ds.mask_rho.isel(eta_rho=-1).data, grid.ds.mask_rho.isel(eta_rho=-2).data
+    )
+    np.testing.assert_array_equal(
+        grid.ds.mask_rho.isel(xi_rho=0).data, grid.ds.mask_rho.isel(xi_rho=1).data
+    )
+    np.testing.assert_array_equal(
+        grid.ds.mask_rho.isel(xi_rho=-1).data, grid.ds.mask_rho.isel(xi_rho=-2).data
+    )