Subnetwork final

install.md

@@ -2,18 +2,21 @@
 
 This file provides more detailed instructions on how to install, configure, and get the project running. Summary instructions are provided in readme.md.
 
-#Step 1 - Clone the repository 
+#Step 1 -- Activate a python environment on your local machine 
+
+#Step 2 - Clone the repository 
 git clone [email protected]:NOAA-OWP/t-route.git
 
-#Step 2 - Compile necessary files
+#Step 3 - Compile necessary files
 cd t-route/src/python_routing_v02/
+
 ./compiler.sh
 
-#Step 3 - Run the model 
+#Step 4 - Run the model 
 
 python3 compute_nhd_routing_SingleSeg_v02.py
 
-#Step 3 - Alternative v01 model (defaults can be overwritten with custom input yaml)
+#Step 4 - Alternative v01 model (defaults can be overwritten with custom input yaml)
 cd .. 
 cd python_routing_v01
 python3 compute_nhd_routing_SingleSeg_v01.py

src/external_connections/next_gen_io.py

@@ -24,7 +24,8 @@ def read_catchment_lateral_flows(path):
           catchment_id = int(file_name_str_list[0][4 :])
           catchment_id_list.append(catchment_id)
           # Read the second column of a csv file and return a series. The index will be an autoincrementing range.
-          catchment_qlats = pd.read_csv(os.path.join(path, file_name), names=[catchment_id], usecols=[1], squeeze=True)
+          catchment_qlats = pd.read_csv(os.path.join(path, file_name), usecols=['Flow'], squeeze=True)
+          catchment_qlats.rename(catchment_id, inplace=True)
           ql.append(catchment_qlats)
 
     qlats = pd.concat(ql, axis='columns').T

src/external_connections/next_gen_network.py

@@ -1,6 +1,6 @@
 #!/usr/bin/env python
 # coding: utf-8
-# example usage: python next_gen_network.py -n sugar_creek_waterbody_data_subset.geojson
+# example usage: python next_gen_network.py -n ./data/coarse/flowpath_data.geojson -s cat-71 cat-42 cat-34 cat-26 cat-68 cat-67 cat-27 cat-52 -i ./data
 
 """This program reads in and modifies lateral flows from the next generation water modeling framework and feeds them into mc_reach.compute_network. 
 This program presupposes that the mc_reach has been compiled prior to running."""
@@ -14,6 +14,7 @@
 from functools import partial
 from itertools import chain, islice
 import next_gen_io
+import json
 
 def _handle_args():
     parser = argparse.ArgumentParser(
@@ -26,73 +27,118 @@ def _handle_args():
         help="Input catchment network",
         dest="supernetwork"
     )
-
+    parser.add_argument(
+        "-s",
+        "--subset",
+        help="Subset of catchments to process",
+        dest="subset",
+        nargs='+'
+    )
+    parser.add_argument(
+        "-i",
+        "--input",
+        help="path to lateral flow inputs",
+        dest="input"
+    )
     return parser.parse_args()
 
-root = pathlib.Path('../..').resolve()
+root = pathlib.Path('.').resolve()
 sys.path.append(r"../python_framework_v02")
 sys.path.append(r"../python_routing_v02")
 sys.path.append(r"../../test/input/next_gen")
 
 test_folder = pathlib.Path(root, "test")
-next_gen_input_folder = test_folder.joinpath("input", "next_gen")
 
-import nhd_network
-import nhd_io
+import troute.nhd_network as nhd_network
+import troute.nhd_io as nhd_io
 import mc_reach
 
 def main():
 
     args = _handle_args()
+
+    next_gen_input_folder = test_folder.joinpath("input", "next_gen")
+    if args.input:
+        next_gen_input_folder = pathlib.Path(args.input)
 
-    #The following 2 values are currently hard coded for this test domain
+    # The following 2 values are currently hard coded for this test domain
     nts = 720  # number of timestep = 1140 * 60(model timestep) = 86400 = day
     dt_mc = 300.0  # time interval for MC
 
-    #Currently tested on the Sugar Creek domain
-    ngen_network_df = nhd_io.read_geopandas(os.path.join(next_gen_input_folder, args.supernetwork))
-
-    #Create dictionary mapping each connection ID
-    ngen_network_dict = dict(zip(ngen_network_df.ID, ngen_network_df.toID))
+    # Currently tested on the Sugar Creek domain
+    ngen_network_df = nhd_io.read_geopandas( args.supernetwork )
+    if args.subset:
+        ngen_network_df = ngen_network_df[ ngen_network_df['realized_catchment'].isin(args.subset) ]
+
+    # Create dictionary mapping each connection ID
+    ngen_network_dict = dict(zip(ngen_network_df.id, ngen_network_df.toid))
+    #ngen_network_dict = dict(zip(ngen_network_df.ID, ngen_network_df.toID))
 
     def node_key_func(x):
-        return int(x[4:])
+        return int(x[3:])
 
-    #Extract the ID integer values
+    # Extract the ID integer values
     waterbody_connections = {node_key_func(k): node_key_func(v) for k, v in ngen_network_dict.items()}
 
-    #Convert dictionary connections to data frame and make ID column the index
+    # Convert dictionary connections to data frame and make ID column the index
     waterbody_df = pd.DataFrame.from_dict(waterbody_connections, orient='index', columns=['to'])
-
-    #Sort ID index column
+    # Sort ID index column
     waterbody_df = waterbody_df.sort_index()
 
     waterbody_df = nhd_io.replace_downstreams(waterbody_df, "to", 0)
 
     connections = nhd_network.extract_connections(waterbody_df, "to")
 
-    #Read and convert catchment lateral flows to format that can be processed by compute_network
+    # Read and convert catchment lateral flows to format that can be processed by compute_network
     qlats = next_gen_io.read_catchment_lateral_flows(next_gen_input_folder)
-
+    print(qlats)
     rconn = nhd_network.reverse_network(connections)
 
     subnets = nhd_network.reachable_network(rconn, check_disjoint=False)
+
+    # read the routelink file
+    nhd_routelink = nhd_io.read_netcdf("data/RouteLink_NHDPLUS.nc")
+    nhd_routelink['dt'] = 300.0
+
+    nhd_routelink.set_index("link", inplace=True)
 
-    waterbody_df['dt'] = 300.0
+    routelink_cols = {
+        "downstream": "to",
+        "dx": "Length",
+        "n": "n",
+        "ncc": "nCC",
+        "s0": "So",
+        "bw": "BtmWdth",
+        "tw": "TopWdth",
+        "twcc": "TopWdthCC",
+        "waterbody": "NHDWaterbodyComID",
+        "musk": "MusK",
+        "musx": "MusX",
+        "cs": "ChSlp",
+    }
 
-    #Setting all below to 1.0 until we can get the appropriate parameters
-    waterbody_df['bw'] = 1.0
-    waterbody_df['tw'] = 1.0
-    waterbody_df['twcc'] = 1.0
-    waterbody_df['dx'] = 1.0
-    waterbody_df['n'] = 1.0
-    waterbody_df['ncc'] = 1.0
-    waterbody_df['cs'] = 1.0
-    waterbody_df['s0'] = 1.0
+    routelink_cols = dict([(value,key) for key, value in routelink_cols.items() ])
+
+    nhd_routelink.rename(columns=routelink_cols, inplace=True)
+
+    with open(next_gen_input_folder/'coarse/crosswalk.json') as f:
+        crosswalk_data = json.load(f)
+    waterbody_df['comid'] = waterbody_df.apply(lambda x: crosswalk_data['cat-' + str(x.name)]['outlet_COMID'], axis=1)
+
+    waterbody_df = waterbody_df.join(nhd_routelink, on='comid', how='left')
+
+    del nhd_routelink
+
+
+    # initial conditions, assume to be zero
+    # TO DO: Allow optional reading of initial conditions from WRF
+    q0 = pd.DataFrame(
+        0, index=waterbody_df.index, columns=["qu0", "qd0", "h0"], dtype="float32"
+    )
 
     #Set types as float32
     waterbody_df = waterbody_df.astype({"dt": "float32", "bw": "float32", "tw": "float32", "twcc": "float32", "dx": "float32", "n": "float32", "ncc": "float32", "cs": "float32", "s0": "float32"})
-
+    
     subreaches = {}
 
     for tw, net in subnets.items():
@@ -104,9 +150,12 @@ def node_key_func(x):
     for twi, (tw, reach) in enumerate(subreaches.items(), 1):
         r = list(chain.from_iterable(reach))
         data_sub = waterbody_df.loc[r, ['dt', 'bw', 'tw', 'twcc', 'dx', 'n', 'ncc', 'cs', 's0']].sort_index()
+        #data_sub = waterbody_df.loc[r, ['dt', 'bw', 'tw', 'twcc', 'dx', 'n', 'ncc', 'cs', 's0']]
         qlat_sub = qlats.loc[r].sort_index()
+        q0_sub = q0.loc[r].sort_index()
+
         results.append(mc_reach.compute_network(
-            nts, reach, subnets[tw], data_sub.index.values, data_sub.columns.values, data_sub.values, qlat_sub.values
+            nts, reach, subnets[tw], data_sub.index.values, data_sub.columns.values, data_sub.values, qlat_sub.values, q0_sub.values
             )
         )