From 9dd7dc4cfbd76ab1b07c7ea5c087a9a472ee2e52 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Mon, 15 Jun 2020 15:02:30 +0200
Subject: [PATCH] [contrasting] add shape plot

---
 .../main_snowfall_article.py                  |  4 ++-
 .../shape_plot.py                             | 17 ++++++++++
 .../snowfall_plot.py                          | 34 ++++++++++++-------
 .../study_visualizer_for_mean_values.py       | 20 ++++++-----
 .../validation_plot.py                        |  1 -
 5 files changed, 53 insertions(+), 23 deletions(-)
 create mode 100644 projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/shape_plot.py

diff --git a/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/main_snowfall_article.py b/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/main_snowfall_article.py
index 484f59b6..c9908b22 100644
--- a/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/main_snowfall_article.py
+++ b/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/main_snowfall_article.py
@@ -4,6 +4,7 @@ from multiprocessing.pool import Pool
 import matplotlib as mpl
 
 from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day
+from projects.contrasting_trends_in_snow_loads.article2_snowfall_versus_time_and_altitude.shape_plot import shape_plot
 from projects.contrasting_trends_in_snow_loads.article2_snowfall_versus_time_and_altitude.snowfall_plot import \
     plot_snowfall_mean, plot_snowfall_change_mean
 from projects.contrasting_trends_in_snow_loads.article2_snowfall_versus_time_and_altitude.study_visualizer_for_mean_values import \
@@ -78,6 +79,7 @@ def intermediate_result(altitudes, massif_names=None,
     plot_snowfall_mean(altitude_to_visualizer)
     plot_selection_curves(altitude_to_visualizer, paper1=False)
     plot_snowfall_change_mean(altitude_to_visualizer)
+    shape_plot(altitude_to_visualizer)
 
 
 def major_result():
@@ -89,7 +91,7 @@ def major_result():
     altitudes = paper_altitudes
     # altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900]
     # altitudes = [900, 1200, 1500, 1800][:2]
-    # altitudes = [2400, 2700][:2]
+    # altitudes = [1800, 2100, 2400, 2700][:2]
     altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000]
     # altitudes = draft_altitudes
 
diff --git a/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/shape_plot.py b/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/shape_plot.py
new file mode 100644
index 00000000..657bf026
--- /dev/null
+++ b/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/shape_plot.py
@@ -0,0 +1,17 @@
+from typing import Dict
+
+import matplotlib
+import matplotlib.pyplot as plt
+
+from projects.contrasting_trends_in_snow_loads.article2_snowfall_versus_time_and_altitude.study_visualizer_for_mean_values import \
+    StudyVisualizerForMeanValues
+
+
+def shape_plot(altitude_to_visualizer: Dict[int, StudyVisualizerForMeanValues]):
+    # Plot map for the repartition of the difference
+    for altitude, visualizer in altitude_to_visualizer.items():
+        label = ' shape parameter'
+        visualizer.plot_abstract_fast(massif_name_to_value=visualizer.massif_name_to_model_shape_last_year,
+                                      label='Model' + label, negative_and_positive_values=True, add_text=True,
+                                      cmap=matplotlib.cm.get_cmap('BrBG_r'), graduation=0.1)
+    plt.close()
diff --git a/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/snowfall_plot.py b/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/snowfall_plot.py
index 8b9b78d9..5399c78f 100644
--- a/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/snowfall_plot.py
+++ b/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/snowfall_plot.py
@@ -28,21 +28,26 @@ def plot_snowfall_change_mean(altitude_to_visualizer: Dict[int, StudyVisualizerF
     # Augmentation every km
     massif_name_to_augmentation_every_km = {m: a * 1000 for m, a in massif_name_to_a.items()}
     visualizer.plot_abstract_fast(massif_name_to_augmentation_every_km,
-                                  label='Augmentation of time derivative of mean annual maxima of {}\n for every km of elevation ({})'.format(SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)], study.variable_unit),
+                                  label='Augmentation of time derivative of mean annual maxima of {}\n for every km of elevation ({})'.format(
+                                      SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)], study.variable_unit),
                                   add_x_label=False)
     # Value at 2000 m
     massif_name_to_mean_at_2000 = {m: a * 2000 + massif_name_to_b[m] for m, a in massif_name_to_a.items()}
