[SCM MODELS DATA] refactor abstract study file

experiment/meteo_france_data/scm_models_data/abstract_study.py

@@ -16,7 +16,8 @@ from matplotlib.colors import Normalize
 from netCDF4 import Dataset
 
 from experiment.meteo_france_data.scm_models_data.abstract_variable import AbstractVariable
-from experiment.meteo_france_data.scm_models_data.scm_constants import ALTITUDES, ZS_INT_23, ZS_INT_MASK, LONGITUDES, LATITUDES
+from experiment.meteo_france_data.scm_models_data.scm_constants import ALTITUDES, ZS_INT_23, ZS_INT_MASK, LONGITUDES, \
+    LATITUDES
 from experiment.meteo_france_data.visualization.utils import get_km_formatter
 from extreme_estimator.extreme_models.margin_model.margin_function.abstract_margin_function import \
     AbstractMarginFunction
@@ -48,84 +49,12 @@ class AbstractStudy(object):
         self.year_max = year_max
         self.multiprocessing = multiprocessing
 
-    def write_to_file(self, df: pd.DataFrame):
-        if not op.exists(self.result_full_path):
-            os.makedirs(self.result_full_path, exist_ok=True)
-        df.to_csv(op.join(self.result_full_path, 'merged_array_{}_altitude.csv'.format(self.altitude)))
-
-    """ Data """
-
-    @property
-    def df_all_daily_time_series_concatenated(self) -> pd.DataFrame:
-        df_list = [pd.DataFrame(time_serie, columns=self.study_massif_names) for time_serie in
-                   self.year_to_daily_time_serie_array.values()]
-        df_concatenated = pd.concat(df_list)
-        return df_concatenated
+    """ Annual maxima """
 
     @property
     def observations_annual_maxima(self) -> AnnualMaxima:
         return AnnualMaxima(df_maxima_gev=pd.DataFrame(self.year_to_annual_maxima, index=self.study_massif_names))
 
-    @property
-    def df_annual_total(self) -> pd.DataFrame:
-        return pd.DataFrame(self.year_to_annual_total, index=self.study_massif_names).transpose()
-
-    def annual_aggregation_function(self, *args, **kwargs):
-        raise NotImplementedError()
-
-    """ Load some attributes only once """
-
-    @property
-    def year_to_dataset_ordered_dict(self) -> OrderedDict:
-        print('This code is quite long... '
-              'You should consider year_to_variable which is way faster when multiprocessing=True')
-        # Map each year to the correspond netCDF4 Dataset
-        path_files, ordered_years = self.ordered_years_and_path_files()
-        datasets = [Dataset(path_file) for path_file in path_files]
-        return OrderedDict(zip(ordered_years, datasets))
-
-    def ordered_years_and_path_files(self):
-        nc_files = [(int(f.split('_')[-2][:4]), f) for f in os.listdir(self.study_full_path) if f.endswith('.nc')]
-        ordered_years, path_files = zip(*[(year, op.join(self.study_full_path, nc_file))
-                                          for year, nc_file in sorted(nc_files, key=lambda t: t[0])
-                                          if self.year_min <= year < self.year_max])
-        return path_files, ordered_years
-
-    @property
-    def ordered_years(self):
-        return self.ordered_years_and_path_files()[1]
-
-    @cached_property
-    def year_to_variable_array(self) -> OrderedDict:
-        # Map each year to the variable array
-        path_files, ordered_years = self.ordered_years_and_path_files()
-        if self.multiprocessing:
-            with Pool(NB_CORES) as p:
-                variables = p.map(self.load_variables, path_files)
-        else:
-            variables = [self.load_variables(path_file) for path_file in path_files]
-        return OrderedDict(zip(ordered_years, variables))
-
-    def load_variables(self, path_file):
-        dataset = Dataset(path_file)
-        keyword = self.load_keyword()
-        if isinstance(keyword, str):
-            return np.array(dataset.variables[keyword])
-        else:
-            return [np.array(dataset.variables[k]) for k in keyword]
-
-    def load_keyword(self):
-        return self.variable_class.keyword()
-
-    @property
-    def start_year_and_stop_year(self) -> Tuple[int, int]:
-        ordered_years = self.ordered_years
-        return ordered_years[0], ordered_years[-1]
-
-    @cached_property
-    def year_to_daily_time_serie_array(self) -> OrderedDict:
-        return self._year_to_daily_time_serie_array
-
     @cached_property
     def year_to_annual_maxima(self) -> OrderedDict:
         # Map each year to an array of size nb_massif
