[SCM] refactor safran study to treat also crocus study

experiment/safran_study/safran.py → experiment/meteo_france_SCM_study/abstract_study.py

 import os
-from typing import List
-
 import os.path as op
 from collections import OrderedDict
+from typing import List
 
 import matplotlib.pyplot as plt
 import pandas as pd
 from netCDF4 import Dataset
 
-from experiment.safran_study.massif import safran_massif_names_from_datasets
-from experiment.safran_study.snowfall_annual_maxima import SafranSnowfall
+from experiment.meteo_france_SCM_study.abstract_variable import AbstractVariable
+from experiment.meteo_france_SCM_study.massif import safran_massif_names_from_datasets
 from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
 from spatio_temporal_dataset.coordinates.spatial_coordinates.abstract_spatial_coordinates import \
     AbstractSpatialCoordinates
@@ -17,12 +16,13 @@ from spatio_temporal_dataset.spatio_temporal_observations.annual_maxima_observat
 from utils import get_full_path, cached_property
 
 
-class Safran(object):
+class AbstractStudy(object):
 
-    def __init__(self, safran_altitude=1800, nb_days_of_snowfall=1):
+    def __init__(self, safran_altitude=1800):
         assert safran_altitude in [1800, 2400]
         self.safran_altitude = safran_altitude
-        self.nb_days_of_snowfall = nb_days_of_snowfall
+        self.model_name = None
+        self.variable_class = None
 
     def write_to_file(self, df):
         if not op.exists(self.result_full_path):
@@ -33,17 +33,14 @@ class Safran(object):
 
     @property
     def df_all_snowfall_concatenated(self):
-        df_list = [pd.DataFrame(snowfall, columns=self.safran_massif_names) for snowfall in self.year_to_snowfall.values()]
+        df_list = [pd.DataFrame(snowfall, columns=self.safran_massif_names) for snowfall in
+                   self.year_to_daily_time_serie.values()]
         df_concatenated = pd.concat(df_list)
         return df_concatenated
 
-    @property
-    def df_annual_maxima(self):
-        return pd.DataFrame(self.year_to_annual_maxima, index=self.safran_massif_names).T
-
     @property
     def observations_annual_maxima(self):
-        return AnnualMaxima(df_maxima_gev=self.df_annual_maxima.T)
+        return AnnualMaxima(df_maxima_gev=pd.DataFrame(self.year_to_annual_maxima, index=self.safran_massif_names))
 
     """ Load some attributes only once """
 
@@ -51,19 +48,33 @@ class Safran(object):
     def year_to_annual_maxima(self) -> OrderedDict:
         # Map each year to an array of size nb_massif
         year_to_annual_maxima = OrderedDict()
-        for year, snowfall in self.year_to_snowfall.items():
-            year_to_annual_maxima[year] = snowfall.max(axis=0)
+        for year, time_serie in self.year_to_daily_time_serie.items():
+            year_to_annual_maxima[year] = time_serie.max(axis=0)
         return year_to_annual_maxima
 
     @cached_property
-    def year_to_snowfall(self) -> OrderedDict:
+    def year_to_dataset_ordered_dict(self) -> OrderedDict:
+        # Map each year to the correspond netCDF4 Dataset
+        year_to_dataset = OrderedDict()
+        nc_files = [(int(f.split('_')[1][:4]), f) for f in os.listdir(self.safran_full_path) if f.endswith('.nc')]
+        for year, nc_file in sorted(nc_files, key=lambda t: t[0]):
+            year_to_dataset[year] = Dataset(op.join(self.safran_full_path, nc_file))
+        return year_to_dataset
+
+    @cached_property
+    def year_to_daily_time_serie(self) -> OrderedDict:
         # Map each year to a matrix of size 365-nb_days_consecutive+1 x nb_massifs
