read_TimeLoop.py 8.37 KiB
import pandas as pd
from itertools import compress


def get_TimeLoop_results(Input_folder, skiprow=[6, 7, 8, 9, 20830]):
    df = pd.read_csv(Input_folder + "/TimeLoop.dat", header=5,
                     sep="\t",
                     skiprows=skiprow)#[6, 7, 8, 9, 18638])#20830])#15717])#13160])
    df['Datetime'] = pd.to_datetime(df['ID'])
    return df


def cut_df_results(df, datetime_begin, datetime_end):
    if(df.index.name != 'Datetime'):
        df = df.set_index('Datetime', drop=True)
    df = df.truncate(before=datetime_begin, after=datetime_end)
    df = df.reset_index(drop=False)
    return df


def get_TimeLoop_dischargecomponents(df, Area_Wat):
    # transform title to get uniformed columns name
    df_heads = list(df.columns)
    BF_Bool = ['outRG1' in df_heads[i] for i in range(0, len(df_heads))]
    RD1_Bool = ['outRD1' in df_heads[i] for i in range(0, len(df_heads))]
    RD2_Bool = ['outRD2' in df_heads[i] for i in range(0, len(df_heads))]
    Runoff_Bool = ['Runoff_' in df_heads[i] for i in range(0, len(df_heads))]
    Precip_Bool = ['precip' in df_heads[i] for i in range(0, len(df_heads))]
    ETR_Bool = ['etact' in df_heads[i] for i in range(0, len(df_heads))]
    Rech_Bool = ['Percolation' in df_heads[i] for i in range(0, len(df_heads))]

    df['BF_mod'] = df[list(compress(df_heads, BF_Bool))[0]] / (1e6 * Area_Wat)
    df['SSF_mod'] = (df[list(compress(df_heads, RD1_Bool))[0]] + df[list(compress(df_heads, RD2_Bool))[0]]) / (
            1e6 * Area_Wat)
    df['Q_mod'] = df[list(compress(df_heads, Runoff_Bool))[0]] / (1e6 * Area_Wat)
    df['P'] = df[list(compress(df_heads, Precip_Bool))[0]]
    df['ETR'] = df[list(compress(df_heads, ETR_Bool))[0]]
    df['Rech_mod'] = df[list(compress(df_heads, Rech_Bool))[0]]
    df = df[['Datetime', 'Q_mod', 'BF_mod', 'SSF_mod', 'P', 'ETR', 'Rech_mod']]
    return df

def get_TimeLoop_dischargecomponents_specific(df, Area_Wat, ID):
    # transform title to get uniformed columns name
    df_heads = list(df.columns)
    BF_Bool = ['outRG1_'+str(ID) in df_heads[i] for i in range(0, len(df_heads))]
    RD1_Bool = ['outRD1_'+str(ID) in df_heads[i] for i in range(0, len(df_heads))]
    RD2_Bool = ['outRD2_'+str(ID) in df_heads[i] for i in range(0, len(df_heads))]
    Runoff_Bool = ['Runoff_'+str(ID) in df_heads[i] for i in range(0, len(df_heads))]
    Precip_Bool = ['precip_'+str(ID) in df_heads[i] for i in range(0, len(df_heads))]
    ETR_Bool = ['etact_'+str(ID) in df_heads[i] for i in range(0, len(df_heads))]
    Rech_Bool = ['Percolation_'+str(ID) in df_heads[i] for i in range(0, len(df_heads))]

    df['BF_mod'] = df[list(compress(df_heads, BF_Bool))[0]] / (1e6 * Area_Wat)
    df['SSF_mod'] = (df[list(compress(df_heads, RD1_Bool))[0]] + df[list(compress(df_heads, RD2_Bool))[0]]) / (
            1e6 * Area_Wat)
    df['Q_mod'] = df[list(compress(df_heads, Runoff_Bool))[0]] / (1e6 * Area_Wat)
    df['P'] = df[list(compress(df_heads, Precip_Bool))[0]]
    df['ETR'] = df[list(compress(df_heads, ETR_Bool))[0]]
    df['Rech_mod'] = df[list(compress(df_heads, Rech_Bool))[0]]
    df = df[['Datetime', 'Q_mod', 'BF_mod', 'SSF_mod', 'P', 'ETR', 'Rech_mod']]
    return df

def get_TimeLoop_dischargecomponents_specific_oldsimu(df, Area_Wat, Code, Reach_number):
    # transform title to get uniformed columns name
    df_heads = list(df.columns)
    BF_Bool = ['reachOutRG1_'+str(Reach_number) in df_heads[i] for i in range(0, len(df_heads))]
    RD1_Bool = ['reachOutRD1_'+str(Reach_number) in df_heads[i] for i in range(0, len(df_heads))]
    RD2_Bool = ['reachOutRD2_'+str(Reach_number) in df_heads[i] for i in range(0, len(df_heads))]
    Precip_Bool = [str(Code)+'precip' in df_heads[i] for i in range(0, len(df_heads))]
    ETR_Bool = [str(Code)+'actET' in df_heads[i] for i in range(0, len(df_heads))]
    Rech_Bool = [str(Code)+'percolation' in df_heads[i] for i in range(0, len(df_heads))]
    
