Bathymetry Tutorial¶

This tutorial covers the XBeachBathy class in rompy-xbeach, which handles bathymetry interpolation and file generation for XBeach simulations.

What You'll Learn¶

1. Loading bathymetry data from various sources
2. Defining a grid aligned with the bathymetry
3. Interpolating bathymetry onto the XBeach grid
4. Extending the domain offshore and laterally
5. Visualising and validating the interpolated bathymetry

Prerequisites¶

• Familiarity with the RegularGrid class (see Grid Tutorial)
• Understanding of coordinate reference systems
• Basic knowledge of data sources (see Source Tutorial for details on SourceGeotiff and other source types)

In [1]:

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

from rompy_xbeach.grid import GeoPoint, RegularGrid
from rompy_xbeach.source import SourceGeotiff
from rompy_xbeach.data import XBeachBathy, SeawardExtensionLinear
from rompy_xbeach.interpolate import RegularGridInterpolator

import warnings
warnings.filterwarnings("ignore")

from pathlib import Path import xarray as xr import rioxarray import matplotlib.pyplot as plt import cartopy.crs as ccrs from rompy_xbeach.grid import GeoPoint, RegularGrid from rompy_xbeach.source import SourceGeotiff from rompy_xbeach.data import XBeachBathy, SeawardExtensionLinear from rompy_xbeach.interpolate import RegularGridInterpolator import warnings warnings.filterwarnings("ignore")

---

1. Loading Source Bathymetry¶

First, let's load and inspect the source bathymetry data. This example uses a GeoTIFF file, but rompy-xbeach supports various data sources.

In [2]:

from pathlib import Path

# Change the data directory to the location in your computer where rompy-xbeach is cloned into
DATADIR = Path("../../../../rompy-xbeach/tests/data")

from pathlib import Path # Change the data directory to the location in your computer where rompy-xbeach is cloned into DATADIR = Path("../../../../rompy-xbeach/tests/data")

In [3]:

# Load bathymetry using rioxarray for inspection
bathy = rioxarray.open_rasterio(DATADIR / "bathy.tif")
bathy = bathy.isel(band=0, drop=True)

print(f"Shape: {bathy.shape}")
print(f"CRS: {bathy.rio.crs}")
print(f"Bounds: x=[{float(bathy.x.min()):.4f}, {float(bathy.x.max()):.4f}], "
      f"y=[{float(bathy.y.min()):.4f}, {float(bathy.y.max()):.4f}]")

# Load bathymetry using rioxarray for inspection bathy = rioxarray.open_rasterio(DATADIR / "bathy.tif") bathy = bathy.isel(band=0, drop=True) print(f"Shape: {bathy.shape}") print(f"CRS: {bathy.rio.crs}") print(f"Bounds: x=[{float(bathy.x.min()):.4f}, {float(bathy.x.max()):.4f}], " f"y=[{float(bathy.y.min()):.4f}, {float(bathy.y.max()):.4f}]")

Shape: (180, 176)
CRS: EPSG:4326
Bounds: x=[115.5915, 115.6289], y=[-32.6483, -32.6101]

In [4]:

# Visualise the source bathymetry
bathy.plot(cmap="terrain", vmin=-20, vmax=20)
plt.title("Source Bathymetry (WGS84)")

# Visualise the source bathymetry bathy.plot(cmap="terrain", vmin=-20, vmax=20) plt.title("Source Bathymetry (WGS84)")

Out[4]:

Text(0.5, 1.0, 'Source Bathymetry (WGS84)')

---

2. Defining the XBeach Grid¶

The grid must be oriented so that:

• The offshore boundary (front) faces the open ocean
• The origin is at the onshore corner
• The x-axis points offshore (cross-shore direction)

The plot() method highlights the origin (red dot) and offshore boundary (red line) to help verify correct orientation.

In [5]:

grid = RegularGrid(
    ori=GeoPoint(
        x=115.594239,
        y=-32.641104,
        crs=bathy.rio.crs,  # Origin in same CRS as bathymetry
    ),
    alfa=347,    # Rotation angle (degrees from East)
    dx=10,       # Cross-shore resolution (m)
    dy=15,       # Alongshore resolution (m)
    nx=230,      # Cross-shore cells
    ny=220,      # Alongshore cells
    crs="28350", # Grid CRS (MGA Zone 50)
)

