to have this notebook working correctly, you just need to install xarray-selafin, numpy and pandas:
pip install xarray-selafin pandas numpy
check how to create Selafins from xarrays also on the wiki at https://github.com/seareport/xarray-selafin/wiki
In [ ]:
import xarray as xr
In [ ]:
HOMETEL = "/home/tomsail/work/opentelemac/"
file = HOMETEL + "examples/python3/converter/mesh.slf"
In [ ]:
slf = xr.open_dataset(file, engine='selafin')
In [ ]:
slf
Out[Â ]:
<xarray.Dataset> Size: 596kB
Dimensions: (time: 1, node: 13541)
Coordinates:
x (node) float32 54kB ...
y (node) float32 54kB ...
* time (time) datetime64[ns] 8B 1900-01-01
Dimensions without coordinates: node
Data variables:
VITESSE U (time, node) float32 54kB ...
VITESSE V (time, node) float32 54kB ...
HAUTEUR D'EAU (time, node) float32 54kB ...
SURFACE LIBRE (time, node) float32 54kB ...
FOND (time, node) float32 54kB ...
FROUDE (time, node) float32 54kB ...
DEBIT SCALAIRE (time, node) float32 54kB ...
DEBIT SUIVANT X (time, node) float32 54kB ...
DEBIT SUIVANT Y (time, node) float32 54kB ...
Attributes:
title: TELEMAC 2D : RUPTURE DE BARRAGE SUR FOND SEC$
language: en
float_size: 4
endian: >
params: (1, 0, 0, 0, 0, 0, 0, 0, 0, 0)
ipobo: [763 765 908 ... 0 0 0]
ikle2: [[ 2193 962 2191]\n [ 1103 2189 648]\n [12311 12298 13...
variables: {'VITESSE U': ('VITESSE U', 'M/S'), 'VITESSE V': ('VITESSE V...
date_start: (1900, 1, 1, 0, 0, 0)In [ ]:
xmin, xmax = slf.x.values.min(), slf.x.values.max()
ymin, ymax = slf.y.values.min(), slf.y.values.max()
xmin, xmax, ymin, ymax
Out[Â ]:
(536.4716, 17763.07, -2343.54, 6837.79)
we can take a 1000m resolution grid that extends accordingly:
In [ ]:
bbox = (0, 20000, -3000, 7000) # xmin xmax ymin ymax
let's generate a wind from scratch
In [ ]:
import xarray as xr
import numpy as np
import pandas as pd
from scipy.spatial import Delaunay
from xarray_selafin.xarray_backend import SelafinAccessor
# Define the mesh grid for the bbox
xmin, xmax, ymin, ymax = 0, 20000, -3000, 7000
resolution = 1000 # in meters
lon = np.arange(xmin, xmax + resolution, resolution)
lat = np.arange(ymin, ymax + resolution, resolution)
nx1d = len(lon)
ny1d = len(lat)
xx = np.tile(lon, ny1d).reshape(ny1d, nx1d).T.ravel()
yy = np.tile(lat, nx1d)
# connectivity
ikle = np.zeros((2*(nx1d-1)*(ny1d-1), 3), dtype=np.int64)
ielem = 0
for i in range(1, nx1d):
for j in range(1, ny1d):
ipoin = (i-1)*ny1d + j - 1
# ~~> first triangle
ikle[ielem][0] = ipoin
ikle[ielem][1] = ipoin + ny1d
ikle[ielem][2] = ipoin + 1
ielem = ielem + 1
# # ~~> second triangle
ikle[ielem][0] = ipoin + ny1d
ikle[ielem][1] = ipoin + ny1d + 1
ikle[ielem][2] = ipoin + 1
ielem = ielem + 1
# Boundaries
ipob3 = np.zeros(2*(nx1d-1)*(ny1d-1), dtype=np.int64)
# ~~> along the x-axis (lon)
for i in range(nx1d):
ipoin = i*ny1d
ipob3[ipoin] = i + 1
ipoin = i*ny1d - 1
ipob3[ipoin] = 2*nx1d+(ny1d-2) - i
# ~~> along the y-axis (alt)
for i in range(1, ny1d):
ipoin = i
ipob3[ipoin] = 2*nx1d + 2*(ny1d-2) - i + 1
ipoin = ny1d*(nx1d-1) + i
ipob3[ipoin] = nx1d + i
In [ ]:
# Define the time range
times = pd.date_range(pd.Timestamp(2023, 1, 1), pd.Timestamp(2023, 2, 1), freq="1D")
In [ ]:
# Function to calculate wind components
def wind_components(speed, direction):
radian_direction = np.deg2rad(direction)
u = speed * np.sin(radian_direction) # U component of the wind (east-west)
v = speed * np.cos(radian_direction) # V component of the wind (north-south)
return u, v
In [ ]:
# Define the initial wind direction and speed
initial_wind_direction = 90 # degrees from north, eastward
wind_speed = 20 # m/s
# Calculate the wind field for each day
wind_u = np.zeros((len(times), len(xx)))
wind_v = np.zeros((len(times), len(xx)))
patm = np.zeros((len(times), len(xx)))
tair = np.zeros((len(times), len(xx)))
for i, time in enumerate(times):
# Calculate the new wind direction
new_direction = (initial_wind_direction + i * 10) % 360
spacial_variability = np.sort(np.random.normal(1.9, 1.02, len(xx)))
# Calculate the wind components
u, v = wind_components(wind_speed, new_direction)
# Set the components for all nodes
wind_u[i, :] = u + 5 * spacial_variability
wind_v[i, :] = v + 5 * spacial_variability
patm[i, :] = 102500 + 100 * spacial_variability
tair[i, :] = 20 + 5 * spacial_variability
In [ ]:
# Create the xarray dataset
ds = xr.Dataset(
{
"WINDX": (("time", "node"), wind_u),
"WINDY": (("time", "node"), wind_v),
"PATM": (("time", "node"), patm),
"TAIR": (("time", "node"), tair),
},
coords={
"x": ("node", xx),
"y": ("node", yy),
"time": times,
}
)
plot the wind field to check
In [ ]:
slf.plot.scatter(x="x", y="y",c=slf.FOND,s=4, cmap="jet", edgecolors='none' ,figsize=(20, 10))
ds.isel(time=-1).plot.quiver(x="x", y="y", u="WINDX", v="WINDX", hue="TAIR", zorder = -1)
Out[Â ]:
<matplotlib.quiver.Quiver at 0x78ba3cbbf0d0>
Important! (otherwise it won't work)
- add connectivity and variables attributes to the dataset
In [ ]:
var_attrs = {
"WINDX": ("WINDX", "M/S"),
"WINDY": ("WINDY", "M/S"),
"PATM": ("PATM", "PASCAL"),
"TAIR": ("TAIR", "DEGREES C"),
}
ds.attrs['ikle2'] = ikle + 1 # convert to 1-based indexing
ds.attrs["variables"] = var_attrs
In [ ]:
# Writing to a SELAFIN file
ds.selafin.write("wind.slf")
(array([], dtype=int64), array([], dtype=int64)) SERAFIN REQUEST ERROR: Unexpected error while determining next boundary node after node 4 The IPOBO table could not be built and set to 0. Check if some nodes are superimposed. SERAFIN REQUEST ERROR: Unexpected error while determining next boundary node after node 4 The IPOBO table could not be built and set to 0. Check if some nodes are superimposed.
ignore the warning about ipobo: we don't need it for this example
Bonus: interactive view with holoviews¶
In [ ]:
import holoviews as hv
hv.extension('bokeh')
In [ ]:
df = slf.isel(time=0).to_dataframe()
scatter_plot = hv.Points(
df, ['x', 'y'],['FOND'],
).opts(
color="FOND",
cmap='jet',
line_width=1,
size=2,
tools=["hover"],
show_legend=True,
hover_fill_color='firebrick',
)
scatter_plot.opts(width=1000, height=600)
Out[Â ]: