From 726a6255b39f94ab8b32e716cbd0fe5d30a58cf3 Mon Sep 17 00:00:00 2001 From: FaugereE Date: Mon, 16 Feb 2026 19:37:47 +0100 Subject: [PATCH 1/2] [EJP] greensurge fix --- bluemath_tk/additive/greensurge.py | 83 ++++++++++++++++++++++++++++-- 1 file changed, 80 insertions(+), 3 deletions(-) diff --git a/bluemath_tk/additive/greensurge.py b/bluemath_tk/additive/greensurge.py index e034924..9041d42 100644 --- a/bluemath_tk/additive/greensurge.py +++ b/bluemath_tk/additive/greensurge.py @@ -13,12 +13,92 @@ from matplotlib.figure import Figure from matplotlib.path import Path from tqdm import tqdm +from shapely.geometry import LineString, Point from ..core.operations import get_degrees_from_uv from ..core.plotting.colors import hex_colors_land, hex_colors_water from ..core.plotting.utils import join_colormaps +def add_forcing_edges_to_computational_domain(poly_clip, vert, tria): + """ + Add edges to the forcing mesh that align with the computational domain boundary. + Parameters + ---------- + poly_clip : shapely.geometry.Polygon + The polygon representing the computational domain boundary. + vert : np.ndarray + Array of shape (n_vertices, 2) containing the coordinates of the vertices in the forcing mesh. + tria : np.ndarray + Array of shape (n_triangles, 3) containing the indices of the vertices for each triangle in the forcing mesh. + Returns + ------- + vert_clean : np.ndarray + Updated array of vertices including new intersection points. + edges_clean : np.ndarray + Updated array of edges including new edges along the computational domain boundary. + """ + + # Extract all edges from forcing mesh triangles + edges_new_full = np.vstack([ + tria[:, [0, 1]], + tria[:, [1, 2]], + tria[:, [2, 0]] + ]) + edges_new = np.unique(np.sort(edges_new_full, axis=1), axis=0) + + # Check which vertices are inside computational domain + mask_node_in = np.array([poly_clip.contains(Point(x, y)) for x, y in vert]) + + # Classify edges by how many endpoints are inside + count_inside = mask_node_in[edges_new].sum(axis=1) + edges_both_inside = edges_new[count_inside == 2] + edges_one_inside = edges_new[count_inside == 1] + + # Clip edges at computational domain boundary + edges_one_inside_new = edges_one_inside.copy() + new_points = [] + + for i, (n1, n2) in enumerate(edges_one_inside): + p1, p2 = vert[n1], vert[n2] + inside1 = mask_node_in[n1] + + line = LineString([p1, p2]) + inter = line.intersection(poly_clip.boundary) + + if inter.is_empty: + continue + + # Handle multiple intersection points + if inter.geom_type == "MultiPoint": + points = list(inter.geoms) + ref = p1 if inside1 else p2 + dists = [Point(ref).distance(pt) for pt in points] + inter_pt = np.array(points[np.argmin(dists)].coords[0]) + else: + inter_pt = np.array(inter.coords[0]) + + new_index = len(vert) + len(new_points) + new_points.append(inter_pt) + + if inside1: + edges_one_inside_new[i, 1] = new_index + else: + edges_one_inside_new[i, 0] = new_index + + # Update vertices with new intersection points + if len(new_points) > 0: + vert_update = np.vstack([vert, np.array(new_points)]) + + # Clean up: keep only used nodes + edges_used = np.vstack([edges_both_inside, edges_one_inside_new]) + used_nodes = np.unique(edges_used.flatten()) + vert_clean = vert_update[used_nodes] + mapping = {old_idx: new_idx for new_idx, old_idx in enumerate(used_nodes)} + edges_clean = np.vectorize(mapping.get)(edges_used) + + return vert_clean, edges_clean + def read_adcirc_grd(grd_file: str) -> Tuple[np.ndarray, np.ndarray, List[str]]: """ Reads the ADCIRC grid file and returns the node and element data. @@ -376,9 +456,6 @@ def plot_greensurge_setup( ax.set_extent([*bnd], crs=ccrs.PlateCarree()) plt.legend(loc="lower left", fontsize=10, markerscale=2.0) ax.set_title("GreenSurge Mesh Setup") - gl = ax.gridlines(draw_labels=True) - gl.top_labels = False - gl.right_labels = False return fig, ax From d271913a8350c0f12818c9f94d37e4024acd993a Mon Sep 17 00:00:00 2001 From: FaugereE Date: Mon, 16 Feb 2026 19:38:09 +0100 Subject: [PATCH 2/2] [EJP] vortex interpolation fix --- bluemath_tk/tcs/vortex.py | 51 ++++++++++++++++++++++++--------------- 1 file changed, 32 insertions(+), 19 deletions(-) diff --git a/bluemath_tk/tcs/vortex.py b/bluemath_tk/tcs/vortex.py index f256d6c..d7ed454 100644 --- a/bluemath_tk/tcs/vortex.py +++ b/bluemath_tk/tcs/vortex.py @@ -216,6 +216,7 @@ def vortex_model_grid( def vortex2delft_3D_FM_nc( mesh: xr.Dataset, ds_vortex: xr.Dataset, + fill_value: float = 0, ) -> xr.Dataset: """ Convert the vortex dataset to a Delft3D FM compatible netCDF forcing file. @@ -232,6 +233,8 @@ def vortex2delft_3D_FM_nc( The name of the output netCDF file, default is "forcing_Tonga_vortex.nc". forcing_ext : str The extension for the forcing file, default is "GreenSurge_GFDcase_wind.ext". + fill_value : float + The fill value to use for missing data, default is 0. Returns ------- xarray.Dataset @@ -245,33 +248,43 @@ def vortex2delft_3D_FM_nc( lat_interp = xr.DataArray(latitude, dims="node") lon_interp = xr.DataArray(longitude, dims="node") - angle = np.deg2rad((270 - ds_vortex.Dir.values) % 360) - W = ds_vortex.W.values + W_node = ds_vortex.W.interp( + lat=lat_interp, + lon=lon_interp, + method="nearest", + ) + + Dir_node = ds_vortex.Dir.interp( + lat=lat_interp, + lon=lon_interp, + method="nearest", + ) + + p_node = ds_vortex.p.interp( + lat=lat_interp, + lon=lon_interp, + method="nearest", + ) + + angle = np.deg2rad((270 - Dir_node.values) % 360) - ds_vortex_interp = xr.Dataset( + windx_node = (W_node.values * np.cos(angle)).astype(np.float32) + windy_node = (W_node.values * np.sin(angle)).astype(np.float32) + + forcing_dataset = xr.Dataset( { - "windx": ( - ("latitude", "longitude", "time"), - (W * np.cos(angle)).astype(np.float32), - ), - "windy": ( - ("latitude", "longitude", "time"), - (W * np.sin(angle)).astype(np.float32), - ), - "airpressure": ( - ("latitude", "longitude", "time"), - ds_vortex.p.values.astype(np.float32), - ), + "windx": (("node", "time"), np.nan_to_num(windx_node, nan=fill_value)), + "windy": (("node", "time"), np.nan_to_num(windy_node, nan=fill_value)), + "airpressure": (("node", "time"), np.nan_to_num(p_node.values.astype(np.float32), nan=fill_value)), }, coords={ - "latitude": ds_vortex.lat.values, - "longitude": ds_vortex.lon.values, + "node": np.arange(lat_interp.size), "time": np.arange(n_time), + "latitude": ("node", latitude), + "longitude": ("node", longitude), }, ) - forcing_dataset = ds_vortex_interp.interp(latitude=lat_interp, longitude=lon_interp) - reference_date_str = ( ds_vortex.time.values[0] .astype("M8[ms]")