Adding window_mask.py for variable unstructured mesh, and revising th…

…e spacing.py and ocn_ww3.py to read the mask file, addressing issue NOAA-EMC#55

unst_msh_gen/README.md

@@ -45,5 +45,27 @@ NOTE: The output mesh will have -180:180 longitude, you can convert this by unis
 	output_file_path: shifted mesh in gmsh format with 0:360 long
 
 
+6- Using variable mode:
+	$python3 window_mask.py --windows '[{"min_lon": -98, "max_lon": -64, "min_lat": 24, "max_lat": 44.5, "hshr": 5}'
+
+NOTE: You can define multiple windows with different shoreline resolution (hshr) in json format.
+NOTE: You can define different background mesh based on lat location in the window_mask.py:
+
+	 # Define the boundaries and scaling values
+   	 upper_bound = 50
+    	 middle_bound = -20
+    	 lower_bound = -90
+    	 scale_north = 9
+    	 scale_middle = 20
+    	 scale_south_upper = 30
+    	 scale_south_lower = 9
+
+NOTE: The output will be wmask.nc file which will have the mesh spacing info.
+
+
+        $python3 ocn_ww3.py --black_sea [option] --mask_file="wmask.nc"
+
+NOTE: This will create the variable mesh based on specified mesh spacing file.
+
 ## Contributing
 This is ongoing effort with the great help of Darren Engwirda, JIGSAW developer.

unst_msh_gen/ocn_ww3.py

@@ -25,6 +25,11 @@
 parser.add_argument('--black_sea', type=int, default=3,
                     help='Set the region mode: 1 no black-sea, 2 black-sea detached, 3 black-sea with connections. Default is 3.')
 
+
+# Add an argument for the mask_file
+parser.add_argument('--mask_file', type=str, default='',
+                    help='Path to the mask file, if any. Leave empty if not used.')
+
 # Parse the command-line arguments
 args = parser.parse_args()
 
@@ -63,10 +68,10 @@ def create_msh():
     # solve |dh/dx| constraints in spacing
     jigsawpy.cmd.marche(opts, spac)
 
-    opts.mesh_file = "mesh_50km_nobc.msh"
+    opts.mesh_file = "uglo_100km_BlkS.msh"
 
     opts.hfun_scal = "absolute"
-    opts.hfun_hmax = +50.           # uniform at 30.km
+    opts.hfun_hmax = 100.           # uniform at 30.km
     opts.mesh_dims = +2             # 2-dim. simplexes
     opts.optm_iter = +64
     opts.optm_cost = "skew-cos"
@@ -78,11 +83,11 @@ def create_siz():
 
 #-- create mesh spacing function for the globe
 
-    hmax = 50.0 # maximum spacing [km]
-    hshr = 50  # shoreline spacing
+    hmax = 100.0 # maximum spacing [km]
+    hshr = 100  # shoreline spacing
     nwav = 400.  # number of cells per sqrt(g*H)
-    hmin = 50.0  # minimum spacing
-    dhdx = 0.1  # allowable spacing gradient
+    hmin = 100.0  # minimum spacing
+    dhdx = 0.05  # allowable spacing gradient
 
     data = nc.Dataset(
         "RTopo_2_0_4_GEBCO_v2023_60sec_pixel.nc", "r")
@@ -112,7 +117,15 @@ def create_siz():
     hmat[high] = hmax
 
     hmat = setup_shoreline_pixels(hmat, land, hshr)
-
+
+#-- apply user-defined scaling: multiply h(x) by mask array
+
+    if hasattr(args, 'mask_file') and args.mask_file:
+       hmat = scale_spacing_via_mask(args, hmat)
+    else:
+    # Handle case where mask_file is not provided
+        print("No mask file provided. Proceeding without scaling...")
+
 #-- and a little nonlinear smoothing
 
     filt = filter_pixels_harmonic(hmat, exp=2)
@@ -139,12 +152,8 @@ def create_siz():
 
     ymid = 41.5 * np.pi / 180.
     xmid = 30.5 * np.pi / 180.
-
-    zoom = +100.0 - 99.0 * np.exp(-(
-        6.75 * (xmat - xmid) ** 2 +
-        12.5 * (ymat - ymid) ** 2) ** 2)
 
-    spac.value = hmat * zoom
+    spac.value = hmat 
     spac.slope = np.array(dhdx)
     spac.value = np.minimum(hmax, spac.value)
 
@@ -518,9 +527,9 @@ def write_gmsh_mesh(filename, node_data, tri):
     point = jigsawpy.R3toS2(geom.radii, point)  # to [lon,lat] in deg
     point*= 180. / np.pi
     depth = np.reshape(-1*mesh.value, (mesh.value.size, 1))
-    #depth[depth < 0] = 50
+    depth[depth <= 0] = 2
     point = np.hstack((point, depth))  # append elev. as 3rd coord.
     cells = [("triangle", mesh.tria3["index"])]
     tri_data=cells[0][1]+1