-        year_to_safran_snowfall = {year: SafranSnowfall(dataset) for year, dataset in
-                                   self.year_to_dataset_ordered_dict.items()}
-        year_to_snowfall = OrderedDict()
+        year_to_variable = {year: self.instantiate_variable_object(dataset) for year, dataset in
+                            self.year_to_dataset_ordered_dict.items()}
+        year_to_daily_time_serie = OrderedDict()
         for year in self.year_to_dataset_ordered_dict.keys():
-            year_to_snowfall[year] = year_to_safran_snowfall[year].annual_snowfall(self.nb_days_of_snowfall)
-        return year_to_snowfall
+            year_to_daily_time_serie[year] = year_to_variable[year].daily_time_serie
+        return year_to_daily_time_serie
+
+    def instantiate_variable_object(self, dataset) -> AbstractVariable:
+        return self.variable_class(dataset)
+
+    ##########
 
     @property
     def safran_massif_names(self) -> List[str]:
@@ -74,15 +85,6 @@ class Safran(object):
     def safran_massif_id_to_massif_name(self):
         return dict(enumerate(self.safran_massif_names))
 
-    @cached_property
-    def year_to_dataset_ordered_dict(self) -> OrderedDict:
-        # Map each year to the correspond netCDF4 Dataset
-        year_to_dataset = OrderedDict()
-        nc_files = [(int(f.split('_')[1][:4]), f) for f in os.listdir(self.safran_full_path) if f.endswith('.nc')]
-        for year, nc_file in sorted(nc_files, key=lambda t: t[0]):
-            year_to_dataset[year] = Dataset(op.join(self.safran_full_path, nc_file))
-        return year_to_dataset
-
     @cached_property
     def massifs_coordinates(self) -> AbstractSpatialCoordinates:
         # Coordinate object that represents the massif coordinates in Lambert extended
@@ -105,7 +107,6 @@ class Safran(object):
         df_centroid = self.load_df_centroid()
         return dict(zip(df_centroid['id'], df_centroid.index))
 
-
     """ Visualization methods """
 
     def visualize(self, ax=None, massif_name_to_fill_kwargs=None, show=True, fill=True):
@@ -142,7 +143,8 @@ class Safran(object):
 
     @property
     def safran_full_path(self) -> str:
-        return op.join(self.full_path, 'safran-crocus_{}'.format(self.safran_altitude), 'Safran')
+        assert self.model_name in ['Safran', 'Crocus']
+        return op.join(self.full_path, 'safran-crocus_{}'.format(self.safran_altitude), self.model_name)
 
     @property
     def map_full_path(self) -> str:
@@ -150,4 +152,4 @@ class Safran(object):
 
     @property
     def result_full_path(self) -> str:
-        return op.join(self.full_path, 'results')
\ No newline at end of file
+        return op.join(self.full_path, 'results')

experiment/meteo_france_SCM_study/abstract_variable.py

+
+
+class AbstractVariable(object):
+
+    def __init__(self, dataset):
+        self.dataset = dataset
+
+    @property
+    def daily_time_serie(self):
+        # Return an array of size length of time series x nb_massif
+        raise NotImplementedError
\ No newline at end of file

experiment/meteo_france_SCM_study/crocus/crocus.py

+from experiment.meteo_france_SCM_study.abstract_study import AbstractStudy
+from experiment.meteo_france_SCM_study.crocus.crocus_variables import CrocusSweVariable, CrocusDepthVariable
+
+
+class Crocus(AbstractStudy):
+
+    def __init__(self, safran_altitude=1800, crocus_variable_class=CrocusSweVariable):
+        super().__init__(safran_altitude)
+        self.model_name = 'Crocus'
+        assert crocus_variable_class in [CrocusSweVariable, CrocusDepthVariable]
+        self.variable_class = crocus_variable_class
+
+if __name__ == '__main__':
+    for variable_class in [CrocusSweVariable, CrocusDepthVariable]:
+        study = Crocus(crocus_variable_class=variable_class)
+        # d = study.year_to_dataset_ordered_dict[1960]
+        # print(d)
+        a = study.year_to_daily_time_serie[1960]
+        print(a.shape)
\ No newline at end of file

