Project scripts: updates on a few scripts on grid manipulation.

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

@@ -65,17 +65,21 @@ def load_chart(hg:schism_grid, sounding_shpfile:Path, region_shpfile:Path, crs_r
     hg_points = gpd.GeoDataFrame(geometry=gpd.points_from_xy(hg.x, hg.y), crs='epsg:4326').to_crs(crs_region)
     joined_gdf = gpd.sjoin(hg_points, channel_polys, how="inner", predicate='within')
 
-    idx = joined_gdf.index.to_numpy()  # get the indices of the points inside the polygons
-    in_channel = np.zeros_like(hg.dp, dtype=bool)
-    in_channel[idx] = True
+    inchannel_idx = joined_gdf.index.to_numpy()  # get the indices of the points inside the polygons
 
     # diagnostic plot
-    plt.figure()
-    hg.plot(value=in_channel.astype(int), fmt=1)
+    # in_channel = np.zeros_like(hg.dp, dtype=bool)
+    # in_channel[inchannel_idx] = True
+    # plt.figure()
+    # hg.plot(value=in_channel.astype(int), fmt=1)
+    # plt.show()
+
+    # find the nearest sounding_xyz point using k-d tree
+    idx = KDTree(sounding_xyz[:, :2]).query(np.c_[hg.x[inchannel_idx], hg.y[inchannel_idx]])[1]
+    inchannel_sounding_z = sounding_xyz[idx, 2]
 
-    # for inpolygon nodes, find the nearest sounding_xyz point using k-d tree
-    idx = KDTree(sounding_xyz[:, :2]).query(np.c_[hg.x, hg.y])[1]
-    dp_sounding = sounding_xyz[idx, 2]
+    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
 
     return dp_sounding
 
@@ -110,8 +114,9 @@ def load_NCF(hg:schism_grid, NCF_shpfile:Path):
     dp_NCF[idx] = joined_gdf['depthmaint'].to_numpy() * 0.3048  # convert from feet to meters
 
     # diagnostic plot
-    plt.figure()
-    hg.plot(value=dp_NCF.astype(int), fmt=1)
+    # plt.figure()
+    # hg.plot(value=dp_NCF.astype(int), fmt=1)
+    # plt.show()
 
     return dp_NCF
 
@@ -133,19 +138,29 @@ def plot_diagnostic(hg:schism_grid, dp:np.ndarray):
 if __name__ == '__main__':
     hg_file = Path('/sciclone/schism10/feiye/Test/RUN02b_JZ/hgrid.gr3')
     hg = read_schism_hgrid_cached(hg_file)
+    dp_orig = hg.dp.copy()
+
+    # 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)
 
     # 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_2_xyz_edited.shp'),
-        region_shpfile=Path('/sciclone/schism10/Hgrid_projects/Charts/Savanna_Cooper/manual_chart_region.shp'), crs_region='esri:102008'
+        sounding_shpfile=Path('/sciclone/schism10/Hgrid_projects/Charts/Savanna_Cooper/savannah_cooper_sounding_3_xyz_edited.shp'),
+        region_shpfile=Path('/sciclone/schism10/Hgrid_projects/Charts/Savanna_Cooper/secofs_chart_loading_zones.shp'), crs_region='esri:102008'
     )
-    hg.dp = np.maximum(hg.dp, dp_from_chart)  # change the depth only if it is deeper than the original depth
+    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"))
 
     # 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(hg.dp, dp_from_NCF)  # change the depth only if it is deeper than the original depth
+    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
+    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"))
 
     pass

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

@@ -0,0 +1,50 @@
+'''
+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
+from shapely.ops import polygonize
+
+wdir = '/sciclone/schism10/feiye/Test/RUN02b_JZ/Dredge/'
+crs = 'esri:102008'
+'''
+coastal = gpd.read_file(f'{wdir}/coastal.shp')
+coastal.set_crs(crs, inplace=True)
+
+lbnd_coastal = gpd.read_file(f'{wdir}/lbnd_coastal.shp')
+lbnd_coastal.set_crs(crs, inplace=True)
+
+raw_watershed = lbnd_coastal.overlay(coastal, how='difference').dissolve()
+
+# make the refined coastal polygons, i.e., coastal minus those intersecting the select_nearshore polygons
+select_nearshore = gpd.read_file(f'/sciclone/data10/feiye/SCHISM_REPOSITORY/schism/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/Grid/select_nearshore_v2.shp')
+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)
+
+levee_buf = gpd.GeoDataFrame(geometry=gpd.read_file(f'{wdir}/levee.shp').buffer(50))
+
+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')
+
+# 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)
+
+total_arcs_polygons = gpd.read_file(f'{wdir}/local_river_polys.shp').to_crs(crs)
+total_arcs_polygons_clipped = total_arcs_polygons.clip(watershed)
+total_arcs_polygons_clipped.to_file(f'{wdir}/total_river_polys_clipped.shp', crs=crs)
+
+pass

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

