Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

BMI prototype for Hyfeature network #598

Open
wants to merge 4 commits into
base: master
Choose a base branch
from
Open

BMI prototype for Hyfeature network #598

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
2afc900
BMI prototype for Hyfeature network
Dec 8, 2022
a89cdc0
BMI prototype for Hyfeature network: missing files
Dec 8, 2022
0e65edb
remove mislocated file
Dec 8, 2022
5a799d9
coastal boundary dataframe from hourly SCHISM data
Dec 16, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
630 changes: 630 additions & 0 deletions src/bmi_troute.py
View file
Open in desktop

Large diffs are not rendered by default.

96 changes: 88 additions & 8 deletions src/troute-network/troute/hyfeature_network_utilities.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@
from joblib import delayed, Parallel
import pyarrow as pa
import pyarrow.parquet as pq
import xarray as xr

import troute.nhd_io as nhd_io

Expand All @@ -23,11 +24,12 @@ def build_forcing_sets(
t0
):

run_sets = forcing_parameters.get("qlat_forcing_sets", None)
nexus_input_folder = forcing_parameters.get("nexus_input_folder", None)
nts = forcing_parameters.get("nts", None)
max_loop_size = forcing_parameters.get("max_loop_size", 12)
dt = forcing_parameters.get("dt", None)
run_sets = forcing_parameters.get("qlat_forcing_sets", None)
nexus_input_folder = forcing_parameters.get("nexus_input_folder", None)
downstream_boundary_input_folder = forcing_parameters.get("downstream_boundary_input_folder", None)
nts = forcing_parameters.get("nts", None)
max_loop_size = forcing_parameters.get("max_loop_size", 12)
dt = forcing_parameters.get("dt", None)

