STOFS3D scripts: added a Utils folder, starting to moving some post-p…

…rocessing functionalities here as a reference for new users to develop their tools. Fixed a bug in generate_adcirc.py, which has not been incorporated in operational forecast. Other minor edits.

src/Utility/Pre-Processing/STOFS-3D-Atl-operation/Utils/post_process_lib.py

@@ -0,0 +1,154 @@
+'''
+Sample functions for post-processing SCHISM output
+'''
+
+from copy import deepcopy
+import numpy as np
+import xarray as xr
+
+
+def split_quads(elements=None):
+    '''
+    Split quad elements to triangles and
+    append additional elements to the end of the element table
+
+    elements is read from SCHISM output file as:
+        elements = ds['SCHISM_hgrid_face_nodes'][:]
+    '''
+
+    if elements is None:
+        raise ValueError('No elements provided')
+
+    if elements.shape[1] == 3:  # already triangles
+        return elements
+    elif elements.shape[1] != 4:
+        raise ValueError('elements should be a numpy array of (n,3) or (n,4)')
+
+    triangles = deepcopy(elements)
+    quad_idx = ~elements[:, -1].mask
+
+    # split quads into two triangles
+    quads = elements[quad_idx]
+    upper_triangle = np.c_[  # last node is -1 (not applicable)
+        quads[:, 0], quads[:, 1], quads[:, 3], -np.ones((quads.shape[0], 1))]
+    lower_triangle = np.c_[  # last node is -1 (not applicable)
+        quads[:, 1], quads[:, 2], quads[:, 3], -np.ones((quads.shape[0], 1))]
+
+    # replace quads with upper triangle
+    triangles[quad_idx, :] = upper_triangle
+    # append lower triangle a the end
+    triangles = np.ma.concatenate([triangles, lower_triangle], axis=0)
+    # mask the last node, because all quads have been changed to triangles
+    triangles.mask[:, -1] = True
+
+    return triangles[:, :3]  # only return the first 3 nodes of each element
+
+
+def cal_maxelev(depth: np.ndarray, elev: np.ndarray, fill_value: float = -99999.0):
+    """
+    Calculate the maximum elevation and the time of the maximum elevation
+
+    - inputs:
+        elev: np.ndarray of shape (ntimes, npoints)
+        depth: np.ndarray of shape (npoints,)
+        fill_value: the value to fill the masked values in elev
+    - outputs:
+        maxelev: np.ndarray of shape (npoints,)
+        time_idx_maxelev: np.ndarray of shape (npoints,)
+    """
+    if fill_value > -1000:
+        raise ValueError("fill_value should be a large negative number")
+
+    # mask large values
+    elev[np.where(elev > 100000)] = fill_value
+    # mask dry nodes
+    elev[elev + depth <= 1e-6] = fill_value  # native schism outputs
+    elev[np.isnan(elev)] = fill_value  # deprecated, adcirc format
+
+    maxelev = np.max(elev, axis=0)
+    time_idx_maxelev = np.argmax(elev, axis=0)
+
+    return maxelev, time_idx_maxelev
+
+
+def cal_disturbance(
+    depth: np.ndarray, elev: np.ndarray,
+    city_node_idx_file: str = None,
+    fillvalue: float = -99999.0
+):
+    """
+    Calculate the maximum disturbance
+
+    - inputs:
+        depth: np.ndarray of shape (npoints,)
+        elev: np.ndarray of shape (ntimes, npoints)
+    - outputs:
+        maxdist: np.ndarray of shape (npoints,)
+    """
+    if fillvalue > -1000:
+        raise ValueError("fillvalue should be a large negative number")
+
+    elev[np.isnan(elev)] = fillvalue  # deprecated, adcirc format
+
+    disturb = elev.copy()  # same as elev in the ocean, so initialize with elev
+
+    # read handle city node indices
+    if city_node_idx_file is not None:
+        city_node_idx = np.loadtxt(city_node_idx_file, encoding='utf-8').astype(bool)
+    else:
+        city_node_idx = np.zeros_like(depth, dtype=bool)
+    # define land nodes, including city nodes
+    land_node_idx = (depth < 0) | city_node_idx
+
+    # Define disturbance:
+    # On land, disturbance is the sum of elevation and depth, i.e., the water depth.
+    # Also, disturbance is zero if the water depth is negative.
+    disturb[:, land_node_idx] = np.maximum(0, elev[:, land_node_idx] + depth[land_node_idx])
+
+    max_disturb = np.max(disturb, axis=0)
+    time_idx_max_disturb = np.argmax(disturb, axis=0)
+
+    # mask small max disturbance (< 0.3 m) on land (including cities)
+    small_dist_on_land = (max_disturb < 0.3) * land_node_idx  # True if both conditions are met
+    max_disturb[small_dist_on_land] = fillvalue
+
+    return max_disturb, time_idx_max_disturb
+
+
+def sample_usage():
+    """
+    Sample usage of cal_maxelev and cal_disturbance
+    """
+
+    # read in multiple SCHISM output files at a time,
+    # useful for operational forecast where each day contains two files
+    ds = xr.open_mfdataset([
+        './outputs/schout_adcirc_20240926.nc',
+        './outputs/schout_adcirc_20240927.nc',
+    ], concat_dim='time', combine='nested')
+
+    # read in the depth and elevation
+    elev = ds['zeta'].values
+    depth = ds['depth'].values
+    # accomodate for a bug which set depth's dimension to (time, x, y)
+    depth = depth[0, :] if depth.ndim == 2 else depth
+
+    # calculate the maximum elevation and the time of the maximum elevation
+    max_elev, time_idx_max_elev = cal_maxelev(depth, elev)
+    # calculate the maximum disturbance and the time of the maximum disturbance
+    max_disturb, time_idx_max_disturb = cal_disturbance(depth, elev, city_node_idx_file=None)
+
+    # Extra detail: mask small disturbance in cities,
+    # which need city node indices from a file (different between operation and shadow forecasts)
+    # max_disturb, time_idx_max_disturb = cal_disturbance(
+    #     depth, elev, city_node_idx_file='./inputs/city_poly.node_id.oper.txt')
+
+    # print some diagnostics
+    max_elev = np.ma.masked_values(max_elev, -99999.0)
+    max_disturb = np.ma.masked_values(max_disturb, -99999.0)
+    print(f'time_idx_max_elev: min={np.min(time_idx_max_elev)}, max={np.max(time_idx_max_elev)}')
+    print(f'max_disturb: min={np.nanmin(max_disturb)}, max={np.nanmax(max_disturb)}')
+
+
+if __name__ == '__main__':
+    sample_usage()

src/Utility/Pre-Processing/STOFS-3D-Atl-operation/pysh/README

@@ -1,3 +1,4 @@
-Changed from v2.1:
-extract_slab_fcst_netcdf4.py
-generate_adcirc.py
+extract_slab_fcst_netcdf4_current.py is used in the latest operational forecast
+extract_slab_fcst_netcdf4.py is under testing
+
+generate_adcirc.py has a bug fix, which has not been incorporated into operational forecast yet

src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/Bathy_edit/bathy_edit.py

@@ -120,8 +120,14 @@ def prepare_dir(wdir: Path, tasks: str):
     the larger files not in the Git repository.
     '''
     script_dir = Path(__file__).parent
-    for task in tasks:
-        shutil.copytree(f'{script_dir}/{task}/', f'{wdir}/{task}/', symlinks=False, dirs_exist_ok=True)
+    if script_dir == wdir:
+        print('The script is already in the working directory; no need to copy.')
+    else:
+        print(f'Copying the script and the subdirectories to {wdir}')
+        for task in tasks:
+            shutil.copytree(
+                f'{script_dir}/{task}/', f'{wdir}/{task}/',
+                symlinks=False, dirs_exist_ok=True)
 
     # copy larger files not in the Git repository
     for task, file_path_list in LARGE_FILES.items():
@@ -150,7 +156,8 @@ def bathy_edit(wdir: Path, hgrid_fname: Path, tasks: list = None):
         from Regional_tweaks.regional_tweaks import tweak_hgrid_depth
         hgrid_obj = tweak_hgrid_depth(
             hgrid=hgrid_obj, regions_dir=f'{wdir}/Regional_tweaks/regions/')
-        print("Finished setting regional tweaks.\n")
+        initial_dp = hgrid_obj.dp.copy()  # treat the regional tweaks as the initial dp
+        print("Finished setting regional tweaks and updating initial dp.\n")
 
     if 'NCF' in tasks:  # load NCF (National Channel Framework)
         from NCF.load_NCF import load_NCF
@@ -232,10 +239,10 @@ def sample_usage():
     '''
     Sample usage of the bathy_edit function.
     '''
-    WDIR = Path('/work/noaa/nosofs/feiye/STOFS-3D-Atl-Example-Setup/Bathy_edit_example/')
+    WDIR = Path('/sciclone/schism10/feiye/STOFS3D-v8/I10/Bathy_edit/')
     HGRID_FNAME = Path(  # Typically, this is the DEM-loaded hgrid
-        '/work/noaa/nosofs/feiye/STOFS-3D-Atl-Example-Setup/'
-        'DEM_loading_example/hgrid.ll.dem_loaded.mpi.gr3'
+        '/sciclone/schism10/feiye/STOFS3D-v8/I10/Bathy_edit/'
+        'DEM_loading/hgrid.ll.dem_loaded.mpi.gr3'
     )
     TASKS = ['Regional_tweaks', 'NCF', 'Levee']
 

src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/Grid/clip_autoarcs.py

@@ -13,7 +13,7 @@
 import geopandas as gpd
 
 # ------------------------- inputs ---------------------------
-WDIR = '/sciclone/schism10/Hgrid_projects/STOFS3D-v8/v22/Clip/'
+WDIR = '/sciclone/schism10/Hgrid_projects/STOFS3D-v8/v23/Clip/'
 CRS = 'esri:102008'
 
 # manual polygons defined in the coastal coverage,

src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/stofs3d_atl_driver.py

@@ -491,40 +491,41 @@ def main():
 
     # -----------------input---------------------
     # hgrid generated by SMS, pre-processed, and converted to *.gr3
-    hgrid_path = ('/work/noaa/nosofs/feiye/STOFS-3D-Atl-Example-Setup/'
-                  '2D_Setup/hgrid_levee_xGEOID20b.gr3')
+    hgrid_path = ('/sciclone/schism10/feiye/STOFS3D-v8/I10a/hgrid.gr3')
 
     # get a configuration preset and make changes if necessary
     # alternatively, you can set the parameters directly on an
     # new instance of ConfigStofs3dAtlantic
-    config = ConfigStofs3dAtlantic.v7_hercules_test()
-    config.rnday = 3
+    config = ConfigStofs3dAtlantic.v7()
+    config.rnday = 36
     config.startdate = datetime(2024, 3, 5)
+    config.nwm_cache_folder = Path(
+        '/sciclone/schism10/feiye/STOFS3D-v8/I09/Source_sink/original_source_sink/20240305/')
 
     # define the project dir, where the run folders are located
-    project_dir = '/work/noaa/nosofs/feiye/STOFS-3D-Atl-Example-Setup/2D_Setup/'
+    project_dir = '/sciclone/schism10/feiye/STOFS3D-v8/'
 
     # run ID. e.g, 00, 01a, 10b, etc.; the run folders are named using this ID as follows:
     # I{runid}: input directory for holding model input files and scripts;
     # R{runid}: run directory, where the run will be submitted to queue;
     # O{runid}: output directory for holding raw outputs and post-processing.
     # under project_dir
-    runid = '00'
+    runid = '10a'
 
-    # swithes to generate differ
+    # swithes to generate different input files
     input_files = {
-        'bctides': True,
+        'bctides': False,
         'vgrid': False,
-        'gr3': True,
+        'gr3': False,
         'nudge_gr3': False,
         'shapiro': True,
         'drag': True,
         'elev_ic': True,
-        'source_sink': True,
+        'source_sink': False,
         'hotstart.nc': False,
         '*D.th.nc': False,
         '*nu.nc': False,
-        '*.prop': True,
+        '*.prop': False,
     }
     # -----------------end input---------------------
 
@@ -780,7 +781,7 @@ def main():
         gen_tvd_prop(hgrid, config.tvd_regions)
 
         os.chdir(run_dir)
-        os.system(f'ln -sf ../I{runid}/{sub_dir}/tvd_prop .')
+        os.system(f'ln -sf ../I{runid}/{sub_dir}/tvd.prop .')
 
     # -----------------hotstart.nc---------------------
     if input_files['hotstart.nc']: