From 77969f35c430a5c37f2d76b554fc681e96285246 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Mon, 25 Jan 2021 15:57:58 +0100
Subject: [PATCH] [slides] fix adamont studies plot. add some plots for the
 slides.

---
 .../adamont_data/adamont_studies.py           |  5 +-
 .../presentation/accumulation_in_winter.py    | 48 ++++++++-----
 .../main_example_snow_depth_total_plot.py     | 67 +++++++++++++++++++
 .../presentation/statistical_model.py         | 13 ++--
 .../projected_data/main_projection.py         | 14 ++--
 5 files changed, 120 insertions(+), 27 deletions(-)
 create mode 100644 projects/exceeding_snow_loads/presentation/main_example_snow_depth_total_plot.py

diff --git a/extreme_data/meteo_france_data/adamont_data/adamont_studies.py b/extreme_data/meteo_france_data/adamont_data/adamont_studies.py
index 55eea8e9..548d8a3c 100644
--- a/extreme_data/meteo_france_data/adamont_data/adamont_studies.py
+++ b/extreme_data/meteo_france_data/adamont_data/adamont_studies.py
@@ -83,9 +83,12 @@ class AdamontStudies(object):
         ax.xaxis.set_ticks(ticks)
         ax.yaxis.grid()
         ax.set_xlim((min(x), max(x)))
+        # Augment the ylim for the legend
+        ylim_min, ylim_max = ax.get_ylim()
+        ax.set_ylim((ylim_min, ylim_max * 1.5))
         ax.tick_params(axis='both', which='major', labelsize=13)
         handles, labels = ax.get_legend_handles_labels()
-        ax.legend(handles[::-1], labels[::-1], ncol=2)
+        ax.legend(handles[::-1], labels[::-1], ncol=2, prop={'size': 7})
         plot_name = 'Annual maxima of {} in {} at {} m'.format(ADAMONT_STUDY_CLASS_TO_ABBREVIATION[self.study_class],
                                                        massif_name.replace('_', ' '),
                                                         self.study.altitude)
diff --git a/projects/exceeding_snow_loads/presentation/accumulation_in_winter.py b/projects/exceeding_snow_loads/presentation/accumulation_in_winter.py
index 4a06f9f2..d73bc1af 100644
--- a/projects/exceeding_snow_loads/presentation/accumulation_in_winter.py
+++ b/projects/exceeding_snow_loads/presentation/accumulation_in_winter.py
@@ -1,23 +1,37 @@
 import matplotlib.pyplot as plt
-from extreme_data.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal
+from extreme_data.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal, CrocusDepth
 from extreme_data.meteo_france_data.scm_models_data.crocus.crocus_variables import AbstractSnowLoadVariable
 
-study = CrocusSnowLoadTotal(altitude=1800)
-year = 1978
-vercors_idx = study.study_massif_names.index('Vercors')
-daily_time_series_vercors_fory_year_of_interest = study.year_to_daily_time_serie_array[year][:, vercors_idx]
-days = [d[5:] for d in study.year_to_days[year]]
 ax = plt.gca()
-x = list(range(len(days)))
-ax.plot(x, daily_time_series_vercors_fory_year_of_interest, linewidth='5')
-ticks_date = ['{:02d}-01'.format(e) for e in [8, 9, 10, 11, 12] + list(range(1, 8))][::2]
-ticks_index = [days.index(d) for d in ticks_date]
-ticks = [x[t] for t in ticks_index]
-labels = ['/'.join(t.split('-')[::-1]) + '/' + (str(year-1) if i < 3 else str(year)) for i, t in enumerate(ticks_date)]
-ax.grid()
-plt.xticks(ticks=ticks_index, labels=labels)
 fontsize = 20
