Skip to content

Commit

Permalink
STOFS-3D scripts: Minor updates on the driver script and tested on He…
Browse files Browse the repository at this point in the history
…rcules.
  • Loading branch information
feiye-vims committed Oct 29, 2024
1 parent c421d36 commit f859616
Show file tree
Hide file tree
Showing 4 changed files with 424 additions and 203 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -13,73 +13,85 @@
import geopandas as gpd

# ------------------------- inputs ---------------------------
wdir = '/sciclone/schism10/Hgrid_projects/STOFS3D-v8/v43s2_RiverMapper/v44/Clip/'
crs = 'esri:102008'
WDIR = '/sciclone/schism10/Hgrid_projects/STOFS3D-v8/v22/Clip/'
CRS = 'esri:102008'

# manual polygons defined in the coastal coverage
coastal_shpfile = f'{wdir}/inputs/coastal.shp' # esri:102008
# manual polygons defined in the coastal coverage,
# may need to delete some polygons so that they are treated as watershed
coastal_shpfile = f'{WDIR}/inputs/coastal.shp' # esri:102008

# land boundary + coastal, i.e., watershed region between the coastline and the land boundary + manual polygons in the coastal coverage
lbnd_coastal_shpfile = f'{wdir}/inputs/lbnd_coastal.shp' # esri:102008
# land boundary + coastal, i.e., watershed region between the coastline and the land boundary
# as well as manual polygons in the coastal coverage
lbnd_coastal_shpfile = f'{WDIR}/inputs/lbnd_coastal.shp' # esri:102008

# use linestring as the levee shapefile
levee_shpfile = f'{wdir}/inputs/levee.shp' # esri:102008
levee_buf_distance = 5
levee_shpfile = f'{WDIR}/inputs/levee.shp' # esri:102008
LEVEE_BUF_DISTANCE = 5

# used to select manual polygons that are supposedly better refined than the auto arcs
select_nearshore_shpfile = f'{wdir}/inputs/select_nearshore_v3.shp' # esri:102008
select_nearshore_shpfile = f'{WDIR}/inputs/select_nearshore_v3.shp' # esri:102008

auto_arcs_file = f'{wdir}/inputs/total_arcs.shp' # lonlat
auto_polys_file = f'{wdir}/inputs/total_river_polys.shp' # lonlat
auto_arcs_file = f'{WDIR}/inputs/total_arcs.shp' # lonlat
auto_polys_file = f'{WDIR}/inputs/total_river_polys.shp' # lonlat
# ------------------------- end inputs ---------------------------

# make outputs directory
output_dir = f'{wdir}/outputs'
output_dir = f'{WDIR}/outputs'
os.makedirs(output_dir, exist_ok=True)
# save a copy of this script to the working directory
os.system(f'cp {__file__} {WDIR}')

coastal = gpd.read_file(coastal_shpfile)
coastal.set_crs(crs, inplace=True)
coastal.set_crs(CRS, inplace=True)

lbnd_coastal = gpd.read_file(lbnd_coastal_shpfile)
lbnd_coastal.set_crs(crs, inplace=True)
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 is a rough estimate of watershed region,
# which excludes manual polygons but may include levee
# This includes most of the watershed region EXCEPT for the manual polygons in the coastal coverage
raw_watershed = lbnd_coastal.overlay(coastal, how='difference').dissolve()

# extract the refined coastal polygons, i.e., coastal minus those intersecting the select_nearshore polygons,
# which are the manual polygons that are supposedly better refined than the auto arcs
# (However this is not always the case because some manualy polygons are not accurate),
# where no intersections with auto arcs are allowed, hence the buffer
# Extract the refined coastal polygons, i.e.,
# coastal minus the un-refined polyogns (those intersecting the select_nearshore),
# which are the manual polygons that are supposedly better configured than the auto arcs.
# It also implies that refined coastal polygons accomodate the connectivity to watershed
# rivers, e.g., a lake shoreline should have denser nodes where it connects to a river.
# However, this is not always the case because some manualy polygons are not accurate,
# No intersections between refined polygons and auto arcs are allowed, hence the buffer.
# Unrefined polygons are allowed to intersect with auto arcs to ensure channel connectivity
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'{output_dir}/coastal_refined.shp', crs=crs)
refined_coastal_buf.to_file(f'{output_dir}/coastal_refined.shp')
print(f'coastal_refined.shp saved to {output_dir}')

