diff --git a/src/troute-network/troute/HYFeaturesNetwork.py b/src/troute-network/troute/HYFeaturesNetwork.py index 1883e3ba5..16edac6aa 100644 --- a/src/troute-network/troute/HYFeaturesNetwork.py +++ b/src/troute-network/troute/HYFeaturesNetwork.py @@ -44,7 +44,7 @@ def node_key_func(x): df = df.loc[lake_id_mask] return df -def read_qlats(forcing_parameters, segment_index, nexus_to_downstream_flowpath_dict): +def read_qlats(forcing_parameters, segment_index, nexus_to_downstream_flowpath_dict, terminal_links): # STEP 5: Read (or set) QLateral Inputs if __showtiming__: start_time = time.time() @@ -92,7 +92,32 @@ def read_qlats(forcing_parameters, segment_index, nexus_to_downstream_flowpath_d qlat_df = pd.concat( (nexuses_flows_df.loc[int(k)].rename(index={int(k):v}) for k,v in nexus_to_downstream_flowpath_dict.items() ), axis=1 ).T - + + #For a terminal nexus, we want to include the lateral flow from the catchment contributing to that nexus + #one way to do that is to cheat and put that lateral flow at the upstream...this is probably the simplest way + #right now. The other is to create a virtual channel segment downstream to "route" i.e accumulate into + #but it isn't clear right now how to do that with flow/velocity/depth requirements + #find the terminal nodes + for tnx, up in terminal_links.items(): + print(up, tnx) + #print(nexuses_flows_df) + #first need to ensure there is an upstream location to dump to + for nex in up: + try: + #FIXME if multiple upstreams exist in this case then a choice is to be made as to which it goes into + #some cases the choice is easy cause the upstream doesn't exist, but in others, it may not be so simple + #in such cases where multiple valid upstream nexuses exist, perhaps the mainstem should be used? + qlat_df.loc[up] += nexuses_flows_df.loc[tnx] + break #flow added, don't add it again! + except KeyError: + #this upstream doesn't actually exist on the network (maybe it is a headwater?) + #or perhaps the output file doesnt exist? If this is the case, this isn't a good trap + #but for now, add the flow to a known good nexus upstream of the terminal + continue + #TODO what happens if can't put the qlat anywhere? Right now this silently ignores the issue... + qlat_df.drop(tnx, inplace=True) + #print(nexuses_flows_df) + # The segment_index has the full network set of segments/flowpaths. # Whereas the set of flowpaths that are downstream of nexuses is a # subset of the segment_index. Therefore, all of the segments/flowpaths @@ -333,8 +358,13 @@ def make_list(s): print("channel initial states complete") if __showtiming__: print("... in %s seconds." % (time.time() - start_time)) - - self._qlateral = read_qlats(forcing_parameters, self._dataframe.index, self.downstream_flowpath_dict) + terminal = self._dataframe[ + ~self._dataframe["downstream"].isin(self._dataframe.index) + ]["downstream"] + upstream_terminal = dict() + for key, value in terminal.items(): + upstream_terminal.setdefault(value, set()).add(key) + self._qlateral = read_qlats(forcing_parameters, self._dataframe.index, self.downstream_flowpath_dict, upstream_terminal) #Mask out all non-simulated waterbodies self._dataframe['waterbody'] = self.waterbody_null #This also remaps the initial NHDComID identity to the HY_Features Waterbody ID for the reservoir...