STOFS-3D scripts: use pyschims' retries for *3D.th.nc and nudging in …

…the stofs_driver script; other minor fixes/edits.

src/Utility/Pre-Processing/STOFS-3D-Atl-operation/Source_sink/v7/gen_sourcesink.py

@@ -41,8 +41,12 @@ def main():
     # startdate = datetime(2024, 3, 5)
 
     # ------------------------- hardwired inputs for operation--------------------------
-    working_dir = '/sciclone/schism10/feiye/STOFS3D-v7/I12w/Source_sink/relocated_source_sink2/'
-    nwm_folder = '/sciclone/schism10/feiye/STOFS3D-v7/I12w/Source_sink/original_source_sink/20240305/'
+    #zy-working_dir = '/sciclone/schism10/feiye/STOFS3D-v7/I12w/Source_sink/relocated_source_sink2/'
+    #zy-nwm_folder = '/sciclone/schism10/feiye/STOFS3D-v7/I12w/Source_sink/original_source_sink/20240305/'
+    working_dir = './'
+    nwm_folder = './'
+
+
     layer = 'conus'
     # ------------------------- end hardwired inputs--------------------------
 

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

@@ -50,7 +50,7 @@
     'Feeder',  # set feeder channel depth, relative, datum doesn't matter
 ]
 
-DEFAULT_TASKS = {'Regional_tweaks', 'NCF', 'Levee'}
+DEFAULT_TASKS = {'Regional_tweaks', 'NCF', 'Levee', 'xGEOID'}
 
 # larger files not included in the Git repository, need to be copied to the working directory
 LARGE_FILES = {
@@ -245,10 +245,9 @@ def sample_usage():
     '''
     Sample usage of the bathy_edit function.
     '''
-    WDIR = Path('/sciclone/schism10/feiye/STOFS3D-v8/I10/Bathy_edit2/')
+    WDIR = Path('/sciclone/schism10/feiye/STOFS3D-v8/I14/Bathy_edit/')
     HGRID_FNAME = Path(  # Typically, this is the DEM-loaded hgrid
-        '/sciclone/schism10/feiye/STOFS3D-v8/I10/Bathy_edit/'
-        'DEM_loading/hgrid.ll.dem_loaded.mpi.gr3'
+        '/sciclone/schism10/feiye/STOFS3D-v8/I14/09b.gr3'
     )
     TASKS = ['xGEOID']
 

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

@@ -32,9 +32,9 @@ def prep_folder(wdir):
     os.chdir(wdir)
 
     # copy over the polygons
-    source_path = ('/sciclone/schism10/hjyoo/task/task6_SECOFS/'
-                   'simulation/Whole_Domain/Grid/Script/xGEOID/VDatum_polygons/')
-    os.system(f'cp -rL {source_path} .')
+    # source_path = ('/sciclone/schism10/hjyoo/task/task6_SECOFS/'
+    #                'simulation/Whole_Domain/Grid/Script/xGEOID/VDatum_polygons/')
+    # os.system(f'cp -rL {source_path} .')
 
 
 def generate_input_txt(hgrid_obj, wdir, n_sub=500000):

src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/Source_sink/NWM/gen_sourcesink_nwm.py

@@ -36,7 +36,7 @@ def gen_sourcesink_nwm(startdate: datetime, rnday: float, cache_folder: Path = N
     output_directory = Path.cwd()
 
     # input directory which saves nc files
-    if cache_folder is not None and cache_folder.exists():  # if cache folder exists, use it
+    if cache_folder is not None and os.path.exists(cache_folder):  # if cache folder exists, use it
         cache = cache_folder
     else:  # if cache folder does not exist, create it
         cache = Path(f'./{startdate.strftime("%Y%m%d")}')

src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/Source_sink/Replace_with_USGS/replace_with_obs.py

@@ -450,10 +450,15 @@ def source_nwm2usgs(
     )
     # build df
     nwm_df = pd.DataFrame({'datetime': my_date_range}).join(pd.DataFrame(stream_flow, columns=fid_subset))
+    # assign utc time zone
+    nwm_df['datetime'] = nwm_df['datetime'].dt.tz_localize('UTC')
     nwm_df = nwm_df.loc[:, ~nwm_df.columns.duplicated()]  # drop duplicate columns
     # nwm_df = nwm_df.set_index('datetime').resample("H").mean()
     # nwm_df = nwm_df.interpolate(method='time')  # fill in missing data
 
+    # set datetime as index
+    nwm_df = nwm_df.set_index('datetime')
+
     print(f'---------------reading nwm data took: {time.time()-t1} s ---------------\n')
 
     # -----------------------------------------------------------------------------------
@@ -505,9 +510,7 @@ def source_nwm2usgs(
             print(f'warning: source {i}_{my_source.hgrid_ie} does not have nearby usgs obs')
             continue
 
-        this_vsource_nwm_array = nwm_df[mysrc_nwm_fid].to_numpy()
-
-        # interpolate usgs to vsource/nwm time
+        # interpolate usgs to vsource time
         this_vs_usgs_array = np.zeros((nrec, len(mysrc_usgs_id)), dtype=float)
         for k, idx in enumerate(mysrc_usgs_idx):
             this_vs_usgs_array[:, k] = np.interp(
@@ -516,10 +519,21 @@ def source_nwm2usgs(
                 my_obs.stations[idx].df.iloc[:, 1].to_numpy()*0.028316847
             )
 
+        # this_vsource_nwm_array = nwm_df[mysrc_nwm_fid].to_numpy()
+        this_vsource_nwm_array = np.zeros((nrec, len(mysrc_nwm_fid)), dtype=float)
+        # interpolate nwm to vsource time
+        for k, idx in enumerate(mysrc_nwm_fid):
+            this_vsource_nwm_array[:, k] = np.interp(
+                (df.index - df.index[0]).total_seconds().to_numpy(),
+                (nwm_df.index - df.index[0]).total_seconds().to_numpy(),
+                nwm_df[mysrc_nwm_fid].to_numpy().squeeze()
+            )
+
+        plt.figure(figsize=(10, 6)) 
         plt.plot(df.index, df['Data'], '.k', label="original vsource")
         for k, usgs_id in enumerate(mysrc_usgs_id):
-            plt.plot(df.index, this_vs_usgs_array[:, k], linestyle1[k], label=f'usgs {usgs_id}')
-            plt.plot(df.index, this_vsource_nwm_array[:, k], linestyle2[k], label=f'NWM near usgs {usgs_id}')
+            plt.plot(df.index, this_vs_usgs_array[:, k], linestyle1[k], linewidth=0.5, label=f'usgs {usgs_id}')
+            plt.plot(df.index, this_vsource_nwm_array[:, k], linestyle2[k], linewidth=0.5, label=f'NWM near usgs {usgs_id}')
 
         # take the diff (including all relavant usgs-nwm pairs found)
         diff_usgs_nwm = this_vs_usgs_array - this_vsource_nwm_array
@@ -543,10 +557,18 @@ def source_nwm2usgs(
 
 if __name__ == "__main__":
     # sample usage (the final product is "adjusted_vsource.th" in the output_dir):
+    # source_nwm2usgs(
+    #     start_time_str="2024-03-05 00:00:00",
+    #     f_shapefile="/sciclone/schism10/Hgrid_projects/STOFS3D-v8/v20p2s2_RiverMapper/shapefiles/LA_nwm_v1p2.shp",
+    #     original_ss_dir='/sciclone/schism10/feiye/STOFS3D-v8/I09/Source_sink/original_source_sink/',
+    #     nwm_data_dir='/sciclone/schism10/feiye/STOFS3D-v8/I09/Source_sink/original_source_sink/20240305/',
+    #     output_dir='/sciclone/schism10/feiye/STOFS3D-v8/I09/Source_sink/USGS_adjusted_sources/',
+    # )
+
     source_nwm2usgs(
-        start_time_str="2024-03-05 00:00:00",
+        start_time_str="2017-12-01 00:00:00",
         f_shapefile="/sciclone/schism10/Hgrid_projects/STOFS3D-v8/v20p2s2_RiverMapper/shapefiles/LA_nwm_v1p2.shp",
-        original_ss_dir='/sciclone/schism10/feiye/STOFS3D-v8/I09/Source_sink/original_source_sink/',
-        nwm_data_dir='/sciclone/schism10/feiye/STOFS3D-v8/I09/Source_sink/original_source_sink/20240305/',
-        output_dir='/sciclone/schism10/feiye/STOFS3D-v8/I09/Source_sink/USGS_adjusted_sources/',
+        original_ss_dir='/sciclone/schism10/feiye/STOFS3D-v8/I14/Source_sink/USGS_adjusted_sources/original/',
+        nwm_data_dir='/sciclone/schism10/feiye/STOFS3D-v8/I13/Source_sink/original_source_sink/20171201/',
+        output_dir='/sciclone/schism10/feiye/STOFS3D-v8/I14/Source_sink/USGS_adjusted_sources/',
     )

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

@@ -21,6 +21,7 @@ def __init__(
         self,
         startdate=datetime(2017, 12, 1),  # start date of the model
         rnday=60,  # number of days to run the model
+        ocean_bnd_ids=[0],  # list of open boundary ids for *D.th.nc
         elev2d_uniform_shift=0.0,  # add a uniform shift to elev2D
         nudging_zone_width=1.5,  # in degrees
         nudging_day=1.0,  # in days
@@ -40,6 +41,7 @@ def __init__(
 
         self.startdate = startdate
         self.rnday = rnday
+        self.ocean_bnd_ids = ocean_bnd_ids
         self.elev2d_uniform_shift = elev2d_uniform_shift
         self.nudging_zone_width = nudging_zone_width
         self.nudging_day = nudging_day
@@ -110,6 +112,7 @@ def v6(cls):
                 '/sciclone/schism10/feiye/STOFS3D-v5/Inputs/v14/Parallel/'
                 'SMS_proj/feeder/feeder.pkl'),
             mandatory_sources_coor=rsf.v19p2_mandatory_sources_coor,
+            ocean_bnd_ids=[0, 1],
             bc_flags=[[5, 5, 4, 4]],
             tvd_regions=[
                 'tvd0_1.reg', 'tvd0_2.reg', 'tvd0_3.reg', 'tvd0_4.reg',
@@ -122,6 +125,7 @@ def v6(cls):
     def v7(cls):
         '''Factory method to create a configuration for STOFS3D-v7'''
         return cls(
+            ocean_bnd_ids=[0, 1],
             elev2d_uniform_shift=-0.42,  # add a uniform shift to elev2D
             nudging_zone_width=7.3,  # default nudging zone
             shapiro_zone_width=11.5,  # default shapiro zone
@@ -157,6 +161,7 @@ def v7(cls):
     def v7_hercules_test(cls):
         '''Factory method to create a configuration for STOFS3D-v7 2D setup'''
         return cls(
+            ocean_bnd_ids=[0, 1],
             elev2d_uniform_shift=-0.42,  # add a uniform shift to elev2D
             nudging_zone_width=7.3,  # default nudging zone
             shapiro_zone_width=11.5,  # default shapiro zone
@@ -183,19 +188,23 @@ def v7_hercules_test(cls):
     def v8_louisianna(cls):
         '''Factory method to create a configuration for STOFS3D-v8's local test in Louisianna'''
         return cls(
+            ocean_bnd_ids=[0],
             nudging_zone_width=0,  # default nudging zone
             shapiro_zone_width=0,  # default shapiro zone
             shapiro_tilt=0,  # default abrupt transition in the shapiro zone
             feeder_info_file='',
             relocate_source=False,
             nwm_cache_folder=Path('/sciclone/schism10/whuang07/schism20/NWM_v2.1/'),
-            bc_flags=[[5, 5, 4, 4]]
+            bc_flags=[[5, 3, 0, 0]],
+            bc_relax=[[None, None, None, None]],
+            bc_const=[[None, None, None, None]],
         )
 
     @classmethod
     def v8(cls):
         '''Factory method to create a configuration for STOFS3D-v8 3D setup'''
         return cls(
+            ocean_bnd_ids=[0, 1],
             elev2d_uniform_shift=-0.42,  # add a uniform shift to elev2D
             nudging_zone_width=7.3,  # default nudging zone
             shapiro_zone_width=11.5,  # default shapiro zone