-    visualizer.plot_abstract_fast(massif_name_to_mean_at_2000, label='Time derivative of mean annual maxima \nof {} at 2000 m ({})'.format(SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)], study.variable_unit),
+    visualizer.plot_abstract_fast(massif_name_to_mean_at_2000,
+                                  label='Time derivative of mean annual maxima \nof {} at 2000 m ({})'.format(
+                                      SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)], study.variable_unit),
                                   add_x_label=False)
     # Altitude for the change of dynamic
     massif_name_to_altitude_change_dynamic = {m: - massif_name_to_b[m] / a for m, a in massif_name_to_a.items()}
     # Keep only those that are in a reasonable range
     massif_name_to_altitude_change_dynamic = {m: d for m, d in massif_name_to_altitude_change_dynamic.items()
                                               if 0 < d < 3000}
-    visualizer.plot_abstract_fast(massif_name_to_altitude_change_dynamic, label='Altitude for the change of dynamic (m)',
+    visualizer.plot_abstract_fast(massif_name_to_altitude_change_dynamic,
+                                  label='Altitude for the change of dynamic (m)',
                                   add_x_label=False, graduation=500)
     # R2 score
-    visualizer.plot_abstract_fast(massif_name_to_r2_score, label='r2 time derivative of the mean', graduation=0.1, add_x_label=False,
+    visualizer.plot_abstract_fast(massif_name_to_r2_score, label='r2 time derivative of the mean', graduation=0.1,
+                                  add_x_label=False,
                                   negative_and_positive_values=False)
 
 
@@ -54,11 +59,13 @@ def plot_snowfall_mean(altitude_to_visualizer: Dict[int, StudyVisualizerForMeanV
     # Augmentation every km
     massif_name_to_augmentation_every_km = {m: a * 1000 for m, a in massif_name_to_a.items()}
     visualizer.plot_abstract_fast(massif_name_to_augmentation_every_km,
-                                  label='Augmentation of mean annual maxima of {} \nfor every km of elevation ({})'.format(SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)], study.variable_unit),
+                                  label='Augmentation of mean annual maxima of {} \nfor every km of elevation ({})'.format(
+                                      SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)], study.variable_unit),
                                   add_x_label=False, negative_and_positive_values=False)
     # Value at 2000 m
     massif_name_to_mean_at_2000 = {m: a * 2000 + massif_name_to_b[m] for m, a in massif_name_to_a.items()}
