
# -------------              Script for IHA assessment          ------------ #
# --Instituto de Pesquisas Hidraulicas - Hidrologia em Grande Escala (HGE)
# --Author: Larissa de Castro Ribeiro (larissa.ribeirocr@gmail.com)
# --Advisor: Dr. Rodrigo Cauduro Dias de Paiva
# --Date: 19/11/2019
# --Updated in: 04/03/2021

# CONSIDERATIONS
# IHA indices generated by this script were developed according to the
# "Indicator of Hydrologic Alteration Version 7.1 User's Manual" by The Nature
# Conservancy (2009) with the following modifications:

# i) the discharges data are divided in hydrological years of maximum and minimum
# for each river segment to calculate the indices;
# ii) date information is transformed to angle in degrees to perform circular
# statistics and obtain maximum and minimum hydrological year dates.
# The alteration of maximum and minimum Julian days is calculated considering
# 182,5 days as a 100% alteration, with standard deviation given in days and the
# coefficient of variation being a ration of standard deviation by 182,5 days;
# c) minimum discharge related indices were calculated considering the hydrological
# year of minimums and maximum discharge related indices considered the
#hydrological year of maximums


# This script calculates the following IHA:
#1 - Mean monthly discharge: M_jan, M_fev, ...
#2 - Mean monthly discharge standard deviation: M_D_jan, M_D_fev,...
#3 - Mean monthly discharge coefficient of variation: M_CV_jan, M_CV_fev, ...

#4 - 1 day, 3 days, 30 days and 90 days minimum discharge: E_1_min1, E_1_min3, ...
#5 - Standard deviation of 1 day minimum discharge: E_1_min1_D, E_1_min3_D
#6 - Coefficient of variation of 1 day minimum discharge: E_1_min1_CV, E_1_min3_CV

#7 - 1 day, 3 days, 30 days and 90 days maximum discharge: E_2_max1, E_2_max3 ...
#8 - Standard deviation of 1 day maximum discharge: E_2_max1_D, E_2_max3_D
#9 - Coefficient of variation of 1 day maximum discharge: E_2_max1_CV, E_2_max3_CV

#10 - Date of minimum discharge: T_1
#11 - Standard deviation of Date of minimum discharge (days): T_1_D
#12 - Coefficient of variation of Date of minimum discharge (days): T_1_CV

#13 - Date of maximum discharge: T_2
#14 - Standard deviation of Date of maximum discharge (days): T_2_D
#15 - Coefficient of variation of Date of maximum discharge (days): T_2_CV

#16 - Ascension rate: G_1
#17 - Standard deviation of ascension rate: G_1_D
#18 - Coefficient of variation of ascension rate: G_1_CV

#19 - Recession rate: G_2
#20 - Standard deviation of recession rate: G_2_D
#21 - Coefficient of variation of recession rate: G_2_CV

#22 - Number of reversions: G_3
#23 - Standard deviation of number of reversions: G_3_D
#24 - Coefficient of variation of number of reversions: G_3_CV


########################### Importing libraries ###########################

import pandas as pd
import datetime
from scipy.stats import variation  
from scipy import stats
from scipy.stats import circstd
from scipy.stats import circvar
import numpy as np
import winsound
from scipy.io import loadmat
import warnings
warnings.filterwarnings("ignore")

######################################################
###                  SECTION 1 - IMPORT DISCHARGE DATA                      ###
######################################################
# Input data are series of discharges in dataframe format;
# rows are time intervals and columns are different river segments

###       CONVERTING INPUT FILE FROM .mat TO DataFrame         ###
##

mat = loadmat("folder/discharge_data.mat")
mdata = mat['...']

df = 0
df = pd.DataFrame(mdata[:,:])

del mat
del mdata

###      CONVERTING INPUT FILE FROM .xls TO DataFrame        ###
##

#df = pd.read_excel(r'../discharge_data.xlsx', engine='openpyxl')

#########################################################    
###                SECTION 2 - CREATING DATE INFORMATION                    ###
#########################################################
### Creating DataFrame in date format and number of days of discharge entries time step ###

# start = INPUT THE DATE OF THE FIRST DISCHARGE ENTRY
start = datetime.datetime.strptime("01/01/1980", "%d/%m/%Y").date()

numdays = len(df)
datelist = []
for x in range(0, numdays):
    datelist.append(start + datetime.timedelta(days=x))

# Creating a DataFrame with days, date, discharges and a river segment
##
# Creating days list (time step)
z = np.arange(len(df))
z = z + 1

df.insert(0, 'IT', z)

