[SCM] add safran temperature & precipitation for comparison with Durand paper

experiment/meteo_france_SCM_study/abstract_study.py

@@ -21,7 +21,7 @@ class AbstractStudy(object):
     ALTITUDES = [1800, 2400]
 
     def __init__(self, variable_class: type, altitude: int = 1800):
-        assert altitude in self.ALTITUDES
+        assert altitude in self.ALTITUDES, altitude
         self.altitude = altitude
         self.model_name = None
         self.variable_class = variable_class
@@ -45,8 +45,11 @@ class AbstractStudy(object):
         return AnnualMaxima(df_maxima_gev=pd.DataFrame(self.year_to_annual_maxima, index=self.safran_massif_names))
 
     @property
-    def df_annual_mean(self) -> pd.DataFrame:
-        return pd.DataFrame(self.year_to_annual_mean, index=self.safran_massif_names).transpose()
+    def df_annual_total(self) -> pd.DataFrame:
+        return pd.DataFrame(self.year_to_annual_total, index=self.safran_massif_names).transpose()
+
+    def annual_aggregation_function(self, *args, **kwargs):
+        raise NotImplementedError()
 
     """ Load some attributes only once """
 
@@ -72,11 +75,11 @@ class AbstractStudy(object):
         return year_to_annual_maxima
 
     @cached_property
-    def year_to_annual_mean(self) -> OrderedDict:
+    def year_to_annual_total(self) -> OrderedDict:
         # 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.items():
-            year_to_annual_mean[year] = time_serie.mean(axis=0)
+            year_to_annual_mean[year] = self.annual_aggregation_function(time_serie, axis=0)
         return year_to_annual_mean
 
     def instantiate_variable_object(self, dataset) -> AbstractVariable:
@@ -151,7 +154,7 @@ class AbstractStudy(object):
                          row_massif[AbstractCoordinates.COORDINATE_Y])
                         for _, row_massif in df_massif.iterrows()]
 
-        for j, coordinate_id in enumerate(set([t[0] for t in coord_tuples])):
+        for _, coordinate_id in enumerate(set([t[0] for t in coord_tuples])):
             # Retrieve the list of coords (x,y) that define the contour of the massif of id coordinate_id
             coords_list = [coords for idx, *coords in coord_tuples if idx == coordinate_id]
             # if j == 0:
@@ -164,8 +167,20 @@ class AbstractStudy(object):
                 massif_name = self.coordinate_id_to_massif_name[coordinate_id]
                 fill_kwargs = massif_name_to_fill_kwargs[massif_name] if massif_name_to_fill_kwargs is not None else {}
                 ax.fill(*coords_list, **fill_kwargs)
+                # x , y = list(self.massifs_coordinates.df_all_coordinates.loc[massif_name])
+                # x , y= coords_list[0][0], coords_list[0][1]
+                # print(x, y)
+                # print(massif_name)
+                # ax.scatter(x, y)
+                # 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)
+        # 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
+        #     ax.text(x, y, massif_name)
+        # print(self.massifs_coordinates.df_all_coordinates.head())
 
 
         if show:

experiment/meteo_france_SCM_study/crocus/crocus.py

@@ -20,6 +20,9 @@ class Crocus(AbstractStudy):
         suffix = '' if self.altitude == 2400 else ' average of data observed every 6 hours'
         return super().variable_name + suffix
 
+    def annual_aggregation_function(self):
+        return np.mean
+
 
 class CrocusSwe(Crocus):
 

experiment/meteo_france_SCM_study/safran/fit_safran.py

 import pandas as pd
 
-from experiment.meteo_france_SCM_study.safran.safran import ExtendedSafran
+from experiment.meteo_france_SCM_study.safran.safran import ExtendedSafranSnowfall
 from utils import VERSION
 
 
@@ -11,7 +11,7 @@ def fit_mle_gev_for_all_safran_and_different_days():
         for nb_day in [1, 3, 7][:]:
             print('alti: {}, nb_day: {}'.format(safran_alti, nb_day))
             # safran = Safran(safran_alti, nb_day)