@@ -134,6 +63,15 @@ class AbstractStudy(object):
             year_to_annual_maxima[year] = time_serie.max(axis=0)
         return year_to_annual_maxima
 
+    """ Annual total """
+
+    @property
+    def df_annual_total(self) -> pd.DataFrame:
+        return pd.DataFrame(self.year_to_annual_total, index=self.study_massif_names).transpose()
+
+    def annual_aggregation_function(self, *args, **kwargs):
+        raise NotImplementedError()
+
     @cached_property
     def year_to_annual_total(self) -> OrderedDict:
         # Map each year to an array of size nb_massif
@@ -145,10 +83,15 @@ class AbstractStudy(object):
     def apply_annual_aggregation(self, time_serie):
         return self.annual_aggregation_function(time_serie, axis=0)
 
-    def instantiate_variable_object(self, variable_array) -> AbstractVariable:
-        return self.variable_class(variable_array)
+    """ Load daily observations """
 
-    """ Private methods to be overwritten """
+    @cached_property
+    def year_to_daily_time_serie_array(self) -> OrderedDict:
+        return self._year_to_daily_time_serie_array
+
+    @property
+    def _year_to_max_daily_time_serie(self) -> OrderedDict:
+        return self._year_to_daily_time_serie_array
 
     @property
     def _year_to_daily_time_serie_array(self) -> OrderedDict:
@@ -167,78 +110,117 @@ class AbstractStudy(object):
             year_to_daily_time_serie_array[year] = daily_time_serie
         return year_to_daily_time_serie_array
 
-    @property
-    def _year_to_max_daily_time_serie(self) -> OrderedDict:
-        return self._year_to_daily_time_serie_array
+    def instantiate_variable_object(self, variable_array) -> AbstractVariable:
+        return self.variable_class(variable_array)
+
+    """ Load Variables and Datasets """
+
+    @cached_property
+    def year_to_variable_array(self) -> OrderedDict:
+        # Map each year to the variable array
+        path_files, ordered_years = self.ordered_years_and_path_files
+        if self.multiprocessing:
+            with Pool(NB_CORES) as p:
+                variables = p.map(self.load_variables, path_files)
+        else:
+            variables = [self.load_variables(path_file) for path_file in path_files]
+        return OrderedDict(zip(ordered_years, variables))
+
+    def load_variables(self, path_file):
+        dataset = Dataset(path_file)
+        keyword = self.load_keyword()
+        if isinstance(keyword, str):
+            return np.array(dataset.variables[keyword])
+        else:
+            return [np.array(dataset.variables[k]) for k in keyword]
 
-    ##########
+    def load_keyword(self):
+        return self.variable_class.keyword()
 
     @property
-    def study_massif_names(self) -> List[str]:
-        return self.altitude_to_massif_names[self.altitude]
+    def year_to_dataset_ordered_dict(self) -> OrderedDict:
+        print('This code is quite long... '
+              'You should consider year_to_variable which is way faster when multiprocessing=True')
+        # Map each year to the correspond netCDF4 Dataset
+        path_files, ordered_years = self.ordered_years_and_path_files
+        datasets = [Dataset(path_file) for path_file in path_files]
+        return OrderedDict(zip(ordered_years, datasets))
 
     @cached_property
