In [38]:
from __future__ import annotations
import hvplot.pandas
import numpy as np
import holoviews as hv
import pandas as pd
import geopandas as gp
import shapely
import panel as pn
from shapely.geometry import Polygon, Point
import typing as T
# from pyposeidon.utils.cfl import parse_hgrid
VOSPOROS_SOUTH = shapely.Polygon([
(26.77069525861703, 40.44054922299234),
(26.563328315257657, 40.4133689586624),
(26.115328982567448, 39.99024357789592),
(26.394107058739323, 39.98182584129004),
(26.785494998192448, 40.342851238592424),
(26.77069525861703, 40.44054922299234),
])
def extract_area(x: np.ndarray, y: np.ndarray, triangles: np.ndarray, lon_range: tuple[float, float], lat_range: tuple[float, float]):
mask = (x >= lon_range[0]) & (x <= lon_range[1]) & \
(y >= lat_range[0]) & (y <= lat_range[1])
node_indices = np.where(mask)[0]
node_map = {old: new for new, old in enumerate(node_indices)}
node_map_reverse = {new: old for new, old in enumerate(node_indices)}
extracted_x = x[node_indices]
extracted_y = y[node_indices]
triangle_indices = np.array([idx for idx, tri in enumerate(triangles) if all(i in node_map for i in tri)])
extracted_triangles = triangles[triangle_indices]
extracted_triangles = np.array([[node_map[i] for i in tri] for tri in extracted_triangles])
return extracted_x, extracted_y, extracted_triangles, node_indices, triangle_indices, node_map_reverse
# 2. Convert triangulation to GeoDataFrame with MultiPolygons
def tri_to_geopandas(points, tris, values):
"""Convert triangulation to GeoDataFrame with values for coloring"""
polygons = [Polygon(points[tri]) for tri in tris]
gdf = gp.GeoDataFrame({
'A': values[tris[:,0]],
'B': values[tris[:,1]],
'C': values[tris[:,2]],
'geometry': polygons
},
crs="EPSG:4326" # using WGS84 for demonstration
)
return gdf
def drop(nodes, elems, drop_indices):
drop_indices = np.unique(drop_indices)
keep_mask = np.ones(len(nodes), dtype=bool)
keep_mask[drop_indices] = False
nodes = nodes[keep_mask]
old_to_new = -np.ones(len(keep_mask), dtype=int)
old_to_new[np.where(keep_mask)[0]] = np.arange(len(nodes))
mask_valid = np.all(np.isin(elems, np.where(keep_mask)[0]), axis=1)
elems = elems[mask_valid]
elems = old_to_new[elems]
return nodes, elems
def render_gdf(gdf, case = "global"):
if case == 'global':
color = 'red'
elif case == 'black sea':
color = "blue"
else:
raise ValueError("case has to be 'global' or 'black sea'!")
if isinstance(gdf, gp.GeoDataFrame):
return gdf.hvplot(geo=True, tiles=True, c=color, label=case)
elif isinstance(gdf, pd.DataFrame):
return gdf.hvplot.points(geo=True, tiles=True, c=color, line_color='k', tools=["box_select"], hover_cols = ['index'], label=case)
def subset_gdf(gdf, poly: shapely.Polygon):
return gdf[gdf.geometry.within(poly)]
def polygon_to_points(gdf):
coords_index = []
for g in gdf.iterrows():
coords_index.append([*g[1].geometry.exterior.xy, np.array([g[1].A,g[1].B,g[1].C,g[1].A])])
pp = np.array(coords_index).transpose(0, 2, 1).reshape(-1, 3)
return pd.DataFrame(pp, columns=['lon', 'lat', 'index'])
global mesh is taken from: https://github.com/seareport/seareport_models
Mesh v3.2 1km resolution at the coastlines
In [ ]:
file = "seareport_models/v3.2/GSHHS_f_0.01_final.gr3"
mesh_dic = parse_hgrid(file)
In [40]:
x, y, depth = mesh_dic['nodes'].T
tris = mesh_dic["elements"]
subset global mesh to vosphoros south strait
In [41]:
x_, y_, tri_, ind_, tri_sub, node_mapping = extract_area(x,y,tris, (25,45), (40,48))
gdf_ = tri_to_geopandas(np.vstack((x_, y_)).T,tri_, ind_)
global_vosphoros_gdf = subset_gdf(gdf_, VOSPOROS_SOUTH)
global_vosphoros_points = polygon_to_points(global_vosphoros_gdf)
In [42]:
points_ = render_gdf(global_vosphoros_points)
global_vosphoros_ = render_gdf(global_vosphoros_gdf)
(global_vosphoros_ * points_).opts(width = 1000, height=800)
Out[42]:
black sea mesh is taken from: https://github.com/tomsail/regional_seamsh/tree/black_sea_bathy
In [43]:
from pyposeidon.mgmsh import read_msh
black_sea_mesh = read_msh("out_blacksea_0.2.msh")
bs_pts = np.vstack((black_sea_mesh.SCHISM_hgrid_node_x, black_sea_mesh.SCHISM_hgrid_node_y)).T
bs_tri = black_sea_mesh.SCHISM_hgrid_face_nodes.data
black_sea_mesh
Info : Reading 'out_blacksea_0.2.msh'... Info : 7718 entities Info : 48945 nodes Info : 97915 elements Info : Storing section $Projection as model attribute 2025-05-13 12:26:33,951 INFO pyposeidon.mgmsh read_msh:76 Analyzing mesh Info : Done reading 'out_blacksea_0.2.msh' 2025-05-13 12:26:33,957 INFO pyposeidon.mgmsh read_msh:99 No of Open boundaries: 6 2025-05-13 12:26:33,962 INFO pyposeidon.mgmsh read_msh:115 No of Land boundaries: 0 2025-05-13 12:26:33,963 INFO pyposeidon.mgmsh read_msh:131 No. of Islands: 0 2025-05-13 12:26:33,966 INFO pyposeidon.mgmsh read_msh:192 Finalizing Dataset
/home/tomsail/miniconda3/envs/pos_test/lib/python3.11/site-packages/numpy/core/fromnumeric.py:59: FutureWarning: 'DataFrame.swapaxes' is deprecated and will be removed in a future version. Please use 'DataFrame.transpose' instead. return bound(*args, **kwds)
Out[43]:
<xarray.Dataset> Size: 4MB
Dimensions: (nSCHISM_hgrid_node: 48945,
nSCHISM_hgrid_face: 91753,
nMaxSCHISM_hgrid_face_nodes: 3, bnodes: 6162)
Coordinates:
* bnodes (bnodes) int64 49kB 0 1 2 3 ... 6284 6287 6288 6289
Dimensions without coordinates: nSCHISM_hgrid_node, nSCHISM_hgrid_face,
nMaxSCHISM_hgrid_face_nodes
Data variables:
SCHISM_hgrid_node_x (nSCHISM_hgrid_node) float64 392kB 28.81 ... 29.07
SCHISM_hgrid_node_y (nSCHISM_hgrid_node) float64 392kB 41.3 ... 40.94
SCHISM_hgrid_face_nodes (nSCHISM_hgrid_face, nMaxSCHISM_hgrid_face_nodes) uint64 2MB ...
depth (nSCHISM_hgrid_node) float64 392kB 0.0 0.0 ... 0.0
node (bnodes) uint64 49kB 145 146 147 ... 521 522 523
type (bnodes) object 49kB 'open' 'open' ... 'open'
id (bnodes) int64 49kB 1 1 1 1 1 2 2 ... 5 5 5 5 6 6 6In [44]:
gdf_ = tri_to_geopandas(bs_pts, bs_tri, np.arange(len(bs_pts)))
black_sea_vosporos_gdf = subset_gdf(gdf_, VOSPOROS_SOUTH)
black_sea_vosporos_points = polygon_to_points(black_sea_vosporos_gdf)
black_sea_vosporos_ = render_gdf(black_sea_vosporos_gdf, "black sea")
In [45]:
selection = hv.streams.Selection1D(source=points_)
# Create a function to print selected points
def print_all_points(df: pd.DataFrame, indices: pd.Index[int], text_box: pn.widgets.TextAreaInput) -> T.Any:
print("Selected indices:", indices)
if indices:
ind = np.unique(indices)
list = df["index"].iloc[ind].astype(int).unique()
value = ",\n".join(map(str, list))
else:
value = "No selection!"
text_box.value = value
points_all = pn.widgets.TextAreaInput(value="", height=200, placeholder="Selected indices will appear here") # type: ignore[no-untyped-call]
selection.add_subscriber(lambda index: print_all_points(df=global_vosphoros_points, indices=index, text_box=points_all))
plot_ = (black_sea_vosporos_ * global_vosphoros_ * points_).opts(width = 700, height=800)
layout = pn.Column(
points_.opts(width=900, height=500),
pn.Row(
pn.Column("## Selected Indices:", points_all),
),
)
out = pn.panel(layout).servable()
out
BokehModel(combine_events=True, render_bundle={'docs_json': {'31c4f898-81d6-43b9-9888-6681c9da68bd': {'version…
Out[45]:
we'll fix the following points and regenerate the meshing for the black sea
In [46]:
pts_ind = [6817898,
6697291,
6549074,
1714532,
1473386]
coords = mesh_dic['nodes'][pts_ind][:,:2]
gdf = gp.GeoDataFrame(
pd.DataFrame(coords, columns=["lon", "lat"]),
geometry=[Point(xy) for xy in coords],
crs="EPSG:4326"
)
gdf['physical'] = "force"
gdf
gdf.to_file("points_4326.shp", driver='ESRI Shapefile')
Out[46]:
| lon | lat | geometry | physical | |
|---|---|---|---|---|
| 0 | 26.249899 | 40.010687 | POINT (26.2499 40.01069) | force |
| 1 | 26.245730 | 40.021553 | POINT (26.24573 40.02155) | force |
| 2 | 26.238869 | 40.034125 | POINT (26.23887 40.03413) | force |
| 3 | 26.230883 | 40.045579 | POINT (26.23088 40.04558) | force |
| 4 | 26.262392 | 40.001939 | POINT (26.26239 40.00194) | force |
all these points need to be suppressed
In [47]:
global_nodes = np.vstack((x,y)).T
In [48]:
discard = [1807639,456464,4792637,2451328,39451,5503929,4015013,6409047,2134299,4386548,4967757,
4856350,7851069,1324320,3643482,7452558,512577,6384066,677206,2528912,7294189,
167240,5582097,7382856,4192983,7134380,1647070,7700285,7907239,2334291,3421735,363349,
7046390,5460790,1145947,1314347,2426245,5698280,3282119,4838108,332392,6336320,1580434,
4864137,2614067,263530,1928581,1756355,4590953,47220,4869691,3173761,7342931,551330,189622,
2312122,7482244,5405589,4401906]
nno, ell = drop(global_nodes, tris, discard)
In [49]:
llo, lla, tri_indices, node_indices, tri_sub, node_mapping = extract_area(nno[:,0],nno[:,1],ell, (25,45), (40,48))
global_black_sea_gdf = tri_to_geopandas(np.vstack((llo,lla)).T, tri_indices, node_indices)
global_vosphoros_gdf = global_black_sea_gdf[global_black_sea_gdf.geometry.within(VOSPOROS_SOUTH)]
global_vosphoros_ = render_gdf(global_vosphoros_gdf)
In [50]:
# now create points from the black sea mesh
points2_ = render_gdf(black_sea_vosporos_points, "black sea")
In [51]:
selection = hv.streams.Selection1D(source=points2_)
points_all = pn.widgets.TextAreaInput(value="", height=200, placeholder="Selected indices will appear here") # type: ignore[no-untyped-call]
def print_all_points(df: pd.DataFrame, indices: pd.Index[int], text_box: pn.widgets.TextAreaInput) -> T.Any:
print("Selected indices:", indices)
if indices:
ind = np.unique(indices)
list = df["index"].iloc[ind].astype(int).unique()
value = ",\n".join(map(str, list))
else:
value = "No selection!"
text_box.value = value
selection.add_subscriber(lambda index: print_all_points(df=black_sea_vosporos_points, indices=index, text_box=points_all))
plot_ = (black_sea_vosporos_ * global_vosphoros_ * points_ * points2_).opts(width = 700, height=800)
layout = pn.Column(
plot_.opts(width=900, height=500),
pn.Row(
pn.Column("## Selected Indices:", points_all),
),
)
out = pn.panel(layout).servable()
out
BokehModel(combine_events=True, render_bundle={'docs_json': {'0914cea4-577c-4250-bbb1-18b5d9c2add2': {'version…
Out[51]:
here are the following node to drop on the back sea mesh
In [52]:
discard2 = [13598,39688,31622,28962,28963,528,529,522,28379,17111,29402,37351,23225,17568,45282,46776,
523,46587,31245,46194,48307,524,47748,44,48416,519,518,42,520,525,521,527,526,35710,43,149, 145]
new_bs_pts, new_bs_tri = drop(bs_pts,bs_tri, discard2)
In [53]:
new_black_sea_gdf = tri_to_geopandas(new_bs_pts, new_bs_tri, np.arange(len(new_bs_pts)))
new_bs_vosphoros_gdf = new_black_sea_gdf[new_black_sea_gdf.geometry.within(VOSPOROS_SOUTH)]
new_black_sea_vosporos_points = polygon_to_points(new_bs_vosphoros_gdf)
new_points2_ = render_gdf(new_black_sea_vosporos_points, "black sea")
In [54]:
plot_ = (render_gdf(global_vosphoros_gdf) *
render_gdf(new_bs_vosphoros_gdf, case = "black sea") * points_ * new_points2_).opts(width = 700, height=800)
In [55]:
plot_
Out[55]:
In [56]:
stitching_map = dict([
[512, 7732209],
[144, 1714532],
[143, 6549074],
[142, 6697291],
[37476, 3418943],
[22941, 1473386]
])
In [57]:
import numpy as np
def stitch_meshes(global_nodes, global_elements, bs_pts, bs_tri, stitching_map):
mapped_bs_indices = set(stitching_map.keys())
unmapped_bs_indices = [i for i in range(len(bs_pts)) if i not in mapped_bs_indices]
index_mapping = {}
new_nodes = global_nodes.tolist()
for i in unmapped_bs_indices:
new_index = len(new_nodes)
index_mapping[i] = new_index
# Insert with dummy z-value (can be interpolated or set to 0)
new_nodes.append([bs_pts[i][0], bs_pts[i][1], 0.0])
for bs_idx, global_idx in stitching_map.items():
index_mapping[bs_idx] = global_idx
new_elements = []
for tri in bs_tri:
new_tri = [index_mapping[idx] for idx in tri]
new_elements.append(new_tri)
all_elements = np.vstack([global_elements, np.array(new_elements, dtype=int)])
return np.array(new_nodes), all_elements
In [58]:
stiched_nodes, stiched_tri = stitch_meshes(mesh_dic['nodes'], mesh_dic['elements'], new_bs_pts, new_bs_tri, stitching_map)
stiched_nodes_fix, stiched_tri_fix = drop(stiched_nodes, stiched_tri, discard)
In [59]:
%matplotlib widget
In [60]:
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize = (10,6))
ax.triplot(stiched_nodes_fix[:,0], stiched_nodes_fix[:,1], stiched_tri_fix, lw=0.2)
ax.set_xlim([26,30])
ax.set_ylim([39.75,41.5])
Out[60]:
[<matplotlib.lines.Line2D at 0x70581ba8a990>, <matplotlib.lines.Line2D at 0x70581ba8a2d0>]
Out[60]:
(26.0, 30.0)
Out[60]:
(39.75, 41.5)
In [61]:
import meshio
points = np.column_stack((stiched_nodes_fix[:,0], stiched_nodes_fix[:,1], np.zeros(len(stiched_nodes_fix)))) # Add z=0
cells = [("triangle", stiched_tri_fix)]
mesh = meshio.Mesh(points=points, cells=cells)
mesh.write("output_mesh.vtk")
In [62]:
from xarray_selafin.xarray_backend import SelafinAccessor
import xarray as xr
ds = xr.Dataset({
"B": (("time", "node"), np.zeros((1,len(stiched_nodes_fix)))),
},
coords = {
"x" : ("node", stiched_nodes_fix[:,0]),
"y" : ("node", stiched_nodes_fix[:,1]),
"time": [pd.Timestamp.now()]
} )
ds.attrs['ikle2'] = stiched_tri_fix + 1
ds
Out[62]:
<xarray.Dataset> Size: 193MB
Dimensions: (time: 1, node: 8039622)
Coordinates:
x (node) float64 64MB -1.011 -27.3 -92.7 104.2 ... 32.13 36.43 29.07
y (node) float64 64MB 45.49 38.64 18.62 1.046 ... 46.11 41.38 40.94
* time (time) datetime64[ns] 8B 2025-05-13T12:27:34.032880
Dimensions without coordinates: node
Data variables:
B (time, node) float64 64MB 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
Attributes:
ikle2: [[3987707 1 7764558]\n [ 1 4523085 7764558]\n [30905...In [ ]:
ds.selafin.write('out.slf')