SWAN BOUNDSPEC SIDE¶

This notebooks shows examples of how to prescribe BOUNDSPEC SIDE spectral boundaries using rompy boundary data objects.

Boundary of type BOUNDSPEC SIDE can be specified from existing spectral data using the rompy_swan.boundary.BoundspecSide object. The boundary location is defined as the middle point of a side of a specified grid object. The boundary spectra are defined by interpolating or nearest-selecting from the source spectra data at the selected boundary location.

This class only supports FILE type boundary since the PAR type does not allow for nonstationary boundary specification (the user can prescribe PAR type boundary by using the SWAN components if required).

We currently support TPAR and SPEC2D file types to write the boundary data.

In [1]:

%load_ext autoreload
%autoreload 2

import os
from pathlib import Path
from tempfile import mkdtemp
import matplotlib.pyplot as plt
import matplotlib
from matplotlib import colors
import cartopy.crs as ccrs
import pandas as pd
import wavespectra

from rompy.core.source import SourceFile
from rompy.core.time import TimeRange
from rompy_swan.grid import SwanGrid
from rompy_swan.boundary import BoundspecSide

import warnings
warnings.filterwarnings('ignore')

%load_ext autoreload %autoreload 2 import os from pathlib import Path from tempfile import mkdtemp import matplotlib.pyplot as plt import matplotlib from matplotlib import colors import cartopy.crs as ccrs import pandas as pd import wavespectra from rompy.core.source import SourceFile from rompy.core.time import TimeRange from rompy_swan.grid import SwanGrid from rompy_swan.boundary import BoundspecSide import warnings warnings.filterwarnings('ignore')

Helper functions¶

In [2]:

def read_tpar(tparfile):
    """Read SWAN TPAR files as a Dataframe."""
    df = pd.read_csv(
        tparfile,
        sep=" ",
        header=None,
        skiprows=1,
        parse_dates=[0],
        date_format="%Y%m%d.%H%M%S",
        index_col=0,
    )
    df.columns = ["hs", "tp", "dpm", "dspr"]
    return df


def get_color(cmap, value, vmin, vmax):
    """Return the HEX color for a value in a colormap."""
    cmap = matplotlib.cm.get_cmap(cmap)
    norm = colors.Normalize(vmin=vmin, vmax=vmax, clip=True)
    return colors.rgb2hex(cmap(norm(value)))


def plot_boundary_side(dset, df, grid, time):
    """Plot side boundary stats.

    This function plots a snapshot of wave statistics from the source spectra dataset
    and the boundary wave stats generated for a grid from those spectra points.

    Parameters
    ----------
    dset: xr.Dataset
        Source spectra dataset used to define boundary from.
    df: pd.DataFrame
        Boundary time series defined from TPAR or SWAN ASCII file.
    grid: SwanGrid
        Grid object defining the SWAN domain.
    time: Datatime | str
        Time to plot.

    """
    # xb, yb = bnd._boundary_points(grid)
    xb = grid.x[:, 0]
    yb = grid.y[:, 0]

    fig, axs = plt.subplots(1, 3, subplot_kw=dict(projection=ccrs.PlateCarree()), figsize=(15, 3.5))

    # Hs
    ax = axs[0]
    p = ax.scatter(dset.lon, dset.lat, s=50, c=dset.sel(time=time).spec.hs(), cmap="turbo", vmin=2.0, vmax=3.5)
    color = get_color("turbo", df.loc[time]["hs"], 2.0, 3.5)
    ax.plot(xb, yb, color=color, linewidth=4, zorder=10)
    plt.colorbar(p, label="Hs (m)")

    # Tp
    ax = axs[1]
    p = ax.scatter(dset.lon, dset.lat, s=50, c=dset.sel(time=time).spec.tp(), cmap="viridis", vmin=13.8, vmax=14.5)
    color = get_color("viridis", df.loc[time]["tp"], 13.8, 14.5)
    ax.plot(xb, yb, color=color, linewidth=4, zorder=10)
    plt.colorbar(p, label="Tp (s)")

    # Dpm
    ax = axs[2]
    p = ax.scatter(dset.lon, dset.lat, s=50, c=dset.sel(time=time).spec.dpm(), cmap="hsv", vmin=220, vmax=225)
    color = get_color("hsv", df.loc[time]["dpm"], 220, 225)
    ax.plot(xb, yb, color=color, linewidth=4, zorder=10)
    # ax.scatter(xb, yb, s=200, c=df.loc[time]["dpm"], marker=">", edgecolor="k", cmap="hsv", vmin=220, vmax=225)
    plt.colorbar(p, label="Dpm (deg)")

    for ax in axs:
        grid.plot(ax=ax, fscale=5)
        ax.set_extent([110.5, 115.5, -29.5, -23.5])
        ax.set_title(time)
        ax.coastlines()

def read_tpar(tparfile): """Read SWAN TPAR files as a Dataframe.""" df = pd.read_csv( tparfile, sep=" ", header=None, skiprows=1, parse_dates=[0], date_format="%Y%m%d.%H%M%S", index_col=0, ) df.columns = ["hs", "tp", "dpm", "dspr"] return df def get_color(cmap, value, vmin, vmax): """Return the HEX color for a value in a colormap.""" cmap = matplotlib.cm.get_cmap(cmap) norm = colors.Normalize(vmin=vmin, vmax=vmax, clip=True) return colors.rgb2hex(cmap(norm(value))) def plot_boundary_side(dset, df, grid, time): """Plot side boundary stats. This function plots a snapshot of wave statistics from the source spectra dataset and the boundary wave stats generated for a grid from those spectra points. Parameters ---------- dset: xr.Dataset Source spectra dataset used to define boundary from. df: pd.DataFrame Boundary time series defined from TPAR or SWAN ASCII file. grid: SwanGrid Grid object defining the SWAN domain. time: Datatime | str Time to plot. """ # xb, yb = bnd._boundary_points(grid) xb = grid.x[:, 0] yb = grid.y[:, 0] fig, axs = plt.subplots(1, 3, subplot_kw=dict(projection=ccrs.PlateCarree()), figsize=(15, 3.5)) # Hs ax = axs[0] p = ax.scatter(dset.lon, dset.lat, s=50, c=dset.sel(time=time).spec.hs(), cmap="turbo", vmin=2.0, vmax=3.5) color = get_color("turbo", df.loc[time]["hs"], 2.0, 3.5) ax.plot(xb, yb, color=color, linewidth=4, zorder=10) plt.colorbar(p, label="Hs (m)") # Tp ax = axs[1] p = ax.scatter(dset.lon, dset.lat, s=50, c=dset.sel(time=time).spec.tp(), cmap="viridis", vmin=13.8, vmax=14.5) color = get_color("viridis", df.loc[time]["tp"], 13.8, 14.5) ax.plot(xb, yb, color=color, linewidth=4, zorder=10) plt.colorbar(p, label="Tp (s)") # Dpm ax = axs[2] p = ax.scatter(dset.lon, dset.lat, s=50, c=dset.sel(time=time).spec.dpm(), cmap="hsv", vmin=220, vmax=225) color = get_color("hsv", df.loc[time]["dpm"], 220, 225) ax.plot(xb, yb, color=color, linewidth=4, zorder=10) # ax.scatter(xb, yb, s=200, c=df.loc[time]["dpm"], marker=">", edgecolor="k", cmap="hsv", vmin=220, vmax=225) plt.colorbar(p, label="Dpm (deg)") for ax in axs: grid.plot(ax=ax, fscale=5) ax.set_extent([110.5, 115.5, -29.5, -23.5]) ax.set_title(time) ax.coastlines()

In [3]:

datadir = Path("../../../../rompy-swan/tests/data")
modeldir = Path("model")
modeldir.mkdir(exist_ok=True)

datadir = Path("../../../../rompy-swan/tests/data") modeldir = Path("model") modeldir.mkdir(exist_ok=True)

Define rompy objects¶

Create instances of source, time and grid objects to use with the boundary classes

In [4]:

# source defines the dataset to be used to create the model boundary

source = SourceFile(
    uri=datadir / "aus-20230101.nc",
    kwargs=dict(engine="netcdf4"),
)
dset = source.open()
dset

# source defines the dataset to be used to create the model boundary source = SourceFile( uri=datadir / "aus-20230101.nc", kwargs=dict(engine="netcdf4"), ) dset = source.open() dset

Out[4]:

<xarray.Dataset> Size: 1MB
Dimensions:  (site: 412, time: 5, freq: 11, dir: 8)
Coordinates:
  * time     (time) datetime64[ns] 40B 2023-01-01 ... 2023-01-02
  * site     (site) int64 3kB 0 4 8 12 16 20 ... 1624 1628 1632 1636 1640 1644
  * freq     (freq) float32 44B 0.05417 0.05959 0.06555 ... 0.1161 0.1277 0.1405
  * dir      (dir) float32 32B 0.0 45.0 90.0 135.0 180.0 225.0 270.0 315.0
Data variables:
    lon      (site) float32 2kB ...
    lat      (site) float32 2kB ...
    efth     (time, site, freq, dir) float64 1MB ...
    dpt      (time, site) float32 8kB ...
    wspd     (time, site) float32 8kB ...
    wdir     (time, site) float32 8kB ...
Attributes: (12/16)
    product_name:           ww3.all_spec.nc
    area:                   Global 0.5 x 0.5 degree
    data_type:              OCO spectra 2D
    format_version:         1.1
    southernmost_latitude:  n/a
    northernmost_latitude:  n/a
    ...                     ...
    minimum_altitude:       n/a
    maximum_altitude:       n/a
    altitude_resolution:    n/a
    start_date:             2023-01-01 00:00:00
    stop_date:              2023-02-01 00:00:00
    field_type:             3-hourly

xarray.Dataset

• Dimensions:
  • site: 412
  • time: 5
  • freq: 11
  • dir: 8
• Coordinates: (4)
  • time
    (time)
    datetime64[ns]
    2023-01-01 ... 2023-01-02

    standard_name :

    time

    array(['2023-01-01T00:00:00.000000000', '2023-01-01T06:00:00.000000000',
           '2023-01-01T12:00:00.000000000', '2023-01-01T18:00:00.000000000',
           '2023-01-02T00:00:00.000000000'], dtype='datetime64[ns]')

  • site
    (site)
    int64
    0 4 8 12 16 ... 1632 1636 1640 1644

    array([   0,    4,    8, ..., 1636, 1640, 1644], shape=(412,))

  • freq
    (freq)
    float32
    0.05417 0.05959 ... 0.1277 0.1405

    standard_name :

    sea_surface_wave_frequency

    units :

    Hz

    array([0.054172, 0.059589, 0.065548, 0.072103, 0.079313, 0.087244, 0.095968,
           0.105565, 0.116122, 0.127734, 0.140507], dtype=float32)

  • dir
    (dir)
    float32
    0.0 45.0 90.0 ... 225.0 270.0 315.0

    standard_name :

    sea_surface_wave_from_direction

    units :

    degree

    array([  0.,  45.,  90., 135., 180., 225., 270., 315.], dtype=float32)

• Data variables: (6)
  • lon
    (site)
    float32
    ...

    standard_name :

    longitude

    units :

    degrees_east

    [412 values with dtype=float32]

  • lat
    (site)
    float32
    ...

    standard_name :

    latitude

    units :

    degrees_north

    [412 values with dtype=float32]

  • efth
    (time, site, freq, dir)
    float64
    ...

    standard_name :

    sea_surface_wave_directional_variance_spectral_density

    units :

    m2 s degree-1

    [181280 values with dtype=float64]

  • dpt
    (time, site)
    float32
    ...

    standard_name :

    sea_floor_depth_below_sea_surface

    units :

    m

    [2060 values with dtype=float32]

  • wspd
    (time, site)
    float32
    ...

    standard_name :

    wind_speed_at_10m_above_ground_level

    units :

    m s-1

    [2060 values with dtype=float32]

  • wdir
    (time, site)
    float32
    ...

    standard_name :

    wind_from_direction_at_10m_above_ground_level

    units :

    degree

    [2060 values with dtype=float32]

• Indexes: (4)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2023-01-01 00:00:00', '2023-01-01 06:00:00',
                   '2023-01-01 12:00:00', '2023-01-01 18:00:00',
                   '2023-01-02 00:00:00'],
                  dtype='datetime64[ns]', name='time', freq=None))

  • site
    PandasIndex

    PandasIndex(Index([   0,    4,    8,   12,   16,   20,   24,   28,   32,   36,
           ...
           1608, 1612, 1616, 1620, 1624, 1628, 1632, 1636, 1640, 1644],
          dtype='int64', name='site', length=412))

  • freq
    PandasIndex

    PandasIndex(Index([0.05417170003056526, 0.05958886817097664, 0.06554775685071945,
           0.07210253179073334, 0.07931278645992279, 0.08724407106637955,
           0.09596847742795944, 0.10556533187627792, 0.11612186580896378,
           0.12773405015468597, 0.14050745964050293],
          dtype='float32', name='freq'))

  • dir
    PandasIndex

    PandasIndex(Index([0.0, 45.0, 90.0, 135.0, 180.0, 225.0, 270.0, 315.0], dtype='float32', name='dir'))

• Attributes: (16)

  product_name :

  ww3.all_spec.nc

  area :

  Global 0.5 x 0.5 degree

  data_type :

  OCO spectra 2D

  format_version :

  1.1

  southernmost_latitude :

  n/a

  northernmost_latitude :

  n/a

  latitude_resolution :

  n/a

  westernmost_longitude :

  n/a

  easternmost_longitude :

  n/a

  longitude_resolution :

  n/a

  minimum_altitude :

  n/a

  maximum_altitude :

  n/a

  altitude_resolution :

  n/a

  start_date :

  2023-01-01 00:00:00

  stop_date :

  2023-02-01 00:00:00

  field_type :

  3-hourly

In [5]:

# times prescribe the times over which the model will run

times = TimeRange(start="2023-01-01T00", end="2023-01-02T00", interval="6h")
print(times)

# times prescribe the times over which the model will run times = TimeRange(start="2023-01-01T00", end="2023-01-02T00", interval="6h") print(times)

	Start: 2023-01-01 00:00:00
	End: 2023-01-02 00:00:00
	Duration: 1 day
	Interval: 6:00:00
	Include End: True

In [6]:

# grid defines the model grid

grid = SwanGrid(x0=112, y0=-29, dx=0.5, dy=0.5, nx=8, ny=9, rot=20)
grid

# grid defines the model grid grid = SwanGrid(x0=112, y0=-29, dx=0.5, dy=0.5, nx=8, ny=9, rot=20) grid

Out[6]:

SwanGrid: REG, 8x9

In [7]:

# Plot model grid and boundary points in the source spectra

fig, ax = grid.plot(fscale=5)
ax.plot(dset.lon, dset.lat, "ok")
ax.set_extent([108, 120, -35.5, -19])

# Plot model grid and boundary points in the source spectra fig, ax = grid.plot(fscale=5) ax.plot(dset.lon, dset.lat, "ok") ax.set_extent([108, 120, -35.5, -19])

Boundary side with IDW and TPAR files¶

Create BOUNDSPEC SIDE type SWAN boundary by interpolating from the source spectra dataset at the centre of the WEST side of the grid.

The boundary is defined as TPAR files by interpolating at the middle point of the side using the Inverse Distance Weighting (IDW) sel method from wavespectra with a tolerance of 2 degrees within which at least three sites in the source spectra must exist (otherwise missing values are returned).

⚠️ When using IDW, missing values are returned for locations that do not have three or more sites withing tolerance in the source dataset. This is okay for land locations but will be problematic for locations expecting an open boundary. Use a larger tolerance when encountering this problem.

In [8]:

workspace = Path(mkdtemp(dir="model"))

bnd = BoundspecSide(
    id="wa",
    source=source,
    sel_method="idw",
    sel_method_kwargs={"tolerance": 2.0},
    location={"side": "west"},
)

filename, cmd = bnd.get(destdir=workspace, grid=grid, time=times)

display(cmd)
os.system(f"cat {filename}");

workspace = Path(mkdtemp(dir="model")) bnd = BoundspecSide( id="wa", source=source, sel_method="idw", sel_method_kwargs={"tolerance": 2.0}, location={"side": "west"}, ) filename, cmd = bnd.get(destdir=workspace, grid=grid, time=times) display(cmd) os.system(f"cat {filename}");

"BOUNDSPEC SIDE WEST CCW CONSTANT FILE fname='wa_tpar_west_000.bnd' seq=1"

TPAR
20230101.000000 2.74 14.39 223.83 30.07
20230101.060000 2.80 14.03 223.88 31.27
20230101.120000 3.07 14.04 223.26 32.15
20230101.180000 3.48 14.17 223.08 33.61
20230102.000000 3.52 14.21 223.68 33.73

In [9]:

# Load the boundary created as a dataframe