-    visualizer.plot_abstract_fast(massif_name_to_mean_at_2000, label='Mean annual maxima of {} at 2000 m ()'.format(SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)], study.variable_unit),
+    visualizer.plot_abstract_fast(massif_name_to_mean_at_2000, label='Mean annual maxima of {} at 2000 m ()'.format(
+        SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)], study.variable_unit),
                                   add_x_label=False, negative_and_positive_values=False)
     # R2 score
     visualizer.plot_abstract_fast(massif_name_to_r2_score, label='r2 mean', graduation=0.1,
@@ -80,11 +87,12 @@ def plot_mean(altitude_to_visualizer: Dict[int, StudyVisualizerForMeanValues], d
             trend_tests = [altitude_to_visualizer[a].massif_name_to_trend_test_that_minimized_aic[massif_name]
                            for a in altitudes_massif]
             if derivative:
-                moment = 'change in 10 years for significant models'
-                values = [t.change_in_mean_for_the_last_x_years(nb_years=10) for t in trend_tests
-                          if t.is_significant]
-                altitudes_values = [a for a in altitudes_massif
-                             if altitude_to_visualizer[a].massif_name_to_trend_test_that_minimized_aic[massif_name].is_significant]
+                nb_years = 10
+                res = [(a, t.change_in_mean_for_the_last_x_years(nb_years=nb_years))
+                       for i, (a, t) in enumerate(zip(altitudes_massif, trend_tests))
+                       if not t.unconstrained_model_is_stationary]
+                altitudes_values, values = zip(*res)
+                moment = 'change in {} years for significant models'.format(nb_years)
             else:
                 moment = 'mean'
                 values = [t.unconstrained_estimator_gev_params_last_year.mean for t in trend_tests]
@@ -92,7 +100,8 @@ def plot_mean(altitude_to_visualizer: Dict[int, StudyVisualizerForMeanValues], d
             # Plot
             if len(altitudes_values) >= 2:
                 massif_name_to_linear_regression_result[massif_name] = fit_linear_regression(altitudes_values, values)
-                plot_values_against_altitudes(ax, altitudes_values, massif_id, massif_name, moment, study, values, visualizer)
+                plot_values_against_altitudes(ax, altitudes_values, massif_id, massif_name, moment, study, values,
+                                              visualizer)
     ax.legend(prop={'size': 7}, ncol=3)
     visualizer.show_or_save_to_file(dpi=500, add_classic_title=False)
     plt.close()
@@ -111,4 +120,3 @@ def plot_values_against_altitudes(ax, altitudes, massif_id, massif_name, moment,
     # ax.set_ylim([lim_down, lim_up])
     ax.tick_params(axis='both', which='major', labelsize=13)
     visualizer.plot_name = plot_name
-
diff --git a/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/study_visualizer_for_mean_values.py b/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/study_visualizer_for_mean_values.py
index a074aad7..c480742b 100644
--- a/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/study_visualizer_for_mean_values.py
+++ b/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/study_visualizer_for_mean_values.py
@@ -47,13 +47,8 @@ class StudyVisualizerForMeanValues(StudyVisualizerForNonStationaryTrends):
     def massif_name_to_text(self):
         d = {}
         for m, t in self.massif_name_to_trend_test_that_minimized_aic.items():
-            if t.is_significant:
-                name = '$\\textbf{'
-                name += t.name
-                name += '}$'
-            else:
-                name = '${}$'.format(t.name)
-            d[m] = name
+            latex_command = 'textbf' if t.is_significant else 'textrm'
+            d[m] = '$\\' + latex_command + '{' + t.name + '}$'
         return d
 
     # Override the main dict massif_name_to_trend_test_that_minimized_aic
@@ -76,11 +71,20 @@ class StudyVisualizerForMeanValues(StudyVisualizerForNonStationaryTrends):
             massif_name_to_empirical_value[massif_name] = np.mean(maxima)
         return massif_name_to_empirical_value
 
+    @cached_property
+    def massif_name_to_model_shape_last_year(self):
+        massif_name_to_model_value_last_year = {}
+        for massif_name, trend_test in self.massif_name_to_trend_test_that_minimized_aic.items():
+            massif_name_to_model_value_last_year[
+                massif_name] = trend_test.unconstrained_estimator_gev_params_last_year.shape
+        return massif_name_to_model_value_last_year
+
     @cached_property
     def massif_name_to_model_mean_last_year(self):
         massif_name_to_model_value_last_year = {}
         for massif_name, trend_test in self.massif_name_to_trend_test_that_minimized_aic.items():
-            massif_name_to_model_value_last_year[massif_name] = trend_test.unconstrained_estimator_gev_params_last_year.mean
+            massif_name_to_model_value_last_year[
+                massif_name] = trend_test.unconstrained_estimator_gev_params_last_year.mean
         return massif_name_to_model_value_last_year
 
     @cached_property
diff --git a/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/validation_plot.py b/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/validation_plot.py
index 971a7c92..7019880d 100644
--- a/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/validation_plot.py
+++ b/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/validation_plot.py
@@ -61,7 +61,6 @@ def plot_relative_difference_map_order_zero(visualizer: StudyVisualizerForMeanVa
     # visualizer.plot_abstract_fast(massif_name_to_value=visualizer.massif_name_to_relative_difference_for_mean,
     #                               label='Relative difference of the model mean w.r.t. the empirical mean \n'
     #                                     'for the ' + label, graduation=1)
-    visualizer.show_or_save_to_file(add_classic_title=False, dpi=500)
     return list(visualizer.massif_name_to_relative_difference_for_mean.values())
 
 
-- 
GitLab