InĀ [1]:
Copied!
%load_ext autoreload
%autoreload 2
%load_ext autoreload
%autoreload 2
InĀ [2]:
Copied!
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
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
import warnings
warnings.filterwarnings("ignore")
Read bathyĀ¶
Read bathy from the test data in the rompy-xbeach repo to inspect it and help it to define the model grid
InĀ [3]:
Copied!
bathy = rioxarray.open_rasterio("../../../rompy-xbeach/tests/data/bathy.tif")
bathy = bathy.isel(band=0, drop=True)
bathy
bathy = rioxarray.open_rasterio("../../../rompy-xbeach/tests/data/bathy.tif")
bathy = bathy.isel(band=0, drop=True)
bathy
Out[3]:
<xarray.DataArray (y: 180, x: 176)> Size: 127kB
[31680 values with dtype=float32]
Coordinates:
* x (x) float64 1kB 115.6 115.6 115.6 115.6 ... 115.6 115.6 115.6
* y (y) float64 1kB -32.65 -32.65 -32.65 ... -32.61 -32.61 -32.61
spatial_ref int64 8B 0
Attributes:
AREA_OR_POINT: Area
scale_factor: 1.0
add_offset: 0.0InĀ [4]:
Copied!
# The bathy is in lat/lon
bathy.plot()
# The bathy is in lat/lon
bathy.plot()
Out[4]:
<matplotlib.collections.QuadMesh at 0x7c65cbe71340>
Set the gridĀ¶
Define the XBeach grid object in EPSG:28350 projection (use latlon to specify the origin)
InĀ [5]:
Copied!
from pyproj import CRS
from pyproj import CRS
InĀ [6]:
Copied!
grid = RegularGrid(
ori=GeoPoint(
x=115.594239,
y=-32.641104,
crs=bathy.rio.crs,
),
alfa=347,
dx=10,
dy=15,
nx=230,
ny=220,
crs="28350",
)
grid = RegularGrid(
ori=GeoPoint(
x=115.594239,
y=-32.641104,
crs=bathy.rio.crs,
),
alfa=347,
dx=10,
dy=15,
nx=230,
ny=220,
crs="28350",
)
Plot grid on bathyĀ¶
The origin and rotation of the grid need to be defined such that the offshore boundary
of the interpolated bathymetry data corresponds to the first row and column. The grid.plot()
method highlights the origin and offshore boundary in red to assist with the grid definition
InĀ [7]:
Copied!
fig, ax = plt.subplots(subplot_kw={"projection": grid.projection}, figsize=(8.5, 8))
# Bathymetry
bathy.plot(
ax=ax,
transform=ccrs.PlateCarree(),
cmap="terrain",
vmin=-20,
vmax=20,
cbar_kwargs=dict(label="Depth (m)"),
)
# Grid
ax = grid.plot(
ax=ax,
buffer=250,
grid_kwargs=dict(facecolor="k", alpha=0.2),
show_offshore=True,
show_origin=True,
)
fig, ax = plt.subplots(subplot_kw={"projection": grid.projection}, figsize=(8.5, 8))
# Bathymetry
bathy.plot(
ax=ax,
transform=ccrs.PlateCarree(),
cmap="terrain",
vmin=-20,
vmax=20,
cbar_kwargs=dict(label="Depth (m)"),
)
# Grid
ax = grid.plot(
ax=ax,
buffer=250,
grid_kwargs=dict(facecolor="k", alpha=0.2),
show_offshore=True,
show_origin=True,
)
ā ļø Here we will define another grid by poorly selecting the origin and rotation
InĀ [8]:
Copied!
# Define origin at the NE corner of the existing grid
x0 = grid.x[-1, -1]
y0 = grid.y[-1, -1]
# Notice we are now specifying the origin in projected coordinates!
grid_bad = RegularGrid(
ori=GeoPoint(x=x0, y=y0, crs=28350),
alfa=347-180,
dx=10,
dy=15,
nx=230,
ny=220,
crs="28350",
)
# Plotting
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,
cbar_kwargs=dict(label="Depth (m)"),
)
ax = grid_bad.plot(ax=ax, buffer=250, grid_kwargs=dict(facecolor="k", alpha=0.2))
# Define origin at the NE corner of the existing grid
x0 = grid.x[-1, -1]
y0 = grid.y[-1, -1]
# Notice we are now specifying the origin in projected coordinates!
grid_bad = RegularGrid(
ori=GeoPoint(x=x0, y=y0, crs=28350),
alfa=347-180,
dx=10,
dy=15,
nx=230,
ny=220,
crs="28350",
)
# Plotting
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,
cbar_kwargs=dict(label="Depth (m)"),
)
ax = grid_bad.plot(ax=ax, buffer=250, grid_kwargs=dict(facecolor="k", alpha=0.2))
Set Data objectĀ¶
InĀ [9]:
Copied!
from rompy_xbeach.source import SourceGeotiff
from rompy_xbeach.data import XBeachBathy, SeawardExtensionLinear
from rompy_xbeach.interpolate import RegularGridInterpolator
from rompy_xbeach.source import SourceGeotiff
from rompy_xbeach.data import XBeachBathy, SeawardExtensionLinear
from rompy_xbeach.interpolate import RegularGridInterpolator
Example 1: no extension of the gridĀ¶
We start off by defining source and interpolator instances
InĀ [10]:
Copied!
# The source instance provides the source dataset with the CRS information
source = SourceGeotiff(filename="../../../rompy-xbeach/tests/data/bathy.tif")
# The interpolator instance specify how the interpolation will be done
interpolator = RegularGridInterpolator(
kwargs=dict(
method="linear",
fill_value=None,
),
)
# The source instance provides the source dataset with the CRS information
source = SourceGeotiff(filename="../../../rompy-xbeach/tests/data/bathy.tif")
# The interpolator instance specify how the interpolation will be done
interpolator = RegularGridInterpolator(
kwargs=dict(
method="linear",
fill_value=None,
),
)
Then we instantiate the bathy data object
InĀ [11]:
Copied!
data = XBeachBathy(
source=source,
posdwn=False,
interpolator=interpolator,
)
data = XBeachBathy(
source=source,
posdwn=False,
interpolator=interpolator,
)
Lastly, we use the previously defined grid instance to interpolate the bathy onto
InĀ [12]:
Copied!
workspace = Path("./examples")
workspace.mkdir(exist_ok=True)
xfile1, yfile1, datafile1, grid1 = data.get(destdir=workspace, grid=grid)
sorted(workspace.glob("*.txt"))
workspace = Path("./examples")
workspace.mkdir(exist_ok=True)
xfile1, yfile1, datafile1, grid1 = data.get(destdir=workspace, grid=grid)
sorted(workspace.glob("*.txt"))
Out[12]:
[PosixPath('examples/bathy.txt'),
PosixPath('examples/xdata.txt'),
PosixPath('examples/ydata.txt')]
The xbeach accessor can be used to read the output data into a xarray Dataset and plot the bathy
InĀ [13]:
Copied!
# Read the data
dset1 = xr.Dataset.xbeach.from_xbeach(datafile1, grid1)
dset1
# Read the data
dset1 = xr.Dataset.xbeach.from_xbeach(datafile1, grid1)
dset1
Out[13]:
<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.53InĀ [14]:
Copied!
# Plot the data
dset1.xbeach.plot_model_bathy(grid1, posdwn=False)
# Plot the data
dset1.xbeach.plot_model_bathy(grid1, posdwn=False)
Example 2: Extend offshore boundaryĀ¶
InĀ [15]:
Copied!
# The offshore extension is prescribed using one of the SeawardExtension subclasses
extension = SeawardExtensionLinear(
depth=25,
slope=0.1,
)
# Instantiate the data object with the new extension
data = XBeachBathy(
source=source,
posdwn=False,
interpolator=interpolator,
extension=extension,
)
# Generate the grid data
xfile2, yfile2, datafile2, grid2 = data.get(destdir=workspace, grid=grid)
# Plot the new bathy
dset2 = xr.Dataset.xbeach.from_xbeach(datafile2, grid2)
dset2.xbeach.plot_model_bathy(grid2, posdwn=False)
# The offshore extension is prescribed using one of the SeawardExtension subclasses
extension = SeawardExtensionLinear(
depth=25,
slope=0.1,
)
# Instantiate the data object with the new extension
data = XBeachBathy(
source=source,
posdwn=False,
interpolator=interpolator,
extension=extension,
)
# Generate the grid data
xfile2, yfile2, datafile2, grid2 = data.get(destdir=workspace, grid=grid)
# Plot the new bathy
dset2 = xr.Dataset.xbeach.from_xbeach(datafile2, grid2)
dset2.xbeach.plot_model_bathy(grid2, posdwn=False)
InĀ [16]:
Copied!
# Compare the shapes of the original and extended grid
print(dset1.dep.shape, grid1.shape)
print(dset2.dep.shape, grid2.shape)
# Compare the shapes of the original and extended grid
print(dset1.dep.shape, grid1.shape)
print(dset2.dep.shape, grid2.shape)
(220, 230) (220, 230) (220, 242) (220, 242)
The lesser the slope the longer the extensionĀ¶
InĀ [17]:
Copied!
extension = SeawardExtensionLinear(
depth=25,
slope=0.02,
)
data = XBeachBathy(
source=source,
posdwn=False,
interpolator=interpolator,
extension=extension,
)
xfile3, yfile3, datafile3, grid3 = data.get(destdir=workspace, grid=grid)
dset3 = xr.Dataset.xbeach.from_xbeach(datafile3, grid3)
dset3.xbeach.plot_model_bathy(grid3, posdwn=False)
extension = SeawardExtensionLinear(
depth=25,
slope=0.02,
)
data = XBeachBathy(
source=source,
posdwn=False,
interpolator=interpolator,
extension=extension,
)
xfile3, yfile3, datafile3, grid3 = data.get(destdir=workspace, grid=grid)
dset3 = xr.Dataset.xbeach.from_xbeach(datafile3, grid3)
dset3.xbeach.plot_model_bathy(grid3, posdwn=False)
Exemple 3: extend offshore and lateral boundaryĀ¶
InĀ [18]:
Copied!
extension = SeawardExtensionLinear(
depth=25,
slope=0.02,
)
data = XBeachBathy(
source=source,
posdwn=False,
interpolator=interpolator,
extension=extension,
left=10,
right=10,
)
xfile4, yfile4, datafile4, grid4 = data.get(destdir=workspace, grid=grid)
dset4 = xr.Dataset.xbeach.from_xbeach(datafile4, grid4)
dset4.xbeach.plot_model_bathy(grid4, posdwn=False)
extension = SeawardExtensionLinear(
depth=25,
slope=0.02,
)
data = XBeachBathy(
source=source,
posdwn=False,
interpolator=interpolator,
extension=extension,
left=10,
right=10,
)
xfile4, yfile4, datafile4, grid4 = data.get(destdir=workspace, grid=grid)
dset4 = xr.Dataset.xbeach.from_xbeach(datafile4, grid4)
dset4.xbeach.plot_model_bathy(grid4, posdwn=False)
InĀ [19]:
Copied!
# Compare the shapes with and without the lateral extension
print(dset3.dep.shape, grid3.shape)
print(dset4.dep.shape, grid4.shape)
# Compare the shapes with and without the lateral extension
print(dset3.dep.shape, grid3.shape)
print(dset4.dep.shape, grid4.shape)
(220, 290) (220, 290) (240, 290) (240, 290)