[contrasting] add study_visualizer_for_mean_values.py for the new study

extreme_data/meteo_france_data/scm_models_data/visualization/plot_utils.py

+import matplotlib.pyplot as plt
+
+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
+
+
+def load_plot(cmap, graduation, label, massif_name_to_value, altitude, fit_method):
+    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)
+    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
+    ax = plt.gca()
+    AbstractStudy.visualize_study(ax=ax,
+                                  massif_name_to_value=massif_name_to_value,
+                                  massif_name_to_color=massif_name_to_color,
+                                  replace_blue_by_white=True,
+                                  axis_off=False,
+                                  cmap=cmap,
+                                  show_label=False,
+                                  add_colorbar=True,
+                                  show=False,
+                                  vmin=min_ratio,
+                                  vmax=max_ratio,
+                                  ticks_values_and_labels=ticks_values_and_labels,
+                                  label=label,
+                                  fontsize_label=10,
+                                  )
+    ax.get_xaxis().set_visible(True)
+    ax.set_xticks([])
+    ax.set_xlabel('Altitude = {}m'.format(altitude), fontsize=15)
+    ax.set_title('Fit method is {}'.format(fit_method))
\ No newline at end of file

extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py

@@ -12,7 +12,9 @@ import numpy as np
 import pandas as pd
 import seaborn as sns
 
+from extreme_data.meteo_france_data.scm_models_data.visualization.plot_utils import load_plot
 from extreme_fit.estimator.margin_estimator.utils import fitted_stationary_gev
+from extreme_fit.model.margin_model.utils import fitmethod_to_str
 from extreme_fit.model.result_from_model_fit.result_from_extremes.eurocode_return_level_uncertainties import \
     EurocodeConfidenceIntervalFromExtremes, compute_eurocode_confidence_interval
 from extreme_data.meteo_france_data.scm_models_data.abstract_extended_study import AbstractExtendedStudy
@@ -320,7 +322,6 @@ class StudyVisualizer(VisualizationParameters):
             zip([massif_name for _, massif_name in massif_ids_and_names], altitudes_and_res))
         return massif_name_to_eurocode_return_level_uncertainty
 
-
     def visualize_all_mean_and_max_graphs(self):
         specified_massif_ids = [self.study.study_massif_names.index(massif_name)
                                 for massif_name in
@@ -409,7 +410,8 @@ class StudyVisualizer(VisualizationParameters):
         tuples_x_y = [(year, np.mean(data[:, massif_id])) for year, data in
                       self.study.year_to_daily_time_serie_array.items()]
         x, y = list(zip(*tuples_x_y))
-        x, y = self.average_smoothing_with_sliding_window(x, y, window_size_for_smoothing=self.window_size_for_smoothing)
+        x, y = self.average_smoothing_with_sliding_window(x, y,
+                                                          window_size_for_smoothing=self.window_size_for_smoothing)
         ax.plot(x, y, color=color_mean)
         ax.set_ylabel('mean'.format(self.window_size_for_smoothing), color=color_mean)
         massif_name = self.study.study_massif_names[massif_id]
@@ -423,7 +425,8 @@ class StudyVisualizer(VisualizationParameters):
             tuples_x_y = [(year, annual_maxima[massif_id]) for year, annual_maxima in
                           self.study.year_to_annual_maxima.items()]
             x, y = list(zip(*tuples_x_y))
-            x, y = self.average_smoothing_with_sliding_window(x, y, window_size_for_smoothing=self.window_size_for_smoothing)
+            x, y = self.average_smoothing_with_sliding_window(x, y,
+                                                              window_size_for_smoothing=self.window_size_for_smoothing)
             self.massif_id_to_smooth_maxima[massif_id] = (x, y)
         return self.massif_id_to_smooth_maxima[massif_id]
 
@@ -701,3 +704,21 @@ class StudyVisualizer(VisualizationParameters):
             x = np.array(x)[nb_to_delete:-nb_to_delete]
         assert len(x) == len(y), "{} vs {}".format(len(x), len(y))
         return x, y
+
+    # PLot functions that should be common
+
+    def plot_map(self, cmap, fit_method, graduation, label, massif_name_to_value, plot_name):
+        load_plot(cmap, graduation, label, massif_name_to_value, self.study.altitude, fitmethod_to_str(fit_method))
+        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):
+        # 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)
+

