From 35620a94b429794376d7abd66994aea6e43af92f Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Thu, 8 Oct 2020 18:04:08 +0200
Subject: [PATCH] [contrasting] improve plots for the first version of the
 paper.

---
 .../altitudes_fit/main_altitudes_studies.py   |  9 +++--
 .../one_fold_analysis/altitude_group.py       | 16 ++++++++
 ...es_visualizer_for_non_stationary_models.py | 30 ++++++++++----
 .../plots/plot_histogram_altitude_studies.py  | 39 +++++++++++++------
 .../preliminary_analysis.py                   |  3 +-
 5 files changed, 73 insertions(+), 24 deletions(-)

diff --git a/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py b/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
index 4b0f672b..64cb1a60 100644
--- a/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
+++ b/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
@@ -55,9 +55,10 @@ def main_loop(altitudes_list, massif_names, seasons, study_classes):
 
 
 def plot_visualizers(massif_names, visualizer_list):
-    plot_histogram_all_trends_against_altitudes(massif_names, visualizer_list)
-    plot_histogram_all_models_against_altitudes(massif_names, visualizer_list)
-    plot_shoe_plot_changes_against_altitude(massif_names, visualizer_list)
+    # plot_histogram_all_trends_against_altitudes(massif_names, visualizer_list)
+    # plot_histogram_all_models_against_altitudes(massif_names, visualizer_list)
+    for relative in [True, False]:
+        plot_shoe_plot_changes_against_altitude(massif_names, visualizer_list, relative=relative)
     # plot_coherence_curves(massif_names, visualizer_list)
     pass
 
@@ -70,7 +71,7 @@ def plot_visualizer(massif_names, visualizer):
     #     for change in [True, False, None]:
     #         studies.plot_mean_maxima_against_altitude(massif_names=massif_names, std=std, change=change)
     # Plot the results for the model that minimizes the individual aic
-    plot_individual_aic(visualizer)
+    # plot_individual_aic(visualizer)
     # Plot the results for the model that minimizes the total aic
     # plot_total_aic(model_classes, visualizer)
     pass
diff --git a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitude_group.py b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitude_group.py
index 4d9a644f..fecd357a 100644
--- a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitude_group.py
+++ b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitude_group.py
@@ -46,6 +46,10 @@ class LowAltitudeGroup(AbstractAltitudeGroup):
     def group_id(self):
         return 1
 
+    @property
+    def graduation_for_return_level(self):
+        return 10
+
     @property
     def name(self):
         return 'low'
@@ -61,6 +65,10 @@ class MidAltitudeGroup(AbstractAltitudeGroup):
     def group_id(self):
         return 2
 
+    @property
+    def graduation_for_return_level(self):
+        return 20
+
     @property
     def name(self):
         return 'mid'
@@ -76,6 +84,10 @@ class HighAltitudeGroup(AbstractAltitudeGroup):
     def group_id(self):
         return 3
 
+    @property
+    def graduation_for_return_level(self):
+        return 40
+
     @property
     def name(self):
         return 'high'
@@ -95,6 +107,10 @@ class VeyHighAltitudeGroup(AbstractAltitudeGroup):
     def name(self):
         return 'very high'
 
+    @property
+    def graduation_for_return_level(self):
+        return 80
+
     @property
     def reference_altitude(self):
         return 3300
diff --git a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
index 6ee246b2..8825fb09 100644
--- a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
+++ b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
@@ -155,11 +155,15 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
         if is_return_level_plot:
             add_colorbar = True
             max_abs_change = None
-            massif_name_to_text = {m: '+- 0' for m in massif_name_to_text.keys()}
+            massif_name_to_text = None
+            graduation = self.altitude_group.graduation_for_return_level
+        else:
+            graduation = 10
 
         negative_and_positive_values = self.moment_names.index(method_name) > 0
         # Plot the map
