[SCM] add permutation which correspond to the sum of rainfall & snowfall. add...

[SCM] add permutation which correspond to the sum of rainfall & snowfall. add test based on Durand paper.

experiment/meteo_france_SCM_study/safran/safran.py

@@ -4,13 +4,14 @@ from experiment.meteo_france_SCM_study.abstract_extended_study import AbstractEx
 from experiment.meteo_france_SCM_study.abstract_study import AbstractStudy
 from experiment.meteo_france_SCM_study.abstract_variable import AbstractVariable
 from experiment.meteo_france_SCM_study.safran.safran_variable import SafranSnowfallVariable, \
-    SafranPrecipitationVariable, SafranTemperatureVariable
+    SafranRainfallVariable, SafranTemperatureVariable, SafranTotalPrecipVariable
 
 
 class Safran(AbstractStudy):
 
     def __init__(self, variable_class: type, *args, **kwargs):
-        assert variable_class in [SafranSnowfallVariable, SafranPrecipitationVariable, SafranTemperatureVariable]
+        assert variable_class in [SafranSnowfallVariable, SafranRainfallVariable, SafranTemperatureVariable,
+                                  SafranTotalPrecipVariable]
         super().__init__(variable_class, *args, **kwargs)
         self.model_name = 'Safran'
 
@@ -42,10 +43,16 @@ class ExtendedSafranSnowfall(AbstractExtendedStudy, SafranSnowfall):
     pass
 
 
-class SafranPrecipitation(SafranFrequency):
+class SafranRainfall(SafranFrequency):
 
     def __init__(self, *args, **kwargs):
-        super().__init__(SafranPrecipitationVariable, *args, **kwargs)
+        super().__init__(SafranRainfallVariable, *args, **kwargs)
+
+
+class SafranTotalPrecip(SafranFrequency):
+
+    def __init__(self, *args, **kwargs):
+        super().__init__(SafranTotalPrecipVariable, *args, **kwargs)
 
 
 class SafranTemperature(Safran):

experiment/meteo_france_SCM_study/safran/safran_variable.py

@@ -28,6 +28,15 @@ class SafranSnowfallVariable(AbstractVariable):
         self.nb_consecutive_days_of_snowfall = nb_consecutive_days_of_snowfall
         # Compute the daily snowfall in kg/m2
         snowfall_rates = np.array(dataset.variables[keyword])
+
+        # Compute the mean snowrate, then multiply it by 60 * 60 * 24
+        # day_duration_in_seconds = 24 * 60 * 60
+        # nb_days = len(snowfall_rates) // 24
+        # print(nb_days)
+        # daily_snowrate = [np.mean(snowfall_rates[24 * i:24 * (i + 1) + 1], axis=0) for i in range(nb_days)]
+        # self.daily_snowfall = day_duration_in_seconds * np.array(daily_snowrate)
+
+        # Compute the hourly snowfall first, then aggregate
         mean_snowfall_rates = 0.5 * (snowfall_rates[:-1] + snowfall_rates[1:])
         hourly_snowfall = 60 * 60 * mean_snowfall_rates
         # Transform the snowfall amount into a dataframe
@@ -40,17 +49,30 @@ class SafranSnowfallVariable(AbstractVariable):
         shifted_list = [self.daily_snowfall[i:] for i in range(self.nb_consecutive_days_of_snowfall)]
         # First element of shifted_list is of length n, Second element of length n-1, Third element n-2....
         # The zip is done with respect to the shortest list
-        snowfall_in_consecutive_days = [sum(e) for e in zip(*shifted_list)]
+        snowfall_in_consecutive_days = np.array([sum(e) for e in zip(*shifted_list)])
         # The returned array is of size n-nb_days+1 x nb_massif
-        return np.array(snowfall_in_consecutive_days)
+        return snowfall_in_consecutive_days
 
 
-class SafranPrecipitationVariable(SafranSnowfallVariable):
+class SafranRainfallVariable(SafranSnowfallVariable):
 
     def __init__(self, dataset, altitude, nb_consecutive_days_of_snowfall=1, keyword='Rainf'):
         super().__init__(dataset, altitude, nb_consecutive_days_of_snowfall, keyword)
 
 
+class SafranTotalPrecipVariable(AbstractVariable):
+
+    def __init__(self, dataset, altitude):
+        super().__init__(dataset, altitude)
+        self.snow_precipitation = SafranSnowfallVariable(dataset=dataset, altitude=altitude)
+        self.rain_precipitation = SafranRainfallVariable(dataset=dataset, altitude=altitude)
+
+    @property
+    def daily_time_serie_array(self) -> np.ndarray:
+        return self.snow_precipitation.daily_time_serie_array + self.rain_precipitation.daily_time_serie_array
+
+
+
 class SafranTemperatureVariable(AbstractVariable):
 
     def __init__(self, dataset, altitude, keyword='Tair'):

experiment/meteo_france_SCM_study/visualization/study_visualization/main_study_visualizer.py

 from experiment.meteo_france_SCM_study.abstract_study import AbstractStudy
 from experiment.meteo_france_SCM_study.crocus.crocus import CrocusDepth, CrocusSwe, ExtendedCrocusDepth, \
     ExtendedCrocusSwe
