diff --git a/experiment/meteo_france_SCM_study/abstract_study.py b/experiment/meteo_france_SCM_study/abstract_study.py
index 264d1fbb1a219c17f58b3e851d770a92d1118f1a..7fc67358c735f921e7a77bdd0bf3c6f08686e02a 100644
--- a/experiment/meteo_france_SCM_study/abstract_study.py
+++ b/experiment/meteo_france_SCM_study/abstract_study.py
@@ -34,7 +34,7 @@ class AbstractStudy(object):
@property
def df_all_snowfall_concatenated(self) -> pd.DataFrame:
- df_list = [pd.DataFrame(snowfall, columns=self.safran_massif_names) for snowfall in
+ df_list = [pd.DataFrame(time_serie, columns=self.safran_massif_names) for time_serie in
self.year_to_daily_time_serie.values()]
df_concatenated = pd.concat(df_list)
return df_concatenated
@@ -111,7 +111,6 @@ class AbstractStudy(object):
""" Visualization methods """
def visualize(self, ax=None, massif_name_to_fill_kwargs=None, show=True, fill=True):
- print("here")
if ax is None:
ax = plt.gca()
df_massif = pd.read_csv(op.join(self.map_full_path, 'massifsalpes.csv'))
@@ -135,6 +134,14 @@ class AbstractStudy(object):
""" Some properties """
+ @property
+ def title(self):
+ return "{} at altitude {}m".format(self.variable_name, self.altitude)
+
+ @property
+ def variable_name(self):
+ return self.variable_class.NAME
+
@property
def relative_path(self) -> str:
return r'local/spatio_temporal_datasets'
diff --git a/experiment/meteo_france_SCM_study/abstract_variable.py b/experiment/meteo_france_SCM_study/abstract_variable.py
index acc653c78c5916291296296cd6af89d0018c03a7..777d7f21f4d9c623e275eb2cfaf3146ba669161d 100644
--- a/experiment/meteo_france_SCM_study/abstract_variable.py
+++ b/experiment/meteo_france_SCM_study/abstract_variable.py
@@ -2,6 +2,8 @@
class AbstractVariable(object):
+ NAME = ''
+
def __init__(self, dataset):
self.dataset = dataset
diff --git a/experiment/meteo_france_SCM_study/crocus/crocus_variables.py b/experiment/meteo_france_SCM_study/crocus/crocus_variables.py
index 1825a43791648c6049cd6e711d0901b238dd4a11..5ba0707ef3925da01d631a1352dd44621c425ff5 100644
--- a/experiment/meteo_france_SCM_study/crocus/crocus_variables.py
+++ b/experiment/meteo_france_SCM_study/crocus/crocus_variables.py
@@ -16,12 +16,14 @@ class CrocusVariable(AbstractVariable):
class CrocusSweVariable(CrocusVariable):
+ NAME = 'Snow Water Equivalent'
def __init__(self, dataset):
super().__init__(dataset, 'SNOWSWE')
class CrocusDepthVariable(CrocusVariable):
+ NAME = 'Snow Depth'
def __init__(self, dataset):
super().__init__(dataset, "SNOWDEPTH")
diff --git a/experiment/meteo_france_SCM_study/main_visualize.py b/experiment/meteo_france_SCM_study/main_visualize.py
index dc249ce215e4791b7f5c5a48a67767915c092bb2..d7203624bd618ca1a4919a5b496cee67879b5ff0 100644
--- a/experiment/meteo_france_SCM_study/main_visualize.py
+++ b/experiment/meteo_france_SCM_study/main_visualize.py
@@ -9,7 +9,7 @@ from experiment.meteo_france_SCM_study.safran.safran_visualizer import StudyVisu
def load_all_studies(study_class, only_first_one=False):
all_studies = []
is_safran_study = study_class == Safran
- nb_days = [1, 5] if is_safran_study else [1]
+ 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)
@@ -20,9 +20,10 @@ def load_all_studies(study_class, only_first_one=False):
if __name__ == '__main__':
- for study_class in [Safran, CrocusSwe, CrocusDepth][:]:
+ for study_class in [Safran, CrocusSwe, CrocusDepth][1:2]:
for study in load_all_studies(study_class, only_first_one=True):
study_visualizer = StudyVisualizer(study)
- # safran_visualizer.visualize_independent_margin_fits(threshold=[None, 20, 40, 60][0])
- study_visualizer.visualize_smooth_margin_fit()
- # safran_visualizer.visualize_full_fit()
+ # study_visualizer.visualize_independent_margin_fits(threshold=[None, 20, 40, 60][0])
+ # study_visualizer.visualize_smooth_margin_fit()
+ study_visualizer.visualize_all_kde_graphs()
+ # study_visualizer.visualize_full_fit()
diff --git a/experiment/meteo_france_SCM_study/safran/safran.py b/experiment/meteo_france_SCM_study/safran/safran.py
index 81b463445c4370b667f7e2aab1c72f725d0a8420..b485bafd1d5eae36a9ce44ab9e0940a8f19175c8 100644
--- a/experiment/meteo_france_SCM_study/safran/safran.py
+++ b/experiment/meteo_france_SCM_study/safran/safran.py
@@ -12,3 +12,8 @@ class Safran(AbstractStudy):
def instantiate_variable_object(self, dataset) -> AbstractVariable:
return self.variable_class(dataset, self.nb_days_of_snowfall)
+
+ @property
+ def variable_name(self):
+ return super().variable_name() + 'cumulated over {} days'.format(self.nb_days_of_snowfall)
+
diff --git a/experiment/meteo_france_SCM_study/safran/safran_snowfall_variable.py b/experiment/meteo_france_SCM_study/safran/safran_snowfall_variable.py
index 04b25737727e049d82fd0900b056260af4a2f47f..491dd3d230db7c386adf35e44d4bba1b6fac72a0 100644
--- a/experiment/meteo_france_SCM_study/safran/safran_snowfall_variable.py
+++ b/experiment/meteo_france_SCM_study/safran/safran_snowfall_variable.py
@@ -19,6 +19,8 @@ class SafranSnowfallVariable(AbstractVariable):
(but here the problem might be that the x_i are not idnependent, they are highly dependent one from another)
"""
+ NAME = 'Snowfall'
+
def __init__(self, dataset, nb_consecutive_days_of_snowfall=1):
super().__init__(dataset)
self.nb_consecutive_days_of_snowfall = nb_consecutive_days_of_snowfall
@@ -46,3 +48,5 @@ class SafranSnowfallVariable(AbstractVariable):
+
+
diff --git a/experiment/meteo_france_SCM_study/safran/safran_visualizer.py b/experiment/meteo_france_SCM_study/safran/safran_visualizer.py
index 517ea10aae1db6e809b7abf6b0f7ca035a755169..cbc0359481023446180de8dad0f062a988d3efca 100644
--- a/experiment/meteo_france_SCM_study/safran/safran_visualizer.py
+++ b/experiment/meteo_france_SCM_study/safran/safran_visualizer.py
@@ -1,3 +1,6 @@
+import math
+import numpy as np
+
import matplotlib.pyplot as plt
import pandas as pd
@@ -11,7 +14,6 @@ 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.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
@@ -21,6 +23,7 @@ class StudyVisualizer(object):
def __init__(self, study: AbstractStudy, show=True):
self.study = study
self.show = show
+ self.window_size_for_smoothing = 21
@property
def observations(self):
@@ -34,6 +37,60 @@ class StudyVisualizer(object):
def dataset(self):
return AbstractDataset(self.observations, self.coordinates)
+ def visualize_all_kde_graphs(self, show=True):
+ 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.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()
+
+ def visualize_kde_graph(self, ax, i, massif_name):
+ self.maxima_plot(ax, i)
+ self.mean_plot(ax, i)
+ ax.set_xlabel('year')
+ ax.set_title(massif_name)
+
+ def mean_plot(self, ax, i):
+ # Display the mean graph
+ # Counting the sum of 3-consecutive days of snowfall does not have any physical meaning,
+ # as we are counting twice some days
+ color_mean = 'g'
+ tuples_x_y = [(year, np.mean(data[:, i])) for year, data in self.study.year_to_daily_time_serie.items()]
+ x, y = list(zip(*tuples_x_y))
+ x, y = self.smooth(x, y)
+ ax.plot(x, y, color=color_mean)
+ ax.set_ylabel('mean', color=color_mean)
+
+ def maxima_plot(self, ax, i):
+ # Display the graph of the max on top
+ color_maxima = 'r'
+ tuples_x_y = [(year, annual_maxima[i]) for year, annual_maxima in self.study.year_to_annual_maxima.items()]
+ x, y = list(zip(*tuples_x_y))
+ ax2 = ax.twinx()
+ ax2.plot(x, y, color=color_maxima)
+ ax2.set_ylabel('maxima', color=color_maxima)
+
+ def smooth(self, x, y):
+ # Average on windows of size 2*M+1 (M elements on each side)
+ filter = np.ones(self.window_size_for_smoothing) / self.window_size_for_smoothing
+ y = np.convolve(y, filter, mode='valid')
+ assert self.window_size_for_smoothing % 2 == 1
+ nb_to_delete = int(self.window_size_for_smoothing // 2)
+ x = np.array(x)[nb_to_delete:-nb_to_delete]
+ assert len(x) == len(y)
+ return x, y
+
+
+
def fit_and_visualize_estimator(self, estimator):
estimator.fit()
axes = estimator.margin_function_fitted.visualize(show=False)