try:
nexus_input_folder = pathlib.Path(nexus_input_folder)
Expand All @@ -38,6 +40,7 @@ def build_forcing_sets(
raise AssertionError("Aborting simulation because the nexus_input_folder:", qlat_input_folder,"does not exist. Please check the the nexus_input_folder variable is correctly entered in the .yaml control file") from None

forcing_glob_filter = forcing_parameters.get("nexus_file_pattern_filter", "*.NEXOUT")
downstream_boundary_glob_filter = forcing_parameters.get("downstream_boundary_file_pattern_filter", "*SCHISM.nc")

if forcing_glob_filter=="nex-*":
print("Reformating qlat nexus files as hourly binary files...")
Expand All @@ -58,7 +61,7 @@ def build_forcing_sets(
nexus_input_folder, forcing_glob_filter = nex_files_to_binary(nexus_files_list, binary_folder)
forcing_parameters["nexus_input_folder"] = nexus_input_folder
forcing_parameters["nexus_file_pattern_filter"] = forcing_glob_filter

# TODO: Throw errors if insufficient input data are available
if run_sets:
#FIXME: Change it for hyfeature
Expand Down Expand Up @@ -136,7 +139,7 @@ def build_forcing_sets(

if k + max_loop_size < len(forcing_filename_list):
run_sets[j]['nexus_files'] = forcing_filename_list[k:k
+ max_loop_size]
+ max_loop_size]
else:
run_sets[j]['nexus_files'] = forcing_filename_list[k:]

Expand All @@ -161,6 +164,81 @@ def build_forcing_sets(
k += max_loop_size
j += 1

# creates downstream boundary forcing file list
if downstream_boundary_input_folder:
# get the first and seconded files from an ordered list of all forcing files
downstream_boundary_input_folder = pathlib.Path(downstream_boundary_input_folder)
all_files = sorted(downstream_boundary_input_folder.glob(downstream_boundary_glob_filter))
first_file = all_files[0]
second_file = all_files[1]

# Deduce the timeinterval of the forcing data from the output timestamps of the first
df = read_file(first_file)
t1_str = pd.to_datetime(df.time.iloc[0]).strftime("%Y-%m-%d_%H:%M:%S")
t1 = datetime.strptime(t1_str,"%Y-%m-%d_%H:%M:%S")
df = read_file(second_file)
t2_str = pd.to_datetime(df.time.iloc[0]).strftime("%Y-%m-%d_%H:%M:%S")
t2 = datetime.strptime(t2_str,"%Y-%m-%d_%H:%M:%S")
dt_downstream_boundary_timedelta = t2 - t1
dt_downstream_boundary = dt_downstream_boundary_timedelta.seconds

bts_subdivisions = dt_downstream_boundary / dt
#forcing_parameters['qts_subdivisions']= qts_subdivisions

# the number of files required for the simulation
nfiles = int(np.ceil(nts / bts_subdivisions))

# list of downstream boundary file datetimes
datetime_list = [t0 + dt_downstream_boundary_timedelta * (n) for n in
range(nfiles)]
datetime_list_str = [datetime.strftime(d, '%Y%m%d%H%M') for d in
datetime_list]
# list of downstream boundary files
downstream_boundary_filename_list = [d_str + downstream_boundary_glob_filter[1:] for d_str in
datetime_list_str]

# check that all forcing files exist
for f in downstream_boundary_filename_list:
try:
J = pathlib.Path(downstream_boundary_input_folder.joinpath(f))
assert J.is_file() == True
except AssertionError:
raise AssertionError("Aborting simulation because downstream boundary file", J, "cannot be not found.") from None

# build run sets list
#run_sets = []
k = 0
j = 0
nts_accum = 0
nts_last = 0
while k < len(downstream_boundary_filename_list):
#run_sets.append({})

if k + max_loop_size < len(downstream_boundary_filename_list):
run_sets[j]['downstream_boundary_files'] = downstream_boundary_filename_list[k:k
+ max_loop_size]
else:
run_sets[j]['downstream_boundary_files'] = downstream_boundary_filename_list[k:]

#nts_accum += len(run_sets[j]['nexus_files']) * qts_subdivisions
#if nts_accum <= nts:
# run_sets[j]['nts'] = int(len(run_sets[j]['nexus_files'])
# * qts_subdivisions)
#else:
# run_sets[j]['nts'] = int(nts - nts_last)

#final_nexout = nexus_input_folder.joinpath(run_sets[j]['nexus_files'
# ][-1])
#df = read_file(final_nexout)
#final_timestamp_str = pd.to_datetime(df.columns[1]).strftime("%Y-%m-%d_%H:%M:%S")

#run_sets[j]['final_timestamp'] = \
# datetime.strptime(final_timestamp_str, '%Y-%m-%d_%H:%M:%S')

#nts_last = nts_accum
k += max_loop_size
j += 1

return run_sets

def build_qlateral_array(
Expand Down Expand Up @@ -296,5 +374,7 @@ def read_file(file_name):
elif extension=='.parquet':
df = pq.read_table(file_name).to_pandas().reset_index()
df.index.name = None

elif extension=='.nc':
df = xr.open_dataset(file_name)
df = df.to_dataframe()
return df
36 changes: 22 additions & 14 deletions src/troute-network/troute/hyfeature_preprocess.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@
def build_hyfeature_network(supernetwork_parameters,
waterbody_parameters,
):

geo_file_path = supernetwork_parameters["geo_file_path"]
cols = supernetwork_parameters["columns"]
terminal_code = supernetwork_parameters.get("terminal_code", 0)
Expand Down Expand Up @@ -692,17 +692,17 @@ def hyfeature_forcing(
-----

"""

# Unpack user-specified forcing parameters
dt = forcing_parameters.get("dt", None)
qts_subdivisions = forcing_parameters.get("qts_subdivisions", None)
nexus_input_folder = forcing_parameters.get("nexus_input_folder", None)
qlat_file_index_col = forcing_parameters.get("qlat_file_index_col", "feature_id")
qlat_file_value_col = forcing_parameters.get("qlat_file_value_col", "q_lateral")
qlat_file_gw_bucket_flux_col = forcing_parameters.get("qlat_file_gw_bucket_flux_col", "qBucket")
qlat_file_terrain_runoff_col = forcing_parameters.get("qlat_file_terrain_runoff_col", "qSfcLatRunoff")

dt = forcing_parameters.get("dt", None)
qts_subdivisions = forcing_parameters.get("qts_subdivisions", None)
nexus_input_folder = forcing_parameters.get("nexus_input_folder", None)
qlat_file_index_col = forcing_parameters.get("qlat_file_index_col", "feature_id")
qlat_file_value_col = forcing_parameters.get("qlat_file_value_col", "q_lateral")
qlat_file_gw_bucket_flux_col = forcing_parameters.get("qlat_file_gw_bucket_flux_col", "qBucket")
qlat_file_terrain_runoff_col = forcing_parameters.get("qlat_file_terrain_runoff_col", "qSfcLatRunoff")
downstream_boundary_input_folder = forcing_parameters.get("downstream_boundary_input_folder", None)

# TODO: find a better way to deal with these defaults and overrides.
run["t0"] = run.get("t0", t0)
run["nts"] = run.get("nts")
Expand All @@ -712,7 +712,8 @@ def hyfeature_forcing(
run["qlat_file_index_col"] = run.get("qlat_file_index_col", qlat_file_index_col)
run["qlat_file_value_col"] = run.get("qlat_file_value_col", qlat_file_value_col)
run["qlat_file_gw_bucket_flux_col"] = run.get("qlat_file_gw_bucket_flux_col", qlat_file_gw_bucket_flux_col)
run["qlat_file_terrain_runoff_col"] = run.get("qlat_file_terrain_runoff_col", qlat_file_terrain_runoff_col)
run["qlat_file_terrain_runoff_col"] = run.get("qlat_file_terrain_runoff_col", qlat_file_terrain_runoff_col)
run["downstream_boundary_input_folder"] = run.get("downstream_boundary_input_folder", downstream_boundary_input_folder)

#---------------------------------------------------------------------------
# Assemble lateral inflow data
Expand Down Expand Up @@ -749,12 +750,19 @@ def hyfeature_forcing(
start_time = time.time()
LOG.info("creating coastal dataframe ...")

coastal_boundary_domain = nhd_io.read_coastal_boundary_domain(coastal_boundary_domain_files)
coastal_boundary_domain = nhd_io.read_coastal_boundary_domain(coastal_boundary_domain_files)
# create dataframe for hourly schism data
coastal_boundary_depth_df = nhd_io.build_coastal_dataframe(
run,
coastal_boundary_domain,
)
# create dataframe for multi hourly schism data
'''
coastal_boundary_depth_df = nhd_io.build_coastal_ncdf_dataframe(
coastal_boundary_elev_files,
coastal_boundary_domain,
)

'''
LOG.debug(
"coastal boundary elevation observation DataFrame creation complete in %s seconds." \
% (time.time() - start_time)
Expand Down
Loading