-            safran = ExtendedSafran(safran_alti, nb_day)
+            safran = ExtendedSafranSnowfall(safran_alti, nb_day)
             df = safran.df_gev_mle_each_massif
             df.index += ' Safran{} with {} days'.format(safran.altitude, safran.nb_days_of_snowfall)
             dfs.append(df)

experiment/meteo_france_SCM_study/safran/safran.py

+import numpy as np
+
 from experiment.meteo_france_SCM_study.abstract_extended_study import AbstractExtendedStudy
 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_snowfall_variable import SafranSnowfallVariable
+from experiment.meteo_france_SCM_study.safran.safran_variable import SafranSnowfallVariable, \
+    SafranPrecipitationVariable, SafranTemperatureVariable
 
 
 class Safran(AbstractStudy):
 
-    def __init__(self, altitude=1800, nb_days_of_snowfall=1):
-        super().__init__(SafranSnowfallVariable, altitude)
-        self.nb_days_of_snowfall = nb_days_of_snowfall
+    def __init__(self, variable_class: type, altitude: int = 1800):
+        super().__init__(variable_class, altitude)
         self.model_name = 'Safran'
 
+
+class SafranFrequency(Safran):
+
+    def __init__(self, variable_class: type, nb_days_of_snowfall=1, altitude: int = 1800, ):
+        super().__init__(variable_class, altitude)
+        self.nb_days_of_snowfall = nb_days_of_snowfall
+
+
     def instantiate_variable_object(self, dataset) -> AbstractVariable:
         return self.variable_class(dataset, self.nb_days_of_snowfall)
 
@@ -18,16 +28,40 @@ class Safran(AbstractStudy):
     def variable_name(self):
         return super().variable_name + ' cumulated over {} days'.format(self.nb_days_of_snowfall)
 
+    def annual_aggregation_function(self, *args, **kwargs):
+        return np.sum(*args, **kwargs)
+
+
+class SafranSnowfall(SafranFrequency):
 
-class ExtendedSafran(AbstractExtendedStudy, Safran):
+    def __init__(self, altitude: int = 1800, nb_days_of_snowfall=1):
+        super().__init__(SafranSnowfallVariable, nb_days_of_snowfall, altitude)
+
+
+class ExtendedSafranSnowfall(AbstractExtendedStudy, SafranSnowfall):
     pass
 
 
+class SafranPrecipitation(SafranFrequency):
+
+    def __init__(self, altitude: int = 1800, nb_days_of_snowfall=1):
+        super().__init__(SafranPrecipitationVariable, nb_days_of_snowfall, altitude)
+
+
+class SafranTemperature(Safran):
+
+    def __init__(self, altitude: int = 1800):
+        super().__init__(SafranTemperatureVariable, altitude)
+
+    def annual_aggregation_function(self, *args, **kwargs):
+        return np.mean(*args, **kwargs)
+
+
 if __name__ == '__main__':
-    study = Safran()
+    study = SafranSnowfall()
     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_mean)
-    print(study.df_annual_mean)
+    print(study.year_to_annual_total)
+    print(study.df_annual_total)

experiment/meteo_france_SCM_study/safran/safran_snowfall_variable.py → experiment/meteo_france_SCM_study/safran/safran_variable.py

@@ -23,16 +23,16 @@ class SafranSnowfallVariable(AbstractVariable):
 
     NAME = 'Snowfall'
 
-    def __init__(self, dataset, altitude, nb_consecutive_days_of_snowfall=1):
+    def __init__(self, dataset, altitude, nb_consecutive_days_of_snowfall=1, keyword='Snowf'):
         super().__init__(dataset, altitude)
         self.nb_consecutive_days_of_snowfall = nb_consecutive_days_of_snowfall
         # Compute the daily snowfall in kg/m2
-        snowfall_rates = np.array(dataset.variables['Snowf'])
+        snowfall_rates = np.array(dataset.variables[keyword])
         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
         nb_days = len(hourly_snowfall) // 24