-from experiment.meteo_france_SCM_study.safran.safran import SafranSnowfall, ExtendedSafranSnowfall, SafranPrecipitation, \
-    SafranTemperature
+from experiment.meteo_france_SCM_study.safran.safran import SafranSnowfall, ExtendedSafranSnowfall, SafranRainfall, \
+    SafranTemperature, SafranTotalPrecip
 
 from experiment.meteo_france_SCM_study.visualization.study_visualization.study_visualizer import StudyVisualizer
 from collections import OrderedDict
@@ -49,7 +49,7 @@ def extended_visualization():
 
 
 def annual_mean_vizu_compare_durand_study():
-    for study_class in [SafranPrecipitation, SafranSnowfall, SafranTemperature][2:]:
+    for study_class in [SafranTotalPrecip, SafranRainfall, SafranSnowfall, SafranTemperature][:1]:
         study = study_class(altitude=1800, year_min=1958, year_max=2002)
         study_visualizer = StudyVisualizer(study)
         study_visualizer.visualize_annual_mean_values()
@@ -59,7 +59,7 @@ def normal_visualization():
     save_to_file = False
     only_first_one = True
     # for study_class in SCM_STUDIES[:1]:
-    for study_class in [SafranPrecipitation, SafranSnowfall, SafranTemperature][1:]:
+    for study_class in [SafranRainfall, SafranSnowfall, SafranTemperature][1:]:
         for study in study_iterator(study_class, only_first_one=only_first_one):
             study_visualizer = StudyVisualizer(study, save_to_file=save_to_file)
             # study_visualizer.visualize_independent_margin_fits(threshold=[None, 20, 40, 60][0])

experiment/meteo_france_SCM_study/visualization/study_visualization/study_visualizer.py

@@ -315,6 +315,7 @@ class StudyVisualizer(object):
         massif_name_to_value = OrderedDict()
         df_annual_total = self.study.df_annual_total
         for massif_id, massif_name in enumerate(self.study.safran_massif_names):
+            # We take the mean over all the annual values
             massif_name_to_value[massif_name] = df_annual_total.loc[:, massif_name].mean()
         self.study.visualize_study(ax=ax, massif_name_to_value=massif_name_to_value, show=self.show)
 

test/test_experiment/test_meteo_france_SCM_study/test_SCM_study.py

@@ -6,7 +6,7 @@ import pandas as pd
 from experiment.meteo_france_SCM_study.crocus.crocus import ExtendedCrocusSwe
 from experiment.meteo_france_SCM_study.visualization.study_visualization.main_study_visualizer import study_iterator
 from experiment.meteo_france_SCM_study.safran.safran import SafranSnowfall, ExtendedSafranSnowfall, SafranTemperature, \
-    SafranPrecipitation
+    SafranRainfall, SafranTotalPrecip
 from experiment.meteo_france_SCM_study.visualization.study_visualization.study_visualizer import StudyVisualizer
 from test.test_utils import load_scm_studies
 
@@ -61,15 +61,16 @@ class TestSCMPrecipitation(TestSCMStudy):
 
     def setUp(self) -> None:
         super().setUp()
-        self.study = SafranPrecipitation(altitude=1800, year_min=1958, year_max=2002)
-
-    # def test_durand(self):
-    #     # Test based on Durand paper
-    #     # Test for the mean temperature between 1958 and 2002
-    #     self.check({
-    #         "Mercantour": 1340,
-    #         'Chablais': 1928,
-    #     })
+        self.study = SafranTotalPrecip(altitude=1800, year_min=1958, year_max=2002)
+
+    def test_durand(self):
+        # Test based on Durand paper
+        # (some small differences probably due to the fact that SAFRAN model has evolved since then)
+        # Test for the mean total precipitation (rainfall + snowfall) between 1958 and 2002
+        self.check({
+            "Mercantour": 1346,
+            'Chablais': 1928,
+        })
 
     def round(self, f):
         return int(f)

test/test_utils.py

@@ -13,8 +13,8 @@ from extreme_estimator.extreme_models.max_stable_model.abstract_max_stable_model
     AbstractMaxStableModelWithCovarianceFunction, CovarianceFunction
 from extreme_estimator.extreme_models.max_stable_model.max_stable_models import Smith, BrownResnick, Schlather, \
     Geometric, ExtremalT, ISchlather
-from experiment.meteo_france_SCM_study.safran.safran import SafranSnowfall, Safran, SafranPrecipitation, \
-    SafranTemperature
+from experiment.meteo_france_SCM_study.safran.safran import SafranSnowfall, Safran, SafranRainfall, \
+    SafranTemperature, SafranTotalPrecip
 from spatio_temporal_dataset.coordinates.spatial_coordinates.alps_station_3D_coordinates import \
     AlpsStation3DCoordinatesWithAnisotropy
 from spatio_temporal_dataset.coordinates.spatial_coordinates.generated_spatial_coordinates import \
@@ -98,8 +98,8 @@ def load_test_spatiotemporal_coordinates(nb_points, nb_steps, train_split_ratio=
 def load_safran_studies(altitudes) -> List[Safran]:
     nb_days_list = [1]
     safran_studies = [safran_class(altitude=safran_altitude, nb_consecutive_days=nb_days)
-            for safran_altitude in altitudes for nb_days in nb_days_list
-            for safran_class in [SafranSnowfall, SafranPrecipitation]]
+                      for safran_altitude in altitudes for nb_days in nb_days_list
+                      for safran_class in [SafranSnowfall, SafranRainfall, SafranTotalPrecip]]
     safran_studies += [SafranTemperature(altitude) for altitude in altitudes]
     return safran_studies