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)
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",
            )