grid = RegularGrid( ori=GeoPoint( x=115.594239, y=-32.641104, crs=bathy.rio.crs, # Origin in same CRS as bathymetry ), alfa=347, # Rotation angle (degrees from East) dx=10, # Cross-shore resolution (m) dy=15, # Alongshore resolution (m) nx=230, # Cross-shore cells ny=220, # Alongshore cells crs="28350", # Grid CRS (MGA Zone 50) )

Visualising Grid on Bathymetry¶

Overlay the grid on the bathymetry to verify alignment:

In [6]:

fig, ax = plt.subplots(subplot_kw={"projection": grid.projection}, figsize=(8.5, 8))

# Plot bathymetry
bathy.plot(
    ax=ax,
    transform=ccrs.PlateCarree(),
    cmap="terrain",
    vmin=-20, vmax=20,
    cbar_kwargs=dict(label="Depth (m)"),
)

# Overlay grid with origin and offshore boundary highlighted
grid.plot(
    ax=ax,
    buffer=250,
    grid_kwargs=dict(facecolor="k", alpha=0.2),
    show_offshore=True,
    show_origin=True,
)

ax.set_title("Grid aligned with bathymetry\n(Red dot = origin, Red line = offshore)")

fig, ax = plt.subplots(subplot_kw={"projection": grid.projection}, figsize=(8.5, 8)) # Plot bathymetry bathy.plot( ax=ax, transform=ccrs.PlateCarree(), cmap="terrain", vmin=-20, vmax=20, cbar_kwargs=dict(label="Depth (m)"), ) # Overlay grid with origin and offshore boundary highlighted grid.plot( ax=ax, buffer=250, grid_kwargs=dict(facecolor="k", alpha=0.2), show_offshore=True, show_origin=True, ) ax.set_title("Grid aligned with bathymetry\n(Red dot = origin, Red line = offshore)")

Out[6]:

Text(0.5, 1.0, 'Grid aligned with bathymetry\n(Red dot = origin, Red line = offshore)')

Common Mistake: Incorrect Grid Orientation¶

If the origin or rotation is wrong, the offshore boundary will face the wrong direction. Here's an example of a poorly-defined grid:

In [7]:

# Define origin at the NE corner (wrong!) with flipped rotation
grid_bad = RegularGrid(
    ori=GeoPoint(x=grid.x[-1, -1], y=grid.y[-1, -1], crs=28350),
    alfa=347 - 180,  # Flipped rotation
    dx=10, dy=15, nx=230, ny=220,
    crs="28350",
)

fig, ax = plt.subplots(subplot_kw={"projection": grid.projection}, figsize=(8.5, 8))
bathy.plot(ax=ax, transform=ccrs.PlateCarree(), cmap="terrain", vmin=-20, vmax=20)
grid_bad.plot(ax=ax, buffer=250, grid_kwargs=dict(facecolor="r", alpha=0.2))
ax.set_title("Incorrectly oriented grid\n(Offshore boundary faces land!)")

# Define origin at the NE corner (wrong!) with flipped rotation grid_bad = RegularGrid( ori=GeoPoint(x=grid.x[-1, -1], y=grid.y[-1, -1], crs=28350), alfa=347 - 180, # Flipped rotation dx=10, dy=15, nx=230, ny=220, crs="28350", ) fig, ax = plt.subplots(subplot_kw={"projection": grid.projection}, figsize=(8.5, 8)) bathy.plot(ax=ax, transform=ccrs.PlateCarree(), cmap="terrain", vmin=-20, vmax=20) grid_bad.plot(ax=ax, buffer=250, grid_kwargs=dict(facecolor="r", alpha=0.2)) ax.set_title("Incorrectly oriented grid\n(Offshore boundary faces land!)")

Out[7]:

Text(0.5, 1.0, 'Incorrectly oriented grid\n(Offshore boundary faces land!)')

---

3. Interpolating Bathymetry¶

The XBeachBathy class handles interpolation from the source data onto the XBeach grid. It requires:

• A source object pointing to the bathymetry data
• An interpolator specifying the interpolation method
• Optional extension for offshore boundary extension

Basic Interpolation (No Extension)¶

In [8]:

# Define the data source
source = SourceGeotiff(filename=DATADIR / "bathy.tif")

# Define the interpolation method
interpolator = RegularGridInterpolator(
    kwargs=dict(
        method="linear",    # Interpolation method
        fill_value=None,    # Extrapolate at boundaries
    ),
)

# Create the bathymetry data object
data = XBeachBathy(
    source=source,
    posdwn=False,           # Negative = deeper (elevation convention)
    interpolator=interpolator,
)

# Define the data source source = SourceGeotiff(filename=DATADIR / "bathy.tif") # Define the interpolation method interpolator = RegularGridInterpolator( kwargs=dict( method="linear", # Interpolation method fill_value=None, # Extrapolate at boundaries ), ) # Create the bathymetry data object data = XBeachBathy( source=source, posdwn=False, # Negative = deeper (elevation convention) interpolator=interpolator, )

In [9]:

# Generate the XBeach files
workspace = Path("./examples")
workspace.mkdir(exist_ok=True)

xfile, yfile, datafile, grid_out = data.get(destdir=workspace, grid=grid)

print("Generated files:")
for f in sorted(workspace.glob("*.txt")):
    print(f"  {f.name}")

# Generate the XBeach files workspace = Path("./examples") workspace.mkdir(exist_ok=True) xfile, yfile, datafile, grid_out = data.get(destdir=workspace, grid=grid) print("Generated files:") for f in sorted(workspace.glob("*.txt")): print(f" {f.name}")

Generated files:
  bathy.txt
  xdata.txt
  ydata.txt

Reading and Visualising the Result¶

The xbeach accessor provides methods to read and plot XBeach bathymetry files:

In [10]:

# Read the interpolated bathymetry
dset = xr.Dataset.xbeach.from_xbeach(datafile, grid_out)
dset

# Read the interpolated bathymetry dset = xr.Dataset.xbeach.from_xbeach(datafile, grid_out) dset

Out[10]:

<xarray.Dataset> Size: 1MB
Dimensions:      (y: 220, x: 230)
Coordinates:
    xc           (y, x) float64 405kB 3.681e+05 3.682e+05 ... 3.711e+05
    yc           (y, x) float64 405kB 6.388e+06 6.388e+06 ... 6.39e+06 6.39e+06
    spatial_ref  int64 8B 0
Dimensions without coordinates: y, x
Data variables:
    dep          (y, x) float64 405kB -13.1 -13.08 -13.08 ... 14.63 14.45 14.53

xarray.Dataset

• Dimensions:
  • y: 220
  • x: 230
• Coordinates: (3)
  • xc
    (y, x)
    float64
    3.681e+05 3.682e+05 ... 3.711e+05

    array([[368145.00235004, 368154.74605069, 368164.48975134, ...,
            370356.82239711, 370366.56609775, 370376.3097984 ],
           [368148.37661586, 368158.12031651, 368167.86401716, ...,
            370360.19666292, 370369.94036357, 370379.68406422],
           [368151.75088167, 368161.49458232, 368171.23828297, ...,
            370363.57092874, 370373.31462939, 370383.05833003],
           ...,
           [368877.21803193, 368886.96173258, 368896.70543323, ...,
            371089.038079  , 371098.78177964, 371108.52548029],
           [368880.59229775, 368890.3359984 , 368900.07969904, ...,
            371092.41234481, 371102.15604546, 371111.89974611],
           [368883.96656356, 368893.71026421, 368903.45396486, ...,
            371095.78661063, 371105.53031127, 371115.27401192]],
          shape=(220, 230))

  • yc
    (y, x)
    float64
    6.388e+06 6.388e+06 ... 6.39e+06

    array([[6387626.58926688, 6387624.33975634, 6387622.0902458 , ...,
            6387115.95037352, 6387113.70086298, 6387111.45135244],
           [6387641.20481785, 6387638.95530731, 6387636.70579677, ...,
            6387130.56592449, 6387128.31641395, 6387126.06690341],
           [6387655.82036883, 6387653.57085828, 6387651.32134774, ...,
            6387145.18147547, 6387142.93196492, 6387140.68245438],
           ...,
           [6390798.16382776, 6390795.91431722, 6390793.66480667, ...,
            6390287.5249344 , 6390285.27542386, 6390283.02591331],
           [6390812.77937873, 6390810.52986819, 6390808.28035764, ...,
            6390302.14048537, 6390299.89097483, 6390297.64146428],
           [6390827.3949297 , 6390825.14541916, 6390822.89590862, ...,
            6390316.75603634, 6390314.5065258 , 6390312.25701526]],
          shape=(220, 230))

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    PROJCS["GDA94 / MGA zone 50",GEOGCS["GDA94",DATUM["Geocentric_Datum_of_Australia_1994",SPHEROID["GRS 1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4283"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","28350"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314140356

    inverse_flattening :

    298.257222101

    reference_ellipsoid_name :

    GRS 1980

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    GDA94

    horizontal_datum_name :

    Geocentric Datum of Australia 1994

    projected_crs_name :

    GDA94 / MGA zone 50

    grid_mapping_name :

    transverse_mercator

    latitude_of_projection_origin :

    0.0

    longitude_of_central_meridian :

    117.0

    false_easting :

    500000.0

    false_northing :

    10000000.0

    scale_factor_at_central_meridian :

    0.9996

    spatial_ref :

    PROJCS["GDA94 / MGA zone 50",GEOGCS["GDA94",DATUM["Geocentric_Datum_of_Australia_1994",SPHEROID["GRS 1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4283"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","28350"]]

    array(0)

• Data variables: (1)
  • dep
    (y, x)
    float64
    -13.1 -13.08 -13.08 ... 14.45 14.53

    array([[-13.10233341, -13.08240353, -13.07586602, ...,  11.43840582,
             14.89170182,  17.03182526],
           [-13.17888252, -13.16692934, -13.06563562, ...,  11.90206947,
             15.78203909,  16.26002547],
           [-13.13399074, -13.10684213, -13.05971619, ...,  12.24710476,
             15.27183094,  15.18587995],
           ...,
           [-13.73772084, -13.58594975, -13.42458863, ...,  13.80981498,
             15.52806884,  15.93475685],
           [-13.66211277, -13.51602798, -13.37380363, ...,  12.05110652,
             12.36159615,  13.02389086],
           [-13.55875508, -13.38831277, -13.27075777, ...,  14.62773471,
             14.44870578,  14.52679177]], shape=(220, 230))

In [11]:

# Plot the model bathymetry
dset.xbeach.plot_model_bathy(grid_out, posdwn=False)

# Plot the model bathymetry dset.xbeach.plot_model_bathy(grid_out, posdwn=False)

---

4. Offshore Extension¶

Often the source bathymetry doesn't extend far enough offshore for the wave boundary with homogeneous depths. The SeawardExtensionLinear class extends the domain with a linear slope to a specified depth.

In [12]:

# Define offshore extension
extension = SeawardExtensionLinear(
    depth=25,    # Target depth at offshore boundary (m)
    slope=0.1,   # Slope of the extension (rise/run)
)

# Create bathymetry with extension
data_extended = XBeachBathy(
    source=source,
    posdwn=False,
    interpolator=interpolator,
    extension=extension,
)

# Generate files
xfile2, yfile2, datafile2, grid2 = data_extended.get(destdir=workspace, grid=grid)

# Define offshore extension extension = SeawardExtensionLinear( depth=25, # Target depth at offshore boundary (m) slope=0.1, # Slope of the extension (rise/run) ) # Create bathymetry with extension data_extended = XBeachBathy( source=source, posdwn=False, interpolator=interpolator, extension=extension, ) # Generate files xfile2, yfile2, datafile2, grid2 = data_extended.get(destdir=workspace, grid=grid)

In [13]:

# Visualise the extended bathymetry
dset2 = xr.Dataset.xbeach.from_xbeach(datafile2, grid2)
dset2.xbeach.plot_model_bathy(grid2, posdwn=False)
plt.title("Bathymetry with Offshore Extension (slope=0.1)")

# Visualise the extended bathymetry dset2 = xr.Dataset.xbeach.from_xbeach(datafile2, grid2) dset2.xbeach.plot_model_bathy(grid2, posdwn=False) plt.title("Bathymetry with Offshore Extension (slope=0.1)")

Out[13]:

Text(0.5, 1.0, 'Bathymetry with Offshore Extension (slope=0.1)')

In [14]:

# Compare grid dimensions
print(f"Original grid:  {grid_out.shape}")
print(f"Extended grid:  {grid2.shape}")
print(f"Extension added {grid2.nx - grid_out.nx} cells offshore")

# Compare grid dimensions print(f"Original grid: {grid_out.shape}") print(f"Extended grid: {grid2.shape}") print(f"Extension added {grid2.nx - grid_out.nx} cells offshore")

Original grid:  (220, 230)
Extended grid:  (220, 242)
Extension added 12 cells offshore

Effect of Slope on Extension Length¶

A gentler slope results in a longer extension to reach the target depth:

In [15]:

# Gentler slope = longer extension
extension_gentle = SeawardExtensionLinear(depth=25, slope=0.02)

data_gentle = XBeachBathy(
    source=source,
    posdwn=False,
    interpolator=interpolator,
    extension=extension_gentle,
)

xfile3, yfile3, datafile3, grid3 = data_gentle.get(destdir=workspace, grid=grid)

dset3 = xr.Dataset.xbeach.from_xbeach(datafile3, grid3)
dset3.xbeach.plot_model_bathy(grid3, posdwn=False)
plt.title("Bathymetry with Gentle Slope Extension (slope=0.02)")

# Gentler slope = longer extension extension_gentle = SeawardExtensionLinear(depth=25, slope=0.02) data_gentle = XBeachBathy( source=source, posdwn=False, interpolator=interpolator, extension=extension_gentle, ) xfile3, yfile3, datafile3, grid3 = data_gentle.get(destdir=workspace, grid=grid) dset3 = xr.Dataset.xbeach.from_xbeach(datafile3, grid3) dset3.xbeach.plot_model_bathy(grid3, posdwn=False) plt.title("Bathymetry with Gentle Slope Extension (slope=0.02)")

Out[15]:

Text(0.5, 1.0, 'Bathymetry with Gentle Slope Extension (slope=0.02)')

In [16]:

print(f"Steep slope (0.1):   {grid2.nx - grid_out.nx} cells added")
print(f"Gentle slope (0.02): {grid3.nx - grid_out.nx} cells added")

print(f"Steep slope (0.1): {grid2.nx - grid_out.nx} cells added") print(f"Gentle slope (0.02): {grid3.nx - grid_out.nx} cells added")

Steep slope (0.1):   12 cells added
Gentle slope (0.02): 60 cells added

---

5. Lateral Extension¶

The domain can also be extended laterally (left and right) to provide buffer zones for boundary conditions. This copies the edge profiles outward.

In [17]:

# Extend offshore AND laterally
data_full = XBeachBathy(
    source=source,
    posdwn=False,
    interpolator=interpolator,
    extension=SeawardExtensionLinear(depth=25, slope=0.02),
    left=10,    # Extend left boundary by 10 cells
    right=10,   # Extend right boundary by 10 cells
)

xfile4, yfile4, datafile4, grid4 = data_full.get(destdir=workspace, grid=grid)

dset4 = xr.Dataset.xbeach.from_xbeach(datafile4, grid4)
dset4.xbeach.plot_model_bathy(grid4, posdwn=False)
plt.title("Bathymetry with Offshore + Lateral Extension")

# Extend offshore AND laterally data_full = XBeachBathy( source=source, posdwn=False, interpolator=interpolator, extension=SeawardExtensionLinear(depth=25, slope=0.02), left=10, # Extend left boundary by 10 cells right=10, # Extend right boundary by 10 cells ) xfile4, yfile4, datafile4, grid4 = data_full.get(destdir=workspace, grid=grid) dset4 = xr.Dataset.xbeach.from_xbeach(datafile4, grid4) dset4.xbeach.plot_model_bathy(grid4, posdwn=False) plt.title("Bathymetry with Offshore + Lateral Extension")

Out[17]:

Text(0.5, 1.0, 'Bathymetry with Offshore + Lateral Extension')

In [18]:

# Compare dimensions
print(f"No extension:        {grid_out.shape}")
print(f"Offshore only:       {grid3.shape}")
print(f"Offshore + lateral:  {grid4.shape}")

# Compare dimensions print(f"No extension: {grid_out.shape}") print(f"Offshore only: {grid3.shape}") print(f"Offshore + lateral: {grid4.shape}")

No extension:        (220, 230)
Offshore only:       (220, 290)
Offshore + lateral:  (240, 290)

---

6. Summary¶

Key Classes¶

Class	Purpose
SourceGeotiff	Load bathymetry from GeoTIFF files
RegularGridInterpolator	Configure interpolation method
XBeachBathy	Interpolate and generate XBeach bathymetry files
SeawardExtensionLinear	Extend domain offshore with linear slope

Key Parameters¶

Parameter	Description
source	Data source object
interpolator	Interpolation configuration
posdwn	Depth convention (True = positive down)
extension	Offshore extension configuration
left, right	Lateral extension (number of cells)

Workflow¶

# 1. Define source and interpolator
source = SourceGeotiff(filename="bathy.tif")
interpolator = RegularGridInterpolator(kwargs=dict(method="linear"))

# 2. Create bathymetry object with optional extensions
bathy = XBeachBathy(
    source=source,
    interpolator=interpolator,
    extension=SeawardExtensionLinear(depth=25, slope=0.05),
    left=10, right=10,
)

# 3. Generate files
xfile, yfile, datafile, grid = bathy.get(destdir=destdir, grid=grid)

Next Steps¶

• See the Source tutorial for other data source types
• See the Wave Boundary tutorial for setting up wave forcing