[SCM] add crocus comparison to Durand paper

experiment/meteo_france_SCM_study/abstract_study.py

@@ -82,9 +82,12 @@ class AbstractStudy(object):
         # Map each year to an array of size nb_massif
         year_to_annual_mean = OrderedDict()
         for year, time_serie in self._year_to_daily_time_serie_array.items():
-            year_to_annual_mean[year] = self.annual_aggregation_function(time_serie, axis=0)
+            year_to_annual_mean[year] = self.apply_annual_aggregation(time_serie)
         return year_to_annual_mean
 
+    def apply_annual_aggregation(self, time_serie):
+        return self.annual_aggregation_function(time_serie, axis=0)
+
     def instantiate_variable_object(self, dataset) -> AbstractVariable:
         return self.variable_class(dataset, self.altitude)
 
@@ -142,11 +145,14 @@ class AbstractStudy(object):
     """ Visualization methods """
 
     def visualize_study(self, ax=None, massif_name_to_value=None, show=True, fill=True, replace_blue_by_white=True,
-                        label=None):
-        massif_names, values = list(zip(*massif_name_to_value.items()))
-        colors = get_color_rbga_shifted(ax, replace_blue_by_white, values, label=label)
-        massif_name_to_fill_kwargs = {massif_name: {'color': color} for massif_name, color in
-                                      zip(massif_names, colors)}
+                        label=None, add_text=False):
+        if massif_name_to_value is None:
+            massif_name_to_fill_kwargs = None
+        else:
+            massif_names, values = list(zip(*massif_name_to_value.items()))
+            colors = get_color_rbga_shifted(ax, replace_blue_by_white, values, label=label)
+            massif_name_to_fill_kwargs = {massif_name: {'color': color} for massif_name, color in
+                                          zip(massif_names, colors)}
 
         if ax is None:
             ax = plt.gca()
@@ -176,12 +182,13 @@ class AbstractStudy(object):
                 # ax.text(x, y, massif_name)
         # Display the center of the massif
         ax.scatter(self.massifs_coordinates.x_coordinates, self.massifs_coordinates.y_coordinates, s=1)
-        # Display the name of the massif
-        for _, row in self.massifs_coordinates.df_all_coordinates.iterrows():
-            x, y = list(row)
-            massif_name = row.name
-            value = massif_name_to_value[massif_name]
-            ax.text(x, y, str(round(value, 1)))
+        # Display the name or value of the massif
+        if add_text:
+            for _, row in self.massifs_coordinates.df_all_coordinates.iterrows():
+                x, y = list(row)
+                massif_name = row.name
+                value = massif_name_to_value[massif_name]
+                ax.text(x, y, str(round(value, 1)))
 
         if show:
             plt.show()

experiment/meteo_france_SCM_study/abstract_variable.py

@@ -2,6 +2,9 @@ import numpy as np
 
 
 class AbstractVariable(object):
+    """
+    All Variable (CROCUS & SAFRAN) are available since 1958-08-01 06:00:00
+    """
 
     NAME = ''
 

experiment/meteo_france_SCM_study/crocus/crocus.py

@@ -17,12 +17,18 @@ class Crocus(AbstractStudy):
 
     @property
     def variable_name(self):
-        suffix = '' if self.altitude == 2400 else ' average of data observed every 6 hours'
+        suffix = '' if self.altitude == 2400 else ' instantaneous data observed sampled every 24 hours'
         return super().variable_name + suffix
 
     def annual_aggregation_function(self, *args, **kwargs):
         return np.mean(*args, **kwargs)
 
+    def apply_annual_aggregation(self, time_serie):
+        # In the Durand paper, we only want the data from November to April
+        # 91 = 30 + 31 + 30 first days of the time serie correspond to the month of August + September + October
+        # 92 = 31 + 30 + 31 last days correspond to the month of May + June + JUly
+        return super().apply_annual_aggregation(time_serie[91:-92, ...])
+
 
 class CrocusSwe(Crocus):
 
@@ -45,8 +51,8 @@ class ExtendedCrocusDepth(AbstractExtendedStudy, CrocusDepth):
 
 
 if __name__ == '__main__':
-    for variable_class in [CrocusSweVariable, CrocusDepthVariable]:
-        study = Crocus(variable_class=variable_class, altitude=2400)
+    for variable_clas in [CrocusSweVariable, CrocusDepthVariable]:
+        study = Crocus(variable_class=variable_clas, altitude=2400)
         d = study.year_to_dataset_ordered_dict[1960]
         time_arr = np.array(d.variables['time'])
         print(time_arr)

