Rompy and core features¶

Rompy is a python package for generating configuration files and forcing data for numerical wave and wave-climate models.

Some Components in rompy¶

1. ModelRun

• Purpose: Combines all components and generates a simulation workspace ready to run
• Key Attributes:
  • run_id
  • period (TimeRange)
  • config
  • output_dir
• Methods:
  • run(): Executes the model simulation

2. Grid

• Purpose: Defines the spatial domain and resolution
• Types:
  • RegularGrid: Uniform spacing in x and y directions
  • UnstructuredGrid: Non-uniform mesh (if supported)
• Key Attributes: x0, y0, nx, ny, dx, dy

3. TimeRange

• Purpose: Specifies the temporal domain of the simulation
• Key Attributes: start, end, interval

4. Forcing

• Purpose: Provides input data for the model
• Types:
  • e.g SwanDataGrid: For SWAN model inputs
  • Other model-specific data types
• Components:
  • e.g Bathymetry
  • e.g Wind
  • e.g Boundary conditions
  • Other environmental variables

5. Config

• Purpose: Holds model-specific configuration settings
• Types:
  • e.g SwanConfig: For SWAN model
  • Other model-specific configs
• Key Components:
  • e.g Physics parameters
  • e.g Numerical scheme settings
  • e.g Output specifications

6. Output

• Purpose: Defines what data to save from the model run
• Components:
  • e.g Variables to output
  • e.g Output locations
  • e.g Output frequency

The base objects contained in the rompy packages are designed to provide scaffolding for building model specific implemntation. However, there is some base functionality that we can demonstrate there using these base objects before looking at how they are used in a model implementation.

In [1]:

from pathlib import Path
import xarray as xr
import cartopy.crs as ccrs
import matplotlib.pyplot as plt

workdir = Path.cwd() / "basedemo"
workdir.mkdir(exist_ok=True)


def print_new_contents(path, old_contents=None):
    """
    Helper function to display directory contents and highlight new files.

    Parameters:
    -----------
    path : Path
        Directory path to inspect
    old_contents : set, optional
        Previously existing file names for comparison

    Returns:
    --------
    set : Set of current file names in the directory
    """
    print(f"\nContents of {path}:")
    for item in path.iterdir():
        if item.is_dir():
            print(f" - {item.name}/")
            print_new_contents(item, old_contents=old_contents)
        if old_contents is None or item.name not in old_contents:
            print(f" - {item.name} (new)")
        else:
            print(f" - {item.name}")
    return set(item.name for item in path.iterdir())


# Initialize workspace contents tracking
contents = print_new_contents(workdir)


def plot_nc(
    ncfile,
    varname,
    timestep=None,
    cmap=None,
    grid=None,
):
    """Helper function to plot a variable from a NetCDF file."""
    # display(ncfile)
    fig = plt.figure(figsize=(10, 5))
    ax = plt.axes(projection=ccrs.PlateCarree())
    ds = xr.open_dataset(ncfile)
    if timestep is not None:
        p = ds[varname].isel(time=timestep).plot(ax=ax, cmap=cmap)
    else:
        p = ds[varname].plot(ax=ax, cmap=cmap)
    ax.coastlines()
    if grid is not None:
        grid.plot(ax=ax)
    plt.title(f"{varname} from {ncfile.name}")
    return p


def plot_spectra(spectra, grid=None, markersize=6):
    """Helper function to plot wave spectra from a NetCDF file."""
    ds = xr.open_dataset(spectra)
    # display(ds)
    fig, ax = plt.subplots(
        figsize=(10, 5), subplot_kw={"projection": ccrs.PlateCarree()}
    )
    if grid:
        grid.plot(ax=ax)
    c = ds.isel(time=[0]).spec.hs()
    p = ax.scatter(
        ds.lon,
        ds.lat,
        s=markersize,
        c=c,
        cmap="turbo",
        vmin=0,
        vmax=4,
        transform=ccrs.PlateCarree(),
    )
    ax.set_title(c.time.to_index().to_pydatetime()[0])
    plt.colorbar(p, label=f"Hs (m)")
    ax.coastlines()
    ax.grid(True)