-ax.set_xlabel("Date".format(year), fontsize=20)
-ax.tick_params(axis='both', which='major', labelsize=15)
-ax.set_ylabel("GSL ({})".format(AbstractSnowLoadVariable.UNIT), fontsize=fontsize)
+altitude = 1800
+studies = [CrocusSnowLoadTotal(altitude=altitude), CrocusDepth(altitude=altitude)]
+colors = ['black', 'grey']
+for i, study in enumerate(studies):
+    color = colors[i]
+    if i == 1:
+        ax = ax.twinx()
+
+    year = 1978
+    vercors_idx = study.study_massif_names.index('Vercors')
+    daily_time_series_vercors_fory_year_of_interest = study.year_to_daily_time_serie_array[year][:, vercors_idx]
+    days = [d[5:] for d in study.year_to_days[year]]
+    x = list(range(len(days)))
+    if i == 0:
+        ylabel = 'ground snow load ({})'.format(AbstractSnowLoadVariable.UNIT)
+    else:
+        ylabel = 'snow depth (m)'
+    ax.set_ylabel(ylabel, fontsize=fontsize)
+    ax.plot(x, daily_time_series_vercors_fory_year_of_interest, linewidth='5', color=color, label=ylabel)
+    ticks_date = ['{:02d}-01'.format(e) for e in [8, 9, 10, 11, 12] + list(range(1, 8))][::2]
+    ticks_index = [days.index(d) for d in ticks_date]
+    ticks = [x[t] for t in ticks_index]
+    labels = ['/'.join(t.split('-')[::-1]) + '/' + (str(year-1) if i < 3 else str(year)) for i, t in enumerate(ticks_date)]
+    ax.tick_params(axis='both', which='major', labelsize=15)
+    ax.set_ylim(bottom=0)
+    if i == 0:
+        ax.grid()
+        plt.xticks(ticks=ticks_index, labels=labels)
+        ax.set_xlabel("Date".format(year), fontsize=fontsize)
+    ax.legend(prop={'size': 14})
 plt.show()
diff --git a/projects/exceeding_snow_loads/presentation/main_example_snow_depth_total_plot.py b/projects/exceeding_snow_loads/presentation/main_example_snow_depth_total_plot.py
new file mode 100644
index 00000000..ebc0d5a8
--- /dev/null
+++ b/projects/exceeding_snow_loads/presentation/main_example_snow_depth_total_plot.py
@@ -0,0 +1,67 @@
+import matplotlib.pyplot as plt
+
+from extreme_data.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal, CrocusDepth
+from extreme_data.meteo_france_data.scm_models_data.visualization.main_study_visualizer import \
+    study_iterator_global, SCM_STUDY_CLASS_TO_ABBREVIATION
+from extreme_data.meteo_france_data.scm_models_data.visualization.study_visualizer import \
+    StudyVisualizer
+from projects.exceeding_snow_loads.utils import dpi_paper1_figure
+
+
+def tuples_for_examples_paper1(examples_for_the_paper=True):
+    if examples_for_the_paper:
+
+        marker_altitude_massif_name_for_paper1 = [
+            # ('magenta', 900, 'Ubaye'),
+            ('darkblue', 1800, 'Vercors'),
+            # ('mediumpurple', 2700, 'Beaufortain'),
+        ]
+    else:
+        marker_altitude_massif_name_for_paper1 = [
+            ('magenta', 600, 'Parpaillon'),
+            # ('darkmagenta', 300, 'Devoluy'),
+            ('mediumpurple', 300, 'Aravis'),
+        ]
+    return marker_altitude_massif_name_for_paper1
+
+
+def max_graph_annual_maxima_poster():
+    """
+    We choose these massif because each represents a different eurocode region
+    we also choose them because they belong to a different climatic area
+    :return:
+    """
+    save_to_file = True
+    study_class = CrocusDepth
+
+    examples_for_the_paper = True
+
+    ax = plt.gca()
+    if examples_for_the_paper:
+        # ax.set_ylim([0, 20])
+        # ax.set_yticks(list(range(0, 21, 2)))
+        linewidth = 5
+    else:
+        linewidth = 3
+
+    marker_altitude_massif_name_for_paper1 = tuples_for_examples_paper1(examples_for_the_paper)
+
+    for color, altitude, massif_name in marker_altitude_massif_name_for_paper1[::-1]:
+        for study in study_iterator_global([study_class], altitudes=[altitude]):
+            study_visualizer = StudyVisualizer(study, save_to_file=save_to_file,
+                                               verbose=True,
+                                               multiprocessing=True)
+            snow_abbreviation = SCM_STUDY_CLASS_TO_ABBREVIATION[study_class]
+            snow_abbreviation = 'snow depth'
+            last_plot = True
+            label = '{} massif at {}m'.format(massif_name, altitude)
+            tight_pad = {'h_pad': 0.2}
+            snow_abbreviation = 'annual maximum of ' + snow_abbreviation
+            study_visualizer.visualize_max_graphs_poster(massif_name, altitude, snow_abbreviation, color, label,
+                                                         last_plot, ax, tight_pad=tight_pad,
+                                                         dpi=dpi_paper1_figure,
+                                                         linewidth=linewidth)
+
+
+if __name__ == '__main__':
+    max_graph_annual_maxima_poster()
diff --git a/projects/exceeding_snow_loads/presentation/statistical_model.py b/projects/exceeding_snow_loads/presentation/statistical_model.py
index 72bdf968..a5aa1576 100644
--- a/projects/exceeding_snow_loads/presentation/statistical_model.py
+++ b/projects/exceeding_snow_loads/presentation/statistical_model.py
@@ -30,11 +30,14 @@ def histogram_for_gev():
     import matplotlib.pyplot as plt
     from extreme_data.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal
     from extreme_data.meteo_france_data.scm_models_data.crocus.crocus_variables import AbstractSnowLoadVariable
