[contrasting] improve display for the snowfall article

extreme_data/meteo_france_data/scm_models_data/visualization/create_shifted_cmap.py

@@ -4,6 +4,7 @@ import matplotlib.cm as cm
 import matplotlib.colorbar as cbar
 import matplotlib.pyplot as plt
 import numpy as np
+from matplotlib.colors import LinearSegmentedColormap
 from mpl_toolkits.axes_grid1 import make_axes_locatable
 
 from extreme_fit.distribution.abstract_params import AbstractParams
@@ -24,6 +25,13 @@ def get_shifted_map(vmin, vmax, cmap=plt.cm.bwr):
     return shifted_cmap
 
 
+def get_half_colormap(cmap):
+    colors = cmap(np.linspace(0.5, 1, cmap.N // 2))
+    # Create a new colormap from those colors
+    cmap2 = LinearSegmentedColormap.from_list('Upper Half', colors)
+    return cmap2
+
+
 def create_colorbase_axis(ax, label, cmap, norm, ticks_values_and_labels=None, fontsize=15):
     divider = make_axes_locatable(ax)
     cax = divider.append_axes('right', size='5%', pad=0.0)
@@ -78,6 +86,15 @@ def ticks_values_and_labels_for_percentages(graduation, max_abs_change):
     ticks_values = [((t / max_abs_change) + 1) / 2 for t in all_ticks_labels]
     return ticks_values, all_ticks_labels
 
+def ticks_values_and_labels_for_positive_value(graduation, max_abs_change):
+    positive_ticks = []
+    tick = 0
+    while tick < max_abs_change:
+        positive_ticks.append(round(tick, 1))
+        tick += graduation
+    ticks_values = [(t / max_abs_change) for t in positive_ticks]
+    return ticks_values, positive_ticks
+
 
 def ticks_and_labels_centered_on_one(max_ratio, min_ratio):
     """When we compute some ratio of two values.

extreme_data/meteo_france_data/scm_models_data/visualization/plot_utils.py

@@ -3,7 +3,8 @@ import matplotlib.colors as mcolors
 
 from extreme_data.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
 from extreme_data.meteo_france_data.scm_models_data.visualization.create_shifted_cmap import \
-    ticks_values_and_labels_for_percentages, get_shifted_map, get_colors
+    ticks_values_and_labels_for_percentages, get_shifted_map, get_colors, ticks_values_and_labels_for_positive_value, \
+    get_half_colormap
 
 
 def plot_against_altitude(altitudes, ax, massif_id, massif_name, values):
@@ -21,12 +22,20 @@ def plot_against_altitude(altitudes, ax, massif_id, massif_name, values):
     ax.plot(altitudes, values, color=color, linewidth=2, label=massif_name_str, linestyle=linestyle)
 
 
-def load_plot(cmap, graduation, label, massif_name_to_value, altitude, fit_method, add_x_label=True):
+def load_plot(cmap, graduation, label, massif_name_to_value, altitude, fit_method, add_x_label=True,
+              negative_and_positive_values=True):
     max_abs_change = max([abs(e) for e in massif_name_to_value.values()])
-    ticks, labels = ticks_values_and_labels_for_percentages(graduation=graduation, max_abs_change=max_abs_change)
-    min_ratio = -max_abs_change
-    max_ratio = max_abs_change
-    cmap = get_shifted_map(min_ratio, max_ratio, cmap)
+    if negative_and_positive_values:
+        ticks, labels = ticks_values_and_labels_for_percentages(graduation=graduation, max_abs_change=max_abs_change)
+        min_ratio = -max_abs_change
+        max_ratio = max_abs_change
+        cmap = get_shifted_map(min_ratio, max_ratio, cmap)
+    else:
+        ticks, labels = ticks_values_and_labels_for_positive_value(graduation=graduation, max_abs_change=max_abs_change)
+        cmap = get_half_colormap(cmap)
+        min_ratio = 0
+        max_ratio = max_abs_change
+
     massif_name_to_color = {m: get_colors([v], cmap, min_ratio, max_ratio)[0]
                             for m, v in massif_name_to_value.items()}
     ticks_values_and_labels = ticks, labels

extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py

@@ -708,19 +708,21 @@ class StudyVisualizer(VisualizationParameters):
 
     # PLot functions that should be common
 
-    def plot_map(self, cmap, fit_method, graduation, label, massif_name_to_value, plot_name, add_x_label=True):
+    def plot_map(self, cmap, fit_method, graduation, label, massif_name_to_value, plot_name, add_x_label=True,
+                 negative_and_positive_values=True):
         load_plot(cmap, graduation, label, massif_name_to_value, self.study.altitude, fitmethod_to_str(fit_method),
-                  add_x_label=add_x_label)
+                  add_x_label=add_x_label, negative_and_positive_values=negative_and_positive_values)
         self.plot_name = plot_name
         # self.show_or_save_to_file(add_classic_title=False, tight_layout=True, no_title=True, dpi=500)
         self.show_or_save_to_file(add_classic_title=False, no_title=True, dpi=500)
         plt.close()
 
-    def plot_abstract(self, massif_name_to_value, label, plot_name, fit_method='', graduation=10.0, cmap=plt.cm.bwr, add_x_label=True):
+    def plot_abstract(self, massif_name_to_value, label, plot_name, fit_method='', graduation=10.0, cmap=plt.cm.bwr,
+                      add_x_label=True, negative_and_positive_values=True):
         # Regroup the plot by altitudes
         plot_name1 = '{}/{}'.format(self.study.altitude, plot_name)
         # Regroup the plot by type of plot also
         plot_name2 = '{}/{}'.format(plot_name.split()[0], plot_name)
         for plot_name in [plot_name1, plot_name2]:
-            self.plot_map(cmap, fit_method, graduation, label, massif_name_to_value, plot_name, add_x_label)
+            self.plot_map(cmap, fit_method, graduation, label, massif_name_to_value, plot_name, add_x_label, negative_and_positive_values)
 

projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/main_snowfall_article.py

@@ -72,8 +72,8 @@ def intermediate_result(altitudes, massif_names=None,
             _ = compute_minimized_aic(visualizer)
 
     # Plots
-    # validation_plot(altitude_to_visualizer)
-    # plot_snowfall_mean(altitude_to_visualizer)
+    validation_plot(altitude_to_visualizer)
+    plot_snowfall_mean(altitude_to_visualizer)
     plot_snowfall_time_derivative_mean(altitude_to_visualizer)
 
 
@@ -84,7 +84,7 @@ def major_result():
     study_classes = [SafranSnowfall1Day]
     model_subsets_for_uncertainty = None
     altitudes = paper_altitudes
-    altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000]
+    altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900]
     # altitudes = [900, 1200]
 
     for study_class in study_classes:

projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/snowfall_plot.py

@@ -22,43 +22,47 @@ def fit_linear_regression(x, y):
 
 def plot_snowfall_time_derivative_mean(altitude_to_visualizer: Dict[int, StudyVisualizerForMeanValues]):
     visualizer = list(altitude_to_visualizer.values())[0]
+    study = visualizer.study
     # Plot the curve for the evolution of the mean
     massif_name_to_a, massif_name_to_b, massif_name_to_r2_score = plot_mean(altitude_to_visualizer, derivative=True)
     # 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 for every km of elevation ()',
+                                  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 of at 2000 m ()',
+    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 ()',
+    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', 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)
 
 
 def plot_snowfall_mean(altitude_to_visualizer: Dict[int, StudyVisualizerForMeanValues]):
     visualizer = list(altitude_to_visualizer.values())[0]
+    study = visualizer.study
     # Plot the curve for the evolution of the mean
     massif_name_to_a, massif_name_to_b, massif_name_to_r2_score = plot_mean(altitude_to_visualizer, derivative=False)
     # 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 for every km of elevation ()',
-                                  add_x_label=False)
+                                  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 ()',
-                                  add_x_label=False)
+    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', graduation=0.1, add_x_label=False)
+    visualizer.plot_abstract_fast(massif_name_to_r2_score, label='r2 mean', graduation=0.1,
+                                  add_x_label=False, negative_and_positive_values=False)
 
 
 def plot_mean(altitude_to_visualizer: Dict[int, StudyVisualizerForMeanValues], derivative=False):

projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/study_visualizer_for_mean_values.py

@@ -35,7 +35,7 @@ class StudyVisualizerForMeanValues(StudyVisualizerForNonStationaryTrends):
         return massif_name_to_empirical_value
 
     @cached_property
-    def massif_name_to_parametric_mean(self):
+    def massif_name_to_model_mean(self):
         massif_name_to_parameter_value = {}
         for massif_name, trend_test in self.massif_name_to_trend_test_that_minimized_aic.items():
             parameter_value = trend_test.unconstrained_average_mean_value
@@ -47,11 +47,12 @@ class StudyVisualizerForMeanValues(StudyVisualizerForNonStationaryTrends):
         massif_name_to_relative_difference = {}
         for massif_name in self.massif_name_to_trend_test_that_minimized_aic.keys():
             e = self.massif_name_to_empirical_mean[massif_name]
-            p = self.massif_name_to_parametric_mean[massif_name]
+            p = self.massif_name_to_model_mean[massif_name]
             relative_diference = 100 * (p-e) / e
             massif_name_to_relative_difference[massif_name] = relative_diference
         return massif_name_to_relative_difference
 
-    def plot_abstract_fast(self, massif_name_to_value, label, graduation=10.0, cmap=plt.cm.coolwarm, add_x_label=True):
-        super().plot_abstract(massif_name_to_value, label, label, self.fit_method, graduation, cmap, add_x_label)
+    def plot_abstract_fast(self, massif_name_to_value, label, graduation=10.0, cmap=plt.cm.coolwarm, add_x_label=True,
+                           negative_and_positive_values=True):
+        super().plot_abstract(massif_name_to_value, label, label, self.fit_method, graduation, cmap, add_x_label, negative_and_positive_values)
 

projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/validation_plot.py

@@ -2,6 +2,8 @@ from typing import Dict
 
 import matplotlib.pyplot as plt
 
+from extreme_data.meteo_france_data.scm_models_data.visualization.main_study_visualizer import \
+    SCM_STUDY_CLASS_TO_ABBREVIATION
 from projects.contrasting_trends_in_snow_loads.snowfall_versus_time_and_altitude.study_visualizer_for_mean_values import \
     StudyVisualizerForMeanValues
 from projects.ogorman.gorman_figures.figure1.study_visualizer_for_double_stationary_fit import \
@@ -28,7 +30,7 @@ def plot_shoe_relative_differences_distribution(altitude_to_relative_differences
     width = 150
     ax.boxplot([altitude_to_relative_differences[a] for a in altitudes], positions=altitudes, widths=width)
     ax.set_xlim([min(altitudes) - width, max(altitudes) + width])
-    ylabel = 'Relative difference between empirical mean and parametric mean (\%)'
+    ylabel = 'Relative difference of the model mean w.r.t. the empirical mean (\%)'
     ax.set_ylabel(ylabel)
     ax.set_xlabel('Altitude (m)')
     ax.legend()
@@ -36,12 +38,13 @@ def plot_shoe_relative_differences_distribution(altitude_to_relative_differences
 
 
 def plot_relative_difference_map(visualizer: StudyVisualizerForMeanValues):
-    label = ' mean annual maxima of {} ({})'.format('', '')
+    study = visualizer.study
+    label = ' mean annual maxima of {} ({})'.format(SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)], study.variable_unit)
     visualizer.plot_abstract_fast(massif_name_to_value=visualizer.massif_name_to_empirical_mean,
-                                  label='Empirical' + label)
-
-    visualizer.plot_abstract_fast(massif_name_to_value=visualizer.massif_name_to_parametric_mean,
-                                  label='Model' + label)
+                                  label='Empirical' + label, negative_and_positive_values=False)
+    visualizer.plot_abstract_fast(massif_name_to_value=visualizer.massif_name_to_model_mean,
+                                  label='Model' + label, negative_and_positive_values=False)
     visualizer.plot_abstract_fast(massif_name_to_value=visualizer.massif_name_to_relative_difference_for_mean,
-                                  label='Relative difference of' + label, graduation=1)
+                                  label='Relative difference of the model mean w.r.t. the empirical mean \n'
+                                        'for the ' + label, graduation=1)
     return list(visualizer.massif_name_to_relative_difference_for_mean.values())