from pathlib import Path import xarray as xr import cartopy.crs as ccrs import matplotlib.pyplot as plt workdir = Path.cwd() / "basedemo" workdir.mkdir(exist_ok=True) def print_new_contents(path, old_contents=None): """ Helper function to display directory contents and highlight new files. Parameters: ----------- path : Path Directory path to inspect old_contents : set, optional Previously existing file names for comparison Returns: -------- set : Set of current file names in the directory """ print(f"\nContents of {path}:") for item in path.iterdir(): if item.is_dir(): print(f" - {item.name}/") print_new_contents(item, old_contents=old_contents) if old_contents is None or item.name not in old_contents: print(f" - {item.name} (new)") else: print(f" - {item.name}") return set(item.name for item in path.iterdir()) # Initialize workspace contents tracking contents = print_new_contents(workdir) def plot_nc( ncfile, varname, timestep=None, cmap=None, grid=None, ): """Helper function to plot a variable from a NetCDF file.""" # display(ncfile) fig = plt.figure(figsize=(10, 5)) ax = plt.axes(projection=ccrs.PlateCarree()) ds = xr.open_dataset(ncfile) if timestep is not None: p = ds[varname].isel(time=timestep).plot(ax=ax, cmap=cmap) else: p = ds[varname].plot(ax=ax, cmap=cmap) ax.coastlines() if grid is not None: grid.plot(ax=ax) plt.title(f"{varname} from {ncfile.name}") return p def plot_spectra(spectra, grid=None, markersize=6): """Helper function to plot wave spectra from a NetCDF file.""" ds = xr.open_dataset(spectra) # display(ds) fig, ax = plt.subplots( figsize=(10, 5), subplot_kw={"projection": ccrs.PlateCarree()} ) if grid: grid.plot(ax=ax) c = ds.isel(time=[0]).spec.hs() p = ax.scatter( ds.lon, ds.lat, s=markersize, c=c, cmap="turbo", vmin=0, vmax=4, transform=ccrs.PlateCarree(), ) ax.set_title(c.time.to_index().to_pydatetime()[0]) plt.colorbar(p, label=f"Hs (m)") ax.coastlines() ax.grid(True)

Contents of /home/tdurrant/source/rompy/rompy-meta/repos/rompy-notebooks/notebooks/common/basedemo:
 - gebco_bathymetry.nc (new)
 - era5_wind.nc (new)
 - ww3_boundary.nc (new)

In [ ]:

TO demonstrate some of the core features of rompy, we will use some sample data files. These files are stored in the tests/data directory of the rompy repository. Here we will be using the following files:

• gebco-1deg.nc: Global bathymetry data from GEBCO
• era5-20230101.nc: A course sample of ERA5 atmospheric data for January 1, 2023
• aus-20230101.nc: Sample wave spectra data from WW3, subset around Australia

In [2]:

DATADIR = Path("../../tests/data/")
# display(sorted(DATADIR.glob("*")))

gebco = DATADIR / "gebco-1deg.nc"
era5 = DATADIR / "era5-20230101.nc"
spectra = DATADIR / "aus-20230101.nc"

for f in [gebco, era5, spectra]:
    assert f.exists(), f"Required data file {f} not found!"

DATADIR = Path("../../tests/data/") # display(sorted(DATADIR.glob("*"))) gebco = DATADIR / "gebco-1deg.nc" era5 = DATADIR / "era5-20230101.nc" spectra = DATADIR / "aus-20230101.nc" for f in [gebco, era5, spectra]: assert f.exists(), f"Required data file {f} not found!"

2. Grid

• Purpose: Defines the spatial domain and resolution
• Types:
  • RegularGrid: Uniform spacing in x and y directions
  • UnstructuredGrid: Non-uniform mesh (if supported)
• Key Attributes: x0, y0, nx, ny, dx, dy

For this demonstration, we will create a RegularGrid off the west coast of Australia.

In [3]:

from rompy.core.grid import RegularGrid