dias_contagem = list(df['IT'])

############################################################
###               SECTION 3 - CREATING EMPTY MATRIX FOR INDICES             ###
############################################################

########################### IHA MATRIX ###########################

z = list(df)
del z[0]

iha = pd.DataFrame(columns=['M_jan','M_fev','M_mar','M_abr',
                            'M_mai','M_jun','M_jul','M_ago',
                            'M_set','M_out','M_nov','M_dez',
                            'M_D_jan','M_D_fev','M_D_mar','M_D_abr',
                            'M_D_mai','M_D_jun','M_D_jul','M_D_ago',
                            'M_D_set','M_D_out','M_D_nov','M_D_dez',
                            'M_CV_jan','M_CV_fev','M_CV_mar',
                            'M_CV_abr','M_CV_mai','M_CV_jun',
                            'M_CV_jul','M_CV_ago','M_CV_set',
                            'M_CV_out','M_CV_nov','M_CV_dez',
                            'E_1_min1','E_1_min3','E_1_min7',
                            'E_1_min30','E_1_min90',
                            'E_1_min1_D','E_1_min3_D','E_1_min7_D',
                            'E_1_min30_D','E_1_min90_D',
                            'E_1_min1_CV','E_1_min3_CV','E_1_min7_CV',
                            'E_1_min30_CV','E_1_min90_CV',
                            'E_2_max1','E_2_max3','E_2_max7',
                            'E_2_max30','E_2_max90',
                            'E_2_max1_D','E_2_max3_D','E_2_max7_D',
                            'E_2_max30_D','E_2_max90_D',
                            'E_2_max1_CV','E_2_max3_CV','E_2_max7_CV',
                            'E_2_max30_CV','E_2_max90_CV',
                            'T_1','T_1_D','T_1_CV',
                            'T_2','T_2_D','T_2_CV','G_1', 
                            'G_2','G_3','G_1_D','G_2_D',
                            'G_3_D','G_1_CV','G_2_CV','G_3_CV','AnoH_Min','AnoH_Max'],index=z)

del x
del numdays
del start


######################################################    
###                   SECTION 4 - CALCULATION OF INDICES                    ###
######################################################
    
########################### LOOP FOR RIVER SEGMENTS ###########################
#Calculates indices for each river segment of the input data
#Converting segment discharge from column to list

before = datetime.datetime.now().isoformat(timespec='minutes')
before = pd.to_datetime(before) 
     
