Adding the mesh refinement based on polygon and adding the config.ini…

…, addressing issue NOAA-EMC#55

unst_msh_gen/README.md

@@ -3,7 +3,7 @@ Unstructured mesh generation for WW3 using JIGSSAW.
 # Description
 Mesh generation script capable of creating unstructured meshes for WW3 global modeling. The tool leverages JIGSAWPY (https://github.com/dengwirda/jigsaw-python) for efficient triangulation.
 Main changes include:
-Implementation of the ocn_ww3.py to create uniform unstructured mesh for global WW3 model.
+Implementation of the ocn_ww3.py to create uniform and variable unstructured mesh for global WW3 model.
 This tool is under active development, with future work focused on variable unstructured mesh generation.
 
 # Installation
@@ -33,35 +33,44 @@ This tool is under active development, with future work focused on variable unst
 
 # Usage
 ## 5- run the script inside of WW3-tools/unst_msh_gen:
-- $python3 ocn_ww3.py --black_sea [option]
-
-		option=  3: default which will have the Black Sea and the connections.
-			 2: will have the Balck sea as a seperate basin.
-			 1: will exclude the Black sea
 
-NOTE: the output will be gmsh format which will be used by WW3.
+## modify config.init
+- Specify the DEM netcdf file using the 'dem_file' in "DataFiles" section.
+- For uniform resolution 'hmax' = 'hmin' = 'hshr' = 'hfun_max' and 'nwave' = 0 in "Spacing" section.
+- To include or exclude Black-Sea change 'black_sea' in "CommandLineArg" section to:
+- 	                  3: will have the Black Sea and the connections.
+-                         2: will have the Balck sea as a seperate basin.
+-                         1: will exclude the Black sea
+-
+- you can specify the mesh name (jigsaw format or ww3 format) using 'mesh_file' and 'ww3_mesh_file'in "MeshSetting" section 
+
+- $python3 ocn_ww3.py
 
-NOTE: for different resolution (uniform) in km, you should change the following:
-- opts.hfun_hmax
-- hmax
-- hshr
-- hmin
+NOTE: the output will be gmsh format which will be used by WW3 (specified by 'ww3_mesh_file')
 
 NOTE: The output mesh will have -180:180 longitude, you can convert this by unisg ShiftMesh.py script, to 0:360 longitude.
 	input_file_path: your jigsaw mesh in gmsh format with -18:180 long
 	output_file_path: shifted mesh in gmsh format with 0:360 long
 
 
 ## 6- Using variable mode:
-- To create a mesh with finer resolution near the US coastlines you can define different region in json format (east coast, west coast and golf od Mexico, Purto Rico, and Hawaii):
+
+## modify config.ini
+- To create a mesh with finer resolution near the US coastlines you can define different region in json format (east coast, west coast and golf od Mexico, Purto Rico, and Hawaii); use the 'window_file' in the "DataFiles" section to specify the windows or parse the windows as;
 
 - $python3 window_mask.py --windows '[{"min_lon": -98, "max_lon": -64, "min_lat": 24, "max_lat": 44.5, "hshr": 5}, {"min_lon": -158, "max_lon": -155, "min_lat": 19, "max_lat": 22, "hshr": 5}, {"min_lon": -128, "max_lon": -64, "min_lat": 34.5, "max_lat": 48.5, "hshr": 5}, {"min_lon": -67.4, "max_lon": -64.1, "min_lat": 17.2, "max_lat": 18.3, "hshr": 5}]'
 
 NOTE: hshr is the shoreline resolution which can be smaller than the hmin which is defined globally.
 
-NOTE: You can define different background mesh based on lat location in the window_mask.py:
+- window_mask.py can read shapefiles in json format as well and assign user defined resolution (scale) to each polygon; use the 'shape_file' in the "DataFiles" to specify the shapefiles or parse the shapefiles as;
+
+- $python3 window_mask.py --shapefiles '[{"path": "./Shapefiles/Arctic.shp", "scale": 10}, {"path": "./Shapefiles/Gulf.shp", "scale": 5}, {"path": "./Shapefiles/Hawaii.shp", "scale": 5}, {"path": "./Shapefiles/PurtoRico.shp", "scale": 10}, {"path": "./Shapefiles/WestCoast.shp", "scale": 5}]'
+
+NOTE: You can create multiple shapefiles using QGIS or ArcGIS and save them in ./Shapefiles directory and assign different resolution (scale) to each polygone.
+
+NOTE: You can define different background mesh based on lat location in the window_mask.py via config.init
 
-NOTE: for the background mesh you can define different resolution (for eaxmple, the globe is divided to three regions based on latitude and corresponding resolutions: 
+NOTE: for the background mesh you can define different resolution (for eaxmple, the globe is divided to three regions based on latitude and corresponding resolutions) in "Scalingsettings" of the config.ini as;
 
 - upper_bound = 50          # Lattitude that starts the upper section (ie lat > upper_bound) 
 - middle_bound = -20        # Determines the boundary between the upper  and middle sections.  The middle section is for upper_bound > lat > middle_bound
@@ -73,13 +82,18 @@ NOTE: for the background mesh you can define different resolution (for eaxmple,
 - scale_south_upper = 30    # km resolution for upper south/lower section 
 - scale_south_lower = 9     # km mesh resolution at lower south/lower section 
 
-NOTE: The output will be wmask.nc file which will have the mesh spacing info.
+NOTE: The output will be "wmask.nc" file which will have the mesh spacing info.
+
+
+NOTE:To create the variable mesh based on mesh spacing file "wmask.nc", in the config.ini change the mask_file = wmask.nc
 
+- $python3 ocn_ww3.py
 
-NOTE:To create the variable mesh based on specified mesh spacing file in "wmask.nc" you can use the following comand:
 
-- $python3 ocn_ww3.py --black_sea [option] --mask_file="wmask.nc"
+## 7- Plotting the mesh info:
 
+- to plot the elements, mesh size, and bathymetry data, use plot_msh.py and change the 'filename' to ww3 format mesh.
+- $python3 plot_msh.py
 
 ## Contributing
-This is ongoing effort with the great help of Darren Engwirda, JIGSAW developer.
+This is ongoing effort by Ali Salimi-Tarazouj with the great help of Darren Engwirda, JIGSAW developer.

unst_msh_gen/ShiftMesh.py

@@ -1,3 +1,13 @@
+"""
+To shift the mesh from -180:180 (lon) to 0:360
+
+input_file_path: -180:180 (lon) mesh in GMSH fomrat
+output_fil_pat: 0:360 (lon) mesh in GMSH format 
+"""
+
+# Author: Ali Salimi-Tarazouj
+
+
 input_file_path = './9km_nobc.ww3'
 output_file_path = './modified_9km_nobc.ww3'
 

unst_msh_gen/config.ini

@@ -0,0 +1,56 @@
+; Spacing configurations for the mesh generation
+[Spacing]
+hmax = 100.0 
+; Maximum spacing in kilometers
+hshr = 100 
+; Spacing at the shoreline in kilometers
+nwav = 400 
+; Number of cells per square root of (gravity * water depth)
+hmin = 100.0 
+; Minimum spacing in kilometers
+dhdx = 0.05 
+; Maximum allowable gradient in spacing
+
+; Scaling settings based on latitude for mesh refinement
+[ScalingSettings]
+upper_bound = 50 
+; Northern latitude above which scale_north is applied
+middle_bound = -20 
+; Middle latitude boundary
+lower_bound = -90 
+; Southern latitude below which scale_south_lower is applied
+scale_north = 9 
+; Scaling factor for northern region
+scale_middle = 20 
+; Scaling factor for middle region
+scale_south_upper = 30 
+; Upper scaling factor for southern region
+scale_south_lower = 9 
+; Lower scaling factor as it approaches the lower bound
+
+; Mesh generation settings
+[MeshSettings]
+hfun_hmax = 100 
+; Global maximum mesh resolution, similar to hmax
+mesh_file = uglo_poly_nBlkS.msh 
+; Filename for the generated mesh in Jigsaw format
+ww3_mesh_file = uglo_poly_nBlkS.ww3 
+; Filename for the mesh in WW3 format
+
+; Command-line arguments and their default values
+[CommandLineArgs]
+black_sea = 3 
+; Black Sea configuration mode: 1=no Black Sea, 2=detached, 3=connected
+mask_file = wmask.nc 
+; Path to the weighting mask file, if any
+
+; Bathymetry and Regional data
+[DataFiles]
+dem_file = /scratch2/NCEPDEV/marine/alisalimi/jigsaw-geo-tutorial/VAR_MSH/CODE/RTopo_2_0_4_GEBCO_v2023_60sec_pixel.nc
+; DEM file
+
+shape_file = [{"path": "/scratch1/NCEPDEV/stmp2/Ali.Salimi/msh-gen/unst_msh_gen/VAR-MESH/Shapefiles/Arctic.shp", "scale": 10}, {"path": "/scratch1/NCEPDEV/stmp2/Ali.Salimi/msh-gen/unst_msh_gen/VAR-MESH/Shapefiles/Gulf.shp", "scale": 5}, {"path": "/scratch1/NCEPDEV/stmp2/Ali.Salimi/msh-gen/unst_msh_gen/VAR-MESH/Shapefiles/Hawaii.shp", "scale": 5}, {"path": "/scratch1/NCEPDEV/stmp2/Ali.Salimi/msh-gen/unst_msh_gen/VAR-MESH/Shapefiles/PurtoRico.shp", "scale": 10}, {"path": "/scratch1/NCEPDEV/stmp2/Ali.Salimi/msh-gen/unst_msh_gen/VAR-MESH/Shapefiles/WestCoast.shp", "scale": 5}]
+; shapefiles for regional mesh tefinement
+
+;window_file = [{"min_lon": -98, "max_lon": -64, "min_lat": 24, "max_lat": 44.5, "hshr": 5}, {"min_lon": -158, "max_lon": -155, "min_lat": 19, "max_lat": 22, "hshr": 5}, {"min_lon": -128, "max_lon": -64, "min_lat": 34.5, "max_lat": 48.5, "hshr": 5}, {"min_lon": -67.4, "max_lon": -64.1, "min_lat": 17.2, "max_lat": 18.3, "hshr": 5}]
+; user-defined windows to refine the mesh regionally.

unst_msh_gen/ocn_ww3.py

@@ -2,11 +2,13 @@
 """ Jigsaw meshes for WW3 with global bathymetry
 """
 
-# Authors: Ali Salimi, Darren
+# Authors: Ali Salimi-Tarazouj, Darren Engwirda
 
 # The DEM file used below can be found at:	
 # https://github.com/dengwirda/dem/releases/tag/v0.1.1
+import os
 import argparse
+import configparser
 import numpy as np
 import netCDF4 as nc
 
@@ -18,18 +20,21 @@
 
 from spacing import *
 
+# Load the configuration file
+config = configparser.ConfigParser()
+config.read('config.ini')
+
 # Create the parser
 parser = argparse.ArgumentParser(description="Run mesh gen with specific region settings.")
 
-# Add an argument for the black_sea option
-parser.add_argument('--black_sea', type=int, default=3,
+# Command-line arguments with defaults from config file
+parser.add_argument('--black_sea', type=int, default=config.getint('CommandLineArgs', 'black_sea'),
                     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='',
+parser.add_argument('--mask_file', type=str, default=config.get('CommandLineArgs', 'mask_file', fallback=''),
                     help='Path to the mask file, if any. Leave empty if not used.')
 
+
 # Parse the command-line arguments
 args = parser.parse_args()
 
@@ -68,10 +73,10 @@ def create_msh():
     # solve |dh/dx| constraints in spacing
     jigsawpy.cmd.marche(opts, spac)
 
-    opts.mesh_file = "uglo_100km_BlkS.msh"  #jigsaw format mesh file
+    opts.mesh_file = config['MeshSettings']['mesh_file']  #jigsaw format mesh file
 
     opts.hfun_scal = "absolute"
-    opts.hfun_hmax = 100.           # global maximum mesh resolution (similar to hmax)
+    opts.hfun_hmax = float(config['MeshSettings']['hfun_hmax'])           # global maximum mesh resolution (similar to hmax)
     opts.mesh_dims = +2             # 2-dim. simplexes
     opts.optm_iter = +64            # number of itereation for the optimization
     opts.optm_cost = "skew-cos"
@@ -83,14 +88,16 @@ def create_siz():
 
     #-- create mesh spacing function for the globe: for uniform mesh hmax = hshr = hmin
 
-    hmax = 100.0 # maximum spacing [km] 
-    hshr = 100  # shoreline spacing
-    nwav = 400.  # number of cells per sqrt(g*H)
-    hmin = 100.0  # minimum spacing
-    dhdx = 0.05  # allowable spacing gradient: for more gradual transition use lower value
+    hmax = float(config['Spacing']['hmax']) # maximum spacing [km] 
+    hshr = float(config['Spacing']['hshr'])  # shoreline spacing
+    nwav = int(config['Spacing']['nwav'])  # number of cells per sqrt(g*H)
+    hmin = float(config['Spacing']['hmin'])  # minimum spacing
+    dhdx = float(config['Spacing']['dhdx'])  # allowable spacing gradient: for more gradual transition use lower value
 
-    data = nc.Dataset(
-        "RTopo_2_0_4_GEBCO_v2023_60sec_pixel.nc", "r")
+    # Load the DEM file from the config
+    dem_file = config['DataFiles']['dem_file']
+
+    data = nc.Dataset(dem_file,"r")
 
     xlon = np.asarray(data["lon"][:])
     ylat = np.asarray(data["lat"][:])
@@ -176,8 +183,10 @@ def inject_dem():
 
     print("*inject-dem...")
 
-    data = nc.Dataset(
-        "RTopo_2_0_4_GEBCO_v2023_60sec_pixel.nc", "r")
+        # Load the DEM file from the config
+    dem_file = config['DataFiles']['dem_file']
+
+    data = nc.Dataset(dem_file,"r")
 
     xlon = np.asarray(data["lon"][:])
     ylat = np.asarray(data["lat"][:])
@@ -531,5 +540,6 @@ def write_gmsh_mesh(filename, node_data, tri):
     point = np.hstack((point, depth))  # append elev. as 3rd coord.
     cells = [("triangle", mesh.tria3["index"])]
     tri_data=cells[0][1]+1
-    write_gmsh_mesh("uglo_100km_BlkS.ww3", point, tri_data)
+    ww3_mesh_file = config['MeshSettings']['ww3_mesh_file']
+    write_gmsh_mesh(ww3_mesh_file, point, tri_data)