extreme_fit/distribution/gev/gev_params.py

@@ -3,7 +3,7 @@ import matplotlib.pyplot as plt
 from typing import List
 
 from cached_property import cached_property
-from mpmath import euler, pi
+from mpmath import euler
 
 from extreme_fit.distribution.abstract_extreme_params import AbstractExtremeParams
 from extreme_fit.distribution.abstract_params import AbstractParams
@@ -86,13 +86,15 @@ class GevParams(AbstractExtremeParams):
     @property
     def mean(self) -> float:
         if self.has_undefined_parameters:
-            return np.nan
+            mean = np.nan
         elif self.shape >= 1:
-            return np.inf
+            mean = np.inf
         elif self.shape == 0:
-            return self.location + self.scale * euler
+            mean = self.location + self.scale * float(euler)
         else:
-            return self.location + self.scale * (self.g(k=1) - 1) / self.shape
+            mean = self.location + self.scale * (self.g(k=1) - 1) / self.shape
+        assert isinstance(mean, float)
+        return mean
 
     @property
     def variance(self) -> float:
@@ -101,7 +103,7 @@ class GevParams(AbstractExtremeParams):
         elif self.shape >= 0.5:
             return np.inf
         elif self.shape == 0.0:
-            return (self.scale * pi) ** 2 / 6
+            return (self.scale * np.pi) ** 2 / 6
         else:
             return ((self.scale / self.shape) ** 2) * (self.g(k=2) - self.g(k=1) ** 2)
 

extreme_fit/model/margin_model/utils.py

@@ -18,5 +18,6 @@ FEVD_MARGIN_FIT_METHOD_TO_ARGUMENT_STR = {
 }
 FEVD_MARGIN_FIT_METHODS = set(FEVD_MARGIN_FIT_METHOD_TO_ARGUMENT_STR.keys())
 
+
 def fitmethod_to_str(fit_method):
-    return str(fit_method).split('.')[-1]
+    return ' '.join(str(fit_method).split('.')[-1].split('_'))

extreme_trend/abstract_gev_trend_test.py

@@ -46,7 +46,7 @@ class AbstractGevTrendTest(object):
     @cached_property
     def unconstrained_estimator(self):
         return fitted_linear_margin_estimator(self.unconstrained_model_class, self.coordinates, self.dataset,
-                                                  self.starting_year, self.fit_method)
+                                              self.starting_year, self.fit_method)
 
     # Likelihood ratio test
 
@@ -108,9 +108,7 @@ class AbstractGevTrendTest(object):
     def relative_change_in_return_level(self, initial_year, final_year):
         return_level_values = []
         for year in [initial_year, final_year]:
-            gev_params = self.unconstrained_estimator.function_from_fit.get_params(
-                coordinate=np.array([year]),
-                is_transformed=False)
+            gev_params = self.get_unconstrained_gev_params(year)
             return_level_values.append(gev_params.quantile(self.quantile_level))
         change_until_final_year = self.time_derivative_times_years(nb_years=final_year - initial_year)
         change_in_between = (return_level_values[1] - return_level_values[0])
@@ -163,8 +161,6 @@ class AbstractGevTrendTest(object):
         ax.plot(standard_gumbel_quantiles, sorted_maxima, linestyle='None',
                 label='Empirical model', marker='o', color='black')
 
-
-
         ax_twiny = ax.twiny()
         return_periods = [10, 25, 50]
         quantiles = self.get_standard_quantiles_for_return_periods(return_periods, psnow)
@@ -182,11 +178,14 @@ class AbstractGevTrendTest(object):
         end_real_proba = 1 - (0.02 / psnow)
         stationary = True
         if stationary:
-            self.plot_model(ax, None, end_proba=end_real_proba, label='Selected model\nwhich is ${}$'.format(self.label),
+            self.plot_model(ax, None, end_proba=end_real_proba,
+                            label='Selected model\nwhich is ${}$'.format(self.label),
                             color='grey')
         else:
-            self.plot_model(ax, 1959, end_proba=end_real_proba, label='Selected model, which is ${}$'.format(self.label))
-            self.plot_model(ax, 2019, end_proba=end_real_proba, label='Selected model, which is ${}$'.format(self.label))
+            self.plot_model(ax, 1959, end_proba=end_real_proba,
+                            label='Selected model, which is ${}$'.format(self.label))
+            self.plot_model(ax, 2019, end_proba=end_real_proba,
+                            label='Selected model, which is ${}$'.format(self.label))
 
         # Plot for the discarded model
         # if 'Verdon' in massif_name and altitude == 300:
@@ -204,9 +203,6 @@ class AbstractGevTrendTest(object):
         #                   + 'with $\zeta_0=0.84$',
         #             color=color)
 
-
-
-
         ax_lim = [-1.5, 4]
         ax.set_xlim(ax_lim)
         ax_twiny.set_xlim(ax_lim)
@@ -231,13 +227,17 @@ class AbstractGevTrendTest(object):
         label = 'Y({})'.format(year) if year is not None else label
         if year is None:
             year = 2019
-        gev_params_year = self.unconstrained_estimator.function_from_fit.get_params(
-                coordinate=np.array([year]),
-                is_transformed=False)
+        gev_params_year = self.get_unconstrained_gev_params(year)
         extended_maxima = [gev_params_year.gumbel_inverse_standardization(q) for q in extended_quantiles]
 
         ax.plot(extended_quantiles, extended_maxima, linestyle='-', label=label, color=color, linewidth=5)
 
+    def get_unconstrained_gev_params(self, year):
+        gev_params_year = self.unconstrained_estimator.function_from_fit.get_params(
+            coordinate=np.array([year]),
+            is_transformed=False)
+        return gev_params_year
+
     def linear_extension(self, ax, q, quantiles, sorted_maxima):
         # Extend the curve linear a bit if the return period 50 is not in the quantiles
         def compute_slope_intercept(x, y):
@@ -368,9 +368,22 @@ class AbstractGevTrendTest(object):
         plt.gca().set_ylim(bottom=0)
 
     def get_gev_params_with_big_shape_and_correct_shape(self):
-        gev_params = self.unconstrained_estimator.function_from_fit.get_params(coordinate=np.array([YEAR_OF_INTEREST_FOR_RETURN_LEVEL]),
-                                                                               is_transformed=False)  # type: GevParams
+        gev_params = self.unconstrained_estimator.function_from_fit.get_params(
+            coordinate=np.array([YEAR_OF_INTEREST_FOR_RETURN_LEVEL]),
+            is_transformed=False)  # type: GevParams
         gev_params_with_corrected_shape = GevParams(loc=gev_params.location,
                                                     scale=gev_params.scale,
                                                     shape=0.5)
         return gev_params, gev_params_with_corrected_shape
+
+    # Mean value
+
+    @property
+    def unconstrained_average_mean_value(self) -> float:
+        mean_values = []
+        for year in self.years:
+            mean = self.get_unconstrained_gev_params(year).mean
+            mean_values.append(mean)
+        average_mean_value = np.mean(mean_values)
+        assert isinstance(average_mean_value, float)
+        return average_mean_value

projects/altitude_spatial_model/altitudes_fit/altitudes_studies.py

@@ -159,7 +159,8 @@ class AltitudesStudies(object):
                 altitudes.append(altitude)
         self.plot_against_altitude(altitudes, ax, massif_id, massif_name, mean_moment)
 
-    def plot_against_altitude(self, altitudes, ax, massif_id, massif_name, mean_moment):
+    @staticmethod
+    def plot_against_altitude(altitudes, ax, massif_id, massif_name, mean_moment):
         di = massif_id // 8
         if di == 0:
             linestyle = '-'

projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/main_snowfall_article.py

+
+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.snowfall_versus_time_and_altitude.study_visualizer_for_mean_values import \
+    StudyVisualizerForMeanValues
+from projects.contrasting_trends_in_snow_loads.snowfall_versus_time_and_altitude.validation_plot import validation_plot
+from projects.exceeding_snow_loads.section_results.plot_trend_curves import plot_trend_map
+
+mpl.rcParams['text.usetex'] = True
+mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
+
+from extreme_data.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal
+from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import \
+    ConfidenceIntervalMethodFromExtremes
+from extreme_trend.visualizers.study_visualizer_for_non_stationary_trends import \
+    StudyVisualizerForNonStationaryTrends
+from extreme_trend.visualizers.utils import load_altitude_to_visualizer
+from projects.exceeding_snow_loads.section_results.plot_uncertainty_curves import plot_uncertainty_massifs
+from projects.exceeding_snow_loads.utils import paper_study_classes, paper_altitudes
+from root_utils import NB_CORES
+
+
+import matplotlib.pyplot as plt
+
+
+def minor_result(altitude):
+    """Plot trends for a single altitude to be fast"""
+    visualizer = StudyVisualizerForNonStationaryTrends(CrocusSnowLoadTotal(altitude=altitude), multiprocessing=True,
+                                                       )
+    visualizer.plot_trends()
+    plt.show()
+
+
+def compute_minimized_aic(visualizer):
+    _ = visualizer.massif_name_to_trend_test_that_minimized_aic
+    return True
+
+
+def intermediate_result(altitudes, massif_names=None,
+                        model_subsets_for_uncertainty=None, uncertainty_methods=None,
+                        study_class=SafranSnowfall1Day,
+                        multiprocessing=False):
+    """
+    Plot all the trends for all altitudes
+    And enable to plot uncertainty plot for some specific massif_names, uncertainty methods to be fast
+    :param altitudes:
+    :param massif_names:
+    :param non_stationary_uncertainty:
+    :param uncertainty_methods:
+    :param study_class:
+    :return:
+    """
+    # Load altitude to visualizer
+    altitude_to_visualizer = load_altitude_to_visualizer(altitudes, massif_names, model_subsets_for_uncertainty,
+                                                         study_class, uncertainty_methods,
+                                                         study_visualizer_class=StudyVisualizerForMeanValues)
+    # Load variable object efficiently
+    for v in altitude_to_visualizer.values():
+        _ = v.study.year_to_variable_object
+    # Compute minimized value efficiently
+    visualizers = list(altitude_to_visualizer.values())
+    if multiprocessing:
+        with Pool(NB_CORES) as p:
+            _ = p.map(compute_minimized_aic, visualizers)
+    else:
+        for visualizer in visualizers:
+            _ = compute_minimized_aic(visualizer)
+
+    # Plots
+    validation_plot(altitude_to_visualizer)
+
+
+def major_result():
+    uncertainty_methods = [ConfidenceIntervalMethodFromExtremes.ci_mle][:]
+    massif_names = ['Beaufortain', 'Vercors']
+    massif_names = None
+    study_classes = [SafranSnowfall1Day]
+    model_subsets_for_uncertainty = None
+    altitudes = paper_altitudes
+    # altitudes = [900, 1200]
+
+    for study_class in study_classes:
+        intermediate_result(altitudes, massif_names, model_subsets_for_uncertainty,
+                            uncertainty_methods, study_class)
+
+
+if __name__ == '__main__':
+    major_result()
+    # intermediate_result(altitudes=[300], massif_names=None,
+    #                     uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
+    #                                          ConfidenceIntervalMethodFromExtremes.ci_mle][1:],
+    #                     multiprocessing=True)

projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/study_visualizer_for_mean_values.py

+import matplotlib.pyplot as plt
+import numpy as np
+from cached_property import cached_property
+
+from extreme_data.eurocode_data.utils import YEAR_OF_INTEREST_FOR_RETURN_LEVEL
+from extreme_data.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
+from extreme_fit.model.margin_model.utils import MarginFitMethod
+from extreme_trend.visualizers.study_visualizer_for_non_stationary_trends import StudyVisualizerForNonStationaryTrends
+
+
+class StudyVisualizerForMeanValues(StudyVisualizerForNonStationaryTrends):
+
+    def __init__(self, study: AbstractStudy, show=True, save_to_file=False, only_one_graph=False, only_first_row=False,
+                 vertical_kde_plot=False, year_for_kde_plot=None, plot_block_maxima_quantiles=False,
+                 temporal_non_stationarity=False, transformation_class=None, verbose=False, multiprocessing=False,
+                 complete_non_stationary_trend_analysis=False, normalization_under_one_observations=True,
+                 uncertainty_methods=None, model_subsets_for_uncertainty=None, uncertainty_massif_names=None,
+                 effective_temporal_covariate=YEAR_OF_INTEREST_FOR_RETURN_LEVEL, relative_change_trend_plot=True,
+                 non_stationary_trend_test_to_marker=None, fit_method=MarginFitMethod.extremes_fevd_mle,
+                 select_only_acceptable_shape_parameter=True, fit_gev_only_on_non_null_maxima=False,
+                 fit_only_time_series_with_ninety_percent_of_non_null_values=True):
+        super().__init__(study, show, save_to_file, only_one_graph, only_first_row, vertical_kde_plot,
+                         year_for_kde_plot, plot_block_maxima_quantiles, temporal_non_stationarity,
+                         transformation_class, verbose, multiprocessing, complete_non_stationary_trend_analysis,
+                         normalization_under_one_observations, uncertainty_methods, model_subsets_for_uncertainty,
+                         uncertainty_massif_names, effective_temporal_covariate, relative_change_trend_plot,
+                         non_stationary_trend_test_to_marker, fit_method, select_only_acceptable_shape_parameter,
+                         fit_gev_only_on_non_null_maxima, fit_only_time_series_with_ninety_percent_of_non_null_values)
+
+    @cached_property
+    def massif_name_to_empirical_mean(self):
+        massif_name_to_empirical_value = {}
+        for massif_name, maxima in self.study.massif_name_to_annual_maxima.items():
+            massif_name_to_empirical_value[massif_name] = np.mean(maxima)
+        return massif_name_to_empirical_value
+
+    @cached_property
+    def massif_name_to_parametric_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
+            print(type(parameter_value), parameter_value)
+            massif_name_to_parameter_value[massif_name] = parameter_value
+        return massif_name_to_parameter_value
+
+    @cached_property
+    def massif_name_to_relative_difference_for_mean(self):
+        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]
+            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):
+        super().plot_abstract(massif_name_to_value, label, label, self.fit_method, graduation, cmap)
+
+

projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/validation_plot.py

+from typing import Dict
+
+import matplotlib.pyplot as plt
+
+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 \
+    StudyVisualizerForReturnLevelChange
+
+
+def validation_plot(altitude_to_visualizer: Dict[int, StudyVisualizerForMeanValues]):
+    # Plot the mean empirical, the mean parametric and the relative difference between the two
+    altitudes = list(altitude_to_visualizer.keys())
+    study_visualizer = list(altitude_to_visualizer.values())[0]
+    altitude_to_relative_differences = {}
+    # Plot map for the repartition of the difference
+    for altitude, visualizer in altitude_to_visualizer.items():
+        altitude_to_relative_differences[altitude] = plot_relative_difference_map(visualizer)
+        study_visualizer.show_or_save_to_file(add_classic_title=False, dpi=500)
+    # Shoe plot with respect to the altitude.
+    plot_shoe_relative_differences_distribution(altitude_to_relative_differences, altitudes)
+    study_visualizer.show_or_save_to_file(add_classic_title=False, dpi=500)
+
+
+def plot_shoe_relative_differences_distribution(altitude_to_relative_differences, altitudes):
+    ax = plt.gca()
+    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 (%)'
+    ax.set_ylabel(ylabel)
+    ax.set_xlabel('Altitude (m)')
+    ax.legend()
+    ax.grid()
+    ax.set_ylim([-8, 5])
+
+
+def plot_relative_difference_map(visualizer: StudyVisualizerForMeanValues):
+    label = ' mean annual maxima of {} ({})'.format('', '')
+    visualizer.plot_abstract_fast(massif_name_to_value=visualizer.massif_name_to_empirical_mean,
+                                  label='Empirical' + label)
+
+    print(visualizer.massif_name_to_parametric_mean)
+    visualizer.plot_abstract_fast(massif_name_to_value=visualizer.massif_name_to_parametric_mean,
+                                  label='Model' + label)
+    print(visualizer.massif_name_to_relative_difference_for_mean)
+    visualizer.plot_abstract_fast(massif_name_to_value=visualizer.massif_name_to_relative_difference_for_mean,
+                                  label='Relative difference of' + label)
+    return list(visualizer.massif_name_to_relative_difference_for_mean.values())

projects/ogorman/gorman_figures/figure1/study_visualizer_for_double_stationary_fit.py

@@ -9,8 +9,7 @@ from extreme_fit.model.margin_model.utils import \
 
 import matplotlib.pyplot as plt
 
-from projects.ogorman.gorman_figures.figure1 import \
-    ResultFromDoubleStationaryFit
+from projects.ogorman.gorman_figures.figure1.result_from_stationary_fit import ResultFromDoubleStationaryFit
 
 
 class StudyVisualizerForReturnLevelChange(StudyVisualizer):
@@ -79,6 +78,7 @@ class StudyVisualizerForReturnLevelChange(StudyVisualizer):
             self.plot_abstract(
                 massif_name_to_value=result.massif_name_to_difference_return_level_and_maxima,
                 label=label, plot_name=plot_name,
+                fit_method=self.fit_method,
                 cmap=plt.cm.coolwarm)
 
     def plot_shape_values(self):
@@ -110,44 +110,5 @@ class StudyVisualizerForReturnLevelChange(StudyVisualizer):
             plot_name = 'Change {} in return levels'.format(b)
             self.plot_abstract(
                 massif_name_to_value=massif_name_to_value,
-                label=label, plot_name=plot_name)
-
-    def plot_abstract(self, massif_name_to_value, label, plot_name, graduation=10.0, cmap=plt.cm.bwr):
-        plot_name1 = '{}/{}'.format(self.study.altitude, plot_name)
-        plot_name2 = '{}/{}'.format(plot_name.split()[0], plot_name)
-        for plot_name in [plot_name1, plot_name2]:
-            self.load_plot(cmap, graduation, label, massif_name_to_value)
-            self.plot_name = plot_name
-            self.show_or_save_to_file(add_classic_title=False, tight_layout=True, no_title=True,
-                                      dpi=500)
-            plt.close()
-
-    def load_plot(self, cmap, graduation, label, massif_name_to_value):
-        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)
-        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
-        ax = plt.gca()
-        AbstractStudy.visualize_study(ax=ax,
-                                      massif_name_to_value=massif_name_to_value,
-                                      massif_name_to_color=massif_name_to_color,
-                                      replace_blue_by_white=True,
-                                      axis_off=False,
-                                      cmap=cmap,
-                                      show_label=False,
-                                      add_colorbar=True,
-                                      show=False,
-                                      vmin=min_ratio,
-                                      vmax=max_ratio,
-                                      ticks_values_and_labels=ticks_values_and_labels,
-                                      label=label,
-                                      fontsize_label=10,
-                                      )
-        ax.get_xaxis().set_visible(True)
-        ax.set_xticks([])
-        ax.set_xlabel('Altitude = {}m'.format(self.study.altitude), fontsize=15)
-        ax.set_title('Fit method is {}'.format(fitmethod_to_str(self.fit_method)))
+                label=label, plot_name=plot_name,
+                fit_method=self.fit_method)