docs: def extractSubNetwork(network_file,
                            out_subset_network_file,
                            outlet_ids,
                            river_id_field,
                            next_down_id_field,
                            river_magnitude_field,
                            safe_mode=True):
          """
          Extracts a subset river network from the main river network based on
          the outlet IDs.
  
          Parameters
          ----------
          network_file: str
              Path to the stream network shapefile.
          out_subset_network_file: str
              Path to the output subset stream network shapefile.
          outlet_ids: list
              List of integers reperesenting the outlet IDs to be included in
              the subset stream network.
          river_id_field: str
              Name of the river ID field in the stream network shapefile.
          next_down_id_field: str
              Name if the field with the river ID of the next downstream
              river segment in the stream network shapefile.
          river_magnitude_field: str
              Name of the river magnitude field in the stream network shapefile.
          safe_mode: bool, optional
              If True, it will kill the simulation early before over taxing
              your computer. If you are confident your computer can handle it,
              set it to False.
  
  
          Here is an example of how to use this:
  
          .. code:: python
  
              import os
              from RAPIDpy.gis.taudem import TauDEM
  
  
              output_directory = '/path/to/output/files'
              network_shp = os.path.join(output_directory,
                                         "stream_reach_file.shp")
              out_shp = os.path.join(output_directory,
                                     "stream_reach_file_subset.shp")
  
              TauDEM.extractSubNetwork(
                  network_file=network_shp,
                  out_subset_network_file=out_shp,
                  outlet_ids=[60830],
                  river_id_field="LINKNO",
                  next_down_id_field="DSLINKNO",
                  river_magnitude_field="Magnitude",
              )
  
          """
          network_shapefile = ogr.Open(network_file)
          network_layer = network_shapefile.GetLayer()
          number_of_features = network_layer.GetFeatureCount()
          network_layer_defn = network_layer.GetLayerDefn()
          rivid_list = np.zeros(number_of_features, dtype=np.int32)
          next_down_rivid_list = np.zeros(number_of_features, dtype=np.int32)
          for feature_idx, drainage_line_feature in enumerate(network_layer):
              rivid_list[feature_idx] = \
                  drainage_line_feature.GetField(river_id_field)
              next_down_rivid_list[feature_idx] = \
                  drainage_line_feature.GetField(next_down_id_field)
  
          def getSubNetworkIDList(outlet_river_id,
                                  _rivid_list,
                                  _next_down_rivid_list):
              """
              Adds ids upstream of the outlet to a list
              """
              sub_network_index_list = []
              try:
                  for feature_ii in \
                          np.where(_next_down_rivid_list == outlet_river_id)[0]:
                      sub_network_index_list.append(feature_ii)
                      sub_network_index_list += \
                          getSubNetworkIDList(_rivid_list[feature_ii],
                                              _rivid_list,
                                              _next_down_rivid_list)
              except IndexError:
                  pass
              return sub_network_index_list
  
          original_recursion_limit = getrecursionlimit()
          try:
              main_sub_network_index_list = []
              for outlet_id in outlet_ids:
                  outlet_index = np.where(rivid_list == outlet_id)[0][0]
                  outlet_feature = network_layer.GetFeature(outlet_index)
                  outlet_magnitude = \
                      outlet_feature.GetField(river_magnitude_field)
                  if outlet_magnitude > original_recursion_limit:
                      if not safe_mode:
                          setrecursionlimit(outlet_magnitude)
                      else:
                          raise Exception("Current recursion limit {0} will not "
                                          "allow extraction for stream magnitude"
                                          " {1}. To override, set safe_mode to "
                                          "False ..."
                                          .format(original_recursion_limit,
                                                  outlet_magnitude))
                  main_sub_network_index_list.append(outlet_index)
                  main_sub_network_index_list += \
                      getSubNetworkIDList(outlet_id,
                                          rivid_list,
                                          next_down_rivid_list)
          except Exception:
              setrecursionlimit(original_recursion_limit)
              raise
  
          setrecursionlimit(original_recursion_limit)
  
          # Write out subset to new shapefile
          shp_drv = ogr.GetDriverByName('ESRI Shapefile')
          # Remove output shapefile if it already exists
          if os.path.exists(out_subset_network_file):
              shp_drv.DeleteDataSource(out_subset_network_file)
  
          network_subset_shp = shp_drv.CreateDataSource(out_subset_network_file)
          network_subset_layer = \
              network_subset_shp.CreateLayer('',
                                             network_layer.GetSpatialRef(),
                                             ogr.wkbLineString)
          # Add input Layer Fields to the output Layer if it is the one we want
          for iii in xrange(network_layer_defn.GetFieldCount()):
              network_subset_layer.CreateField(
                  network_layer_defn.GetFieldDefn(iii))
          network_subset_layer_defn = network_subset_layer.GetLayerDefn()
  
          for feature_index in main_sub_network_index_list:
              subset_feature = network_layer.GetFeature(feature_index)
              # add to list
              new_feat = ogr.Feature(network_subset_layer_defn)
  
              # Add field values from input Layer
              for iii in xrange(network_layer_defn.GetFieldCount()):
                  new_feat.SetField(
                      network_subset_layer_defn.GetFieldDefn(iii).GetNameRef(),
                      subset_feature.GetField(iii))
  
              # Set geometry as centroid
              geom = subset_feature.GetGeometryRef()
              new_feat.SetGeometry(geom.Clone())
              # Add new feature to output Layer
              network_subset_layer.CreateFeature(new_feat)