global TOMAWAC simulations¶
We'll generate a world contiguous mesh for global wave simulations.
install base environment¶
To be able to run this notebook, first install the pre-requisites for seamsh:
mamba create -n tomawac_mesh gdal gmsh scipy python=3.12
Optionally¶
create then a virtual environment on top of conda:
python -mvenv .venv
source .venv/bin/activate
install libraries¶
pip install -e .
You're ready to run the notebook, and generate different version of the mesh.
This notebook contains to following sections:
- gather data:
- A. GEBCO bathy through
seareport_data - B. coastlines from natural earth coasltines
- A. GEBCO bathy through
- Prepare the distance fields from the above data
- only coastline is used to constrain with the Distance from the shore, but you can activate more contraints if necessary
- Mesh
- Interpolate bathy onto the mesh
- Convert to Selafin (using
xarray-selafinavailable viapip) - Generate the CLI file
- Interpolate ERA5 wind data onto the mesh
- Convert to wind Selafin binary
- Visualise TOMAWAC results directly in the notebook:
1 - Generate the mesh¶
Particularities of the meshes available in the mesh/ folder:¶
- For TOMAWAC simulations, I needed to add a hole in the pole. Only because of the triangle(s) that contain the pole (bug in routine
GEOELT.f). - The mesh is contiguous, e.g. it wrap around the globe and the triangles at the dateline (+/- 180°) are connected alltogether.
- The only constraint for the mesh size is the distance from the coast. You can activate the wavelength contrainst (proportional to
sqrt(depth)) or the bathy gradient constraint.
2 - Run the model¶
You will need to compile and run TELEMAC from the main branch, because it contains the latest version of 1D time series exports, very useful for global model analysis (and comparison against buoys).
All the modified routines to run TOMAWAC global are in wac/princi/ folder.
Look for SEB or TOM for the modified parts of the source code.
Get the whole project:¶
available on Google drive : https://drive.google.com/drive/folders/1SMNkHymAeuRfKh27pwPf7TgdIlOhTBvI?usp=sharing
In [1]:
import seamsh
import numpy as np
import xarray as xr
import seareport_data as D
import utils
import hvplot.pandas
import hvplot.xarray
Mesh settings¶
here are the main parameters used for the mesher settings:
In [2]:
# Meshing parameters for contributions
OPTS = {
"max": 150000, # max resolution
"min": 10000, # min resolution
"dist": 0.15, # The rate of expansion in decimal percent from the shoreline
"wl": 50, # (not used here) number of element to resolve WL
"m2": 12.42 * 3600, # (not used here) for wavelength: M2 period in seconds
"g": 9.81, # (not used here) for wavelength: m/s^2
"slope": 5, # (not used here) number of element to resolve bathy gradient
"grade": 1.5, # the rate of growth in decimal percent
}
REMOVE_POLE = True # for tomawac simulations
Coastlines¶
In [3]:
SHP_IN = "ne_10m_coastline/ne_10m_coastline.shp"
SHPOUT = "ne_10m_coastline/ne_10m_coastline_poly.shp"
POLE= "ne_10m_coastline/pole_ring.shp"
utils.get_ne_coastline(SHP_IN)
utils.remove_small_islands(SHP_IN, SHPOUT, OPTS['min'])
utils.add_hole_in_pole(POLE)
define the CRS for the meshing¶
In [ ]:
from osgeo import osr
domain_srs = osr.SpatialReference()
domain_srs.ImportFromProj4("+ellps=WGS84 +proj=stere +lat_0=90")
# "Cartesian" projection (i.e. no projection) is used to compute the distance to the coast.
cart_srs = osr.SpatialReference()
cart_srs.ImportFromProj4("+ellps=WGS84 +proj=cart +units=m +x_0=0 +y_0=0")
wgs84_srs = osr.SpatialReference()
wgs84_srs.ImportFromProj4("+ellps=WGS84 +proj=longlat +datum=WGS84 +no_defs")
Distance field¶
In [ ]:
domain = seamsh.geometry.Domain(domain_srs)
domain.add_boundary_curves_shp(SHPOUT, "featurecla", seamsh.geometry.CurveType.POLYLINE)
dist_coast = seamsh.field.Distance(domain, OPTS["min"]/2, projection=cart_srs)
Bathy field (not used here)¶
In [6]:
bathy = D.gebco_ds("sub_ice")
bathy_subset = bathy.isel(lon=slice(0,-1,10), lat=slice(0,-1,10))
bathy_grad = utils.calc_gradient(bathy_subset, 'elevation')
bathy_grad
Out[6]:
<xarray.Dataset> Size: 373MB
Dimensions: (lat: 4320, lon: 8640)
Coordinates:
* lon (lon) float64 69kB -180.0 -180.0 -179.9 ... 179.9 179.9 180.0
* lat (lat) float64 35kB -90.0 -89.96 -89.91 ... 89.88 89.92 89.96
Data variables:
crs |S1 1B ...
elevation (lat, lon) int16 75MB -18 -18 -18 -18 ... -4212 -4212 -4212 -4212
gradient (lat, lon) float64 299MB 0.02388 0.02409 ... 2.912e-08 2.912e-08
Attributes: (12/36)
title: The GEBCO_2025 Grid - a continuous terra...
summary: The GEBCO_2025 Grid is a continuous, glo...
keywords: BATHYMETRY/SEAFLOOR TOPOGRAPHY, DIGITAL ...
Conventions: CF-1.6, ACDD-1.3
id: DOI: 10.5285/37c52e96-24ea-67ce-e063-708...
naming_authority: https://dx.doi.org
... ...
geospatial_vertical_units: meters
geospatial_vertical_resolution: 1.0
geospatial_vertical_positive: up
identifier_product_doi: DOI: 10.5285/37c52e96-24ea-67ce-e063-708...
references: DOI: 10.5285/37c52e96-24ea-67ce-e063-708...
node_offset: 1.0In [ ]:
! mkdir -p data
In [ ]:
BATHY = "data/gebco_2024_1000_4k.tif"
BATHY_GRADIENT = "data/gebco_2024_grad_smooth_1000_4k.tif"
utils.to_raster(bathy_subset.elevation, BATHY)
utils.to_raster(bathy_grad.gradient, BATHY_GRADIENT)
In [ ]:
bath_field = seamsh.field.Raster(BATHY)
bath_field._projection = osr.SpatialReference()
bath_field._projection.ImportFromProj4("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
grad_field = seamsh.field.Raster(BATHY_GRADIENT)
grad_field._projection = osr.SpatialReference()
grad_field._projection.ImportFromProj4("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
# %%
MESH¶
Mesh size function¶
In [10]:
# it is necessary to convert the mesh size to stereographic coordinates
def mesh_size(x, projection,R = 6371000):
depth = bath_field(x, projection)
depth[depth>0] = 0
grad = grad_field(x, projection)
lon, lat = utils.stereo_to_wgs84(x, source_epsg=projection)
s_wl = OPTS["m2"] / OPTS["wl"] * np.sqrt(9.81 * - depth)
s_wl = np.clip(s_wl, OPTS["min"]*2,None)
s_grad = abs(depth) / (grad + 1e-10)*(2*np.pi / OPTS["slope"])
s_grad = np.clip(s_grad, OPTS["min"]*3,None)
s_coast = dist_coast(x, projection)* OPTS["grade"] + OPTS["min"]
# s_final = np.c_[s_coast, s_wl, s_grad].min(axis=1)
s_final = s_coast
stereo_factor = 2/(1+x[:,0]**2/R**2+x[:,1]**2/R**2) * (3/4*np.cos(np.deg2rad(lat+90))+5/4)
return np.clip(s_final,OPTS["min"],OPTS["max"])/stereo_factor
Coarsen boundaries
In [ ]:
coarse = seamsh.geometry.coarsen_boundaries(domain, (0, 0), domain_srs, mesh_size)
if REMOVE_POLE:
coarse.add_boundary_curves_shp(POLE,"featurecla", seamsh.geometry.CurveType.POLYLINE)
else:
x = [[0., 0.]]
coarse.add_interior_points(x,"featurecla",domain_srs)
finally mesh
In [12]:
seamsh.gmsh.mesh(
coarse,
"natural_earth.msh",
mesh_size,
output_srs=domain_srs,
intermediate_file_name="-",
# binary=True
)
seamsh.gmsh.reproject("natural_earth.msh", domain_srs, "natural_earth_wgs84.msh", wgs84_srs)
* Generate mesh * Build topology Build gmsh model Build mesh size field Mesh with gmsh Write "natural_earth.msh" (msh version 4.1)
Interpolate Bathy¶
In [13]:
import gmsh
gmsh.open("natural_earth_wgs84.msh")
tri_i, tri_n = gmsh.model.mesh.getElementsByType(2)
tri_n = tri_n.reshape([-1, 3])
node_i, nodes, _ = gmsh.model.mesh.getNodes()
nodes = nodes.reshape([-1, 3])
x = nodes[:, 0]
y = nodes[:, 1]
z = nodes[:, 2]
element = np.subtract(tri_n, 1)
In [14]:
subset = bathy.isel(lon=slice(0, -1, 10), lat=slice(0, -1, 10)).pad(lon=(1,1),lat=(1,1),mode="edge")
lon_values = subset.lon.values
lat_values = subset.lat.values
lon_values[-1] = 180.0
lon_values[0] = -180.0
lat_values[-1] = 90.0
lat_values[0] = -90.0
subset = subset.assign_coords(lon=lon_values,lat=lat_values)
subset
Out[14]:
<xarray.Dataset> Size: 75MB
Dimensions: (lat: 4322, lon: 8642)
Coordinates:
* lon (lon) float64 69kB -180.0 -180.0 -180.0 ... 179.9 180.0 180.0
* lat (lat) float64 35kB -90.0 -90.0 -89.96 -89.91 ... 89.92 89.96 90.0
Data variables:
crs |S1 1B ...
elevation (lat, lon) int16 75MB -18 -18 -18 -18 ... -4212 -4212 -4212 -4212
Attributes: (12/36)
title: The GEBCO_2025 Grid - a continuous terra...
summary: The GEBCO_2025 Grid is a continuous, glo...
keywords: BATHYMETRY/SEAFLOOR TOPOGRAPHY, DIGITAL ...
Conventions: CF-1.6, ACDD-1.3
id: DOI: 10.5285/37c52e96-24ea-67ce-e063-708...
naming_authority: https://dx.doi.org
... ...
geospatial_vertical_units: meters
geospatial_vertical_resolution: 1.0
geospatial_vertical_positive: up
identifier_product_doi: DOI: 10.5285/37c52e96-24ea-67ce-e063-708...
references: DOI: 10.5285/37c52e96-24ea-67ce-e063-708...
node_offset: 1.0In [15]:
from scipy.interpolate import RegularGridInterpolator as RGI
RG = RGI((subset.lon, subset.lat), subset.elevation.data.T, method="linear")
export to Selafin¶
In [16]:
!mkdir -p mesh
In [17]:
from xarray_selafin.xarray_backend import SelafinAccessor
import pandas as pd
slf_ds = xr.Dataset(
coords={
"time": [pd.Timestamp.now()],
"x": ("node", x),
"y": ("node", y),
},
data_vars={
"B" : (("time", "node"), [RG((x,y))])
}
)
slf_ds.attrs['ikle2'] = element + 1
filebase = f"mesh/global_{OPTS["min"]/1000}-{int(OPTS["max"]/1000)}km-smo{OPTS["grade"]-1}-{len(slf_ds.x)}nodes"
slf_ds.selafin.write(filebase+".slf")
# export CLI file
_=utils.export_cli(slf_ds, filebase+".cli")
In [18]:
slf_ds
Out[18]:
<xarray.Dataset> Size: 3MB
Dimensions: (time: 1, node: 133059)
Coordinates:
* time (time) datetime64[ns] 8B 2025-09-26T18:09:23.790500
x (node) float64 1MB -26.24 -156.1 -66.72 ... -2.335 110.9 -60.49
y (node) float64 1MB -58.48 -85.13 -78.44 -80.3 ... 49.1 79.26 9.148
Dimensions without coordinates: node
Data variables:
B (time, node) float64 1MB -155.3 -871.3 -361.1 ... -1.926e+03 -5.747
Attributes:
ikle2: [[106330 113322 86664]\n [ 1957 124004 18862]\n [ 64490 1239...In [19]:
filebase
Out[19]:
'mesh/global_10.0-150km-smo0.5-133059nodes'
In [ ]:
slf_ds.B.hvplot.scatter(
x="x",
y="y",
c='B',
cmap="fire",
s=1,
clim=(-5000, 0),
).opts(
width=800,
height=500,
)
Out[ ]:
BokehModel(combine_events=True, render_bundle={'docs_json': {'77a0288a-0fd2-480d-92b0-03ec038329b6': {'version…
Interpolate wind¶
In [21]:
atm_nc = xr.open_dataset("data/era5_202307_uvp.nc")
atm_nc
Out[21]:
<xarray.Dataset> Size: 19GB
Dimensions: (time: 744, latitude: 721, longitude: 1440)
Coordinates:
* longitude (longitude) float32 6kB 0.0 0.25 0.5 0.75 ... 359.2 359.5 359.8
* latitude (latitude) float32 3kB 90.0 89.75 89.5 ... -89.5 -89.75 -90.0
* time (time) datetime64[ns] 6kB 2023-07-01 ... 2023-07-31T23:00:00
Data variables:
msl (time, latitude, longitude) float64 6GB ...
u10 (time, latitude, longitude) float64 6GB ...
v10 (time, latitude, longitude) float64 6GB ...
Attributes:
Conventions: CF-1.6
history: 2024-04-04 09:42:07 GMT by grib_to_netcdf-2.25.1: /opt/ecmw...In [22]:
lon_values = atm_nc.longitude.values
lat_values = atm_nc.latitude.values
lon_values[-1] = 360.0
atm_nc = atm_nc.assign_coords(lon=lon_values)
In [23]:
atm_data = dict()
for var in ["msl", "u10", "v10"]:
RG = RGI((atm_nc.longitude, atm_nc.latitude), atm_nc[var].data.T, method="linear")
atm_data[var] = RG((slf_ds.x.data % 360, slf_ds.y.data))
slf_wind =xr.Dataset(coords = {
"x": ("node", slf_ds.x.data),
"y": ("node", slf_ds.y.data),
"time": atm_nc.time.data
},
data_vars={
"PATM": (("node", "time"), atm_data["msl"]),
"WINDX": (("node", "time"), atm_data["u10"]),
"WINDY": (("node", "time"), atm_data["v10"]),
}
)
slf_wind.attrs["ikle2"] = slf_ds.ikle2
slf_wind.attrs["variables"] = {'PATM': ('PATM', 'PASCAL'), 'WINDX': ('WINDX', 'M/S'), 'WINDY': ('WINDY', 'M/S')}
In [24]:
slf_wind.selafin.write(filebase+"_wind.slf")
Inspect results¶
In [25]:
res2d = xr.open_dataset("wac/r2d_global_tom.slf")
res1d = xr.open_dataset("wac/r1d_global_tom.slf")
In [26]:
res2d.max(dim="time").WH.hvplot.scatter(x="x",
y="y",
c='WH',
cmap="rainbow4",
s=1,
).opts(
clim=(0, 10),
width=800,
height=500,
title="Max Hm0 over July 2023"
)
Out[26]:
In [27]:
res1d.max(dim="time").WH.hvplot.scatter(x="x",
y="y",
c='WH',
cmap="rainbow4",
hover_cols=["node"]
).opts(
clim=(0, 10),
width=800,
height=500,
title="Max Hm0 over July 2023"
)
Out[27]:
In [28]:
res1d.isel(node=389).WH.hvplot(
).opts(
width=800,
height=500,
title=f"Hm0 over July 2023, for station located at {res1d.isel(node=389).x.values:0.2f}°E, {res1d.isel(node=389).y.values:0.2f}°N"
)
Out[28]: