diff --git a/src/Utility/Pre-Processing/STOFS-3D-Atl-operation/Utils/post_process_lib.py b/src/Utility/Pre-Processing/STOFS-3D-Atl-operation/Utils/post_process_lib.py new file mode 100644 index 00000000..824fb638 --- /dev/null +++ b/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() diff --git a/src/Utility/Pre-Processing/STOFS-3D-Atl-operation/pysh/README b/src/Utility/Pre-Processing/STOFS-3D-Atl-operation/pysh/README index e41cbdda..935b2efe 100644 --- a/src/Utility/Pre-Processing/STOFS-3D-Atl-operation/pysh/README +++ b/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 diff --git a/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/Bathy_edit/bathy_edit.py b/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/Bathy_edit/bathy_edit.py index 9c3d8148..b4c690e0 100755 --- a/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/Bathy_edit/bathy_edit.py +++ b/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'] diff --git a/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/Grid/clip_autoarcs.py b/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/Grid/clip_autoarcs.py index 586225e4..1c465b56 100755 --- a/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/Grid/clip_autoarcs.py +++ b/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, diff --git a/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/stofs3d_atl_driver.py b/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/stofs3d_atl_driver.py index 75135f73..32f86ac2 100755 --- a/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/stofs3d_atl_driver.py +++ b/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']: