diff --git a/extreme_data/meteo_france_data/scm_models_data/visualization/main_study_visualizer.py b/extreme_data/meteo_france_data/scm_models_data/visualization/main_study_visualizer.py
index cfc583396558fa2617a60d15657bbf84bf3809ac..271753085c518ba9a2aeb21e86cb751b4c8787ce 100644
--- a/extreme_data/meteo_france_data/scm_models_data/visualization/main_study_visualizer.py
+++ b/extreme_data/meteo_france_data/scm_models_data/visualization/main_study_visualizer.py
@@ -9,7 +9,9 @@ from extreme_data.meteo_france_data.scm_models_data.crocus.crocus import CrocusD
 from extreme_data.meteo_france_data.scm_models_data.crocus.crocus_variables import CrocusDensityVariable
 from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall, ExtendedSafranSnowfall, \
     SafranRainfall, \
-    SafranTemperature, SafranPrecipitation, SafranSnowfall1Day
+    SafranTemperature, SafranPrecipitation, SafranSnowfall1Day, SafranSnowfall3Days, SafranSnowfall5Days, \
+    SafranSnowfall7Days, SafranPrecipitation1Day, SafranPrecipitation3Days, SafranPrecipitation5Days, \
+    SafranPrecipitation7Days
 from extreme_data.meteo_france_data.scm_models_data.visualization.study_visualizer import \
     StudyVisualizer
 from projects.exceeding_snow_loads.section_discussion.crocus_study_comparison_with_eurocode import \
@@ -30,6 +32,13 @@ SCM_STUDY_NAME_TO_SCM_STUDY = dict(zip(SCM_STUDIES_NAMES, SCM_STUDIES))
 SCM_STUDY_CLASS_TO_ABBREVIATION = {
     SafranSnowfall: 'SF3',
     SafranSnowfall1Day: 'SF1',
+    SafranSnowfall3Days: 'SF3',
+    SafranSnowfall5Days: 'SF5',
+    SafranSnowfall7Days: 'SF7',
+    SafranPrecipitation1Day: 'PR1',
+    SafranPrecipitation3Days: 'PR3',
+    SafranPrecipitation5Days: 'PR5',
+    SafranPrecipitation7Days: 'PR7',
     CrocusSweTotal: 'SWE',
     CrocusSwe3Days: 'SWE3',
     CrocusSnowLoadEurocode: 'GSL from annual maximum of HS \nand {}'.format(eurocode_snow_density),
diff --git a/projects/contrasting_trends_in_snow_loads/altitudes_fit/altitudes_studies.py b/projects/contrasting_trends_in_snow_loads/altitudes_fit/altitudes_studies.py
index d78cafe84dcf36e07b5fef0d5822902ff97ec635..84a6c8f640b19e134850a46dcd466b36319f8334 100644
--- a/projects/contrasting_trends_in_snow_loads/altitudes_fit/altitudes_studies.py
+++ b/projects/contrasting_trends_in_snow_loads/altitudes_fit/altitudes_studies.py
@@ -1,4 +1,5 @@
 import pandas as pd
+import numpy as np
 from collections import OrderedDict
 
 from cached_property import cached_property
@@ -17,6 +18,7 @@ from spatio_temporal_dataset.coordinates.temporal_coordinates.generated_temporal
 from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
 from spatio_temporal_dataset.spatio_temporal_observations.annual_maxima_observations import AnnualMaxima
 import matplotlib.pyplot as plt
+import matplotlib.colors as mcolors
 
 
 class AltitudesStudies(object):
@@ -91,13 +93,16 @@ class AltitudesStudies(object):
         study_visualizer.plot_name = plot_name
         study_visualizer.show_or_save_to_file(add_classic_title=False, dpi=500)
 
-    def plot_maxima_time_series(self, massif_names=None, show=False):
+    def run_for_each_massif(self, function, massif_names, **kwargs):
         massif_names = massif_names if massif_names is not None else self.study.all_massif_names()
         assert isinstance(massif_names, list)
-        for massif_name in massif_names:
-            self._plot_maxima_time_series(massif_name, show=show)
+        for i, massif_name in enumerate(massif_names):
+            function(massif_name, massif_id=i, **kwargs)
+
+    def plot_maxima_time_series(self, massif_names=None, show=False):
+        self.run_for_each_massif(self._plot_maxima_time_series, massif_names, show=show)
 
-    def _plot_maxima_time_series(self, massif_name, show=False):
+    def _plot_maxima_time_series(self, massif_name, massif_id, show=False):
         ax = plt.gca()
         x = self.study.ordered_years
         for altitude, study in list(self.altitude_to_study.items())[::-1]:
@@ -113,3 +118,55 @@ class AltitudesStudies(object):
         ax.set_xlabel('years', fontsize=15)
         self.show_or_save_to_file(plot_name=plot_name, show=show)
         ax.clear()
+
+    def plot_mean_maxima_against_altitude(self, massif_names=None, show=False, std=False, change=False):
+        ax = plt.gca()
+        self.run_for_each_massif(self._plot_mean_maxima_against_altitude, massif_names, ax=ax, std=std, change=change)
+        ax.legend(prop={'size': 7}, ncol=3)
+        moment = 'Mean' if not std else 'Std'
+        if change is None:
+            moment += ' relative change for'
+        elif change:
+            moment += ' change for'
+        plot_name = '{} annual maxima of {}'.format(moment, SCM_STUDY_CLASS_TO_ABBREVIATION[self.study_class])
+        ax.set_ylabel('{} ({})'.format(plot_name, self.study.variable_unit), fontsize=15)
+        ax.set_xlabel('altitudes', fontsize=15)
+        lim_down, lim_up = ax.get_ylim()
+        lim_up += (lim_up - lim_down) / 3
+        ax.set_ylim([lim_down, lim_up])
+        ax.tick_params(axis='both', which='major', labelsize=13)
+        self.show_or_save_to_file(plot_name=plot_name, show=show)
+        ax.clear()
+
+    def _plot_mean_maxima_against_altitude(self, massif_name, massif_id, ax=None, std=False, change=False):
+        assert ax is not None
+        altitudes = []
+        mean_moment = []
+        for altitude, study in self.altitude_to_study.items():
+            if massif_name in study.massif_name_to_annual_maxima:
+                annual_maxima = study.massif_name_to_annual_maxima[massif_name]
+                function = np.std if std else np.mean
+                if change in [True, None]:
+                    after = function(annual_maxima[31:])
+                    before = function(annual_maxima[:31])
+                    moment = after - before
+                    if change is None:
+                        moment /= before
+                else:
+                    moment = function(annual_maxima)
+                mean_moment.append(moment)
+                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):
+        di = massif_id // 8
+        if di == 0:
+            linestyle = '-'
+        elif di == 1:
+            linestyle = 'dotted'
+        else:
+            linestyle = '--'
+        colors = list(mcolors.TABLEAU_COLORS)
+        colors[-3:-1] = []  # remove gray and olive
+        color = colors[massif_id % 8]
+        ax.plot(altitudes, mean_moment, color=color, linewidth=2, label=massif_name, linestyle=linestyle)
diff --git a/projects/contrasting_trends_in_snow_loads/altitudes_fit/main_altitudes_studies.py b/projects/contrasting_trends_in_snow_loads/altitudes_fit/main_altitudes_studies.py
index f7af11a9c12b821b0f959795b67ca6de027cc8a2..759b35c8a02b6fa236afff734adcda746389f19a 100644
--- a/projects/contrasting_trends_in_snow_loads/altitudes_fit/main_altitudes_studies.py
+++ b/projects/contrasting_trends_in_snow_loads/altitudes_fit/main_altitudes_studies.py
@@ -1,15 +1,26 @@
-from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day
+from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day, SafranSnowfall3Days, \
+    SafranSnowfall5Days, SafranSnowfall7Days, SafranPrecipitation1Day, SafranPrecipitation3Days, \
+    SafranPrecipitation5Days, SafranPrecipitation7Days
 from projects.contrasting_trends_in_snow_loads.altitudes_fit.altitudes_studies import AltitudesStudies
 
 
-def main_plots():
-    altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000]
-    study_class = SafranSnowfall1Day
-    studies = AltitudesStudies(study_class, altitudes)
-    # massifs_names = ['Vercors']
-    # studies.plot_maxima_time_series(massif_names=massifs_names)
-    studies.plot_maxima_time_series()
+def main_plots_moments():
+    altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000]
+    study_classes = [SafranSnowfall1Day, SafranSnowfall3Days, SafranSnowfall5Days, SafranSnowfall7Days][:2]
+    study_classes = [SafranPrecipitation1Day, SafranPrecipitation3Days, SafranPrecipitation5Days,
+                     SafranPrecipitation7Days][:]
+    study_classes = [SafranSnowfall1Day, SafranPrecipitation1Day]
+
+    for study_class in study_classes:
+
+        studies = AltitudesStudies(study_class, altitudes)
+        # massifs_names = ['Vercors', 'Chartreuse', 'Belledonne']
+        # studies.plot_mean_maxima_against_altitude(massif_names=massifs_names, show=True)
+        # studies.plot_maxima_time_series()
+        for std in [True, False]:
+            for change in [True, False, None]:
+                studies.plot_mean_maxima_against_altitude(std=std, change=change)
 
 
 if __name__ == '__main__':
-    main_plots()
+    main_plots_moments()
diff --git a/projects/contrasting_trends_in_snow_loads/spatial trends/main_constrasting_result.py b/projects/contrasting_trends_in_snow_loads/spatial trends/main_constrasting_result.py
index 407958e6797b5b59dca648e43175de007959cacf..57d3cbea1a0522ac5fed6c52922a3abf0c4b53c4 100644
--- a/projects/contrasting_trends_in_snow_loads/spatial trends/main_constrasting_result.py	
+++ b/projects/contrasting_trends_in_snow_loads/spatial trends/main_constrasting_result.py	
@@ -70,8 +70,7 @@ def major_result():
     altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][:]
     # altitudes = [1200, 1500, 1800, 2100, 2400, 2700][:]
     snow_load_classes = [CrocusSnowLoad1Day, CrocusSnowLoad3Days, CrocusSnowLoad5Days, CrocusSnowLoad7Days][:]
-    precipitation_classes = [SafranPrecipitation1Day, SafranPrecipitation3Days, SafranPrecipitation5Days,
-                             SafranPrecipitation7Days][:]
+    precipitation_classes = [SafranPrecipitation1Day, SafranPrecipitation3Days, SafranPrecipitation5Days, SafranPrecipitation7Days][:]
     snowfall_classes = [SafranSnowfall1Day, SafranSnowfall3Days, SafranSnowfall5Days, SafranSnowfall7Days]
     rainfall_classes = [SafranRainfall1Day, SafranRainfall3Days, SafranRainfall5Days, SafranRainfall7Days]
     study_classes = precipitation_classes + snow_load_classes