df = read_tpar(filename)
df

# Load the boundary created as a dataframe df = read_tpar(filename) df

Out[9]:

	hs	tp	dpm	dspr
0
2023-01-01 00:00:00	2.74	14.39	223.83	30.07
2023-01-01 06:00:00	2.80	14.03	223.88	31.27
2023-01-01 12:00:00	3.07	14.04	223.26	32.15
2023-01-01 18:00:00	3.48	14.17	223.08	33.61
2023-01-02 00:00:00	3.52	14.21	223.68	33.73

Plot boundary stats, note the values weighted by the four nearest neighbours at the centre of the WEST boundary

In [10]:

for time in df.index[::2]:
    plot_boundary_side(dset, df, grid, time)

for time in df.index[::2]: plot_boundary_side(dset, df, grid, time)

Boundary side with nearest and SPEC2D files¶

Create BOUNDSPEC SIDE type SWAN boundary by nearest-selecting from the source spectra dataset at the centre of the WEST side of the grid.

The boundary is defined as SPEC2D files by nearest-selecting at the middle point of the side using the Nearest neighbour (nearest) sel method from wavespectra with a tolerance of 2 degrees within which at least one site in the source spectra must exist (otherwise an exception is raised).

⚠️ When using nearest, exceptions are raised for locations that do not have at least one site placed withing tolerance in the source dataset. Use a larger tolerance when encountering this problem.

In [11]:

workspace = Path(mkdtemp(dir="model"))

bnd = BoundspecSide(
    id="wa",
    source=source,
    sel_method="nearest",
    sel_method_kwargs={"tolerance": 2.0},
    location={"side": "west"},
    file_type="spec2d",
)

filename, cmd = bnd.get(destdir=workspace, grid=grid, time=times)

display(cmd)
os.system(f"head -n 50 {filename}");

workspace = Path(mkdtemp(dir="model")) bnd = BoundspecSide( id="wa", source=source, sel_method="nearest", sel_method_kwargs={"tolerance": 2.0}, location={"side": "west"}, file_type="spec2d", ) filename, cmd = bnd.get(destdir=workspace, grid=grid, time=times) display(cmd) os.system(f"head -n 50 {filename}");

"BOUNDSPEC SIDE WEST CCW CONSTANT FILE fname='wa_spec2d_west_000.bnd' seq=1"

SWAN   1                                Swan standard spectral file
$   Created by wavespectra
$   
TIME                                    time-dependent data
     1                                  time coding option
LONLAT                                  locations in spherical coordinates
     1                                  number of locations
  111.000000  -27.000000
AFREQ                                   absolute frequencies in Hz
    11                                  number of frequencies
    0.05417
    0.05959
    0.06555
    0.07210
    0.07931
    0.08724
    0.09597
    0.10557
    0.11612
    0.12773
    0.14051
NDIR                                    spectral nautical directions in degr
     8                                  number of directions
     0.0000
    45.0000
    90.0000
   135.0000
   180.0000
   225.0000
   270.0000
   315.0000
QUANT
     1                                  number of quantities in table
VaDens                                  variance densities in m2/Hz/degr
m2/Hz/degr                              unit
   -99                                  exception value
20230101.000000                         date and time
FACTOR
    2.51854154E-05
    0    0    0    0  560 1009    0    0
    0    0    0    0  500 2847    0    0
    0    0    0    0  310 9450    0    0
    0    0    0    0  286 9998    0    0
    0    0    0    0  298 5650    0    0
    0    0    0    4  211 2795    4    8
    4    0    4   48  354 1481    4   24
    4    0   12  238  890 1056    4   16
    0    0   36  596 1660  802    0   16
    0    0   36 1302 2041  580    0    8
    0    0   40 1882 1815  369    0    0

In [12]:

# Load the boundary created as a dataset

ds = wavespectra.read_swan(filename, as_site=True)
ds

# Load the boundary created as a dataset ds = wavespectra.read_swan(filename, as_site=True) ds

Out[12]:

<xarray.Dataset> Size: 4kB
Dimensions:  (time: 5, site: 1, freq: 11, dir: 8)
Coordinates:
  * time     (time) datetime64[ns] 40B 2023-01-01 ... 2023-01-02
  * site     (site) <U18 72B 'wa_spec2d_west_000'
  * freq     (freq) float64 88B 0.05417 0.05959 0.06555 ... 0.1161 0.1277 0.1405
  * dir      (dir) float64 64B 0.0 45.0 90.0 135.0 180.0 225.0 270.0 315.0