# get final watershed by subtracting polygons that cannot be intersected by auto arcs from the raw watershed
# get final watershed
print('subtracting refined coastal from watershed')
watershed = raw_watershed.overlay(refined_coastal_buf, how='difference').dissolve()

# add additional watershed region (manually specified) to the watershed
pass

# subtract levee from watershed
levee_buf = gpd.GeoDataFrame(geometry=gpd.read_file(levee_shpfile).buffer(levee_buf_distance))
print('subtracting levee from watershed')
levee_buf = gpd.GeoDataFrame(geometry=gpd.read_file(levee_shpfile).buffer(LEVEE_BUF_DISTANCE))
levee_buf.to_file(f'{output_dir}/levee_buf.shp')
watershed = watershed.overlay(levee_buf, how='difference').dissolve()
watershed.to_file(f'{output_dir}/watershed.shp', crs=crs)
watershed.to_file(f'{output_dir}/watershed.shp')

# add additional watershed region (manually specified) to the watershed
print('break here if using qgis for clipping, which is faster')
watershed = gpd.read_file(f'{output_dir}/watershed.shp')

# clip the auto arcs to the watershed boundary
total_arcs = gpd.read_file(auto_arcs_file).to_crs(crs)
total_arcs = gpd.read_file(auto_arcs_file).to_crs(CRS)
total_arcs_clipped = total_arcs.clip(watershed)
total_arcs_clipped.to_file(f'{output_dir}/total_arcs_clipped.shp', crs=crs)
total_arcs_clipped.to_file(f'{output_dir}/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(auto_polys_file).to_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(auto_polys_file).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'{output_dir}/total_river_polys_clipped_test.shp', crs=crs)

pass
total_arcs_polygons_clipped.to_file(f'{output_dir}/total_river_polys_clipped_test.shp', crs=CRS)
print('done')
Original file line number Diff line number Diff line change
@@ -1,37 +1,59 @@
import os
'''
Generate tvd.prop file for the high-order transport schemes in SCHISM
'''
import numpy as np

from pylib import schism_grid, read_schism_bpfile, inside_polygon
import socket
import pylib
from pylib import read_schism_bpfile, inside_polygon
if 'sciclone' in socket.gethostname():
from pylib_experimental.schism_file import cread_schism_hgrid as schism_read
else:
from pylib import schism_grid as schism_read

if __name__ == '__main__':

#Read hgrid
hgrid_name='hgrid.gr3'
gr3_pkl = 'hgrid.pkl'
if os.path.exists(gr3_pkl):
print('Reading hgrid.pkl')
gd = schism_grid(gr3_pkl)
else:
print('Reading hgrid.gr3')
gd = schism_grid(hgrid_name)
gd.save(gr3_pkl)

def gen_tvd_prop(hg: pylib.schism_grid, regions: list):
'''
Generate tvd.prop file for the high-order transport schemes in SCHISM
#gd = read_schism_hgrid('hgrid.gr3')
gd.compute_ctr()
Parameters
----------
hg : schism_grid
SCHISM grid object
regions : list
list of region files,
in the final output tvd.prop: 1 for inside region, 0 for outside region
regions = ['tvd0_1.reg', 'tvd0_2.reg', 'tvd0_3.reg', 'tvd0_4.reg', \
'tvd0_5.reg', 'tvd0_6.reg', 'tvd0_7.reg', \
'upwind_east_Carribbean.rgn', 'upwind_west_Carribbean.rgn']
Returns
-------
None
#write constant (=1) *prop
ab = np.ones(gd.ne)
'''
hg.compute_ctr()
tvd = np.zeros(hg.ne)

#reset values (=0) within region
for region in regions:
bp = read_schism_bpfile(region, fmt=1)
sind = inside_polygon(np.c_[gd.xctr,gd.yctr], bp.x, bp.y)
ab[sind==1] = 0
idx = inside_polygon(np.c_[hg.xctr, hg.yctr], bp.x, bp.y)
tvd[idx == 1] = 1

hg.write_prop('tvd.prop', value=tvd, fmt='{:d}')


gd.write_prop('tvd.prop',value=ab, fmt='{:d}')
def sample_usage():
'''
Sample usage of gen_tvd_prop
'''
gd = schism_read('hgrid.gr3')

regions = [
'tvd0_1.reg', 'tvd0_2.reg', 'tvd0_3.reg', 'tvd0_4.reg',
'tvd0_5.reg', 'tvd0_6.reg', 'tvd0_7.reg',
'upwind_east_Carribbean.rgn', 'upwind_west_Carribbean.rgn'
]

gen_tvd_prop(gd, regions)


if __name__ == '__main__':
sample_usage()
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,8 @@
import numpy as np
from scipy import spatial

from pylib_essentials.schism_file import source_sink, TimeHistory, read_schism_hgrid_cached, schism_bpfile
from pylib import schism_grid
from pylib_essentials.schism_file import source_sink, TimeHistory, schism_bpfile
from pylib_essentials.utility_functions import inside_polygon

# Global Var
Expand Down Expand Up @@ -161,7 +162,7 @@ def relocate_sources2(
# -------------------------------------info from old source_sink-------------------------------------
# Read original source/sink generated by pyschism based on NWM
old_source_sink = source_sink.from_files(source_dir=old_ss_dir)
old_gd = read_schism_hgrid_cached(f'{old_ss_dir}/hgrid.gr3', overwrite_cache=False)
old_gd = schism_grid(f'{old_ss_dir}/hgrid.gr3')
old_gd.compute_ctr() # computer element center coordinates

# read source to fid mapping from json
Expand Down Expand Up @@ -197,7 +198,7 @@ def relocate_sources2(

# -------------------------------------new hgrid-------------------------------------
# get new hgrid (with feeders)
new_gd = read_schism_hgrid_cached(hgrid_fname, overwrite_cache=False)
new_gd = schism_grid(hgrid_fname)
new_gd.compute_ctr()

# -------------------- process nan values in mandatory_sources_coor ------------------
Expand Down Expand Up @@ -399,7 +400,7 @@ def relocate_sources(

# -------------------------------------info from old source_sink-------------------------------------
# get old hgrid (without feeders); hgrid with feeders may also be used
old_gd = read_schism_hgrid_cached(f'{old_ss_dir}/hgrid.gr3', overwrite_cache=False)
old_gd = schism_grid(f'{old_ss_dir}/hgrid.gr3')

# get old source coordinates
old_gd.compute_ctr()
Expand All @@ -424,7 +425,7 @@ def relocate_sources(

# -------------------------------------new hgrid-------------------------------------
# get new hgrid (with feeders)
new_gd = read_schism_hgrid_cached(hgrid_fname, overwrite_cache=False)
new_gd = schism_grid(hgrid_fname)
new_gd.compute_ctr()

# -------------------------------------manual intervetion as pre-requisits -------------------------------------
Expand Down Expand Up @@ -609,7 +610,7 @@ def samples():
# copy data files from git to the working directory for record if needed

# read hgrid
hgrid = read_schism_hgrid_cached(hgrid_fname, overwrite_cache=False)
hgrid = schism_grid(hgrid_fname)

# read region
region = schism_bpfile()
Expand Down
Loading

0 comments on commit f859616

Please sign in to comment.