experiment/meteo_france_SCM_study/crocus/crocus_variables.py

+import numpy as np
+
+from experiment.meteo_france_SCM_study.abstract_variable import AbstractVariable
+
+
+class CrocusVariable(AbstractVariable):
+
+    def __init__(self, dataset, variable_name):
+        super().__init__(dataset)
+        self.variable_name = variable_name
+
+    @property
+    def daily_time_serie(self):
+        # So far the dimension of the time serie is 1460 x 23
+        return np.array(self.dataset.variables[self.variable_name])[:, 0, :]
+
+
+class CrocusSweVariable(CrocusVariable):
+
+    def __init__(self, dataset):
+        super().__init__(dataset, 'SNOWSWE')
+
+
+class CrocusDepthVariable(CrocusVariable):
+
+    def __init__(self, dataset):
+        super().__init__(dataset, "SNOWDEPTH")

experiment/safran_study/main_visualize_safran.py → experiment/meteo_france_SCM_study/main_visualize.py

-from experiment.safran_study.safran import Safran
+from experiment.meteo_france_SCM_study.safran.safran import Safran
 from itertools import product
 
-from experiment.safran_study.safran_visualizer import SafranVisualizer
+from experiment.meteo_france_SCM_study.safran.safran_visualizer import SafranVisualizer
 
 
 def load_all_safran(only_first_one=False):

experiment/safran_study/fit_safran.py → experiment/meteo_france_SCM_study/safran/fit_safran.py

 import pandas as pd
 
-from experiment.safran_study.safran_extended import ExtendedSafran
+from experiment.meteo_france_SCM_study.safran.safran_extended import ExtendedSafran
 from utils import VERSION
 
 

experiment/meteo_france_SCM_study/safran/safran.py

+from typing import List
+
+import os.path as op
+from collections import OrderedDict
+
+import matplotlib.pyplot as plt
+import pandas as pd
+
+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.massif import safran_massif_names_from_datasets
+from experiment.meteo_france_SCM_study.safran.safran_snowfall_variable import SafranSnowfallVariable
+from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
+from spatio_temporal_dataset.coordinates.spatial_coordinates.abstract_spatial_coordinates import \
+    AbstractSpatialCoordinates
+from spatio_temporal_dataset.spatio_temporal_observations.annual_maxima_observations import AnnualMaxima
+from utils import cached_property
+
+
+class Safran(AbstractStudy):
+
+    def __init__(self, safran_altitude=1800, nb_days_of_snowfall=1):
+        super().__init__(safran_altitude)
+        self.nb_days_of_snowfall = nb_days_of_snowfall
+        self.model_name = 'Safran'
+        self.variable_class = SafranSnowfallVariable
+
+
+    def instantiate_variable_object(self, dataset) -> AbstractVariable:
+        return self.variable_class(dataset, self.nb_days_of_snowfall)
+
+

experiment/safran_study/safran_extended.py → experiment/meteo_france_SCM_study/safran/safran_extended.py

@@ -2,7 +2,7 @@ from collections import OrderedDict
 
 import pandas as pd
 
-from experiment.safran_study.safran import Safran
+from experiment.meteo_france_SCM_study.safran.safran import Safran
 from utils import cached_property
 
 

experiment/safran_study/snowfall_annual_maxima.py → experiment/meteo_france_SCM_study/safran/safran_snowfall_variable.py

 import numpy as np
 
+from experiment.meteo_france_SCM_study.abstract_variable import AbstractVariable
 
-class SafranSnowfall(object):
+
+class SafranSnowfallVariable(AbstractVariable):
     """"
     Hypothesis:
 
@@ -17,12 +19,9 @@ class SafranSnowfall(object):
         (but here the problem might be that the x_i are not idnependent, they are highly dependent one from another)
     """
 
-    def __init__(self, dataset):
-        self.dataset = dataset
-        # Corresponding starting times
-        # start_datetime = datetime(year=year, month=8, day=1, hour=6)
-        # start_times = np.array(dataset.variables['time'])[:-1]
-
+    def __init__(self, dataset, nb_consecutive_days_of_snowfall=1):
+        super().__init__(dataset)
+        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'])
         mean_snowfall_rates = 0.5 * (snowfall_rates[:-1] + snowfall_rates[1:])
@@ -31,13 +30,13 @@ class SafranSnowfall(object):
         nb_days = len(hourly_snowfall) // 24
         self.daily_snowfall = [sum(hourly_snowfall[24 * i:24 * (i+1)]) for i in range(nb_days)]
 
-    def annual_snowfall(self, nb_days_of_snowfall=1):
+    @property
+    def daily_time_serie(self):
         # Aggregate the daily snowfall by the number of consecutive days
-        shifted_list = [self.daily_snowfall[i:] for i in range(nb_days_of_snowfall)]
+        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)]
-        # Return the maximum of the aggregated list
         # The returned array is of size n-nb_days+1 x nb_massif
         return np.array(snowfall_in_consecutive_days)
 

experiment/safran_study/safran_visualizer.py → experiment/meteo_france_SCM_study/safran/safran_visualizer.py

-import matplotlib as mpl
-import matplotlib.colorbar as cbar
 import matplotlib.pyplot as plt
 import pandas as pd
-from matplotlib import cm
-from mpl_toolkits.axes_grid1 import make_axes_locatable
 
 from extreme_estimator.estimator.full_estimator.abstract_full_estimator import \
     FullEstimatorInASingleStepWithSmoothMargin
 from extreme_estimator.estimator.margin_estimator.abstract_margin_estimator import SmoothMarginEstimator
 from extreme_estimator.extreme_models.margin_model.smooth_margin_model import LinearAllParametersAllDimsMarginModel
 from extreme_estimator.extreme_models.max_stable_model.max_stable_models import Smith
-from extreme_estimator.margin_fits.extreme_params import ExtremeParams
 from extreme_estimator.margin_fits.gev.gev_params import GevParams
 from extreme_estimator.margin_fits.gev.gevmle_fit import GevMleFit
 from extreme_estimator.margin_fits.gpd.gpd_params import GpdParams
 from extreme_estimator.margin_fits.gpd.gpdmle_fit import GpdMleFit
-from experiment.safran_study.safran import Safran
-from extreme_estimator.margin_fits.plot.create_shifted_cmap import plot_extreme_param, get_color_rbga
+from experiment.meteo_france_SCM_study.safran.safran import Safran
+from extreme_estimator.margin_fits.plot.create_shifted_cmap import get_color_rbga_shifted
 from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
 
 
@@ -76,7 +71,7 @@ class SafranVisualizer(object):
             massif_name_to_value = df.loc[gev_param_name, :].to_dict()
             # Compute the middle point of the values for the color map
             values = list(massif_name_to_value.values())
-            colors = get_color_rbga(ax, gev_param_name, values)
+            colors = get_color_rbga_shifted(ax, gev_param_name, values)
             massif_name_to_fill_kwargs = {massif_name: {'color': color} for massif_name, color in
                                           zip(self.safran.safran_massif_names, colors)}
             self.safran.visualize(ax=ax, massif_name_to_fill_kwargs=massif_name_to_fill_kwargs, show=False)
@@ -98,7 +93,7 @@ class SafranVisualizer(object):
     @property
     def df_gev_mle_each_massif(self):
         # Fit a margin_fits on each massif