-    write_gmsh_mesh("50km_nobc.ww3", point, tri_data)
+    write_gmsh_mesh("uglo_100km_BlkS.ww3", point, tri_data)
 

unst_msh_gen/spacing.py

@@ -2,13 +2,14 @@
 """ Mesh spacing utilities, for barotropic flows
 """
 
-# Authors: Darren Engwirda
+# Authors: Darren Engwirda, Ali Salimi
 
 # routines to compute mesh spacing functions that scale
 # with shallow-water wave lengths, elev. gradients, etc 
 
 import numpy as np
 import jigsawpy
+import netCDF4 as nc
 
 from skimage.filters import gaussian, median
 from skimage.measure import label, regionprops_table
@@ -174,4 +175,14 @@ def elev_sharpness_spacing(
     vals = np.asarray(vals, dtype=np.float32)
 
     return vals
-
+
+
+def scale_spacing_via_mask(args, vals):
+    print("User-defined h(x) scaling...")
+    if hasattr(args, 'mask_file') and args.mask_file:  # Check if mask_file exists and is not empty
+        data = nc.Dataset(args.mask_file, "r")
+        vals = np.asarray(data["val"][:], dtype=np.float32)
+    else:
+        print("No mask file provided. Scaling will not be applied.")
+    return vals
+

unst_msh_gen/window_mask.py

@@ -0,0 +1,83 @@
+""" Jigsaw meshes for WW3 with global bathymetry
+"""
+
+# Authors: Ali Salimi, Darren
+
+# This script will create mesh spacing based on user defined windows in json format.
+
+import numpy as np
+import netCDF4 as nc
+import argparse
+import json
+
+def create_mask_file(windows):
+    # Open the data file
+    data = nc.Dataset("RTopo_2_0_4_GEBCO_v2023_60sec_pixel.nc", "r")
+
+    # Extract longitude and latitude arrays
+    xlon = np.asarray(data["lon"][:])
+    ylat = np.asarray(data["lat"][:])
+    elev = np.asarray(data["bed_elevation"][:], dtype=np.float32) + np.asarray(data["ice_thickness"][:], dtype=np.float32)
+
+    # Compute midpoints for longitude and latitude
+    xmid = 0.5 * (xlon[:-1] + xlon[+1:])
+    ymid = 0.5 * (ylat[:-1] + ylat[+1:])
+
+    # Create a meshgrid of midpoints
+    xmat, ymat = np.meshgrid(xmid, ymid)
+
+    # Initialize scal with a default condition
+    #scal = np.where(ymat > 50, 9, np.where(ymat < -20, 30, 20))
+
+    # Define the boundaries and scaling values
+    upper_bound = 50
+    middle_bound = -20
+    lower_bound = -90
+    scale_north = 9
+    scale_middle = 20
+    scale_south_upper = 30
+    scale_south_lower = 9
+
+    # Calculate the scaling using conditions
+    scal = np.where(ymat > upper_bound,
+                    scale_north,  # Use 9 km for ymat > 50
+                    np.where(ymat > middle_bound,
+                             scale_middle,  # Use 20 km for -20 < ymat <= 50
+                             # Gradual decrease from 30 km to 9 km as latitude decreases from -20 to -90
+                             scale_south_upper + (scale_south_lower - scale_south_upper) * 
+                             (ymat - middle_bound) / (lower_bound - middle_bound)
+                            ))
+
+    # Process each window
+    for window in windows:
+        # Apply the 'hshr' value for the current window where conditions are met
+        window_mask_lon = np.logical_and(xmat >= window["min_lon"], xmat <= window["max_lon"])
+        window_mask_lat = np.logical_and(ymat >= window["min_lat"], ymat <= window["max_lat"])
+        window_mask = np.logical_and(window_mask_lon, window_mask_lat)
+
+        mask_elevation = np.logical_and(elev >= -4., elev <= +8.)
+        mask_final = np.logical_and(mask_elevation, window_mask)
+        scal[mask_final] = window["hshr"]
+
+    # Create a new NetCDF file to store the mask
+    data_out = nc.Dataset("wmask.nc", "w")
+    data_out.createDimension("nlon", xmid.size)
+    data_out.createDimension("nlat", ymid.size)
+    if "val" not in data_out.variables.keys():
+        data_out.createVariable("val", "f4", ("nlat", "nlon"))
+    data_out["val"][:, :] = scal
+    data_out.close()
+
+def parse_windows_from_args():
+    parser = argparse.ArgumentParser(description='Create a mask file with multiple windows.')
+    parser.add_argument('--windows', type=str, required=True,
+                        help='JSON string with window definitions. Example: \'[{"min_lon": -130, "max_lon": -64, "min_lat": 24.5, "max_lat": 47.5, "hshr": 5}]\'')
+
+    args = parser.parse_args()
+    windows = json.loads(args.windows)
+    return windows
+
+if __name__ == "__main__":
+    windows = parse_windows_from_args()
+    create_mask_file(windows)
+