Creating the boundary forcing

Creating the boundary forcing#

[1]:

from roms_tools import Grid, BoundaryForcing

As always, the first step is to create our grid. In this example, we choose a grid that contains Greenland’s south east coast and the Irminger Sea. The grid has a horizontal resolution of 2 km and 100 vertical layers.

[2]:

grid = Grid(
    nx=500,
    ny=500,
    size_x=1000,
    size_y=1000,
    center_lon=-37,
    center_lat=63.5,
    rot=25,
    N=100,
    topography_source={
        "name": "SRTM15",
        "path": "/anvil/projects/x-ees250129/Datasets/SRTM15/SRTM15_V2.6.nc",
    },
)

[3]:

grid.plot()

Note that the northern boundary lies entirely on land, so we can exclude it from the boundary forcing generation.

Next, we specify the temporal range that we want to make the surface forcing for.

[3]:

from datetime import datetime

[4]:

start_time = datetime(2012, 1, 30)
end_time = datetime(2012, 2, 4)

ROMS-Tools can create two types of boundary forcing:

physical boundary forcing like temperature, salinity, velocities, and sea surface height
biogeochemical (BGC) boundary forcing like alkalinity, dissolved inorganic phosphate, etc.

As with surface forcing, ROMS accepts multiple boundary forcing files, so we create these two types separately.

Saving as NetCDF or YAML file#

We can now save our boundary forcings as a NetCDF files.

We need to specify the desired target path.

[24]:

filepath = "/anvil/scratch/x-kthyng/forcing/my_bgc_boundary_forcing.nc"

[25]:

%time unified_bgc_boundary_forcing.save(filepath)

2026-06-29 18:22:44 - INFO - Writing the following NetCDF files:
/anvil/scratch/x-kthyng/forcing/my_bgc_boundary_forcing_clim.nc