+    from extreme_data.meteo_france_data.scm_models_data.crocus.crocus import CrocusDepth
     ax = plt.gca()
-    study = CrocusSnowLoadTotal(altitude=1800)
+    study_class = CrocusDepth
+    study = study_class(altitude=1800)
     s = study.observations_annual_maxima.df_maxima_gev.loc['Vercors']
     x_gev = s.values
     gev_params = fitted_stationary_gev(x_gev)
+    print(gev_params.return_level(return_period=50))
     samples = gev_params.sample(10000)
     nb = 10
     epsilon = 0.0
@@ -46,11 +49,11 @@ def histogram_for_gev():
     # x = np.linspace(0.0, 10, 1000)
     # y = gev_params.density(x)
     # ax.plot(x, y, linewidth=5)
-    ax.set_xlabel('Annual maxima of GSL ({})'.format(AbstractSnowLoadVariable.UNIT), fontsize=15)
+    ax.set_xlabel('Annual maximum of snow depth (m)', fontsize=15)
     ax.set_ylabel('Probability', fontsize=15)
     ax.tick_params(axis='both', which='major', labelsize=15)
     ax.set_yticks([0, 0.1, 0.2, 0.3])
-    ax.set_xlim([0, 10])
+    ax.set_xlim([0, 2.5])
     ax.set_ylim([0, 0.3])
 
 
@@ -71,6 +74,6 @@ def histogram_for_normal():
 
 if __name__ == '__main__':
     # binomial_observation()
-    # histogram_for_gev()
-    histogram_for_normal()
+    histogram_for_gev()
+    # histogram_for_normal()
     plt.show()
diff --git a/projects/projected_snowfall/projected_data/main_projection.py b/projects/projected_snowfall/projected_data/main_projection.py
index 8498017a..71d33c7e 100644
--- a/projects/projected_snowfall/projected_data/main_projection.py
+++ b/projects/projected_snowfall/projected_data/main_projection.py
@@ -18,14 +18,20 @@ from projects.projected_snowfall.comparison_with_scm.comparison_historical_visua
 
 
 def main():
-    fast = True
-    adamont_scenario = AdamontScenario.rcp85
+    fast = None
+    adamont_scenario = AdamontScenario.rcp85_extended
     year_min = 1982 if adamont_scenario is AdamontScenario.rcp85_extended else 2006
     # Set the year_min and year_max for the comparison
-    if fast:
+    if fast is True:
         year_max = [2030][0]
         massif_names = ['Vanoise']
         altitudes = [1800]
+    elif fast is None:
+        # year_min = [1951][0]
+        # year_min = [1951][0]
+        year_max = [2005][0]
+        massif_names = ['Vercors']
+        altitudes =  [900, 1200, 1500, 1800, 2100, 2400]
     else:
         year_max = [2100][0]
         massif_names = None
@@ -41,7 +47,7 @@ def main():
                                          altitude=altitude, year_min=year_min,
                                          year_max=year_max, season=season,
                                          scenario=adamont_scenario)
-        adamont_studies.plot_maxima_time_series(massif_names)
+        adamont_studies.plot_maxima_time_series_adamont(massif_names)
 
 
 if __name__ == '__main__':
-- 
GitLab