-        self.daily_snowfall = [sum(hourly_snowfall[24 * i:24 * (i+1)]) for i in range(nb_days)]
+        self.daily_snowfall = [sum(hourly_snowfall[24 * i:24 * (i + 1)]) for i in range(nb_days)]
 
     @property
     def daily_time_serie(self):
@@ -45,10 +45,23 @@ class SafranSnowfallVariable(AbstractVariable):
         return np.array(snowfall_in_consecutive_days)
 
 
+class SafranPrecipitationVariable(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 SafranTemperatureVariable(AbstractVariable):
 
+    def __init__(self, dataset, altitude, keyword='Tair'):
+        super().__init__(dataset, altitude)
+        # Temperature are in K, I transform them as celsius
+        self.hourly_temperature = np.array(dataset.variables[keyword]) - 273.15
+        nb_days = len(self.hourly_temperature) // 24
+        self.daily_temperature = [np.mean(self.hourly_temperature[24 * i:24 * (i + 1)]) for i in range(nb_days)]
 
 
+    @property
+    def daily_time_serie(self):
+        return self.daily_temperature
 

experiment/meteo_france_SCM_study/visualization/studies_visualization/main_studies_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 Safran, ExtendedSafran
+from experiment.meteo_france_SCM_study.safran.safran import SafranSnowfall, ExtendedSafranSnowfall
 from experiment.meteo_france_SCM_study.visualization.studies_visualization.studies import Studies
 from experiment.meteo_france_SCM_study.visualization.studies_visualization.studies_visualizer import StudiesVisualizer
 
 from experiment.meteo_france_SCM_study.visualization.study_visualization.study_visualizer import StudyVisualizer
 from collections import OrderedDict
 
-SCM_STUDIES = [Safran, CrocusSwe, CrocusDepth]
-SCM_EXTENDED_STUDIES = [ExtendedSafran, ExtendedCrocusSwe, ExtendedCrocusDepth]
+SCM_STUDIES = [SafranSnowfall, CrocusSwe, CrocusDepth]
+SCM_EXTENDED_STUDIES = [ExtendedSafranSnowfall, ExtendedCrocusSwe, ExtendedCrocusDepth]
 SCM_STUDY_TO_EXTENDED_STUDY = OrderedDict(zip(SCM_STUDIES, SCM_EXTENDED_STUDIES))
 
 

experiment/meteo_france_SCM_study/visualization/studies_visualization/studies_visualizer.py

@@ -27,7 +27,7 @@ class StudiesVisualizer(object):
         # Load the dictionary that maps each massif_name to its corresponding time series
         mean_series = []
         for study in self.studies.altitude_to_study.values():
-            mean_serie = study.df_annual_mean.loc[:, massif_names].mean(axis=0)
+            mean_serie = study.df_annual_total.loc[:, massif_names].mean(axis=0)
             mean_series.append(mean_serie)
         df_mean = pd.concat(mean_series, axis=1) # type: pd.DataFrame
         df_mean.columns = self.studies.altitude_list

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 Safran, ExtendedSafran
+from experiment.meteo_france_SCM_study.safran.safran import SafranSnowfall, ExtendedSafranSnowfall, SafranPrecipitation, \
+    SafranTemperature
 
 from experiment.meteo_france_SCM_study.visualization.study_visualization.study_visualizer import StudyVisualizer
 from collections import OrderedDict
 
-SCM_STUDIES = [Safran, CrocusSwe, CrocusDepth]
-SCM_EXTENDED_STUDIES = [ExtendedSafran, ExtendedCrocusSwe, ExtendedCrocusDepth]
+SCM_STUDIES = [SafranSnowfall, CrocusSwe, CrocusDepth]
+SCM_EXTENDED_STUDIES = [ExtendedSafranSnowfall, ExtendedCrocusSwe, ExtendedCrocusDepth]
 SCM_STUDY_TO_EXTENDED_STUDY = OrderedDict(zip(SCM_STUDIES, SCM_EXTENDED_STUDIES))
 
 
 def study_iterator(study_class, only_first_one=False, both_altitude=False, verbose=True):
     all_studies = []
-    is_safran_study = study_class in [Safran, ExtendedSafran]
+    is_safran_study = study_class in [SafranSnowfall, ExtendedSafranSnowfall]
     nb_days = [1] if is_safran_study else [1]
     if verbose:
         print('Loading studies....')
@@ -49,11 +50,12 @@ def extended_visualization():
 def normal_visualization():
     save_to_file = False
     only_first_one = True
-    for study_class in SCM_STUDIES[:1]:
+    # for study_class in SCM_STUDIES[:1]:
+    for study_class in [SafranPrecipitation, 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_mean_daily_values()
+            study_visualizer.visualize_annual_mean_values()
             # study_visualizer.visualize_linear_margin_fit(only_first_max_stable=True)
 
 

experiment/meteo_france_SCM_study/visualization/study_visualization/study_visualizer.py

@@ -308,19 +308,15 @@ class StudyVisualizer(object):
             # plot_name = 'Full Likelihood with Linear marginals and max stable dependency structure'
             plt.show()
 
-    def visualize_mean_daily_values(self, ax=None):
+    def visualize_annual_mean_values(self, ax=None):
         if ax is None:
             _, ax = plt.subplots(1, 1, figsize=self.figsize)
 
-        massif_name_to_value = {
-            massif_name: np.mean(self.get_all_massif_data(massif_id))
-            for massif_id, massif_name in enumerate(self.study.safran_massif_names)
-        }
-
-        self.study.visualize_study(ax=ax, massif_name_to_value=massif_name_to_value, show=False)
-
-        if self.show:
-            plt.show()
+        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):
+            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)
 
     """ Statistics methods """
 