Data variables:
    efth     (time, site, freq, dir) float64 4kB 0.0 0.0 0.0 ... 0.01061 0.0 0.0
    lat      (site) float64 8B -27.0
    lon      (site) float64 8B 111.0

xarray.Dataset

• Dimensions:
  • time: 5
  • site: 1
  • freq: 11
  • dir: 8
• Coordinates: (4)
  • time
    (time)
    datetime64[ns]
    2023-01-01 ... 2023-01-02

    array(['2023-01-01T00:00:00.000000000', '2023-01-01T06:00:00.000000000',
           '2023-01-01T12:00:00.000000000', '2023-01-01T18:00:00.000000000',
           '2023-01-02T00:00:00.000000000'], dtype='datetime64[ns]')

  • site
    (site)
    <U18
    'wa_spec2d_west_000'

    array(['wa_spec2d_west_000'], dtype='<U18')

  • freq
    (freq)
    float64
    0.05417 0.05959 ... 0.1277 0.1405

    array([0.05417, 0.05959, 0.06555, 0.0721 , 0.07931, 0.08724, 0.09597, 0.10557,
           0.11612, 0.12773, 0.14051])

  • dir
    (dir)
    float64
    0.0 45.0 90.0 ... 225.0 270.0 315.0

    array([  0.,  45.,  90., 135., 180., 225., 270., 315.])

• Data variables: (3)
  • efth
    (time, site, freq, dir)
    float64
    0.0 0.0 0.0 0.0 ... 0.01061 0.0 0.0

    standard_name :

    sea_surface_wave_directional_variance_spectral_density

    units :

    m2 s degree-1

    _units :

    m^{2}.s.degree^{-1}

    _variable_name :

    VaDens

    array([[[[0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
              0.00000000e+00, 1.41038326e-02, 2.54120841e-02,
              0.00000000e+00, 0.00000000e+00],
             [0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
              0.00000000e+00, 1.25927077e-02, 7.17028776e-02,
              0.00000000e+00, 0.00000000e+00],
             [0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
              0.00000000e+00, 7.80747877e-03, 2.38002176e-01,
              0.00000000e+00, 0.00000000e+00],
             [0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
              0.00000000e+00, 7.20302880e-03, 2.51803783e-01,
              0.00000000e+00, 0.00000000e+00],
             [0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
              0.00000000e+00, 7.50525379e-03, 1.42297597e-01,
              0.00000000e+00, 0.00000000e+00],
             [0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
              1.00741662e-04, 5.31412265e-03, 7.03932360e-02,
              1.00741662e-04, 2.01483323e-04],
             [1.00741662e-04, 0.00000000e+00, 1.00741662e-04,
              1.20889994e-03, 8.91563705e-03, 3.72996002e-02,
    ...
              8.97669742e-04, 1.20021769e-02, 1.75411317e-01,
              0.00000000e+00, 0.00000000e+00],
             [0.00000000e+00, 0.00000000e+00, 2.99223248e-04,
              7.31434605e-03, 2.06131571e-02, 9.00994445e-02,
              0.00000000e+00, 0.00000000e+00],
             [0.00000000e+00, 0.00000000e+00, 1.09715191e-03,
              2.58994344e-02, 7.55039995e-02, 5.09011991e-02,
              9.97410825e-05, 1.99482165e-04],
             [0.00000000e+00, 0.00000000e+00, 1.39637516e-03,
              5.26965386e-02, 1.12906905e-01, 3.62060129e-02,
              9.97410825e-05, 2.99223248e-04],
             [0.00000000e+00, 0.00000000e+00, 1.09715191e-03,
              9.82117192e-02, 1.09814932e-01, 2.48022825e-02,
              0.00000000e+00, 2.99223248e-04],
             [0.00000000e+00, 0.00000000e+00, 1.89508057e-03,
              9.40890878e-02, 8.09897590e-02, 1.42962218e-02,
              0.00000000e+00, 1.99482165e-04],
             [0.00000000e+00, 0.00000000e+00, 4.08938438e-03,
              4.89063775e-02, 5.28960208e-02, 1.06058018e-02,
              0.00000000e+00, 0.00000000e+00]]]])

  • lat
    (site)
    float64
    -27.0

    standard_name :

    latitude

    units :

    degrees_north

    array([-27.])

  • lon
    (site)
    float64
    111.0

    standard_name :

    longitude

    units :

    degrees_east

    array([111.])