-    def all_massif_names(self) -> List[str]:
-        """
-        Pour l'identification des massifs, le numéro de la variable massif_num correspond à celui de l'attribut num_opp
-        """
-        metadata_path = op.join(self.full_path, self.REANALYSIS_FOLDER, 'metadata')
-        dbf = Dbf5(op.join(metadata_path, 'massifs_alpes.dbf'))
-        df = dbf.to_dataframe().copy()  # type: pd.DataFrame
-        dbf.f.close()
-        df.sort_values(by='num_opp', inplace=True)
-        all_massif_names = list(df['nom'])
-        # Correct a massif name
-        all_massif_names[all_massif_names.index('Beaufortin')] = 'Beaufortain'
-        return all_massif_names
+    def ordered_years_and_path_files(self):
+        nc_files = [(int(f.split('_')[-2][:4]), f) for f in os.listdir(self.study_full_path) if f.endswith('.nc')]
+        ordered_years, path_files = zip(*[(year, op.join(self.study_full_path, nc_file))
+                                          for year, nc_file in sorted(nc_files, key=lambda t: t[0])
+                                          if self.year_min <= year < self.year_max])
+        return path_files, ordered_years
+
+    """ Temporal properties """
 
     @property
-    def original_safran_massif_id_to_massif_name(self) -> Dict[int, str]:
-        return {massif_id: massif_name for massif_id, massif_name in enumerate(self.all_massif_names)}
+    def ordered_years(self):
+        return self.ordered_years_and_path_files[1]
 
-    @cached_property
-    def massifs_coordinates_for_display(self) -> AbstractSpatialCoordinates:
-        # Coordinate object that represents the massif coordinates in Lambert extended
-        # extracted for a csv file, and used only for display purposes
-        df = self.load_df_centroid()
-        # Filter, keep massifs present at the altitude of interest
-        df = df.loc[self.study_massif_names]
-        # Build coordinate object from df_centroid
-        return AbstractSpatialCoordinates.from_df(df)
+    @property
+    def start_year_and_stop_year(self) -> Tuple[int, int]:
+        ordered_years = self.ordered_years
+        return ordered_years[0], ordered_years[-1]
+
+    """ Spatial properties """
+
+    @property
+    def study_massif_names(self) -> List[str]:
+        return self.altitude_to_massif_names[self.altitude]
 
     @property
     def df_massifs_longitude_and_latitude(self) -> pd.DataFrame:
         # DataFrame object that represents the massif coordinates in degrees extracted from the SCM data
+        # Another way of getting the latitudes and longitudes could have been the following:
         # any_ordered_dict = list(self.year_to_dataset_ordered_dict.values())[0]
         # longitude = np.array(any_ordered_dict.variables['longitude'])
         # latitude = np.array(any_ordered_dict.variables['latitude'])
         longitude = np.array(LONGITUDES)
         latitude = np.array(LATITUDES)
-        index = self.altitude_to_massif_names[self.altitude]
         columns = [AbstractSpatialCoordinates.COORDINATE_X, AbstractSpatialCoordinates.COORDINATE_Y]
         data = dict(zip(columns, [longitude[self.altitude_mask], latitude[self.altitude_mask]]))
-        return pd.DataFrame(data=data, index=index, columns=columns)
+        return pd.DataFrame(data=data, index=self.study_massif_names, columns=columns)
 
-    def load_df_centroid(self) -> pd.DataFrame:
-        # Load df_centroid containing all the massif names
-        df_centroid = pd.read_csv(op.join(self.map_full_path, 'coordonnees_massifs_alpes.csv'))
-        df_centroid.set_index('NOM', inplace=True)
-        # Check that the names of massifs are the same
-        symmetric_difference = set(df_centroid.index).symmetric_difference(self.all_massif_names)
-        assert len(symmetric_difference) == 0, symmetric_difference
-        # Sort the column in the order of the SAFRAN dataset
-        df_centroid = df_centroid.reindex(self.all_massif_names, axis=0)
-        for coord_column in [AbstractCoordinates.COORDINATE_X, AbstractCoordinates.COORDINATE_Y]:
-            df_centroid.loc[:, coord_column] = df_centroid[coord_column].str.replace(',', '.').astype(float)
-        return df_centroid
+    @property
+    def missing_massif_name(self):
+        return set(self.all_massif_names) - set(self.altitude_to_massif_names[self.altitude])
+
+    @cached_property
+    def altitude_mask(self):
+        altitude_mask = ZS_INT_MASK == self.altitude
+        assert np.sum(altitude_mask) == len(self.altitude_to_massif_names[self.altitude])
+        return altitude_mask
+
+    """ Path properties """
 
     @property
