Merge pull request #144 from GeoOcean/WMS

Wms

Author: FaugereE

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
 

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]")