##Finding Overlapping Catchments
import geopandas as gpd
import matplotlib.pyplot as plt
import os
import fnmatch
import pandas as pd
import numpy as np
import itertools
import sys

data_dir = r"data_dir"
result_dir = r"result_dir"
data_dir = str(sys.argv[1])
result_dir = str(sys.argv[2])
result_filename = str(sys.argv[3])

c_shp=[] 
list_files = os.listdir(data_dir) 
for each_file in list_files :  
  if fnmatch.fnmatch(each_file, "*.shp"): 
        print(each_file)
        c_shp.append(each_file) 

c_combos = list(itertools.combinations(c_shp, 2))
no_overlap_result=[]
overlap_result=[]
for c in c_combos: 
   
    file_name_0 = c[0] 
    file_name_1 = c[1] 
    data_0 = gpd.read_file(os.path.join(data_dir, file_name_0)) 
    
    data_1 = gpd.read_file(os.path.join(data_dir, file_name_1))
    data_0 = data_0.to_crs(crs=3857)
    data_1 = data_1.to_crs(crs=3857)
    result = gpd.overlay(data_0, data_1, how='intersection', keep_geom_type=False) 
    result['area'] = result['geometry'].area
    if result.shape[0]> 0:
           overlap_result.append([c[0],c[1],result.area.values[0]]) 
    else:
           no_overlap_result.append([c[0],c[1]]) 
    
df= pd.DataFrame(overlap_result, columns = ['catchment1', 'catchment2', 'area'])
print(df)  
results_full_filename = os.path.join(result_dir, result_filename) 
df.to_csv(results_full_filename, sep=';', index=False)
####END HERE####

#%%
##Boxplot for all variables (Here we considered T as an example)
import matplotlib.pyplot as plt
import pandas as pd
import os

data1 = pd.read_csv("data_dir_Obs")
data2 = pd.read_csv("data_dir_Mixed")
data3 = pd.read_csv("data_dir_GLDAS")
data4 = pd.read_csv("data_dir_ERA5")


datasets = [data1["Mean T (C)"], data2["Mean T (C)"], data3["Mean T (C)"], data4["Mean T (C)"]]
box_colors = ["blue", "brown", "green", "purple"]
median_color = "black"
# Create the boxplot
plt.figure(figsize=(8, 6))
boxplot = plt.boxplot(datasets, patch_artist=True, showfliers=False)  

for patch, color in zip(boxplot['boxes'], box_colors):
    patch.set_facecolor(color)
for median in boxplot['medians']:
    median.set(color=median_color)

plt.ylabel("Mean Temperature (C)", fontsize=20)
labels = ["Obs", "Mixed", "GLDAS", "ERA5"]
plt.xticks(range(1, len(datasets) + 1), labels, fontsize=20)
plt.yticks(fontsize=20)

desktop_path = os.path.expanduser("data_dir")  
file_path = os.path.join(desktop_path, "Global-T.tiff")  
plt.savefig(file_path, dpi=300)
plt.show()
####END HERE####
#%%

##Budyko Plot (Here we consider Obs dataset as an example)
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import os
data =pd.read_csv("data_dir")
PETPratio = data['PET/P']
ETPratio = data['ET/P']
# Plot the data points
plt.scatter(PET/P, ET/P, color='purple', s=8, label='Data Points')
# Define a range of values for potential evapotranspiration (PET) fraction (ET/P)
budyko_curve_x = np.arange(0, 8, 0.05)
energy_limit_x = np.arange(0, 1.0001, 0.05)
x = np.arange(0, 1.0001, 0.05)
budyko_curve_y = np.power((budyko_curve_x*np.tanh(1/budyko_curve_x)*(1-np.exp(-budyko_curve_x))),0.5)
water_limit_y = 1+budyko_curve_x*0
energy_limit_y = energy_limit_x
y = 1 + x*0
plt.plot(budyko_curve_x,budyko_curve_y,  linestyle='--')
plt.plot(energy_limit_y,energy_limit_x, c='black')
plt.plot(budyko_curve_x,water_limit_y,c='black')
plt.plot(y,x,linestyle='-', c='red',label='_nolegend_')
plt.ylabel("ET/P")
plt.xlabel("PET/P")
plt.xticks(fontsize=16)
plt.yticks(fontsize=16)
desktop_path = os.path.expanduser("data_dir")  
file_path = os.path.join(desktop_path, "Global_budyko_Obs.tiff")  
plt.savefig(file_path, dpi=300)
plt.show()
####END HERE####
#%%

##Correlation (Here we consider Obs dataset as an example)
##T-ET
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import os
from scipy.stats import linregress
plt.rcParams.update({'font.size': 18})
file_path = "data_dir"
df = pd.read_csv(file_path)
df = df.dropna(subset=['Mean T (C)', 'Mean ET (mm/year)'])
correlation = df['Mean T (C)'].corr(df['Mean ET (mm/year)'])
print(f"Correlation between Mean T (C) and Mean ET (mm/year): {correlation}")
sns.regplot(x='Mean T (C)', y='Mean ET (mm/year)', data=df, color='blue', line_kws={'color':'black', 'linestyle':'--', 'linewidth':2})
slope, intercept, r_value, p_value, std_err = linregress(df['Mean T (C)'], df['Mean ET (mm/year)'])
line_eq = f"y = {slope:.2f}x + ({intercept:.2f})"
plt.text(0.02, 0.8, f"r2 = {r_value**2:.3f}", transform=plt.gca().transAxes)
plt.text(0.02, 0.9, line_eq, transform=plt.gca().transAxes)
plt.text(0.02, 0.7, f"p = {p_value:.3f}", transform=plt.gca().transAxes)
plt.xlim(-15,30)
plt.ylim(-50,2500)
desktop_path = os.path.expanduser("data_dir")  
file_path = os.path.join(desktop_path, "Obs.T-ET.Global.tiff")  
plt.savefig(file_path, dpi=300)
plt.show()

##T-R
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import os
from scipy.stats import linregress
plt.rcParams.update({'font.size': 18})
file_path = "data_dir"
df = pd.read_csv(file_path)
df = df.dropna(subset=['Mean T (C)', 'Mean R (mm/year)'])
correlation = df['Mean T (C)'].corr(df['Mean R (mm/year)'])
print(f"Correlation between Mean T (C) and Mean R (mm/year): {correlation}")
sns.regplot(x='Mean T (C)', y='Mean R (mm/year)', data=df, color='blue', line_kws={'color':'black', 'linestyle':'--', 'linewidth':2})
slope, intercept, r_value, p_value, std_err = linregress(df['Mean T (C)'], df['Mean R (mm/year)'])
line_eq = f"y = {slope:.2f}x + ({intercept:.2f})"
plt.text(0.02, 0.8, f"r2 = {r_value**2:.3f}", transform=plt.gca().transAxes)
plt.text(0.03, 0.9, line_eq, transform=plt.gca().transAxes)
plt.text(0.03, 0.7, f"p = {p_value:.3f}", transform=plt.gca().transAxes)
plt.xlim(-15,30)
plt.ylim(-50,3500)
desktop_path = os.path.expanduser("data_dir")  
file_path = os.path.join(desktop_path, "Obs.T-R.Global.tiff")  
plt.savefig(file_path, dpi=300)
plt.show()

##P-R
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import os
from scipy.stats import linregress
plt.rcParams.update({'font.size': 18})
file_path = "data_dir"
df = pd.read_csv(file_path)
df = df.dropna(subset=['Mean P (mm/year)', 'Mean R (mm/year)'])
correlation = df['Mean P (mm/year)'].corr(df['Mean R (mm/year)'])
print(f"Correlation between Mean P (mm/year) and Mean R (mm/year): {correlation}")
sns.regplot(x='Mean P (mm/year)', y='Mean R (mm/year)', data=df, color='blue', line_kws={'color':'black', 'linestyle':'--', 'linewidth':2})
slope, intercept, r_value, p_value, std_err = linregress(df['Mean P (mm/year)'], df['Mean R (mm/year)'])
line_eq = f"y = {slope:.2f}x + ({intercept:.2f})"
plt.text(0.02, 0.8, f"r2 = {r_value**2:.3f}", transform=plt.gca().transAxes)
plt.text(0.03, 0.9, line_eq, transform=plt.gca().transAxes)
plt.text(0.03, 0.7, f"p = {p_value:.3f}", transform=plt.gca().transAxes)
plt.xlim(0,4000)
plt.ylim(-50,3500)
desktop_path = os.path.expanduser("data_dir")  
file_path = os.path.join(desktop_path, "Obs.P-R.Global.tiff")  
plt.savefig(file_path, dpi=300)
plt.show()

##P-ET
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import os
from scipy.stats import linregress
plt.rcParams.update({'font.size': 18})
file_path = "data_dir"
df = pd.read_csv(file_path)
df = df.dropna(subset=['Mean P (mm/year)', 'Mean ET (mm/year)'])
correlation = df['Mean P (mm/year)'].corr(df['Mean ET (mm/year)'])
print(f"Correlation between Mean P (mm/year) and Mean ET (mm/year): {correlation}")
sns.regplot(x='Mean P (mm/year)', y='Mean ET (mm/year)', data=df, color='blue', line_kws={'color':'black', 'linestyle':'--', 'linewidth':2})
slope, intercept, r_value, p_value, std_err = linregress(df['Mean P (mm/year)'], df['Mean ET (mm/year)'])
line_eq = f"y = {slope:.2f}x + ({intercept:.2f})"
plt.text(0.02, 0.8, f"r2 = {r_value**2:.3f}", transform=plt.gca().transAxes)
plt.text(0.03, 0.9, line_eq, transform=plt.gca().transAxes)
plt.text(0.03, 0.7, f"p = {p_value:.3f}", transform=plt.gca().transAxes)
plt.xlim(0,4000)
plt.ylim(-50,2500)
desktop_path = os.path.expanduser("data_dir")  
file_path = os.path.join(desktop_path, "Obs.P-ET.Global.tiff")  
plt.savefig(file_path, dpi=300)
plt.show()