[########################################] | 100% Completed | 7.36 ss
CPU times: user 6.94 s, sys: 2.39 s, total: 9.32 s
Wall time: 8.05 s

[25]:

[PosixPath('/anvil/scratch/x-kthyng/forcing/my_bgc_boundary_forcing_clim.nc')]

Next, let’s save the physical boundary forcing. We will use the group = True flag (default) which ensures that the data is split into multiple NetCDFs, grouped into months.

[26]:

filepath = "/anvil/scratch/x-kthyng/forcing/my_boundary_forcing.nc"

[27]:

%time boundary_forcing.save(filepath, group=True)

2026-06-29 18:22:52 - INFO - Writing the following NetCDF files:
/anvil/scratch/x-kthyng/forcing/my_boundary_forcing_201201.nc
/anvil/scratch/x-kthyng/forcing/my_boundary_forcing_201202.nc

[########################################] | 100% Completed | 7.06 ss
[########################################] | 100% Completed | 19.01 s
CPU times: user 3.79 s, sys: 4.21 s, total: 8 s
Wall time: 27.4 s

[27]:

[PosixPath('/anvil/scratch/x-kthyng/forcing/my_boundary_forcing_201201.nc'),
 PosixPath('/anvil/scratch/x-kthyng/forcing/my_boundary_forcing_201202.nc')]

From the file paths printed to the screen, you will notice that the year and month information has been appended to the specified path.

Note

The way data is grouped depends on its temporal frequency. For datasets with a daily or higher frequency, the data is grouped by month. For datasets with lower temporal frequency (e.g., monthly data), the grouping is done by year.

We can also export the parameters of our BoundaryForcing object to a YAML file.

[31]:

yaml_filepath = "/anvil/scratch/x-kthyng/forcing/my_boundary_forcing.yaml"

[32]:

boundary_forcing.to_yaml(yaml_filepath)

This is the YAML file that was created.

[33]:

# Open and read the YAML file
with open(yaml_filepath, "r") as file:
    file_contents = file.read()

# Print the contents
print(file_contents)

---
roms_tools_version: 10000.dev377+g7f6dd1f62
---
Grid:
  nx: 500
  ny: 500
  size_x: 1000
  size_y: 1000
  center_lon: -37
  center_lat: 63.5
  rot: 25
  N: 100
  theta_s: 5.0
  theta_b: 2.0
  hc: 300.0
  topography_source:
    name: SRTM15
    path: /anvil/projects/x-ees250129/Datasets/SRTM15/SRTM15_V2.6.nc
  mask_shapefile: null
  close_narrow_channels: false
  hmin: 5.0
  filename: null
BoundaryForcing:
  start_time: '2012-01-30T00:00:00'
  end_time: '2012-02-04T00:00:00'
  boundaries:
    south: true
    east: true
    north: false
    west: true
  source:
    name: GLORYS
    path:
    - /anvil/projects/x-ees250129/cstar-forge-data/source-data/GLORYS_GLOBAL/cmems_mod_glo_phy_my_0.083deg_P1D-m_GLOBAL_20120129.nc
    - /anvil/projects/x-ees250129/cstar-forge-data/source-data/GLORYS_GLOBAL/cmems_mod_glo_phy_my_0.083deg_P1D-m_GLOBAL_20120130.nc
    - /anvil/projects/x-ees250129/cstar-forge-data/source-data/GLORYS_GLOBAL/cmems_mod_glo_phy_my_0.083deg_P1D-m_GLOBAL_20120131.nc
    - /anvil/projects/x-ees250129/cstar-forge-data/source-data/GLORYS_GLOBAL/cmems_mod_glo_phy_my_0.083deg_P1D-m_GLOBAL_20120201.nc
    - /anvil/projects/x-ees250129/cstar-forge-data/source-data/GLORYS_GLOBAL/cmems_mod_glo_phy_my_0.083deg_P1D-m_GLOBAL_20120202.nc
    - /anvil/projects/x-ees250129/cstar-forge-data/source-data/GLORYS_GLOBAL/cmems_mod_glo_phy_my_0.083deg_P1D-m_GLOBAL_20120203.nc
    - /anvil/projects/x-ees250129/cstar-forge-data/source-data/GLORYS_GLOBAL/cmems_mod_glo_phy_my_0.083deg_P1D-m_GLOBAL_20120204.nc
    - /anvil/projects/x-ees250129/cstar-forge-data/source-data/GLORYS_GLOBAL/cmems_mod_glo_phy_my_0.083deg_P1D-m_GLOBAL_20120205.nc
    climatology: false
  type: physics
  prefill: null
  prefill_kwargs: null
  regrid_method: null
  extrap_method: null
  extrap_kwargs: null
  model_reference_date: '2000-01-01T00:00:00'
  chunks: null
  initial_slice_bounds: null
  bypass_validation: false
  bgc_interpolation_method: depth

Horizontal regridding near land and coarse source data#

The source data (such as GLORYS) and the ROMS grid generally have different land masks, often because they have different resolutions. A naive linear interpolation would propagate NaNs into the ROMS boundary wherever any of the surrounding source cells is land, leaving holes at ocean points that ROMS needs.

ROMS-Tools avoids this with a NaN-aware horizontal regridding that needs no separate fill step:

When xESMF is available (the default for conda installs), boundary data is regridded with masked bilinear interpolation: the interpolation weights are renormalized over only the valid (ocean) source cells, so a target point surrounded by mostly land still receives a sensible ocean value. Target points whose source neighbors are all land are then filled by inverse-distance-weighted extrapolation (configurable via the extrap_method argument). Together these guarantee NaN-free boundaries.
When xESMF is not available (e.g. Windows/pip installs), ROMS-Tools falls back to a nearest-neighbor pre-fill of the source data followed by scipy linear interpolation. Because every masked source cell is first filled with its nearest valid ocean value, the subsequent interpolation cannot propagate NaNs. This is numerically close to, though not identical to, the masked-bilinear result.

This replaces the older, slower 2D Poisson (AMG) fill that was previously required to plug NaN holes before regridding, giving a substantial speedup and memory reduction on large, multi-year domains.

You can also choose the regrid engine explicitly with the regrid_method argument ("auto" (default; xESMF if installed, otherwise scipy), "xesmf", or "scipy"), independently of any source pre-fill, and tune the extrapolation with extrap_kwargs (e.g. num_src_pnts and dist_exponent for "inverse_dist").

Pre-filling the source data#

The default NaN-aware regridding above interpolates directly from the valid ocean cells of the source. In some cases you may instead want to fill the source’s land/void cells first and then regrid the now NaN-free source. This is controlled by the prefill argument (with method-specific options in prefill_kwargs):

Strategy	`prefill` value	Needs xESMF	Cost	Character
Masked bilinear + extrapolation (default)	`None`	yes*	low	Interpolates over real ocean only; `extrap_method` fills boundary points with no valid source neighbor. *Falls back to a nearest-neighbor pre-fill + scipy when xESMF is unavailable.
AMG Poisson fill (legacy)	`"2d_lateral_fill"`	no	high	Smoothest; a global harmonic fill of the whole source. Slower and more memory-hungry on large domains.
Inverse-distance fill	`"inverse_dist"`	yes	medium	Smooth inverse-distance-weighted fill of the source; tune with `num_src_pnts` / `dist_exponent`.
Nearest-source fill	`"nearest_s2d"`	yes	low	Blocky (single nearest source cell).
Nearest-neighbor fill	`"nearest_neighbor"`	no	lowest	Blocky and connectivity-blind; also the automatic fallback when xESMF is unavailable.

There are two independent stages here: ``prefill`` fills the source before regridding (after which the regrid is a plain bilinear interpolation, so extrap_method no longer applies), while ``extrap_method`` fills destination boundary points during the default masked-bilinear regrid. The regrid engine itself is selected with regrid_method ("auto", "xesmf", "scipy").

Let’s visualize how each method fills the coarse CESM source over the Wales domain. We crop the source to the domain (plus a buffer) and apply each fill in turn:

[42]:

from roms_tools.datasets.lat_lon_datasets import CESMBGCDataset
from roms_tools.setup.utils import get_target_coords

target_coords = get_target_coords(grid)
source = CESMBGCDataset(
    filename=cesm_bgc_path,
    start_time=start_time,
    end_time=end_time,
    climatology=True,
    use_dask=True,
)

lon_name = source.dim_names["longitude"]
lat_name = source.dim_names["latitude"]

methods = ["2d_lateral_fill", "inverse_dist", "nearest_s2d", "nearest_neighbor"]
fig, axs = plt.subplots(2, 2, figsize=(11, 9), constrained_layout=True)

for method, ax in zip(methods, axs.ravel()):
    # A fresh copy per method so each fill starts from the raw (NaN) source.
    sub = source.choose_subdomain(
        target_coords, return_copy=True, unchunk_lateral_dims=True
    )
    sub.convert_to_float64()
    sub.extrapolate_deepest_to_bottom()   # fill deep NaNs down each column first
    sub.apply_prefill(method)

    field = sub.ds["ALK"]
    # reduce to a single (lat, lon) surface slice for plotting
    extra_dims = [d for d in field.dims if d not in (lon_name, lat_name)]
    field.isel({d: 0 for d in extra_dims}).plot(ax=ax, x=lon_name, y=lat_name)
    ax.set_title(f"prefill='{method}'")

The AMG (2d_lateral_fill) and inverse-distance fills produce smooth fields, while the two nearest-style fills are blockier. For inverse_dist you can trade smoothness against locality via prefill_kwargs (more source points and a smaller distance exponent give a smoother fill):

[43]:

sub = source.choose_subdomain(
    target_coords, return_copy=True, unchunk_lateral_dims=True
)
sub.convert_to_float64()
sub.apply_prefill(
    "inverse_dist", prefill_kwargs={"num_src_pnts": 16, "dist_exponent": 1.0}
)

field = sub.ds["ALK"]
extra_dims = [d for d in field.dims if d not in (lon_name, lat_name)]
field.isel({d: 0 for d in extra_dims}).plot(x=lon_name, y=lat_name)

[43]:

<matplotlib.collections.QuadMesh at 0x14717f99c910>

To use a pre-fill when building the physics boundary forcing itself, pass prefill (and optionally prefill_kwargs and regrid_method) to BoundaryForcing:

bf = BoundaryForcing(
    grid=grid,
    start_time=start_time,
    end_time=end_time,
    boundaries={
        "south": True,
        "east": True,
        "north": False,  # northern boundary excluded
        "west": True,
    },
    source={"name": "GLORYS", "path": glorys_path},
    type="physics",  # "physics" or "bgc"; default is "physics"
    prefill="inverse_dist",
    prefill_kwargs={"num_src_pnts": 16},
    model_reference_date=datetime(2000, 1, 1), # this is the default
    use_dask=True,
)

Pass regrid_method="scipy" to reproduce pre-v4 outputs byte-for-byte (the legacy apply_2d_horizontal_fill=True is now equivalent to prefill="2d_lateral_fill", regrid_method="scipy").

[ ]:

Creating the boundary forcing

Contents

Creating the boundary forcing#

Physical boundary forcing#

Biogeochemical (BGC) boundary forcing#

Interpolate BGC variables by density or mixed-layer depth instead of plain depth#

Saving as NetCDF or YAML file#

Creating boundary forcing from an existing YAML file#

Horizontal regridding near land and coarse source data#

Coarse resolution of the source data#

Pre-filling the source data#