diff --git a/budyko_equation.py b/budyko_equation.py
index 1a0917014a2b85712f08727c5f3a132b08e45a62..7d5bdce41ca73c8237b02ad2fc4bdf86a8331746 100644
--- a/budyko_equation.py
+++ b/budyko_equation.py
@@ -17,7 +17,6 @@ import seaborn as sns
 
 
 # Load hydrological signatures df
-# df_hydro_sig_path = 'C:/Users/laura.lindeperg/Documents/INRAE_2021/616_stations_hydrosig_df.csv'
 df_hydro_sig_path = 'C:/Users/laura.lindeperg/Documents/INRAE_2021/616_stations_hydrosig_dim_df.csv'
 df_hydro_sig = pd.read_csv(df_hydro_sig_path)
 
@@ -36,6 +35,10 @@ rrse_241 = rrse_241.loc[:, 'Code']
 df_geol_path = 'C:/Users/laura.lindeperg/Documents/INRAE_2021/CODE/fhysa/616_stations_geol_df.csv'
 df_geol = pd.read_csv(df_geol_path)
 
+# Geomorpho
+df_geomorpho_path = 'C:/Users/laura.lindeperg/Documents/INRAE_2021/CODE/fhysa/616_stations_geomorpho_df.csv'
+df_geomorpho = pd.read_csv(df_geomorpho_path)
+
 
 # Border catchments
 border_catchments = ['A3792010', 'A3902010', 'B4224310', 'B4601010', 'B4631010', 'D0156510', 'O1900010', 'O2620010', 'O6140010', 'U2512010', 'U2542010', 'U2722010','V1000010', 'V1020010']
@@ -53,6 +56,9 @@ weird_catchments=weird_catchments.loc[:,'Code']
 geol = df_geol.loc[:, ['code', 'maingeol_description']]
 my_data = pd.merge(df_hydro_sig, geol, on='code')
 
+elevation = df_geomorpho.loc[:, ['code', 'mean_elevation']]
+my_data = pd.merge(df_hydro_sig, elevation, on='code')
+
 # Filtering workflow
 filtered_watersheds = my_data.loc[df_hydro_sig.loc[:, 'code'].isin(rrse_241)==True]
 filtered_watersheds = filtered_watersheds.loc[filtered_watersheds.loc[:, 'code'].isin(border_catchments)==False]
@@ -68,21 +74,23 @@ def budyko_equation(x):
     return y
 
 x = np.linspace(0, 1.2, 100)
+
 ET0_P = np.array(1/filtered_watersheds.aridity_ratio)
 E_P = np.array(1-filtered_watersheds.runoff_ratio)
 
 ET0_P = np.array(1/df_hydro_sig.aridity_ratio)
 E_P = np.array(1-df_hydro_sig.runoff_ratio)
 
-fig, ax = plt.subplots()  # Create a figure containing a single axes.
+fig, ax = plt.subplots(figsize=(8,5))  # Create a figure containing a single axes.
 ax.plot([0, 1, 1.2], [0, 1, 1])  # Plot some data on the axes.
 ax.plot(x, budyko_equation(x))
-sns.scatterplot(ax=ax, x=1/filtered_watersheds.aridity_ratio, y=1-filtered_watersheds.runoff_ratio, hue = 'maingeol_description', data = filtered_watersheds)
+sns.scatterplot(ax=ax, x=1/filtered_watersheds.aridity_ratio, y=1-filtered_watersheds.runoff_ratio, hue = 'mean_elevation', palette = 'viridis', data = filtered_watersheds).set(xlabel = 'ET0/P = 1/aridity_ratio', ylabel = 'E/P = 1-runoff_ratio')
+plt.legend(title = 'mean elevation [m]')
 for i in filtered_watersheds.code:
     watershed_i = filtered_watersheds.loc[filtered_watersheds.code == i]
     if float(1-watershed_i.runoff_ratio) > float(1/watershed_i.aridity_ratio):
         plt.text(1/watershed_i.aridity_ratio, 1-watershed_i.runoff_ratio, i, horizontalalignment='left', size='medium', color='black')
-
+        print(i)
 # ax.scatter(ET0_P, E_P)
 # plt.legend(loc='center left', bbox_to_anchor=(1.04, 0.5))