-    def coordinate_id_to_massif_name(self) -> Dict[int, str]:
-        df_centroid = self.load_df_centroid()
-        return dict(zip(df_centroid['id'], df_centroid.index))
+    def title(self):
+        return "{}/at altitude {}m ({} mountain chains)".format(self.variable_name, self.altitude,
+                                                                len(self.study_massif_names))
+
+    @property
+    def variable_name(self):
+        return self.variable_class.NAME + ' (in {})'.format(self.variable_unit)
+
+    @property
+    def variable_unit(self):
+        return self.variable_class.UNIT
 
     """ Visualization methods """
 
+    @cached_property
+    def massifs_coordinates_for_display(self) -> AbstractSpatialCoordinates:
+        # Coordinate object that represents the massif coordinates in Lambert extended
+        # extracted for a csv file, and used only for display purposes
+        df = self.load_df_centroid()
+        # Filter, keep massifs present at the altitude of interest
+        df = df.loc[self.study_massif_names]
+        # Build coordinate object from df_centroid
+        return AbstractSpatialCoordinates.from_df(df)
+
     def visualize_study(self, ax=None, massif_name_to_value=None, show=True, fill=True, replace_blue_by_white=True,
                         label=None, add_text=False, cmap=None, vmax=100, vmin=0):
         if massif_name_to_value is None:
@@ -300,49 +282,68 @@ class AbstractStudy(object):
         if show:
             plt.show()
 
+    """ 
+    CLASS ATTRIBUTES COMMON TO ALL OBJECTS 
+    (written as object attributes/methods for simplicity)
+    """
+
+    """ Path properties """
+
     @property
-    def idx_to_coords_list(self):
-        df_massif = pd.read_csv(op.join(self.map_full_path, 'massifsalpes.csv'))
-        coord_tuples = [(row_massif['idx'], row_massif[AbstractCoordinates.COORDINATE_X],
-                         row_massif[AbstractCoordinates.COORDINATE_Y])
-                        for _, row_massif in df_massif.iterrows()]
-        all_idxs = set([t[0] for t in coord_tuples])
-        return {idx: [coords for idx_loop, *coords in coord_tuples if idx == idx_loop] for idx in all_idxs}
+    def relative_path(self) -> str:
+        return r'local/spatio_temporal_datasets'
 
     @property
-    def all_coords_list(self):
-        all_values = []
-        for e in self.idx_to_coords_list.values():
-            all_values.extend(e)
-        return list(zip(*all_values))
+    def full_path(self) -> str:
+        return get_full_path(relative_path=self.relative_path)
 
     @property
-    def visualization_x_limits(self):
-        return min(self.all_coords_list[0]), max(self.all_coords_list[0])
+    def map_full_path(self) -> str:
+        return op.join(self.full_path, 'map')
 
     @property
-    def visualization_y_limits(self):
-        return min(self.all_coords_list[1]), max(self.all_coords_list[1])
+    def result_full_path(self) -> str:
+        return op.join(self.full_path, 'results')
+
+    @property
+    def study_full_path(self) -> str:
+        assert self.model_name in ['Safran', 'Crocus']
+        study_folder = 'meteo' if self.model_name is 'Safran' else 'pro'
+        return op.join(self.full_path, self.REANALYSIS_FOLDER, study_folder)
+
+    """  Spatial properties """
+
+    @property
+    def original_safran_massif_id_to_massif_name(self) -> Dict[int, str]:
+        return {massif_id: massif_name for massif_id, massif_name in enumerate(self.all_massif_names)}
 
     @cached_property
