diff --git a/extreme_data/meteo_france_data/scm_models_data/crocus/crocus_snow_density.py b/extreme_data/meteo_france_data/scm_models_data/crocus/crocus_snow_density.py
new file mode 100644
index 0000000000000000000000000000000000000000..6035221e395eb48ea354242b40325779faa58064
--- /dev/null
+++ b/extreme_data/meteo_france_data/scm_models_data/crocus/crocus_snow_density.py
@@ -0,0 +1,30 @@
+from collections import OrderedDict
+
+from cached_property import cached_property
+
+from extreme_data.meteo_france_data.scm_models_data.crocus.crocus import Crocus, CrocusSweTotal
+from extreme_data.meteo_france_data.scm_models_data.crocus.crocus_variables import CrocusDensityVariable
+
+
+class CrocusSnowDensity(Crocus):
+
+    def __init__(self, *args, **kwargs):
+        super().__init__(CrocusDensityVariable, *args, **kwargs)
+        study_swe = CrocusSweTotal(*args, **kwargs)
+        self._year_to_annual_maxima = OrderedDict()
+        for year in study_swe.ordered_years:
+            daily_time_series_swe = study_swe.year_to_daily_time_serie_array[year]
+            daily_time_series_snow_depth = self.year_to_daily_time_serie_array[year]
+            daily_time_series_density = daily_time_series_swe / daily_time_series_snow_depth
+            ind_to_exclude = daily_time_series_snow_depth < 0.1
+            daily_time_series_density[ind_to_exclude] = 0
+            self._year_to_annual_maxima[year] = daily_time_series_density.max(axis=0)
+
+    @cached_property
+    def year_to_annual_maxima(self) -> OrderedDict:
+        return self._year_to_annual_maxima
+
+
+if __name__ == '__main__':
+    study = CrocusSnowDensity()
+    print(study.year_to_annual_maxima[1959])
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 fb55246b692b1a3dc91b1f671a2532bc0e2162b6..e8d68e3bda9c71da78e1bed1c3358c70b3356e5d 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
@@ -8,6 +8,7 @@ from extreme_data.meteo_france_data.scm_models_data.crocus.crocus import CrocusD
     ExtendedCrocusDepth, \
     ExtendedCrocusSweTotal, CrocusDaysWithSnowOnGround, CrocusSwe3Days, CrocusSnowLoad3Days, CrocusSnowLoadTotal, \
     CrocusSnowLoadEurocode, CrocusSnowLoad5Days, CrocusSnowLoad7Days
+from extreme_data.meteo_france_data.scm_models_data.crocus.crocus_snow_density import CrocusSnowDensity
 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, \
@@ -55,7 +56,8 @@ SCM_STUDY_CLASS_TO_ABBREVIATION = {
     CrocusDifferenceSnowLoad: ('max GSL - GSL from max HS \n with {}'.format(eurocode_snow_density)),
     CrocusSnowDepthDifference: 'max HS - HS at max of GSL',
     CrocusSnowDepthAtMaxofSwe: 'HS at max of GSL',
-    SafranSnowfallSimulationRCP85: 'SF1 RCP85 projections'
+    CrocusSnowDensity: 'Density',
+    SafranSnowfallSimulationRCP85: 'SF1 RCP85 projections',
 }
 
 altitude_massif_name_and_study_class_for_poster = [