-        massif_to_gev_mle = {massif_name: GevMleFit(self.safran.df_annual_maxima[massif_name]).gev_params.summary_serie
+        massif_to_gev_mle = {massif_name: GevMleFit(self.safran.observations_annual_maxima.loc[massif_name]).gev_params.summary_serie
                              for massif_name in self.safran.safran_massif_names}
         return pd.DataFrame(massif_to_gev_mle, columns=self.safran.safran_massif_names)
 

extreme_estimator/extreme_models/margin_model/margin_function/abstract_margin_function.py

@@ -6,7 +6,7 @@ import numpy as np
 import pandas as pd
 
 from extreme_estimator.margin_fits.gev.gev_params import GevParams
-from extreme_estimator.margin_fits.plot.create_shifted_cmap import plot_extreme_param
+from extreme_estimator.margin_fits.plot.create_shifted_cmap import plot_extreme_param, imshow_shifted
 from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
 from spatio_temporal_dataset.slicer.split import Split
 
@@ -47,7 +47,8 @@ class AbstractMarginFunction(object):
 
     # Visualization function
 
-    def set_datapoint_display_parameters(self, spatio_temporal_split=Split.all, datapoint_marker=None, filter=None, color=None,
+    def set_datapoint_display_parameters(self, spatio_temporal_split=Split.all, datapoint_marker=None, filter=None,
+                                         color=None,
                                          linewidth=1, datapoint_display=False):
         self.datapoint_display = datapoint_display
         self.spatio_temporal_split = spatio_temporal_split
@@ -132,13 +133,10 @@ class AbstractMarginFunction(object):
         grid = self.grid_2D(x, y)
         if ax is None:
             ax = plt.gca()
-        imshow_method = ax.imshow
-        values = grid[gev_param_name]
 
-        norm, shifted_cmap = plot_extreme_param(ax, gev_param_name, values)
+        # Special display
+        imshow_shifted(ax, gev_param_name, grid[gev_param_name], x, y)
 
-        imshow_method(values, extent=(x.min(), x.max(), y.min(), y.max()),
-                      interpolation='nearest', cmap=shifted_cmap)
         # X axis
         ax.set_xlabel('coordinate X')
         plt.setp(ax.get_xticklabels(), visible=True)
@@ -151,7 +149,6 @@ class AbstractMarginFunction(object):
         if show:
             plt.show()
 
-
     def grid_2D(self, x, y):
         resolution = 100
         grid = []

extreme_estimator/margin_fits/plot/create_shifted_cmap.py

@@ -32,7 +32,7 @@ def plot_extreme_param(ax, gev_param_name, values):
     return norm, shifted_cmap
 
 
-def get_color_rbga(ax, gev_param_name, values):
+def get_color_rbga_shifted(ax, gev_param_name, values):
     """
     For some display it was necessary to transform dark blue values into white values
     """
@@ -40,5 +40,20 @@ def get_color_rbga(ax, gev_param_name, values):
     m = cm.ScalarMappable(norm=norm, cmap=shifted_cmap)
     colors = [m.to_rgba(value) for value in values]
     if gev_param_name != ExtremeParams.SHAPE:
-        colors = [color if color != (0, 0, 1, 1) else (1, 1, 1, 1) for color in colors]
+        colors = [color if color[2] == 1 else (1, 1, 1, 1) for color in colors]
     return colors
+
+
+def imshow_shifted(ax, gev_param_name, values, x, y):
+    norm, shifted_cmap = plot_extreme_param(ax, gev_param_name, values)
+    shifted_cmap.set_bad(color='white')
+
+    masked_array = values
+    if gev_param_name != ExtremeParams.SHAPE:
+        epsilon = 1.0
+        value = np.min(values)
+        # The right blue corner will be blue (but most of the time, another display will be on top)
+        masked_array[-1, -1] = value - epsilon
+
+    ax.imshow(masked_array, extent=(x.min(), x.max(), y.min(), y.max()), cmap=shifted_cmap)
+

test/test_experiment/test_safran/test_safran_visualizer.py

 # import unittest
 #
-# from experiment.safran_study.safran import Safran
-# from experiment.safran_study.safran_visualizer import SafranVisualizer
+# from experiment.meteo_france_SCM_study.safran import Safran
+# from experiment.meteo_france_SCM_study.safran_visualizer import SafranVisualizer
 #
 #
 # class TestSafranVisualizer(unittest.TestCase):

test/test_utils.py

@@ -7,7 +7,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.safran_study.safran import Safran
+from experiment.meteo_france_SCM_study.safran.safran import 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 \