-    def mask_french_alps(self):
-        resolution = AbstractMarginFunction.VISUALIZATION_RESOLUTION
-        mask_french_alps = np.zeros([resolution, resolution])
-        for polygon in self.idx_to_coords_list.values():
-            xy_values = list(zip(*polygon))
-            normalized_polygon = []
-            for values, (minlim, max_lim) in zip(xy_values, [self.visualization_x_limits, self.visualization_y_limits]):
-                values -= minlim
-                values *= resolution / (max_lim - minlim)
-                normalized_polygon.append(values)
-            normalized_polygon = list(zip(*normalized_polygon))
-            img = Image.new('L', (resolution, resolution), 0)
-            ImageDraw.Draw(img).polygon(normalized_polygon, outline=1, fill=1)
-            mask_massif = np.array(img)
-            mask_french_alps += mask_massif
-        return ~np.array(mask_french_alps, dtype=bool)
+    def all_massif_names(self) -> List[str]:
+        """
+        Pour l'identification des massifs, le numéro de la variable massif_num correspond à celui de l'attribut num_opp
+        """
+        metadata_path = op.join(self.full_path, self.REANALYSIS_FOLDER, 'metadata')
+        dbf = Dbf5(op.join(metadata_path, 'massifs_alpes.dbf'))
+        df = dbf.to_dataframe().copy()  # type: pd.DataFrame
+        dbf.f.close()
+        df.sort_values(by='num_opp', inplace=True)
+        all_massif_names = list(df['nom'])
+        # Correct a massif name
+        all_massif_names[all_massif_names.index('Beaufortin')] = 'Beaufortain'
+        return all_massif_names
 
-    """ Some properties """
+    def load_df_centroid(self) -> pd.DataFrame:
+        # Load df_centroid containing all the massif names
+        df_centroid = pd.read_csv(op.join(self.map_full_path, 'coordonnees_massifs_alpes.csv'))
+        df_centroid.set_index('NOM', inplace=True)
+        # Check that the names of massifs are the same
+        symmetric_difference = set(df_centroid.index).symmetric_difference(self.all_massif_names)
+        assert len(symmetric_difference) == 0, symmetric_difference
+        # Sort the column in the order of the SAFRAN dataset
+        df_centroid = df_centroid.reindex(self.all_massif_names, axis=0)
+        for coord_column in [AbstractCoordinates.COORDINATE_X, AbstractCoordinates.COORDINATE_Y]:
+            df_centroid.loc[:, coord_column] = df_centroid[coord_column].str.replace(',', '.').astype(float)
+        return df_centroid
 
     @cached_property
     def massif_name_to_altitudes(self) -> Dict[str, List[int]]:
@@ -364,49 +365,51 @@ class AbstractStudy(object):
                 altitude_to_massif_names[altitude].append(massif_name)
         return altitude_to_massif_names
 
-    @property
-    def missing_massif_name(self):
-        return set(self.all_massif_names) - set(self.altitude_to_massif_names[self.altitude])
-
-    @cached_property
-    def altitude_mask(self):
-        altitude_mask = ZS_INT_MASK == self.altitude
-        assert np.sum(altitude_mask) == len(self.altitude_to_massif_names[self.altitude])
-        return altitude_mask
-
-    @property
-    def title(self):
-        return "{}/at altitude {}m ({} mountain chains)".format(self.variable_name, self.altitude,
-                                                                len(self.study_massif_names))
-
-    @property
-    def variable_name(self):
-        return self.variable_class.NAME + ' (in {})'.format(self.variable_unit)
+    """ Visualization methods """
 
     @property
-    def variable_unit(self):
-        return self.variable_class.UNIT
-
-    """ Some path properties """
+    def coordinate_id_to_massif_name(self) -> Dict[int, str]:
+        df_centroid = self.load_df_centroid()
+        return dict(zip(df_centroid['id'], df_centroid.index))
 
     @property
-    def relative_path(self) -> str:
-        return r'local/spatio_temporal_datasets'
+    def idx_to_coords_list(self):
+        df_massif = pd.read_csv(op.join(self.map_full_path, 'massifsalpes.csv'))
+        coord_tuples = [(row_massif['idx'], row_massif[AbstractCoordinates.COORDINATE_X],
+                         row_massif[AbstractCoordinates.COORDINATE_Y])
+                        for _, row_massif in df_massif.iterrows()]
+        all_idxs = set([t[0] for t in coord_tuples])
+        return {idx: [coords for idx_loop, *coords in coord_tuples if idx == idx_loop] for idx in all_idxs}
 
     @property
-    def full_path(self) -> str:
-        return get_full_path(relative_path=self.relative_path)
+    def all_coords_list(self):
+        all_values = []
+        for e in self.idx_to_coords_list.values():
+            all_values.extend(e)
+        return list(zip(*all_values))
 
     @property
