diff --git a/experiment/meteo_france_SCM_study/abstract_study.py b/experiment/meteo_france_SCM_study/abstract_study.py
index dcb669986930a37e7d31a96c2719f0db28b868de..7b60a4c4af2526b152ae240dc466322f25cb97d4 100644
--- a/experiment/meteo_france_SCM_study/abstract_study.py
+++ b/experiment/meteo_france_SCM_study/abstract_study.py
@@ -49,7 +49,7 @@ class AbstractStudy(object):
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')]
+ nc_files = [(int(f.split('_')[-2][: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
@@ -94,7 +94,7 @@ class AbstractStudy(object):
@property
def original_safran_massif_names(self):
# Load the names of the massif as defined by SAFRAN
- return safran_massif_names_from_datasets(list(self.year_to_dataset_ordered_dict.values()))
+ return safran_massif_names_from_datasets(list(self.year_to_dataset_ordered_dict.values()), self.altitude)
@property
def original_safran_massif_id_to_massif_name(self) -> Dict[int, str]:
diff --git a/experiment/meteo_france_SCM_study/main_visualize.py b/experiment/meteo_france_SCM_study/main_visualize.py
index 07d8a01033a6f035cc53bc09cda6cd1d33be02f3..9110ac96520621533c500d425dd78a9d399fb30e 100644
--- a/experiment/meteo_france_SCM_study/main_visualize.py
+++ b/experiment/meteo_france_SCM_study/main_visualize.py
@@ -5,37 +5,63 @@ from experiment.meteo_france_SCM_study.safran.safran import Safran, ExtendedSafr
from itertools import product
from experiment.meteo_france_SCM_study.safran.safran_visualizer import StudyVisualizer
+from collections import OrderedDict
+SCM_STUDIES = [Safran, CrocusSwe, CrocusDepth]
+SCM_EXTENDED_STUDIES = [ExtendedSafran, ExtendedCrocusSwe, ExtendedCrocusDepth]
+SCM_STUDY_TO_EXTENDED_STUDY = OrderedDict(zip(SCM_STUDIES, SCM_EXTENDED_STUDIES))
-def load_all_studies(study_class, only_first_one=False):
+
+def study_iterator(study_class, only_first_one=False, both_altitude=False, verbose=True):
all_studies = []
- is_safran_study = study_class == Safran
+ is_safran_study = study_class in [Safran, ExtendedSafran]
nb_days = [3, 1] if is_safran_study else [1]
- for alti, nb_day in product(AbstractStudy.ALTITUDES, nb_days):
- print('alti: {}, nb_day: {}'.format(alti, nb_day))
- study = Safran(alti, nb_day) if is_safran_study else study_class(alti)
- all_studies.append(study)
+ if verbose:
+ print('Loading studies....')
+ for nb_day in nb_days:
+ for alti in AbstractStudy.ALTITUDES[::1]:
+ if verbose:
+ print('alti: {}, nb_day: {}'.format(alti, nb_day))
+ study = study_class(alti, nb_day) if is_safran_study else study_class(alti)
+ yield study
+ if only_first_one and not both_altitude:
+ break
if only_first_one:
break
+
return all_studies
def extended_visualization():
- for study_class in [ExtendedSafran, ExtendedCrocusSwe, ExtendedCrocusDepth][:1]:
- for study in load_all_studies(study_class, only_first_one=True):
+ for study_class in SCM_EXTENDED_STUDIES[1:2]:
+ for study in study_iterator(study_class, only_first_one=True):
study_visualizer = StudyVisualizer(study)
study_visualizer.visualize_all_kde_graphs()
def normal_visualization():
- for study_class in [Safran, CrocusSwe, CrocusDepth][:1]:
- for study in load_all_studies(study_class, only_first_one=True):
+ for study_class in SCM_STUDIES[:1]:
+ for study in study_iterator(study_class, only_first_one=True):
study_visualizer = StudyVisualizer(study)
# study_visualizer.visualize_independent_margin_fits(threshold=[None, 20, 40, 60][0])
- # study_visualizer.visualize_smooth_margin_fit()
- study_visualizer.visualize_full_fit()
+ study_visualizer.visualize_linear_margin_fit()
+
+
+def complete_analysis(only_first_one=False):
+ """An overview of everything that is possible with study OR extended study"""
+ for study_class, extended_study_class in list(SCM_STUDY_TO_EXTENDED_STUDY.items())[:]:
+ # First explore everything you can do with the extended study class
+ print('Extended study')
+ for extended_study in study_iterator(extended_study_class, only_first_one=only_first_one):
+ study_visualizer = StudyVisualizer(extended_study, save_to_file=True)
+ study_visualizer.visualize_all_kde_graphs()
+ print('Study normal')
+ for study in study_iterator(study_class, only_first_one=only_first_one):
+ study_visualizer = StudyVisualizer(study, save_to_file=True)
+ study_visualizer.visualize_linear_margin_fit()
if __name__ == '__main__':
- normal_visualization()
- # extended_visualization()
\ No newline at end of file
+ # normal_visualization()
+ extended_visualization()
+ # complete_analysis()
diff --git a/experiment/meteo_france_SCM_study/massif.py b/experiment/meteo_france_SCM_study/massif.py
index eb65d7baf6db95804e60f031a8cc5f6f8ef3901c..1fad20c9442bcf9679c118bdf65d716e4f0c2bda 100644
--- a/experiment/meteo_france_SCM_study/massif.py
+++ b/experiment/meteo_france_SCM_study/massif.py
@@ -1,8 +1,14 @@
from utils import first
-MASSIF_NAMES = ['Pelvoux', 'Queyras', 'Mont-Blanc', 'Aravis', 'Haute-Tarentaise', 'Vercors', 'Alpes-Azur', 'Oisans',
- 'Mercantour', 'Chartreuse', 'Haute-Maurienne', 'Belledonne', 'Thabor', 'Parpaillon', 'Bauges',
- 'Chablais', 'Ubaye', 'Grandes-Rousses', 'Devoluy', 'Champsaur', 'Vanoise', 'Beaufortain', 'Maurienne']
+MASSIF_NAMES_1800 = ['Pelvoux', 'Queyras', 'Mont-Blanc', 'Aravis', 'Haute-Tarentaise', 'Vercors', 'Alpes-Azur',
+ 'Oisans',
+ 'Mercantour', 'Chartreuse', 'Haute-Maurienne', 'Belledonne', 'Thabor', 'Parpaillon', 'Bauges',
+ 'Chablais', 'Ubaye', 'Grandes-Rousses', 'Devoluy', 'Champsaur', 'Vanoise', 'Beaufortain',
+ 'Maurienne']
+# Some massif like Chartreuse do not have massif whose altitude is higher or equal to 2400
+MASSIF_NAMES_2400 = ['Pelvoux', 'Queyras', 'Mont-Blanc', 'Aravis', 'Haute-Tarentaise', 'Vercors', 'Alpes-Azur', 'Oisans'
+ , 'Mercantour', 'Haute-Maurienne', 'Belledonne', 'Thabor', 'Parpaillon', 'Chablais', 'Ubaye', 'Grandes-Rousses',
+ 'Devoluy', 'Champsaur', 'Vanoise', 'Beaufortain', 'Maurienne']
class Massif(object):
@@ -19,12 +25,13 @@ class Massif(object):
return cls(name.strip(), int(id), float(lat), float(lon))
-def safran_massif_names_from_datasets(datasets):
+def safran_massif_names_from_datasets(datasets, altitude):
# Massifs names are extracted from SAFRAN dataset
+ reference_massif_list = MASSIF_NAMES_1800 if altitude == 1800 else MASSIF_NAMES_2400
if hasattr(datasets[0], 'massifsList'):
# Assert the all the datasets have the same indexing for the massif
assert len(set([dataset.massifsList for dataset in datasets])) == 1
# List of the name of the massif used by all the SAFRAN datasets
safran_names = [Massif.from_str(massif_str).name for massif_str in first(datasets).massifsList.split('/')]
- assert MASSIF_NAMES == safran_names
- return MASSIF_NAMES
+ assert reference_massif_list == safran_names, '{} \n{}'.format(reference_massif_list, safran_names)
+ return reference_massif_list
diff --git a/experiment/meteo_france_SCM_study/safran/safran_visualizer.py b/experiment/meteo_france_SCM_study/safran/safran_visualizer.py
index d211fb387490ed0cf85696b6fdd9d43088e52e65..113558758e4a6f27fcfeb0aad0e5b971d03176aa 100644
--- a/experiment/meteo_france_SCM_study/safran/safran_visualizer.py
+++ b/experiment/meteo_france_SCM_study/safran/safran_visualizer.py
@@ -1,4 +1,6 @@
import math
+import os
+import os.path as op
import numpy as np
import matplotlib.pyplot as plt
@@ -17,15 +19,17 @@ from extreme_estimator.margin_fits.gpd.gpdmle_fit import GpdMleFit
from extreme_estimator.margin_fits.plot.create_shifted_cmap import get_color_rbga_shifted
from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
from test.test_utils import load_test_max_stable_models
-from utils import get_display_name_from_object_type
+from utils import get_display_name_from_object_type, VERSION, VERSION_TIME
class StudyVisualizer(object):
- def __init__(self, study: AbstractStudy, show=True):
+ def __init__(self, study: AbstractStudy, show=True, save_to_file=False):
+ self.save_to_file = save_to_file
self.study = study
- self.show = show
+ self.show = False if self.save_to_file else show
self.window_size_for_smoothing = 21
+ self.figsize=(16.0, 10.0)
@property
def observations(self):
@@ -39,21 +43,18 @@ class StudyVisualizer(object):
def dataset(self):
return AbstractDataset(self.observations, self.coordinates)
- def visualize_all_kde_graphs(self, show=True):
+ def visualize_all_kde_graphs(self):
massif_names = self.study.safran_massif_names
nb_columns = 5
nb_rows = math.ceil(len(massif_names) / nb_columns)
- fig, axes = plt.subplots(nb_rows, nb_columns)
+ fig, axes = plt.subplots(nb_rows, nb_columns, figsize=self.figsize)
fig.subplots_adjust(hspace=1.0, wspace=1.0)
for i, massif_name in enumerate(massif_names):
row_id, column_id = i // nb_columns, i % nb_columns
ax = axes[row_id, column_id]
self.visualize_kde_graph(ax, i, massif_name)
- title = self.study.title
- title += ' (mean computed with a sliding window of size {})'.format(self.window_size_for_smoothing)
- fig.suptitle(title)
- if show:
- plt.show()
+ plot_name = ' mean with sliding window of size {}'.format(self.window_size_for_smoothing)
+ self.show_or_save_to_file(plot_name)
def visualize_kde_graph(self, ax, i, massif_name):
self.maxima_plot(ax, i)
@@ -91,34 +92,43 @@ class StudyVisualizer(object):
assert len(x) == len(y)
return x, y
- def fit_and_visualize_estimator(self, estimator, axes=None, show=True, title=None):
+ def visualize_linear_margin_fit(self):
+ plot_name = 'Full Likelihood with Linear marginals and max stable dependency structure'
+ max_stable_models = load_test_max_stable_models(only_one_covariance_function=True)[:1]
+ fig, axes = plt.subplots(len(max_stable_models) + 1, len(GevParams.SUMMARY_NAMES), figsize=self.figsize)
+ fig.subplots_adjust(hspace=1.0, wspace=1.0)
+ margin_class = LinearAllParametersAllDimsMarginModel
+ # Plot the smooth margin only
+ margin_model = margin_class(coordinates=self.coordinates)
+ estimator = SmoothMarginEstimator(dataset=self.dataset, margin_model=margin_model)
+ self.fit_and_visualize_estimator(estimator, axes[0], title='without max stable')
+ # Plot the smooth margin fitted with a max stable
+ for i, max_stable_model in enumerate(max_stable_models, 1):
+ margin_model = margin_class(coordinates=self.coordinates)
+ estimator = FullEstimatorInASingleStepWithSmoothMargin(self.dataset, margin_model, max_stable_model)
+ title = get_display_name_from_object_type(type(max_stable_model))
+ self.fit_and_visualize_estimator(estimator, axes[i], title=title)
+ self.show_or_save_to_file(plot_name)
+
+ def fit_and_visualize_estimator(self, estimator, axes=None, title=None):
estimator.fit()
axes = estimator.margin_function_fitted.visualize_function(show=False, axes=axes, title=title)
for ax in axes:
self.study.visualize(ax, fill=False, show=False)
- if show:
- plt.suptitle(self.study.title)
- plt.show()
- def visualize_smooth_margin_fit(self):
- margin_model = LinearAllParametersAllDimsMarginModel(coordinates=self.coordinates)
- estimator = SmoothMarginEstimator(dataset=self.dataset, margin_model=margin_model)
- self.fit_and_visualize_estimator(estimator)
-
- def visualize_full_fit(self):
- max_stable_models = load_test_max_stable_models()
- fig, axes = plt.subplots(len(max_stable_models), len(GevParams.SUMMARY_NAMES))
- fig.subplots_adjust(hspace=1.0, wspace=1.0)
- for i, max_stable_model in enumerate(max_stable_models):
- margin_model = LinearAllParametersAllDimsMarginModel(coordinates=self.coordinates)
- estimator = FullEstimatorInASingleStepWithSmoothMargin(self.dataset, margin_model, max_stable_model)
- title = get_display_name_from_object_type(type(max_stable_model))
- print(title)
- self.fit_and_visualize_estimator(estimator, axes[i], show=False, title=title)
+ def show_or_save_to_file(self, plot_name):
title = self.study.title
- title += '\nMethod: Full Likelihood with Linear marginals and max stable dependency structure'
+ title += '\n' + plot_name
plt.suptitle(title)
- plt.show()
+ if self.show:
+ plt.show()
+ if self.save_to_file:
+ filename = "{}/{}/{}".format(VERSION_TIME, '_'.join(self.study.title.split()), '_'.join(plot_name.split()))
+ filepath = op.join(self.study.result_full_path, filename + '.png')
+ dir = op.dirname(filepath)
+ if not op.exists(dir):
+ os.makedirs(dir, exist_ok=True)
+ plt.savefig(filepath)
def visualize_independent_margin_fits(self, threshold=None, axes=None):
if threshold is None:
@@ -162,8 +172,8 @@ class StudyVisualizer(object):
def df_gev_mle_each_massif(self):
# Fit a margin_fits on each massif
massif_to_gev_mle = {
- massif_name: GevMleFit(self.study.observations_annual_maxima.loc[massif_name]).gev_params.summary_serie
- for massif_name in self.study.safran_massif_names}
+ massif_name: GevMleFit(self.study.observations_annual_maxima.loc[massif_name]).gev_params.summary_serie
+ for massif_name in self.study.safran_massif_names}
return pd.DataFrame(massif_to_gev_mle, columns=self.study.safran_massif_names)
def df_gpd_mle_each_massif(self, threshold):
diff --git a/test/test_meteo_france_SCM_study/__init__.py b/test/test_meteo_france_SCM_study/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/test/test_meteo_france_SCM_study/test_SCM_study.py b/test/test_meteo_france_SCM_study/test_SCM_study.py
index cd4d7f5578a64298579f660d155fc343691fe92a..ec4765c956e62566aa65b1faea02da862495ad1f 100644
--- a/test/test_meteo_france_SCM_study/test_SCM_study.py
+++ b/test/test_meteo_france_SCM_study/test_SCM_study.py
@@ -1,23 +1,12 @@
-import os
import os.path as op
-from collections import OrderedDict
-from typing import List, Dict
+import unittest
-import matplotlib.pyplot as plt
import pandas as pd
-from netCDF4 import Dataset
-
-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 import Safran
-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 get_full_path, cached_property
-import unittest
-from test.test_utils import load_safran_objects
+from experiment.meteo_france_SCM_study.crocus.crocus import ExtendedCrocusSwe
+from experiment.meteo_france_SCM_study.main_visualize import study_iterator
+from experiment.meteo_france_SCM_study.safran.safran import Safran, ExtendedSafran
+from experiment.meteo_france_SCM_study.safran.safran_visualizer import StudyVisualizer
class TestSCMStudy(unittest.TestCase):
@@ -31,6 +20,13 @@ class TestSCMStudy(unittest.TestCase):
# Assert that the massif names are the same between SAFRAN and the coordinate file
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 in study_iterator(study_class, only_first_one=True, both_altitude=True, verbose=False):
+ study_visualizer = StudyVisualizer(study, show=False, save_to_file=False)
+ study_visualizer.visualize_all_kde_graphs()
+ self.assertTrue(True)
+
if __name__ == '__main__':
unittest.main()
diff --git a/test/test_utils.py b/test/test_utils.py
index 2bb92efefa07f84a633241d0ac9c4624381170d4..3cf6be931db34ef29e3539455b7662c49f51c6e9 100644
--- a/test/test_utils.py
+++ b/test/test_utils.py
@@ -46,13 +46,17 @@ def load_smooth_margin_models(coordinates):
return [margin_class(coordinates=coordinates) for margin_class in TEST_MARGIN_TYPES]
-def load_test_max_stable_models():
+def load_test_max_stable_models(only_one_covariance_function=False):
+ default_covariance_function = CovarianceFunction.cauchy
# Load all max stable model
max_stable_models = []
for max_stable_class in TEST_MAX_STABLE_MODEL:
if issubclass(max_stable_class, AbstractMaxStableModelWithCovarianceFunction):
- max_stable_models.extend([max_stable_class(covariance_function=covariance_function)
- for covariance_function in CovarianceFunction])
+ if only_one_covariance_function:
+ max_stable_models.append(max_stable_class(covariance_function=default_covariance_function))
+ else:
+ max_stable_models.extend([max_stable_class(covariance_function=covariance_function)
+ for covariance_function in CovarianceFunction])
else:
max_stable_models.append(max_stable_class())
return max_stable_models
diff --git a/utils.py b/utils.py
index d3e9071ef04240e3a1d9a6a826ee5a2b829faa60..92daa944111b97e75ddbee7c72c7403285d48b55 100644
--- a/utils.py
+++ b/utils.py
@@ -2,7 +2,9 @@ import datetime
import os.path as op
VERSION = datetime.datetime.now()
-
+VERSION_TIME = str(VERSION).split('.')[0]
+for c in [' ', ':', '-']:
+ VERSION_TIME = VERSION_TIME.replace(c, '_')
def get_root_path() -> str:
return op.dirname(op.abspath(__file__))
@@ -62,3 +64,4 @@ if __name__ == '__main__':
e = Example()
print(e.big_attribute)
print(e.big_attribute)
+ print(VERSION_TIME)