-        self.plot_map(cmap=plt.cm.coolwarm, fit_method=self.fit_method, graduation=10,
+
+        self.plot_map(cmap=plt.cm.coolwarm, fit_method=self.fit_method, graduation=graduation,
                       label=ylabel, massif_name_to_value=massif_name_to_value,
                       plot_name=plot_name, add_x_label=True,
                       negative_and_positive_values=negative_and_positive_values,
@@ -172,8 +176,8 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
 
     @property
     def add_colorbar(self):
-        # return isinstance(self.altitude_group, (VeyHighAltitudeGroup))
-        return isinstance(self.altitude_group, (VeyHighAltitudeGroup, MidAltitudeGroup))
+        return isinstance(self.altitude_group, (VeyHighAltitudeGroup))
+        # return isinstance(self.altitude_group, (VeyHighAltitudeGroup, MidAltitudeGroup))
 
     def plot_against_years(self, method_name, order):
         ax = plt.gca()
@@ -406,16 +410,26 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
         percents = 100 * nbs / nb_valid_massif_names
         return [nb_valid_massif_names] + list(percents)
 
-    def all_changes(self, massif_names):
+    def all_changes(self, massif_names, relative=False):
         """return percents which contain decrease, significant decrease, increase, significant increase percentages"""
         valid_massif_names = self.get_valid_names(massif_names)
         changes = []
+        non_stationary_changes = []
+        non_stationary_significant_changes = []
         for one_fold in [one_fold for m, one_fold in self.massif_name_to_one_fold_fit.items()
                          if m in valid_massif_names]:
             # Compute changes
-            changes.append(one_fold.change_in_return_level_for_reference_altitude)
-
-        return changes
+            if relative:
+                change = one_fold.relative_change_in_return_level_for_reference_altitude
+            else:
+                change = one_fold.change_in_return_level_for_reference_altitude
+            changes.append(change)
+            if change != 0:
+                non_stationary_changes.append(change)
+                if one_fold.is_significant:
+                    non_stationary_significant_changes.append(change)
+
+        return changes, non_stationary_changes, non_stationary_significant_changes
 
     def get_valid_names(self, massif_names):
         valid_massif_names = set(self.massif_name_to_one_fold_fit.keys())
diff --git a/projects/altitude_spatial_model/altitudes_fit/plots/plot_histogram_altitude_studies.py b/projects/altitude_spatial_model/altitudes_fit/plots/plot_histogram_altitude_studies.py
index 58a56bb3..8c04dd74 100644
--- a/projects/altitude_spatial_model/altitudes_fit/plots/plot_histogram_altitude_studies.py
+++ b/projects/altitude_spatial_model/altitudes_fit/plots/plot_histogram_altitude_studies.py
@@ -4,6 +4,7 @@ from typing import List
 import numpy as np
 
 import matplotlib.pyplot as plt
+from matplotlib.lines import Line2D
 
 from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitudes_studies_visualizer_for_non_stationary_models import \
     AltitudesStudiesVisualizerForNonStationaryModels
@@ -114,10 +115,14 @@ def plot_histogram_all_trends_against_altitudes(massif_names, visualizer_list: L
 
 
 def plot_shoe_plot_changes_against_altitude(massif_names, visualizer_list: List[
-    AltitudesStudiesVisualizerForNonStationaryModels]):
+    AltitudesStudiesVisualizerForNonStationaryModels],
+                                            relative=False):
     visualizer = visualizer_list[0]
 
-    all_changes = [v.all_changes(massif_names) for v in visualizer_list]
+    all_changes = [v.all_changes(massif_names, relative=relative) for v in visualizer_list]
+    all_changes = list(zip(*all_changes))
+    labels = ['All models', 'Non-stationary models', 'Non-stationary and significant models']
+    colors = ['darkgreen', 'forestgreen', 'limegreen']
     nb_massifs = [len(v.get_valid_names(massif_names)) for v in visualizer_list]
 
     plt.close()
@@ -128,27 +133,39 @@ def plot_shoe_plot_changes_against_altitude(massif_names, visualizer_list: List[
     labelsize = 10
     linewidth = 3
 
-    x = np.array([3 * width * (i + 1) for i in range(len(nb_massifs))])
 
-    ax.boxplot(all_changes, positions=x, widths=width)
+    x = np.array([4 * width * (i + 1) for i in range(len(nb_massifs))])
 
-    ax.legend(loc='upper left', prop={'size': size})
-    ax.set_ylabel('Changes between 1969 and 2019 in {}-year return levels ({})'.format(OneFoldFit.return_period,
-                                                                                       visualizer.study.variable_unit),
+    for changes in all_changes:
+        print(changes)
+    for j, (changes, label, color) in enumerate(list(zip(all_changes, labels, colors)), -1):
+        bplot = ax.boxplot(list(changes), positions=x + j * width, widths=width, patch_artist=True)
+        for patch in bplot['boxes']:
+            patch.set_facecolor(color)
+
+    custom_lines = [Line2D([0], [0], color=color, lw=4) for color in colors]
+    ax.legend(custom_lines, labels, loc='lower right')
+
+    start = 'Relative changes' if relative else 'Changes'
+    unit = '\%' if relative else visualizer.study.variable_unit
+    ax.set_ylabel('{} between 1969 and 2019 in {}-year return levels ({})'.format(start, OneFoldFit.return_period,
+                                                                                  unit),
                   fontsize=legend_fontsize)
     ax.set_xlabel('Altitudes', fontsize=legend_fontsize)
     ax.tick_params(axis='both', which='major', labelsize=labelsize)
     ax.set_xticks(x)
     ax.yaxis.grid()
+
     altitudes = [v.altitude_group.reference_altitude for v in visualizer_list]
     ax.set_xticklabels([str(a) for a in altitudes])
 
-    shift = 300
-    ax.set_xlim((min(altitudes) - shift, max(altitudes) + shift))
+    shift = 2 * width
+    ax.set_xlim((min(x) - shift, max(x) + shift))
 
-    plot_nb_massif_on_upper_axis(ax, labelsize, legend_fontsize, nb_massifs, x)
+    # I could display the number of massif used to build each box plot.
+    # plot_nb_massif_on_upper_axis(ax, labelsize, legend_fontsize, nb_massifs, x)
 
-    visualizer.plot_name = 'All changes'
+    visualizer.plot_name = 'All ' + start
     visualizer.show_or_save_to_file(add_classic_title=False, no_title=True)
 
     plt.close()
diff --git a/projects/altitude_spatial_model/preliminary_analysis.py b/projects/altitude_spatial_model/preliminary_analysis.py
index 581cabdb..b2942820 100644
--- a/projects/altitude_spatial_model/preliminary_analysis.py
+++ b/projects/altitude_spatial_model/preliminary_analysis.py
@@ -33,6 +33,7 @@ class PointwiseGevStudyVisualizer(AltitudesStudies):
                     linear_coef, _, r2 = self._plot_gev_params_against_altitude_one_massif(ax, massif_name, param_name)
                     massif_name_to_linear_coef[massif_name] = 100 * linear_coef[0]
                     massif_name_to_r2_score[massif_name] = str(round(r2, 2))
+            print(param_name, np.mean([c for c in massif_name_to_linear_coef.values()]))
 
             # Display x label
             xticks = [1000 * i for i in range(1, 5)]
@@ -252,7 +253,7 @@ if __name__ == '__main__':
     altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300]
     altitudes = [600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900]
     # altitudes = paper_altitudes
-    altitudes = [1800, 2100]
+    # altitudes = [1800, 2100]
     visualizer = PointwiseGevStudyVisualizer(SafranSnowfall1Day, altitudes=altitudes)
     visualizer.plot_gev_params_against_altitude()
     # visualizer.plot_gev_params_against_time_for_all_altitudes()
-- 
GitLab