experiment/meteo_france_SCM_study/crocus/crocus_variables.py

@@ -4,6 +4,8 @@ from experiment.meteo_france_SCM_study.abstract_variable import AbstractVariable
 
 
 class CrocusVariable(AbstractVariable):
+    """Crocus data is every 6 hours. To obtain daily data, we select one data out of 4
+    (in order to have data that will still be comparable to an instantaneous variable"""
 
     def __init__(self, dataset, altitude, variable_name):
         super().__init__(dataset, altitude)
@@ -34,8 +36,6 @@ class CrocusSweVariable(CrocusVariable):
 
 
 class CrocusDepthVariable(CrocusVariable):
-    """Crocus Depth  data is every 6 hours
-    To obtain daily data, we take the average over the 4 slots of 6 hours that compose a full day """
     NAME = 'Snow Depth'
 
     def __init__(self, dataset, altitude):

experiment/meteo_france_SCM_study/safran/safran.py

@@ -73,8 +73,8 @@ if __name__ == '__main__':
     for year, dataset in study.year_to_dataset_ordered_dict.items():
         print('{}: {}'.format(year, dataset.massifsList))
     d = study.year_to_dataset_ordered_dict[1958]
-    # print(d.variables['time'])
+    print(d.variables['time'])
     # print(study.year_to_daily_time_serie[1958].shape)
     # print(len(d.variables['time']))
     # print(study.year_to_annual_total)
-    print(study.df_annual_total.columns)
+    # print(study.df_annual_total.columns)

experiment/meteo_france_SCM_study/safran/safran_variable.py

@@ -72,7 +72,6 @@ class SafranTotalPrecipVariable(AbstractVariable):
         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'):
@@ -82,8 +81,6 @@ class SafranTemperatureVariable(AbstractVariable):
         nb_days = len(self.hourly_temperature) // 24
         self.daily_temperature = [np.mean(self.hourly_temperature[24 * i:24 * (i + 1)], axis=0) for i in range(nb_days)]
 
-
     @property
     def daily_time_serie_array(self):
         return np.array(self.daily_temperature)
-

experiment/meteo_france_SCM_study/visualization/study_visualization/main_study_visualizer.py

@@ -48,23 +48,29 @@ def extended_visualization():
     #         study_visualizer.visualize_all_experimental_law()
 
 
-def annual_mean_vizu_compare_durand_study():
-    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()
+def annual_mean_vizu_compare_durand_study(safran=True, take_mean_value=True):
+    if safran:
+        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(take_mean_value=True)
+    else:
+        for study_class in [CrocusSwe, CrocusDepth][1:]:
+            study = study_class(altitude=1800, year_min=1958, year_max=2005)
+            study_visualizer = StudyVisualizer(study)
+            study_visualizer.visualize_annual_mean_values(take_mean_value=take_mean_value)
 
 
 def normal_visualization():
     save_to_file = False
     only_first_one = True
     # for study_class in SCM_STUDIES[:1]:
-    for study_class in [SafranRainfall, 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])
-            study_visualizer.visualize_annual_mean_values()
-            # study_visualizer.visualize_linear_margin_fit(only_first_max_stable=True)
+            # study_visualizer.visualize_annual_mean_values()
+            study_visualizer.visualize_linear_margin_fit(only_first_max_stable=only_first_one)
 
 
 def complete_analysis(only_first_one=False):
@@ -83,7 +89,7 @@ def complete_analysis(only_first_one=False):
 
 
 if __name__ == '__main__':
-    annual_mean_vizu_compare_durand_study()
+    annual_mean_vizu_compare_durand_study(safran=False, take_mean_value=True)
     # normal_visualization()
     # extended_visualization()
     # complete_analysis()

experiment/meteo_france_SCM_study/visualization/study_visualization/study_visualizer.py

@@ -308,16 +308,19 @@ class StudyVisualizer(object):
             # plot_name = 'Full Likelihood with Linear marginals and max stable dependency structure'
             plt.show()
 
-    def visualize_annual_mean_values(self, ax=None):
+    def visualize_annual_mean_values(self, ax=None, take_mean_value=True):
         if ax is None:
             _, ax = plt.subplots(1, 1, figsize=self.figsize)
 
         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)
+            # We take the mean over all the annual values, otherwise we take the max
+            value = df_annual_total.loc[:, massif_name]
+            value = value.mean() if take_mean_value else value.max()
+            massif_name_to_value[massif_name] = value
+        self.study.visualize_study(ax=ax, massif_name_to_value=massif_name_to_value, show=self.show, add_text=True,
+                                   label=self.study.variable_name)
 
     """ Statistics methods """