From 8dda4ed6a341afdbf7697a822c033392d9edf022 Mon Sep 17 00:00:00 2001
From: Poulard Christine <christine.poulard@irstea.fr>
Date: Sat, 8 May 2021 23:47:37 +0000
Subject: [PATCH] Upload New File, Atelier G
---
GenerateurCruesMaxAnnuelles.py | 378 +++++++++++++++++++++++++++++++++
1 file changed, 378 insertions(+)
create mode 100644 GenerateurCruesMaxAnnuelles.py
diff --git a/GenerateurCruesMaxAnnuelles.py b/GenerateurCruesMaxAnnuelles.py
new file mode 100644
index 0000000..8c95379
--- /dev/null
+++ b/GenerateurCruesMaxAnnuelles.py
@@ -0,0 +1,378 @@
+"""
+CPoulard ; mai 2021
+côté "métier" :
+ montrer le principe des formules Tchégodaiev / Hazen utilisées en hydrologie (TD 1ere année ENTPE) et la relation fréquence/période de retour
+ démo du caractère aléatoire des crues : on constitue progressivement une série d' "observations" issues d'une loi de Gumbel
+coté "Python" : mettre en oeuvre deux widgets de matplotlib (Checkbutton et Slider), faciles d'emploi mais au rendu pas très abouti ; définir les fonctions associées à des événements sur ces boutons ; manipuler les fonctions en tant qu'objet (fonction_en_cours = une des 3 fonctions disponibles...)
+"""
+
+from matplotlib import pyplot as plt
+import matplotlib.gridspec as gridspec
+from matplotlib.widgets import Slider, CheckButtons, Button # widgets du module matplotlib !
+import numpy as np
+from numpy.random import gumbel
+
+
+# FONCTIONS
+#
+# trois fonctions correspondant chacune à une manière de définir un échantillon de nb valeurs
+def echantillon_uniforme(nb):
+ return sorted(range(1, nb+1), reverse=True)
+
+def echantillon_constant(nb):
+ # non encore testé
+ return [1]*nb
+
+def echantillon_gumbel(nb, x0=196, gr=102):
+ # non encore testé
+
+ return gumbel(loc=x0, scale=gr, size=nb)
+
+class SerieGumbel:
+ T_Gumbel_plot = [1.1, 2, 5, 10, 20, 50, 100]
+
+ def __init__(self, gr, x0):
+ self.gr = gr
+ self.x0 = x0
+ self.serie_chrono = []
+ self.serie_triee = []
+ self.serie_couleurs = []
+
+ self.QT10 = self.QdeT(10)
+ self.QT50 = self.QdeT(50)
+ self.QT100 = self.QdeT(100)
+ self.QGumbel_plot = [self.QdeT(T) for T in SerieGumbel.T_Gumbel_plot]
+ self.gr_echantillon = None
+ self.x0_echantillon = None
+ print("quantiles ", self.QT10, self.QT50, self.QT100)
+
+ def attribution_couleur(self, Q):
+
+ if Q < self.QT10:
+ couleur = "black"
+ elif Q < self.QT50:
+ couleur = "green"
+ elif Q < self.QT100:
+ couleur = "blue"
+ else:
+ couleur = 'red'
+ return couleur
+
+ def QdeT(self, T):
+ """Donne la valeur estimée pour une période de retour donnée selon une loi de Gumbel
+ Args:
+ T (float) : période de retour
+ gr, x0 (float) : paramètres d'échelle et de position de la loi de Gumbel
+ Returns :
+ QdeT(float) : valeur estimée
+ """
+ Fnondepassement = 1 - (1 / T) # fréquence au non dépassement = 1 - f
+ u = -np.log(-np.log(Fnondepassement)) # log x s'obtient par np.log(x)
+ Q = self.x0 + self.gr * u
+ return Q
+
+ def TdeQ(self, Q):
+ """Donne la valeur estimée pour une période de retour donnée selon une loi de Gumbel
+ Args:
+ Q (float) : valeur
+ gr, x0 (float) : paramètres d'échelle et de position de la loi de Gumbel
+ Returns :
+ TdeQ(float) : période de retour estimée
+ """
+ u = (Q - self.x0) / self.gr # log x s'obtient par math.log(x)
+ Fnondepassement = np.exp(-1 * np.exp(-u))
+
+ T = 1 / (1 - Fnondepassement) # fréquence au non dépassement = 1 - f
+
+ return T
+
+ def ajout_obs(self, nb):
+ nouvelles_obs = gumbel(loc=self.x0, scale=self.gr, size=nb)
+ self.serie_chrono.extend(nouvelles_obs)
+ self.serie_couleurs.extend([self.attribution_couleur(Q) for Q in nouvelles_obs])
+ self.serie_triee.extend(nouvelles_obs)
+ self.serie_triee.sort(reverse=True)
+
+ self.gr_echantillon = 0.78 * np.std(self.serie_chrono) # gradex = paramètre d'échelle
+ self.x0_echantillon = np.mean(self.serie_chrono) - 0.577 * self.gr_echantillon # "ordonnée à l'origine" = paramètre de position
+
+
+# Loi de Gumbel (max annuels)
+#définition des lois d'ajustement
+def AjustementParametresDeGumbel(moyenne, ecart_type): # 2 arguments à calculer par ailleurs
+ """Ajuste une loi de Gumbel, usuelle en hydrologie pour ajuster des valeurs max annuelles
+ à partir "seulement" de la moyenne et l'écart-type de l'échantillon
+ Elle renvoie ses 2 paramètres gr et x0 calculés par la méthode des moments"
+ # ils sont donc moins sensibles aux points exceptionnels que la régression linéaire (distances au CARRE)
+ Args:
+ moyenne (float) : moyenne de l'échantillon
+ ecart_type (float) : ecart-type de l'échantillon
+ Returns :
+ gr(float) : paramètre d'échelle, gradex
+ x0(float) : paramètre de position
+
+ """
+
+ gr = 0.78 * ecart_type # gradex = paramètre d'échelle
+ x0 = moyenne - 0.577 * gr # "ordonnée à l'origine" = paramètre de position
+
+ return gr, x0
+
+
+def QdeT_Gumbel(T, gr, x0):
+ """Donne la valeur estimée pour une période de retour donnée selon une loi de Gumbel
+ Args:
+ T (float) : période de retour
+ gr, x0 (float) : paramètres d'échelle et de position de la loi de Gumbel
+ Returns :
+ QdeT(float) : valeur estimée
+ """
+ Fnondepassement = 1 - (1 / T) # fréquence au non dépassement = 1 - f
+ u = -np.log(-np.log(Fnondepassement)) # log x s'obtient par np.log(x)
+
+ QdeT = x0 + gr * u
+
+ return QdeT
+
+def TdeQ_Gumbel(Q, gr, x0):
+ """Donne la valeur estimée pour une période de retour donnée selon une loi de Gumbel
+ Args:
+ Q (float) : valeur
+ gr, x0 (float) : paramètres d'échelle et de position de la loi de Gumbel
+ Returns :
+ TdeQ(float) : période de retour estimée
+ """
+ u = (Q - x0) / gr # log x s'obtient par math.log(x)
+ Fnondepassement = np.exp(-1 * np.exp(-u))
+
+ TdeQ = 1 / (1 - Fnondepassement) # fréquence au non dépassement = 1 - f
+
+ return TdeQ
+
+
+# fonction paramétrée calculant la fréquence par une fonction paramétrée par a et b
+# puis 2 fonctions calculant la période de retour par la même fonction paramétrée et par une fonction avec valeurs fixes
+# EVIDEMMENT IL Y A REDONDANCE pour l'exemple ; la première fonction suffirait !
+
+def freq_emp_parametre(n_annees, a, b):
+ # pour la formule de Chegodayev, a = 0.3 et b=0.4
+
+ plotting_positions_freq = [ 1 - (((indice + 1) - a) / ((n_annees + b) )) for indice in range(n_annees)]
+ # au lieu d'envoyer une liste brute pour l'affichage,
+ # on crée une liste de freq formattées (:.2f = flottants 2 décimales) séparées par des "/"
+ print("Fréquences empiriques : ", " / ".join([f"{freq:.2f}" for freq in plotting_positions_freq]))
+ return plotting_positions_freq
+
+def T_emp_parametre(n_annees, a , b):
+ # T=1/freq de la fonction freq_emp_parametre, cas général de la précédente
+ plotting_positions_param = [((n_annees + b) / ((indice + 1) - a)) for indice in range(n_annees)]
+ print(f"Périodes de retour empiriques, a={a:.1f}, b={b:.1f}, N={n_annees} : ", " / ".join([f"{freq:.2f}" for freq in plotting_positions_param]))
+ return plotting_positions_param
+
+
+def update_vlines(courbe_en_vlines, x, ymin=None, ymax=None):
+ # mise à jour de courbes de type vlines
+ # https://stackoverflow.com/questions/29331401/updating-pyplot-vlines-in-interactive-plot
+ seg_old = courbe_en_vlines.get_segments()
+ if ymin is None:
+ ymin = seg_old[0][0, 1]
+ if ymax is None:
+ ymax = seg_old[0][1, 1]
+
+ seg_new = [np.array([[xx, ymin],
+ [xx, ymax]]) for xx in x]
+
+ courbe_en_vlines.set_segments(seg_new)
+
+def retracer_parametree(echantillon, a, b):
+ y_max = max(echantillon)
+ if en_periode_de_retour:
+ T_emp = T_emp_parametre(len(echantillon), a, b)
+ update_vlines(courbe_estim_T_vlines, T_emp, ymin=0, ymax=y_max)
+ ax_T.set_title(f"T empirique pour a={a:.2f} et b={b:.2f} ")
+ else:
+ #freq = [ 1 - (1 / T) for T in T_emp]
+ freq = freq_emp_parametre(len(echantillon), a , b)
+ update_vlines(courbe_estim_T_vlines, freq, ymin=0, ymax=y_max)
+ ax_T.set_title(f"fréquence empirique pour a={a:.2f} et b={b:.2f} ")
+
+ #ax.set_xbounds(lower=min(T_emp), upper=max(T_emp))
+ # on retrace
+ fig_pp.canvas.draw_idle()
+
+def retracer_Gumbel(serie_Gumbel):
+
+ if en_periode_de_retour:
+ courbe_Gumbel.set_xdata(SerieGumbel.T_Gumbel_plot)
+ else:
+ #freq = [ 1 - (1 / T) for T in T_emp]
+ freq = [ 1 - (1/T) for T in SerieGumbel.T_Gumbel_plot]
+ courbe_Gumbel.set_xdata(freq)
+
+ # on retrace
+ fig_pp.canvas.draw_idle()
+
+def retracer_estim_Gumbel(serie_Gumbel):
+ gr, x0 = AjustementParametresDeGumbel(np.mean(serie_Gumbel), np.std(serie_Gumbel))
+ if en_periode_de_retour:
+ T_Gumbel = [TdeQ_Gumbel(Q, gr, x0) for Q in serie_Gumbel]
+ courbe_estim_Gumbel.set_data(T_Gumbel, serie_Gumbel)
+ else:
+ # freq = [ 1 - (1 / T) for T in T_emp]
+ freq = freq_emp_parametre(len(serie_Gumbel), a, b)
+ courbe_estim_Gumbel.set_data(freq, serie_Gumbel)
+
+ # on retrace
+ fig_pp.canvas.draw_idle()
+
+
+# https://matplotlib.org/stable/gallery/widgets/slider_demo.html
+# https://stackoverflow.com/questions/13656387/can-i-make-matplotlib-sliders-more-discrete
+def update_a(val):
+ # val = valeur renvoyée par le widget
+ global a
+ a = val # comme on a précisé "global a", changer a dans la fonction = changer a du prog principal
+ retracer_parametree(serie_Gumbel.serie_triee, a, b)
+
+def update_b(val):
+ # val = valeur renvoyée par le widget
+ global b
+ b = val # comme on a précisé "global b", changer b dans la fonction = changer b du prog principal
+ retracer_parametree(serie_Gumbel.serie_triee, a, b)
+
+def update_series(event):
+ # bouton "ajouter des obs" activé
+
+ # comme on a précisé "global N_total, echantillon", sa modifications sera propagée jusque dans le prog principal
+ global N_total
+ N_total += N_serie
+ serie_Gumbel.ajout_obs(N_serie)
+
+ ax_chrono.clear()
+ ax_chrono.set_title(f"Qmax annuels observés sur {N_total} années")
+ courbe_Gumbel1, = ax_chrono.plot(range(len(serie_Gumbel.serie_chrono)), serie_Gumbel.serie_chrono, ls='--')
+ courbe_Gumbel = ax_chrono.scatter(range(len(serie_Gumbel.serie_chrono)), serie_Gumbel.serie_chrono,
+ c=serie_Gumbel.serie_couleurs, alpha=0.7, marker="*", s=20)
+ retracer_estim_Gumbel(serie_Gumbel.serie_triee)
+ fig_pp.canvas.draw_idle()
+ retracer_parametree(serie_Gumbel.serie_triee, a, b)
+
+def switch_freq(label):
+ global en_periode_de_retour
+
+ en_periode_de_retour = not en_periode_de_retour
+ retracer_parametree(serie_Gumbel.serie_triee, a, b)
+ retracer_Gumbel(serie_Gumbel.serie_triee)
+ retracer_estim_Gumbel(serie_Gumbel.serie_triee)
+ if en_periode_de_retour:
+ ax_T.set_xlim(left=0, right=N_total * 1.1)
+ ax_T.set_xlabel("Période de retour, années")
+ else:
+ # freq = [ 1 - (1 / T) for T in T_emp]
+ ax_T.set_xlim(left=0, right=1.1)
+ ax_T.set_xlabel("fréquence de non-dépassement")
+
+ fig_pp.canvas.draw_idle()
+
+#plt.ion()
+
+a_Tchego, b_Tchego = 0.3 , 0.4 # affectation de plusieurs valeurs simultanées par tuples
+# initialisation des paramètres de départ (modifiables par curseur)
+a = a_Tchego
+b = b_Tchego
+
+# on va constituer une série d' "observations" par ajouts de N_series obs à la fois
+N_serie = 10
+serie_Gumbel = SerieGumbel(x0=196, gr=102)
+serie_Gumbel.ajout_obs(N_serie)
+
+
+rangs = echantillon_uniforme(N_serie)
+
+N_total = N_serie
+en_periode_de_retour = True
+
+print(serie_Gumbel.serie_chrono)
+
+
+fig_pp = plt.figure()
+gs = gridspec.GridSpec(16, 2)
+ax_chrono = plt.subplot(gs[0:6, :])
+ax_T = plt.subplot(gs[6:12, :])
+ax_slide_a = plt.subplot(gs[15, 0])
+ax_slide_b = plt.subplot(gs[15, 1])
+ax_chb = plt.subplot(gs[14, 0])
+ax_bouton = plt.subplot(gs[14, 1])
+
+# ax_chb = fig.add_subplot()
+plt.subplots_adjust(wspace=1, hspace=0.5,left=0.1,top=0.85,right=0.9,bottom=0.1)
+fig_pp.canvas.set_window_title("ScE - Hydrologie -Démo Gumbel")
+
+fig_pp.suptitle(f"Démo Gumbel : tirages et leurs périodes de retour 'empiriques")
+
+ax_T.set_xlabel("Période de retour (années)")
+
+# LES OBS
+courbe_Gumbel1, = ax_chrono.plot(range(len(serie_Gumbel.serie_chrono)), serie_Gumbel.serie_chrono, ls='--')
+courbe_Gumbel = ax_chrono.scatter(range(len(serie_Gumbel.serie_chrono)), serie_Gumbel.serie_chrono, c=serie_Gumbel.serie_couleurs, alpha=0.7, label= f"Q(t) observé sur {N_total} années, ", marker="*", s=20)
+
+#LA COURBE DE GUMBEL, inchangée (sauf si tracé en fréquence...)
+courbe_Gumbel, = ax_T.plot(SerieGumbel.T_Gumbel_plot, serie_Gumbel.QGumbel_plot, color='red', alpha=0.7, linewidth=5, solid_capstyle='round',
+ zorder=10, label = "Q(T) du bassin versant", marker="None", ls='--', markersize=7)
+
+gr_est, x0_est = serie_Gumbel.gr_echantillon, serie_Gumbel.x0_echantillon
+
+courbe_estim_T_vlines = ax_T.vlines(x=T_emp_parametre(N_total, a, b), ymin=0, ymax=serie_Gumbel.serie_triee, color="purple", alpha=0.7, linewidth=2,
+ zorder=2, label= f"T ou f empirique, fonction de a et b", ls="--")
+courbe_estim_Gumbel, = ax_T.plot([TdeQ_Gumbel(Q, gr_est, x0_est) for Q in serie_Gumbel.serie_triee], serie_Gumbel.serie_triee, color='blue', alpha=0.7, label=" ajustement Gumbel sur obs",
+ zorder=10, marker="x", markersize=7)
+
+
+for ax_s in [ax_slide_a, ax_slide_b, ax_bouton, ax_chb]:
+ ax_s.xaxis.set_visible(False)
+ ax_s.yaxis.set_visible(False)
+
+for pos in ['right', 'top', 'bottom', 'left']:
+ ax_chb.spines[pos].set_visible(False)
+
+print("ici")
+
+# nom connu même hors de la fonction pour éviter le GC ?
+# a
+slider_a = Slider(
+ ax_slide_a, "paramètre a ", valmin=0.01, valmax = 0.99, valfmt='%0.2f', valinit=a, color="green")
+
+slider_a.on_changed(update_a)
+
+# b
+slider_b = Slider(
+ ax_slide_b, "paramètre b ", valmin=0.01, valmax = 0.99, valfmt='%0.2f', valinit=b, color="blue")
+
+slider_b.on_changed(update_b)
+
+# bouton pour ajouter N_serie valeurs
+valeurs_possibles_N = np.linspace(0, 100, 1)
+bouton_obs = Button(
+ ax_bouton, f" Ajouter {N_serie} obs ", color="purple")
+
+bouton_obs.on_clicked(update_series)
+
+# switch T / freq
+
+chb_enT = CheckButtons(ax_chb, [ "axe des x en fréquence (de 0 à 1)"], [False])
+"""
+for rect in chb_enT.rectangles:
+ rect.set_width(0.1)
+ rect.set_height(1)
+ rect.xy = (0.05,0)
+"""
+
+chb_enT.on_clicked(switch_freq)
+
+
+ax_T.legend(title="formules", bbox_to_anchor=(1.05, 0.10), loc='lower right')
+#plt.tight_layout()
+#fig_pp.savefig("ma_figure.png")
+#fig_pp.canvas.draw()
+plt.show()
+#plt.ioff()
--
GitLab