grid = RegularGrid(
    x0=110.0,  # Longitude of the lower-left corner
    y0=-35.2,  # Latitude of the lower-left corner
    rot=4.0,  # Rotation angle in degrees
    dx=0.5,  # Grid spacing in the x-direction (degrees)
    dy=0.5,  # Grid spacing in the y-direction (degrees)
    nx=15,  # Number of grid points in the x-direction
    ny=25,  # Number of grid points in the y-direction
)

# Each grid now has a bbox and a plot method
print(grid.bbox())  # Get the bounding box of the grid

from rompy.core.grid import RegularGrid grid = RegularGrid( x0=110.0, # Longitude of the lower-left corner y0=-35.2, # Latitude of the lower-left corner rot=4.0, # Rotation angle in degrees dx=0.5, # Grid spacing in the x-direction (degrees) dy=0.5, # Grid spacing in the y-direction (degrees) nx=15, # Number of grid points in the x-direction ny=25, # Number of grid points in the y-direction ) # Each grid now has a bbox and a plot method print(grid.bbox()) # Get the bounding box of the grid

[np.float64(109.1629223150705), np.float64(-35.2), np.float64(116.98294835181876), np.float64(-22.740936080673237)]

In [4]:

grid.plot()

grid.plot()

Out[4]:

(<Figure size 633.54x1000 with 1 Axes>, <GeoAxes: >)

In [5]:

plot_nc(gebco, "elevation", cmap="terrain")

plot_nc(gebco, "elevation", cmap="terrain")

Out[5]:

<cartopy.mpl.geocollection.GeoQuadMesh at 0x151c7eb95be0>

In [6]:

from rompy.core.data import DataGrid
from rompy.core.source import SourceFile

gebco_grid = DataGrid(
    id="gebco_bathymetry",
    source=SourceFile(uri=gebco),
    variables=["elevation"],
    coords={"y": "lat", "x": "lon"},
)
output = gebco_grid.get(grid=grid, destdir=workdir)
print_new_contents(workdir, old_contents=contents)

from rompy.core.data import DataGrid from rompy.core.source import SourceFile gebco_grid = DataGrid( id="gebco_bathymetry", source=SourceFile(uri=gebco), variables=["elevation"], coords={"y": "lat", "x": "lon"}, ) output = gebco_grid.get(grid=grid, destdir=workdir) print_new_contents(workdir, old_contents=contents)

Contents of /home/tdurrant/source/rompy/rompy-meta/repos/rompy-notebooks/notebooks/common/basedemo:
 - gebco_bathymetry.nc
 - era5_wind.nc
 - ww3_boundary.nc

Out[6]:

{'era5_wind.nc', 'gebco_bathymetry.nc', 'ww3_boundary.nc'}

In [7]:

plot_nc(output, "elevation", cmap="terrain")

plot_nc(output, "elevation", cmap="terrain")

Out[7]:

<cartopy.mpl.geocollection.GeoQuadMesh at 0x151c78e5fb30>

So what we have done here is:

1. Created a RegularGrid object defining our area of interest
2. Created a DataGrid object pointing to the GEBCO bathymetry file
3. Used the DataGrid's get method to extract the bathymetry for our grid
4. The extracted data is saved to a new NetCDF file in our working directory

For forcing data, which a time component, we can use a TimeRange object to define the time period of interest. Here we will create a TimeRange of 1 day with hourly intervals, and use it to extract wind data from the ERA5 file. By passite a time and grid to the DataGrid's get method, we can extract the data for our grid and time range.

Lets first inspect the sample ERA5 data file

In [8]:

plot_nc(era5, "u10", cmap="viridis", timestep=0)
plot_nc(era5, "v10", cmap="viridis", timestep=0)

plot_nc(era5, "u10", cmap="viridis", timestep=0) plot_nc(era5, "v10", cmap="viridis", timestep=0)

Out[8]:

<cartopy.mpl.geocollection.GeoQuadMesh at 0x151c782f2900>

Now lets create a datagrid object and use the Timerange and grid to extract the wind data. Note also that ERA5 has Latitude in reverse order, use filter to reverse it

In [9]:

from rompy.core.filters import Filter
from rompy.core.time import TimeRange

timerange = TimeRange(
    start="2023-01-01T00:00:00",
    end="2023-01-02T00:00:00",
    interval="1H",
)

era5_forcig = DataGrid(
    id="era5_wind",
    source=SourceFile(uri=era5),
    variables=["u10", "v10"],
    coords={"t": "time", "y": "latitude", "x": "longitude"},
    filter=Filter(sort=dict(coords=["latitude"])),
)
output = era5_forcig.get(grid=grid, time=timerange, destdir=workdir)
print_new_contents(workdir, old_contents=contents)

from rompy.core.filters import Filter from rompy.core.time import TimeRange timerange = TimeRange( start="2023-01-01T00:00:00", end="2023-01-02T00:00:00", interval="1H", ) era5_forcig = DataGrid( id="era5_wind", source=SourceFile(uri=era5), variables=["u10", "v10"], coords={"t": "time", "y": "latitude", "x": "longitude"}, filter=Filter(sort=dict(coords=["latitude"])), ) output = era5_forcig.get(grid=grid, time=timerange, destdir=workdir) print_new_contents(workdir, old_contents=contents)

Contents of /home/tdurrant/source/rompy/rompy-meta/repos/rompy-notebooks/notebooks/common/basedemo:
 - gebco_bathymetry.nc
 - era5_wind.nc
 - ww3_boundary.nc

Out[9]:

{'era5_wind.nc', 'gebco_bathymetry.nc', 'ww3_boundary.nc'}

Now lets have a look at the extracted wind data

In [10]:

# Plot u10 at the first timestep
plot_nc(output, "u10", cmap="viridis", timestep=0)
# Plot v10 at the first timestep
plot_nc(output, "v10", cmap="viridis", timestep=0)

# Plot u10 at the first timestep plot_nc(output, "u10", cmap="viridis", timestep=0) # Plot v10 at the first timestep plot_nc(output, "v10", cmap="viridis", timestep=0)

Out[10]:

<cartopy.mpl.geocollection.GeoQuadMesh at 0x151c7823eb70>

We can also extract wave spectra data in a similar way. Here we will extract wave spectra data from the sample WW3 spectra file. Spectra is slightly different in that we are extracting and interpolating points to specific locations along the grid boundary, rather than simply extracting a hyperslab from a gridded dataset. For this we use the BoundaryDataGrid object.

First, lets inspect the sample WW3 spectra data file

In [11]:

plot_spectra(spectra, grid=grid)

plot_spectra(spectra, grid=grid)

Now we can create a BoundaryWaveStation object and use it to extract the wave spectra data along the grid boundary.

In [12]:

from rompy.core.boundary import BoundaryWaveStation

boundary_stations = BoundaryWaveStation(
    id="ww3_boundary",
    source=SourceFile(uri=spectra),
    sel_method="idw",
    sel_method_kwargs={
        "tolerance": 4
    },  # points are sparse around the offshore boundary so make sure tolerance is appropriate
)
output = boundary_stations.get(grid=grid, time=timerange, destdir=workdir)
print_new_contents(workdir, old_contents=contents)

plot_spectra(output, grid=grid, markersize=20)


plt.show()

from rompy.core.boundary import BoundaryWaveStation boundary_stations = BoundaryWaveStation( id="ww3_boundary", source=SourceFile(uri=spectra), sel_method="idw", sel_method_kwargs={ "tolerance": 4 }, # points are sparse around the offshore boundary so make sure tolerance is appropriate ) output = boundary_stations.get(grid=grid, time=timerange, destdir=workdir) print_new_contents(workdir, old_contents=contents) plot_spectra(output, grid=grid, markersize=20) plt.show()

Contents of /home/tdurrant/source/rompy/rompy-meta/repos/rompy-notebooks/notebooks/common/basedemo:
 - gebco_bathymetry.nc
 - era5_wind.nc
 - ww3_boundary.nc

So what we have done here is:

1. Taken an oddly spaced spectral output dataset
2. Used the BoundaryWaveStations's get method interpolated those spectra to the boundary using (in this case), inverse distance weighting for the period of our run
3. The extracted and interpoated data is saved to a new NetCDF file in our working directory, ready to be ingested