<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 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)