diff --git a/src/Utility/Pre-Processing/SECOFS/Bathy_edit/Chart_NCF/load_chart_NCF.py b/src/Utility/Pre-Processing/SECOFS/Bathy_edit/Chart/load_chart.py old mode 100644 new mode 100755 similarity index 71% rename from src/Utility/Pre-Processing/SECOFS/Bathy_edit/Chart_NCF/load_chart_NCF.py rename to src/Utility/Pre-Processing/SECOFS/Bathy_edit/Chart/load_chart.py index ef882d5a0..e4b460dcf --- a/src/Utility/Pre-Processing/SECOFS/Bathy_edit/Chart_NCF/load_chart_NCF.py +++ b/src/Utility/Pre-Processing/SECOFS/Bathy_edit/Chart/load_chart.py @@ -6,7 +6,6 @@ from matplotlib import pyplot as plt from pylib_essentials.schism_file import read_schism_hgrid_cached, grd2sms, schism_grid -from pylib_essentials.utility_functions import inside_polygon def prep_chart_data(sounding_shpfile:Path): ''' @@ -85,43 +84,6 @@ def load_chart(hg:schism_grid, sounding_shpfile:Path, region_shpfile:Path, crs_r return dp_sounding - -def load_NCF(hg:schism_grid, NCF_shpfile:Path): - ''' - Load the maintained depth from National Channel Framework data into the hgrid. - The original NCF polygons are enlarged to accommodate the mismatch between the hgrid and the NCF data. - The maintained depth is converted from feet to meters. - ''' - - # get the nodes in the NCF (National Channel Framework) polygons - NCF_shpfile = Path('/sciclone/schism10/Hgrid_projects/Charts/Mississippi/channel_quarter_NCF.shp') - NCF_data = gpd.read_file(NCF_shpfile) - - # remove the polygons that are not in the bounding box of the hgrid - NCF_data = NCF_data.cx[hg.x.min():hg.x.max(), hg.y.min():hg.y.max()] - # convert any multipolygons to polygons - NCF_data = NCF_data.explode() - - # enlarge the polygons by 4 m - NCF_data['geometry'] = NCF_data['geometry'].to_crs('esri:102008').buffer(4).to_crs('epsg:4326') - # extract the polygons to a list of nx2 numpy arrays - # NCF_polygons = [np.array(poly.exterior.coords) for poly in NCF_data['geometry']] - - # put hgrid points into a Point GeoDataFrame - hg_points = gpd.GeoDataFrame(geometry=gpd.points_from_xy(hg.x, hg.y), crs='epsg:4326') - # determine which points are inside the polygons and get the maintained depth - joined_gdf = gpd.sjoin(hg_points, NCF_data, how="inner", op='within') - idx = joined_gdf.index.to_numpy() # get the indices of the points inside the polygons - dp_NCF = np.ones_like(hg.dp, dtype=float) * -9999 # initialize with a large negative number - 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.show() - - return dp_NCF - def plot_diagnostic(hg:schism_grid, dp:np.ndarray): # # plot to check the changed nodes # plt.figure() @@ -138,7 +100,10 @@ def plot_diagnostic(hg:schism_grid, dp:np.ndarray): pass if __name__ == '__main__': - hg_file = Path('/sciclone/schism10/feiye/SECOFS/Inputs/I02e/Edit_bathy/Load_charts_NCF/hgrid_levee_loaded.gr3') + # -----------------inputs----------------------------- + hg_file = Path('/sciclone/schism10/feiye/SECOFS/Inputs/I02e/Edit_bathy/Load_charts_NCF/hgrid_dem_levee_loaded.gr3') + # -----------------end inputs----------------------------- + hg = schism_grid(str(hg_file)) dp_orig = hg.dp.copy() @@ -170,16 +135,8 @@ def plot_diagnostic(hg:schism_grid, dp:np.ndarray): 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 + hg.save(f"{hg_file.parent}/{hg_file.stem}_chart_loaded.gr3", fmt=1) 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")) - - 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")) + pass - pass \ No newline at end of file diff --git a/src/Utility/Pre-Processing/SECOFS/Bathy_edit/Levee/set_levees.py b/src/Utility/Pre-Processing/SECOFS/Bathy_edit/Levee/set_levees.py index 67d80545a..07779a293 100755 --- a/src/Utility/Pre-Processing/SECOFS/Bathy_edit/Levee/set_levees.py +++ b/src/Utility/Pre-Processing/SECOFS/Bathy_edit/Levee/set_levees.py @@ -179,8 +179,8 @@ def set_levees(hgrid_name='', gd:schism_grid=None): if __name__ == "__main__": wdir = './' - gd = set_levees(f'{wdir}/hgrid.ll.new') # input is the bathy-loaded hgrid - gd.save(f'{wdir}/hgrid_levee_loaded.gr3') + gd = set_levees(f'{wdir}/hgrid.ll.dem_loaded.gr3') # input is the bathy-loaded hgrid + gd.save(f'{wdir}/hgrid_dem_levee_loaded.gr3') diff --git a/src/Utility/Pre-Processing/SECOFS/Bathy_edit/NCF/load_NCF.py b/src/Utility/Pre-Processing/SECOFS/Bathy_edit/NCF/load_NCF.py new file mode 100755 index 000000000..0f73b816b --- /dev/null +++ b/src/Utility/Pre-Processing/SECOFS/Bathy_edit/NCF/load_NCF.py @@ -0,0 +1,75 @@ +import geopandas as gpd +import numpy as np +from pathlib import Path + +from matplotlib import pyplot as plt + +from pylib_essentials.schism_file import read_schism_hgrid_cached, grd2sms, schism_grid + + +def load_NCF(hg:schism_grid, NCF_shpfile:Path): + ''' + Load the maintained depth from National Channel Framework data into the hgrid. + The original NCF polygons are enlarged to accommodate the mismatch between the hgrid and the NCF data. + The maintained depth is converted from feet to meters. + ''' + + # get the nodes in the NCF (National Channel Framework) polygons + NCF_shpfile = Path('/sciclone/schism10/Hgrid_projects/Charts/Mississippi/channel_quarter_NCF.shp') + NCF_data = gpd.read_file(NCF_shpfile) + + # remove the polygons that are not in the bounding box of the hgrid + NCF_data = NCF_data.cx[hg.x.min():hg.x.max(), hg.y.min():hg.y.max()] + # convert any multipolygons to polygons + NCF_data = NCF_data.explode() + + # enlarge the polygons by 4 m + NCF_data['geometry'] = NCF_data['geometry'].to_crs('esri:102008').buffer(4).to_crs('epsg:4326') + # extract the polygons to a list of nx2 numpy arrays + # NCF_polygons = [np.array(poly.exterior.coords) for poly in NCF_data['geometry']] + + # put hgrid points into a Point GeoDataFrame + hg_points = gpd.GeoDataFrame(geometry=gpd.points_from_xy(hg.x, hg.y), crs='epsg:4326') + # determine which points are inside the polygons and get the maintained depth + joined_gdf = gpd.sjoin(hg_points, NCF_data, how="inner", predicate='within') + idx = joined_gdf.index.to_numpy() # get the indices of the points inside the polygons + dp_NCF = np.ones_like(hg.dp, dtype=float) * -9999 # initialize with a large negative number + 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.show() + + return dp_NCF + +def plot_diagnostic(hg:schism_grid, dp:np.ndarray): + # # plot to check the changed nodes + # plt.figure() + # idx = abs(dp_orig - hg.dp) > 1e-5 + # # add a 1:1 line for reference from the lower left to upper right + # # Determine the range for the 1:1 line + # min_val = min(min(dp_orig[idx]), min(hg.dp[idx])) + # max_val = max(max(dp_orig[idx]), max(hg.dp[idx])) + # # Add a 1:1 line for reference + # plt.plot([min_val, max_val], [min_val, max_val], 'k--') + # plt.scatter(dp_orig[idx], hg.dp[idx]) + # plt.axis('equal') + # plt.show() + pass + +if __name__ == '__main__': + # -----------------inputs----------------------------- + hg_file = Path('./hgrid_dem_levee_loaded.gr3') + # -----------------end inputs----------------------------- + + hg = schism_grid(str(hg_file)) + dp_orig = hg.dp.copy() + + # 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.save(f"{hg_file.parent}/{hg_file.stem}_NCF_loaded.gr3", fmt=1) + + pass diff --git a/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/Bathy_edit/DEM_loading/DEM_info.json b/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/Bathy_edit/DEM_loading/DEM_info.json index 0340035ac..8ed8c9529 100644 --- a/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/Bathy_edit/DEM_loading/DEM_info.json +++ b/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/Bathy_edit/DEM_loading/DEM_info.json @@ -113,5 +113,10 @@ "Annotation": "CuDEM (Continously updated DEM): Updated Florida", "Note": "https://vims0-my.sharepoint.com/personal/feiye_vims_edu/_layouts/15/Doc.aspx?sourcedoc={d6aa1259-6a28-4b12-9caf-58f05600b3b5}&action=edit&wd=target%28DEM.one%7Cae181e84-4618-4faf-8af2-32fe3a6c3d29%2FCuDEM%7C7ed05443-5395-4dee-81b2-a8f7b0d530f2%2F%29&wdorigin=703&wdpreservelink=1", "File": "/sciclone/schism10/lcui01/schism20/DEM/CUDEM/FL/" + }, + "BlueTopo": { + "Annotation": "BlueTopo", + "Note": "", + "File": "/sciclone/schism10/Hgrid_projects/DEMs/BlueTopo/" } } diff --git a/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/Bathy_edit/DEM_loading/pload_depth.py b/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/Bathy_edit/DEM_loading/pload_depth.py index 542e7ed2b..fc5601d6f 100755 --- a/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/Bathy_edit/DEM_loading/pload_depth.py +++ b/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/Bathy_edit/DEM_loading/pload_depth.py @@ -14,7 +14,7 @@ #Input #----------------------------------------------------------------------------- grd='hgrid.ll' #grid name -grdout='hgrid.ll.new' #grid name with depth loaded +grdout='hgrid.ll.dem_loaded.gr3' #grid name with depth loaded #parameter regions=[