-    def map_full_path(self) -> str:
-        return op.join(self.full_path, 'map')
+    def visualization_x_limits(self):
+        return min(self.all_coords_list[0]), max(self.all_coords_list[0])
 
     @property
-    def result_full_path(self) -> str:
-        return op.join(self.full_path, 'results')
+    def visualization_y_limits(self):
+        return min(self.all_coords_list[1]), max(self.all_coords_list[1])
 
-    @property
-    def study_full_path(self) -> str:
-        assert self.model_name in ['Safran', 'Crocus']
-        study_folder = 'meteo' if self.model_name is 'Safran' else 'pro'
-        return op.join(self.full_path, self.REANALYSIS_FOLDER, study_folder)
+    @cached_property
+    def mask_french_alps(self):
+        resolution = AbstractMarginFunction.VISUALIZATION_RESOLUTION
+        mask_french_alps = np.zeros([resolution, resolution])
+        for polygon in self.idx_to_coords_list.values():
+            xy_values = list(zip(*polygon))
+            normalized_polygon = []
+            for values, (minlim, max_lim) in zip(xy_values, [self.visualization_x_limits, self.visualization_y_limits]):
+                values -= minlim
+                values *= resolution / (max_lim - minlim)
+                normalized_polygon.append(values)
+            normalized_polygon = list(zip(*normalized_polygon))
+            img = Image.new('L', (resolution, resolution), 0)
+            ImageDraw.Draw(img).polygon(normalized_polygon, outline=1, fill=1)
+            mask_massif = np.array(img)
+            mask_french_alps += mask_massif
+        return ~np.array(mask_french_alps, dtype=bool)

experiment/meteo_france_data/scm_models_data/safran/safran.py

@@ -67,7 +67,7 @@ if __name__ == '__main__':
     # print(study.all_massif_names)
     # print(study.massif_name_to_altitudes)
     # print(study.year_to_daily_time_serie_array[1958].shape)
-    print(study.missing_massif_name)
+    # print(study.missing_massif_name)
 
     # print(len(d.variables['time']))
     # print(study.year_to_annual_total)

experiment/meteo_france_data/visualization/study_visualization/study_visualizer.py

@@ -749,9 +749,16 @@ class StudyVisualizer(object):
                              for massif_name in self.study.study_massif_names}
         return pd.DataFrame(massif_to_gev_mle, columns=self.study.study_massif_names)
 
+    @property
+    def df_all_daily_time_series_concatenated(self) -> pd.DataFrame:
+        df_list = [pd.DataFrame(time_serie, columns=self.study.study_massif_names) for time_serie in
+                   self.study.year_to_daily_time_serie_array.values()]
+        df_concatenated = pd.concat(df_list)
+        return df_concatenated
+
     def df_gpd_mle(self, threshold) -> pd.DataFrame:
         # Fit a margin fit on each massif
-        massif_to_gev_mle = {massif_name: GpdMleFit(self.study.df_all_daily_time_series_concatenated[massif_name],
+        massif_to_gev_mle = {massif_name: GpdMleFit(self.df_all_daily_time_series_concatenated[massif_name],
                                                     threshold=threshold).gpd_params.summary_serie
                              for massif_name in self.study.study_massif_names}
         return pd.DataFrame(massif_to_gev_mle, columns=self.study.study_massif_names)

test/test_experiment/test_SCM_study.py

@@ -32,7 +32,7 @@ class TestSCMAllStudy(unittest.TestCase):
                                                only_first_one=False, verbose=False,
                                                altitudes=sample(set(ALL_ALTITUDES), k=nb_sample), nb_days=nb_days):
                 self.assertTrue('day' in study.variable_name)
-                first_path_file = study.ordered_years_and_path_files()[0][0]
+                first_path_file = study.ordered_years_and_path_files[0][0]
                 variable_array = study.load_variables(path_file=first_path_file)
                 variable_object = study.instantiate_variable_object(variable_array)
                 self.assertEqual((365, 263), variable_object.daily_time_serie_array.shape,