From 7f54cc8048487ef4b82b3f106edc6af492e8ff4c Mon Sep 17 00:00:00 2001 From: sclaw Date: Mon, 10 Jun 2024 16:30:20 -0400 Subject: [PATCH 01/43] first draft moving water around --- src/muskingumcunge/network.py | 60 +++++++++++++++++++++++++++++++++++ 1 file changed, 60 insertions(+) create mode 100644 src/muskingumcunge/network.py diff --git a/src/muskingumcunge/network.py b/src/muskingumcunge/network.py new file mode 100644 index 0000000..ff1c72c --- /dev/null +++ b/src/muskingumcunge/network.py @@ -0,0 +1,60 @@ +from collections import defaultdict +import pandas as pd +import numpy as np + + +class Network: + + def __init__(self, edge_dict): + self.edge_dict = edge_dict # key=u/s, value=d/s + self.chlid_dict = defaultdict(list) + for u, v in edge_dict.items(): + self.chlid_dict[v].append(u) + self.headwaters = list(set(edge_dict.keys()).difference(set(edge_dict.values()))) + self.root = set(edge_dict.values()).difference(set(edge_dict.keys())).pop() + self.post_order = list() + q = [self.root] + working = True + while working: + cur = q.pop() + children = self.chlid_dict[cur] + if len(children) == 0: + self.post_order.append(cur) + elif all([c in self.post_order for c in children]): + self.post_order.append(cur) + else: + q.append(cur) + q.extend(children) + if len(q) == 0: + working = False + self.post_order = self.post_order[:-1] + + self.reach_dict = None + self.forcing_df = None + self.channel_outflows = {k: None for k in self.edge_dict.keys()} + + def export_post_order(self, path): + out_df = pd.DataFrame(self.post_order, columns=['comid']) + out_df['order'] = range(1, len(out_df) + 1) + out_df.to_csv(path, index=False) + + def load_reaches(self, reach_dict): + self.reach_dict = reach_dict + + def load_forcings(self, forcing_path): + self.forcing_df = pd.read_csv(forcing_path, index_col=0, parse_dates=True) + for n in self.headwaters: + self.channel_outflows[n] = np.zeros(self.forcing_df.shape[0]) + + def run_event(self): + for node in self.post_order: + reach = self.reach_dict[node] + children = self.chlid_dict[node] + us_hydro = [self.channel_outflows[c] for c in children] + us_hydro = np.sum(us_hydro, axis=0) + # routed = reach.route_hydrograph(us_hydro) + routed = us_hydro + outflow = routed + self.forcing_df[node] + self.channel_outflows[node] = outflow + + From bbf2aa9357f575a7e541781588ec6e2cd5fd8f46 Mon Sep 17 00:00:00 2001 From: sclaw Date: Tue, 11 Jun 2024 14:32:56 -0400 Subject: [PATCH 02/43] functional --- src/muskingumcunge/network.py | 41 +++++++++++++++++++++++++++++------ 1 file changed, 34 insertions(+), 7 deletions(-) diff --git a/src/muskingumcunge/network.py b/src/muskingumcunge/network.py index ff1c72c..ec0b966 100644 --- a/src/muskingumcunge/network.py +++ b/src/muskingumcunge/network.py @@ -11,6 +11,9 @@ def __init__(self, edge_dict): for u, v in edge_dict.items(): self.chlid_dict[v].append(u) self.headwaters = list(set(edge_dict.keys()).difference(set(edge_dict.values()))) + roots = set(edge_dict.values()).difference(set(edge_dict.keys())) + if len(roots) > 1: + raise ValueError(f'Multiple roots detected: {roots}') self.root = set(edge_dict.values()).difference(set(edge_dict.keys())).pop() self.post_order = list() q = [self.root] @@ -32,6 +35,7 @@ def __init__(self, edge_dict): self.reach_dict = None self.forcing_df = None self.channel_outflows = {k: None for k in self.edge_dict.keys()} + self.out_df = None def export_post_order(self, path): out_df = pd.DataFrame(self.post_order, columns=['comid']) @@ -41,20 +45,43 @@ def export_post_order(self, path): def load_reaches(self, reach_dict): self.reach_dict = reach_dict - def load_forcings(self, forcing_path): - self.forcing_df = pd.read_csv(forcing_path, index_col=0, parse_dates=True) + def load_forcings(self, forcing_df): + self.forcing_df = forcing_df + self.forcing_df = self.forcing_df[self.forcing_df.columns.intersection(self.post_order)].copy() for n in self.headwaters: self.channel_outflows[n] = np.zeros(self.forcing_df.shape[0]) + + def load_headwater_forcings(self, headwater_forcings): + for n in self.headwaters: + self.channel_outflows[n] = headwater_forcings[n].to_numpy() def run_event(self): + # calculate total inflow + volume_in = self.forcing_df.sum().sum() + np.sum([self.channel_outflows[c].sum() for c in self.headwaters]) + dt = (self.forcing_df.index[1] - self.forcing_df.index[0]).seconds / 3600 for node in self.post_order: reach = self.reach_dict[node] - children = self.chlid_dict[node] - us_hydro = [self.channel_outflows[c] for c in children] - us_hydro = np.sum(us_hydro, axis=0) - # routed = reach.route_hydrograph(us_hydro) - routed = us_hydro + if node in self.headwaters: + us_hydro = self.channel_outflows[node] + else: + children = self.chlid_dict[node] + us_hydro = [self.channel_outflows[c] for c in children] + us_hydro = np.sum(us_hydro, axis=0) + routed = reach.route_hydrograph(us_hydro, dt) + # routed = us_hydro outflow = routed + self.forcing_df[node] self.channel_outflows[node] = outflow + + # calculate total outflow + us_root = self.chlid_dict[self.root] + us_hydro = [self.channel_outflows[c] for c in us_root] + self.channel_outflows[self.root] = np.sum(us_hydro, axis=0) + volume_out = self.channel_outflows[self.root].sum() + print(f"Volume conservation = {round(100 *volume_out / volume_in, 2)} %") + + # Post-processing + out_df = pd.DataFrame(self.channel_outflows) + out_df['dt'] = self.forcing_df.index + self.out_df = out_df.set_index('dt') From b69ca4cc0a81fe6fac4bd9d26a474b5493653f06 Mon Sep 17 00:00:00 2001 From: sclaw Date: Wed, 12 Jun 2024 11:58:32 -0400 Subject: [PATCH 03/43] add routing instead of dry run --- src/muskingumcunge/network.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/src/muskingumcunge/network.py b/src/muskingumcunge/network.py index ec0b966..f6f6ef7 100644 --- a/src/muskingumcunge/network.py +++ b/src/muskingumcunge/network.py @@ -67,8 +67,16 @@ def run_event(self): children = self.chlid_dict[node] us_hydro = [self.channel_outflows[c] for c in children] us_hydro = np.sum(us_hydro, axis=0) - routed = reach.route_hydrograph(us_hydro, dt) - # routed = us_hydro + + dx, subreaches = reach.optimize_route_params(us_hydro, dt) + if subreaches > 1: + routed = us_hydro.copy() + reach.reach_length = dx + for i in range(subreaches): + routed = reach.route_hydrograph(routed, dt) + else: + routed = reach.route_hydrograph(us_hydro, dt) + outflow = routed + self.forcing_df[node] self.channel_outflows[node] = outflow From 065dc3d166880be31ec76c93d46aa7b4109dc231 Mon Sep 17 00:00:00 2001 From: sclaw Date: Wed, 12 Jun 2024 12:00:06 -0400 Subject: [PATCH 04/43] add method for auto dx, fixed dt --- src/muskingumcunge/reach.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/src/muskingumcunge/reach.py b/src/muskingumcunge/reach.py index 66f18ef..d6b46b9 100644 --- a/src/muskingumcunge/reach.py +++ b/src/muskingumcunge/reach.py @@ -146,6 +146,23 @@ def route_hydrograph_mct(self, inflows, dt, max_iter=1000): reynoldref1 = reynoldref2 return np.array(outflows) + + def optimize_route_params(self, inflows, dt): + # Ponce method + qref = (inflows.max() + inflows.min()) / 2 + cref = np.interp(qref, self.geometry['discharge'], self.geometry['celerity']) + twref = np.interp(qref, self.geometry['discharge'], self.geometry['top_width']) + dxc = dt * 60 * 60 * cref # Courant length + dxd = (qref / twref) / (self.slope * cref) # characteristic reach length + dxmax = 0.5 * (dxc + dxd) + peak_loc = np.argmax(self.geometry['discharge'] > inflows.max()) + peak_loc = max([peak_loc, 1]) + cmax = self.geometry['celerity'][:peak_loc].max() # I think the fact that this is the max should cover us for making all stable + dxmin = cmax * (dt * 60 * 60) + dx = max([dxmin, dxmax]) + subreaches = int(np.ceil(self.reach_length / dx)) + dx = self.reach_length / subreaches + return dx, subreaches class TrapezoidalReach(BaseReach): From ae43e984e562e52dd97b1054d2e7f82401a00654 Mon Sep 17 00:00:00 2001 From: sclaw Date: Wed, 12 Jun 2024 16:26:33 -0400 Subject: [PATCH 05/43] send laterals to MC routing. update method flow for readability. add conservation of mass option --- src/muskingumcunge/network.py | 25 ++++++++++++++++--------- 1 file changed, 16 insertions(+), 9 deletions(-) diff --git a/src/muskingumcunge/network.py b/src/muskingumcunge/network.py index f6f6ef7..6686c36 100644 --- a/src/muskingumcunge/network.py +++ b/src/muskingumcunge/network.py @@ -55,7 +55,7 @@ def load_headwater_forcings(self, headwater_forcings): for n in self.headwaters: self.channel_outflows[n] = headwater_forcings[n].to_numpy() - def run_event(self): + def run_event(self, optimize_dx=True, conserve_mass=False): # calculate total inflow volume_in = self.forcing_df.sum().sum() + np.sum([self.channel_outflows[c].sum() for c in self.headwaters]) dt = (self.forcing_df.index[1] - self.forcing_df.index[0]).seconds / 3600 @@ -68,17 +68,24 @@ def run_event(self): us_hydro = [self.channel_outflows[c] for c in children] us_hydro = np.sum(us_hydro, axis=0) - dx, subreaches = reach.optimize_route_params(us_hydro, dt) - if subreaches > 1: - routed = us_hydro.copy() + laterals = self.forcing_df[node].to_numpy() + + if optimize_dx: + dx, subreaches = reach.optimize_route_params(us_hydro, dt) reach.reach_length = dx - for i in range(subreaches): - routed = reach.route_hydrograph(routed, dt) else: - routed = reach.route_hydrograph(us_hydro, dt) + subreaches = 1 + + routed = us_hydro + for i in range(subreaches): + routed = reach.route_hydrograph(routed, dt, lateral=laterals) - outflow = routed + self.forcing_df[node] - self.channel_outflows[node] = outflow + # enforce mass conservation + if conserve_mass: + in_flow = us_hydro.sum() + laterals.sum() + out_flow = routed.sum() + routed = routed * (in_flow / out_flow) + self.channel_outflows[node] = routed # calculate total outflow us_root = self.chlid_dict[self.root] From d66f9867bc933c96cae8fd06f13dc322b533b007 Mon Sep 17 00:00:00 2001 From: sclaw Date: Wed, 12 Jun 2024 16:27:13 -0400 Subject: [PATCH 06/43] add correct processing of laterals --- src/muskingumcunge/reach.py | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/src/muskingumcunge/reach.py b/src/muskingumcunge/reach.py index d6b46b9..5416a15 100644 --- a/src/muskingumcunge/reach.py +++ b/src/muskingumcunge/reach.py @@ -56,17 +56,21 @@ def route_hydrograph_c(self, inflows, dt): return croute(inflows, dt, reach_length, slope, geometry) - def route_hydrograph(self, inflows, dt, max_iter=1000): + def route_hydrograph(self, inflows, dt, lateral=None, max_iter=1000): outflows = list() outflows.append((inflows[0])) assert max(inflows) < max(self.geometry['discharge']), 'Rating Curve does not cover range of flowrates in hydrograph' + + if lateral is None: + lateral = np.zeros_like(inflows) + lateral = (lateral[:-1] + lateral[1:]) / 2 + lateral = np.append(lateral, 0) for i in range(len(inflows) - 1): - q_guess = sum([inflows[i], inflows[i + 1], outflows[i]]) / 3 + q_guess = sum([inflows[i], inflows[i + 1], outflows[i], lateral[i]]) / 3 last_guess = q_guess * 2 counter = 1 - # while abs(last_guess - q_guess) > 0.003: # from handbook of hydrology page 328 - while counter < 2: + while abs(last_guess - q_guess) > 0.003: # from handbook of hydrology page 328 counter += 1 last_guess = q_guess.copy() reach_q = sum([inflows[i], inflows[i + 1], outflows[i], q_guess]) / 4 @@ -82,11 +86,9 @@ def route_hydrograph(self, inflows, dt, max_iter=1000): c0 = (-1 + courant + reynold) / (1 + courant + reynold) c1 = (1 + courant - reynold) / (1 + courant + reynold) c2 = (1 - courant + reynold) / (1 + courant + reynold) + c3 = (2 * courant) / (1 + courant + reynold) - k = self.reach_length / c_tmp - x = 0.5 * (1 - reynold) - - q_guess = (c0 * inflows[i + 1]) + (c1 * inflows[i]) + (c2 * outflows[i]) + q_guess = (c0 * inflows[i + 1]) + (c1 * inflows[i]) + (c2 * outflows[i]) + (c3 * lateral[i]) q_guess = max(min(inflows), q_guess) if counter == max_iter: last_guess = q_guess From 54781c14f96e5207dbd5e9c5169b55ce287736c5 Mon Sep 17 00:00:00 2001 From: sclaw Date: Tue, 18 Jun 2024 16:21:26 -0400 Subject: [PATCH 07/43] add option to choose where laterals enter --- src/muskingumcunge/network.py | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/src/muskingumcunge/network.py b/src/muskingumcunge/network.py index 6686c36..6e1591f 100644 --- a/src/muskingumcunge/network.py +++ b/src/muskingumcunge/network.py @@ -55,7 +55,7 @@ def load_headwater_forcings(self, headwater_forcings): for n in self.headwaters: self.channel_outflows[n] = headwater_forcings[n].to_numpy() - def run_event(self, optimize_dx=True, conserve_mass=False): + def run_event(self, optimize_dx=True, conserve_mass=False, lat_addition='middle'): # calculate total inflow volume_in = self.forcing_df.sum().sum() + np.sum([self.channel_outflows[c].sum() for c in self.headwaters]) dt = (self.forcing_df.index[1] - self.forcing_df.index[0]).seconds / 3600 @@ -69,6 +69,8 @@ def run_event(self, optimize_dx=True, conserve_mass=False): us_hydro = np.sum(us_hydro, axis=0) laterals = self.forcing_df[node].to_numpy() + if lat_addition == 'top': + us_hydro = us_hydro + laterals if optimize_dx: dx, subreaches = reach.optimize_route_params(us_hydro, dt) @@ -78,7 +80,17 @@ def run_event(self, optimize_dx=True, conserve_mass=False): routed = us_hydro for i in range(subreaches): - routed = reach.route_hydrograph(routed, dt, lateral=laterals) + try: + if lat_addition == 'middle': + l = laterals + else: + l = None + routed = reach.route_hydrograph(routed, dt, lateral=l) + except AssertionError as e: + print(f"Error routing {node}: {e}") + + if lat_addition == 'bottom': + routed = routed + laterals # enforce mass conservation if conserve_mass: From 3acf156bdba95f23a40947b63ce382c6b8dc28d4 Mon Sep 17 00:00:00 2001 From: sclaw Date: Thu, 20 Jun 2024 10:23:44 -0400 Subject: [PATCH 08/43] incremental --- src/muskingumcunge/network.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/muskingumcunge/network.py b/src/muskingumcunge/network.py index 6e1591f..f841d8c 100644 --- a/src/muskingumcunge/network.py +++ b/src/muskingumcunge/network.py @@ -30,7 +30,7 @@ def __init__(self, edge_dict): q.extend(children) if len(q) == 0: working = False - self.post_order = self.post_order[:-1] + self.post_order = self.post_order[:-1] # remove root self.reach_dict = None self.forcing_df = None @@ -48,18 +48,21 @@ def load_reaches(self, reach_dict): def load_forcings(self, forcing_df): self.forcing_df = forcing_df self.forcing_df = self.forcing_df[self.forcing_df.columns.intersection(self.post_order)].copy() - for n in self.headwaters: - self.channel_outflows[n] = np.zeros(self.forcing_df.shape[0]) def load_headwater_forcings(self, headwater_forcings): for n in self.headwaters: self.channel_outflows[n] = headwater_forcings[n].to_numpy() + if n in self.post_order: + self.post_order.remove(n) def run_event(self, optimize_dx=True, conserve_mass=False, lat_addition='middle'): # calculate total inflow volume_in = self.forcing_df.sum().sum() + np.sum([self.channel_outflows[c].sum() for c in self.headwaters]) dt = (self.forcing_df.index[1] - self.forcing_df.index[0]).seconds / 3600 + counter = 1 for node in self.post_order: + if counter % 100 == 0: + print(f"{counter} / {len(self.post_order)}") reach = self.reach_dict[node] if node in self.headwaters: us_hydro = self.channel_outflows[node] From 9cbf3a4cf43b04d67fdf08a61f06e4aa18a3c5e8 Mon Sep 17 00:00:00 2001 From: sclaw Date: Thu, 20 Jun 2024 11:26:00 -0400 Subject: [PATCH 09/43] update headwater handling --- src/muskingumcunge/network.py | 31 ++++++++++++++++++++----------- 1 file changed, 20 insertions(+), 11 deletions(-) diff --git a/src/muskingumcunge/network.py b/src/muskingumcunge/network.py index f841d8c..fc37563 100644 --- a/src/muskingumcunge/network.py +++ b/src/muskingumcunge/network.py @@ -10,11 +10,21 @@ def __init__(self, edge_dict): self.chlid_dict = defaultdict(list) for u, v in edge_dict.items(): self.chlid_dict[v].append(u) - self.headwaters = list(set(edge_dict.keys()).difference(set(edge_dict.values()))) + roots = set(edge_dict.values()).difference(set(edge_dict.keys())) if len(roots) > 1: raise ValueError(f'Multiple roots detected: {roots}') self.root = set(edge_dict.values()).difference(set(edge_dict.keys())).pop() + self.headwaters = list() + self.post_order = list() + self.calculate_post_order() + self.reach_dict = None + self.forcing_df = None + self.channel_outflows = {k: None for k in self.edge_dict.keys()} + self.out_df = None + + def calculate_post_order(self): + self.headwaters = list() self.post_order = list() q = [self.root] working = True @@ -32,28 +42,27 @@ def __init__(self, edge_dict): working = False self.post_order = self.post_order[:-1] # remove root - self.reach_dict = None - self.forcing_df = None - self.channel_outflows = {k: None for k in self.edge_dict.keys()} - self.out_df = None - def export_post_order(self, path): out_df = pd.DataFrame(self.post_order, columns=['comid']) out_df['order'] = range(1, len(out_df) + 1) out_df.to_csv(path, index=False) def load_reaches(self, reach_dict): - self.reach_dict = reach_dict + self.reach_dict = reach_dict # Todo: add informational note about whether all reaches in post order have been loaded def load_forcings(self, forcing_df): self.forcing_df = forcing_df self.forcing_df = self.forcing_df[self.forcing_df.columns.intersection(self.post_order)].copy() def load_headwater_forcings(self, headwater_forcings): - for n in self.headwaters: + headwater_forcings = headwater_forcings.fillna(0) + missing = list(set(self.headwaters).difference(set(headwater_forcings.columns))) + print(f"Missing headwater forcings for: {missing}") + for n in headwater_forcings.columns: self.channel_outflows[n] = headwater_forcings[n].to_numpy() - if n in self.post_order: - self.post_order.remove(n) + self.headwaters.append(n) + for n in missing: + self.channel_outflows[n] = np.zeros(len(self.forcing_df)) def run_event(self, optimize_dx=True, conserve_mass=False, lat_addition='middle'): # calculate total inflow @@ -65,7 +74,7 @@ def run_event(self, optimize_dx=True, conserve_mass=False, lat_addition='middle' print(f"{counter} / {len(self.post_order)}") reach = self.reach_dict[node] if node in self.headwaters: - us_hydro = self.channel_outflows[node] + continue else: children = self.chlid_dict[node] us_hydro = [self.channel_outflows[c] for c in children] From d966e3d2423b8d1edd598395e7976910d4c8f550 Mon Sep 17 00:00:00 2001 From: sclaw Date: Thu, 20 Jun 2024 11:47:21 -0400 Subject: [PATCH 10/43] add lake forcings --- src/muskingumcunge/network.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/muskingumcunge/network.py b/src/muskingumcunge/network.py index fc37563..b7a3893 100644 --- a/src/muskingumcunge/network.py +++ b/src/muskingumcunge/network.py @@ -64,6 +64,12 @@ def load_headwater_forcings(self, headwater_forcings): for n in missing: self.channel_outflows[n] = np.zeros(len(self.forcing_df)) + def load_lake_forcings(self, lake_forcings): + lake_forcings = lake_forcings.fillna(0) + for n in lake_forcings.columns: + self.channel_outflows[n] = lake_forcings[n].to_numpy() + self.headwaters.append(n) + def run_event(self, optimize_dx=True, conserve_mass=False, lat_addition='middle'): # calculate total inflow volume_in = self.forcing_df.sum().sum() + np.sum([self.channel_outflows[c].sum() for c in self.headwaters]) From cf567e56c9c2f15fa4bb59ff7d6931bc53fc28f1 Mon Sep 17 00:00:00 2001 From: sclaw Date: Thu, 20 Jun 2024 17:55:40 -0400 Subject: [PATCH 11/43] match wrf compound to nwm wrf hydro github --- src/muskingumcunge/reach.py | 83 ++++++++++++++++++++----------------- 1 file changed, 46 insertions(+), 37 deletions(-) diff --git a/src/muskingumcunge/reach.py b/src/muskingumcunge/reach.py index 5416a15..25b11f6 100644 --- a/src/muskingumcunge/reach.py +++ b/src/muskingumcunge/reach.py @@ -91,6 +91,7 @@ def route_hydrograph(self, inflows, dt, lateral=None, max_iter=1000): q_guess = (c0 * inflows[i + 1]) + (c1 * inflows[i]) + (c2 * outflows[i]) + (c3 * lateral[i]) q_guess = max(min(inflows), q_guess) if counter == max_iter: + print('hit max iter') last_guess = q_guess outflows.append(q_guess) @@ -241,6 +242,13 @@ def generate_geometry(self): da = geom['area'][1:] - geom['area'][:-1] dq_da = dq / da geom['celerity'] = np.append(dq_da, dq_da[-1]) + + # dpdy = (geom['wetted_perimeter'][1:] - geom['wetted_perimeter'][:-1]) / (geom['stage'][1:] - geom['stage'][:-1]) + # dpdy = np.append(dpdy, dpdy[-1]) + # dpdy[dpdy < 0.0001] = 0.0001 + # k_prime = (5 / 3) - ((2 / 3) * (geom['area'] / (geom['top_width'] * geom['wetted_perimeter'])) * dpdy) + # geom['celerity'] = k_prime * (geom['discharge'] / geom['area']) + geom['celerity'][geom['celerity'] < 0.0001] = 0.0001 geom['log_q'] = np.log(geom['discharge']) @@ -340,7 +348,7 @@ def __init__(self, k, x): self.k = k self.x = x - def route_hydrograph(self, inflows, dt): + def route_hydrograph(self, inflows, dt, lateral=None): outflows = list() outflows.append((inflows[0])) @@ -359,13 +367,13 @@ class WRFCompound(BaseReach): """ def __init__(self, bottom_width, side_slope, bankfull_depth, floodplain_width, channel_n, floodplain_n, slope, reach_length, max_stage=10, stage_resolution=50): - self.bottom_width = bottom_width - self.side_slope = side_slope - self.bankfull_depth = bankfull_depth - assert floodplain_width >= self.bottom_width + (self.side_slope * self.bankfull_depth), 'floodplain width smaller than channel width at bankfull depth' - self.floodplain_width = floodplain_width - self.channel_n = channel_n - self.floodplain_n = floodplain_n + self.Bw = bottom_width + self.z = side_slope # Should be rise/run per side + self.bfd = bankfull_depth + self.TwCC = floodplain_width + self.n = channel_n + self.nCC = floodplain_n + self.So = slope self.slope = slope self.reach_length = reach_length self.max_stage = max_stage @@ -377,37 +385,38 @@ def generate_geometry(self): geom = self.geometry geom['stage'] = np.linspace(0, self.max_stage, self.resolution) - geom['top_width'] = self.bottom_width + (2 * geom['stage'] * self.side_slope) - geom['area'] = ((geom['top_width'] + self.bottom_width) / 2) * geom['stage'] - geom['wetted_perimeter'] = (((((geom['stage'] * self.side_slope) ** 2) + (geom['stage'] ** 2)) ** 0.5) * 2) + self.bottom_width - geom['hydraulic_radius'] = geom['area'] / geom['wetted_perimeter'] - geom['hydraulic_radius'][0] = geom['hydraulic_radius'][1] - geom['discharge'] = (1 / self.channel_n) * geom['area'] * (geom['hydraulic_radius'] ** (2 / 3)) * (self.slope ** 0.5) - geom['dpdy'] = np.repeat(np.sqrt((self.side_slope ** 2) + 1) * 2, geom['wetted_perimeter'].shape[0]) - - # Add compound channel - bankfull_indices = (geom['stage'] > self.bankfull_depth) - tw_at_bkf = self.bottom_width + (2 * self.bankfull_depth * self.side_slope) - a_at_bkf = ((self.bottom_width + tw_at_bkf) / 2) * self.bankfull_depth - p_at_bkf = (((((self.bankfull_depth * self.side_slope) ** 2) + (self.bankfull_depth ** 2)) ** 0.5) * 2) + self.bottom_width - area = a_at_bkf + (self.floodplain_width * (geom['stage'] - self.bankfull_depth)) - perimeter = p_at_bkf + (self.floodplain_width - tw_at_bkf) + (geom['stage'] - self.bankfull_depth) - radius = area / perimeter - radius[~bankfull_indices] = 1 - n_adj = ((self.channel_n * p_at_bkf) + ((perimeter - p_at_bkf) * self.floodplain_n)) / perimeter - discharge = (1 / n_adj) * area * (radius ** (2 / 3)) * (self.slope ** 0.5) - - geom['wetted_perimeter'][bankfull_indices] = perimeter[bankfull_indices] - geom['area'][bankfull_indices] = area[bankfull_indices] - geom['top_width'][bankfull_indices] = self.floodplain_width - geom['discharge'][bankfull_indices] = discharge[bankfull_indices] - geom['dpdy'][bankfull_indices] = 2 - k_prime = (5 / 3) - ((2 / 3) * (geom['area'] / (geom['top_width'] * geom['wetted_perimeter'])) * geom['dpdy']) - geom['celerity'] = k_prime * (geom['discharge'] / geom['area']) - geom['celerity'][0] = geom['celerity'][1] + # Trapezoidal Channel + geom['area'] = (self.Bw + geom['stage'] * self.z) * geom['stage'] + geom['wetted_perimeter'] = (self.Bw + 2.0 * geom['stage'] * np.sqrt(1.0 + self.z*self.z)) + geom['hydraulic_radius'] = geom['area'] / geom['wetted_perimeter'] + read_celerity = lambda y, r, n: (np.sqrt(self.So)/n)*((5./3.)*r**(2./3.) - ((2./3.)*r**(5./3.)*(2.0*np.sqrt(1.0 + self.z*self.z)/(self.Bw+2.0*y*self.z)))) + geom['celerity'] = read_celerity(geom['stage'], geom['hydraulic_radius'], self.n) + geom['mannings_n'] = np.repeat(self.n, self.resolution) + geom['top_width'] = self.Bw + (2 * self.z * geom['stage']) + + # Compound Channel + mask = (geom['stage'] > self.bfd) + area_trap = (self.Bw + self.bfd * self.z) * self.bfd + wetted_trap = (self.Bw + 2.0 * self.bfd * np.sqrt(1.0 + self.z*self.z)) + hr_trap = area_trap / wetted_trap + celerity_trap = read_celerity(self.bfd, hr_trap, self.n) + area_CC = (geom['stage'] - self.bfd) * self.TwCC + wetted_CC = (self.TwCC + (2.0 * (geom['stage'] - self.bfd))) + radius = (area_CC + area_trap) / (wetted_CC + wetted_trap) + celerity_cc = read_celerity(geom['stage'], radius, self.nCC) + celerity_CC = ((celerity_trap * area_trap) + (celerity_cc * area_CC)) / (area_trap + area_CC) + n_cc = ((self.n * wetted_trap) + (wetted_CC * self.nCC)) / (wetted_trap + wetted_CC) + + geom['area'][mask] = area_CC[mask] + geom['wetted_perimeter'][mask] = wetted_CC[mask] + geom['hydraulic_radius'][mask] = radius[mask] + geom['celerity'][mask] = celerity_CC[mask] + geom['mannings_n'][mask] = n_cc[mask] + geom['top_width'][mask] = self.TwCC + + geom['discharge'] = (1.0 / geom['mannings_n']) * geom['area'] * (geom['hydraulic_radius'] ** (2.0 / 3.0)) * (self.So ** 0.5) geom['log_q'] = np.log(geom['discharge']) - geom['log_width'] = np.log(geom['top_width']) From 2285d78bac9270bff910efe899d4a9ecd7251890 Mon Sep 17 00:00:00 2001 From: sclaw Date: Thu, 20 Jun 2024 17:55:58 -0400 Subject: [PATCH 12/43] make stage hydrographs --- src/muskingumcunge/network.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/muskingumcunge/network.py b/src/muskingumcunge/network.py index b7a3893..6a56dc5 100644 --- a/src/muskingumcunge/network.py +++ b/src/muskingumcunge/network.py @@ -21,6 +21,7 @@ def __init__(self, edge_dict): self.reach_dict = None self.forcing_df = None self.channel_outflows = {k: None for k in self.edge_dict.keys()} + self.channel_outflows_stage = {k: None for k in self.edge_dict.keys()} self.out_df = None def calculate_post_order(self): @@ -116,6 +117,7 @@ def run_event(self, optimize_dx=True, conserve_mass=False, lat_addition='middle' out_flow = routed.sum() routed = routed * (in_flow / out_flow) self.channel_outflows[node] = routed + self.channel_outflows_stage[node] = np.interp(routed, reach.geometry['discharge'], reach.geometry['stage']) # calculate total outflow us_root = self.chlid_dict[self.root] @@ -128,5 +130,8 @@ def run_event(self, optimize_dx=True, conserve_mass=False, lat_addition='middle' out_df = pd.DataFrame(self.channel_outflows) out_df['dt'] = self.forcing_df.index self.out_df = out_df.set_index('dt') + stage_df = pd.DataFrame(self.channel_outflows_stage) + stage_df['dt'] = self.forcing_df.index + self.stage_df = stage_df.set_index('dt') From 4abff740d426e189ff155a4ac37d6881c79c96d2 Mon Sep 17 00:00:00 2001 From: sclaw Date: Fri, 21 Jun 2024 09:04:20 -0400 Subject: [PATCH 13/43] catchup --- src/muskingumcunge/reach.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/muskingumcunge/reach.py b/src/muskingumcunge/reach.py index 25b11f6..4cfd254 100644 --- a/src/muskingumcunge/reach.py +++ b/src/muskingumcunge/reach.py @@ -65,7 +65,6 @@ def route_hydrograph(self, inflows, dt, lateral=None, max_iter=1000): lateral = np.zeros_like(inflows) lateral = (lateral[:-1] + lateral[1:]) / 2 lateral = np.append(lateral, 0) - for i in range(len(inflows) - 1): q_guess = sum([inflows[i], inflows[i + 1], outflows[i], lateral[i]]) / 3 last_guess = q_guess * 2 @@ -73,7 +72,8 @@ def route_hydrograph(self, inflows, dt, lateral=None, max_iter=1000): while abs(last_guess - q_guess) > 0.003: # from handbook of hydrology page 328 counter += 1 last_guess = q_guess.copy() - reach_q = sum([inflows[i], inflows[i + 1], outflows[i], q_guess]) / 4 + # reach_q = sum([inflows[i], inflows[i + 1], outflows[i], q_guess]) / 4 + reach_q = sum([inflows[i], inflows[i], outflows[i], q_guess]) / 4 # Interpolate log_reach_q = np.log(reach_q) @@ -91,10 +91,10 @@ def route_hydrograph(self, inflows, dt, lateral=None, max_iter=1000): q_guess = (c0 * inflows[i + 1]) + (c1 * inflows[i]) + (c2 * outflows[i]) + (c3 * lateral[i]) q_guess = max(min(inflows), q_guess) if counter == max_iter: - print('hit max iter') last_guess = q_guess outflows.append(q_guess) + return np.array(outflows) def route_hydrograph_mct(self, inflows, dt, max_iter=1000): @@ -412,6 +412,8 @@ def generate_geometry(self): geom['wetted_perimeter'][mask] = wetted_CC[mask] geom['hydraulic_radius'][mask] = radius[mask] geom['celerity'][mask] = celerity_CC[mask] + if np.any(geom['celerity'] < 0) or np.any(np.isnan(geom['celerity'])) or np.any(np.isinf(geom['celerity'])): + print('bad celerity') geom['mannings_n'][mask] = n_cc[mask] geom['top_width'][mask] = self.TwCC From 6a60bbc9a302cc3a4792db4522edf243e9080ad8 Mon Sep 17 00:00:00 2001 From: sclaw Date: Fri, 21 Jun 2024 09:04:33 -0400 Subject: [PATCH 14/43] controller for network runs --- src/muskingumcunge/controller.py | 186 +++++++++++++++++++++++++++++++ 1 file changed, 186 insertions(+) create mode 100644 src/muskingumcunge/controller.py diff --git a/src/muskingumcunge/controller.py b/src/muskingumcunge/controller.py new file mode 100644 index 0000000..8d2bb23 --- /dev/null +++ b/src/muskingumcunge/controller.py @@ -0,0 +1,186 @@ +import os +import sys +import json +import pandas as pd +import numpy as np +from muskingumcunge.reach import CompoundReach, CustomReach, MuskingumReach, WRFCompound +from muskingumcunge.network import Network + + +def create_project(path): + template = { + "base_dir": os.path.dirname(path), + "meta_path": os.path.join(os.path.dirname(path), "reach_meta.csv"), + "lateral_path": os.path.join(os.path.dirname(path), "laterals.csv"), + "head_path": os.path.join(os.path.dirname(path), "headwaters.csv"), + "lake_path": os.path.join(os.path.dirname(path), "lake.csv"), + "geometry_dir": "", + "geometry_source": "NWM", + "optimize_dx": True, + "conserve_mass": False, + "lat_addition": "middle", + "id_field": "comid", + "to_id_field": "ds_comid", + "mode": 'basin', + "out_dir": os.path.join(os.path.dirname(path), 'outputs'), + } + os.makedirs(template['out_dir'], exist_ok=True) + with open(path, 'w') as f: + json.dump(template, f, indent=4) + +def clean_ts(df): + if 'Unnamed: 0' in df.columns: + df = df.drop(columns='Unnamed: 0') + if 'datetime' in df.columns: + df["datetime"] = df["datetime"].map(lambda x: x + " 00:00:00" if len(x) == 10 else x) + df["datetime"] = pd.to_datetime(df["datetime"]) + df = df.set_index("datetime") + return df + +def load_geom(meta_df, source='NWM', geom_dir=None): + reaches = dict() + if source == 'HAND': + hand_tw = pd.read_csv(os.path.join(geom_dir, 'area.csv')) + hand_el = pd.read_csv(os.path.join(geom_dir, 'el.csv')) + hand_wp = pd.read_csv(os.path.join(geom_dir, 'p.csv')) + hand_ar = pd.read_csv(os.path.join(geom_dir, 'vol.csv')) + for r in list(meta_df.index): + ch_n = meta_df.loc[r, 'n'] + fp_n = meta_df.loc[r, 'nCC'] + slope = meta_df.loc[r, 'slope'] + length = meta_df.loc[r, 'length'] + + tw = meta_df.loc[r, 'TopWdth'] + bw = meta_df.loc[r, 'BtmWdth'] + z = meta_df.loc[r, 'ChSlp'] + tw_cc = meta_df.loc[r, 'TopWdth'] + bf = (tw - bw) / (2 * z) + + if source == 'NWM': + reaches[r] = WRFCompound(bw, z, bf, tw_cc, ch_n, fp_n, slope, length, max_stage=10*bf, stage_resolution=3000) + elif source == 'HAND': + try: + tw = hand_tw[r].to_numpy() / length + wp = hand_wp[r].to_numpy() / length + ar = hand_ar[r].to_numpy() / length + el = hand_el[r].to_numpy() + n = np.ones(el.shape) * fp_n + n[el < bf] = ch_n + n[el >= bf] = fp_n + reaches[r] = CustomReach(n, slope, length, el, tw, ar, wp) + except: + # Error catching. Default to NWM channel + print(f'Error loading HAND data for reach {r}. Defaulting to NWM channel') + reaches[r] = WRFCompound(bw, z, bf, tw_cc, ch_n, fp_n, slope, length, max_stage=10*bf, stage_resolution=3000) + elif source == 'MUSKINGUM': + x = meta_df.loc[r, 'MusX'] + k = meta_df.loc[r, 'MusK'] + reaches[r] = MuskingumReach(k, x) + + # reaches[r].geometry['celerity'] = reaches[r].geometry['celerity'] * 0.1 # debugging + return reaches + +def execute_by_reach(meta_path): + # Load data + with open(meta_path, 'r') as f: + paths = json.load(f) + + meta_df = pd.read_csv(paths['meta_path']) + meta_df['comid'] = meta_df['comid'].astype(int).astype(str) + meta_df['ds_comid'] = meta_df['ds_comid'].astype(int).astype(str) + meta_df = meta_df.set_index("comid") + + upstream = pd.read_csv(os.path.join(paths['base_dir'], 'upstream_hydrographs.csv')) + upstream = clean_ts(upstream) + upstream = upstream.reindex(sorted(upstream.columns), axis=1) + + # Set up reaches + reaches = load_geom(meta_df, source=paths['geometry_source'], geom_dir=paths['geometry_dir']) + + # Run + out_df = pd.DataFrame(index=upstream.index, columns=upstream.columns) + stage_df = pd.DataFrame(index=upstream.index, columns=upstream.columns) + dt = (upstream.index[1] - upstream.index[0]).seconds / 3600 + counter = 1 + for r in upstream.columns: + if counter % 100 == 0: + print(f"{counter} / {len(upstream.columns)}") + try: + reach = reaches[r] + outflow = reach.route_hydrograph(upstream[r].to_numpy(), dt) + except Exception as e: + print(f"Error routing reach {r}: ({type(e).__name__}) {e}") + outflow = np.zeros(len(upstream)) + out_df[r] = outflow + stage_df[r] = np.interp(outflow, reach.geometry['discharge'], reach.geometry['stage']) + + os.makedirs(paths['out_dir'], exist_ok=True) + out_df.to_csv(os.path.join(paths['out_dir'], "model_outflows.csv")) + stage_df.to_csv(os.path.join(paths['out_dir'], "model_outflows_stage.csv")) + +def execute(meta_path): + # Load data + with open(meta_path, 'r') as f: + paths = json.load(f) + + meta_df = pd.read_csv(paths['meta_path']) + meta_df['comid'] = meta_df['comid'].astype(int).astype(str) + meta_df['ds_comid'] = meta_df['ds_comid'].astype(int).astype(str) + meta_df = meta_df.set_index("comid") + edge_dict = meta_df['ds_comid'].to_dict() + + lateral_df = pd.read_csv(paths['lateral_path']) + lateral_df = clean_ts(lateral_df) + lateral_df = lateral_df.reindex(sorted(lateral_df.columns), axis=1) + + head_df = pd.read_csv(paths['head_path']) + head_df = clean_ts(head_df) + head_df = head_df.reindex(sorted(head_df.columns), axis=1) + + lake_df = pd.read_csv(paths['lake_path']) + lake_df = clean_ts(lake_df) + lake_df = lake_df.reindex(sorted(lake_df.columns), axis=1) + + # # resample to 5 minute increments + # lateral_df = lateral_df.resample('5T').ffill() + # head_df = head_df.resample('5T').ffill() + # lake_df = lake_df.resample('5T').ffill() + + # Error checking + missing_forcings = list(set(meta_df.index).difference(set(lateral_df.columns))) + print(f'Missing forcing data for: {missing_forcings}') + + network = Network(edge_dict) + + # Set up reaches + reaches = load_geom(meta_df, source=paths['geometry_source'], geom_dir=paths['geometry_dir']) + + # Set up run + network = Network(edge_dict) + network.load_forcings(lateral_df) + network.load_headwater_forcings(head_df) + network.load_lake_forcings(lake_df) + network.load_reaches(reaches) + + # Run model + network.run_event(optimize_dx=paths['optimize_dx'], conserve_mass=paths['conserve_mass'], lat_addition=paths['lat_addition']) + os.makedirs(paths['out_dir'], exist_ok=True) + out_path = os.path.join(paths['out_dir'], "model_outflows.csv") + network.out_df.to_csv(out_path) + out_path = os.path.join(paths['out_dir'], "model_outflows_stage.csv") + network.stage_df.to_csv(out_path) + + return network + +def handler(path): + with open(path, 'r') as f: + dict = json.load(f) + + if dict['mode'] == 'reach': + execute_by_reach(path) + elif dict['mode'] == 'basin': + execute(path) + +if __name__ == '__main__': + path = sys.argv[1] + handler(path) \ No newline at end of file From f7059e8f8d52f5d9d1ca403ee071d6909971de18 Mon Sep 17 00:00:00 2001 From: sclaw Date: Fri, 21 Jun 2024 09:04:51 -0400 Subject: [PATCH 15/43] add conda yml --- environment.yml | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) create mode 100644 environment.yml diff --git a/environment.yml b/environment.yml new file mode 100644 index 0000000..4e44082 --- /dev/null +++ b/environment.yml @@ -0,0 +1,16 @@ +name: floodplains_muskingum +channels: + - conda-forge +dependencies: + - python + - matplotlib + - ca-certificates + - openssl + - scipy + - certifi + - anaconda::cython + - pandas + - anaconda::pillow + - statsmodels + - anaconda::seaborn +prefix: /users/k/l/klawson1/.conda/envs/floodplains_muskingum From 164acf808c52247276fb455409fe6b1964d931f9 Mon Sep 17 00:00:00 2001 From: sclaw Date: Thu, 27 Jun 2024 12:12:21 -0400 Subject: [PATCH 16/43] export celerity and discharge curves --- src/muskingumcunge/controller.py | 39 +++++++++++++++++++++++++++----- 1 file changed, 33 insertions(+), 6 deletions(-) diff --git a/src/muskingumcunge/controller.py b/src/muskingumcunge/controller.py index 8d2bb23..83fa80a 100644 --- a/src/muskingumcunge/controller.py +++ b/src/muskingumcunge/controller.py @@ -39,6 +39,7 @@ def clean_ts(df): def load_geom(meta_df, source='NWM', geom_dir=None): reaches = dict() + geom_error_count = 0 if source == 'HAND': hand_tw = pd.read_csv(os.path.join(geom_dir, 'area.csv')) hand_el = pd.read_csv(os.path.join(geom_dir, 'el.csv')) @@ -57,10 +58,12 @@ def load_geom(meta_df, source='NWM', geom_dir=None): bf = (tw - bw) / (2 * z) if source == 'NWM': - reaches[r] = WRFCompound(bw, z, bf, tw_cc, ch_n, fp_n, slope, length, max_stage=10*bf, stage_resolution=3000) + reaches[r] = WRFCompound(bw, z, bf, tw_cc, ch_n, fp_n, slope, length, max_stage=10*bf, stage_resolution=999) elif source == 'HAND': try: tw = hand_tw[r].to_numpy() / length + if np.all(tw == 0): + raise RuntimeError('All top-width values are zero') wp = hand_wp[r].to_numpy() / length ar = hand_ar[r].to_numpy() / length el = hand_el[r].to_numpy() @@ -68,16 +71,18 @@ def load_geom(meta_df, source='NWM', geom_dir=None): n[el < bf] = ch_n n[el >= bf] = fp_n reaches[r] = CustomReach(n, slope, length, el, tw, ar, wp) - except: + except Exception as e: # Error catching. Default to NWM channel print(f'Error loading HAND data for reach {r}. Defaulting to NWM channel') - reaches[r] = WRFCompound(bw, z, bf, tw_cc, ch_n, fp_n, slope, length, max_stage=10*bf, stage_resolution=3000) + print(f'Error: {e}') + geom_error_count += 1 + reaches[r] = WRFCompound(bw, z, bf, tw_cc, ch_n, fp_n, slope, length, max_stage=10*bf, stage_resolution=len(hand_tw)) elif source == 'MUSKINGUM': x = meta_df.loc[r, 'MusX'] k = meta_df.loc[r, 'MusK'] reaches[r] = MuskingumReach(k, x) - - # reaches[r].geometry['celerity'] = reaches[r].geometry['celerity'] * 0.1 # debugging + + print(f'Error loading {geom_error_count} / {len(reaches)} reaches') return reaches def execute_by_reach(meta_path): @@ -97,6 +102,17 @@ def execute_by_reach(meta_path): # Set up reaches reaches = load_geom(meta_df, source=paths['geometry_source'], geom_dir=paths['geometry_dir']) + # export celerity curves + stages = {r: reaches[r].geometry['stage'] for r in reaches} + stage_df = pd.DataFrame().from_dict(stages, orient='columns') + stage_df.to_csv(os.path.join(paths['out_dir'], "stage.csv")) + celerities = {r: reaches[r].geometry['celerity'] for r in reaches} + celerity_df = pd.DataFrame().from_dict(celerities, orient='columns') + celerity_df.to_csv(os.path.join(paths['out_dir'], "celerity.csv")) + discharge = {r: reaches[r].geometry['discharge'] for r in reaches} + discharge_df = pd.DataFrame().from_dict(discharge, orient='columns') + discharge_df.to_csv(os.path.join(paths['out_dir'], "discharge.csv")) + # Run out_df = pd.DataFrame(index=upstream.index, columns=upstream.columns) stage_df = pd.DataFrame(index=upstream.index, columns=upstream.columns) @@ -154,7 +170,18 @@ def execute(meta_path): # Set up reaches reaches = load_geom(meta_df, source=paths['geometry_source'], geom_dir=paths['geometry_dir']) - + + # export celerity curves + stages = {r: reaches[r].geometry['stage'] for r in reaches} + stage_df = pd.DataFrame().from_dict(stages, orient='columns') + stage_df.to_csv(os.path.join(paths['out_dir'], "stage.csv")) + celerities = {r: reaches[r].geometry['celerity'] for r in reaches} + celerity_df = pd.DataFrame().from_dict(celerities, orient='columns') + celerity_df.to_csv(os.path.join(paths['out_dir'], "celerity.csv")) + discharge = {r: reaches[r].geometry['discharge'] for r in reaches} + discharge_df = pd.DataFrame().from_dict(discharge, orient='columns') + discharge_df.to_csv(os.path.join(paths['out_dir'], "discharge.csv")) + # Set up run network = Network(edge_dict) network.load_forcings(lateral_df) From 04c7a76fb74afde06d6020fe6d196f314099673f Mon Sep 17 00:00:00 2001 From: sclaw Date: Thu, 27 Jun 2024 12:12:51 -0400 Subject: [PATCH 17/43] fix issue where gemetry was pooled across multiplpe reaches (for compound and custom) --- src/muskingumcunge/reach.py | 27 +++++++++++++++++++++------ 1 file changed, 21 insertions(+), 6 deletions(-) diff --git a/src/muskingumcunge/reach.py b/src/muskingumcunge/reach.py index 4cfd254..dc79cb0 100644 --- a/src/muskingumcunge/reach.py +++ b/src/muskingumcunge/reach.py @@ -301,7 +301,13 @@ def __init__(self, mannings_n, slope, reach_length, stages, top_widths, areas, p self.generate_geometry(stages, top_widths, areas, perimeters) def generate_geometry(self, stages, top_widths, areas, perimeters): - geom = self.geometry + geom = {'stage': None, + 'top_width': None, + 'area': None, + 'wetted_perimeter': None, + 'hydraulic_radius': None, + 'mannings_n': None, + 'discharge': None} geom['stage'] = stages geom['top_width'] = top_widths geom['log_width'] = np.log(geom['top_width']) @@ -315,7 +321,7 @@ def generate_geometry(self, stages, top_widths, areas, perimeters): geom['mannings_n'] = np.repeat(self.mannings_n, len(stages)) geom['discharge'] = (1 / geom['mannings_n']) * geom['area'] * (geom['hydraulic_radius'] ** (2 / 3)) * (self.slope ** 0.5) geom['discharge'][geom['discharge'] <= 0] = geom['discharge'][geom['discharge'] > 0].min() - self.clean_looped_rating_curve() + geom['discharge'] = self.clean_looped_rating_curve(geom['discharge']) geom['log_q'] = np.log(geom['discharge']) dp = geom['wetted_perimeter'][1:] - geom['wetted_perimeter'][:-1] @@ -331,16 +337,18 @@ def generate_geometry(self, stages, top_widths, areas, perimeters): geom['celerity'][0] = geom['celerity'][1] geom['celerity'] = np.nan_to_num(geom['celerity']) ## TODO: geom['celerity'][<0] = 0 + self.geometry = geom - def clean_looped_rating_curve(self): + def clean_looped_rating_curve(self, discharge): # Enforce monotonic discharge inreases q_last = np.nan - for ind, q in enumerate(self.geometry['discharge']): + for ind, q in enumerate(discharge): if q < q_last: - self.geometry['discharge'][ind] = q_last + discharge[ind] = q_last else: q_last = q + return discharge class MuskingumReach: @@ -382,7 +390,13 @@ def __init__(self, bottom_width, side_slope, bankfull_depth, floodplain_width, c self.generate_geometry() def generate_geometry(self): - geom = self.geometry + geom = {'stage': None, + 'top_width': None, + 'area': None, + 'wetted_perimeter': None, + 'hydraulic_radius': None, + 'mannings_n': None, + 'discharge': None} geom['stage'] = np.linspace(0, self.max_stage, self.resolution) @@ -420,6 +434,7 @@ def generate_geometry(self): geom['discharge'] = (1.0 / geom['mannings_n']) * geom['area'] * (geom['hydraulic_radius'] ** (2.0 / 3.0)) * (self.So ** 0.5) geom['log_q'] = np.log(geom['discharge']) geom['log_width'] = np.log(geom['top_width']) + self.geometry = geom class SigmoidalReach(BaseReach): From f6fa51bf985b9596757a5f848bd929dd13fa7f48 Mon Sep 17 00:00:00 2001 From: sclaw Date: Fri, 12 Jul 2024 16:55:58 -0400 Subject: [PATCH 18/43] add option to use regression instead of hydrofabric two-stage params --- src/muskingumcunge/controller.py | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/src/muskingumcunge/controller.py b/src/muskingumcunge/controller.py index 83fa80a..4395722 100644 --- a/src/muskingumcunge/controller.py +++ b/src/muskingumcunge/controller.py @@ -14,7 +14,7 @@ def create_project(path): "lateral_path": os.path.join(os.path.dirname(path), "laterals.csv"), "head_path": os.path.join(os.path.dirname(path), "headwaters.csv"), "lake_path": os.path.join(os.path.dirname(path), "lake.csv"), - "geometry_dir": "", + "geometry_dir": "/users/k/l/klawson1/netfiles/ciroh/floodplainsData/runs/NWM/geometry", "geometry_source": "NWM", "optimize_dx": True, "conserve_mass": False, @@ -50,15 +50,24 @@ def load_geom(meta_df, source='NWM', geom_dir=None): fp_n = meta_df.loc[r, 'nCC'] slope = meta_df.loc[r, 'slope'] length = meta_df.loc[r, 'length'] + da = meta_df.loc[r, 'TotDASqKm'] tw = meta_df.loc[r, 'TopWdth'] bw = meta_df.loc[r, 'BtmWdth'] z = meta_df.loc[r, 'ChSlp'] - tw_cc = meta_df.loc[r, 'TopWdth'] + tw_cc = meta_df.loc[r, 'TopWdthCC'] bf = (tw - bw) / (2 * z) if source == 'NWM': - reaches[r] = WRFCompound(bw, z, bf, tw_cc, ch_n, fp_n, slope, length, max_stage=10*bf, stage_resolution=999) + reaches[r] = WRFCompound(bw, z, bf, tw_cc, ch_n, fp_n, slope, length, max_stage=20*bf, stage_resolution=999) + elif source == 'NWM_Regression': + tw = 2.44 * (da ** 0.34) + a_ch = 0.75 * (da ** 0.53) + bf = (a_ch / tw) * 1.25 + bw = ((2 * a_ch) / bf) - tw + z = (tw - bw) / bf + tw_cc = 3 * tw + reaches[r] = WRFCompound(bw, z, bf, tw_cc, ch_n, fp_n, slope, length, max_stage=20*bf, stage_resolution=999) elif source == 'HAND': try: tw = hand_tw[r].to_numpy() / length From e55d12e4787e127d295ad7b4a40dc3a16b7e0c69 Mon Sep 17 00:00:00 2001 From: sclaw Date: Fri, 12 Jul 2024 17:22:53 -0400 Subject: [PATCH 19/43] dirty working commit --- src/muskingumcunge/reach.py | 62 ++++++++++++++++++++++++++----------- 1 file changed, 44 insertions(+), 18 deletions(-) diff --git a/src/muskingumcunge/reach.py b/src/muskingumcunge/reach.py index dc79cb0..83a4b8b 100644 --- a/src/muskingumcunge/reach.py +++ b/src/muskingumcunge/reach.py @@ -57,8 +57,10 @@ def route_hydrograph_c(self, inflows, dt): return croute(inflows, dt, reach_length, slope, geometry) def route_hydrograph(self, inflows, dt, lateral=None, max_iter=1000): + dt = dt * 60 * 60 # MUST DELETE AFTER DEBUGGING outflows = list() - outflows.append((inflows[0])) + # outflows.append((inflows[0])) + outflows.append(0.9 * inflows[0]) assert max(inflows) < max(self.geometry['discharge']), 'Rating Curve does not cover range of flowrates in hydrograph' if lateral is None: @@ -66,7 +68,8 @@ def route_hydrograph(self, inflows, dt, lateral=None, max_iter=1000): lateral = (lateral[:-1] + lateral[1:]) / 2 lateral = np.append(lateral, 0) for i in range(len(inflows) - 1): - q_guess = sum([inflows[i], inflows[i + 1], outflows[i], lateral[i]]) / 3 + # q_guess = sum([inflows[i], inflows[i + 1], outflows[i], lateral[i]]) / 3 + q_guess = sum([inflows[i], inflows[i], outflows[i], lateral[i]]) / 3 last_guess = q_guess * 2 counter = 1 while abs(last_guess - q_guess) > 0.003: # from handbook of hydrology page 328 @@ -80,16 +83,34 @@ def route_hydrograph(self, inflows, dt, lateral=None, max_iter=1000): b_tmp = np.exp(np.interp(log_reach_q, self.geometry['log_q'], self.geometry['log_width'])) c_tmp = np.interp(log_reach_q, self.geometry['log_q'], self.geometry['celerity']) - courant = c_tmp * dt * 60 * 60 / self.reach_length - reynold = reach_q / (self.slope * c_tmp * self.reach_length * b_tmp) + # courant = c_tmp * dt * 60 * 60 / self.reach_length + # reynold = reach_q / (self.slope * c_tmp * self.reach_length * b_tmp) + + # c0 = (-1 + courant + reynold) / (1 + courant + reynold) + # c1 = (1 + courant - reynold) / (1 + courant + reynold) + # c2 = (1 - courant + reynold) / (1 + courant + reynold) + # c3 = (2 * courant) / (1 + courant + reynold) + + # q_guess = (c0 * inflows[i + 1]) + (c1 * inflows[i]) + (c2 * outflows[i]) + (c3 * lateral[i]) + # q_guess = max(min(inflows), q_guess) + + # WRF way + if c_tmp > 0: + Km = self.reach_length / c_tmp + Km = max(dt, Km) + else: + Km = dt + X = 0.5*(1-(reach_q/(2.0*b_tmp*self.slope*c_tmp*self.reach_length))) + X = min(0.5, max(0.0, X)) + D = (Km*(1.000 - X) + (dt)/2.0000) + + C1 = (Km*X + dt/2.000000)/D + C2 = (dt/2.0000 - Km*X)/D + C3 = (Km*(1.00000000-X)-dt/2.000000)/D + C4 = (lateral[i]*dt)/D + # q_guess = ((C1*inflows[i]) + (C2*inflows[i + 1]) + (C3*outflows[i])) + q_guess = ((C1*inflows[i]) + (C2*inflows[i]) + (C3*outflows[i])) - c0 = (-1 + courant + reynold) / (1 + courant + reynold) - c1 = (1 + courant - reynold) / (1 + courant + reynold) - c2 = (1 - courant + reynold) / (1 + courant + reynold) - c3 = (2 * courant) / (1 + courant + reynold) - - q_guess = (c0 * inflows[i + 1]) + (c1 * inflows[i]) + (c2 * outflows[i]) + (c3 * lateral[i]) - q_guess = max(min(inflows), q_guess) if counter == max_iter: last_guess = q_guess outflows.append(q_guess) @@ -376,7 +397,7 @@ class WRFCompound(BaseReach): def __init__(self, bottom_width, side_slope, bankfull_depth, floodplain_width, channel_n, floodplain_n, slope, reach_length, max_stage=10, stage_resolution=50): self.Bw = bottom_width - self.z = side_slope # Should be rise/run per side + self.z = side_slope / 2 # Should be rise/run per side self.bfd = bankfull_depth self.TwCC = floodplain_width self.n = channel_n @@ -401,13 +422,14 @@ def generate_geometry(self): geom['stage'] = np.linspace(0, self.max_stage, self.resolution) # Trapezoidal Channel + geom['top_width'] = self.Bw + (self.z * geom['stage'] * 2) geom['area'] = (self.Bw + geom['stage'] * self.z) * geom['stage'] geom['wetted_perimeter'] = (self.Bw + 2.0 * geom['stage'] * np.sqrt(1.0 + self.z*self.z)) geom['hydraulic_radius'] = geom['area'] / geom['wetted_perimeter'] read_celerity = lambda y, r, n: (np.sqrt(self.So)/n)*((5./3.)*r**(2./3.) - ((2./3.)*r**(5./3.)*(2.0*np.sqrt(1.0 + self.z*self.z)/(self.Bw+2.0*y*self.z)))) geom['celerity'] = read_celerity(geom['stage'], geom['hydraulic_radius'], self.n) geom['mannings_n'] = np.repeat(self.n, self.resolution) - geom['top_width'] = self.Bw + (2 * self.z * geom['stage']) + geom['discharge'] = (1.0 / geom['mannings_n']) * geom['area'] * (geom['hydraulic_radius'] ** (2.0 / 3.0)) * (self.So ** 0.5) # Compound Channel mask = (geom['stage'] > self.bfd) @@ -418,20 +440,24 @@ def generate_geometry(self): area_CC = (geom['stage'] - self.bfd) * self.TwCC wetted_CC = (self.TwCC + (2.0 * (geom['stage'] - self.bfd))) radius = (area_CC + area_trap) / (wetted_CC + wetted_trap) - celerity_cc = read_celerity(geom['stage'], radius, self.nCC) - celerity_CC = ((celerity_trap * area_trap) + (celerity_cc * area_CC)) / (area_trap + area_CC) + # celerity_cc = read_celerity(geom['stage'], radius, self.nCC) + celerity_cc = (np.sqrt(self.So)/self.nCC)*((5./3.)*(geom['stage'] - self.bfd)**(2./3.)) + celerity_cc = ((celerity_trap * area_trap) + (celerity_cc * area_CC)) / (area_trap + area_CC) n_cc = ((self.n * wetted_trap) + (wetted_CC * self.nCC)) / (wetted_trap + wetted_CC) + discharge_cc = (1.0 / self.nCC) * area_CC * ((geom['stage'] - self.bfd) ** (2.0 / 3.0)) * (self.So ** 0.5) + discharge_trap = (1.0 / self.n) * area_trap * ((area_trap / wetted_trap) ** (2.0 / 3.0)) * (self.So ** 0.5) - geom['area'][mask] = area_CC[mask] + geom['area'][mask] = area_CC[mask] + area_trap geom['wetted_perimeter'][mask] = wetted_CC[mask] geom['hydraulic_radius'][mask] = radius[mask] - geom['celerity'][mask] = celerity_CC[mask] + geom['celerity'][mask] = celerity_cc[mask] if np.any(geom['celerity'] < 0) or np.any(np.isnan(geom['celerity'])) or np.any(np.isinf(geom['celerity'])): print('bad celerity') geom['mannings_n'][mask] = n_cc[mask] geom['top_width'][mask] = self.TwCC + geom['discharge'][mask] = discharge_cc[mask] + discharge_trap - geom['discharge'] = (1.0 / geom['mannings_n']) * geom['area'] * (geom['hydraulic_radius'] ** (2.0 / 3.0)) * (self.So ** 0.5) + # geom['discharge'] = (1.0 / geom['mannings_n']) * geom['area'] * (geom['hydraulic_radius'] ** (2.0 / 3.0)) * (self.So ** 0.5) geom['log_q'] = np.log(geom['discharge']) geom['log_width'] = np.log(geom['top_width']) self.geometry = geom From 1e4d2bbb694f377ad34d425fc8459b2db397235a Mon Sep 17 00:00:00 2001 From: sclaw Date: Tue, 16 Jul 2024 13:32:12 -0400 Subject: [PATCH 20/43] catchup --- src/muskingumcunge/reach.py | 25 ++++++++++++++++++------- 1 file changed, 18 insertions(+), 7 deletions(-) diff --git a/src/muskingumcunge/reach.py b/src/muskingumcunge/reach.py index 83a4b8b..3a66ae7 100644 --- a/src/muskingumcunge/reach.py +++ b/src/muskingumcunge/reach.py @@ -56,11 +56,13 @@ def route_hydrograph_c(self, inflows, dt): return croute(inflows, dt, reach_length, slope, geometry) - def route_hydrograph(self, inflows, dt, lateral=None, max_iter=1000): + def route_hydrograph(self, inflows, dt, lateral=None, max_iter=1000, initial_outflow=None): dt = dt * 60 * 60 # MUST DELETE AFTER DEBUGGING outflows = list() - # outflows.append((inflows[0])) - outflows.append(0.9 * inflows[0]) + if initial_outflow is not None: + outflows.append(initial_outflow) + else: + outflows.append(inflows[0]) assert max(inflows) < max(self.geometry['discharge']), 'Rating Curve does not cover range of flowrates in hydrograph' if lateral is None: @@ -118,6 +120,12 @@ def route_hydrograph(self, inflows, dt, lateral=None, max_iter=1000): return np.array(outflows) + def route_hydrograph_wrf_hydro(self, inflows, dt, lateral=None, max_iter=1000, initial_outflow=None): + pass + + + + def route_hydrograph_mct(self, inflows, dt, max_iter=1000): outflows = list() outflows.append(inflows[0]) @@ -436,16 +444,18 @@ def generate_geometry(self): area_trap = (self.Bw + self.bfd * self.z) * self.bfd wetted_trap = (self.Bw + 2.0 * self.bfd * np.sqrt(1.0 + self.z*self.z)) hr_trap = area_trap / wetted_trap - celerity_trap = read_celerity(self.bfd, hr_trap, self.n) + # celerity_trap = read_celerity(self.bfd, hr_trap, self.n) area_CC = (geom['stage'] - self.bfd) * self.TwCC wetted_CC = (self.TwCC + (2.0 * (geom['stage'] - self.bfd))) radius = (area_CC + area_trap) / (wetted_CC + wetted_trap) + celerity_trap = read_celerity(self.bfd, radius, self.n) # Is this an error? I feel like radius should be capped at the bankfull radius? # celerity_cc = read_celerity(geom['stage'], radius, self.nCC) celerity_cc = (np.sqrt(self.So)/self.nCC)*((5./3.)*(geom['stage'] - self.bfd)**(2./3.)) celerity_cc = ((celerity_trap * area_trap) + (celerity_cc * area_CC)) / (area_trap + area_CC) n_cc = ((self.n * wetted_trap) + (wetted_CC * self.nCC)) / (wetted_trap + wetted_CC) - discharge_cc = (1.0 / self.nCC) * area_CC * ((geom['stage'] - self.bfd) ** (2.0 / 3.0)) * (self.So ** 0.5) - discharge_trap = (1.0 / self.n) * area_trap * ((area_trap / wetted_trap) ** (2.0 / 3.0)) * (self.So ** 0.5) + # discharge_cc = (1.0 / self.nCC) * area_CC * ((geom['stage'] - self.bfd) ** (2.0 / 3.0)) * (self.So ** 0.5) + # discharge_trap = (1.0 / self.n) * area_trap * ((area_trap / wetted_trap) ** (2.0 / 3.0)) * (self.So ** 0.5) + discharge_cc = (1.0 / self.nCC) * (area_CC + area_trap) * (radius ** (2.0 / 3.0)) * (self.So ** 0.5) geom['area'][mask] = area_CC[mask] + area_trap geom['wetted_perimeter'][mask] = wetted_CC[mask] @@ -455,7 +465,8 @@ def generate_geometry(self): print('bad celerity') geom['mannings_n'][mask] = n_cc[mask] geom['top_width'][mask] = self.TwCC - geom['discharge'][mask] = discharge_cc[mask] + discharge_trap + # geom['discharge'][mask] = discharge_cc[mask] + discharge_trap + geom['discharge'][mask] = discharge_cc[mask] # geom['discharge'] = (1.0 / geom['mannings_n']) * geom['area'] * (geom['hydraulic_radius'] ** (2.0 / 3.0)) * (self.So ** 0.5) geom['log_q'] = np.log(geom['discharge']) From 0d9fbb32afb0fed48db1fe1c46954b873a8f13e6 Mon Sep 17 00:00:00 2001 From: sclaw Date: Tue, 16 Jul 2024 19:23:12 -0400 Subject: [PATCH 21/43] fix lake handling. fix NWM geometry formats. Add initial outflow boundary condition --- src/muskingumcunge/controller.py | 38 ++++++++++++++++++++++++++------ 1 file changed, 31 insertions(+), 7 deletions(-) diff --git a/src/muskingumcunge/controller.py b/src/muskingumcunge/controller.py index 4395722..89fe47b 100644 --- a/src/muskingumcunge/controller.py +++ b/src/muskingumcunge/controller.py @@ -54,9 +54,9 @@ def load_geom(meta_df, source='NWM', geom_dir=None): tw = meta_df.loc[r, 'TopWdth'] bw = meta_df.loc[r, 'BtmWdth'] - z = meta_df.loc[r, 'ChSlp'] + z = 2 / (meta_df.loc[r, 'ChSlp']) tw_cc = meta_df.loc[r, 'TopWdthCC'] - bf = (tw - bw) / (2 * z) + bf = (tw - bw) / (z) if source == 'NWM': reaches[r] = WRFCompound(bw, z, bf, tw_cc, ch_n, fp_n, slope, length, max_stage=20*bf, stage_resolution=999) @@ -166,10 +166,27 @@ def execute(meta_path): lake_df = clean_ts(lake_df) lake_df = lake_df.reindex(sorted(lake_df.columns), axis=1) - # # resample to 5 minute increments - # lateral_df = lateral_df.resample('5T').ffill() - # head_df = head_df.resample('5T').ffill() - # lake_df = lake_df.resample('5T').ffill() + downstream = pd.read_csv(os.path.join(paths['base_dir'], 'downstream_hydrographs.csv')) + downstream = clean_ts(downstream) + downstream = downstream.reindex(sorted(downstream.columns), axis=1) + + # pre process lakes + tmp_cols = list(lake_df.columns) + lakes = [float(c) for c in tmp_cols] + for l in lakes: + # Find root of each lake + lake_reaches = list(meta_df[meta_df['NHDWaterbodyComID'] == l].index) + lake_root = lake_reaches[0] + outflow = meta_df.loc[lake_root, 'ds_comid'] + while outflow in lake_reaches: + lake_root = meta_df.loc[lake_root, 'ds_comid'] + outflow = meta_df.loc[lake_root, 'ds_comid'] + head_df[lake_root] = downstream[lake_root] + + # resample to 5 minute increments + lateral_df = lateral_df.resample('5T').ffill() + head_df = head_df.resample('5T').ffill() + lake_df = lake_df.resample('5T').ffill() # Error checking missing_forcings = list(set(meta_df.index).difference(set(lateral_df.columns))) @@ -195,9 +212,16 @@ def execute(meta_path): network = Network(edge_dict) network.load_forcings(lateral_df) network.load_headwater_forcings(head_df) - network.load_lake_forcings(lake_df) network.load_reaches(reaches) + # Set up initial outflows + for n in network.post_order: + obs_out = downstream[n].to_numpy()[0] + if np.isnan(obs_out): + obs_out = None + network.init_outflows[n] = obs_out + network.load_initial_conditions(network.init_outflows) # kind of unnecessary + # Run model network.run_event(optimize_dx=paths['optimize_dx'], conserve_mass=paths['conserve_mass'], lat_addition=paths['lat_addition']) os.makedirs(paths['out_dir'], exist_ok=True) From 58c364d164be4a5412f4c120ffa84a4791a836f5 Mon Sep 17 00:00:00 2001 From: sclaw Date: Tue, 16 Jul 2024 19:24:16 -0400 Subject: [PATCH 22/43] add ability to have initial outflow boundary condition --- src/muskingumcunge/network.py | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/src/muskingumcunge/network.py b/src/muskingumcunge/network.py index 6a56dc5..96a3805 100644 --- a/src/muskingumcunge/network.py +++ b/src/muskingumcunge/network.py @@ -22,6 +22,7 @@ def __init__(self, edge_dict): self.forcing_df = None self.channel_outflows = {k: None for k in self.edge_dict.keys()} self.channel_outflows_stage = {k: None for k in self.edge_dict.keys()} + self.init_outflows = {k: None for k in self.edge_dict.keys()} self.out_df = None def calculate_post_order(self): @@ -34,6 +35,7 @@ def calculate_post_order(self): children = self.chlid_dict[cur] if len(children) == 0: self.post_order.append(cur) + self.headwaters.append(cur) elif all([c in self.post_order for c in children]): self.post_order.append(cur) else: @@ -65,11 +67,9 @@ def load_headwater_forcings(self, headwater_forcings): for n in missing: self.channel_outflows[n] = np.zeros(len(self.forcing_df)) - def load_lake_forcings(self, lake_forcings): - lake_forcings = lake_forcings.fillna(0) - for n in lake_forcings.columns: - self.channel_outflows[n] = lake_forcings[n].to_numpy() - self.headwaters.append(n) + def load_initial_conditions(self, initial_conditions): + # loads initial outflows for each reach + self.init_outflows = initial_conditions def run_event(self, optimize_dx=True, conserve_mass=False, lat_addition='middle'): # calculate total inflow @@ -104,7 +104,11 @@ def run_event(self, optimize_dx=True, conserve_mass=False, lat_addition='middle' l = laterals else: l = None - routed = reach.route_hydrograph(routed, dt, lateral=l) + if i == 0 and node in self.init_outflows: + init_out = self.init_outflows[node] + else: + init_out = None + routed = reach.route_hydrograph(routed, dt, lateral=l, initial_outflow=init_out) except AssertionError as e: print(f"Error routing {node}: {e}") From 65564ad6e502778e436676751f998d4c5b3027e6 Mon Sep 17 00:00:00 2001 From: sclaw Date: Wed, 17 Jul 2024 09:06:47 -0400 Subject: [PATCH 23/43] stable matches NWM3.0 --- src/muskingumcunge/reach.py | 343 +++++++++++++++++++++++++++++++----- 1 file changed, 298 insertions(+), 45 deletions(-) diff --git a/src/muskingumcunge/reach.py b/src/muskingumcunge/reach.py index 3a66ae7..50600bb 100644 --- a/src/muskingumcunge/reach.py +++ b/src/muskingumcunge/reach.py @@ -119,12 +119,6 @@ def route_hydrograph(self, inflows, dt, lateral=None, max_iter=1000, initial_out return np.array(outflows) - - def route_hydrograph_wrf_hydro(self, inflows, dt, lateral=None, max_iter=1000, initial_outflow=None): - pass - - - def route_hydrograph_mct(self, inflows, dt, max_iter=1000): outflows = list() @@ -417,8 +411,40 @@ def __init__(self, bottom_width, side_slope, bankfull_depth, floodplain_width, c self.resolution = stage_resolution self.generate_geometry() - + def generate_geometry(self): + geom = dict() + geom = {'stage': None, + 'top_width': None, + 'area': None, + 'wetted_perimeter': None, + 'hydraulic_radius': None, + 'mannings_n': None, + 'discharge': None} + + geom['stage'] = np.linspace(0, self.max_stage, self.resolution) + for param in ['top_width', 'area', 'wetted_perimeter', 'hydraulic_radius', 'mannings_n', 'discharge', 'celerity']: + geom[param] = np.zeros(geom['stage'].shape) + + for ind, h_0 in enumerate(geom['stage']): + Ck, WP, WPC, n, nCC, AREA, AREAC, R, So = self.calc_params(h_0) + if h_0 > self.bfd: + Twl = self.TwCC + else: + Twl = self.Bw + 2.0*self.z*h_0 + q = ((1/(((WP*n)+(WPC*nCC))/(WP+WPC))) * (AREA+AREAC) * (R**(2./3.)) * np.sqrt(So)) + geom['top_width'][ind] = Twl + geom['area'][ind] = AREA + AREAC + geom['wetted_perimeter'][ind] = WP + WPC + geom['hydraulic_radius'][ind] = R + geom['mannings_n'][ind] = ((n * WP) + (nCC * WPC)) / (WP + WPC) + geom['celerity'][ind] = Ck + geom['discharge'][ind] = q + + self.geometry = geom + + + def generate_geometry_old(self): geom = {'stage': None, 'top_width': None, 'area': None, @@ -429,49 +455,276 @@ def generate_geometry(self): geom['stage'] = np.linspace(0, self.max_stage, self.resolution) + # rename variables to match WRF code + Bw = self.Bw + Tw = self.Bw + (2 * self.bfd * self.z) + TwCC = self.TwCC + nCC = self.nCC + So = self.So + n = self.n + z = self.z + + AREA, AREAC = None, None + Z = None # trapezoid distance + R, RC = None, None + WP, WPC = None, None + h_0 = geom['stage'] + bfd = self.bfd + + Ck = None + Twl = Bw + 2.0*z*h_0 # Top width at simulated flow + + compound_mask = (geom['stage'] > self.bfd) + # Trapezoidal Channel - geom['top_width'] = self.Bw + (self.z * geom['stage'] * 2) - geom['area'] = (self.Bw + geom['stage'] * self.z) * geom['stage'] - geom['wetted_perimeter'] = (self.Bw + 2.0 * geom['stage'] * np.sqrt(1.0 + self.z*self.z)) - geom['hydraulic_radius'] = geom['area'] / geom['wetted_perimeter'] - read_celerity = lambda y, r, n: (np.sqrt(self.So)/n)*((5./3.)*r**(2./3.) - ((2./3.)*r**(5./3.)*(2.0*np.sqrt(1.0 + self.z*self.z)/(self.Bw+2.0*y*self.z)))) - geom['celerity'] = read_celerity(geom['stage'], geom['hydraulic_radius'], self.n) - geom['mannings_n'] = np.repeat(self.n, self.resolution) - geom['discharge'] = (1.0 / geom['mannings_n']) * geom['area'] * (geom['hydraulic_radius'] ** (2.0 / 3.0)) * (self.So ** 0.5) + AREA = (Bw + h_0 * z ) * h_0 + WP = (Bw + 2.0 * h_0 * np.sqrt(1.0 + z*z)) + R = AREA / WP # Compound Channel - mask = (geom['stage'] > self.bfd) - area_trap = (self.Bw + self.bfd * self.z) * self.bfd - wetted_trap = (self.Bw + 2.0 * self.bfd * np.sqrt(1.0 + self.z*self.z)) - hr_trap = area_trap / wetted_trap - # celerity_trap = read_celerity(self.bfd, hr_trap, self.n) - area_CC = (geom['stage'] - self.bfd) * self.TwCC - wetted_CC = (self.TwCC + (2.0 * (geom['stage'] - self.bfd))) - radius = (area_CC + area_trap) / (wetted_CC + wetted_trap) - celerity_trap = read_celerity(self.bfd, radius, self.n) # Is this an error? I feel like radius should be capped at the bankfull radius? - # celerity_cc = read_celerity(geom['stage'], radius, self.nCC) - celerity_cc = (np.sqrt(self.So)/self.nCC)*((5./3.)*(geom['stage'] - self.bfd)**(2./3.)) - celerity_cc = ((celerity_trap * area_trap) + (celerity_cc * area_CC)) / (area_trap + area_CC) - n_cc = ((self.n * wetted_trap) + (wetted_CC * self.nCC)) / (wetted_trap + wetted_CC) - # discharge_cc = (1.0 / self.nCC) * area_CC * ((geom['stage'] - self.bfd) ** (2.0 / 3.0)) * (self.So ** 0.5) - # discharge_trap = (1.0 / self.n) * area_trap * ((area_trap / wetted_trap) ** (2.0 / 3.0)) * (self.So ** 0.5) - discharge_cc = (1.0 / self.nCC) * (area_CC + area_trap) * (radius ** (2.0 / 3.0)) * (self.So ** 0.5) - - geom['area'][mask] = area_CC[mask] + area_trap - geom['wetted_perimeter'][mask] = wetted_CC[mask] - geom['hydraulic_radius'][mask] = radius[mask] - geom['celerity'][mask] = celerity_cc[mask] - if np.any(geom['celerity'] < 0) or np.any(np.isnan(geom['celerity'])) or np.any(np.isinf(geom['celerity'])): - print('bad celerity') - geom['mannings_n'][mask] = n_cc[mask] - geom['top_width'][mask] = self.TwCC - # geom['discharge'][mask] = discharge_cc[mask] + discharge_trap - geom['discharge'][mask] = discharge_cc[mask] + AREA[compound_mask] = ((Bw + bfd * z) * bfd) + AREAC = np.zeros(AREA.shape) + AREAC[compound_mask] = (TwCC * (h_0 -bfd))[compound_mask] + WP[compound_mask] = (Bw + 2.0 * bfd * np.sqrt(1.0 + z*z)) + WPC = np.zeros(WP.shape) + WPC[compound_mask] = (TwCC + (2.0 * (h_0-bfd)))[compound_mask] + R[compound_mask] = ((AREA + AREAC)/(WP +WPC))[compound_mask] + R[np.isinf(R)] = 0.0 + + # Trapazoidal Channel + Ck = (np.sqrt(So)/n)*((5./3.)*R**(2./3.) - ((2./3.)*R**(5./3.)*(2.0*np.sqrt(1.0 + z*z)/(Bw+2.0*h_0*z)))) + Ck[h_0 == 0] = 0.0 + + # Compound Channel + Ck[compound_mask] = (((np.sqrt(So)/n)*((5./3.)*R**(2./3.) - ((2./3.)*R**(5./3.)*(2.0*np.sqrt(1.0 + z*z)/(Bw+2.0*bfd*z))))*AREA + ((np.sqrt(So)/(nCC))*(5./3.)*(h_0-bfd)**(2./3.))*AREAC)/(AREA+AREAC))[compound_mask] + Ck[Ck < 0] = 0.0 + + discharge = ((1/(((WP*n)+(WPC*nCC))/(WP+WPC))) * (AREA+AREAC) * (R**(2./3.)) * np.sqrt(So)) + + # log to class variables + self.geometry['stage'] = geom['stage'] + self.geometry['top_width'] = Twl + self.geometry['celerity'] = Ck + self.geometry['discharge'] = discharge + self.geometry['log_q'] = np.log(self.geometry['discharge']) + self.geometry['log_width'] = np.log(self.geometry['top_width']) + + + # geom['top_width'] = self.Bw + (self.z * geom['stage'] * 2) + # geom['area'] = (self.Bw + geom['stage'] * self.z) * geom['stage'] + # geom['wetted_perimeter'] = (self.Bw + 2.0 * geom['stage'] * np.sqrt(1.0 + self.z*self.z)) + # geom['hydraulic_radius'] = geom['area'] / geom['wetted_perimeter'] + # read_celerity = lambda y, r, n: (np.sqrt(self.So)/n)*((5./3.)*r**(2./3.) - ((2./3.)*r**(5./3.)*(2.0*np.sqrt(1.0 + self.z*self.z)/(self.Bw+2.0*y*self.z)))) + # geom['celerity'] = read_celerity(geom['stage'], geom['hydraulic_radius'], self.n) + # geom['mannings_n'] = np.repeat(self.n, self.resolution) # geom['discharge'] = (1.0 / geom['mannings_n']) * geom['area'] * (geom['hydraulic_radius'] ** (2.0 / 3.0)) * (self.So ** 0.5) - geom['log_q'] = np.log(geom['discharge']) - geom['log_width'] = np.log(geom['top_width']) - self.geometry = geom + + # # Compound Channel + # mask = (geom['stage'] > self.bfd) + # area_trap = (self.Bw + self.bfd * self.z) * self.bfd + # wetted_trap = (self.Bw + 2.0 * self.bfd * np.sqrt(1.0 + self.z*self.z)) + # hr_trap = area_trap / wetted_trap + # # celerity_trap = read_celerity(self.bfd, hr_trap, self.n) + # area_CC = (geom['stage'] - self.bfd) * self.TwCC + # wetted_CC = (self.TwCC + (2.0 * (geom['stage'] - self.bfd))) + # radius = (area_CC + area_trap) / (wetted_CC + wetted_trap) + # celerity_trap = read_celerity(self.bfd, radius, self.n) # Is this an error? I feel like radius should be capped at the bankfull radius? + # # celerity_cc = read_celerity(geom['stage'], radius, self.nCC) + # celerity_cc = (np.sqrt(self.So)/self.nCC)*((5./3.)*(geom['stage'] - self.bfd)**(2./3.)) + # celerity_cc = ((celerity_trap * area_trap) + (celerity_cc * area_CC)) / (area_trap + area_CC) + # n_cc = ((self.n * wetted_trap) + (wetted_CC * self.nCC)) / (wetted_trap + wetted_CC) + # # discharge_cc = (1.0 / self.nCC) * area_CC * ((geom['stage'] - self.bfd) ** (2.0 / 3.0)) * (self.So ** 0.5) + # # discharge_trap = (1.0 / self.n) * area_trap * ((area_trap / wetted_trap) ** (2.0 / 3.0)) * (self.So ** 0.5) + # discharge_cc = (1.0 / self.nCC) * (area_CC + area_trap) * (radius ** (2.0 / 3.0)) * (self.So ** 0.5) + + # geom['area'][mask] = area_CC[mask] + area_trap + # geom['wetted_perimeter'][mask] = wetted_CC[mask] + # geom['hydraulic_radius'][mask] = radius[mask] + # geom['celerity'][mask] = celerity_cc[mask] + # if np.any(geom['celerity'] < 0) or np.any(np.isnan(geom['celerity'])) or np.any(np.isinf(geom['celerity'])): + # print('bad celerity') + # geom['mannings_n'][mask] = n_cc[mask] + # geom['top_width'][mask] = self.TwCC + # # geom['discharge'][mask] = discharge_cc[mask] + discharge_trap + # geom['discharge'][mask] = discharge_cc[mask] + + # geom['discharge'] = (1.0 / geom['mannings_n']) * geom['area'] * (geom['hydraulic_radius'] ** (2.0 / 3.0)) * (self.So ** 0.5) + # geom['log_q'] = np.log(geom['discharge']) + # geom['log_width'] = np.log(geom['top_width']) + # self.geometry = geom + + def calc_params(self, h_0): + Bw = self.Bw + bfd = self.bfd + z = self.z + TwCC = self.TwCC + So = self.So + n = self.n + nCC = self.nCC + + WPC = 0 + AREAC = 0 + + if h_0 > bfd: + AREA = (Bw + bfd * z) * bfd + AREAC = TwCC * (h_0 - bfd) + WP = (Bw + 2.0 * bfd * np.sqrt(1.0 + z*z)) + WPC = TwCC + (2.0 * (h_0 - bfd)) + R = (AREA + AREAC) / (WP + WPC) + else: + AREA = (Bw + h_0 * z) * h_0 + WP = (Bw + 2.0 * h_0 * np.sqrt(1.0 + z*z)) + if WP > 0: + R = AREA / WP + else: + R = 0.0 + + if h_0 > bfd: + Ck = ((np.sqrt(So)/n)*((5./3.)*R**(2./3.) - ((2./3.)*R**(5./3.)*(2.0*np.sqrt(1.0 + z*z)/(Bw+2.0*bfd*z))))*AREA + ((np.sqrt(So)/(nCC))*(5./3.)*(h_0-bfd)**(2./3.))*AREAC)/(AREA+AREAC) + Ck = max(0.0, Ck) + else: + if h_0 > 0.0: + Ck = max(0.0,(np.sqrt(So)/n)*((5./3.)*R**(2./3.) - ((2./3.)*R**(5./3.)*(2.0*np.sqrt(1.0 + z*z)/(Bw+2.0*h_0*z))))) + else: + Ck = 0.0 + return Ck, WP, WPC, n, nCC, AREA, AREAC, R, So + + def route_hydrograph(self, inflows, dt, lateral=None, max_iter=1000, initial_outflow=None): + dt = dt * 60 * 60 # MUST DELETE AFTER DEBUGGING + dx = self.reach_length + outflows = list() + if initial_outflow is not None: + outflows.append(initial_outflow) + else: + outflows.append(inflows[0]) + assert max(inflows) < max(self.geometry['discharge']), 'Rating Curve does not cover range of flowrates in hydrograph' + + if lateral is None: + lateral = np.zeros_like(inflows) + lateral = (lateral[:-1] + lateral[1:]) / 2 + lateral = np.append(lateral, 0) + + # local variables + mindepth = 0.01 + max_iter = 100 + short_ts = True + + # initialize secant method + depth = np.interp(outflows[0], self.geometry['discharge'], self.geometry['stage']) + h = (depth * 1.33) + mindepth + h_0 = (depth * 0.67) + + for i in range(len(inflows) - 1): + if not ((inflows[i] > 0.0) or (inflows[i + 1] > 0.0) or (outflows[i] > 0.0)): + outflows.append(0.0) + continue + Qj_0 = 0.0 + iter = 0 + rerror = 1.0 + aerror = 0.01 + tries = 0 + + while rerror > 0.01 and aerror >= mindepth and iter <= max_iter: + qup = inflows[i] + if short_ts: + quc = inflows[i] + else: + quc = inflows[i + 1] + qdp = outflows[i] + + # Lower Interval + # Ck, WP, WPC, n, nCC, AREA, AREAC, R, So = self.calc_params(h_0) + # if h_0 > self.bfd: + # Twl = self.TwCC + # else: + # Twl = self.Bw + 2.0*self.z*h_0 + # q = ((1/(((WP*n)+(WPC*nCC))/(WP+WPC))) * (AREA+AREAC) * (R**(2./3.)) * np.sqrt(So)) + Ck = np.interp(h_0, self.geometry['stage'], self.geometry['celerity']) + Twl = np.interp(h_0, self.geometry['stage'], self.geometry['top_width']) + q = np.interp(h_0, self.geometry['stage'], self.geometry['discharge']) + + if Ck > 0.000000: + Km = max(dt, dx / Ck) + else: + Km = dt + + X = 0.5*(1-(Qj_0/(2.0*Twl*self.So*Ck*dx))) + X = min(0.5, max(0.0, X)) + + D = (Km*(1.000 - X) + dt/2.0000) + + C1 = (Km*X + dt/2.000000)/D + C2 = (dt/2.0000 - Km*X)/D + C3 = (Km*(1.00000000-X)-dt/2.000000)/D + C4 = (lateral[i]*dt)/D + Qj_0 = ((C1*quc) + (C2*qup) + (C3*qdp) + C4) - q + + # Upper Interval + # Ck, WP, WPC, n, nCC, AREA, AREAC, R, So = self.calc_params(h) + # if h > self.bfd: + # Twl = self.TwCC + # else: + # Twl = self.Bw + 2.0*self.z*h + # q = ((1/(((WP*n)+(WPC*nCC))/(WP+WPC))) * (AREA+AREAC) * (R**(2./3.)) * np.sqrt(So)) + Ck = np.interp(h, self.geometry['stage'], self.geometry['celerity']) + Twl = np.interp(h, self.geometry['stage'], self.geometry['top_width']) + q = np.interp(h, self.geometry['stage'], self.geometry['discharge']) + + if Ck > 0.000000: + Km = max(dt, dx / Ck) + else: + Km = dt + + if Ck > 0.0: + X = 0.5*(1-(((C1*qup)+(C2*quc)+(C3*qdp) + C4)/(2.0*Twl*self.slope*Ck*dx))) + X = min(0.5, max(0.25, X)) + else: + X = 0.5 + + D = (Km*(1.000 - X) + dt/2.0000) + + C1 = (Km*X + dt/2.000000)/D + C2 = (dt/2.0000 - Km*X)/D + C3 = (Km*(1.00000000-X)-dt/2.000000)/D + C4 = (lateral[i]*dt)/D + Qj = ((C1*quc) + (C2*qup) + (C3*qdp) + C4) - q + + # Secant Method + if Qj - Qj_0 != 0: + h_1 = h - ((Qj * (h_0 - h))/(Qj_0 - Qj)) + if h_1 < 0.0: + h_1 = h + else: + h_1 = h + + if h > 0.0: + rerror = abs((h_1 - h) / h) + aerror = abs(h_1 - h) + else: + rerror = 0.0 + aerror = 0.9 + + h_0 = max(0.0, h) + h = max(0.0, h_1) + iter += 1 + + if h < mindepth: + break + + if iter >= max_iter: + tries += 1 + if tries <= 4: + h = h * 1.33 + h_0 = h_0 * 0.67 + max_iter += 25 + + qdc = ((C1*qup)+(C2*quc)+(C3*qdp) + C4) + if qdc < 0.0: + qdc = 0.0 + outflows.append(qdc) + + return np.array(outflows) class SigmoidalReach(BaseReach): From 62fd8ca054e5ad26e972a2ca93546c67036feb0f Mon Sep 17 00:00:00 2001 From: sclaw Date: Wed, 17 Jul 2024 09:10:22 -0400 Subject: [PATCH 24/43] move methods and cleanup --- src/muskingumcunge/reach.py | 376 ++++++++++++------------------------ 1 file changed, 127 insertions(+), 249 deletions(-) diff --git a/src/muskingumcunge/reach.py b/src/muskingumcunge/reach.py index 50600bb..f23beab 100644 --- a/src/muskingumcunge/reach.py +++ b/src/muskingumcunge/reach.py @@ -173,6 +173,131 @@ def route_hydrograph_mct(self, inflows, dt, max_iter=1000): return np.array(outflows) + def route_hydrograph_wrf_secant(self, inflows, dt, lateral=None, max_iter=1000, initial_outflow=None): + # Secant method solver from WRF-HYDRO + dt = dt * 60 * 60 # MUST DELETE AFTER DEBUGGING + dx = self.reach_length + outflows = list() + if initial_outflow is not None: + outflows.append(initial_outflow) + else: + outflows.append(inflows[0]) + assert max(inflows) < max(self.geometry['discharge']), 'Rating Curve does not cover range of flowrates in hydrograph' + + if lateral is None: + lateral = np.zeros_like(inflows) + lateral = (lateral[:-1] + lateral[1:]) / 2 + lateral = np.append(lateral, 0) + + # local variables + mindepth = 0.01 + max_iter = 100 + short_ts = True + + # initialize secant method + depth = np.interp(outflows[0], self.geometry['discharge'], self.geometry['stage']) + h = (depth * 1.33) + mindepth + h_0 = (depth * 0.67) + + for i in range(len(inflows) - 1): + if not ((inflows[i] > 0.0) or (inflows[i + 1] > 0.0) or (outflows[i] > 0.0)): + outflows.append(0.0) + continue + Qj_0 = 0.0 + iter = 0 + rerror = 1.0 + aerror = 0.01 + tries = 0 + + while rerror > 0.01 and aerror >= mindepth and iter <= max_iter: + qup = inflows[i] + if short_ts: + quc = inflows[i] + else: + quc = inflows[i + 1] + qdp = outflows[i] + + # Lower Interval + Ck = np.interp(h_0, self.geometry['stage'], self.geometry['celerity']) + Twl = np.interp(h_0, self.geometry['stage'], self.geometry['top_width']) + q = np.interp(h_0, self.geometry['stage'], self.geometry['discharge']) + + if Ck > 0.000000: + Km = max(dt, dx / Ck) + else: + Km = dt + + X = 0.5*(1-(Qj_0/(2.0*Twl*self.So*Ck*dx))) + X = min(0.5, max(0.0, X)) + + D = (Km*(1.000 - X) + dt/2.0000) + + C1 = (Km*X + dt/2.000000)/D + C2 = (dt/2.0000 - Km*X)/D + C3 = (Km*(1.00000000-X)-dt/2.000000)/D + C4 = (lateral[i]*dt)/D + Qj_0 = ((C1*quc) + (C2*qup) + (C3*qdp) + C4) - q + + # Upper Interval + Ck = np.interp(h, self.geometry['stage'], self.geometry['celerity']) + Twl = np.interp(h, self.geometry['stage'], self.geometry['top_width']) + q = np.interp(h, self.geometry['stage'], self.geometry['discharge']) + + if Ck > 0.000000: + Km = max(dt, dx / Ck) + else: + Km = dt + + if Ck > 0.0: + X = 0.5*(1-(((C1*qup)+(C2*quc)+(C3*qdp) + C4)/(2.0*Twl*self.slope*Ck*dx))) + X = min(0.5, max(0.25, X)) + else: + X = 0.5 + + D = (Km*(1.000 - X) + dt/2.0000) + + C1 = (Km*X + dt/2.000000)/D + C2 = (dt/2.0000 - Km*X)/D + C3 = (Km*(1.00000000-X)-dt/2.000000)/D + C4 = (lateral[i]*dt)/D + Qj = ((C1*quc) + (C2*qup) + (C3*qdp) + C4) - q + + # Secant Method + if Qj - Qj_0 != 0: + h_1 = h - ((Qj * (h_0 - h))/(Qj_0 - Qj)) + if h_1 < 0.0: + h_1 = h + else: + h_1 = h + + if h > 0.0: + rerror = abs((h_1 - h) / h) + aerror = abs(h_1 - h) + else: + rerror = 0.0 + aerror = 0.9 + + h_0 = max(0.0, h) + h = max(0.0, h_1) + iter += 1 + + if h < mindepth: + break + + if iter >= max_iter: + tries += 1 + if tries <= 4: + h = h * 1.33 + h_0 = h_0 * 0.67 + max_iter += 25 + + qdc = ((C1*qup)+(C2*quc)+(C3*qdp) + C4) + if qdc < 0.0: + qdc = 0.0 + outflows.append(qdc) + + return np.array(outflows) + def optimize_route_params(self, inflows, dt): # Ponce method qref = (inflows.max() + inflows.min()) / 2 @@ -427,7 +552,7 @@ def generate_geometry(self): geom[param] = np.zeros(geom['stage'].shape) for ind, h_0 in enumerate(geom['stage']): - Ck, WP, WPC, n, nCC, AREA, AREAC, R, So = self.calc_params(h_0) + Ck, WP, WPC, n, nCC, AREA, AREAC, R, So = self.get_geom_at_stage(h_0) if h_0 > self.bfd: Twl = self.TwCC else: @@ -443,118 +568,7 @@ def generate_geometry(self): self.geometry = geom - - def generate_geometry_old(self): - geom = {'stage': None, - 'top_width': None, - 'area': None, - 'wetted_perimeter': None, - 'hydraulic_radius': None, - 'mannings_n': None, - 'discharge': None} - - geom['stage'] = np.linspace(0, self.max_stage, self.resolution) - - # rename variables to match WRF code - Bw = self.Bw - Tw = self.Bw + (2 * self.bfd * self.z) - TwCC = self.TwCC - nCC = self.nCC - So = self.So - n = self.n - z = self.z - - AREA, AREAC = None, None - Z = None # trapezoid distance - R, RC = None, None - WP, WPC = None, None - h_0 = geom['stage'] - bfd = self.bfd - - Ck = None - Twl = Bw + 2.0*z*h_0 # Top width at simulated flow - - compound_mask = (geom['stage'] > self.bfd) - - # Trapezoidal Channel - AREA = (Bw + h_0 * z ) * h_0 - WP = (Bw + 2.0 * h_0 * np.sqrt(1.0 + z*z)) - R = AREA / WP - - # Compound Channel - AREA[compound_mask] = ((Bw + bfd * z) * bfd) - AREAC = np.zeros(AREA.shape) - AREAC[compound_mask] = (TwCC * (h_0 -bfd))[compound_mask] - WP[compound_mask] = (Bw + 2.0 * bfd * np.sqrt(1.0 + z*z)) - WPC = np.zeros(WP.shape) - WPC[compound_mask] = (TwCC + (2.0 * (h_0-bfd)))[compound_mask] - R[compound_mask] = ((AREA + AREAC)/(WP +WPC))[compound_mask] - - R[np.isinf(R)] = 0.0 - - # Trapazoidal Channel - Ck = (np.sqrt(So)/n)*((5./3.)*R**(2./3.) - ((2./3.)*R**(5./3.)*(2.0*np.sqrt(1.0 + z*z)/(Bw+2.0*h_0*z)))) - Ck[h_0 == 0] = 0.0 - - # Compound Channel - Ck[compound_mask] = (((np.sqrt(So)/n)*((5./3.)*R**(2./3.) - ((2./3.)*R**(5./3.)*(2.0*np.sqrt(1.0 + z*z)/(Bw+2.0*bfd*z))))*AREA + ((np.sqrt(So)/(nCC))*(5./3.)*(h_0-bfd)**(2./3.))*AREAC)/(AREA+AREAC))[compound_mask] - Ck[Ck < 0] = 0.0 - - discharge = ((1/(((WP*n)+(WPC*nCC))/(WP+WPC))) * (AREA+AREAC) * (R**(2./3.)) * np.sqrt(So)) - - # log to class variables - self.geometry['stage'] = geom['stage'] - self.geometry['top_width'] = Twl - self.geometry['celerity'] = Ck - self.geometry['discharge'] = discharge - self.geometry['log_q'] = np.log(self.geometry['discharge']) - self.geometry['log_width'] = np.log(self.geometry['top_width']) - - - # geom['top_width'] = self.Bw + (self.z * geom['stage'] * 2) - # geom['area'] = (self.Bw + geom['stage'] * self.z) * geom['stage'] - # geom['wetted_perimeter'] = (self.Bw + 2.0 * geom['stage'] * np.sqrt(1.0 + self.z*self.z)) - # geom['hydraulic_radius'] = geom['area'] / geom['wetted_perimeter'] - # read_celerity = lambda y, r, n: (np.sqrt(self.So)/n)*((5./3.)*r**(2./3.) - ((2./3.)*r**(5./3.)*(2.0*np.sqrt(1.0 + self.z*self.z)/(self.Bw+2.0*y*self.z)))) - # geom['celerity'] = read_celerity(geom['stage'], geom['hydraulic_radius'], self.n) - # geom['mannings_n'] = np.repeat(self.n, self.resolution) - # geom['discharge'] = (1.0 / geom['mannings_n']) * geom['area'] * (geom['hydraulic_radius'] ** (2.0 / 3.0)) * (self.So ** 0.5) - - # # Compound Channel - # mask = (geom['stage'] > self.bfd) - # area_trap = (self.Bw + self.bfd * self.z) * self.bfd - # wetted_trap = (self.Bw + 2.0 * self.bfd * np.sqrt(1.0 + self.z*self.z)) - # hr_trap = area_trap / wetted_trap - # # celerity_trap = read_celerity(self.bfd, hr_trap, self.n) - # area_CC = (geom['stage'] - self.bfd) * self.TwCC - # wetted_CC = (self.TwCC + (2.0 * (geom['stage'] - self.bfd))) - # radius = (area_CC + area_trap) / (wetted_CC + wetted_trap) - # celerity_trap = read_celerity(self.bfd, radius, self.n) # Is this an error? I feel like radius should be capped at the bankfull radius? - # # celerity_cc = read_celerity(geom['stage'], radius, self.nCC) - # celerity_cc = (np.sqrt(self.So)/self.nCC)*((5./3.)*(geom['stage'] - self.bfd)**(2./3.)) - # celerity_cc = ((celerity_trap * area_trap) + (celerity_cc * area_CC)) / (area_trap + area_CC) - # n_cc = ((self.n * wetted_trap) + (wetted_CC * self.nCC)) / (wetted_trap + wetted_CC) - # # discharge_cc = (1.0 / self.nCC) * area_CC * ((geom['stage'] - self.bfd) ** (2.0 / 3.0)) * (self.So ** 0.5) - # # discharge_trap = (1.0 / self.n) * area_trap * ((area_trap / wetted_trap) ** (2.0 / 3.0)) * (self.So ** 0.5) - # discharge_cc = (1.0 / self.nCC) * (area_CC + area_trap) * (radius ** (2.0 / 3.0)) * (self.So ** 0.5) - - # geom['area'][mask] = area_CC[mask] + area_trap - # geom['wetted_perimeter'][mask] = wetted_CC[mask] - # geom['hydraulic_radius'][mask] = radius[mask] - # geom['celerity'][mask] = celerity_cc[mask] - # if np.any(geom['celerity'] < 0) or np.any(np.isnan(geom['celerity'])) or np.any(np.isinf(geom['celerity'])): - # print('bad celerity') - # geom['mannings_n'][mask] = n_cc[mask] - # geom['top_width'][mask] = self.TwCC - # # geom['discharge'][mask] = discharge_cc[mask] + discharge_trap - # geom['discharge'][mask] = discharge_cc[mask] - - # geom['discharge'] = (1.0 / geom['mannings_n']) * geom['area'] * (geom['hydraulic_radius'] ** (2.0 / 3.0)) * (self.So ** 0.5) - # geom['log_q'] = np.log(geom['discharge']) - # geom['log_width'] = np.log(geom['top_width']) - # self.geometry = geom - - def calc_params(self, h_0): + def get_geom_at_stage(self, h_0): Bw = self.Bw bfd = self.bfd z = self.z @@ -590,142 +604,6 @@ def calc_params(self, h_0): Ck = 0.0 return Ck, WP, WPC, n, nCC, AREA, AREAC, R, So - def route_hydrograph(self, inflows, dt, lateral=None, max_iter=1000, initial_outflow=None): - dt = dt * 60 * 60 # MUST DELETE AFTER DEBUGGING - dx = self.reach_length - outflows = list() - if initial_outflow is not None: - outflows.append(initial_outflow) - else: - outflows.append(inflows[0]) - assert max(inflows) < max(self.geometry['discharge']), 'Rating Curve does not cover range of flowrates in hydrograph' - - if lateral is None: - lateral = np.zeros_like(inflows) - lateral = (lateral[:-1] + lateral[1:]) / 2 - lateral = np.append(lateral, 0) - - # local variables - mindepth = 0.01 - max_iter = 100 - short_ts = True - - # initialize secant method - depth = np.interp(outflows[0], self.geometry['discharge'], self.geometry['stage']) - h = (depth * 1.33) + mindepth - h_0 = (depth * 0.67) - - for i in range(len(inflows) - 1): - if not ((inflows[i] > 0.0) or (inflows[i + 1] > 0.0) or (outflows[i] > 0.0)): - outflows.append(0.0) - continue - Qj_0 = 0.0 - iter = 0 - rerror = 1.0 - aerror = 0.01 - tries = 0 - - while rerror > 0.01 and aerror >= mindepth and iter <= max_iter: - qup = inflows[i] - if short_ts: - quc = inflows[i] - else: - quc = inflows[i + 1] - qdp = outflows[i] - - # Lower Interval - # Ck, WP, WPC, n, nCC, AREA, AREAC, R, So = self.calc_params(h_0) - # if h_0 > self.bfd: - # Twl = self.TwCC - # else: - # Twl = self.Bw + 2.0*self.z*h_0 - # q = ((1/(((WP*n)+(WPC*nCC))/(WP+WPC))) * (AREA+AREAC) * (R**(2./3.)) * np.sqrt(So)) - Ck = np.interp(h_0, self.geometry['stage'], self.geometry['celerity']) - Twl = np.interp(h_0, self.geometry['stage'], self.geometry['top_width']) - q = np.interp(h_0, self.geometry['stage'], self.geometry['discharge']) - - if Ck > 0.000000: - Km = max(dt, dx / Ck) - else: - Km = dt - - X = 0.5*(1-(Qj_0/(2.0*Twl*self.So*Ck*dx))) - X = min(0.5, max(0.0, X)) - - D = (Km*(1.000 - X) + dt/2.0000) - - C1 = (Km*X + dt/2.000000)/D - C2 = (dt/2.0000 - Km*X)/D - C3 = (Km*(1.00000000-X)-dt/2.000000)/D - C4 = (lateral[i]*dt)/D - Qj_0 = ((C1*quc) + (C2*qup) + (C3*qdp) + C4) - q - - # Upper Interval - # Ck, WP, WPC, n, nCC, AREA, AREAC, R, So = self.calc_params(h) - # if h > self.bfd: - # Twl = self.TwCC - # else: - # Twl = self.Bw + 2.0*self.z*h - # q = ((1/(((WP*n)+(WPC*nCC))/(WP+WPC))) * (AREA+AREAC) * (R**(2./3.)) * np.sqrt(So)) - Ck = np.interp(h, self.geometry['stage'], self.geometry['celerity']) - Twl = np.interp(h, self.geometry['stage'], self.geometry['top_width']) - q = np.interp(h, self.geometry['stage'], self.geometry['discharge']) - - if Ck > 0.000000: - Km = max(dt, dx / Ck) - else: - Km = dt - - if Ck > 0.0: - X = 0.5*(1-(((C1*qup)+(C2*quc)+(C3*qdp) + C4)/(2.0*Twl*self.slope*Ck*dx))) - X = min(0.5, max(0.25, X)) - else: - X = 0.5 - - D = (Km*(1.000 - X) + dt/2.0000) - - C1 = (Km*X + dt/2.000000)/D - C2 = (dt/2.0000 - Km*X)/D - C3 = (Km*(1.00000000-X)-dt/2.000000)/D - C4 = (lateral[i]*dt)/D - Qj = ((C1*quc) + (C2*qup) + (C3*qdp) + C4) - q - - # Secant Method - if Qj - Qj_0 != 0: - h_1 = h - ((Qj * (h_0 - h))/(Qj_0 - Qj)) - if h_1 < 0.0: - h_1 = h - else: - h_1 = h - - if h > 0.0: - rerror = abs((h_1 - h) / h) - aerror = abs(h_1 - h) - else: - rerror = 0.0 - aerror = 0.9 - - h_0 = max(0.0, h) - h = max(0.0, h_1) - iter += 1 - - if h < mindepth: - break - - if iter >= max_iter: - tries += 1 - if tries <= 4: - h = h * 1.33 - h_0 = h_0 * 0.67 - max_iter += 25 - - qdc = ((C1*qup)+(C2*quc)+(C3*qdp) + C4) - if qdc < 0.0: - qdc = 0.0 - outflows.append(qdc) - - return np.array(outflows) - class SigmoidalReach(BaseReach): geometry = {'stage': None, From a378d2eab07a860e8ca8299fbb76e44dee203b6b Mon Sep 17 00:00:00 2001 From: sclaw Date: Thu, 18 Jul 2024 10:49:43 -0400 Subject: [PATCH 25/43] worshopping above max stage extrapolation --- src/muskingumcunge/reach.py | 61 +++++++++++++++++++++++++++++++++++-- 1 file changed, 59 insertions(+), 2 deletions(-) diff --git a/src/muskingumcunge/reach.py b/src/muskingumcunge/reach.py index f23beab..e45cd4b 100644 --- a/src/muskingumcunge/reach.py +++ b/src/muskingumcunge/reach.py @@ -182,7 +182,64 @@ def route_hydrograph_wrf_secant(self, inflows, dt, lateral=None, max_iter=1000, outflows.append(initial_outflow) else: outflows.append(inflows[0]) - assert max(inflows) < max(self.geometry['discharge']), 'Rating Curve does not cover range of flowrates in hydrograph' + + if max(inflows) > max(self.geometry['discharge']): # putting in for loop for debugging + # for param in ['top_width', 'celerity', 'discharge', 'stage']: + # self.geometry[param] = np.append(self.geometry[param], np.linspace(self.geometry[param][-1] + 1, 9999, 100)) + # extrapolate out to max inflow + 1 + add_stages = np.linspace(0, 9 * self.geometry['stage'][-1], 200) + # define trend range for extrapolation as last meter of stages or 0.5 bankfull depth + start_trend = np.argmax(self.geometry['stage'] > self.geometry['stage'][-1] - 0.25) + # start_trend = int(len(self.geometry['stage']) - (len(self.geometry['stage']) / 40)) + # start_trend = int(len(self.geometry['stage']) - 6) + + params = ['top_width', 'celerity', 'discharge'] + scales = ['linear', 'loglog', 'linear'] + import matplotlib.pyplot as plt + import os + for param, scale in zip(params, scales): + fig, ax = plt.subplots() + ax.plot(self.geometry['stage'], self.geometry[param], c='k') + if scale == 'log': + end_trend = (np.log(self.geometry[param][-1]) - np.log(self.geometry[param][start_trend])) / (self.geometry['stage'][-1] - self.geometry['stage'][start_trend]) + # add_vals = np.exp(add_stages * end_trend) + self.geometry[param][-1] + add_vals = np.exp((add_stages * end_trend) + np.log(self.geometry[param][-1])) + elif scale == 'linear': + end_trend = (self.geometry[param][-1] - self.geometry[param][start_trend]) / (self.geometry['stage'][-1] - self.geometry['stage'][start_trend]) + add_vals = (add_stages * end_trend) + self.geometry[param][-1] + elif scale == 'logx': + end_trend = (self.geometry[param][-1] - self.geometry[param][start_trend]) / (np.log(self.geometry['stage'][-1]) - np.log(self.geometry['stage'][start_trend])) + tmp_s = np.log(add_stages + self.geometry['stage'][-1]) - np.log(self.geometry['stage'][-1]) + add_vals = (tmp_s * end_trend) + self.geometry[param][-1] + elif scale == 'loglog': + end_trend = (np.log(self.geometry[param][-1]) - np.log(self.geometry[param][start_trend])) / (np.log(self.geometry['stage'][-1]) - np.log(self.geometry['stage'][start_trend])) + tmp_s = np.log(add_stages + self.geometry['stage'][-1]) - np.log(self.geometry['stage'][-1]) + add_vals = np.exp((tmp_s * end_trend) + np.log(self.geometry[param][-1])) + + self.geometry[param] = np.append(self.geometry[param], add_vals) + ax.plot(np.append(self.geometry['stage'], add_stages + self.geometry['stage'][-1]), self.geometry[param], c='k', ls='--') + + max_s = add_stages[-1] + self.geometry['stage'][-1] + c_tmp, WP_tmp, WPC_tmp, n_tmp, nCC_tmp, AREA_tmp, AREAC_tmp, R_tmp, So_tmp = self.get_geom_at_stage(max_s) + q_tmp = ((1/(((WP_tmp*n_tmp)+(WPC_tmp*nCC_tmp))/(WP_tmp+WPC_tmp))) * (AREA_tmp+AREAC_tmp) * (R_tmp**(2./3.)) * np.sqrt(So_tmp)) + if param == 'discharge': + ax.scatter(max_s, q_tmp, ec='r', fc='none', s=20) + elif param == 'celerity': + ax.scatter(max_s, c_tmp, ec='r', fc='none', s=20) + + ax.set(xlabel='Stage (m)', ylabel=param) + out_dir = r"/users/k/l/klawson1/netfiles/ciroh/floodplainsData/retrospective/otter/2011/scott_working/extrapolation" + os.makedirs(out_dir, exist_ok=True) + existing_files = [int(f.split('.')[0]) for f in os.listdir(out_dir) if f.endswith('.png')] + if len(existing_files) == 0: + existing_files = [0] + existing_files = np.max(existing_files) + fig.savefig(os.path.join(out_dir, f"{existing_files + 1}.png"), dpi=300) + plt.close(fig) + + add_stages = add_stages + self.geometry['stage'][-1] + self.geometry['stage'] = np.append(self.geometry['stage'], add_stages) + assert max(inflows) < max(self.geometry['discharge']), 'Rating Curve does not cover range of flowrates in hydrograph' if lateral is None: lateral = np.zeros_like(inflows) @@ -227,7 +284,7 @@ def route_hydrograph_wrf_secant(self, inflows, dt, lateral=None, max_iter=1000, else: Km = dt - X = 0.5*(1-(Qj_0/(2.0*Twl*self.So*Ck*dx))) + X = 0.5*(1-(Qj_0/(2.0*Twl*self.slope*Ck*dx))) X = min(0.5, max(0.0, X)) D = (Km*(1.000 - X) + dt/2.0000) From 7b3a352f8a602fab3a6d6b67df4aafc6f32a9009 Mon Sep 17 00:00:00 2001 From: sclaw Date: Thu, 18 Jul 2024 12:04:01 -0400 Subject: [PATCH 26/43] update status outputs --- src/muskingumcunge/network.py | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/src/muskingumcunge/network.py b/src/muskingumcunge/network.py index 96a3805..99f2e71 100644 --- a/src/muskingumcunge/network.py +++ b/src/muskingumcunge/network.py @@ -1,6 +1,7 @@ from collections import defaultdict import pandas as pd import numpy as np +import time class Network: @@ -73,12 +74,16 @@ def load_initial_conditions(self, initial_conditions): def run_event(self, optimize_dx=True, conserve_mass=False, lat_addition='middle'): # calculate total inflow - volume_in = self.forcing_df.sum().sum() + np.sum([self.channel_outflows[c].sum() for c in self.headwaters]) + t_start = time.perf_counter() dt = (self.forcing_df.index[1] - self.forcing_df.index[0]).seconds / 3600 - counter = 1 + counter = 0 + pct5 = int(len(self.post_order) / 20) for node in self.post_order: - if counter % 100 == 0: - print(f"{counter} / {len(self.post_order)}") + counter += 1 + if counter % pct5 == 0: + pct_done = int(100 * counter / len(self.post_order)) + print(f'{pct_done}% done') + reach = self.reach_dict[node] if node in self.headwaters: continue @@ -108,7 +113,7 @@ def run_event(self, optimize_dx=True, conserve_mass=False, lat_addition='middle' init_out = self.init_outflows[node] else: init_out = None - routed = reach.route_hydrograph(routed, dt, lateral=l, initial_outflow=init_out) + routed = reach.route_hydrograph_wrf_secant(routed, dt, lateral=l, initial_outflow=init_out) except AssertionError as e: print(f"Error routing {node}: {e}") @@ -123,12 +128,7 @@ def run_event(self, optimize_dx=True, conserve_mass=False, lat_addition='middle' self.channel_outflows[node] = routed self.channel_outflows_stage[node] = np.interp(routed, reach.geometry['discharge'], reach.geometry['stage']) - # calculate total outflow - us_root = self.chlid_dict[self.root] - us_hydro = [self.channel_outflows[c] for c in us_root] - self.channel_outflows[self.root] = np.sum(us_hydro, axis=0) - volume_out = self.channel_outflows[self.root].sum() - print(f"Volume conservation = {round(100 *volume_out / volume_in, 2)} %") + print(f'Routing complete in {round(time.perf_counter() - t_start, 1)} seconds') # Post-processing out_df = pd.DataFrame(self.channel_outflows) From d1e3f9d2354a699f3accd9fad86a2d89307490e3 Mon Sep 17 00:00:00 2001 From: sclaw Date: Thu, 18 Jul 2024 12:04:41 -0400 Subject: [PATCH 27/43] adjust time resampling and stage resolution --- src/muskingumcunge/controller.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/muskingumcunge/controller.py b/src/muskingumcunge/controller.py index 89fe47b..61e270f 100644 --- a/src/muskingumcunge/controller.py +++ b/src/muskingumcunge/controller.py @@ -59,7 +59,7 @@ def load_geom(meta_df, source='NWM', geom_dir=None): bf = (tw - bw) / (z) if source == 'NWM': - reaches[r] = WRFCompound(bw, z, bf, tw_cc, ch_n, fp_n, slope, length, max_stage=20*bf, stage_resolution=999) + reaches[r] = WRFCompound(bw, z, bf, tw_cc, ch_n, fp_n, slope, length, max_stage=50*bf, stage_resolution=1000) elif source == 'NWM_Regression': tw = 2.44 * (da ** 0.34) a_ch = 0.75 * (da ** 0.53) @@ -67,7 +67,7 @@ def load_geom(meta_df, source='NWM', geom_dir=None): bw = ((2 * a_ch) / bf) - tw z = (tw - bw) / bf tw_cc = 3 * tw - reaches[r] = WRFCompound(bw, z, bf, tw_cc, ch_n, fp_n, slope, length, max_stage=20*bf, stage_resolution=999) + reaches[r] = WRFCompound(bw, z, bf, tw_cc, ch_n, fp_n, slope, length, max_stage=50*bf, stage_resolution=1000) elif source == 'HAND': try: tw = hand_tw[r].to_numpy() / length @@ -85,7 +85,7 @@ def load_geom(meta_df, source='NWM', geom_dir=None): print(f'Error loading HAND data for reach {r}. Defaulting to NWM channel') print(f'Error: {e}') geom_error_count += 1 - reaches[r] = WRFCompound(bw, z, bf, tw_cc, ch_n, fp_n, slope, length, max_stage=10*bf, stage_resolution=len(hand_tw)) + reaches[r] = WRFCompound(bw, z, bf, tw_cc, ch_n, fp_n, slope, length, max_stage=12*bf, stage_resolution=len(hand_tw)) elif source == 'MUSKINGUM': x = meta_df.loc[r, 'MusX'] k = meta_df.loc[r, 'MusK'] @@ -184,9 +184,10 @@ def execute(meta_path): head_df[lake_root] = downstream[lake_root] # resample to 5 minute increments - lateral_df = lateral_df.resample('5T').ffill() - head_df = head_df.resample('5T').ffill() - lake_df = lake_df.resample('5T').ffill() + resample_dt = pd.Timedelta('5m') + lateral_df = lateral_df.resample(resample_dt).ffill() + head_df = head_df.resample(resample_dt).ffill() + lake_df = lake_df.resample(resample_dt).ffill() # Error checking missing_forcings = list(set(meta_df.index).difference(set(lateral_df.columns))) From 88409a229c8953ad78594cd7512a2b55798c5596 Mon Sep 17 00:00:00 2001 From: sclaw Date: Thu, 18 Jul 2024 12:05:30 -0400 Subject: [PATCH 28/43] handle stages above max stage as vertical sidewalls --- src/muskingumcunge/reach.py | 75 +++++++++++-------------------------- 1 file changed, 21 insertions(+), 54 deletions(-) diff --git a/src/muskingumcunge/reach.py b/src/muskingumcunge/reach.py index e45cd4b..032f593 100644 --- a/src/muskingumcunge/reach.py +++ b/src/muskingumcunge/reach.py @@ -184,61 +184,28 @@ def route_hydrograph_wrf_secant(self, inflows, dt, lateral=None, max_iter=1000, outflows.append(inflows[0]) if max(inflows) > max(self.geometry['discharge']): # putting in for loop for debugging - # for param in ['top_width', 'celerity', 'discharge', 'stage']: - # self.geometry[param] = np.append(self.geometry[param], np.linspace(self.geometry[param][-1] + 1, 9999, 100)) - # extrapolate out to max inflow + 1 - add_stages = np.linspace(0, 9 * self.geometry['stage'][-1], 200) - # define trend range for extrapolation as last meter of stages or 0.5 bankfull depth - start_trend = np.argmax(self.geometry['stage'] > self.geometry['stage'][-1] - 0.25) - # start_trend = int(len(self.geometry['stage']) - (len(self.geometry['stage']) / 40)) - # start_trend = int(len(self.geometry['stage']) - 6) - - params = ['top_width', 'celerity', 'discharge'] - scales = ['linear', 'loglog', 'linear'] - import matplotlib.pyplot as plt - import os - for param, scale in zip(params, scales): - fig, ax = plt.subplots() - ax.plot(self.geometry['stage'], self.geometry[param], c='k') - if scale == 'log': - end_trend = (np.log(self.geometry[param][-1]) - np.log(self.geometry[param][start_trend])) / (self.geometry['stage'][-1] - self.geometry['stage'][start_trend]) - # add_vals = np.exp(add_stages * end_trend) + self.geometry[param][-1] - add_vals = np.exp((add_stages * end_trend) + np.log(self.geometry[param][-1])) - elif scale == 'linear': - end_trend = (self.geometry[param][-1] - self.geometry[param][start_trend]) / (self.geometry['stage'][-1] - self.geometry['stage'][start_trend]) - add_vals = (add_stages * end_trend) + self.geometry[param][-1] - elif scale == 'logx': - end_trend = (self.geometry[param][-1] - self.geometry[param][start_trend]) / (np.log(self.geometry['stage'][-1]) - np.log(self.geometry['stage'][start_trend])) - tmp_s = np.log(add_stages + self.geometry['stage'][-1]) - np.log(self.geometry['stage'][-1]) - add_vals = (tmp_s * end_trend) + self.geometry[param][-1] - elif scale == 'loglog': - end_trend = (np.log(self.geometry[param][-1]) - np.log(self.geometry[param][start_trend])) / (np.log(self.geometry['stage'][-1]) - np.log(self.geometry['stage'][start_trend])) - tmp_s = np.log(add_stages + self.geometry['stage'][-1]) - np.log(self.geometry['stage'][-1]) - add_vals = np.exp((tmp_s * end_trend) + np.log(self.geometry[param][-1])) - - self.geometry[param] = np.append(self.geometry[param], add_vals) - ax.plot(np.append(self.geometry['stage'], add_stages + self.geometry['stage'][-1]), self.geometry[param], c='k', ls='--') - - max_s = add_stages[-1] + self.geometry['stage'][-1] - c_tmp, WP_tmp, WPC_tmp, n_tmp, nCC_tmp, AREA_tmp, AREAC_tmp, R_tmp, So_tmp = self.get_geom_at_stage(max_s) - q_tmp = ((1/(((WP_tmp*n_tmp)+(WPC_tmp*nCC_tmp))/(WP_tmp+WPC_tmp))) * (AREA_tmp+AREAC_tmp) * (R_tmp**(2./3.)) * np.sqrt(So_tmp)) - if param == 'discharge': - ax.scatter(max_s, q_tmp, ec='r', fc='none', s=20) - elif param == 'celerity': - ax.scatter(max_s, c_tmp, ec='r', fc='none', s=20) - - ax.set(xlabel='Stage (m)', ylabel=param) - out_dir = r"/users/k/l/klawson1/netfiles/ciroh/floodplainsData/retrospective/otter/2011/scott_working/extrapolation" - os.makedirs(out_dir, exist_ok=True) - existing_files = [int(f.split('.')[0]) for f in os.listdir(out_dir) if f.endswith('.png')] - if len(existing_files) == 0: - existing_files = [0] - existing_files = np.max(existing_files) - fig.savefig(os.path.join(out_dir, f"{existing_files + 1}.png"), dpi=300) - plt.close(fig) + old_max_stage = self.geometry['stage'][-1] + add_res = 200 + add_stages = np.linspace(0, 9 * self.geometry['stage'][-1], add_res) + add_widths = np.repeat(self.geometry['top_width'][-1], add_res) + add_areas = (add_stages * add_widths) + self.geometry['area'][-1] + add_wp = (2 * add_stages) + self.geometry['wetted_perimeter'][-1] + add_hr = add_areas / add_wp + add_n = np.repeat(self.geometry['mannings_n'][-1], add_res) + add_q = (1 / add_n) * add_areas * (add_hr ** (2 / 3)) * (self.slope ** 0.5) + dpdy = np.repeat(2, add_res) + k_prime = (5 / 3) - ((2 / 3) * (add_areas / (add_widths * add_wp)) * dpdy) + add_c = k_prime * (add_q / add_areas) + + self.geometry['stage'] = np.append(self.geometry['stage'], add_stages + self.geometry['stage'][-1]) + self.geometry['top_width'] = np.append(self.geometry['top_width'], add_widths) + self.geometry['area'] = np.append(self.geometry['area'], add_areas) + self.geometry['wetted_perimeter'] = np.append(self.geometry['wetted_perimeter'], add_wp) + self.geometry['hydraulic_radius'] = np.append(self.geometry['hydraulic_radius'], add_hr) + self.geometry['mannings_n'] = np.append(self.geometry['mannings_n'], add_n) + self.geometry['discharge'] = np.append(self.geometry['discharge'], add_q) + self.geometry['celerity'] = np.append(self.geometry['celerity'], add_c) - add_stages = add_stages + self.geometry['stage'][-1] - self.geometry['stage'] = np.append(self.geometry['stage'], add_stages) assert max(inflows) < max(self.geometry['discharge']), 'Rating Curve does not cover range of flowrates in hydrograph' if lateral is None: From 762b1b36cb0a4bfea37798feeb53a7619829cd8b Mon Sep 17 00:00:00 2001 From: sclaw Date: Thu, 18 Jul 2024 12:06:41 -0400 Subject: [PATCH 29/43] cleanup --- src/muskingumcunge/reach.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/muskingumcunge/reach.py b/src/muskingumcunge/reach.py index 032f593..f8156ff 100644 --- a/src/muskingumcunge/reach.py +++ b/src/muskingumcunge/reach.py @@ -183,8 +183,7 @@ def route_hydrograph_wrf_secant(self, inflows, dt, lateral=None, max_iter=1000, else: outflows.append(inflows[0]) - if max(inflows) > max(self.geometry['discharge']): # putting in for loop for debugging - old_max_stage = self.geometry['stage'][-1] + if max(inflows) > max(self.geometry['discharge']): add_res = 200 add_stages = np.linspace(0, 9 * self.geometry['stage'][-1], add_res) add_widths = np.repeat(self.geometry['top_width'][-1], add_res) From 73a60438697f19c6e6a20dd55fd01ab608120394 Mon Sep 17 00:00:00 2001 From: sclaw Date: Wed, 24 Jul 2024 08:55:58 -0400 Subject: [PATCH 30/43] stable initial commit --- samples/basin_analysis/preprocessing.py | 214 ++++++++++++++++++++++++ 1 file changed, 214 insertions(+) create mode 100644 samples/basin_analysis/preprocessing.py diff --git a/samples/basin_analysis/preprocessing.py b/samples/basin_analysis/preprocessing.py new file mode 100644 index 0000000..a462d34 --- /dev/null +++ b/samples/basin_analysis/preprocessing.py @@ -0,0 +1,214 @@ +import pandas as pd +import numpy as np +import json +import os +from muskingumcunge.controller import create_project +import netCDF4 as nc + + +# Load data +# paths = r"/users/k/l/klawson1/netfiles/ciroh/floodplainsData/retrospective/lamoille/2002/mc_run.json" +# paths = r"/users/k/l/klawson1/netfiles/ciroh/floodplainsData/retrospective/lamoille/2019/mc_run.json" +paths = r"/users/k/l/klawson1/netfiles/ciroh/floodplainsData/retrospective/otter/2011/mc_run.json" +# paths = r"/users/k/l/klawson1/netfiles/ciroh/floodplainsData/retrospective/otter/2019/mc_run.json" +# paths = r"/users/k/l/klawson1/netfiles/ciroh/floodplainsData/retrospective/lewis/1996/mc_run.json" + + +def clean(df): + if 'Unnamed: 0' in df.columns: + df = df.drop(columns='Unnamed: 0') + if 'datetime' in df.columns: + df = df[df['datetime'].notnull()] + df["datetime"] = df["datetime"].map(lambda x: x + " 00:00:00" if len(x) == 10 else x) + df["datetime"] = pd.to_datetime(df["datetime"]) + df = df.set_index("datetime") + return df + +def find_root(reaches): + dataset = nc.Dataset(param_path) + comids = np.array(dataset.variables['link'][:]) + to = np.array(dataset.variables['to'][:]) + outflow_dict = {us: ds for us, ds in zip(comids, to)} + r1 = reaches[0] # just need to trace one reach to outlet + outflow = outflow_dict[r1] + while outflow in reaches: + r1 = outflow + outflow = outflow_dict[r1] + return r1 + +def extract_from_root(root): + print(f'Extracting data from root {root}') + print(f'Loading parameters from {param_path}') + dataset = nc.Dataset(param_path) + comids = np.array(dataset.variables['link'][:]) + to = np.array(dataset.variables['to'][:]) + + print('traversing network...') + keep = list() + q = [root] + while len(q) > 0: + cur = q.pop() + children = comids[to == cur] + q.extend(children) + keep.extend(children) + + print('extracting subset...') + indices = np.where(np.isin(comids, keep)) + data = {} + for var_name in ['link', 'to', 'n', 'nCC', 'order', 'Length', 'So', 'ChSlp', 'BtmWdth', 'TopWdthCC', 'TopWdth', 'NHDWaterbodyComID', 'MusK', 'MusX', 'Kchan']: + var_data = dataset.variables[var_name][:] + if var_data.ndim == 0: + continue + elif var_data.ndim > 1: + data[var_name] = var_data[indices, :] + else: + data[var_name] = var_data[indices] + + # Close the dataset + dataset.close() + + print('creating DataFrame...') + # Create a DataFrame from the dictionary + df = pd.DataFrame(data) + return df + + +run_dir = r'/netfiles/ciroh/floodplainsData/retrospective/basin_run_files' +all_paths = [os.path.join(run_dir, f) for f in os.listdir(run_dir) if f.endswith('.json')] + +for paths in all_paths: + print(f'Processing {paths}') + + # load data + if not os.path.exists(paths): + create_project(paths) # Optional make template + with open(paths, 'r') as f: + paths = json.load(f) + lateral_path = paths['lateral_path'] + head_path = paths['head_path'] + lake_path = paths['lake_path'] + out_path = paths['meta_path'] + if os.path.exists(out_path): + continue + meta_path = r"/users/k/l/klawson1/netfiles/ciroh/floodplainsData/runs/NWM/network/reach_data.csv" + param_path = r"/users/k/l/klawson1/netfiles/ciroh/floodplainsData/runs/NWM/working/RouteLink_CONUS.nc" + da_path = r"/users/k/l/klawson1/netfiles/ciroh/floodplainsData/runs/NWM/working/drain_areas.csv" + geom_dir = paths['geometry_dir'] + ds_hydrographs = os.path.join(paths['base_dir'], 'downstream_hydrographs.csv') + + ds = pd.read_csv(ds_hydrographs) + ds = clean(ds) + reaches = [int(i) for i in ds.columns] + root = find_root(reaches) + meta_df = extract_from_root(root) + meta_df = meta_df.sort_values("link") + meta_df = meta_df.rename(columns={"link": paths['id_field'], "to": paths['to_id_field'], 'Length': 'length', 'So': 'slope'}) + meta_df = meta_df.drop_duplicates() # remove completely duplicated rows + meta_df = meta_df.set_index("comid") + + # load DAs + das = pd.read_csv(da_path) + das = das.set_index("ID") + meta_df = meta_df.join(das, how='left') + replace_nan = lambda x: np.nanmean(meta_df['TotDASqKm']) if np.isnan(x) else x + meta_df['TotDASqKm'] = meta_df['TotDASqKm'].map(replace_nan) # maybe should update later + + # Remove divergences by keeping largest connection first + count = len(meta_df) + meta_df = meta_df.sort_values('order', ascending=False) + meta_df = meta_df[~meta_df.index.duplicated(keep='first')].copy().sort_index() + print(f'Removed {count - len(meta_df)} divergences') + + # Check for single root + roots = list(set(meta_df['ds_comid']).difference(set(meta_df.index))) + if len(roots) > 1: + raise RuntimeError(f'Multiple roots detected: {roots}') + + # Handle missing data by grabbing next ds reach + missing = meta_df[meta_df[['n', 'nCC', 'slope', 'ChSlp', 'BtmWdth', 'TopWdthCC', 'TopWdth', 'TotDASqKm']].isnull().values.any(axis=1)].index + print(f'Handling missing data for {len(missing)} reaches') + print(missing) + for col in missing: + isna = True + ds = col + while isna: + ds = meta_df.loc[ds, 'ds_comid'] + if ds in missing: + continue + else: + isna = False + if ds == roots[0]: + meta_df.loc[col, ['n', 'ncc', 'TotDASqKm', 'slope']] = [0.06, 0.12, meta_df['TotDASqKm'].max(), np.exp(np.mean(np.log(meta_df['slope'])))] + else: + meta_df.loc[col, ['n', 'ncc', 'TotDASqKm', 'slope']] = meta_df.loc[ds, ['n', 'ncc', 'TotDASqKm', 'slope']] + meta_df.loc[col, 'length'] = meta_df['length'].mean() + + # Add lakes to headwater forcings + try: + lake_df = pd.read_csv(lake_path) + tmp_cols = list(lake_df.columns) + tmp_cols.remove('Unnamed: 0') + tmp_cols.remove('datetime') + lakes = [float(c) for c in tmp_cols] + force_lakes = list() + for l in lakes: + # Find root of each lake + lake_reaches = list(meta_df[meta_df['NHDWaterbodyComID'] == l].index) + lake_root = lake_reaches[0] + outflow = meta_df.loc[lake_root, 'ds_comid'] + while outflow in lake_reaches: + lake_root = meta_df.loc[lake_root, 'ds_comid'] + outflow = meta_df.loc[lake_root, 'ds_comid'] + force_lakes.append(str(int(lake_root))) + # ds[force_lakes].to_csv(os.path.join(paths['base_dir'], 'lake2.csv'), index=False) + except: + force_lakes = list() + + # Prune upstream of headwater forcings + head_df = pd.read_csv(head_path) # reaches we don't want to route through + tmp_cols = list(head_df.columns) + tmp_cols.remove('Unnamed: 0') + tmp_cols.remove('datetime') + heads = [float(c) for c in tmp_cols] + small = list(meta_df[meta_df['TotDASqKm'] < 5.2].index) # also treat small reaches as headwaters + heads.extend(small) + heads.extend([float(i) for i in force_lakes]) + prune_counter = 0 + for h in heads: + us = list(meta_df[meta_df['ds_comid'] == h].index) + while len(us) > 0: + tmp = us.pop() + us.extend(list(meta_df[meta_df['ds_comid'] == tmp].index)) + meta_df = meta_df.drop(index=tmp).copy() + prune_counter += 1 + print(f'Removed {prune_counter} reaches to align with headwater forcings') + + # # Prune upstream of lake forcings + # lake_df = pd.read_csv(lake_path) + # tmp_cols = list(lake_df.columns) + # tmp_cols.remove('Unnamed: 0') + # tmp_cols.remove('datetime') + # lakes = [float(c) for c in tmp_cols] + # prune_counter = 0 + # for l in lakes: + # # Find root of each lake + # lake_reaches = list(meta_df[meta_df['NHDWaterbodyComID'] == l].index) + # lake_root = lake_reaches[0] + # outflow = meta_df.loc[lake_root, 'ds_comid'] + # while outflow in lake_reaches: + # lake_root = meta_df.loc[lake_root, 'ds_comid'] + # outflow = meta_df.loc[lake_root, 'ds_comid'] + # lake_root_cache = meta_df.loc[lake_root] + # us = list(meta_df[meta_df['ds_comid'] == lake_root].index) + # while len(us) > 0: + # tmp = us.pop() + # us.extend(list(meta_df[meta_df['ds_comid'] == tmp].index)) + # meta_df = meta_df.drop(index=tmp).copy() + # prune_counter += 1 + # meta_df.loc[lake_root] = lake_root_cache + # print(f'Removed {prune_counter} reaches to align with lake forcings') + + # Export metadata + meta_df['comid'] = meta_df.index + meta_df = meta_df[['comid', 'ds_comid', 'n', 'nCC', 'order', 'TotDASqKm', 'length', 'slope', 'ChSlp', 'BtmWdth', 'TopWdthCC', 'TopWdth', 'NHDWaterbodyComID', 'MusK', 'MusX', 'Kchan']] + meta_df.to_csv(out_path, index=False) From aed1f4fbfdd4a5d12bfc8e97190dd98f981868fa Mon Sep 17 00:00:00 2001 From: sclaw Date: Wed, 24 Jul 2024 09:01:27 -0400 Subject: [PATCH 31/43] stable initial commit --- samples/basin_analysis/make_runs.py | 39 +++++++++++++++++++++ samples/basin_analysis/run_all_basinwide.py | 8 +++++ 2 files changed, 47 insertions(+) create mode 100644 samples/basin_analysis/make_runs.py create mode 100644 samples/basin_analysis/run_all_basinwide.py diff --git a/samples/basin_analysis/make_runs.py b/samples/basin_analysis/make_runs.py new file mode 100644 index 0000000..00c494c --- /dev/null +++ b/samples/basin_analysis/make_runs.py @@ -0,0 +1,39 @@ +import os +import json + + +def create_project(data_path, run_path, geom_type): + template = { + "base_dir": data_path, + "meta_path": os.path.join(data_path, "reach_meta.csv"), + "lateral_path": os.path.join(data_path, "laterals.csv"), + "head_path": os.path.join(data_path, "headwaters.csv"), + "lake_path": os.path.join(data_path, "lake.csv"), + "geometry_dir": "/netfiles/ciroh/floodplainsData/runs/NWM/geometry", + "geometry_source": geom_type.upper(), + "optimize_dx": False, + "conserve_mass": False, + "lat_addition": "top", + "id_field": "comid", + "to_id_field": "ds_comid", + "mode": 'basin', + "out_dir": os.path.join(data_path, f'{geom_type.lower()}_outputs'), + } + os.makedirs(template['out_dir'], exist_ok=True) + with open(run_path, 'w') as f: + json.dump(template, f, indent=4) + +base_dir = r"/netfiles/ciroh/floodplainsData/retrospective" +runs_path = os.path.join(base_dir, 'basin_run_files') +catchments = ['lamoille', 'lewis', 'mettawee', 'missisquoi', 'otter', 'winooski'] +for catchment in catchments: + tmp_dir = os.path.join(base_dir, catchment) + folders = [f for f in os.listdir(tmp_dir) if os.path.isdir(os.path.join(tmp_dir, f))] + for f in folders: + data_path = os.path.join(tmp_dir, f) + geom_type = 'nwm' + run_path = os.path.join(runs_path, f'{catchment}_{f}_{geom_type}.json') + create_project(data_path, run_path, geom_type) + geom_type = 'hand' + run_path = os.path.join(runs_path, f'{catchment}_{f}_{geom_type}.json') + create_project(data_path, run_path, geom_type) diff --git a/samples/basin_analysis/run_all_basinwide.py b/samples/basin_analysis/run_all_basinwide.py new file mode 100644 index 0000000..12b4c91 --- /dev/null +++ b/samples/basin_analysis/run_all_basinwide.py @@ -0,0 +1,8 @@ +import os +from muskingumcunge.controller import handler + +run_dir = r'/netfiles/ciroh/floodplainsData/retrospective/basin_run_files' +all_paths = [os.path.join(run_dir, f) for f in os.listdir(run_dir) if f.endswith('.json')] + +for path in all_paths: + handler(path) \ No newline at end of file From 6d8b7442ebe702252ba0127298e98150e32cdec6 Mon Sep 17 00:00:00 2001 From: sclaw Date: Wed, 24 Jul 2024 12:38:03 -0400 Subject: [PATCH 32/43] add hail mary stages for the NWM reach --- src/muskingumcunge/reach.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/muskingumcunge/reach.py b/src/muskingumcunge/reach.py index f8156ff..e5e8d6e 100644 --- a/src/muskingumcunge/reach.py +++ b/src/muskingumcunge/reach.py @@ -571,6 +571,8 @@ def generate_geometry(self): 'discharge': None} geom['stage'] = np.linspace(0, self.max_stage, self.resolution) + far_out_stages = np.exp(np.linspace(np.log(self.max_stage), np.log(10 * self.max_stage), 50)) # add some hail mary stages + geom['stage'] = np.append(geom['stage'], far_out_stages) for param in ['top_width', 'area', 'wetted_perimeter', 'hydraulic_radius', 'mannings_n', 'discharge', 'celerity']: geom[param] = np.zeros(geom['stage'].shape) From f08fc747a3c9afdc7a2cb40d5392d94ca5a1ddef Mon Sep 17 00:00:00 2001 From: sclaw Date: Wed, 24 Jul 2024 12:38:37 -0400 Subject: [PATCH 33/43] try to make progress prints at 5% increments --- src/muskingumcunge/network.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/muskingumcunge/network.py b/src/muskingumcunge/network.py index 99f2e71..e8d5b20 100644 --- a/src/muskingumcunge/network.py +++ b/src/muskingumcunge/network.py @@ -77,7 +77,7 @@ def run_event(self, optimize_dx=True, conserve_mass=False, lat_addition='middle' t_start = time.perf_counter() dt = (self.forcing_df.index[1] - self.forcing_df.index[0]).seconds / 3600 counter = 0 - pct5 = int(len(self.post_order) / 20) + pct5 = int(len(self.post_order) / 20) + 1 for node in self.post_order: counter += 1 if counter % pct5 == 0: From 9f1dea6b989d0140f4dec361f481a2375f2037dd Mon Sep 17 00:00:00 2001 From: sclaw Date: Wed, 24 Jul 2024 12:39:55 -0400 Subject: [PATCH 34/43] change the way forcings are put into model --- samples/basin_analysis/preprocessing.py | 85 ++++++++++++++----------- src/muskingumcunge/controller.py | 55 +++++----------- 2 files changed, 63 insertions(+), 77 deletions(-) diff --git a/samples/basin_analysis/preprocessing.py b/samples/basin_analysis/preprocessing.py index a462d34..e808dab 100644 --- a/samples/basin_analysis/preprocessing.py +++ b/samples/basin_analysis/preprocessing.py @@ -22,6 +22,8 @@ def clean(df): df["datetime"] = df["datetime"].map(lambda x: x + " 00:00:00" if len(x) == 10 else x) df["datetime"] = pd.to_datetime(df["datetime"]) df = df.set_index("datetime") + if 'NA' in df.columns: + df = df.drop(columns='NA') return df def find_root(reaches): @@ -79,24 +81,29 @@ def extract_from_root(root): for paths in all_paths: print(f'Processing {paths}') - # load data + # Establish file structure and paths if not os.path.exists(paths): create_project(paths) # Optional make template with open(paths, 'r') as f: paths = json.load(f) lateral_path = paths['lateral_path'] - head_path = paths['head_path'] - lake_path = paths['lake_path'] + inflow_path = paths['inflow_path'] out_path = paths['meta_path'] - if os.path.exists(out_path): - continue - meta_path = r"/users/k/l/klawson1/netfiles/ciroh/floodplainsData/runs/NWM/network/reach_data.csv" - param_path = r"/users/k/l/klawson1/netfiles/ciroh/floodplainsData/runs/NWM/working/RouteLink_CONUS.nc" - da_path = r"/users/k/l/klawson1/netfiles/ciroh/floodplainsData/runs/NWM/working/drain_areas.csv" geom_dir = paths['geometry_dir'] - ds_hydrographs = os.path.join(paths['base_dir'], 'downstream_hydrographs.csv') - ds = pd.read_csv(ds_hydrographs) + ds_hydrographs = os.path.join(paths['base_dir'], 'downstream_hydrographs.csv') + head_path = os.path.join(paths['base_dir'], 'headwaters.csv') + lake_path = os.path.join(paths['base_dir'], 'lake.csv') + + meta_path = r"/netfiles/ciroh/floodplainsData/runs/NWM/network/reach_data.csv" + param_path = r"/netfiles/ciroh/floodplainsData/runs/NWM/working/RouteLink_CONUS.nc" + da_path = r"/netfiles/ciroh/floodplainsData/runs/NWM/working/drain_areas.csv" + if os.path.exists(out_path): + print(f'{out_path} already exists, skipping') + continue + + # Load data + ds = pd.read_csv(ds_hydrographs) # treating this as the model domain ds = clean(ds) reaches = [int(i) for i in ds.columns] root = find_root(reaches) @@ -154,13 +161,12 @@ def extract_from_root(root): for l in lakes: # Find root of each lake lake_reaches = list(meta_df[meta_df['NHDWaterbodyComID'] == l].index) - lake_root = lake_reaches[0] + lake_root = lake_reaches[0] # initialize outflow = meta_df.loc[lake_root, 'ds_comid'] while outflow in lake_reaches: - lake_root = meta_df.loc[lake_root, 'ds_comid'] + lake_root = meta_df.loc[lake_root, 'ds_comid'] # update lake root outflow = meta_df.loc[lake_root, 'ds_comid'] force_lakes.append(str(int(lake_root))) - # ds[force_lakes].to_csv(os.path.join(paths['base_dir'], 'lake2.csv'), index=False) except: force_lakes = list() @@ -183,32 +189,35 @@ def extract_from_root(root): prune_counter += 1 print(f'Removed {prune_counter} reaches to align with headwater forcings') - # # Prune upstream of lake forcings - # lake_df = pd.read_csv(lake_path) - # tmp_cols = list(lake_df.columns) - # tmp_cols.remove('Unnamed: 0') - # tmp_cols.remove('datetime') - # lakes = [float(c) for c in tmp_cols] - # prune_counter = 0 - # for l in lakes: - # # Find root of each lake - # lake_reaches = list(meta_df[meta_df['NHDWaterbodyComID'] == l].index) - # lake_root = lake_reaches[0] - # outflow = meta_df.loc[lake_root, 'ds_comid'] - # while outflow in lake_reaches: - # lake_root = meta_df.loc[lake_root, 'ds_comid'] - # outflow = meta_df.loc[lake_root, 'ds_comid'] - # lake_root_cache = meta_df.loc[lake_root] - # us = list(meta_df[meta_df['ds_comid'] == lake_root].index) - # while len(us) > 0: - # tmp = us.pop() - # us.extend(list(meta_df[meta_df['ds_comid'] == tmp].index)) - # meta_df = meta_df.drop(index=tmp).copy() - # prune_counter += 1 - # meta_df.loc[lake_root] = lake_root_cache - # print(f'Removed {prune_counter} reaches to align with lake forcings') + # Check that all final network have valid geometry + geom_area_path = os.path.join(geom_dir, 'area.csv') + geom_area = pd.read_csv(geom_area_path) + valid_geom = geom_area.columns[~(geom_area == 0).all(axis=0)].astype(int) + missing_geom = list(set(meta_df.index).difference(set(valid_geom))) + if len(missing_geom) > 0: + print(f'Missing geometry for {len(missing_geom)} reaches') + print(missing_geom) + meta_df['valid_geom'] = meta_df.index.isin(valid_geom) + + # make an inflows file + headwaters = list(set(meta_df.index).difference(set(meta_df['ds_comid']))) + headwaters = [str(h) for h in headwaters] + try: + inflow_df = ds[headwaters].copy() + except KeyError: + missing = list(set(headwaters).difference(set(ds.columns))) + missing_ds = pd.DataFrame(index=missing, columns=['missing_ds']) + missing_ds['missing_ds'] = True + missing_ds.to_csv(os.path.join(paths['base_dir'], 'missing_ds.csv')) + raise RuntimeError(f'Missing downstream reaches in {ds_hydrographs}') + inflow_df.to_csv(inflow_path) + + # Make an initial conditions file + init_flow = pd.DataFrame(ds.loc[ds.index[0], meta_df.index.astype(str)]) + init_flow = init_flow.rename(columns={ds.index[0]: 'init_discharge'}) + init_flow.to_csv(paths['init_conditions_path']) # Export metadata meta_df['comid'] = meta_df.index - meta_df = meta_df[['comid', 'ds_comid', 'n', 'nCC', 'order', 'TotDASqKm', 'length', 'slope', 'ChSlp', 'BtmWdth', 'TopWdthCC', 'TopWdth', 'NHDWaterbodyComID', 'MusK', 'MusX', 'Kchan']] + meta_df = meta_df[['comid', 'ds_comid', 'n', 'nCC', 'order', 'TotDASqKm', 'length', 'slope', 'ChSlp', 'BtmWdth', 'TopWdthCC', 'TopWdth', 'NHDWaterbodyComID', 'MusK', 'MusX', 'Kchan', 'valid_geom']] meta_df.to_csv(out_path, index=False) diff --git a/src/muskingumcunge/controller.py b/src/muskingumcunge/controller.py index 61e270f..53b8b91 100644 --- a/src/muskingumcunge/controller.py +++ b/src/muskingumcunge/controller.py @@ -32,6 +32,7 @@ def clean_ts(df): if 'Unnamed: 0' in df.columns: df = df.drop(columns='Unnamed: 0') if 'datetime' in df.columns: + df = df[df['datetime'].notnull()] df["datetime"] = df["datetime"].map(lambda x: x + " 00:00:00" if len(x) == 10 else x) df["datetime"] = pd.to_datetime(df["datetime"]) df = df.set_index("datetime") @@ -59,7 +60,7 @@ def load_geom(meta_df, source='NWM', geom_dir=None): bf = (tw - bw) / (z) if source == 'NWM': - reaches[r] = WRFCompound(bw, z, bf, tw_cc, ch_n, fp_n, slope, length, max_stage=50*bf, stage_resolution=1000) + reaches[r] = WRFCompound(bw, z, bf, tw_cc, ch_n, fp_n, slope, length, max_stage=12*bf, stage_resolution=1000) elif source == 'NWM_Regression': tw = 2.44 * (da ** 0.34) a_ch = 0.75 * (da ** 0.53) @@ -67,7 +68,7 @@ def load_geom(meta_df, source='NWM', geom_dir=None): bw = ((2 * a_ch) / bf) - tw z = (tw - bw) / bf tw_cc = 3 * tw - reaches[r] = WRFCompound(bw, z, bf, tw_cc, ch_n, fp_n, slope, length, max_stage=50*bf, stage_resolution=1000) + reaches[r] = WRFCompound(bw, z, bf, tw_cc, ch_n, fp_n, slope, length, max_stage=12*bf, stage_resolution=1000) elif source == 'HAND': try: tw = hand_tw[r].to_numpy() / length @@ -158,36 +159,14 @@ def execute(meta_path): lateral_df = clean_ts(lateral_df) lateral_df = lateral_df.reindex(sorted(lateral_df.columns), axis=1) - head_df = pd.read_csv(paths['head_path']) - head_df = clean_ts(head_df) - head_df = head_df.reindex(sorted(head_df.columns), axis=1) - - lake_df = pd.read_csv(paths['lake_path']) - lake_df = clean_ts(lake_df) - lake_df = lake_df.reindex(sorted(lake_df.columns), axis=1) - - downstream = pd.read_csv(os.path.join(paths['base_dir'], 'downstream_hydrographs.csv')) - downstream = clean_ts(downstream) - downstream = downstream.reindex(sorted(downstream.columns), axis=1) - - # pre process lakes - tmp_cols = list(lake_df.columns) - lakes = [float(c) for c in tmp_cols] - for l in lakes: - # Find root of each lake - lake_reaches = list(meta_df[meta_df['NHDWaterbodyComID'] == l].index) - lake_root = lake_reaches[0] - outflow = meta_df.loc[lake_root, 'ds_comid'] - while outflow in lake_reaches: - lake_root = meta_df.loc[lake_root, 'ds_comid'] - outflow = meta_df.loc[lake_root, 'ds_comid'] - head_df[lake_root] = downstream[lake_root] + inflow_df = pd.read_csv(paths['inflow_path']) + inflow_df = clean_ts(inflow_df) + inflow_df = inflow_df.reindex(sorted(inflow_df.columns), axis=1) # resample to 5 minute increments - resample_dt = pd.Timedelta('5m') + resample_dt = pd.Timedelta(paths['dt']) lateral_df = lateral_df.resample(resample_dt).ffill() - head_df = head_df.resample(resample_dt).ffill() - lake_df = lake_df.resample(resample_dt).ffill() + inflow_df = inflow_df.resample(resample_dt).ffill() # Error checking missing_forcings = list(set(meta_df.index).difference(set(lateral_df.columns))) @@ -199,29 +178,27 @@ def execute(meta_path): reaches = load_geom(meta_df, source=paths['geometry_source'], geom_dir=paths['geometry_dir']) # export celerity curves - stages = {r: reaches[r].geometry['stage'] for r in reaches} + stages = {r: reaches[r].geometry['stage'] for r in reaches if type(reaches[r]) == CustomReach} stage_df = pd.DataFrame().from_dict(stages, orient='columns') stage_df.to_csv(os.path.join(paths['out_dir'], "stage.csv")) - celerities = {r: reaches[r].geometry['celerity'] for r in reaches} + celerities = {r: reaches[r].geometry['celerity'] for r in reaches if type(reaches[r]) == CustomReach} celerity_df = pd.DataFrame().from_dict(celerities, orient='columns') celerity_df.to_csv(os.path.join(paths['out_dir'], "celerity.csv")) - discharge = {r: reaches[r].geometry['discharge'] for r in reaches} + discharge = {r: reaches[r].geometry['discharge'] for r in reaches if type(reaches[r]) == CustomReach} discharge_df = pd.DataFrame().from_dict(discharge, orient='columns') discharge_df.to_csv(os.path.join(paths['out_dir'], "discharge.csv")) # Set up run network = Network(edge_dict) network.load_forcings(lateral_df) - network.load_headwater_forcings(head_df) + network.load_headwater_forcings(inflow_df) network.load_reaches(reaches) # Set up initial outflows - for n in network.post_order: - obs_out = downstream[n].to_numpy()[0] - if np.isnan(obs_out): - obs_out = None - network.init_outflows[n] = obs_out - network.load_initial_conditions(network.init_outflows) # kind of unnecessary + init_outflows = pd.read_csv(paths['init_conditions_path']) + init_outflows = init_outflows.set_index('comid') + init_outflows = {str(r): init_outflows.loc[r, 'init_discharge'] for r in init_outflows.index} + network.load_initial_conditions(init_outflows) # kind of unnecessary # Run model network.run_event(optimize_dx=paths['optimize_dx'], conserve_mass=paths['conserve_mass'], lat_addition=paths['lat_addition']) From 1ccea597bceb0b83942b8d22f3fb0a60829f0d09 Mon Sep 17 00:00:00 2001 From: sclaw Date: Thu, 25 Jul 2024 14:57:14 -0400 Subject: [PATCH 35/43] add time formatting and update forcing format --- samples/basin_analysis/make_runs.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/samples/basin_analysis/make_runs.py b/samples/basin_analysis/make_runs.py index 00c494c..4f36860 100644 --- a/samples/basin_analysis/make_runs.py +++ b/samples/basin_analysis/make_runs.py @@ -7,10 +7,11 @@ def create_project(data_path, run_path, geom_type): "base_dir": data_path, "meta_path": os.path.join(data_path, "reach_meta.csv"), "lateral_path": os.path.join(data_path, "laterals.csv"), - "head_path": os.path.join(data_path, "headwaters.csv"), - "lake_path": os.path.join(data_path, "lake.csv"), + "inflow_path": os.path.join(data_path, "inflows.csv"), + "init_conditions_path": os.path.join(data_path, "init_conditions.csv"), "geometry_dir": "/netfiles/ciroh/floodplainsData/runs/NWM/geometry", "geometry_source": geom_type.upper(), + "dt": "5m", "optimize_dx": False, "conserve_mass": False, "lat_addition": "top", From ec40b397bd1b5f50482ed1dfe08238475857f83d Mon Sep 17 00:00:00 2001 From: sclaw Date: Thu, 25 Jul 2024 14:58:43 -0400 Subject: [PATCH 36/43] update error geom handling and correct format for nan initial inflows --- src/muskingumcunge/controller.py | 31 ++++++++++++++++++++++--------- 1 file changed, 22 insertions(+), 9 deletions(-) diff --git a/src/muskingumcunge/controller.py b/src/muskingumcunge/controller.py index 53b8b91..79505e7 100644 --- a/src/muskingumcunge/controller.py +++ b/src/muskingumcunge/controller.py @@ -36,11 +36,14 @@ def clean_ts(df): df["datetime"] = df["datetime"].map(lambda x: x + " 00:00:00" if len(x) == 10 else x) df["datetime"] = pd.to_datetime(df["datetime"]) df = df.set_index("datetime") + df = df.astype(np.float64) + if 'NA' in df.columns: + df = df.drop(columns='NA') return df def load_geom(meta_df, source='NWM', geom_dir=None): reaches = dict() - geom_error_count = 0 + missing_geom = list() if source == 'HAND': hand_tw = pd.read_csv(os.path.join(geom_dir, 'area.csv')) hand_el = pd.read_csv(os.path.join(geom_dir, 'el.csv')) @@ -60,7 +63,7 @@ def load_geom(meta_df, source='NWM', geom_dir=None): bf = (tw - bw) / (z) if source == 'NWM': - reaches[r] = WRFCompound(bw, z, bf, tw_cc, ch_n, fp_n, slope, length, max_stage=12*bf, stage_resolution=1000) + reaches[r] = WRFCompound(bw, z, bf, tw_cc, ch_n, fp_n, slope, length, max_stage=100*bf, stage_resolution=10000) elif source == 'NWM_Regression': tw = 2.44 * (da ** 0.34) a_ch = 0.75 * (da ** 0.53) @@ -85,15 +88,18 @@ def load_geom(meta_df, source='NWM', geom_dir=None): # Error catching. Default to NWM channel print(f'Error loading HAND data for reach {r}. Defaulting to NWM channel') print(f'Error: {e}') - geom_error_count += 1 + missing_geom.append(r) reaches[r] = WRFCompound(bw, z, bf, tw_cc, ch_n, fp_n, slope, length, max_stage=12*bf, stage_resolution=len(hand_tw)) elif source == 'MUSKINGUM': x = meta_df.loc[r, 'MusX'] k = meta_df.loc[r, 'MusK'] reaches[r] = MuskingumReach(k, x) + + if len(missing_geom) > 0: + print(f'Error loading {len(missing_geom)} / {len(reaches)} reaches') + print(missing_geom) - print(f'Error loading {geom_error_count} / {len(reaches)} reaches') - return reaches + return reaches, missing_geom def execute_by_reach(meta_path): # Load data @@ -172,11 +178,15 @@ def execute(meta_path): missing_forcings = list(set(meta_df.index).difference(set(lateral_df.columns))) print(f'Missing forcing data for: {missing_forcings}') - network = Network(edge_dict) - # Set up reaches - reaches = load_geom(meta_df, source=paths['geometry_source'], geom_dir=paths['geometry_dir']) - + reaches, missing_geom = load_geom(meta_df, source=paths['geometry_source'], geom_dir=paths['geometry_dir']) + geom_status_df = pd.DataFrame(index=meta_df.index) + geom_status_df['valid_geom'] = ~meta_df.index.isin(missing_geom) + try: + geom_status_df.to_csv(os.path.join(paths['out_dir'], "geom_status.csv")) + except PermissionError: + print("Error writing geom_status.csv. Check permissions") + # export celerity curves stages = {r: reaches[r].geometry['stage'] for r in reaches if type(reaches[r]) == CustomReach} stage_df = pd.DataFrame().from_dict(stages, orient='columns') @@ -198,6 +208,9 @@ def execute(meta_path): init_outflows = pd.read_csv(paths['init_conditions_path']) init_outflows = init_outflows.set_index('comid') init_outflows = {str(r): init_outflows.loc[r, 'init_discharge'] for r in init_outflows.index} + for r in init_outflows: + if np.isnan(init_outflows[r]): + init_outflows[r] = None network.load_initial_conditions(init_outflows) # kind of unnecessary # Run model From 5f1ab2bcb5759d476b00fb68bffdef1933f251ba Mon Sep 17 00:00:00 2001 From: sclaw Date: Thu, 25 Jul 2024 15:03:00 -0400 Subject: [PATCH 37/43] update depths in wrf_secant method. correct qup and quc placement. more congruence with WRF-HYDRO fortran code. option to vectorize geometry. --- src/muskingumcunge/reach.py | 96 +++++++++++++++++++++++++++---------- 1 file changed, 72 insertions(+), 24 deletions(-) diff --git a/src/muskingumcunge/reach.py b/src/muskingumcunge/reach.py index e5e8d6e..02e5945 100644 --- a/src/muskingumcunge/reach.py +++ b/src/muskingumcunge/reach.py @@ -219,12 +219,12 @@ def route_hydrograph_wrf_secant(self, inflows, dt, lateral=None, max_iter=1000, # initialize secant method depth = np.interp(outflows[0], self.geometry['discharge'], self.geometry['stage']) - h = (depth * 1.33) + mindepth - h_0 = (depth * 0.67) for i in range(len(inflows) - 1): if not ((inflows[i] > 0.0) or (inflows[i + 1] > 0.0) or (outflows[i] > 0.0)): - outflows.append(0.0) + depth = 0.0 + qdc = 0.0 + outflows.append(qdc) continue Qj_0 = 0.0 iter = 0 @@ -232,6 +232,9 @@ def route_hydrograph_wrf_secant(self, inflows, dt, lateral=None, max_iter=1000, aerror = 0.01 tries = 0 + h = (depth * 1.33) + mindepth + h_0 = (depth * 0.67) + while rerror > 0.01 and aerror >= mindepth and iter <= max_iter: qup = inflows[i] if short_ts: @@ -244,14 +247,17 @@ def route_hydrograph_wrf_secant(self, inflows, dt, lateral=None, max_iter=1000, Ck = np.interp(h_0, self.geometry['stage'], self.geometry['celerity']) Twl = np.interp(h_0, self.geometry['stage'], self.geometry['top_width']) q = np.interp(h_0, self.geometry['stage'], self.geometry['discharge']) - + if Ck > 0.000000: Km = max(dt, dx / Ck) else: Km = dt - X = 0.5*(1-(Qj_0/(2.0*Twl*self.slope*Ck*dx))) - X = min(0.5, max(0.0, X)) + if Ck > 0.0: + X = 0.5*(1-(Qj_0/(2.0*Twl*self.slope*Ck*dx))) + X = min(0.5, max(0.0, X)) + else: + X = 0.5 D = (Km*(1.000 - X) + dt/2.0000) @@ -259,7 +265,7 @@ def route_hydrograph_wrf_secant(self, inflows, dt, lateral=None, max_iter=1000, C2 = (dt/2.0000 - Km*X)/D C3 = (Km*(1.00000000-X)-dt/2.000000)/D C4 = (lateral[i]*dt)/D - Qj_0 = ((C1*quc) + (C2*qup) + (C3*qdp) + C4) - q + Qj_0 = ((C1*qup) + (C2*quc) + (C3*qdp) + C4) - q # Upper Interval Ck = np.interp(h, self.geometry['stage'], self.geometry['celerity']) @@ -283,7 +289,7 @@ def route_hydrograph_wrf_secant(self, inflows, dt, lateral=None, max_iter=1000, C2 = (dt/2.0000 - Km*X)/D C3 = (Km*(1.00000000-X)-dt/2.000000)/D C4 = (lateral[i]*dt)/D - Qj = ((C1*quc) + (C2*qup) + (C3*qdp) + C4) - q + Qj = ((C1*qup) + (C2*quc) + (C3*qdp) + C4) - q # Secant Method if Qj - Qj_0 != 0: @@ -314,10 +320,17 @@ def route_hydrograph_wrf_secant(self, inflows, dt, lateral=None, max_iter=1000, h_0 = h_0 * 0.67 max_iter += 25 - qdc = ((C1*qup)+(C2*quc)+(C3*qdp) + C4) + if ((C1*qup)+(C2*quc)+(C3*qdp) + C4) < 0: + if C4 < 0 and abs(C4) > ((C1*qup)+(C2*quc)+(C3*qdp)): + qdc = 0.0 + else: + qdc = max(((C1*qup)+(C2*quc)+C4), ((C1*qup)+(C3*qdp)+C4)) + else: + qdc = ((C1*qup)+(C2*quc)+(C3*qdp) + C4) if qdc < 0.0: qdc = 0.0 outflows.append(qdc) + depth = h return np.array(outflows) @@ -576,23 +589,58 @@ def generate_geometry(self): for param in ['top_width', 'area', 'wetted_perimeter', 'hydraulic_radius', 'mannings_n', 'discharge', 'celerity']: geom[param] = np.zeros(geom['stage'].shape) - for ind, h_0 in enumerate(geom['stage']): - Ck, WP, WPC, n, nCC, AREA, AREAC, R, So = self.get_geom_at_stage(h_0) - if h_0 > self.bfd: - Twl = self.TwCC - else: - Twl = self.Bw + 2.0*self.z*h_0 - q = ((1/(((WP*n)+(WPC*nCC))/(WP+WPC))) * (AREA+AREAC) * (R**(2./3.)) * np.sqrt(So)) - geom['top_width'][ind] = Twl - geom['area'][ind] = AREA + AREAC - geom['wetted_perimeter'][ind] = WP + WPC - geom['hydraulic_radius'][ind] = R - geom['mannings_n'][ind] = ((n * WP) + (nCC * WPC)) / (WP + WPC) - geom['celerity'][ind] = Ck - geom['discharge'][ind] = q - + mask = geom['stage'] > self.bfd + vec_geom = self.vectorized_get_geom_at_stage(geom['stage']) # Ck, WP, WPC, n, nCC, AREA, AREAC, R, So + geom['celerity'] = vec_geom[0] + geom['wetted_perimeter'] = vec_geom[1] + vec_geom[2] + geom['area'] = vec_geom[5] + vec_geom[6] + geom['hydraulic_radius'] = vec_geom[7] + geom['mannings_n'] = ((vec_geom[3] * vec_geom[1]) + (vec_geom[4] * vec_geom[2])) / (vec_geom[1] + vec_geom[2]) + Twl = self.Bw + 2.0*self.z*geom['stage'] + Twl[mask] = self.TwCC + geom['top_width'] = Twl + geom['discharge'] = ((1/geom['mannings_n']) * geom['area'] * (geom['hydraulic_radius']**(2./3.)) * np.sqrt(self.So)) + np.any(np.isnan(geom['discharge'])) + self.geometry = geom + def vectorized_get_geom_at_stage(self, h_0): + Bw = self.Bw + bfd = self.bfd + z = self.z + TwCC = self.TwCC + So = self.So + n = self.n + nCC = self.nCC + mask = h_0 > bfd + + # get areas + AREA = (Bw + h_0 * z) * h_0 + AREA[mask] = (Bw + bfd * z) * bfd + AREAC = np.zeros_like(AREA) + AREAC[mask] = (TwCC * (h_0 - bfd))[mask] + + # get wetted perimeters + WP = (Bw + 2.0 * h_0 * np.sqrt(1.0 + z*z)) + WP[mask] = (Bw + 2.0 * bfd * np.sqrt(1.0 + z*z)) + WPC = np.zeros_like(WP) + WPC[mask] = (TwCC + (2.0 * (h_0 - bfd)))[mask] + + # get hydraulic radius + R = (AREA + AREAC) / (WP + WPC) + R[np.isnan(R)] = 0 + + # get overbank celerity + Ck = (np.sqrt(So)/n)*((5./3.)*R**(2./3.) - ((2./3.)*R**(5./3.)*(2.0*np.sqrt(1.0 + z*z)/(Bw+2.0*h_0*z)))) + tmp_a = AREA+AREAC + tmp_a[tmp_a < 0.001] = 0.001 + tmp_d = h_0 - bfd + tmp_d[tmp_d < 0] = 0 + Ck[mask] = (((np.sqrt(So)/n)*((5./3.)*R**(2./3.) - ((2./3.)*R**(5./3.)*(2.0*np.sqrt(1.0 + z*z)/(Bw+2.0*bfd*z))))*AREA + ((np.sqrt(So)/(nCC))*(5./3.)*(tmp_d)**(2./3.))*AREAC)/(tmp_a))[mask] + Ck = np.maximum(0.0, Ck) + + return Ck, WP, WPC, n, nCC, AREA, AREAC, R, So + def get_geom_at_stage(self, h_0): Bw = self.Bw bfd = self.bfd From 3e22a5acc29b2b9db78a0a2f176921e93773911d Mon Sep 17 00:00:00 2001 From: sclaw Date: Mon, 29 Jul 2024 11:07:18 -0400 Subject: [PATCH 38/43] proper handling of lakes --- samples/basin_analysis/preprocessing.py | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/samples/basin_analysis/preprocessing.py b/samples/basin_analysis/preprocessing.py index e808dab..2540c1d 100644 --- a/samples/basin_analysis/preprocessing.py +++ b/samples/basin_analysis/preprocessing.py @@ -75,8 +75,10 @@ def extract_from_root(root): return df -run_dir = r'/netfiles/ciroh/floodplainsData/retrospective/basin_run_files' +run_dir = r'/netfiles/ciroh/floodplainsData/retrospective/basin_run_files/MSQ' all_paths = [os.path.join(run_dir, f) for f in os.listdir(run_dir) if f.endswith('.json')] +# all_paths = [r'/netfiles/ciroh/floodplainsData/retrospective/basin_run_files/lamoille_2019_hand.json'] + for paths in all_paths: print(f'Processing {paths}') @@ -98,9 +100,9 @@ def extract_from_root(root): meta_path = r"/netfiles/ciroh/floodplainsData/runs/NWM/network/reach_data.csv" param_path = r"/netfiles/ciroh/floodplainsData/runs/NWM/working/RouteLink_CONUS.nc" da_path = r"/netfiles/ciroh/floodplainsData/runs/NWM/working/drain_areas.csv" - if os.path.exists(out_path): - print(f'{out_path} already exists, skipping') - continue + # if os.path.exists(out_path): + # print(f'{out_path} already exists, skipping') + # continue # Load data ds = pd.read_csv(ds_hydrographs) # treating this as the model domain @@ -152,11 +154,8 @@ def extract_from_root(root): # Add lakes to headwater forcings try: - lake_df = pd.read_csv(lake_path) - tmp_cols = list(lake_df.columns) - tmp_cols.remove('Unnamed: 0') - tmp_cols.remove('datetime') - lakes = [float(c) for c in tmp_cols] + lakes = meta_df['NHDWaterbodyComID'].unique() + lakes = lakes[lakes != -9999] force_lakes = list() for l in lakes: # Find root of each lake From cb3c509c1f35dac592d8dea2fbe54c3e2d39eea8 Mon Sep 17 00:00:00 2001 From: sclaw Date: Mon, 29 Jul 2024 16:33:07 -0400 Subject: [PATCH 39/43] refactor to put all solvers under the same method --- src/muskingumcunge/controller.py | 10 +- src/muskingumcunge/network.py | 6 +- src/muskingumcunge/reach.py | 258 +++++++++---------------------- src/muskingumcunge/route.pyx | 192 +++++++++++++++++++++-- 4 files changed, 260 insertions(+), 206 deletions(-) diff --git a/src/muskingumcunge/controller.py b/src/muskingumcunge/controller.py index 79505e7..ed35f95 100644 --- a/src/muskingumcunge/controller.py +++ b/src/muskingumcunge/controller.py @@ -14,11 +14,13 @@ def create_project(path): "lateral_path": os.path.join(os.path.dirname(path), "laterals.csv"), "head_path": os.path.join(os.path.dirname(path), "headwaters.csv"), "lake_path": os.path.join(os.path.dirname(path), "lake.csv"), - "geometry_dir": "/users/k/l/klawson1/netfiles/ciroh/floodplainsData/runs/NWM/geometry", + "geometry_dir": "", # expecting the format of the geometry folder generated by https://github.com/CIROH-UVM/floodplain-routing "geometry_source": "NWM", - "optimize_dx": True, + "optimize_dx": False, "conserve_mass": False, - "lat_addition": "middle", + "short_ts": False, + "dt": "5m", + "lat_addition": "top", "id_field": "comid", "to_id_field": "ds_comid", "mode": 'basin', @@ -214,7 +216,7 @@ def execute(meta_path): network.load_initial_conditions(init_outflows) # kind of unnecessary # Run model - network.run_event(optimize_dx=paths['optimize_dx'], conserve_mass=paths['conserve_mass'], lat_addition=paths['lat_addition']) + network.run_event(optimize_dx=paths['optimize_dx'], conserve_mass=paths['conserve_mass'], lat_addition=paths['lat_addition'], short_ts=paths['short_ts']) os.makedirs(paths['out_dir'], exist_ok=True) out_path = os.path.join(paths['out_dir'], "model_outflows.csv") network.out_df.to_csv(out_path) diff --git a/src/muskingumcunge/network.py b/src/muskingumcunge/network.py index e8d5b20..ce08449 100644 --- a/src/muskingumcunge/network.py +++ b/src/muskingumcunge/network.py @@ -72,10 +72,10 @@ def load_initial_conditions(self, initial_conditions): # loads initial outflows for each reach self.init_outflows = initial_conditions - def run_event(self, optimize_dx=True, conserve_mass=False, lat_addition='middle'): + def run_event(self, optimize_dx=True, conserve_mass=False, lat_addition='middle', short_ts=False): # calculate total inflow t_start = time.perf_counter() - dt = (self.forcing_df.index[1] - self.forcing_df.index[0]).seconds / 3600 + dt = (self.forcing_df.index[1] - self.forcing_df.index[0]).seconds counter = 0 pct5 = int(len(self.post_order) / 20) + 1 for node in self.post_order: @@ -113,7 +113,7 @@ def run_event(self, optimize_dx=True, conserve_mass=False, lat_addition='middle' init_out = self.init_outflows[node] else: init_out = None - routed = reach.route_hydrograph_wrf_secant(routed, dt, lateral=l, initial_outflow=init_out) + routed = reach.route_hydrograph(routed, dt, lateral=l, initial_outflow=init_out, short_ts=short_ts, solver='wrf-c') except AssertionError as e: print(f"Error routing {node}: {e}") diff --git a/src/muskingumcunge/reach.py b/src/muskingumcunge/reach.py index 02e5945..ff6ad78 100644 --- a/src/muskingumcunge/reach.py +++ b/src/muskingumcunge/reach.py @@ -1,19 +1,8 @@ import numpy as np from scipy.ndimage import gaussian_filter1d -from .route import croute +from .route import croute_fread, croute_wrf class BaseReach: - geometry = {'stage': None, - 'top_width': None, - 'area': None, - 'wetted_perimeter': None, - 'hydraulic_radius': None, - 'mannings_n': None} - rating_curve = {'stage': None, - 'discharge': None} - muskingum_params = {'stage': None, - 'k': None, - 'x': None} def __init__(self, width, mannings_n, slope, reach_length, max_stage=10, stage_resolution=50): self.width = width @@ -48,141 +37,22 @@ def calculate_x_ref(self, inflows, dt): c_ref = np.interp(np.log(q_ref), self.geometry['log_q'], self.geometry['celerity']) x_ref = 0.5 * ((c_ref * dt * 60 * 60) + (q_ref / (b_ref * self.slope * c_ref))) return x_ref - - def route_hydrograph_c(self, inflows, dt): - reach_length = self.reach_length - slope = self.slope - geometry = self.geometry - - return croute(inflows, dt, reach_length, slope, geometry) - - def route_hydrograph(self, inflows, dt, lateral=None, max_iter=1000, initial_outflow=None): - dt = dt * 60 * 60 # MUST DELETE AFTER DEBUGGING - outflows = list() + + def route_hydrograph(self, inflows, dt, lateral=None, initial_outflow=None, short_ts=False, solver='wrf-c'): + assert np.all(inflows >= 0), 'Negative inflows detected' + assert np.all(~np.isnan(inflows)), 'NaN inflows detected' + # initialize outflow array + outflows = np.zeros(inflows.shape) if initial_outflow is not None: - outflows.append(initial_outflow) + outflows[0] = initial_outflow else: - outflows.append(inflows[0]) - assert max(inflows) < max(self.geometry['discharge']), 'Rating Curve does not cover range of flowrates in hydrograph' + outflows[0] = inflows[0] + # initialize lateraals if none if lateral is None: lateral = np.zeros_like(inflows) - lateral = (lateral[:-1] + lateral[1:]) / 2 - lateral = np.append(lateral, 0) - for i in range(len(inflows) - 1): - # q_guess = sum([inflows[i], inflows[i + 1], outflows[i], lateral[i]]) / 3 - q_guess = sum([inflows[i], inflows[i], outflows[i], lateral[i]]) / 3 - last_guess = q_guess * 2 - counter = 1 - while abs(last_guess - q_guess) > 0.003: # from handbook of hydrology page 328 - counter += 1 - last_guess = q_guess.copy() - # reach_q = sum([inflows[i], inflows[i + 1], outflows[i], q_guess]) / 4 - reach_q = sum([inflows[i], inflows[i], outflows[i], q_guess]) / 4 - - # Interpolate - log_reach_q = np.log(reach_q) - b_tmp = np.exp(np.interp(log_reach_q, self.geometry['log_q'], self.geometry['log_width'])) - c_tmp = np.interp(log_reach_q, self.geometry['log_q'], self.geometry['celerity']) - - # courant = c_tmp * dt * 60 * 60 / self.reach_length - # reynold = reach_q / (self.slope * c_tmp * self.reach_length * b_tmp) - - # c0 = (-1 + courant + reynold) / (1 + courant + reynold) - # c1 = (1 + courant - reynold) / (1 + courant + reynold) - # c2 = (1 - courant + reynold) / (1 + courant + reynold) - # c3 = (2 * courant) / (1 + courant + reynold) - - # q_guess = (c0 * inflows[i + 1]) + (c1 * inflows[i]) + (c2 * outflows[i]) + (c3 * lateral[i]) - # q_guess = max(min(inflows), q_guess) - - # WRF way - if c_tmp > 0: - Km = self.reach_length / c_tmp - Km = max(dt, Km) - else: - Km = dt - X = 0.5*(1-(reach_q/(2.0*b_tmp*self.slope*c_tmp*self.reach_length))) - X = min(0.5, max(0.0, X)) - D = (Km*(1.000 - X) + (dt)/2.0000) - - C1 = (Km*X + dt/2.000000)/D - C2 = (dt/2.0000 - Km*X)/D - C3 = (Km*(1.00000000-X)-dt/2.000000)/D - C4 = (lateral[i]*dt)/D - # q_guess = ((C1*inflows[i]) + (C2*inflows[i + 1]) + (C3*outflows[i])) - q_guess = ((C1*inflows[i]) + (C2*inflows[i]) + (C3*outflows[i])) - - if counter == max_iter: - last_guess = q_guess - outflows.append(q_guess) - - - return np.array(outflows) - - def route_hydrograph_mct(self, inflows, dt, max_iter=1000): - outflows = list() - outflows.append(inflows[0]) - assert max(inflows) < max(self.geometry['discharge']), 'Rating Curve does not cover range of flowrates in hydrograph' - - # warm up parameters - qref1 = (inflows[0] + outflows[0]) / 2 - aref1 = np.exp(np.interp(np.log(qref1), self.geometry['log_q'], self.geometry['log_area'])) - bref1 = np.exp(np.interp(np.log(qref1), self.geometry['log_q'], self.geometry['log_width'])) - cref1 = np.interp(np.log(qref1), self.geometry['log_q'], self.geometry['celerity']) - betaref1 = (cref1 * aref1) / qref1 - courantref1 = (cref1 / betaref1) * (dt * 60 * 60 / self.reach_length) - reynoldref1 = qref1 / (betaref1 * bref1 * self.slope * cref1 * self.reach_length) - for i in range(len(inflows) - 1): - q_guess = outflows[i] + (inflows[i + 1] - inflows[i]) - last_guess = q_guess * 2 - counter = 1 - while abs(last_guess - q_guess) > 0.003: # from handbook of hydrology page 328 - counter += 1 - last_guess = q_guess.copy() - - qref2 = (inflows[i + 1] + q_guess) / 2 - # interpolate - aref2 = np.exp(np.interp(np.log(qref2), self.geometry['log_q'], self.geometry['log_area'])) - bref2 = np.exp(np.interp(np.log(qref2), self.geometry['log_q'], self.geometry['log_width'])) - cref2 = np.interp(np.log(qref2), self.geometry['log_q'], self.geometry['celerity']) - betaref2 = (cref2 * aref2) / qref2 - courantref2 = (cref2 / betaref2) * (dt * 60 * 60 / self.reach_length) - reynoldref2 = qref2 / (betaref2 * bref2 * self.slope * cref2 * self.reach_length) - - # MCT parameters - c0 = (-1 + courantref1 + reynoldref1) / (1 + courantref2 + reynoldref2) - c1 = ((1 + courantref1 - reynoldref1) / (1 + courantref2 + reynoldref2)) * (courantref2 / courantref1) - c2 = ((1 - courantref1 + reynoldref1) / (1 + courantref2 + reynoldref2)) * (courantref2 / courantref1) - - # Estimate outflow - q_guess = (c0 * inflows[i + 1]) + (c1 * inflows[i]) + (c2 * outflows[i]) - q_guess = max(min(inflows), q_guess) - if counter == max_iter: - last_guess = q_guess - outflows.append(q_guess) - if (cref2 / betaref2) > 1: - print((cref2 / betaref2)) - qref1 = qref2 - aref1 = aref2 - bref1 = bref2 - cref1 = bref2 - betaref1 = betaref2 - courantref1 = courantref2 - reynoldref1 = reynoldref2 - - return np.array(outflows) - - def route_hydrograph_wrf_secant(self, inflows, dt, lateral=None, max_iter=1000, initial_outflow=None): - # Secant method solver from WRF-HYDRO - dt = dt * 60 * 60 # MUST DELETE AFTER DEBUGGING - dx = self.reach_length - outflows = list() - if initial_outflow is not None: - outflows.append(initial_outflow) - else: - outflows.append(inflows[0]) + # check if rating curve covers range of inflows. If not, add more points if max(inflows) > max(self.geometry['discharge']): add_res = 200 add_stages = np.linspace(0, 9 * self.geometry['stage'][-1], add_res) @@ -204,27 +74,74 @@ def route_hydrograph_wrf_secant(self, inflows, dt, lateral=None, max_iter=1000, self.geometry['mannings_n'] = np.append(self.geometry['mannings_n'], add_n) self.geometry['discharge'] = np.append(self.geometry['discharge'], add_q) self.geometry['celerity'] = np.append(self.geometry['celerity'], add_c) + self.geometry['log_q'] = np.log(self.geometry['discharge']) + self.geometry['log_width'] = np.log(self.geometry['top_width']) assert max(inflows) < max(self.geometry['discharge']), 'Rating Curve does not cover range of flowrates in hydrograph' - if lateral is None: - lateral = np.zeros_like(inflows) - lateral = (lateral[:-1] + lateral[1:]) / 2 - lateral = np.append(lateral, 0) + # initialize secant method + depth = np.interp(outflows[0], self.geometry['discharge'], self.geometry['stage']) + + if solver == 'wrf': + outflows = self.route_hydrograph_wrf(inflows, outflows, lateral, depth, self.reach_length, dt, short_ts) + elif solver == 'wrf-c': + outflows = croute_wrf(inflows, outflows, lateral, depth, self.geometry['stage'], self.geometry['celerity'], self.geometry['top_width'], self.geometry['discharge'], self.slope, self.reach_length, dt, short_ts) + elif solver == 'fread': + outflows = self.route_hydrograph_fread(inflows, outflows, lateral, dt, short_ts) + elif solver == 'fread-c': + inflows = inflows + lateral # did not program middle addition for fread + outflows = croute_fread(inflows, outflows, dt, self.reach_length, self.slope, self.geometry, short_ts) + + return np.array(outflows) + + def route_hydrograph_fread(self, inflows, outflows, lateral, dt, short_ts): + max_iter = 1000 + for i in range(len(inflows) - 1): + if short_ts: + us_cur = np.copy(inflows[i]) + else: + us_cur = np.copy(inflows[i + 1]) + q_guess = sum([inflows[i], us_cur, outflows[i], lateral[i]]) / 3 + last_guess = q_guess * 2 + counter = 1 + while abs(last_guess - q_guess) > 0.003: # from handbook of hydrology page 328 + counter += 1 + last_guess = q_guess.copy() + reach_q = sum([inflows[i], us_cur, outflows[i], q_guess]) / 4 + + # Interpolate + log_reach_q = np.log(reach_q) + b_tmp = np.exp(np.interp(log_reach_q, self.geometry['log_q'], self.geometry['log_width'])) + c_tmp = np.interp(log_reach_q, self.geometry['log_q'], self.geometry['celerity']) + + courant = c_tmp * dt/ self.reach_length + reynold = reach_q / (self.slope * c_tmp * self.reach_length * b_tmp) + + c0 = (-1 + courant + reynold) / (1 + courant + reynold) + c1 = (1 + courant - reynold) / (1 + courant + reynold) + c2 = (1 - courant + reynold) / (1 + courant + reynold) + c3 = (2 * courant) / (1 + courant + reynold) + + q_guess = (c0 * us_cur) + (c1 * inflows[i]) + (c2 * outflows[i]) + (c3 * lateral[i]) + q_guess = max(min(inflows), q_guess) + + if counter == max_iter: + last_guess = q_guess + outflows[i + 1] = q_guess + return np.array(outflows) + + + def route_hydrograph_wrf(self, inflows, outflows, lateral, depth, dx, dt, short_ts): # local variables mindepth = 0.01 max_iter = 100 - short_ts = True - - # initialize secant method - depth = np.interp(outflows[0], self.geometry['discharge'], self.geometry['stage']) for i in range(len(inflows) - 1): if not ((inflows[i] > 0.0) or (inflows[i + 1] > 0.0) or (outflows[i] > 0.0)): depth = 0.0 qdc = 0.0 - outflows.append(qdc) + outflows[i + 1] = qdc continue Qj_0 = 0.0 iter = 0 @@ -329,23 +246,24 @@ def route_hydrograph_wrf_secant(self, inflows, dt, lateral=None, max_iter=1000, qdc = ((C1*qup)+(C2*quc)+(C3*qdp) + C4) if qdc < 0.0: qdc = 0.0 - outflows.append(qdc) + outflows[i + 1] = qdc depth = h - return np.array(outflows) + return outflows + def optimize_route_params(self, inflows, dt): # Ponce method qref = (inflows.max() + inflows.min()) / 2 cref = np.interp(qref, self.geometry['discharge'], self.geometry['celerity']) twref = np.interp(qref, self.geometry['discharge'], self.geometry['top_width']) - dxc = dt * 60 * 60 * cref # Courant length + dxc = dt * cref # Courant length dxd = (qref / twref) / (self.slope * cref) # characteristic reach length dxmax = 0.5 * (dxc + dxd) peak_loc = np.argmax(self.geometry['discharge'] > inflows.max()) peak_loc = max([peak_loc, 1]) cmax = self.geometry['celerity'][:peak_loc].max() # I think the fact that this is the max should cover us for making all stable - dxmin = cmax * (dt * 60 * 60) + dxmin = cmax * (dt) dx = max([dxmin, dxmax]) subreaches = int(np.ceil(self.reach_length / dx)) dx = self.reach_length / subreaches @@ -427,12 +345,6 @@ def generate_geometry(self): dq_da = dq / da geom['celerity'] = np.append(dq_da, dq_da[-1]) - # dpdy = (geom['wetted_perimeter'][1:] - geom['wetted_perimeter'][:-1]) / (geom['stage'][1:] - geom['stage'][:-1]) - # dpdy = np.append(dpdy, dpdy[-1]) - # dpdy[dpdy < 0.0001] = 0.0001 - # k_prime = (5 / 3) - ((2 / 3) * (geom['area'] / (geom['top_width'] * geom['wetted_perimeter'])) * dpdy) - # geom['celerity'] = k_prime * (geom['discharge'] / geom['area']) - geom['celerity'][geom['celerity'] < 0.0001] = 0.0001 geom['log_q'] = np.log(geom['discharge']) @@ -460,22 +372,10 @@ def generate_muskingum_params(self): c[0] = c[1] self.muskingum_params['k'] = self.reach_length / c - self.muskingum_params['k'] /= (60 * 60) # Convert seconds to hours self.muskingum_params['x'] = (1 / 2) - (self.rating_curve['discharge'] / (2 * c * self.width * self.slope * self.reach_length)) self.muskingum_params['x'][self.muskingum_params['x'] < 0] = 0 class CustomReach(BaseReach): - geometry = {'stage': None, - 'top_width': None, - 'area': None, - 'wetted_perimeter': None, - 'hydraulic_radius': None, - 'mannings_n': None} - rating_curve = {'stage': None, - 'discharge': None} - muskingum_params = {'stage': None, - 'k': None, - 'x': None} def __init__(self, mannings_n, slope, reach_length, stages, top_widths, areas, perimeters): self.mannings_n = mannings_n @@ -522,7 +422,6 @@ def generate_geometry(self, stages, top_widths, areas, perimeters): geom['celerity'] = np.nan_to_num(geom['celerity']) ## TODO: geom['celerity'][<0] = 0 self.geometry = geom - def clean_looped_rating_curve(self, discharge): # Enforce monotonic discharge inreases @@ -600,6 +499,8 @@ def generate_geometry(self): Twl[mask] = self.TwCC geom['top_width'] = Twl geom['discharge'] = ((1/geom['mannings_n']) * geom['area'] * (geom['hydraulic_radius']**(2./3.)) * np.sqrt(self.So)) + geom['log_q'] = np.log(geom['discharge']) + geom['log_width'] = np.log(geom['top_width']) np.any(np.isnan(geom['discharge'])) self.geometry = geom @@ -679,17 +580,6 @@ def get_geom_at_stage(self, h_0): class SigmoidalReach(BaseReach): - geometry = {'stage': None, - 'top_width': None, - 'area': None, - 'wetted_perimeter': None, - 'hydraulic_radius': None, - 'mannings_n': None} - rating_curve = {'stage': None, - 'discharge': None} - muskingum_params = {'stage': None, - 'k': None, - 'x': None} def __init__(self, mannings_n, slope, reach_length, ch_w, fp_w, bkf_el, fp_s, max_stage=10, stage_resolution=50): self.mannings_n = mannings_n diff --git a/src/muskingumcunge/route.pyx b/src/muskingumcunge/route.pyx index 57699a9..9c1998b 100644 --- a/src/muskingumcunge/route.pyx +++ b/src/muskingumcunge/route.pyx @@ -2,14 +2,12 @@ import numpy as np import cython from libc.math cimport log, exp -def croute(py_inflows, dt, reach_length, slope, geometry, max_iter=1000): - inflows: cython.double[:] = py_inflows - outflows: cython.double[:] = np.empty_like(py_inflows) +def croute_fread(inflows: cython.double[:], outflows: cython.double[:], dt: cython.float, reach_length: cython.float, slope: cython.float, geometry, short_ts: cython.bint): geom_log_q: cython.double[:] = geometry['log_q'] geom_log_w: cython.double[:] = geometry['log_width'] geom_q: cython.double[:] = geometry['discharge'] geom_c: cython.double[:] = geometry['celerity'] - series_length: cython.int = len(py_inflows) + series_length: cython.int = len(inflows) i: cython.int q_guess: cython.float last_guess: cython.float @@ -23,41 +21,205 @@ def croute(py_inflows, dt, reach_length, slope, geometry, max_iter=1000): c0: cython.float c1: cython.float c2: cython.float - min_flow: cython.float = py_inflows.min() - max_iter: cython.int = max_iter - reach_length: cython.float = reach_length - dt: cython.float = dt - slope: cython.float = slope + min_flow: cython.float = np.min(inflows) + max_iter: cython.int = 1000 assert np.log(max(inflows)) < max(geom_log_q), 'Rating Curve does not cover range of flowrates in hydrograph' - - outflows[0] = inflows[0] for i in range(series_length - 1): - q_guess = (inflows[i] + inflows[i + 1] + outflows[i]) / 3 + if short_ts: + us_cur = inflows[i] + else: + us_cur = inflows[i + 1] + q_guess = (inflows[i] + us_cur + outflows[i]) / 3 last_guess = q_guess * 2 counter = 1 while abs(last_guess - q_guess) > 0.003: counter += 1 last_guess = q_guess - reach_q = (inflows[i] + inflows[i + 1] + outflows[i] + q_guess) / 4 + reach_q = (inflows[i] + us_cur + outflows[i] + q_guess) / 4 # Interpolate log_reach_q = log(reach_q) b_tmp = exp(np.interp(log_reach_q, geom_log_q, geom_log_w)) c_tmp = np.interp(log_reach_q, geom_log_q, geom_c) - courant = c_tmp * dt * 60 * 60 / reach_length + courant = c_tmp * dt / reach_length reynold = reach_q / (slope * c_tmp * reach_length * b_tmp) c0 = (-1 + courant + reynold) / (1 + courant + reynold) c1 = (1 + courant - reynold) / (1 + courant + reynold) c2 = (1 - courant + reynold) / (1 + courant + reynold) - q_guess = (c0 * inflows[i + 1]) + (c1 * inflows[i]) + (c2 * outflows[i]) + q_guess = (c0 * us_cur) + (c1 * inflows[i]) + (c2 * outflows[i]) q_guess = max(min_flow, q_guess) if counter == max_iter: last_guess = q_guess outflows[i + 1] = q_guess + return np.array(outflows) + +def interpolate(x: cython.double, xp: cython.double[:], yp: cython.double[:]): + i: cython.int + for i in range(len(xp) - 1): + if xp[i] <= x <= xp[i + 1]: + return yp[i] + ((x - xp[i]) * ((yp[i + 1] - yp[i]) / (xp[i + 1] - xp[i]))) + return 0.0 + +def croute_wrf( + inflows: cython.double[:], + outflows: cython.double[:], + lateral: cython.double[:], + depth: cython.double, + stages: cython.double[:], + cks: cython.double[:], + twls: cython.double[:], + qs: cython.double[:], + slope: cython.double, + dx: cython.double, + dt: cython.double, + short_ts: cython.bint): + + max_iter: cython.int = 100 + mindepth: cython.double = 0.01 + i: cython.int + + Qj_0: cython.double + Qj: cython.double + iters: cython.int + rerror: cython.double + aerror: cython.double + tries: cython.int + + h: cython.double + h_0: cython.double + qup: cython.double = 0 + quc: cython.double = 0 + qdp: cython.double = 0 + qdc: cython.double = 0 + depth: cython.double + Ck: cython.double + Twl: cython.double + q: cython.double + Km: cython.double + X: cython.double + D: cython.double + C1: cython.double = 0 + C2: cython.double = 0 + C3: cython.double = 0 + C4: cython.double = 0 + h_1: cython.double + + for i in range(len(inflows) - 1): + if not ((inflows[i] > 0.0) or (inflows[i + 1] > 0.0) or (outflows[i] > 0.0)): + depth = 0.0 + qdc = 0.0 + outflows[i + 1] = qdc + continue + Qj_0 = 0.0 + iters = 0 + rerror = 1.0 + aerror = 0.01 + tries = 0 + + h = (depth * 1.33) + mindepth + h_0 = (depth * 0.67) + + while rerror > 0.01 and aerror >= mindepth and iters <= max_iter: + qup = inflows[i] + if short_ts: + quc = inflows[i] + else: + quc = inflows[i + 1] + qdp = outflows[i] + + # Lower Interval + Ck = interpolate(h_0, stages, cks) + Twl = interpolate(h_0, stages, twls) + q = interpolate(h_0, stages, qs) + + if Ck > 0.000000: + Km = max(dt, dx / Ck) + else: + Km = dt + + if Ck > 0.0: + X = 0.5*(1-(Qj_0/(2.0*Twl*slope*Ck*dx))) + X = min(0.5, max(0.0, X)) + else: + X = 0.5 + + D = (Km*(1.000 - X) + dt/2.0000) + + C1 = (Km*X + dt/2.000000)/D + C2 = (dt/2.0000 - Km*X)/D + C3 = (Km*(1.00000000-X)-dt/2.000000)/D + C4 = (lateral[i]*dt)/D + Qj_0 = ((C1*qup) + (C2*quc) + (C3*qdp) + C4) - q + + # Upper Interval + Ck = interpolate(h, stages, cks) + Twl = interpolate(h, stages, twls) + q = interpolate(h, stages, qs) + + if Ck > 0.000000: + Km = max(dt, dx / Ck) + else: + Km = dt + + if Ck > 0.0: + X = 0.5*(1-(((C1*qup)+(C2*quc)+(C3*qdp) + C4)/(2.0*Twl*slope*Ck*dx))) + X = min(0.5, max(0.25, X)) + else: + X = 0.5 + + D = (Km*(1.000 - X) + dt/2.0000) + + C1 = (Km*X + dt/2.000000)/D + C2 = (dt/2.0000 - Km*X)/D + C3 = (Km*(1.00000000-X)-dt/2.000000)/D + C4 = (lateral[i]*dt)/D + Qj = ((C1*qup) + (C2*quc) + (C3*qdp) + C4) - q + + # Secant Method + if Qj - Qj_0 != 0: + h_1 = h - ((Qj * (h_0 - h))/(Qj_0 - Qj)) + if h_1 < 0.0: + h_1 = h + else: + h_1 = h + + if h > 0.0: + rerror = abs((h_1 - h) / h) + aerror = abs(h_1 - h) + else: + rerror = 0.0 + aerror = 0.9 + + h_0 = max(0.0, h) + h = max(0.0, h_1) + iters += 1 + + if h < mindepth: + break + + if iters >= max_iter: + tries += 1 + if tries <= 4: + h = h * 1.33 + h_0 = h_0 * 0.67 + max_iter += 25 + + if ((C1*qup)+(C2*quc)+(C3*qdp) + C4) < 0: + if C4 < 0 and abs(C4) > ((C1*qup)+(C2*quc)+(C3*qdp)): + qdc = 0.0 + else: + qdc = max(((C1*qup)+(C2*quc)+C4), ((C1*qup)+(C3*qdp)+C4)) + else: + qdc = ((C1*qup)+(C2*quc)+(C3*qdp) + C4) + if qdc < 0.0: + qdc = 0.0 + outflows[i + 1] = qdc + depth = h + return np.array(outflows) \ No newline at end of file From 5c2936955d8253ba298904631ec610ebb99430e4 Mon Sep 17 00:00:00 2001 From: sclaw Date: Mon, 29 Jul 2024 16:34:59 -0400 Subject: [PATCH 40/43] update time format --- samples/simple_routing/simple_routing.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/samples/simple_routing/simple_routing.py b/samples/simple_routing/simple_routing.py index d0b22e7..b5cefe1 100644 --- a/samples/simple_routing/simple_routing.py +++ b/samples/simple_routing/simple_routing.py @@ -16,7 +16,7 @@ reach = TrapezoidalReach(100, 100, 0.1, 0.001, 2000) # Route flows -outflows = reach.route_hydrograph(inflows, 0.1) +outflows = reach.route_hydrograph(inflows, 360, solver='fread') # Plot results fig, ax = plt.subplots() @@ -26,4 +26,5 @@ ax.set_ylabel('Flowrate(cms)') plt.legend() plt.tight_layout() -plt.show() \ No newline at end of file +plt.savefig('simple_routing.png', dpi=300) +# plt.show() \ No newline at end of file From b59cf2419d307b3fe84785c6d7c9f97774ea4b74 Mon Sep 17 00:00:00 2001 From: sclaw Date: Mon, 29 Jul 2024 16:36:14 -0400 Subject: [PATCH 41/43] update mannings definition. change solver --- samples/CIROH/muskingum_cunge.py | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/samples/CIROH/muskingum_cunge.py b/samples/CIROH/muskingum_cunge.py index 251c9e1..32b0d52 100644 --- a/samples/CIROH/muskingum_cunge.py +++ b/samples/CIROH/muskingum_cunge.py @@ -7,6 +7,7 @@ from scipy.ndimage import gaussian_filter1d from statsmodels.nonparametric.smoothers_lowess import lowess import time +import sys ### Static Data ### @@ -75,6 +76,7 @@ def execute(meta_path, debug_plots=False): reach_data = pd.read_csv(run_dict['reach_meta_path']) reach_data[run_dict['id_field']] = reach_data[run_dict['id_field']].astype(np.int64).astype(str) reach_data = reach_data.set_index(run_dict['id_field']) + mannings_n = 0.095 # Setup Hydrographs hydrographs = ['Q2_Short', 'Q2_Medium', 'Q2_Long', 'Q10_Short', 'Q10_Medium', 'Q10_Long', 'Q50_Short', 'Q50_Medium', 'Q50_Long', 'Q100_Short', 'Q100_Medium', 'Q100_Long'] @@ -97,7 +99,7 @@ def execute(meta_path, debug_plots=False): counter = 1 t_start = time.perf_counter() reaches = reach_data.index.to_list() - for reach in reaches[625:]: + for reach in reaches: print(f'{counter} / {len(reaches)} | {round((len(reaches) - counter) * ((time.perf_counter() - t_start) / counter), 1)} seconds left | {(len(reaches) - counter)} left | | {((time.perf_counter() - t_start) / counter)} rate') counter += 1 @@ -136,7 +138,7 @@ def execute(meta_path, debug_plots=False): continue # Create reach - mc_reach = CustomReach(0.035, slope, 1000, tmp_geom['el'], tmp_geom['area'], tmp_geom['vol'], tmp_geom['p']) + mc_reach = CustomReach(mannings_n, slope, 1000, tmp_geom['el'], tmp_geom['area'], tmp_geom['vol'], tmp_geom['p']) # Route hydrographs results_dict[run_dict['id_field']].append(reach) @@ -234,7 +236,7 @@ def execute(meta_path, debug_plots=False): outflows = inflows.copy() for iter in range(subreaches): try: - outflows = mc_reach.route_hydrograph_c(outflows, dt) + outflows = mc_reach.route_hydrograph(outflows, (dt * 60 * 60), lateral=None, initial_outflow=None, short_ts=False, solver='fread-c') except AssertionError: outflows = np.repeat(np.nan, outflows.shape[0]) break @@ -244,7 +246,7 @@ def execute(meta_path, debug_plots=False): attenuation_per_km = raw_attenuation / (dx * subreaches / 1000) pct_attenuation = raw_attenuation / inflows.max() pct_attenuation_km = (1 - ((1 - pct_attenuation) ** (1000 / (dx * subreaches)))) - tmp_diff_number = ((9 * np.pi) / 50) * (((0.035 ** (6 / 5)) * (max(inflows) ** (1 / 5))) / (((slope * 100) ** (8 / 5)) * (dt * 20))) + tmp_diff_number = ((9 * np.pi) / 50) * (((mannings_n ** (6 / 5)) * (max(inflows) ** (1 / 5))) / (((slope * 100) ** (8 / 5)) * (dt * 20))) peak_loc = np.argmax(outflows) dqs = outflows[1:] - outflows[:-1] rising_dqs = np.abs(dqs[:peak_loc]) @@ -303,8 +305,8 @@ def execute(meta_path, debug_plots=False): out_data = pd.DataFrame(results_dict) out_data = out_data.set_index(run_dict['id_field']) os.makedirs(os.path.dirname(run_dict['muskingum_path']), exist_ok=True) - # out_data.to_csv(run_dict['muskingum_path']) + out_data.to_csv(run_dict['muskingum_path']) if __name__ == '__main__': - run_path = r"/users/k/l/klawson1/netfiles/ciroh/floodplainsData/runs/8/run_metadata.json" - execute(run_path, debug_plots=True) + run_path = sys.argv[1] + execute(run_path, debug_plots=False) From 48246ffd64cd4c02089cdb49ddf981602b49942d Mon Sep 17 00:00:00 2001 From: sclaw Date: Mon, 29 Jul 2024 16:37:27 -0400 Subject: [PATCH 42/43] update config format --- samples/basin_analysis/make_runs.py | 1 + 1 file changed, 1 insertion(+) diff --git a/samples/basin_analysis/make_runs.py b/samples/basin_analysis/make_runs.py index 4f36860..f644065 100644 --- a/samples/basin_analysis/make_runs.py +++ b/samples/basin_analysis/make_runs.py @@ -14,6 +14,7 @@ def create_project(data_path, run_path, geom_type): "dt": "5m", "optimize_dx": False, "conserve_mass": False, + "short_ts": True, "lat_addition": "top", "id_field": "comid", "to_id_field": "ds_comid", From 8bfdb8be12012b58566ce3589e4207d3fe7f7c13 Mon Sep 17 00:00:00 2001 From: sclaw Date: Mon, 29 Jul 2024 16:37:42 -0400 Subject: [PATCH 43/43] catchup --- samples/basin_analysis/run_all_basinwide.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/samples/basin_analysis/run_all_basinwide.py b/samples/basin_analysis/run_all_basinwide.py index 12b4c91..045dd12 100644 --- a/samples/basin_analysis/run_all_basinwide.py +++ b/samples/basin_analysis/run_all_basinwide.py @@ -1,8 +1,13 @@ import os +import json from muskingumcunge.controller import handler +import numpy as np run_dir = r'/netfiles/ciroh/floodplainsData/retrospective/basin_run_files' all_paths = [os.path.join(run_dir, f) for f in os.listdir(run_dir) if f.endswith('.json')] for path in all_paths: + with open(path, 'r') as f: + dict = json.load(f) + handler(path) \ No newline at end of file