######################################################################
#             / ____|  ____|  __ \| |  | |  ____|  ____|             #
#            | |    | |__  | |__) | |__| | |__  | |__                #
#            | |    |  __| |  ___/|  __  |  __| |  __|               #
#            | |____| |____| |    | |  | | |____| |____              #
#             \_____|______|_|    |_|  |_|______|______|             #
######################################################################
#
# CEPHEE
# Copyright (C) 2024 Toulouse INP
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
# See the GNU General Public License for more details :
# <http://www.gnu.org/licenses/>.
#
######################################################################

# global
import rasterio.mask
import geopandas as gpd
from shapely.geometry import linestring
from shapely.ops import nearest_points
from scipy.spatial import cKDTree
from geopandas import GeoDataFrame
from collections import OrderedDict

# local
from core.Tools import *
from core.Reach import Reach

# optional for QGIS information display
try:
    from qgis.core import QgsMessageLog
except:
    pass


class ModelCatchment:
    """ TODO

    """
    def __init__(self, qgis_mode):
        self.init_mode = 0
        self.reach = []
        self.ordered_network = gpd.GeoDataFrame()
        self.hydro_network = None
        self.id_outlet = 0
        self.list_of_outlet = []
        self.outlet_point = None
        self.junction = None # geodataframe
        self.DEM_stack = None
        self.crs = None
        self.globalGrid = None
        self.globalDEM = None
        self.global_DEM_path = None
        self.Map = None
        self.outline = None
        self.qgis_mode = qgis_mode

    def display(self, message, level=0):
        if self.qgis_mode:
            QgsMessageLog.logMessage(message=message, tag="CEPHEE", level=level)
        else:
            if level == 0:
                print("Info: " + message)
            elif level == 1:
                print("Warning: " + message)
            elif level == 2:
                print("Error: " + message)

    def clipRiverFromDEM(self):
        """ Clip the hydrographic network as a function of boundary limit obtained from DEM file.
        Multilines are converted in linestring and added as a new reach
        """
        network = self.hydro_network
        network.reset_index()
        remove_index = []
        Xmin, Ymin, Xmax, Ymax =self.DEM_stack['global_extent']
        coords = ((Xmin,Ymin),(Xmax,Ymin),(Xmax,Ymax),(Xmin,Ymax))
        dem_polygon = Polygon(coords)
        rows_out = []  # Liste des lignes de sortie

        for _, row in network.iterrows():
            geom = row['geometry'].intersection(dem_polygon)
            if geom.is_empty:
                continue

            # Cas 1 : LineString simple
            if isinstance(geom, LineString):
                new_row = row.copy()
                new_row['geometry'] = geom
                rows_out.append(new_row)

            # Cas 2 : MultiLineString → on fait une ligne par sous-géométrie
            elif isinstance(geom, MultiLineString):
                for line in geom.geoms:
                    new_row = row.copy()
                    new_row['geometry'] = line
                    rows_out.append(new_row)

        # Reconstruction du GeoDataFrame complet avec les mêmes colonnes et CRS
        new_network = gpd.GeoDataFrame(rows_out, columns=network.columns, crs=self.crs)

        self.hydro_network = new_network

    def projectionOnDEM(self, verbose =True):
        """ Project the river reach (in the hydrographic network) on the DEM. The projection is made sequentially for
        each DEM file.

        The variable type_projection gives the method for projection:
            interpolation : use of interpolation 2D from scipy for each point on the DEM file
            raster : transform line in binary map (fast but only one measurement by DEM cell size)
        Metadata from BD carthage are saved and adapted


        :param interpolation_method:  method for 2D interpolation (nearest, cubic, linear)
        :type interpolation_method: string:
        """
        self.display("Hydro_network projection on DEMs", 0)
        network = self.hydro_network
        network.reset_index()

        pd_temp = pd.DataFrame(columns=['River', 'Reach', 'CdEntiteHy', 'Classe', 'NomEntiteH','geometry'])
        projected_network = GeoDataFrame(pd_temp,crs =self.crs)

        allX,allY,allid,allnpoint =[],[],[],[]
        for i, row in network.iterrows():
            line = row['geometry']
            allX += [coord[0] for coord in line.coords]
            allY += [coord[1] for coord in line.coords]
            allid += [i for _ in line.coords]
            allnpoint += [j for j in range(len(line.coords))]

        list_of_point_tot = []
        list_of_z = []
        list_of_id_tot = []
        list_of_dist_point_tot = []

        for idx, data_DEMi in enumerate(self.DEM_stack['data_DEM']):
            if verbose:
                self.display('Hydro_network projection on  DEM n°'+str(idx+1)+'/'
                             + str(len(self.DEM_stack['file_list'])), 0)

            DEM_polygon = Polygon(data_DEMi['polygon_coords'])
            current_file = self.DEM_stack['file_list'][idx]
            DEMi = rasterio.open(current_file, 'r', crs=self.crs)
            nodata = DEMi.nodata

            list_of_point = []
            list_of_id = []
            list_of_dist_point = []

            for x,y,id,npoint in zip(allX, allY, allid, allnpoint):

                if DEM_polygon.contains(Point(x, y)):
                    list_of_point.append(Point(x, y))
                    list_of_dist_point.append(npoint)
                    list_of_id.append(id)

            X = [pt.x for pt in list_of_point]
            Y = [pt.y for pt in list_of_point]
            Z = projOnDEM(X, Y, DEMi)
            Z= [float(z[0]) for z in Z]

            ind_nodata = [i for i,z in enumerate(Z) if z == nodata or z == float(-32768) or np.isnan(z)]
            for i in sorted(ind_nodata, reverse=True):
                del list_of_point[i]
                del list_of_dist_point[i]
                del list_of_id[i]
                del Z[i]

            list_of_z = list_of_z + Z
            list_of_point_tot = list_of_point_tot + list_of_point
            list_of_id_tot = list_of_id_tot + list_of_id
            list_of_dist_point_tot = list_of_dist_point_tot + list_of_dist_point

        n_lines_inc =0
        for n_lines, row in network.iterrows():
            id_selected_point = [ii for ii, d in enumerate(list_of_id_tot) if d == n_lines]
            list_point_selected = []
            list_of_dist_point_selected = []
            list_z_selected = []

            if len(id_selected_point) > 1:
                for id_point in id_selected_point:
                    if not list_of_z[id_point] == data_DEMi['no_data'] or list_of_z[id_point] > -99:
                        list_point_selected.append(list_of_point_tot[id_point])
                        list_z_selected.append(list_of_z[id_point])
                        list_of_dist_point_selected.append(list_of_dist_point_tot[id_point])
                list_x = [pt.x for pt in list_point_selected]
                list_y = [pt.y for pt in list_point_selected]
                list_indpoint = [npoint for npoint in list_of_dist_point_selected]
                # trie des points dans le sens initial
                sort_list_x = [i for _, i in sorted(zip(list_indpoint, list_x))]
                sort_list_y = [i for _, i in sorted(zip(list_indpoint, list_y))]
                sort_list_z_selected = [float(i) for _, i in sorted(zip(list_indpoint, list_z_selected))]

                line_proj = LineString(
                    [(x, y, z) for x, y, z in zip(sort_list_x, sort_list_y, sort_list_z_selected )])


                if not line_proj.is_empty:
                    try :
                        projected_network.loc[n_lines_inc, 'River'] = row['River']
                    except:
                        projected_network.loc[n_lines_inc, 'River'] = row['gid']

                    projected_network.loc[n_lines_inc, 'Reach'] = 0
                    projected_network.loc[n_lines_inc, 'CdEntiteHy'] = row['CdEntiteHy']
                    try:
                        projected_network.loc[n_lines_inc, 'Classe'] = row['Classe']
                    except:
                        projected_network.loc[n_lines_inc, 'Classe'] = row['Reach']
                    projected_network.loc[n_lines_inc, 'NomEntiteH'] = row['NomEntiteH']
                    projected_network.loc[n_lines_inc, 'geometry'] = line_proj
                    n_lines_inc += 1

        self.hydro_network = projected_network

    def order_network(self,min_dist_junction):
        """ Sort river reach from upstream to downstream considering averaged slope. Consecutive reaches without
        confluence are gathered (a reach can belong to only one river)

        :param min_dist_junction: Minimal distance between junctions
        :type min_dist_junction: float
        """
        self.display("Network ordering", 0)
        n_river = 0
        pd_temp=pd.DataFrame(columns=['River','Reach','CdEntiteHy','Classe','NomEntiteH'])
        self.ordered_network = GeoDataFrame(pd_temp, geometry=[],crs =self.crs)

        if self.DEM_stack:
            # Construct coords_min_max
            x_coords = []
            y_coords = []
            z_max = []
            indexes = []
            for i, row in self.hydro_network.iterrows():

                if not row['geometry'].is_empty:
                    if row['geometry'].coords[-1][2] > row['geometry'].coords[0][2]:
                        # invert LineString if z_end > z_start
                        row['geometry'] = LineString(row['geometry'].coords[::-1])
                        self.hydro_network.loc[i,'geometry'] = row['geometry']
                    x_coords.append((row['geometry'].coords[0][0], row['geometry'].coords[-1][0]))
                    y_coords.append((row['geometry'].coords[0][1], row['geometry'].coords[-1][1]))
                    z_max.append(row['geometry'].coords[0][2])
                    indexes.append(i)

            # Connecting reaches with connection upstream/downstream smaller min_dist_junction
            dist = np.zeros(len(z_max))
            while np.max(z_max) > 0:
                # Looking for the highest point not yet processed
                current_ind = np.argmax(z_max)
                self.ordered_network.loc[n_river] = self.hydro_network.loc[indexes[current_ind]].copy()
                current_linestring = self.hydro_network.loc[indexes[current_ind], 'geometry']
                z_max[current_ind] = - np.inf # z_max set to -inf when processed
                # Connecting all successive reaches
                while current_ind > -1:
                    # computing distance between downstream of current_ind and other available reaches
                    for j in range(len(z_max)):
                        if z_max[j] > 0:
                            dist[j] = ((x_coords[current_ind][1] - x_coords[j][0])**2
                                + (y_coords[current_ind][1] - y_coords[j][0])**2 ) **0.5
                        else:
                            dist[j] = np.inf

                    closest_ind = np.argmin(dist)
                    if dist[closest_ind] < min_dist_junction:
                        total_coords = (list(current_linestring.coords) +
                                        list(self.hydro_network.loc[indexes[closest_ind], 'geometry'].coords))
                        current_linestring = LineString(total_coords)
                        z_max[closest_ind] = - np.inf
                        current_ind = closest_ind
                    else:
                        current_ind = -1
                # Re-setting line string of the current reach of the ordered network
                self.ordered_network.loc[n_river, 'geometry'] = current_linestring
                n_river += 1
        else:
            if len(self.hydro_network) > 1:
                print('more than 1 reach, DEM must be provided')
            else:
                self.ordered_network = self.hydro_network
                self.ordered_network['River'] = self.ordered_network['gid']
                self.ordered_network.loc[0, 'Reach'] = 0
                self.list_of_outlet = [[self.ordered_network.loc[0, 'River'], [self.ordered_network.loc[0, 'gid']], 0]]


    def find_junctionAndOutlet(self, dist_min):
        """ Find the location of junctions between 2 reaches and find outlet on the DEM boundary for catchments.

        :param dist_min: minimum distance between 2 reaches to consider a junction
        :type dist_min: float
        """
        self.display("Finding junction and outlet", 0)
        gdf = self.ordered_network
        junction = GeoDataFrame(columns=['River1','River2','geometry'],crs =self.crs)
        outlet = []
        dist = []
        new_index = [i for i in range(len(gdf))]
        for i in range(len(gdf)):
            dist.append(self.outlet_point.distance(gdf.loc[i,'geometry']))

        id_outlet = np.argmin(dist) #indice du reach le plus proche de l'outlet
        gdf.index = new_index
        PointNearOutlet = nearest_points( gdf.loc[id_outlet,'geometry'], self.outlet_point)[0]
        c0 = np.transpose(np.array([ gdf.loc[id_outlet,'geometry'].xy[:][0], gdf['geometry'].iloc[id_outlet].xy[:][1]]))
        t0 = cKDTree(c0)
        c1 = np.transpose(np.array([PointNearOutlet.x,PointNearOutlet.y]))
        distance, neighbours = t0.query(c1)
        z_outlet = gdf.loc[id_outlet,'geometry'].coords[neighbours][2]
        if z_outlet <= gdf.loc[id_outlet,'geometry'].coords[0][2]: # on garde du point 1 jusqu'à l'exutoire
            if neighbours>1:
                gdf.loc[id_outlet, 'geometry'] = LineString([ gdf.loc[id_outlet,'geometry'].coords[jj] for jj in range(neighbours)])

        else: #on garde de l'exutoire jusqu'au dernier point
            if len(gdf['geometry'].iloc[id_outlet].coords)-neighbours>1:
                gdf.loc[id_outlet,'geometry'] = LineString([gdf.loc[id_outlet,'geometry'].coords[jj] for jj in range(neighbours,len(gdf.loc[id_outlet,'geometry'].coords))])
        self.outlet_point = PointNearOutlet #ajustement de l'outlet sur le reach
        # Suppression des reach à l'aval de l'exutoire
        remove_index = []
        for i in range(len(gdf)):  # boucle sur les reaches
            Zjj = [gdf.loc[i,'geometry'].coords[ii][2] for ii in range(len(gdf.loc[i,'geometry'].coords))]
            if np.max(Zjj) < z_outlet:
                self.display('Some internal reaches are lower than outlet')
                #remove_index.append(i)
        gdf.drop(remove_index, inplace=True)
        gdf.reset_index()


        minLinei = []
        ind_junction = 0
        for i in range(len(gdf)): #boucle sur les reaches
            minLinei.append(gdf.loc[i,'geometry'].coords[-1]) #stockage du point le plus à l'aval
            # calcul des distances entre les points de 2 reaches
            for j in range(i): #boucle sur les reaches avant celui considéré (symétrie des distances entre points)
                line1 = LineString([(x, y) for (x, y, z) in gdf.loc[i,'geometry'].coords])
                line2 = LineString([(x, y) for (x, y, z) in gdf.loc[j,'geometry'].coords])
                c0 = np.transpose(np.array([line1.xy[:][0],line1.xy[:][1]]))
                t0 = cKDTree(c0)
                c1 = np.transpose(np.array([line2.xy[:][0],line2.xy[:][1]]))
                distance, neighbours = t0.query(c1)
                ind_c1 = int(np.argmin(distance))
                ind_c0 = int(neighbours[ind_c1])
                if distance[ind_c1]<=dist_min:
                    Zjunction = gdf.loc[i,'geometry'].coords[ind_c0][2]
                    junction.loc[ind_junction, 'River1'] = gdf['River'].loc[i]#[[i, 0], [j, 0]]
                    junction.loc[ind_junction, 'River2'] = gdf['River'].loc[j]
                    junction.loc[ind_junction, 'geometry'] = Point(c0[ind_c0][0],c0[ind_c0][1],Zjunction)
                    ind_junction +=1
        self.junction = junction

        zmin= [minLinei[jj][2] for jj in range(len(minLinei))] # list des cotes de tous les points à l'aval des reachs
        ind_river = [gdf.loc[i,'River'] for i in range(len(gdf))]
        ind_reach = []
        ind_min = np.argmin(zmin)  #exutoire le plus bas
        ind_reach.append(ind_river[ind_min])
        for coord_min in minLinei:
            outlet.append([Point(coord_min), [], [], []])
        self.list_of_outlet += outlet

    def renameReachFromJunction(self,params):
        """ Rename reach considering junction.

        River is partitioned to get reach with a consistent discharge for each reach. For each river , reaches have
        a increasing id from the downstream to upstream. Id river is also corrected from 0 to Nriver. The name of reach
        refers to river and number of reach from downstream.

        :param compute_global: whether to compute global reach considering junction.
        :type compute_global: boolean
        """
        junction = self.junction
        gdf = self.ordered_network.copy()
        #on cherche la jonction la plus basse
        zjunction = []
        n_reach = 0
        junction_temp = GeoDataFrame(columns=['geometry', 'River1','Reach1','River2','Reach2','River3','Reach3','area'],crs =self.crs)
        junction_temp2 = GeoDataFrame(columns=['geometry', 'River1', 'Reach1', 'River2', 'Reach2', 'River3', 'Reach3', 'area'], crs=self.crs)
        pd_temp = pd.DataFrame(columns=['River', 'Reach', 'CdEntiteHy', 'Classe', 'NomEntiteH'])
        gdf_temp = GeoDataFrame(pd_temp, geometry=[],crs =self.crs)
        gdf_temp2  = GeoDataFrame(pd_temp, geometry=[],crs =self.crs)
        #decoupage des rivières en reach
        for i in range(len(gdf)):
            line = LineString([(x, y, z) for (x, y, z) in gdf.loc[i,'geometry'].coords])
            ind_river = gdf.loc[i,'River']
            points_to_cut = []
            for ind_junction in range(len(self.junction)):
                if self.junction.loc[ind_junction,'River1'] == ind_river or \
                        self.junction.loc[ind_junction,'River2'] == ind_river:
                    if self.junction.loc[ind_junction,'geometry'].distance(line) < params['N']['minDistJunction']:
                        points_to_cut.append(self.junction.loc[ind_junction,'geometry'])
                        zjunction.append(self.junction.loc[ind_junction,'geometry'].z)

            points_to_cut.append(Point(gdf['geometry'].loc[i].coords[0][0],gdf['geometry'].loc[i].coords[0][1]))
            points_to_cut.append(Point(gdf['geometry'].loc[i].coords[-1][0],gdf['geometry'].loc[i].coords[-1][1]))
            zjunction.append(gdf['geometry'].loc[i].coords[0][2])
            zjunction.append(gdf['geometry'].loc[i].coords[-1][2])
            line = LineString([(x, y,z) for (x, y, z) in gdf['geometry'].loc[i].coords])
            # Obtenir les coordonnées de la LineString
            coords = list(line.coords)
            # Projeter les points sur la ligne pour trouver les points les plus proches
            projected_coords = []
            for point,z in zip (points_to_cut,zjunction):
                # Calculer la position projetée le long de la ligne
                projected_distance = line.project(point)
                # Trouver les coordonnées réelles sur la ligne à cette distance projetée
                projected_point = line.interpolate(projected_distance)
                projected_coords.append((projected_point.x, projected_point.y,z))
            # Ajouter les points projetés à la liste des coordonnées de la ligne
            for projected_point in projected_coords:
                coords.append(projected_point)
            # Trier les points (coordonnées) dans l'ordre de la ligne en fonction de leur distance sur la ligne
            sorted_coords = sorted(coords, key=lambda coord: line.project(Point(coord)))
            segments = []
            current_segment = []
            for coord in sorted_coords:
                current_segment.append(coord)
                # Si le point actuel est un point projeté, terminer le segment ici
                if coord in projected_coords:
                    if len(current_segment) > 1 and current_segment[0][0] != current_segment[-1][0] \
                            and current_segment[0][1] != current_segment[-1][1]:
                        segments.append(LineString(current_segment))
                    # Démarrer un nouveau segment à partir de ce point
                    current_segment = [coord]

            for ind_reach, segment in enumerate(segments):
                gdf_temp.loc[n_reach,'River'] = ind_river
                gdf_temp.loc[n_reach,'Reach'] = ind_reach
                gdf_temp.loc[n_reach,'CdEntiteHy'] = gdf.loc[i,'CdEntiteHy']
                gdf_temp.loc[n_reach,'Classe'] = gdf.loc[i,'Classe']
                gdf_temp.loc[n_reach,'NomEntiteH'] = gdf.loc[i,'NomEntiteH']
                gdf_temp.loc[n_reach, 'geometry'] = segment
                n_reach +=1

        #modification des noms des reachs de l'aval vers l'amont
        ind_river_all = []
        for j in range(len(gdf_temp)): #cherche les index des rivières
            ind_river_all.append(gdf_temp['River'].loc[j])
        ind_river_all = np.unique(ind_river_all)
        for i in ind_river_all:
            Nreach = 0
            for j in range(len(gdf_temp)):
                if gdf_temp['River'].loc[j] == i:
                    Nreach+=1
            for j in range(len(gdf_temp)):
                if gdf_temp.loc[j,'River'] == i:
                    gdf_temp.loc[j,'Reach'] = Nreach-gdf_temp.loc[j,'Reach']-1
        #modification du nom des reachs dans la variable junction  découlant de la présence des jonctions
        for i in range(len(junction)):
            dist_end = []
            for j in range(len(gdf_temp)): #pour chaque jonction on cherche les 3 biefs les plus proches
                line1 = LineString([(gdf_temp['geometry'].loc[j].coords[jj][0], gdf_temp['geometry'].loc[j].coords[jj][1])
                                    for jj  in range(len(gdf_temp['geometry'].loc[j].coords))])
                p = junction.loc[i, 'geometry']
                dist_end.append(p.distance(line1))
            #on suppose uniquement 3 tronçons sur une jonction
            list_of_reach =[]
            list_of_river =[]
            list_of_z =[]

            for _ in range(3):
                ind_min=np.argmin(dist_end)
                z = [coord[2] for coord in gdf_temp.loc[ind_min,'geometry'].coords]
                list_of_reach.append(gdf_temp['Reach'].loc[ind_min])
                list_of_river.append(gdf_temp['River'].loc[ind_min])
                list_of_z.append(np.nanmin(z[:-1]))
                dist_end[ind_min] = np.inf

            sorted_reach = [[i, j] for _, i, j in sorted(zip(list_of_z,list_of_river, list_of_reach))]
            if params['C']['computeGlobal']:
            #utilisation de l'aire drainée et ajout de la valeur pour chaque jonction
                neighborhood = find_pixels_in_neighborhood(self.Map['acc'],self.Map['X'],self.Map['Y'],
                                                         junction.loc[i, 'geometry'].x, junction.loc[i, 'geometry'].y,
                                                         params['C']['window_size'])
                Area = np.max(neighborhood)
            else:
                Area =0

            #check consistency  (upstream and downstream reach can be inverted if they are short)
            groupes = OrderedDict()
            for river, reach in sorted_reach:
                groupes.setdefault(river, []).append(reach)

            sorted_reach2 = [[river, reach] for river in groupes for reach in sorted(groupes[river])]
            sorted_reach = [row[:] for row in sorted_reach2]

            if sorted_reach2[1][0] == sorted_reach2[2][0]: #la rivière commune doit être à l'aval
                sorted_reach[1][0] = sorted_reach2[0][0]
                sorted_reach[1][1] = sorted_reach2[0][1]
                sorted_reach[0][0] = sorted_reach2[1][0]
                sorted_reach[0][1] = sorted_reach2[1][1]



            junction_temp.loc[i, 'geometry'] = junction.loc[i, 'geometry']
            junction_temp.loc[i, 'River1'] = sorted_reach[0][0]
            junction_temp.loc[i, 'Reach1'] = sorted_reach[0][1]
            junction_temp.loc[i, 'River2'] = sorted_reach[1][0]
            junction_temp.loc[i, 'Reach2'] = sorted_reach[1][1]
            junction_temp.loc[i, 'River3'] = sorted_reach[2][0]
            junction_temp.loc[i, 'Reach3'] = sorted_reach[2][1]
            junction_temp.loc[i, 'area'] = Area

        junction_temp.sort_values(by='area')


        #ajout des noms de rivières dans outlet

        #recherche de la rivière exutoire
        dist_end = []
        for j in range(len(gdf_temp)):
            line1 = LineString([(x, y) for (x, y, z) in gdf_temp['geometry'].loc[j].coords])
            dist_line = self.outlet_point.distance(line1)
            if not np.isnan(dist_line):
                dist_end.append(self.outlet_point.distance(line1))
            else:
                dist_end.append(np.Inf)
        ind_min = np.argmin(dist_end)
        river1 = gdf_temp['River'].loc[ind_min ] #indice de la rivière connecté à l'exutoire
        reach1 = gdf_temp['Reach'].loc[ind_min] #indice du tronçon connecté à l'exutoire
        #recherche de toutes les rivières connectées à la rivière exutoire
        ind_river = [river1]
        Nind_prev = 0
        Nind = len(ind_river)

        while Nind > Nind_prev: #tant que l'on trouve une nouvelle rivière connecté au précédent
            Nind_prev = len(ind_river)
            #recherche des rivières liées à cet exutoire
            for j in range(len(junction_temp)):
                nriver = [junction_temp.loc[j,'River1'],
                        junction_temp.loc[j, 'River2'],
                        junction_temp.loc[j, 'River3']
                        ]
                for t in range(3):
                    if nriver[t] in ind_river:
                       #on ajoute les nouvelles rivières et on supprime celle en double
                        for ii in range(len(nriver)):
                            ind_river.append(nriver[ii])
                        ind_river = list(np.unique(ind_river))
            Nind=len(ind_river)
        # s'il ya eu calcul de l'accumulation, il y a un seul DEM d'où l'indice 0
        X =self.outlet_point.x
        Y = self.outlet_point.y

        if params['C']['computeGlobal']:
            neighborhood = find_pixels_in_neighborhood(self.Map['acc'], self.Map['X'], self.Map['Y'], X, Y, params['C']['window_size'])
            area = np.max(neighborhood)
        else:
            area = 0

        dist_point=[]
        for j in range(len(self.list_of_outlet)):
            dist_point.append(self.outlet_point.distance(self.list_of_outlet[j][0] ))

        self.id_outlet = np.argmin(dist_point)
        self.list_of_outlet[self.id_outlet ][3] = area
        self.list_of_outlet[self.id_outlet ][1] = ind_river
        self.list_of_outlet[self.id_outlet ][2] = [river1, reach1]

        nriver =0
        #keep only river connected to outlet
        for i, row in gdf_temp.iterrows():
            if  row['River'] in ind_river:
                gdf_temp2.loc[nriver] = row
                nriver +=1

        nriver = 0
        #keep only junction connected to outlet
        for i, row in junction_temp.iterrows():
            if row['River1'] in ind_river:
                junction_temp2.loc[nriver] = row
                nriver +=1

        gdf_temp2.set_geometry(col='geometry', inplace=True)
        junction_temp2.set_geometry(col='geometry', inplace=True)
        self.ordered_network = gdf_temp2
        self.junction = junction_temp2


    def interpolateReach(self, step):
        """ Change the distance between point of the center line of the reach. These points are those where sections
         will be plotted

        :param step: distance in meter between two floowing point on centerline of reach
        :type step: int
        """
        reach = self.reach
        ind_river = self.list_of_outlet[self.id_outlet ][1]
        # loop over all model reaches

        for j in range(len(reach)):
            # Select reach linked to the current outlet
            if reach[j].geodata['River'] in ind_river:
                geom = reach[j].geodata['geometry']
                num_vert = int(max(round(geom.length / step), 1))
                if geom.length > 1:
                    line_int = LineString([geom.interpolate(float(n) / num_vert, normalized=True) for n in range(num_vert + 1)])
                    reach[j].line_int=line_int
                    reach[j].Xinterp=distance_along_line(geom,line_int)
                else:
                    reach[j].line_int=LineString()

    def createReach(self):
        """Create the reach variable from data in the hydrographic network

        """
        # Create ordered_network if necessary
        if self.ordered_network.empty:
            self.ordered_network = GeoDataFrame(
                {'River': 0, 'Reach': 0, 'CdEntiteHy': 0, 'Classe': 1,
                 'NomEntiteH': 'river0', 'geometry': self.hydro_network.loc[0,'geometry']},crs =self.crs)
        # Create reaches connected to outlet
        network = self.ordered_network
        ind_river = self.list_of_outlet[self.id_outlet ][1]
        self.reach = []
        for _, row in self.ordered_network.iterrows():
            if row['River'] in ind_river:
                reach = Reach(str(row['NomEntiteH']) + "_" + str(row['Reach']))
                reach.geodata = row.copy()
                self.reach.append(reach)

    def compute_slope(self,npoly =None,data ='bed'):
        for reach in self.reach:
            if len(reach.section_list)>1:
                if reach.section_list[0].geom.has_z:
                    reach.compute_slope(npoly,data)
