From f85ace18ff9a9b947f61bbf6420b0af40e20961b Mon Sep 17 00:00:00 2001
From: Laura LINDEPERG <laura.lindeperg@irstea.priv>
Date: Tue, 25 May 2021 11:07:02 +0200
Subject: [PATCH] Budyko Analysis

---
 budyko_equation.py | 92 ++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 92 insertions(+)
 create mode 100644 budyko_equation.py

diff --git a/budyko_equation.py b/budyko_equation.py
new file mode 100644
index 0000000..1a09170
--- /dev/null
+++ b/budyko_equation.py
@@ -0,0 +1,92 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Fri May 21 15:45:10 2021
+
+@author: laura.lindeperg
+"""
+
+import geopandas as gpd
+import pandas as pd
+import matplotlib.pyplot as plt
+import numpy as np
+import seaborn as sns
+
+
+
+# **************************** Data *****************************
+
+
+# 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)
+
+
+# Watersheds
+shp_watersheds_path = 'C:/Users/laura.lindeperg/Documents/DonneesLaura/Watersheds/616_Catchments.shp'
+shp_watersheds = gpd.read_file(shp_watersheds_path)
+
+rrse_241_path = 'C:/Users/laura.lindeperg/Documents/DonneesLaura/STATIONS/Stations241_RRSE.csv'
+rrse_241 = pd.read_csv(rrse_241_path)
+rrse_241 = rrse_241.rename(columns = {'levels.RRSE.':'Code'})
+rrse_241 = rrse_241.loc[:, 'Code']
+
+
+# Geol
+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)
+
+
+# Border catchments
+border_catchments = ['A3792010', 'A3902010', 'B4224310', 'B4601010', 'B4631010', 'D0156510', 'O1900010', 'O2620010', 'O6140010', 'U2512010', 'U2542010', 'U2722010','V1000010', 'V1020010']
+
+# Catchments size filter
+size = shp_watersheds.loc[shp_watersheds.loc[:, 'S_km2'] < 5000]
+
+# Main geol proportion filter
+proportion_geol = df_geol.loc[df_geol.loc[:, 'maingeol_proportion'] > 0.70]
+
+# Weird hydrographs filter
+weird_catchments = pd.DataFrame(np.array(['H6201010', 'L4010710', 'Y1232010']), columns=['Code'])
+weird_catchments=weird_catchments.loc[:,'Code']
+
+geol = df_geol.loc[:, ['code', 'maingeol_description']]
+my_data = pd.merge(df_hydro_sig, geol, 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]
+filtered_watersheds = filtered_watersheds.loc[filtered_watersheds.loc[:, 'code'].isin(size.loc[:,'Code'])==True]
+filtered_watersheds = filtered_watersheds.loc[filtered_watersheds.loc[:, 'code'].isin(proportion_geol.loc[:,'code'])==True]
+filtered_watersheds = filtered_watersheds.loc[filtered_watersheds.loc[:, 'code'].isin(weird_catchments)==False]
+
+
+
+def budyko_equation(x):
+    # x=np.array(x)
+    y = np.sqrt(x*np.tanh(1/x)*(1-np.exp(-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.
+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)
+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')
+
+# ax.scatter(ET0_P, E_P)
+# plt.legend(loc='center left', bbox_to_anchor=(1.04, 0.5))
+
+sns.relplot(x=1/filtered_watersheds.aridity_ratio, y=1-filtered_watersheds.runoff_ratio)
+
+sns.relplot(x=1/df_hydro_sig.aridity_ratio, y=1-df_hydro_sig.runoff_ratio)
+
-- 
GitLab