k=0
for k in range (54,len(z)):
    
    vazoes_trecho = list(df.iloc[:,k+1])
    
    matriz = 0
    matriz = pd.DataFrame({"dia": dias_contagem,"data": datelist,"vazao":vazoes_trecho})
    matriz['data'] = pd.to_datetime(matriz['data'])
          
    print("Calculando trecho " + str(k+1) + " de " + str(len(df.columns)))


    ################## Determining the beginning of the hydrological year ##################
    # For minimum
    ##
    q = pd.DataFrame()
    q['data'] = matriz['data']
    q['data'] = q['data'].dt.year
    q= q.groupby('data').size()
    
    yearlist = 0
    
    # Here you must input the series of discharge starting date

    yearlist = [datetime.datetime.strptime("1980", "%Y").year + index for index 
                in range(len(q))] ##TODO colocar o ano inicial
    yearlist_df = pd.DataFrame({"ano": yearlist})
    yearlist_df['ano'] = pd.to_datetime(yearlist_df['ano'])
    
    anos_min = 0
    anos_min = pd.DataFrame(index=list(yearlist),
                            columns = ['Data da vazao minima','Mes da vazao minima'])
    anos_min['Data da vazao minima'] = pd.to_datetime(anos_min['Data da vazao minima'])
      
    i = 0
    for i in range (0,len(yearlist)):
        ano1 = matriz.copy()
        ano1 = ano1.loc[(ano1['data'].dt.year==yearlist[i])]
        mini = ano1.copy()
        mini = mini.loc[(mini['vazao']==mini['vazao'].min())]
        
        if mini.size > 0:
            
            mini2 = mini.iloc[0] 
            anos_min.at[yearlist[i],'Data da vazao minima'] = mini2['data']
            anos_min['Mes da vazao minima'] = anos_min['Data da vazao minima'].dt.month
        else:
            pass
    
    anos_min = anos_min.groupby('Mes da vazao minima').size().sort_values(ascending=False)
    mes_ano_min = anos_min.index[0] - 5
    
    if mes_ano_min <= 0:
        mes_ano_min = mes_ano_min+12
    mes_anterior = mes_ano_min - 1
    if mes_anterior == 0:
        mes_anterior = 12
        
    iha.iloc[k,81] = mes_ano_min
       
    # Deleting data outside of the hydrological year range
    
    matriz_min = 0
    matriz_min = matriz.copy()
    
    i = 0
    for i in range (0,365):
        while matriz_min['data'][matriz_min.index[i]].month!=mes_ano_min:
            matriz_min = matriz_min.drop(index=[matriz_min.index[i]],axis=1)
        break
  
    i=0           
    matriz_min2 = matriz_min.iloc[::-1]      
    for i in range (0,365):
        while matriz_min2['data'][matriz_min2.index[i]].month!=mes_anterior:
            matriz_min2 = matriz_min2.drop(index=[matriz_min2.index[i]],axis=1)
        break
    
    # DataFrame with hydrological years of minimums
    matriz_min = matriz_min2.iloc[::-1]
    del(matriz_min2)
                       
    # Hydrological year of maximums
    ##   
    
    anos_max = 0
    anos_max = pd.DataFrame(index=list(yearlist),
                             columns = ['Data da vazao maxima','Mes da vazao maxima'])
    anos_max['Data da vazao maxima'] = pd.to_datetime(anos_max['Data da vazao maxima'])
      
    i = 0
    for i in range (0,len(yearlist)):
        ano1 = matriz.copy()
        ano1 = ano1.loc[(ano1['data'].dt.year==yearlist[i])]
        maxi = ano1.copy()
        maxi = maxi.loc[(maxi['vazao']==maxi['vazao'].max())]
        
        if maxi.size > 0:
            
            maxi2 = maxi.iloc[0] 
            anos_max.at[yearlist[i],'Data da vazao maxima'] = maxi2['data']
            anos_max['Mes da vazao maxima'] = anos_max['Data da vazao maxima'].dt.month
        else:
            pass  
         
    anos_max = anos_max.groupby('Mes da vazao maxima').size().sort_values(ascending=False)
    mes_ano_max = anos_max.index[0]-5
    
    if mes_ano_max <= 0:
        mes_ano_max = mes_ano_max+12
    mes_anterior = mes_ano_max - 1
    if mes_anterior == 0:
        mes_anterior = 12
        
    iha.iloc[k,82] = mes_ano_max
       
    # Deleting data outside of the hydrological year range
    
    matriz_max = 0
    matriz_max = matriz.copy()
    
    i = 0
    for i in range (0,365):
        while matriz_max['data'][matriz_max.index[i]].month!=mes_ano_max:
            matriz_max = matriz_max.drop(index=[matriz_max.index[i]],axis=1)
        break
             
    matriz_max2 = matriz_max.iloc[::-1]  
    i = 0   
    for i in range (0,365):
        while matriz_max2['data'][matriz_max2.index[i]].month!=mes_anterior:
            matriz_max2 = matriz_max2.drop(index=[matriz_max2.index[i]],axis=1)
        break
    
    # DataFrama with hydrological years of maximums
    matriz_max = matriz_max2.iloc[::-1]
    del(matriz_max2)
    
    ################## G1 - MEAN MAGNITUDE METRICS  ##################
    ##
    # mean monthly discharge for each year
    
    mag_medias_AH = 0
    mag_medias_AH = pd.DataFrame(index = np.arange(len(yearlist)),
                                  columns=['jan','fev','mar','abr',
                                          'mai','jun','jul','ago',
                                            'set','out','nov','dez'])
    
    meses = 12
    i = 0
    j = 0
    for i in range (0,len(yearlist)):
        mes = matriz.loc[
                    matriz['data'].dt.year == matriz.iloc[0,1].year+i
                    ]
        for j in range (0,meses):
            mes_i = mes.loc[(mes['data'].dt.month==j+1)]
            mag_medias_AH.iloc[i,j] = mes_i['vazao'].mean()
    
    meses2 = 12
    j = 0
    for j in range (0,meses2):
        iha.iloc[k,j] = mag_medias_AH.iloc[:,j].mean()
        iha.iloc[k,j+12] = mag_medias_AH.iloc[:,j].std(skipna=True)
        iha.iloc[k,j+24] =  mag_medias_AH.iloc[:,j].std(skipna=True
                )/mag_medias_AH.iloc[:,j].mean()
    
    ################## G2 - EXTREME MAGNITUDE METRICS ##################
    ##
    # for hydrological year of minimum
    
    yearlist2 = 0
    yearlist2 = matriz_min['data'].dt.year
    yearlist2 = pd.DataFrame(yearlist2)
    yearlist2 = yearlist2.groupby('data').size()
    
    mag_medias_AH_min = 0
    mag_medias_AH_min = pd.DataFrame(index = np.arange(len(yearlist2)-1),
                                      columns=['jan','fev','mar','abr',
                                              'mai','jun','jul','ago',
                                                'set','out','nov','dez'])
    
    # Calculating metrics for consecutive minimum discharges for a hydrological year of minimums

    mag_ext_AH_min = 0
    mag_ext_AH_min = pd.DataFrame(index = np.arange(len(yearlist2)-1),
                                  columns=['vazao minima de 1 dia',
                                        'vazao minima de 3 dias',
                                        'vazao minima de 7 dias',
                                        'vazao minima de 30 dias',
                                        'vazao minima de 90 dias'])
    
    ## Calculating timing for extreme minimum discharges
    ex_min = 0
    ex_min = pd.DataFrame(index = np.arange(len(yearlist2)-1),
                                  columns=['data', 'vazao', 'graus', 'dias no ano'])
    
    dias = 0
    dias = np.arange(0, len(ex_min))
    
    
    meses = 12
    
    try:
        i = 0
        j = 0
        for i in range (0,len(yearlist2)):
            mes = matriz_min.loc[
                    matriz_min['data'].dt.year == matriz_min.iloc[0,1].year+i
                    ]
            mes = mes.loc[mes['data'].dt.month >= mes_ano_min]
            mes2 = matriz_min.loc[
                    matriz_min['data'].dt.year == matriz_min.iloc[0,1].year+i+1
                    ]
            mes2 = mes2.loc[mes2['data'].dt.month < mes_ano_min]
            mes = mes.append(mes2)
            
            mag_ext_AH_min.iloc[i,0] = mes['vazao'].min()
            med3_mov_min = mes.copy()
            med3_mov_min['vazao'] = mes.vazao.rolling(window=3).mean()
            mag_ext_AH_min.iloc[i,1] = med3_mov_min['vazao'].min()
            med7_mov_min = mes.copy()
            med7_mov_min['vazao'] = mes.vazao.rolling(window=7).mean()
            mag_ext_AH_min.iloc[i,2] = med7_mov_min['vazao'].min()
            med30_mov_min = mes.copy()
            med30_mov_min['vazao'] = mes.vazao.rolling(window=30).mean()
            mag_ext_AH_min.iloc[i,3] = med30_mov_min['vazao'].min()
            med90_mov_min = mes.copy()
            med90_mov_min['vazao'] = mes.vazao.rolling(window=90).mean()
            mag_ext_AH_min.iloc[i,4] = med90_mov_min['vazao'].min()
            
            dias[i] = len(mes)
                                    
            for j in range (0,meses):
                if mes['vazao'].isna().sum().sum() < 364:
                    mes_i = mes.loc[(mes['data'].dt.month==j+1)]
                    mag_medias_AH_min.iloc[i, j] = mes_i['vazao'].mean()
                    mini2 = mes.loc[(mes['vazao']==mes['vazao'].min())]
                    mini2 = mini2.iloc[0]
                    ex_min['data'][i] = mini2['data']
                    ex_min['vazao'][i] = mini2['vazao']
                else:
                    ex_min['data'][i] = np.nan
                    ex_min['vazao'][i] = np.nan
                    
                
    except IndexError:
        pass
    
    ex_min['dias no ano'] = 0
    ex_min['dias no ano'] = pd.DataFrame(dias)
    extrem = 5
    
    mag_ext_AH_min = mag_ext_AH_min.dropna()
    
    j=0
    for j in range (0,extrem):
        iha.iloc[k,j+36] = mag_ext_AH_min.iloc[:,j].mean()
        iha.iloc[k,j+41] = stats.tstd(mag_ext_AH_min.iloc[:,j])
        if (mag_ext_AH_min.iloc[:,j]).sum() != 0:
            iha.iloc[k,j+46] = variation(mag_ext_AH_min.iloc[:,j], axis = 0)
        else:
            iha.iloc[k,j+46] = np.nan
    
    # For hydrological year of maximums
          
    yearlist3 = 0
    yearlist3 = matriz_max['data'].dt.year
    yearlist3 = pd.DataFrame(yearlist3)
    yearlist3 = yearlist3.groupby('data').size()
    
    mag_medias_AH_max = 0
    mag_medias_AH_max = pd.DataFrame(index = np.arange(len(yearlist3)-1),
                                  columns=['jan','fev','mar','abr',
                                          'mai','jun','jul','ago',
                                            'set','out','nov','dez'])
    
    # Calculating metrics of consecutive maximum discharges for hydrological year of maximums
    mag_ext_AH_max = 0
    mag_ext_AH_max = pd.DataFrame(index = np.arange(len(yearlist3)-1),
                                  columns=['vazao maxima de 1 dia',
                                        'vazao maxima de 3 dias',
                                        'vazao maxima de 7 dias',
                                        'vazao maxima de 30 dias',
                                        'vazao maxima de 90 dias'])
    
    # Calculating timing of extreme for maximum discharge
    
    ex_max = 0
    ex_max = pd.DataFrame(index = np.arange(len(yearlist)-1),
                                  columns=['data', 'vazao', 'graus', 'dias no ano'])
    dias = 0
    dias = np.arange(0, len(ex_max))
    
    meses = 12
    
    try:
        j = 0
        i = 0
        for i in range (0,len(yearlist3)):
            mes = matriz_max.loc[
                    matriz_max['data'].dt.year == matriz_max.iloc[0,1].year+i
                    ]
            mes = mes.loc[mes['data'].dt.month >= mes_ano_max]
            mes2 = matriz_max.loc[
                    matriz_max['data'].dt.year == matriz_max.iloc[0,1].year+i+1
                    ]
            mes2 = mes2.loc[mes2['data'].dt.month < mes_ano_max]
            mes = mes.append(mes2)
            
            mag_ext_AH_max.iloc[i,0] = mes['vazao'].max()
            med3_mov_max = mes.copy()
            med3_mov_max['vazao'] = mes.vazao.rolling(window=3).mean()
            mag_ext_AH_max.iloc[i,1] = med3_mov_max['vazao'].max()
            med7_mov_max = mes.copy()
            med7_mov_max['vazao'] = mes.vazao.rolling(window=7).mean()
            mag_ext_AH_max.iloc[i,2] = med7_mov_max['vazao'].max()
            med30_mov_max = mes.copy()
            med30_mov_max['vazao'] = mes.vazao.rolling(window=30).mean()
            mag_ext_AH_max.iloc[i,3] = med30_mov_max['vazao'].max()
            med90_mov_max = mes.copy()
            med90_mov_max['vazao'] = mes.vazao.rolling(window=90).mean()
            mag_ext_AH_max.iloc[i,4] = med90_mov_max['vazao'].max()
            dias[i] = len(mes)
                        
            for j in range (0,meses):               
                if mes['vazao'].isna().sum().sum() < 364:
                    mes_i = mes.loc[(mes['data'].dt.month==j+1)]
                    mag_medias_AH_max.iloc[i, j] = mes_i['vazao'].mean()
                    maxi2 = mes.loc[(mes['vazao']==mes['vazao'].min())]
                    maxi2 = maxi2.iloc[0]
                    ex_max['data'][i] = maxi2['data']
                    ex_max['vazao'][i] = maxi2['vazao']
                else:
                    ex_max['data'][i] = np.nan
                    ex_max['vazao'][i] = np.nan
                     
                
    except IndexError:
        pass
    
    ex_max['dias no ano'] = 0
    ex_max['dias no ano'] = pd.DataFrame(dias)
    extrem = 5
    
    mag_ext_AH_max = mag_ext_AH_max.dropna()
    
    j=0
    for j in range (0,extrem):
        iha.iloc[k,j+51] = mag_ext_AH_max.iloc[:,j].mean()
        iha.iloc[k,j+56] = stats.tstd(mag_ext_AH_max.iloc[:,j])
        iha.iloc[k,j+61] = variation(mag_ext_AH_max.iloc[:,j], axis = 0)
    
    
    ######################### G3 - TIMING OF EXTREMES #########################
    ##Calculating timing of extremes for minimum discharge
    
    #dates to degrees
    graus = 0
    graus = np.arange(0, len(ex_min)).astype(float)
    
       
    try:
        i = 0
        for i in range(0, len(ex_min)):
            if np.isnan(ex_min['vazao'][i]) == False:
                if ex_min['dias no ano'][i] == 365:
                    if ex_min['data'][i].month == 1:
                        graus[i] = (360/365)*ex_min['data'][i].day
                    elif ex_min['data'][i].month == 2:
                        graus[i] = (360/365)*(ex_min['data'][i].day+31)
                    elif ex_min['data'][i].month == 3:
                        graus[i] = (360/365)*(ex_min['data'][i].day+59)
                    elif ex_min['data'][i].month == 4:
                        graus[i] = (360/365)*(ex_min['data'][i].day+90)
                    elif ex_min['data'][i].month == 5:
                        graus[i] = (360/365)*(ex_min['data'][i].day+120)
                    elif ex_min['data'][i].month == 6:
                        graus[i] = (360/365)*(ex_min['data'][i].day+151)
                    elif ex_min['data'][i].month == 7:
                        graus[i] = (360/365)*(ex_min['data'][i].day+181)
                    elif ex_min['data'][i].month == 8:
                        graus[i] = (360/365)*(ex_min['data'][i].day+212)
                    elif ex_min['data'][i].month == 9:
                        graus[i] = (360/365)*(ex_min['data'][i].day+243)
                    elif ex_min['data'][i].month == 10:
                        graus[i] = (360/365)*(ex_min['data'][i].day+273)
                    elif ex_min['data'][i].month == 11:
                        graus[i] = (360/365)*(ex_min['data'][i].day+304)
                    elif ex_min['data'][i].month == 12:
                        graus[i] = (360/365)*(ex_min['data'][i].day+334)
        
                elif ex_min['dias no ano'][i] == 366:
                    if ex_min['data'][i].month == 1:
                        graus[i] = (360/366)*ex_min['data'][i].day
                    elif ex_min['data'][i].month == 2:
                        graus[i] = (360/366)*(ex_min['data'][i].day+31)
                    elif ex_min['data'][i].month == 3:
                        graus[i] = (360/366)*(ex_min['data'][i].day+60)
                    elif ex_min['data'][i].month == 4:
                        graus[i] = (360/366)*(ex_min['data'][i].day+91)
                    elif ex_min['data'][i].month == 5:
                        graus[i] = (360/366)*(ex_min['data'][i].day+121)
                    elif ex_min['data'][i].month == 6:
                        graus[i] = (360/366)*(ex_min['data'][i].day+152)
                    elif ex_min['data'][i].month == 7:
                        graus[i] = (360/366)*(ex_min['data'][i].day+182)
                    elif ex_min['data'][i].month == 8:
                        graus[i] = (360/366)*(ex_min['data'][i].day+213)
                    elif ex_min['data'][i].month == 9:
                        graus[i] = (360/366)*(ex_min['data'][i].day+244)
                    elif ex_min['data'][i].month == 10:
                        graus[i] = (360/366)*(ex_min['data'][i].day+274)
                    elif ex_min['data'][i].month == 11:
                        graus[i] = (360/366)*(ex_min['data'][i].day+305)
                    elif ex_min['data'][i].month == 12:
                        graus[i] = (360/366)*(ex_min['data'][i].day+335)
        
                else:
                    print("One of the years has less than 365 days, check input data")
                    
            else:
                graus[i] = np.nan
        
    except IndexError:
        pass

    ex_min['graus'] = 0
    ex_min['graus'] = pd.DataFrame(graus)
    
    sin_mean = 0
    sin_mean = np.sum(np.sin(np.radians(ex_min['graus'])))/len(ex_min)
    cos_mean = 0
    cos_mean = np.sum(np.cos(np.radians(ex_min['graus'])))/len(ex_min)
    
    timing_min = 0
    timing_min = np.degrees(np.arctan2(sin_mean, cos_mean))
    if timing_min < 0:
        timing_min = timing_min + 360
    if timing_min == 0:
        timing_min = 360
    timing_min = int(np.round(timing_min*365/360))
    graus = graus[np.logical_not(np.isnan(graus))]
    timing_min_std = int(round(np.degrees(circstd(graus*np.pi/180))*365/360))
    timing_min_var = circvar(graus*np.pi/180)
    
    # transforming degrees to month-day
    if timing_min > 0 and timing_min < 32:
        timing_min = "01-"+str(timing_min)
    elif timing_min >= 32 and timing_min< 60:
        timing_min = "02-"+str(timing_min-31)
    elif timing_min >= 60 and timing_min< 91:
        timing_min = "03-"+str(timing_min-59)
    elif timing_min >= 91 and timing_min< 121:
        timing_min = "04-"+str(timing_min-90)
    elif timing_min >= 121 and timing_min< 152:
        timing_min = "05-"+str(timing_min-120)
    elif timing_min >= 152 and timing_min< 182:
        timing_min = "06-"+str(timing_min-151)
    elif timing_min >= 182 and timing_min< 213:
        timing_min = "07-"+str(timing_min-181)
    elif timing_min >= 213 and timing_min< 244:
        timing_min = "08-"+str(timing_min-212)
    elif timing_min >= 244 and timing_min< 274:
        timing_min = "09-"+str(timing_min-243)
    elif timing_min >= 274 and timing_min< 305:
        timing_min = "10-"+str(timing_min-273)
    elif timing_min >= 305 and timing_min< 335:
        timing_min = "11-"+str(timing_min-304)
    elif timing_min >= 335 and timing_min<= 365:
        timing_min = "12-"+str(timing_min-334)
    
    
    
    iha.iloc[k,66] = timing_min
    iha.iloc[k,67] = timing_min_std
    iha.iloc[k,68] = timing_min_var
    
    # Calculating timing of extreme for maximum discharge
           
    # dates to degrees
    graus = 0
    graus = np.arange(0, len(ex_max)).astype(float)
    
    
    try:
        i = 0
        for i in range(0, len(ex_max)):
            if np.isnan(ex_max['vazao'][i]) == False:
                if ex_max['dias no ano'][i] == 365:
                    if ex_max['data'][i].month == 1:
                        graus[i] = (360/365)*ex_max['data'][i].day
                    elif ex_max['data'][i].month == 2:
                        graus[i] = (360/365)*(ex_max['data'][i].day+31)
                    elif ex_max['data'][i].month == 3:
                        graus[i] = (360/365)*(ex_max['data'][i].day+59)
                    elif ex_max['data'][i].month == 4:
                        graus[i] = (360/365)*(ex_max['data'][i].day+90)
                    elif ex_max['data'][i].month == 5:
                        graus[i] = (360/365)*(ex_max['data'][i].day+120)
                    elif ex_max['data'][i].month == 6:
                        graus[i] = (360/365)*(ex_max['data'][i].day+151)
                    elif ex_max['data'][i].month == 7:
                        graus[i] = (360/365)*(ex_max['data'][i].day+181)
                    elif ex_max['data'][i].month == 8:
                        graus[i] = (360/365)*(ex_max['data'][i].day+212)
                    elif ex_max['data'][i].month == 9:
                        graus[i] = (360/365)*(ex_max['data'][i].day+243)
                    elif ex_max['data'][i].month == 10:
                        graus[i] = (360/365)*(ex_max['data'][i].day+273)
                    elif ex_max['data'][i].month == 11:
                        graus[i] = (360/365)*(ex_max['data'][i].day+304)
                    elif ex_max['data'][i].month == 12:
                        graus[i] = (360/365)*(ex_max['data'][i].day+334)
        
                elif ex_max['dias no ano'][i] == 366:
                    if ex_max['data'][i].month == 1:
                        graus[i] = (360/366)*ex_max['data'][i].day
                    elif ex_max['data'][i].month == 2:
                        graus[i] = (360/366)*(ex_max['data'][i].day+31)
                    elif ex_max['data'][i].month == 3:
                        graus[i] = (360/366)*(ex_max['data'][i].day+60)
                    elif ex_max['data'][i].month == 4:
                        graus[i] = (360/366)*(ex_max['data'][i].day+91)
                    elif ex_max['data'][i].month == 5:
                        graus[i] = (360/366)*(ex_max['data'][i].day+121)
                    elif ex_max['data'][i].month == 6:
                        graus[i] = (360/366)*(ex_max['data'][i].day+152)
                    elif ex_max['data'][i].month == 7:
                        graus[i] = (360/366)*(ex_max['data'][i].day+182)
                    elif ex_max['data'][i].month == 8:
                        graus[i] = (360/366)*(ex_max['data'][i].day+213)
                    elif ex_max['data'][i].month == 9:
                        graus[i] = (360/366)*(ex_max['data'][i].day+244)
                    elif ex_max['data'][i].month == 10:
                        graus[i] = (360/366)*(ex_max['data'][i].day+274)
                    elif ex_max['data'][i].month == 11:
                        graus[i] = (360/366)*(ex_max['data'][i].day+305)
                    elif ex_max['data'][i].month == 12:
                        graus[i] = (360/366)*(ex_max['data'][i].day+335)
        
                else:
                    print("One of the years has less than 365 days, check input data")
                    
            else:
                graus[i] = np.nan
    
    except IndexError:
        pass
        
            
    ex_max['graus'] = 0
    ex_max['graus'] = pd.DataFrame(graus)
    
    sin_mean = 0
    sin_mean = np.sum(np.sin(np.radians(ex_max['graus'])))/len(ex_max)
    cos_mean = 0
    cos_mean = np.sum(np.cos(np.radians(ex_max['graus'])))/len(ex_max)
    
    timing_max = 0
    timing_max = np.degrees(np.arctan2(sin_mean, cos_mean))
    if timing_max < 0:
        timing_max = timing_max + 360
    if timing_max == 0:
        timing_max = 360
    timing_max = int(np.round(timing_max*365/360))
    graus = graus[np.logical_not(np.isnan(graus))]
    timing_max_std = int(round(np.degrees(circstd(graus*np.pi/180))*365/360))
    timing_max_var = circvar(graus*np.pi/180)
    
    # transforming degrees to month-day
    if timing_max > 0 and timing_max < 32:
        timing_max = "01-"+str(timing_max)
    elif timing_max >= 32 and timing_max< 60:
        timing_max = "02-"+str(timing_max-31)
    elif timing_max >= 60 and timing_max< 91:
        timing_max = "03-"+str(timing_max-59)
    elif timing_max >= 91 and timing_max< 121:
        timing_max = "04-"+str(timing_max-90)
    elif timing_max >= 121 and timing_max< 152:
        timing_max = "05-"+str(timing_max-120)
    elif timing_max >= 152 and timing_max< 182:
        timing_max = "06-"+str(timing_max-151)
    elif timing_max >= 182 and timing_max< 213:
        timing_max = "07-"+str(timing_max-181)
    elif timing_max >= 213 and timing_max< 244:
        timing_max = "08-"+str(timing_max-212)
    elif timing_max >= 244 and timing_max< 274:
        timing_max = "09-"+str(timing_max-243)
    elif timing_max >= 274 and timing_max< 305:
        timing_max = "10-"+str(timing_max-273)
    elif timing_max >= 305 and timing_max< 335:
        timing_max = "11-"+str(timing_max-304)
    elif timing_max >= 335 and timing_max<= 365:
        timing_max = "12-"+str(timing_max-334)
    
           
    iha.iloc[k,69] = timing_max
    iha.iloc[k,70] = timing_max_std
    iha.iloc[k,71] = timing_max_var
    
    ################## G5 - GRADIENT OF VARIABLES IN TIME ##################
    # Calculating ascension rate (GA), descent rate (GD) and reversion rate (RV)
        
    mag_G = 0
    mag_G = pd.DataFrame(index = np.arange(len(yearlist)), columns=['GA','GD','RV'])
    
    try:
        i = 0
    
        for i in range (0,len(mag_G)):
            j = 0
           
            ano_g = matriz.loc[
                    matriz['data'].dt.year == matriz.iloc[0,1].year+i]
            ano_g['GA'] = 0
            ano_g['GD'] = 0
            ano_g['GAD'] = 0
            ano_g['rev'] = 0
            
                               
            ano_g['GAD'] = ano_g['vazao']-ano_g['vazao'].shift(1)
            ano_g['GD'] = np.where(ano_g['GAD']<0,ano_g['GAD'],0)
            ano_g['GA'] = np.where(ano_g['GAD']>0,ano_g['GAD'],0)
            ano_g['rev'] = np.where(
                (ano_g['GAD']*ano_g['GAD'].shift(-1)) < 0,1,0
                )
    
            ano_g.iloc[0,3] = 0
            ano_g.iloc[0,4] = 0
            mag_G.iloc[i,0] = sum(ano_g['GA'])/(len(ano_g) -
            ano_g['GA'].isin([0]).sum())
            mag_G.iloc[i,1] = sum(-ano_g['GD'])/(len(ano_g) -
            ano_g['GD'].isin([0]).sum())
            mag_G.iloc[i,2] = sum(ano_g['rev'])
    
    except IndexError:
        pass
    
    # Filling the IHA matrix with rates values
    mag_G = mag_G.dropna()
    
    iha.iloc[k,72] = mag_G['GA'].mean()
    iha.iloc[k,73] = mag_G['GD'].mean()
    iha.iloc[k,74] = mag_G['RV'].mean()
    
    iha.iloc[k,75] = stats.tstd(mag_G['GA'])
    iha.iloc[k,76] = stats.tstd(mag_G['GD'])
    iha.iloc[k,77] = stats.tstd(mag_G['RV'])
    
    iha.iloc[k,78] = variation(mag_G['GA'])
    iha.iloc[k,79] = variation(mag_G['GD'])
    iha.iloc[k,80] = variation(mag_G['RV'])
    
   
##############################################################################
    
after = datetime.datetime.now().isoformat(timespec='minutes')
after = pd.to_datetime(after) 
tempo_execucao = after-before
print(tempo_execucao)

duration = 1000  # milliseconds
freq = 440  # Hz
winsound.Beep(freq, duration)


##########################         END          ###########################
 