test/test_experiment/test_meteo_france_SCM_study/test_SCM_study.py

@@ -5,7 +5,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 Safran, ExtendedSafran
+from experiment.meteo_france_SCM_study.safran.safran import SafranSnowfall, ExtendedSafranSnowfall
 from experiment.meteo_france_SCM_study.visualization.study_visualization.study_visualizer import StudyVisualizer
 from test.test_utils import load_scm_studies
 
@@ -14,7 +14,7 @@ class TestSCMStudy(unittest.TestCase):
 
     def setUp(self) -> None:
         super().setUp()
-        self.study = Safran()
+        self.study = SafranSnowfall()
 
     def test_massif_safran(self):
         df_centroid = pd.read_csv(op.join(self.study.map_full_path, 'coordonnees_massifs_alpes.csv'))
@@ -22,7 +22,7 @@ class TestSCMStudy(unittest.TestCase):
         assert not set(self.study.safran_massif_names).symmetric_difference(set(df_centroid['NOM']))
 
     def test_extended_run(self):
-        for study_class in [ExtendedSafran, ExtendedCrocusSwe]:
+        for study_class in [ExtendedSafranSnowfall, ExtendedCrocusSwe]:
             for study in study_iterator(study_class, only_first_one=True, both_altitude=False, verbose=False):
                 study_visualizer = StudyVisualizer(study, show=False, save_to_file=False)
                 study_visualizer.visualize_all_mean_and_max_graphs()

test/test_utils.py

@@ -13,7 +13,7 @@ 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 Safran
+from experiment.meteo_france_SCM_study.safran.safran import SafranSnowfall, Safran
 from spatio_temporal_dataset.coordinates.spatial_coordinates.alps_station_3D_coordinates import \
     AlpsStation3DCoordinatesWithAnisotropy
 from spatio_temporal_dataset.coordinates.spatial_coordinates.generated_spatial_coordinates import \
@@ -96,7 +96,7 @@ def load_test_spatiotemporal_coordinates(nb_points, nb_steps, train_split_ratio=
 
 def load_safran_studies(altitudes) -> List[Safran]:
     nb_days_list = [1]
-    return [Safran(safran_altitude, nb_days) for safran_altitude in altitudes for nb_days in nb_days_list]
+    return [SafranSnowfall(safran_altitude, nb_days) for safran_altitude in altitudes for nb_days in nb_days_list]
 
 
 def load_crocus_studies(altitudes) -> List[Crocus]: