docs: def extractSubsetFromWatershed(subset_network_file,
                                     subset_network_river_id_field,
                                     watershed_file,
                                     watershed_network_river_id_field,
                                     out_watershed_subset_file):
          """
          Extract catchment by using subset network file.
          Use this after using either
          :func:`~RAPIDpy.gis.taudem.TauDEM.extractSubNetwork()`
          or :func:`~RAPIDpy.gis.taudem.TauDEM.extractLargestSubNetwork()`.
  
          Parameters
          ----------
          subset_network_file: str
              Path to the pre-subsetted stream network shapefile.
          subset_network_river_id_field: str
              The field name with the river ID in the stream network shapefile.
          watershed_file: str
              Path to the watershed shapefile.
          watershed_network_river_id_field: str
              Name of the field with the river ID in the watershed shapefile.
          out_watershed_subset_file: str
              The path to output the subset watershed shapefile.
  
  
          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")
              water_shp = os.path.join(output_directory,
                                      "watershed_shapefile.shp")
              out_shp = os.path.join(output_directory,
                                     "watershed_shapefile_subset.shp")
              TauDEM.extractSubsetFromWatershed(
                  subset_network_filenetwork_shp,
                  subset_network_river_id_field="LINKNO",
                  watershed_file=water_shp,
                  watershed_network_river_id_field="LINKNO",
                  out_watershed_subset_file=out_shp)
  
          """
          subset_network_shapefile = ogr.Open(subset_network_file)
          subset_network_layer = subset_network_shapefile.GetLayer()
  
          ogr_watershed_shapefile = ogr.Open(watershed_file)
          ogr_watershed_shapefile_lyr = ogr_watershed_shapefile.GetLayer()
          ogr_watershed_shapefile_lyr_defn = \
              ogr_watershed_shapefile_lyr.GetLayerDefn()
  
          number_of_features = ogr_watershed_shapefile_lyr.GetFeatureCount()
          watershed_rivid_list = np.zeros(number_of_features, dtype=np.int32)
          for feature_idx, watershed_feature in \
                  enumerate(ogr_watershed_shapefile_lyr):
              watershed_rivid_list[feature_idx] = \
                  watershed_feature.GetField(watershed_network_river_id_field)
  
          shp_drv = ogr.GetDriverByName('ESRI Shapefile')
          # Remove output shapefile if it already exists
          if os.path.exists(out_watershed_subset_file):
              shp_drv.DeleteDataSource(out_watershed_subset_file)
  
          subset_watershed_shapefile = \
              shp_drv.CreateDataSource(out_watershed_subset_file)
          subset_watershed_layer = \
              subset_watershed_shapefile.CreateLayer(
                  '',
                  ogr_watershed_shapefile_lyr.GetSpatialRef(),
                  ogr.wkbPolygon)
          # Add input Layer Fields to the output Layer if it is the one we want
          for iii in xrange(ogr_watershed_shapefile_lyr_defn.GetFieldCount()):
              subset_watershed_layer.CreateField(
                  ogr_watershed_shapefile_lyr_defn.GetFieldDefn(iii))
          subset_watershed_layer_defn = subset_watershed_layer.GetLayerDefn()
  
          for drainage_line_feature in subset_network_layer:
              try:
                  watershed_feature_index = \
                      np.where(watershed_rivid_list ==
                               drainage_line_feature.GetField(
                                   subset_network_river_id_field))[0][0]
              except IndexError:
                  log("{0} {1} not found ...".format(
                      subset_network_river_id_field,
                      drainage_line_feature.GetField(
                          subset_network_river_id_field)))
                  continue
  
              subset_feature = \
                  ogr_watershed_shapefile_lyr.GetFeature(watershed_feature_index)
              # add to list
              new_feat = ogr.Feature(subset_watershed_layer_defn)
  
              # Add field values from input Layer
              for iii in \
                      xrange(ogr_watershed_shapefile_lyr_defn.GetFieldCount()):
                  new_feat.SetField(
                      subset_watershed_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
              subset_watershed_layer.CreateFeature(new_feat)