• Indexes: (4)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2023-01-01 00:00:00', '2023-01-01 06:00:00',
                   '2023-01-01 12:00:00', '2023-01-01 18:00:00',
                   '2023-01-02 00:00:00'],
                  dtype='datetime64[ns]', name='time', freq=None))

  • site
    PandasIndex

    PandasIndex(Index(['wa_spec2d_west_000'], dtype='object', name='site'))

  • freq
    PandasIndex

    PandasIndex(Index([0.05417, 0.05959, 0.06555,  0.0721, 0.07931, 0.08724, 0.09597, 0.10557,
           0.11612, 0.12773, 0.14051],
          dtype='float64', name='freq'))

  • dir
    PandasIndex

    PandasIndex(Index([0.0, 45.0, 90.0, 135.0, 180.0, 225.0, 270.0, 315.0], dtype='float64', name='dir'))

• Attributes: (0)

In [13]:

# Calculate boundary stats from dataset for plotting

df = ds.isel(site=0, drop=True).spec.stats(["hs", "tp", "dpm", "dspr"]).to_pandas()
df

# Calculate boundary stats from dataset for plotting df = ds.isel(site=0, drop=True).spec.stats(["hs", "tp", "dpm", "dspr"]).to_pandas() df

Out[13]:

	hs	tp	dpm	dspr
time
2023-01-01 00:00:00	2.709378	14.361055	223.864182	31.245212
2023-01-01 06:00:00	2.799166	13.995838	223.720566	32.982870
2023-01-01 12:00:00	3.009826	14.025472	223.009888	33.819244
2023-01-01 18:00:00	3.453083	14.141555	222.878220	35.180063
2023-01-02 00:00:00	3.510868	14.194441	223.638062	35.336880

Plot boundary stats, note the values are identical to the nearest neighbour at the centre of the WEST boundary

In [14]:

for time in df.index[::2]:
    plot_boundary_side(dset, df, grid, time)

for time in df.index[::2]: plot_boundary_side(dset, df, grid, time)

No nearest neighbours¶

When the nearest neighbour sel method is choosen there must be at least one site in the source spectra within the radius defined by tolerance otherwise an exception is raised

In [15]:

workspace = Path(mkdtemp(dir="model"))

bnd = BoundspecSide(
    id="wa",
    source=source,
    sel_method="nearest",
    sel_method_kwargs={"tolerance": 0.1},
    location={"side": "west"},
)

try:
    filename, cmd = bnd.get(destdir=workspace, grid=grid, time=times)
except AssertionError as e:
    print(e)

workspace = Path(mkdtemp(dir="model")) bnd = BoundspecSide( id="wa", source=source, sel_method="nearest", sel_method_kwargs={"tolerance": 0.1}, location={"side": "west"}, ) try: filename, cmd = bnd.get(destdir=workspace, grid=grid, time=times) except AssertionError as e: print(e)

Nearest site in dataset from ((np.float64(111.31595971334866), np.float64(-27.120614758428182))) is 0.338199 deg away but tolerance is 0.1 deg

The idw sel method requires at least two neighbours within the radius defined by tolerance however it does not raise an exception if those are not found but return missing values instead

In [16]:

bnd = BoundspecSide(
    id="wa",
    source=source,
    sel_method="idw",
    sel_method_kwargs={"tolerance": 0.1},
    location={"side": "west"},
)

filename, cmd = bnd.get(destdir=workspace, grid=grid, time=times)

os.system(f"cat {filename}");

bnd = BoundspecSide( id="wa", source=source, sel_method="idw", sel_method_kwargs={"tolerance": 0.1}, location={"side": "west"}, ) filename, cmd = bnd.get(destdir=workspace, grid=grid, time=times) os.system(f"cat {filename}");

TPAR
20230101.000000 0.00 0.0 0.0 0.0
20230101.060000 0.00 0.0 0.0 0.0
20230101.120000 0.00 0.0 0.0 0.0
20230101.180000 0.00 0.0 0.0 0.0
20230102.000000 0.00 0.0 0.0 0.0