@@ -9,26 +9,41 @@
 '''
 
 import geopandas as gpd
+from shapely.ops import polygonize
 
-wdir = '/sciclone/schism10/Hgrid_projects/STOFS3D-v7/new19/auto_arc_process/'
+wdir = '/sciclone/schism10/feiye/Test/RUN02b_JZ/Dredge/'
 crs = 'esri:102008'
 
-coastal = gpd.read_file(f'{wdir}/coastal.shp', crs=crs)
-lbnd_coastal = gpd.read_file(f'{wdir}/lbnd_coast.shp', crs=crs)
+coastal = gpd.read_file(f'{wdir}/coastal.shp')
+coastal.set_crs(crs, inplace=True)
+
+lbnd_coastal = gpd.read_file(f'{wdir}/lbnd_coastal.shp')
+lbnd_coastal.set_crs(crs, inplace=True)
 
 raw_watershed = lbnd_coastal.overlay(coastal, how='difference').dissolve()
 
-coastal_refined_buf = gpd.GeoDataFrame(geometry=gpd.read_file(f'{wdir}/coastal_refined.shp', crs=crs).buffer(50))
-levee_buf = gpd.GeoDataFrame(geometry=gpd.read_file(f'{wdir}/levee.shp', crs=crs).buffer(50))
+# make the refined coastal polygons, i.e., coastal minus those intersecting the select_nearshore polygons
+select_nearshore = gpd.read_file(f'/sciclone/data10/feiye/SCHISM_REPOSITORY/schism/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/Grid/select_nearshore_v2.shp')
+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)
+
+levee_buf = gpd.GeoDataFrame(geometry=gpd.read_file(f'{wdir}/levee.shp').buffer(50))
 
-watershed = raw_watershed.overlay(coastal_refined_buf, how='difference').dissolve()
+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', crs=crs)
+watershed = gpd.read_file(f'{wdir}/watershed.shp')
 
 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)
 
+# 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)
+# total_arcs_polygons_clipped.to_file(f'{wdir}/total_river_polys_clipped.shp', crs=crs)
+
 pass

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

@@ -0,0 +1,31 @@
+from pylib_essentials.schism_file import read_schism_hgrid_cached, grd2sms
+from pathlib import Path
+import geopandas as gpd
+import numpy as np
+
+wdir = Path('/sciclone/schism10/feiye/Test/RUN02b_JZ/Dredge/')
+dredge_depth = 5
+
+hg = read_schism_hgrid_cached(f'{wdir}/hgrid_chart_NCF_loaded.gr3')
+
+# load river polygons
+river_polys = gpd.read_file(wdir / 'total_river_polys_clipped_dissolved.shp')
+# 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')
+
+# determine in-channel nodes
+hg_points = gpd.GeoDataFrame(geometry=gpd.points_from_xy(hg.x, hg.y), crs='epsg:4326')
+joined_gdf = gpd.sjoin(hg_points, river_polys, how="inner", op='within')
+idx = joined_gdf.index.to_numpy() 
+
+in_channel = np.zeros_like(hg.dp, dtype=bool)
+in_channel[idx] = True
+hg.plot(value=in_channel.astype(int), fmt=1)
+
+# 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)
+
+pass

src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/Grid/select_nearshore_v2.map

@@ -0,0 +1,91 @@
+MAP VERSION 8
+BEGCOV
+COVFLDR "select_nearshore_v2"
+COVNAME "select_nearshore_v2"
+COVELEV 328.000000
+COVID 26020
+COVGUID 3f94fda7-550c-4ce5-ba6c-a6c96e98770d
+COVATTS VISIBLE 0
+COVATTS ACTIVECOVERAGE select_nearshore_v2
+COV_WKT PROJCS["North_America_Albers_Equal_Area_Conic",GEOGCS["GCS_North_American_1983",DATUM["North_American_Datum_1983",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["False_Easting",0],PARAMETER["False_Northing",0],PARAMETER["longitude_of_center",-96],PARAMETER["Standard_Parallel_1",20],PARAMETER["Standard_Parallel_2",60],PARAMETER["latitude_of_center",40],UNIT["Meter",1],AUTHORITY["EPSG","102008"]]END_COV_WKT 
+COV_VERT_DATUM 0
+COV_VERT_UNITS 0
+COVATTS PROPERTIES MESH_GENERATION
+NODE
+XY 1165395.1907756 -1184142.815209 328.0
+ID 2
+END
+ARC
+ID 1
+ARCELEVATION 316.285714
+NODES        2        2
+ARCVERTICES 26
+1223700.000000000000000 -1296000.000000000000000 328.000000000000000
+1252700.000000000000000 -1384200.000000000000000 328.000000000000000
+1342100.000000000000000 -1520700.000000000000000 328.000000000000000
+1377240.000000000000000 -1667400.000000000000000 0.000000000000000
+1479499.649472640594468 -1653706.501646494958550 328.000000000000000
+1519310.000000000000000 -1616980.000000000000000 328.000000000000000
+1544340.192105188500136 -1563126.583915059454739 328.000000000000000
+1545018.906015781220049 -1554071.209364138310775 328.000000000000000
+1559500.000000000000000 -1502100.000000000000000 328.000000000000000
+1541600.000000000000000 -1382700.000000000000000 328.000000000000000
+1479800.000000000000000 -1204200.000000000000000 328.000000000000000
+1429600.000000000000000 -1108600.000000000000000 328.000000000000000
+1441700.000000000000000 -830000.000000000000000 328.000000000000000
+1688900.000000000000000 -537100.000000000000000 328.000000000000000
+1785500.000000000000000 -355300.000000000000000 328.000000000000000
+1763600.000000000000000 -183900.000000000000000 328.000000000000000
+1773000.000000000000000 5100.000000000000000 328.000000000000000
+1984900.000000000000000 -65500.000000000000000 328.000000000000000
+1794800.000000000000000 -461300.000000000000000 328.000000000000000
+1548600.000000000000000 -879800.000000000000000 328.000000000000000
+1588100.000000000000000 -1260000.000000000000000 328.000000000000000
+1600600.000000000000000 -1612100.000000000000000 328.000000000000000
+1520600.000000000000000 -1721100.000000000000000 328.000000000000000
+1262000.000000000000000 -1760600.000000000000000 328.000000000000000
+834100.000000000000000 -1315000.000000000000000 328.000000000000000
+838200.000000000000000 -1234000.000000000000000 328.000000000000000
+DISTNODE 0
+NODETYPE 0
+ARCBIAS 1
+MERGED 0 0 0 0 
+END
+POLYGON
+ARCS 0
+HARCS 1
+1
+ID 0
+POLYMIN 1.90838e-08 1.90838e-08
+POLYMAX 3.40282e+38 3.40282e+38
+POLYINTERPTYPE 0 0
+POLYEXTRAPTYPE 1 1
+POLYEXTRAP 0 0
+POLYTRUN 0 0
+MAT
+MN	1	material 01
+MC	1	255	0		0
+MS	1	0
+PAVE 0.3
+CONSTANTBATH 0
+END
+POLYGON
+ARCS 1
+1  0 0
+ID 1
+POLYMIN 1.90838e-08 1.90838e-08
+POLYMAX 3.40282e+38 3.40282e+38
+POLYINTERPTYPE 0 0
+POLYEXTRAPTYPE 1 1
+POLYEXTRAP 0 0
+POLYTRUN 0 0
+MAT
+MN	1	material 01
+MC	1	255	0		0
+MS	1	0
+PAVE 0.3
+CONSTANTBATH 0
+END
+ENDCOV
+BEGTS
+LEND

src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/Grid/select_nearshore_v2.prj

@@ -0,0 +1 @@
+PROJCS["North_America_Albers_Equal_Area_Conic",GEOGCS["GCS_North_American_1983",DATUM["North_American_Datum_1983",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["False_Easting",0],PARAMETER["False_Northing",0],PARAMETER["longitude_of_center",-96],PARAMETER["Standard_Parallel_1",20],PARAMETER["Standard_Parallel_2",60],PARAMETER["latitude_of_center",40],UNIT["Meter",1],AUTHORITY["EPSG","102008"]]