Compound flooding scripts: Minor edits on the meshing script for SECOFS.

src/Utility/Pre-Processing/SECOFS/Bathy_edit/load_chart_NCF.py

@@ -81,6 +81,8 @@ def load_chart(hg:schism_grid, sounding_shpfile:Path, region_shpfile:Path, crs_r
     dp_sounding = np.ones_like(hg.dp, dtype=float) * -9999  # initialize with a large negative number
     dp_sounding[inchannel_idx] = inchannel_sounding_z  # only update the nodes in the channel
 
+    # hg.plot(value=dp_sounding, fmt=1)
+
     return dp_sounding
 
 
@@ -136,30 +138,47 @@ def plot_diagnostic(hg:schism_grid, dp:np.ndarray):
     pass
 
 if __name__ == '__main__':
-    hg_file = Path('/sciclone/schism10/feiye/Test/RUN02b_JZ/hgrid.gr3')
-    hg = read_schism_hgrid_cached(hg_file)
+    hg_file = Path('/sciclone/schism10/feiye/SECOFS/Inputs/I02e/Edit_bathy/Load_charts_NCF/hgrid_levee_loaded.gr3')
+    hg = schism_grid(str(hg_file))
     dp_orig = hg.dp.copy()
 
+    ''' prepare the chart data
     # extract xyz coordinates into an array
     sounding = gpd.read_file('/sciclone/schism10/Hgrid_projects/Charts/Savanna_Cooper/savannah_cooper_sounding_3_xyz_edited.shp')
     sounding_xyz = np.c_[np.array([point.coords[0] for point in sounding['geometry']]), sounding['z'].values]
     np.savetxt('/sciclone/schism10/Hgrid_projects/Charts/Savanna_Cooper/savannah_cooper_sounding_3_xyz_edited.txt', sounding_xyz)
 
+    # load txt file
+    xyz_xgeoid = np.loadtxt('/sciclone/schism10/Hgrid_projects/DEMs/vdatum/vdatum/result/reversed_input.txt')
+    xyz_xgeoid = - xyz_xgeoid  # convert from elevation to bathymetry
+    # convert from MLLW to xgeoid
+    idx = np.argwhere(abs(xyz_xgeoid[:, 2]) < 9e5).flatten()
+    mean_datum_shift = np.mean(xyz_xgeoid[idx, 2] - sounding_xyz[idx, 2])  # use average shift for invalid values
+    datum_shift = np.ones_like(xyz_xgeoid[:, 2]) * mean_datum_shift
+    datum_shift[idx] = xyz_xgeoid[idx, 2] - sounding_xyz[idx, 2]
+    sounding_xyz[:, 2] = sounding_xyz[:, 2] + datum_shift
+    # replace 'z' column in the shapefile
+    sounding['z'] = sounding_xyz[:, 2]
+    sounding.to_file('/sciclone/schism10/Hgrid_projects/Charts/Savanna_Cooper/savannah_cooper_sounding_3_xyz_edited_xgeoid.shp')
+    '''
+
+
     # the chart data needs a manual region file to limit the nodes to be loaded
     dp_from_chart = load_chart(
         hg=hg,
-        sounding_shpfile=Path('/sciclone/schism10/Hgrid_projects/Charts/Savanna_Cooper/savannah_cooper_sounding_3_xyz_edited.shp'),
+        sounding_shpfile=Path('/sciclone/schism10/Hgrid_projects/Charts/Savanna_Cooper/savannah_cooper_sounding_3_xyz_edited_xgeoid.shp'),
         region_shpfile=Path('/sciclone/schism10/Hgrid_projects/Charts/Savanna_Cooper/secofs_chart_loading_zones.shp'), crs_region='esri:102008'
     )
     hg.dp = np.maximum(dp_orig, dp_from_chart)  # change the depth only if it is deeper than the original depth
-    # grd2sms(hg, (f"{hg_file.parent}/{hg_file.stem}_chart_loaded.2dm"))
+    grd2sms(hg, (f"{hg_file.parent}/{hg_file.stem}_chart_loaded.2dm"))
 
     # the NCF data already defines the region
     dp_from_NCF = load_NCF(hg=hg, NCF_shpfile=Path('/sciclone/schism10/Hgrid_projects/Charts/Mississippi/channel_quarter_NCF.shp'))
     hg.dp = np.maximum(dp_orig, dp_from_NCF)  # change the depth only if it is deeper than the original depth
     # grd2sms(hg, (f"{hg_file.parent}/{hg_file.stem}_NCF_loaded.2dm"))
 
-    hg.dp = np.maximum(dp_orig, dp_from_NCF, dp_from_chart)  # change the depth only if it is deeper than the original depth
+    dp_new = np.maximum(dp_from_chart, dp_from_NCF)  # take the maximum of the two
+    hg.dp = np.maximum(dp_orig, dp_new)  # change the depth only if it is deeper than the original depth
     hg.save(f"{hg_file.parent}/{hg_file.stem}_chart_NCF_loaded.gr3", fmt=1)
     grd2sms(hg, (f"{hg_file.parent}/{hg_file.stem}_chart_NCF_loaded.2dm"))
 

src/Utility/Pre-Processing/SECOFS/Grid/clip_autoarcs.py

@@ -0,0 +1,57 @@
+'''
+This script clips the auto arcs to the watershed boundary,
+which is the part of domain that is not ocean, levee, or manually refined coastal polygons.
+A 50-m buffer is applied to the levee and refined coastal polygons to ensure that the auto arcs
+does not directly intersect with the levee and refined coastal polygons.
+However, no buffer is applied to the interface between watershed and ocean (i.e., shoreline)
+, allowing intersections between the auto arcs and shoreline to provide channels where
+there is no manual refinement.
+'''
+
+import geopandas as gpd
+
+# ------------------------- inputs ---------------------------
+wdir = '/sciclone/schism10/Hgrid_projects/SECOFS/new20_JZ/'
+crs = 'esri:102008'
+coastal_shpfile = f'{wdir}/coastal.shp'  # esri:102008
+lbnd_coastal_shpfile = f'{wdir}/lbnd_coastal.shp'  # esri:102008
+levee_shpfile = '/sciclone/schism10/feiye/SECOFS/Inputs/I02e/SMS/levee.shp'
+select_nearshore_shpfile = '/sciclone/data10/feiye/SCHISM_REPOSITORY/schism/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/Grid/select_nearshore_v2.shp'  
+# ------------------------- end inputs ---------------------------
+
+coastal = gpd.read_file(coastal_shpfile)
+coastal.set_crs(crs, inplace=True)
+
+lbnd_coastal = gpd.read_file(lbnd_coastal_shpfile)
+lbnd_coastal.set_crs(crs, inplace=True)
+# raw_watershed is a rough estimate of watershed region, which excludes manual polygons but may include levee
+raw_watershed = lbnd_coastal.overlay(coastal, how='difference').dissolve()
+
+# extract the refined coastal polygons, i.e., coastal minus those intersecting the select_nearshore polygons, where no intersections with auto arcs are allowed, hence the buffer
+select_nearshore = gpd.read_file(select_nearshore_shpfile)
+nearshore_coastal = gpd.sjoin(coastal, select_nearshore, how="inner", predicate='intersects')
+refined_coastal = coastal.overlay(nearshore_coastal, how='difference').dissolve()
+refined_coastal_buf = gpd.GeoDataFrame(geometry=refined_coastal.buffer(50))
+refined_coastal_buf.to_file(f'{wdir}/coastal_refined.shp', crs=crs)
+# auto-arcs cannot intersect with the levee, hence the buffer
+levee_buf = gpd.GeoDataFrame(geometry=gpd.read_file(levee_shpfile).buffer(50))
+# get final watershed by subtracting polygons that cannot be intersected by auto arcs from the raw watershed
+watershed = raw_watershed.overlay(refined_coastal_buf, how='difference').dissolve()
+watershed = watershed.overlay(levee_buf, how='difference').dissolve()
+watershed.to_file(f'{wdir}/watershed.shp', crs=crs)
+watershed = gpd.read_file(f'{wdir}/watershed.shp')
+
+# clip the auto arcs to the watershed boundary
+total_arcs = gpd.read_file(f'{wdir}/total_arcs.shp').to_crs(crs)
+total_arcs_clipped = total_arcs.clip(watershed)
+total_arcs_clipped.to_file(f'{wdir}/total_arcs_clipped.shp', crs=crs)
+
+# clip the polygons formed by auto arcs to the watershed boundary; this step may fail, in that case do it manually in qgis
+total_arcs_polygons = gpd.read_file(f'{wdir}/total_river_polys.shp').to_crs(crs)
+total_arcs_polygons_clipped = total_arcs_polygons.buffer(0).clip(watershed)
+# put total_arcs_polygons_clipped back to gdf and dissolve
+total_arcs_polygons_clipped = gpd.GeoDataFrame(geometry=total_arcs_polygons_clipped).dissolve()
+
+total_arcs_polygons_clipped.to_file(f'{wdir}/total_river_polys_clipped_test.shp', crs=crs)
+
+pass

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

@@ -1,15 +1,16 @@
-from pylib_essentials.schism_file import read_schism_hgrid_cached, grd2sms
+from pylib_essentials.schism_file import read_schism_hgrid_cached, grd2sms, schism_grid
 from pathlib import Path
 import geopandas as gpd
 import numpy as np
 
-wdir = Path('/sciclone/schism10/feiye/Test/RUN02b_JZ/Dredge/')
-dredge_depth = 5
+wdir = Path('/sciclone/schism10/feiye/SECOFS/Inputs/I02e/Edit_bathy/Dredge')
+hg_file = Path(f'{wdir}/hgrid_levee_loaded_chart_NCF_loaded.gr3')
+dredge_depth = 2
 
-hg = read_schism_hgrid_cached(f'{wdir}/hgrid_chart_NCF_loaded.gr3')
+hg = schism_grid(str(hg_file))  # epsg:4326
 
 # load river polygons
-river_polys = gpd.read_file(wdir / 'total_river_polys_clipped_dissolved.shp')
+river_polys = gpd.read_file(wdir / 'total_river_polys_clipped_dissolved.shp')  # epsg:4326
 # shrink the polygons by 1 m to exclude bank nodes
 river_polys.geometry = river_polys.geometry.to_crs('esri:102008').buffer(-1).to_crs('epsg:4326')
 
@@ -25,7 +26,7 @@
 # dredge the in-channel nodes
 hg.dp[idx] += dredge_depth
 
-grd2sms(hg, f'{wdir}/hgrid_chart_NCF_loaded_dredged_{dredge_depth}m.2dm')
-hg.save(f'{wdir}/hgrid_chart_NCF_loaded_dredged_{dredge_depth}m.gr3', fmt=1)
+grd2sms(hg, f'{wdir}/{hg_file.stem}_dredged_{dredge_depth}m.2dm')
+hg.save(f'{wdir}/{hg_file.stem}_dredged_{dredge_depth}m.gr3', fmt=1)
 
 pass