#     print(str(sum(BF_Bool))+ ' & '+ str(sum(RD1_Bool))+ ' & '+ str(sum(RD2_Bool))+ ' & '+ str(sum(Precip_Bool))+ ' & '+ str(sum(ETR_Bool))+ ' & '+ str(sum(Rech_Bool))+'\n')
    
    df['BF_mod'] = df[list(compress(df_heads, BF_Bool))[0]] / (1e6 * Area_Wat)
    df['SSF_mod'] = (df[list(compress(df_heads, RD1_Bool))[0]] + df[list(compress(df_heads, RD2_Bool))[0]]) / (
            1e6 * Area_Wat)
    df['Q_mod'] = df['BF_mod'] + df['SSF_mod']
    df['P'] = df[list(compress(df_heads, Precip_Bool))[0]]
    df['ETR'] = df[list(compress(df_heads, ETR_Bool))[0]]
    df['Rech_mod'] = df[list(compress(df_heads, Rech_Bool))[0]]
    df = df[['Datetime', 'Q_mod', 'BF_mod', 'SSF_mod', 'P', 'ETR', 'Rech_mod']]
    return df

def check_mass_balance(df):
#     df = get_TimeLoop_results(Input_folder)
#     df = get_TimeLoop_dischargecomponents(df, Area_wat)
    P_mean = df['P'].mean()
    ETR_mod_mean = df['ETR'].mean()
    Q_mod_mean = df['Q_mod'].mean()
    Rech_mod_mean = df['Rech_mod'].mean()
    absolute_error = P_mean - ETR_mod_mean - Q_mod_mean
    relative_error = absolute_error / P_mean
    print('Abs error is ' + str(absolute_error) + ' mm/j \n Relative error is ' + str(relative_error) + '%\n')
    print('P_mean = ' + str(round(P_mean*365.25,0)) + ' mm/yr ' + 'ETR_mean = ' + str(round(ETR_mod_mean*365.25,0)) + ' mm/yr ' + 'Q_mean = ' + str(round(Q_mod_mean*365.25,0)) + ' mm/yr ' + 'Rech_mean = ' + str(round(Rech_mod_mean*365.25,0)) + ' mm/yr ' + '\n')
    return absolute_error, relative_error


def get_daily_interannual_average(df):
    # compute daily average
    if(df.index.name != 'Datetime'):
        df = df.set_index('Datetime', drop=True)
    df_dailyinterannualav = df.groupby(df.index.dayofyear).mean()
    if(max(df_dailyinterannualav.index)==366):
        df_dailyinterannualav = df_dailyinterannualav.drop(index=366)
    df = df.reset_index(drop=False)
    return df_dailyinterannualav

def get_daily_data(df):
    if(df.index.name != 'Datetime'):
        df = df.set_index('Datetime', drop=True)
    df = df.resample('D').mean()
    df = df.reset_index(drop=False)
    return df


def get_monthly_data(df):
    if(df.index.name != 'Datetime'):
        df = df.set_index('Datetime', drop=True)
    df = df.resample('M').mean()  # ,loffset='15d'
    df = df.reset_index(drop=False)
    return df


def merge_df(df_obs, df_mod):
    df_merged = pd.merge(df_obs, df_mod, on='Datetime')
    return df_merged


def get_performance(df_merged):
    from numpy import log10, corrcoef, sqrt
    df_merged = df_merged.dropna()
    NSE = 1 - sum((df_merged['Q_mod'] - df_merged['Q']) ** 2) / sum((df_merged['Q'] - df_merged['Q'].mean()) ** 2)
    NSE_log = 1 - sum((log10(df_merged['Q_mod']) - log10(df_merged['Q'])) ** 2) / sum(
        (log10(df_merged['Q']) - log10(df_merged['Q'].mean())) ** 2)
    alpha = df_merged['Q_mod'].std() / df_merged['Q'].std()
    beta = df_merged['Q_mod'].mean() / df_merged['Q'].mean()
    r = corrcoef(df_merged['Q_mod'], df_merged['Q'])
    KGE = 1 - sqrt((r[1, 0] - 1) ** 2 + (alpha - 1) ** 2 + (beta - 1) ** 2)
    return NSE, NSE_log, KGE

def get_bias(df_merged):
    Bias = ((df_merged['Q']-df_merged['Q_mod'])/df_merged['Q_mod']).mean()
    return Bias

def get_monthly_interannual_average(df):
    # data_set=pd.merge(data_fontainedutheil,data_surfex,on='Datetime')
    # compute monthly average
    df_monthlyinterannualav = get_monthly_data(df)
    df_monthlyinterannualav = df_monthlyinterannualav.groupby(df_monthlyinterannualav.index.month).mean()
    return df_monthlyinterannualav


def unit_test_function():
    Area_J2K = 232
    Input_folder = '/home/jean.marcais/Modeles/JAMS/jamsmodeldata/J2K_Gallaure/output/current/'
    df = get_TimeLoop_results(Input_folder, Area_J2K)
    df_dailyinterannualav = get_monthly_interannual_average(df)
    return df, df_dailyinterannualav


def plot_some_results(df):
    from matplotlib import pyplot as plt
    import matplotlib as mp
    mp.use('TkAgg')
    plt.plot(df['Datetime'], df['BF_mod'])
    plt.plot(df['Datetime'], df['Q_mod'])
    plt.plot(Q_real_df['Datetime'], Q_real_df['Q'])
    plt.show()

    plt.figure(1)
    plt.plot(df.index, df['BF_mod'])
    plt.plot(df.index, df['Q_mod'])
    plt.plot(Q_real_df.index, Q_real_df['Q'])
    plt.show()

    plt.figure(1)
    plt.plot(df_dailyinterannualav.index, df_dailyinterannualav['BF_mod'])
    plt.plot(df_dailyinterannualav.index, df_dailyinterannualav['Q_mod'])
    plt.plot(Q_real_df_dailyinterannualav.index, Q_real_df_dailyinterannualav['Q'])
    plt.show()

# Q_real_df.index=Q_real_df['Datetime']
# Q_real_df_dailyinterannualav=Q_real_df.groupby(Q_real_df.index.dayofyear).mean()