STOFS3D scripts: remove leading zeros in *source.th for operation. Ot…

…her minor changes.

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

@@ -84,7 +84,7 @@ def main():
     idx = np.argsort(np.array([int(key) for key in source2fid_dict.keys()]))
     # write source time history file
     np.savetxt(f'{working_dir}/vsource.th',
-               np.c_[times, sources[:, idx]], fmt='%10.4f', delimiter=' ')
+               np.c_[times, sources[:, idx]], fmt='%.4f', delimiter=' ')
 
 
 if __name__ == '__main__':

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

@@ -18,6 +18,12 @@
 
 Some larger files are not included in the Git repository,
 you need to specify the paths in the "LARGE_FILES" dictionary below.
+
+Note that xGEOID is just a wrapper around vdatum.jar,
+which takes a lot of memory and quite slow.
+Running xGEOID on a large grid may take hours.
+Use viz.sciclone.wm.edu; other SciClone subclusters often run out of memory.
+This is deprecated, use Felicio's workflow after testing.
 '''
 
 import os
@@ -29,7 +35,7 @@
 
 from pylib import grd2sms, sms2grd
 try:  # c++ function to speed up the grid reading
-    from pylib_experimental.schism_file import cread_schism_hgrid as schism_read
+    from pylib_experimental.schism_file import xread_schism_hgrid as schism_read
 except ImportError:
     from pylib import schism_grid as schism_read
 
@@ -217,11 +223,11 @@ def bathy_edit(wdir: Path, hgrid_fname: Path, tasks: list = None):
         # A grid without feeder is needed to identify which feeder points are outside and should be deepened
         # Only the boundary matters, the interior of the grid doesn't matter,
         # so if you don't have a grid without feeders, you can just generate a simplified grid with the lbnd_ocean map
-        gd_no_feeder = sms2grd('/sciclone/schism10/Hgrid_projects/STOFS3D-v7/v20.0/no_feeder.2dm')
+        gd_no_feeder = sms2grd('/sciclone/schism10/Hgrid_projects/STOFS3D-v8/v23.1/hgrid.2dm')
         gd_no_feeder.proj(prj0='esri:102008', prj1='epsg:4326')
         initial_dp = hgrid_obj.dp.copy()
         hgrid_obj = set_feeder_dp(
-            feeder_info_dir='/sciclone/schism10/Hgrid_projects/STOFS3D-v7/v20.0/Feeder/',
+            feeder_info_dir='/sciclone/schism10/Hgrid_projects/STOFS3D-v7/v23.3/Feeder/',
             hgrid_obj=hgrid_obj, hgrid_obj_no_feeder=gd_no_feeder
         )
         dp_diff = initial_dp - hgrid_obj.dp
@@ -239,12 +245,12 @@ def sample_usage():
     '''
     Sample usage of the bathy_edit function.
     '''
-    WDIR = Path('/sciclone/schism10/feiye/STOFS3D-v8/I10/Bathy_edit/')
+    WDIR = Path('/sciclone/schism10/feiye/STOFS3D-v8/I10/Bathy_edit2/')
     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'
     )
-    TASKS = ['Regional_tweaks', 'NCF', 'Levee']
+    TASKS = ['xGEOID']
 
     bathy_edit(wdir=WDIR, hgrid_fname=HGRID_FNAME, tasks=TASKS)
 

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

@@ -223,7 +223,7 @@ def convert2xgeoid(wdir, hgrid_obj, diag_output=None):
 
     # this generates the input files for vdatum.jar, including
     # hgrid_stofs3d_inland_?.txt and hgrid_stofs3d_inland_ches_del.txt
-    generated_input_files = generate_input_txt(hgrid_obj=hgrid_obj, wdir=wdir, n_sub=100000)
+    generated_input_files = generate_input_txt(hgrid_obj=hgrid_obj, wdir=wdir, n_sub=500000)
 
     # see if the input files are complete
 
@@ -241,7 +241,7 @@ def convert2xgeoid(wdir, hgrid_obj, diag_output=None):
 
     # the first group should have no failed files, since they are strictly in region 4
     # this may not be true for other domains, so manually go over the workflow first before using the script
-    input_fnames = glob("*_[0-9].txt")
+    input_fnames = glob("hgrid*.txt")
     # Starting the processes
     processes = [subprocess.Popen(
             f"java -jar vdatum.jar  ihorz:NAD83_2011 ivert:navd88:m:{z_convention} "

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

@@ -103,6 +103,45 @@
     [-72.06195833333334, 41.525600000000004, np.nan, np.nan],  # Quinebaug River, CT
 ]).reshape(-1, 4)
 
+v23p3_mandatory_sources_coor = np.array([
+    [-73.90869, 42.13509, np.nan, np.nan],  # Hudson River, NY
+    [-74.94442, 40.34478, np.nan, np.nan],  # Delaware River, NJ
+    [-76.1011, 39.5802, np.nan, np.nan],  # Susquehanna River, VA
+    [-78.425288, 34.508177, np.nan, np.nan],  # Cape Fear River, NC
+    [-91.72306, 31.04462, np.nan, np.nan],  # Red River (upstream of Atchafalaya River), LA
+    [-80.10808, 33.50005, np.nan, np.nan],  # Santee River, SC
+    [-79.81703, 33.59694, np.nan, np.nan],  # Black River, SC
+    [-79.57210, 33.71223, np.nan, np.nan],  # Black Mingo Creek, SC
+    [-79.49997, 33.84686, np.nan, np.nan],  # Lynches River, SC
+    [-79.48467, 33.93939, np.nan, np.nan],  # Pee Dee River, SC
+    [-79.33247, 33.98196, np.nan, np.nan],  # Little Pee Dee River, SC
+    [-77.917829, 34.749979, np.nan, np.nan],  # Northeast Cape Fear River, NC
+    [-87.9523, 30.8472, np.nan, np.nan],  # Mobile River, AL
+    [-96.695401, 28.968284, -96.69652166667, 28.990345],  # Lavaca River, TX
+    [-96.548436, 28.999706, -96.554498, 29.024612666667],  # Lake Texana, TX
+    [-93.83342666667, 30.355123333333, -93.83342666667, 30.355123333333],  # Cypress Creek, TX
+    [-89.764476, 30.551926, -89.76781133333, 30.538070666667],  # Lotts Creek, LA
+    [-87.219805, 30.567296, -87.24471466667, 30.601442333333],  # Escambia River, FL
+    [-83.987035, 30.331327, np.nan, np.nan],  # Horsehead Creek and Little River, FL
+    [-83.928038, 30.30404, np.nan, np.nan],  # Bailey Mill Creek, FL
+    [-82.950913, 29.958097, -82.99605566667, 30.007415],  # Suwannee River, FL
+    [-81.02370433333333, 27.315079666666666, np.nan, np.nan],  # Kissimmee River, FL
+    [-81.997572, 30.786870, -82.040457, 30.74494233333333],  # St Marys River, FL
+    [-91.56184, 31.05043, np.nan, np.nan],  # Mississippi River
+    [-79.43425, 33.84487, -79.50974266666667, 33.85385866666667],  # Lyches River, SC
+    [-74.74868, 39.47915, -74.75470666666668, 39.485390333333335],  # Great Egg Harbor River, NJ
+    [-73.94009733333333, 42.06972966666667, np.nan, np.nan],  # Saugeties Creek, NY
+    [-73.971293, 41.920595999999996, np.nan, np.nan],  # Hudson River branch, NY
+    [-73.92918633333333, 41.592421333333334, np.nan, np.nan],  # Hudson River branch, NY
+    [-73.07229533333333, 41.303546000000004, np.nan, np.nan],  # Housatonic River, CT
+    [-72.625735, 41.656137666666666, np.nan, np.nan],  # Connecticut River, CT
+    [-72.64970633333333, 41.572111666666665, np.nan, np.nan],  # Mattabesset River, CT
+    [-72.470818, 41.47020933333334, np.nan, np.nan],  # Salmon River, CT
+    [-72.11158266666666, 41.455657333333335, np.nan, np.nan],  # Stony Brook, CT
+    [-72.090553, 41.535118000000004, np.nan, np.nan],  # Yantic River, CT
+    [-72.06195833333334, 41.525600000000004, np.nan, np.nan],  # Quinebaug River, CT
+]).reshape(-1, 4)
+
 
 def nearest_neighbour(points_a, points_b):
     '''A wrapper for scipy.spatial.cKDTree to find the nearest neighbour of points_a in points_b'''
@@ -130,6 +169,17 @@ def lonlat2cpp(lon, lat, lon0=0, lat0=0):
     return [xout, yout]
 
 
+def cpp2lonlat(x, y, lon0=0, lat0=0):
+    """
+    Convert cpp in meters to lon/lat
+    """
+    R_EARTH = 6378206.4
+    lon0_radian, lat0_radian = lon0/180*np.pi, lat0/180*np.pi
+    lon_radian = x / R_EARTH / np.cos(lat0_radian) + lon0_radian
+    lat0_radian = y / R_EARTH
+    return [lon_radian/np.pi*180, lat0_radian/np.pi*180]
+
+
 def relocate_sources2(
     old_ss_dir=None, outdir=None,
     feeder_info_file=None, hgrid_fname=None, allow_neglection=True,
@@ -262,6 +312,8 @@ def relocate_sources2(
             relocation_distance[i] = distance[target_old_source]
         else:
             if i < len(mandatory_sources_coor):
+                # mandatory sources are the first len(manadatory_sources_coor) new_sources,
+                # which must be relocated
                 raise ValueError(f'mandatory new source {i}: {x, y} cannot be mapped to an old source')
 
     # -------------------------------------clean up odd cases-------------------------------------
@@ -358,6 +410,7 @@ def relocate_sources2(
 
 def relocate_sources(
     old_ss_dir=None, outdir=None, relocate_map=None,
+    no_feeder=False,
     allow_neglection=True, feeder_info_file=None, hgrid_fname=None,
     max_search_radius=1000.0, mandatory_sources_coor=np.empty((0, 4)),
     region_list=None,
@@ -373,8 +426,10 @@ def relocate_sources(
         which is useful for the operational forecast to save time,
         in this case no need to provide the arguments below;
         otherwise, the relocation map will be generated.
+    : no_feeder: if true, feeder_bases (which should still be present in a mesh without feeders)
+        will be used; otherwise, feeder_heads will be used.
     :allow_neglection: Allow neglecting channels that cannot be matched to a new source location.
-        This is useful when the relocation is only for a specific region.
+        This is useful when the relocation is restricted to a specific region.
         If False, any unassigned old sources will be linked to the closest new source.
     :feeder_info_file:
         Feeder channel information generated by make_feeder_channel.py after RiverMapper
@@ -453,7 +508,13 @@ def relocate_sources(
         # -------------------------------------relocation-------------------------------------
         # Find matching source point at mandatory_sources_coor and feeders
         # These are the desired new source locations
-        new_sources_coor = np.r_[mandatory_sources_coor[:, :2], feeder_heads[:, :2]]
+        if no_feeder:
+            # If feeders are not implemented in the mesh, use "feeder_bases" instead of "feeder_heads",
+            # because "feeder_bases" are the locations where the feeders are to be placed.
+            # The relocation still works in the sense that only resolved rivers receive new sources.
+            new_sources_coor = np.r_[mandatory_sources_coor[:, :2], feeder_bases[:, :2]]
+        else:
+            new_sources_coor = np.r_[mandatory_sources_coor[:, :2], feeder_heads[:, :2]]
         # These are the locations used to search for the closest old source
         # ; in other words, to link the new source to the old source.
         # For the mandatory sources, these may be different from the new source locations
@@ -650,7 +711,7 @@ def test():
         old_ss_dir='/sciclone/schism10/feiye/STOFS3D-v7/Inputs/I12z/Source_sink/original_source_sink/',
         feeder_info_file='/sciclone/schism10/Hgrid_projects/STOFS3D-v7/v20.0/Feeder/feeder_heads_bases.xy',
         hgrid_fname=f'{wdir}/hgrid.gr3', outdir=wdir,
-        max_search_radius=2100, mandatory_sources_coor=v19p2_mandatory_sources_coor,
+        max_search_radius=2100, mandatory_sources_coor=v23p3_mandatory_sources_coor,
         allow_neglection=False
     )