diff --git a/extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py b/extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py
index d7caa81ae9e855799fb2b0140df9ba73dd63823d..71b0e606139925936521a27fc21b361f4f29f785 100644
--- a/extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py
+++ b/extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py
@@ -13,6 +13,8 @@ import numpy as np
import pandas as pd
import seaborn as sns
+from extreme_data.meteo_france_data.adamont_data.abstract_adamont_study import AbstractAdamontStudy
+from extreme_data.meteo_france_data.adamont_data.adamont_scenario import gcm_rcm_couple_to_str
from extreme_data.meteo_france_data.scm_models_data.visualization.plot_utils import load_plot
from extreme_fit.estimator.margin_estimator.utils import fitted_stationary_gev
from extreme_fit.model.margin_model.utils import fitmethod_to_str
@@ -541,6 +543,10 @@ class StudyVisualizer(VisualizationParameters):
def show_or_save_to_file(self, add_classic_title=True, no_title=False, tight_layout=False, tight_pad=None,
dpi=None, folder_for_variable=True):
+ if isinstance(self.study, AbstractAdamontStudy):
+ prefix = gcm_rcm_couple_to_str(self.study.gcm_rcm_couple)
+ self.plot_name = prefix + ' ' + self.plot_name
+
assert self.plot_name is not None
if add_classic_title:
title = self.study.title
diff --git a/projects/altitude_spatial_model/altitudes_fit/altitudes_studies.py b/projects/altitude_spatial_model/altitudes_fit/altitudes_studies.py
index b0566196a1a68992f641a9cb443f06211418aca8..e5162707dce67d36bfad0a19337e26adff32cf99 100644
--- a/projects/altitude_spatial_model/altitudes_fit/altitudes_studies.py
+++ b/projects/altitude_spatial_model/altitudes_fit/altitudes_studies.py
@@ -5,7 +5,7 @@ from collections import OrderedDict
from cached_property import cached_property
from extreme_data.meteo_france_data.adamont_data.abstract_adamont_study import AbstractAdamontStudy
-from extreme_data.meteo_france_data.adamont_data.adamont_scenario import scenario_to_str
+from extreme_data.meteo_france_data.adamont_data.adamont_scenario import scenario_to_str, gcm_rcm_couple_to_str
from extreme_data.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
from extreme_data.meteo_france_data.scm_models_data.visualization.main_study_visualizer import \
SCM_STUDY_CLASS_TO_ABBREVIATION, STUDY_CLASS_TO_ABBREVIATION
diff --git a/projects/projected_snowfall/elevation_temporal_model_for_projections/ensemble_fit/abstract_ensemble_fit.py b/projects/projected_snowfall/elevation_temporal_model_for_projections/ensemble_fit/abstract_ensemble_fit.py
index d2c9cf659961692b111cd34a62f710983cd6644a..b378302178c941c3a84b651b235905ccab44912c 100644
--- a/projects/projected_snowfall/elevation_temporal_model_for_projections/ensemble_fit/abstract_ensemble_fit.py
+++ b/projects/projected_snowfall/elevation_temporal_model_for_projections/ensemble_fit/abstract_ensemble_fit.py
@@ -1,4 +1,27 @@
+from typing import Dict, Tuple, List
+
+from extreme_fit.model.margin_model.utils import MarginFitMethod
+from projects.altitude_spatial_model.altitudes_fit.altitudes_studies import AltitudesStudies
+from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitude_group import DefaultAltitudeGroup
class AbstractEnsembleFit(object):
- pass
\ No newline at end of file
+
+ def __init__(self, massif_names, gcm_rcm_couple_to_altitude_studies: Dict[Tuple[str, str], AltitudesStudies],
+ models_classes,
+ fit_method=MarginFitMethod.extremes_fevd_mle,
+ temporal_covariate_for_fit=None,
+ only_models_that_pass_goodness_of_fit_test=True,
+ confidence_interval_based_on_delta_method=False,
+ ):
+ self.massif_names = massif_names
+ self.models_classes = models_classes
+ self.gcm_rcm_couple_to_altitude_studies = gcm_rcm_couple_to_altitude_studies
+ self.fit_method = fit_method
+ self.temporal_covariate_for_fit = temporal_covariate_for_fit
+ self.only_models_that_pass_goodness_of_fit_test = only_models_that_pass_goodness_of_fit_test
+ self.confidence_interval_based_on_delta_method = confidence_interval_based_on_delta_method
+
+
+ def plot(self):
+ raise NotImplementedError
\ No newline at end of file
diff --git a/projects/projected_snowfall/elevation_temporal_model_for_projections/ensemble_fit/independent_ensemble_fit.py b/projects/projected_snowfall/elevation_temporal_model_for_projections/ensemble_fit/independent_ensemble_fit.py
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..9de456cc17bb4685981ef3a36d4ab5b72fa6814e 100644
--- a/projects/projected_snowfall/elevation_temporal_model_for_projections/ensemble_fit/independent_ensemble_fit.py
+++ b/projects/projected_snowfall/elevation_temporal_model_for_projections/ensemble_fit/independent_ensemble_fit.py
@@ -0,0 +1,48 @@
+from typing import Dict, Tuple, List
+
+from extreme_fit.model.margin_model.utils import MarginFitMethod
+from projects.altitude_spatial_model.altitudes_fit.altitudes_studies import AltitudesStudies
+from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitude_group import DefaultAltitudeGroup
+from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitudes_studies_visualizer_for_non_stationary_models import \
+ AltitudesStudiesVisualizerForNonStationaryModels
+from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.one_fold_fit import OneFoldFit
+from projects.altitude_spatial_model.altitudes_fit.utils_altitude_studies_visualizer import compute_and_assign_max_abs
+from projects.projected_snowfall.elevation_temporal_model_for_projections.ensemble_fit.abstract_ensemble_fit import \
+ AbstractEnsembleFit
+
+
+class IndependentEnsembleFit(AbstractEnsembleFit):
+ """For each gcm_rcm_couple, we create a OneFoldFit"""
+
+ def __init__(self, massif_names, gcm_rcm_couple_to_altitude_studies: Dict[Tuple[str, str], AltitudesStudies], models_classes,
+ fit_method=MarginFitMethod.extremes_fevd_mle, temporal_covariate_for_fit=None, only_models_that_pass_goodness_of_fit_test=True,
+ confidence_interval_based_on_delta_method=False):
+ super().__init__(massif_names, gcm_rcm_couple_to_altitude_studies, models_classes, fit_method, temporal_covariate_for_fit, only_models_that_pass_goodness_of_fit_test,
+ confidence_interval_based_on_delta_method)
+
+ # Load a classical visualizer
+ self.gcm_rcm_couple_to_visualizer = {}
+ for gcm_rcm_couple, studies in gcm_rcm_couple_to_altitude_studies.items():
+ visualizer = AltitudesStudiesVisualizerForNonStationaryModels(studies, self.models_classes,
+ False,
+ self.massif_names, self.fit_method,
+ self.temporal_covariate_for_fit,
+ self.only_models_that_pass_goodness_of_fit_test,
+ self.confidence_interval_based_on_delta_method)
+ self.gcm_rcm_couple_to_visualizer[gcm_rcm_couple] = visualizer
+
+ # Assign max
+ visualizer_list = list(self.gcm_rcm_couple_to_visualizer.values())
+ compute_and_assign_max_abs(visualizer_list)
+
+ def plot(self):
+ for v in self.gcm_rcm_couple_to_visualizer.values():
+ self.plot_one_visualizer(v)
+
+ @staticmethod
+ def plot_one_visualizer(visualizer):
+ # visualizer.studies.plot_maxima_time_series(['Vanoise'])
+ visualizer.plot_shape_map()
+ visualizer.plot_moments()
+ # visualizer.plot_qqplots()
+
diff --git a/projects/projected_snowfall/elevation_temporal_model_for_projections/main_elevation_temporal_for_projections_ensemble.py b/projects/projected_snowfall/elevation_temporal_model_for_projections/main_elevation_temporal_for_projections_ensemble.py
index f3f05a462abbd40b01d52f3088aef06e52721f55..33baed2b2ba1335d1e9bca5215491f0acf9acd0e 100644
--- a/projects/projected_snowfall/elevation_temporal_model_for_projections/main_elevation_temporal_for_projections_ensemble.py
+++ b/projects/projected_snowfall/elevation_temporal_model_for_projections/main_elevation_temporal_for_projections_ensemble.py
@@ -2,12 +2,25 @@ import datetime
import time
from typing import List
+import matplotlib as mpl
+
+
+mpl.rcParams['text.usetex'] = True
+mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
+
import matplotlib
+from extreme_fit.model.utils import set_seed_for_test
+from extreme_data.meteo_france_data.adamont_data.adamont.adamont_snowfall import AdamontSnowfall
+from extreme_data.meteo_france_data.adamont_data.adamont_scenario import AdamontScenario, get_gcm_rcm_couples
+from projects.projected_snowfall.elevation_temporal_model_for_projections.ensemble_fit.independent_ensemble_fit import \
+ IndependentEnsembleFit
from projects.projected_snowfall.elevation_temporal_model_for_projections.utils_projected_visualizer import \
load_projected_visualizer_list
from projects.projected_snowfall.elevation_temporal_model_for_projections.visualizer_for_projection_ensemble import \
VisualizerForProjectionEnsemble
+from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_covariate_for_fit import \
+ AnomalyTemperatureTemporalCovariate
matplotlib.use('Agg')
@@ -17,96 +30,60 @@ from projects.altitude_spatial_model.altitudes_fit.plots.plot_histogram_altitude
from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_extract_eurocode_return_level import \
AbstractExtractEurocodeReturnLevel
-import matplotlib as mpl
-from extreme_fit.model.utils import set_seed_for_test
-
-mpl.rcParams['text.usetex'] = True
-mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitude_group import altitudes_for_groups
-from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day, SafranSnowfall3Days, \
- SafranSnowfall5Days, SafranSnowfall7Days
from extreme_data.meteo_france_data.scm_models_data.utils import Season
def main():
- study_classes = [SafranSnowfall1Day
- , SafranSnowfall3Days,
- SafranSnowfall5Days, SafranSnowfall7Days][:1]
- seasons = [Season.annual, Season.winter, Season.spring, Season.automn][:1]
-
+ study_classes = [AdamontSnowfall][:1]
+ scenario = AdamontScenario.rcp85
+ gcm_rcm_couples = get_gcm_rcm_couples(scenario)[:2]
+ ensemble_fit_class = [IndependentEnsembleFit]
+ temporal_covariate_for_fit = [None, AnomalyTemperatureTemporalCovariate][0]
set_seed_for_test()
- model_must_pass_the_test = False
AbstractExtractEurocodeReturnLevel.ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY = 0.2
- fast = False
+ fast = True
if fast is None:
massif_names = None
AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP = 10
altitudes_list = altitudes_for_groups[2:3]
elif fast:
AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP = 10
- massif_names = ['Vanoise', 'Haute-Maurienne', 'Vercors'][:1]
+ massif_names = ['Vanoise'][:]
altitudes_list = altitudes_for_groups[1:2]
else:
massif_names = None
altitudes_list = altitudes_for_groups[:]
start = time.time()
- main_loop(altitudes_list, massif_names, seasons, study_classes, model_must_pass_the_test)
+ main_loop(gcm_rcm_couples, altitudes_list, massif_names, study_classes, ensemble_fit_class, scenario, temporal_covariate_for_fit)
end = time.time()
duration = str(datetime.timedelta(seconds=end - start))
print('Total duration', duration)
-def main_loop(altitudes_list, massif_names, seasons, study_classes, model_must_pass_the_test):
+def main_loop(gcm_rcm_couples, altitudes_list, massif_names, study_classes, ensemble_fit_classes, scenario, temporal_covariate_for_fit):
assert isinstance(altitudes_list, List)
assert isinstance(altitudes_list[0], List)
- for season in seasons:
- for study_class in study_classes:
- print('Inner loop', season, study_class)
- visualizer_list = load_projected_visualizer_list(season, study_class, altitudes_list, massif_names,
- model_must_pass_the_test
- )
- plot_visualizers(massif_names, visualizer_list)
- for visualizer in visualizer_list:
- plot_visualizer(massif_names, visualizer)
- del visualizer_list
- time.sleep(2)
-
-
-def plot_visualizers(massif_names, visualizer_list):
- # plot_histogram_all_models_against_altitudes(massif_names, visualizer_list)
- plot_histogram_all_trends_against_altitudes(massif_names, visualizer_list)
- # plot_shoe_plot_ratio_interval_size_against_altitude(massif_names, visualizer_list)
- for relative in [True, False]:
- plot_shoe_plot_changes_against_altitude(massif_names, visualizer_list, relative=relative)
- # plot_coherence_curves(massif_names, visualizer_list)
- # plot_coherence_curves(['Vanoise'], visualizer_list)
- pass
-
-
-def plot_visualizer(massif_names, visualizer):
- # Plot time series
- # visualizer.studies.plot_maxima_time_series(massif_names)
- # visualizer.studies.plot_maxima_time_series(['Vanoise'])
-
- # Plot the results for the model that minimizes the individual aic
- plots(visualizer)
-
- # Plot the results for the model that minimizes the total aic
- # plot_total_aic(model_classes, visualizer)
- pass
-
-
-def plots(visualizer: VisualizerForProjectionEnsemble):
- # visualizer.plot_shape_map()
- visualizer.plot_moments()
- # visualizer.plot_qqplots()
- # for std in [True, False]:
- # visualizer.studies.plot_mean_maxima_against_altitude(std=std)
+ for study_class in study_classes:
+ print('Inner loop', study_class)
+ visualizer_list = load_projected_visualizer_list(gcm_rcm_couples=gcm_rcm_couples, ensemble_fit_classes=ensemble_fit_classes,
+ season=Season.annual, study_class=study_class,
+ altitudes_list=altitudes_list, massif_names=massif_names,
+ scenario=scenario,
+ temporal_covariate_for_fit=temporal_covariate_for_fit)
+ for v in visualizer_list:
+ v.plot()
+
+ del visualizer_list
+ time.sleep(2)
+
+
+
if __name__ == '__main__':
diff --git a/projects/projected_snowfall/elevation_temporal_model_for_projections/utils_projected_visualizer.py b/projects/projected_snowfall/elevation_temporal_model_for_projections/utils_projected_visualizer.py
index f3bd8d953e60b4beb0f3bc2dab3af272d85fc919..a7f70c95ac4c40fc00426890cd4fec7c71e8fc1a 100644
--- a/projects/projected_snowfall/elevation_temporal_model_for_projections/utils_projected_visualizer.py
+++ b/projects/projected_snowfall/elevation_temporal_model_for_projections/utils_projected_visualizer.py
@@ -1,44 +1,38 @@
+from extreme_data.meteo_france_data.adamont_data.adamont_scenario import AdamontScenario, rcp_scenarios
from extreme_fit.model.margin_model.polynomial_margin_model.utils import \
ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS
from projects.altitude_spatial_model.altitudes_fit.altitudes_studies import AltitudesStudies
from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitudes_studies_visualizer_for_non_stationary_models import \
AltitudesStudiesVisualizerForNonStationaryModels
+from projects.projected_snowfall.elevation_temporal_model_for_projections.visualizer_for_projection_ensemble import \
+ VisualizerForProjectionEnsemble
-def load_projected_visualizer_list(season, study_class, altitudes_list, massif_names, model_must_pass_the_test=True, **kwargs_study):
+def load_projected_visualizer_list(gcm_rcm_couples, ensemble_fit_classes,
+ season, study_class, altitudes_list, massif_names, model_must_pass_the_test=False,
+ scenario=AdamontScenario.rcp85,
+ temporal_covariate_for_fit=None, **kwargs_study):
model_classes = ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS
+ assert scenario in rcp_scenarios
visualizer_list = []
# Load all studies
for altitudes in altitudes_list:
- studies = AltitudesStudies(study_class, altitudes, season=season, **kwargs_study)
- visualizer = AltitudesStudiesVisualizerForNonStationaryModels(studies=studies,
- model_classes=model_classes,
- massif_names=massif_names,
- show=False,
- temporal_covariate_for_fit=None,
- confidence_interval_based_on_delta_method=False,
- display_only_model_that_pass_anderson_test=model_must_pass_the_test
- )
+ gcm_rcm_couple_to_altitude_studies = {}
+ for gcm_rcm_couple in gcm_rcm_couples:
+ studies = AltitudesStudies(study_class, altitudes, season=season,
+ scenario=scenario, gcm_rcm_couple=gcm_rcm_couple, **kwargs_study)
+ gcm_rcm_couple_to_altitude_studies[gcm_rcm_couple] = studies
+ visualizer = VisualizerForProjectionEnsemble(gcm_rcm_couple_to_altitude_studies=gcm_rcm_couple_to_altitude_studies,
+ model_classes=model_classes,
+ ensemble_fit_classes=ensemble_fit_classes,
+ massif_names=massif_names,
+ show=False,
+ temporal_covariate_for_fit=temporal_covariate_for_fit,
+ confidence_interval_based_on_delta_method=False,
+ display_only_model_that_pass_gof_test=model_must_pass_the_test
+ )
visualizer_list.append(visualizer)
- compute_and_assign_max_abs(visualizer_list)
return visualizer_list
-def compute_and_assign_max_abs(visualizer_list):
- # Compute the max abs for all metrics
- d = {}
- for method_name in AltitudesStudiesVisualizerForNonStationaryModels.moment_names:
- for order in AltitudesStudiesVisualizerForNonStationaryModels.orders:
- c = (method_name, order)
- max_abs = max([
- max([abs(e) for e in v.method_name_and_order_to_d(method_name, order).values()
- ]) for v in visualizer_list])
- d[c] = max_abs
- # Assign the max abs dictionary
- for v in visualizer_list:
- v._method_name_and_order_to_max_abs = d
- # Compute the max abs for the shape parameter
- max_abs_for_shape = max([max([abs(e) for e in v.massif_name_to_shape.values()]) for v in visualizer_list])
- for v in visualizer_list:
- v._max_abs_for_shape = max_abs_for_shape
diff --git a/projects/projected_snowfall/elevation_temporal_model_for_projections/visualizer_for_projection_ensemble.py b/projects/projected_snowfall/elevation_temporal_model_for_projections/visualizer_for_projection_ensemble.py
index a88cc1a9faf5ad4fe4a7c286f47f946222e02c6b..86ffe40664346da73db1eea2a49cc4f1c28a9438 100644
--- a/projects/projected_snowfall/elevation_temporal_model_for_projections/visualizer_for_projection_ensemble.py
+++ b/projects/projected_snowfall/elevation_temporal_model_for_projections/visualizer_for_projection_ensemble.py
@@ -33,7 +33,7 @@ from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
class VisualizerForProjectionEnsemble(StudyVisualizer):
- def __init__(self, gcm_rcm_couple_to_altitude_studies: Dict[str,AltitudesStudies],
+ def __init__(self, gcm_rcm_couple_to_altitude_studies: Dict[str, AltitudesStudies],
model_classes: List[AbstractSpatioTemporalPolynomialModel],
show=False,
ensemble_fit_classes=None,
@@ -46,513 +46,23 @@ class VisualizerForProjectionEnsemble(StudyVisualizer):
studies = list(gcm_rcm_couple_to_altitude_studies.values())[0]
study = studies.study
super().__init__(study, show=show, save_to_file=not show)
- self.ensemble_fit_classes = ensemble_fit_classes
- self.gcm_rcm_couple_to_altitude_studies = gcm_rcm_couple_to_altitude_studies
- self.non_stationary_models = model_classes
- self.fit_method = fit_method
- self.temporal_covariate_for_fit = temporal_covariate_for_fit
- self.display_only_model_that_pass_test = display_only_model_that_pass_gof_test
- self.massif_names = massif_names if massif_names is not None else study.all_massif_names()
- self.massif_name_to_massif_id = {m: i for i, m in enumerate(self.massif_names)}
- self.altitude_group = get_altitude_group_from_altitudes(studies.altitudes)
- self.confidence_interval_based_on_delta_method = confidence_interval_based_on_delta_method
-
- # Load ensemble_fit mapper
- ensemble_fit_class_to_object = {}
- for ensemble_fit_class in ensemble_fit_classes:
- pass
-
# Load one fold fit
self.massif_name_to_massif_altitudes = {}
- self._massif_name_to_one_fold_fit = {}
- for massif_name in self.massif_names:
- # Load valid massif altitudes
- massif_altitudes = self.get_massif_altitudes(massif_name)
- if self.load_condition(massif_altitudes):
- # Save the massif altitudes only for those who pass the condition
- self.massif_name_to_massif_altitudes[massif_name] = massif_altitudes
- # Load dataset
- dataset = studies.spatio_temporal_dataset(massif_name=massif_name, massif_altitudes=massif_altitudes)
- old_fold_fit = OneFoldFit(massif_name, dataset, model_classes, self.fit_method,
- self.temporal_covariate_for_fit,
- type(self.altitude_group),
- self.display_only_model_that_pass_test,
- self.confidence_interval_based_on_delta_method)
- self._massif_name_to_one_fold_fit[massif_name] = old_fold_fit
- # Print number of massif without any validated fit
- massifs_without_any_validated_fit = [massif_name
- for massif_name, old_fold_fit in self._massif_name_to_one_fold_fit.items()
- if not old_fold_fit.has_at_least_one_valid_model]
- print('Not validated:', len(massifs_without_any_validated_fit), massifs_without_any_validated_fit)
- # Cache
- self._method_name_and_order_to_massif_name_to_value = {}
- self._method_name_and_order_to_max_abs = {}
- self._max_abs_for_shape = None
-
- moment_names = ['moment', 'changes_of_moment', 'relative_changes_of_moment'][:]
- orders = [1, 2, None][2:]
-
- def get_massif_altitudes(self, massif_name):
- valid_altitudes = [altitude for altitude, study in self.studies.altitude_to_study.items()
- if massif_name in study.study_massif_names]
- massif_altitudes = []
- for altitude in valid_altitudes:
- study = self.studies.altitude_to_study[altitude]
- annual_maxima = study.massif_name_to_annual_maxima[massif_name]
- percentage_of_non_zeros = 100 * np.count_nonzero(annual_maxima) / len(annual_maxima)
- if percentage_of_non_zeros > 90:
- massif_altitudes.append(altitude)
- # else:
- # print(massif_name, altitude, percentage_of_non_zeros)
- return massif_altitudes
-
- def load_condition(self, massif_altitudes):
- # At least two altitudes for the estimated
- # reference_altitude_is_in_altitudes = (self.altitude_group.reference_altitude in massif_altitudes)
- at_least_two_altitudes = (len(massif_altitudes) >= 2)
- # return reference_altitude_is_in_altitudes and at_least_two_altitudes
- return at_least_two_altitudes
-
- @property
- def massif_name_to_one_fold_fit(self) -> Dict[str, OneFoldFit]:
- return {massif_name: old_fold_fit for massif_name, old_fold_fit in self._massif_name_to_one_fold_fit.items()
- if not self.display_only_model_that_pass_test
- or old_fold_fit.has_at_least_one_valid_model}
-
- def plot_moments(self):
- for method_name in self.moment_names[:2]:
- for order in self.orders:
- # self.plot_against_years(method_name, order)
- self.plot_map_moment(method_name, order)
-
- def method_name_and_order_to_max_abs(self, method_name, order):
- c = (method_name, order)
- if c not in self._method_name_and_order_to_max_abs:
- return None
- else:
- return self._method_name_and_order_to_max_abs[c]
-
- def method_name_and_order_to_d(self, method_name, order):
- c = (method_name, order)
- if c not in self._method_name_and_order_to_massif_name_to_value:
- # Compute values
- massif_name_to_value = {}
- for massif_name, one_fold_fit in self.massif_name_to_one_fold_fit.items():
- value = \
- one_fold_fit.__getattribute__(method_name)([self.altitude_group.reference_altitude], order=order)[0]
- massif_name_to_value[massif_name] = value
- # Remove values
- if any([np.isinf(v) for v in massif_name_to_value.values()]):
- print("shape to large > 0.5, thus removing std that are infinite")
- massif_name_to_value = {m: v for m, v in massif_name_to_value.items()
- if not np.isinf(v)}
- # Store it
- self._method_name_and_order_to_massif_name_to_value[c] = massif_name_to_value
- return self._method_name_and_order_to_massif_name_to_value[c]
-
- def ratio_groups(self):
- return [self.ratio_uncertainty_interval_size(altitude, 2019) for altitude in self.studies.altitudes]
-
- def ratio_uncertainty_interval_size(self, altitude, year):
- study = self.studies.altitude_to_study[altitude]
- massif_name_to_interval = study.massif_name_to_stationary_gev_params_and_confidence(OneFoldFit.quantile_level,
- self.confidence_interval_based_on_delta_method)[
- 1]
- massif_names_with_pointwise_interval = set(massif_name_to_interval)
- valid_massif_names = set(self.massif_name_to_one_fold_fit.keys())
- intersection_massif_names = valid_massif_names.intersection(massif_names_with_pointwise_interval)
- ratios = []
- for massif_name in intersection_massif_names:
- one_fold_fit = self.massif_name_to_one_fold_fit[massif_name]
- new_interval_size = one_fold_fit.best_confidence_interval(altitude, year).interval_size
- old_interval_size = massif_name_to_interval[massif_name].interval_size
- ratio = new_interval_size / old_interval_size
- ratios.append(ratio)
- return ratios
-
- def plot_map_moment(self, method_name, order):
- massif_name_to_value = self.method_name_and_order_to_d(method_name, order)
- # Plot settings
- moment = ' '.join(method_name.split('_'))
- str_for_last_year = ' in 2019'
- moment = moment.replace('moment', '{}{}'.format(OneFoldFit.get_moment_str(order=order), str_for_last_year))
- plot_name = '{}{} '.format(OneFoldFit.folder_for_plots, moment)
-
-
- if 'change' in method_name:
- plot_name = plot_name.replace(str_for_last_year, '')
- plot_name += ' between {} and {}'.format(2019 - OneFoldFit.nb_years, 2019)
- if 'relative' not in method_name:
- # Put the relative score as text on the plot for the change.
- massif_name_to_text = {m: ('+' if v > 0 else '') + str(int(v)) + '\%' for m, v in
- self.method_name_and_order_to_d(self.moment_names[2], order).items()}
-
- parenthesis = self.study.variable_unit if 'relative' not in method_name else '\%'
- ylabel = '{} ({})'.format(plot_name, parenthesis)
-
- max_abs_change = self.method_name_and_order_to_max_abs(method_name, order)
- add_colorbar = self.add_colorbar
-
- is_return_level_plot = (self.moment_names.index(method_name) == 0) and (order is None)
- if is_return_level_plot:
-
- cmap = plt.cm.Spectral
- cmap = remove_the_extreme_colors(cmap, epsilon=0.25)
- cmap = get_inverse_colormap(cmap)
- add_colorbar = True
- max_abs_change = None
- massif_name_to_text = {m: round(v) for m, v in massif_name_to_value.items()}
- graduation = self.altitude_group.graduation_for_return_level
- fontsize_label = 17
- else:
- # cmap = plt.cm.RdYlGn
- cmap = [plt.cm.coolwarm, plt.cm.bwr, plt.cm.seismic][1]
- # cmap = get_inverse_colormap(cmap)
- # cmap = get_cmap_with_inverted_blue_and_green_channels(cmap)
- cmap = remove_the_extreme_colors(cmap)
- graduation = 10
- fontsize_label = 10
-
- negative_and_positive_values = self.moment_names.index(method_name) > 0
- # Plot the map
-
- self.plot_map(cmap=cmap, graduation=graduation,
- label=ylabel, massif_name_to_value=massif_name_to_value,
- plot_name=plot_name, add_x_label=True,
- negative_and_positive_values=negative_and_positive_values,
- altitude=self.altitude_group.reference_altitude,
- add_colorbar=add_colorbar,
- max_abs_change=max_abs_change,
- massif_name_to_text=massif_name_to_text,
- xlabel=self.altitude_group.xlabel,
- fontsize_label=fontsize_label,
- )
-
-
- @property
- def add_colorbar(self):
- # return isinstance(self.altitude_group, (VeyHighAltitudeGroup))
- return isinstance(self.altitude_group, (VeyHighAltitudeGroup, MidAltitudeGroup))
-
- def plot_against_years(self, method_name, order):
- ax = plt.gca()
- min_altitude, *_, max_altitude = self.studies.altitudes
- altitudes_plot = np.linspace(min_altitude, max_altitude, num=50)
- for massif_name, one_fold_fit in self.massif_name_to_one_fold_fit.items():
- massif_altitudes = self.studies.massif_name_to_altitudes[massif_name]
- ind = (min(massif_altitudes) <= altitudes_plot) & (altitudes_plot <= max(massif_altitudes))
- massif_altitudes_plot = altitudes_plot[ind]
- values = one_fold_fit.__getattribute__(method_name)(massif_altitudes_plot, order=order)
- massif_id = self.massif_name_to_massif_id[massif_name]
- plot_against_altitude(massif_altitudes_plot, ax, massif_id, massif_name, values)
- # Plot settings
- ax.legend(prop={'size': 7}, ncol=3)
- moment = ' '.join(method_name.split('_'))
- moment = moment.replace('moment', '{} in 2019'.format(OneFoldFit.get_moment_str(order=order)))
- plot_name = '{}Model {} annual maxima of {}'.format(OneFoldFit.folder_for_plots, moment,
- SCM_STUDY_CLASS_TO_ABBREVIATION[self.studies.study_class])
- ax.set_ylabel('{} ({})'.format(plot_name, self.study.variable_unit), fontsize=15)
- ax.set_xlabel('altitudes', fontsize=15)
- ax.tick_params(axis='both', which='major', labelsize=13)
- self.studies.show_or_save_to_file(plot_name=plot_name, show=self.show, no_title=True)
- ax.clear()
-
- # def plot_abstract_fast(self, massif_name_to_value, label, graduation=10.0, cmap=plt.cm.coolwarm, add_x_label=True,
- # negative_and_positive_values=True, massif_name_to_text=None):
- # plot_name = '{}{}'.format(OneFoldFit.folder_for_plots, label)
- # self.plot_map(cmap, self.fit_method, graduation, label, massif_name_to_value, plot_name, add_x_label,
- # negative_and_positive_values,
- # massif_name_to_text)
-
- @property
- def massif_name_to_shape(self):
- return {massif_name: one_fold_fit.best_shape
- for massif_name, one_fold_fit in self.massif_name_to_one_fold_fit.items()}
-
- @property
- def massif_name_to_best_name(self):
- return {massif_name: one_fold_fit.best_name
- for massif_name, one_fold_fit in self.massif_name_to_one_fold_fit.items()}
-
- def plot_best_coef_maps(self):
- for param_name in GevParams.PARAM_NAMES:
- coordinate_names = [AbstractCoordinates.COORDINATE_X, AbstractCoordinates.COORDINATE_T]
- dim_to_coordinate_name = dict(zip([0, 1], coordinate_names))
- for dim in [0, 1, (0, 1)]:
- coordinate_name = LinearCoef.coefficient_name(dim, dim_to_coordinate_name)
- for degree in range(4):
- coef_name = ' '.join([param_name + coordinate_name + str(degree)])
- massif_name_to_best_coef = {}
- for massif_name, one_fold_fit in self.massif_name_to_one_fold_fit.items():
- best_coef = one_fold_fit.best_coef(param_name, dim, degree)
- if best_coef is not None:
- massif_name_to_best_coef[massif_name] = best_coef
-
- if len(massif_name_to_best_coef) > 0:
- for evaluate_coordinate in [False, True][:1]:
- if evaluate_coordinate:
- coef_name += 'evaluated at coordinates'
- for massif_name in massif_name_to_best_coef.keys():
- if AbstractCoordinates.COORDINATE_X in coordinate_name:
- massif_name_to_best_coef[massif_name] *= np.power(1000, degree)
- if AbstractCoordinates.COORDINATE_T in coordinate_name:
- massif_name_to_best_coef[massif_name] *= np.power(2019, degree)
- self.plot_best_coef_map(coef_name.replace('_', ''), massif_name_to_best_coef)
-
- def plot_best_coef_map(self, coef_name, massif_name_to_best_coef):
- values = list(massif_name_to_best_coef.values())
- graduation = (max(values) - min(values)) / 6
- print(coef_name, graduation, max(values), min(values))
- negative_and_positive_values = (max(values) > 0) and (min(values) < 0)
- raise NotImplementedError
- self.plot_map(massif_name_to_value=massif_name_to_best_coef,
- label='{}Coef/{} plot for {} {}'.format(OneFoldFit.folder_for_plots,
- coef_name,
- SCM_STUDY_CLASS_TO_ABBREVIATION[
- type(self.study)],
- self.study.variable_unit),
- add_x_label=False, graduation=graduation, massif_name_to_text=self.massif_name_to_best_name,
- negative_and_positive_values=negative_and_positive_values)
-
- def plot_shape_map(self):
+ self.ensemble_class_to_ensemble_fit = {}
+ for ensemble_fit_class in ensemble_fit_classes:
+ ensemble_fit = ensemble_fit_class(massif_names, gcm_rcm_couple_to_altitude_studies, model_classes,
+ fit_method, temporal_covariate_for_fit,
+ display_only_model_that_pass_gof_test,
+ confidence_interval_based_on_delta_method)
+ self.ensemble_class_to_ensemble_fit[ensemble_fit_class] = ensemble_fit
- label = 'Shape parameter in 2019 (no unit)'
- max_abs_change = self._max_abs_for_shape + 0.05
- self.plot_map(massif_name_to_value=self.massif_name_to_shape,
- label=label,
- plot_name=label,
- fontsize_label=15,
- add_x_label=True, graduation=0.1,
- massif_name_to_text=self.massif_name_to_best_name,
- cmap=matplotlib.cm.get_cmap('BrBG_r'),
- altitude=self.altitude_group.reference_altitude,
- add_colorbar=self.add_colorbar,
- max_abs_change=max_abs_change,
- xlabel=self.altitude_group.xlabel,
- )
+ def plot(self):
+ self.plot_independent()
+ self.plot_dependent()
- def plot_altitude_for_the_peak(self):
+ def plot_dependent(self):
pass
- def plot_year_for_the_peak(self, plot_mean=True):
- t_list = self.study.ordered_years
- return_period = 50
- for massif_name, one_fold_fit in self.massif_name_to_one_fold_fit.items():
- ax = plt.gca()
- # One plot for each altitude
- altitudes = np.arange(500, min(3000, max(self.studies.altitudes)), 500)
- for altitude in altitudes:
- i = 0
- while self.studies.altitudes[i] < altitude:
- i += 1
- nearest_altitude = self.studies.altitudes[i]
- nearest_study = self.studies.altitude_to_study[nearest_altitude]
- if massif_name in nearest_study.study_massif_names:
- y_list = []
- for t in t_list:
- coordinate = np.array([altitude, t])
- gev_params = one_fold_fit.best_function_from_fit.get_params(coordinate, is_transformed=False)
- if plot_mean:
- y = gev_params.mean
- else:
- y = gev_params.return_level(return_period=return_period)
- y_list.append(y)
- label = '{} m'.format(altitude)
- ax.plot(t_list, y_list, label=label)
- ax.legend()
- # Modify the limits of the y axis
- lim_down, lim_up = ax.get_ylim()
- ax_lim = (0, lim_up)
- ax.set_ylim(ax_lim)
- ax.set_xlabel('Year')
- if plot_mean:
- ylabel = 'Mean {} maxima'.format(self.study.season_name)
- else:
- ylabel = '{}-year return level'.format(return_period)
- ax.set_ylabel('{} of {} in {} ({})'.format(ylabel, SCM_STUDY_CLASS_TO_ABBREVIATION[type(self.study)],
- massif_name.replace('_', ' '), self.study.variable_unit))
- peak_year_folder = 'Peak year ' + ylabel
- plot_name = '{}{}/Peak year for {}'.format(OneFoldFit.folder_for_plots, peak_year_folder,
- massif_name.replace('_', ''))
- self.studies.show_or_save_to_file(plot_name=plot_name, show=self.show, no_title=True, tight_layout=True)
- plt.close()
-
- # Plots "altitude switch" and "peak year"
-
- @property
- def massif_name_to_is_decreasing_parabol(self):
- # For the test we only activate the Mont-Blanc massif
- d = {massif_name: False for massif_name in self.massif_name_to_one_fold_fit.keys()}
- if max(self.study.ordered_years) < 2030:
- for massif_name in ['Vanoise', 'Aravis', 'Beaufortain', 'Chablais']:
- d[massif_name] = True
- return d
-
- @property
- def massif_name_to_altitudes_switch_and_peak_years(self):
- return {massif_name: self.compute_couple_peak_year_and_altitude_switch(massif_name)
- for massif_name, is_decreasing_parabol in self.massif_name_to_is_decreasing_parabol.items()
- if is_decreasing_parabol}
-
- def compute_couple_peak_year_and_altitude_switch(self, massif_name):
- # Get the altitude limits
- altitudes = self.study.massif_name_to_altitudes[massif_name]
- # use a step of 100m for instance
- step = 10
- altitudes = list(np.arange(min(altitudes), max(altitudes) + step, step))
- # Get all the correspond peak years
- margin_function = self.massif_name_to_one_fold_fit[massif_name].best_function_from_fit
- peak_years = []
- year_left = 1900
- switch_altitudes = []
- for altitude in altitudes:
- year_left = self.compute_peak_year(margin_function, altitude, year_left)
- if year_left > 2020:
- break
- peak_years.append(year_left)
- switch_altitudes.append(altitude)
- print(switch_altitudes)
- print(peak_years)
- return switch_altitudes, peak_years
-
- def compute_peak_year(self, margin_function: AbstractMarginFunction, altitude, year_left):
- year_right = year_left + 0.1
- mean_left = margin_function.get_params(np.array([altitude, year_left])).mean
- mean_right = margin_function.get_params(np.array([altitude, year_right])).mean
- print(year_left, year_right, mean_left, mean_right)
- if mean_right < mean_left:
- return year_left
- else:
- return self.compute_peak_year(margin_function, altitude, year_right)
-
- def plot_peak_year_against_altitude(self):
- ax = plt.gca()
- for massif_name, (altitudes, peak_years) in self.massif_name_to_altitudes_switch_and_peak_years.items():
- ax.plot(altitudes, peak_years, label=massif_name)
- ax.legend()
- ax.set_xlabel('Altitude')
- ax.set_ylabel('Peak years')
- plot_name = 'Peak Years'
- self.studies.show_or_save_to_file(plot_name=plot_name, show=self.show)
- plt.close()
-
- def plot_altitude_switch_against_peak_year(self):
- ax = plt.gca()
- for massif_name, (altitudes, peak_years) in self.massif_name_to_altitudes_switch_and_peak_years.items():
- ax.plot(peak_years, altitudes, label=massif_name)
- ax.legend()
- ax.set_xlabel('Peak years')
- ax.set_ylabel('Altitude')
- plot_name = 'Switch altitude'
- self.studies.show_or_save_to_file(plot_name=plot_name, show=self.show)
- plt.close()
-
- def all_trends(self, massif_names):
- """return percents which contain decrease, significant decrease, increase, significant increase percentages"""
- valid_massif_names = self.get_valid_names(massif_names)
-
- nb_valid_massif_names = len(valid_massif_names)
- nbs = np.zeros(4)
- for one_fold in [one_fold for m, one_fold in self.massif_name_to_one_fold_fit.items()
- if m in valid_massif_names]:
- # Compute nb of non stationary models
- if one_fold.change_in_return_level_for_reference_altitude == 0:
- continue
- # Compute nbs
- idx = 0 if one_fold.change_in_return_level_for_reference_altitude < 0 else 2
- nbs[idx] += 1
- if one_fold.is_significant:
- nbs[idx + 1] += 1
-
- percents = 100 * nbs / nb_valid_massif_names
- return [nb_valid_massif_names] + list(percents)
-
- def all_changes(self, massif_names, relative=False):
- """return percents which contain decrease, significant decrease, increase, significant increase percentages"""
- valid_massif_names = self.get_valid_names(massif_names)
- changes = []
- non_stationary_changes = []
- non_stationary_significant_changes = []
- for one_fold in [one_fold for m, one_fold in self.massif_name_to_one_fold_fit.items()
- if m in valid_massif_names]:
- # Compute changes
- if relative:
- change = one_fold.relative_change_in_return_level_for_reference_altitude
- else:
- change = one_fold.change_in_return_level_for_reference_altitude
- changes.append(change)
- if change != 0:
- non_stationary_changes.append(change)
- if one_fold.is_significant:
- non_stationary_significant_changes.append(change)
-
- moment = 'relative mean' if relative else 'Mean'
- print('{} for {}m'.format(moment, self.altitude_group.reference_altitude), np.mean(changes))
- return changes, non_stationary_changes, non_stationary_significant_changes
-
- def get_valid_names(self, massif_names):
- valid_massif_names = set(self.massif_name_to_one_fold_fit.keys())
- if massif_names is not None:
- valid_massif_names = valid_massif_names.intersection(set(massif_names))
- return valid_massif_names
-
- def model_name_to_percentages(self, massif_names, only_significant=False):
- valid_massif_names = self.get_valid_names(massif_names)
- nb_valid_massif_names = len(valid_massif_names)
- best_names = [one_fold_fit.best_estimator.margin_model.name_str
- for m, one_fold_fit in self.massif_name_to_one_fold_fit.items()
- if m in valid_massif_names and (not only_significant or one_fold_fit.is_significant)]
- counter = Counter(best_names)
- d = {name: 100 * c / nb_valid_massif_names for name, c in counter.items()}
- # Add 0 for the name not present
- for name in self.model_names:
- if name not in d:
- d[name] = 0
- return d
-
- @property
- def model_names(self):
- massif_name = list(self.massif_name_to_one_fold_fit.keys())[0]
- return self.massif_name_to_one_fold_fit[massif_name].model_names
-
- def plot_qqplots(self):
- for massif_name, one_fold_fit in self.massif_name_to_one_fold_fit.items():
- ax = plt.gca()
- altitudes = self.massif_name_to_massif_altitudes[massif_name]
- massif_name_corrected = massif_name.replace('_', ' ')
-
- all_quantiles = []
-
- for altitude in self.studies.altitudes:
- coordinate_for_filter = (altitude, None)
- unconstrained_empirical_quantiles = one_fold_fit.best_estimator.sorted_empirical_standard_gumbel_quantiles(coordinate_for_filter=coordinate_for_filter)
- n = len(unconstrained_empirical_quantiles)
- if n > 0:
- assert n == 61
- standard_gumbel_quantiles = one_fold_fit.standard_gumbel_quantiles(n=n)
- ax.plot(standard_gumbel_quantiles, unconstrained_empirical_quantiles, linestyle='None',
- label='{} m'.format(altitude), marker='o')
-
- all_quantiles.extend(standard_gumbel_quantiles)
- all_quantiles.extend(unconstrained_empirical_quantiles)
-
- size_label = 20
- ax.set_xlabel("Theoretical quantile", fontsize=size_label)
- ax.set_ylabel("Empirical quantile", fontsize=size_label)
-
- epsilon = 0.1
- ax_lim = [min(all_quantiles) - epsilon, max(all_quantiles) + epsilon]
- ax.set_xlim(ax_lim)
- ax.set_ylim(ax_lim)
-
- ax.plot(ax_lim, ax_lim, color='k')
- ticks = [i for i in range(ceil(ax_lim[0]), floor(ax_lim[1]) + 1)]
- ax.set_xticks(ticks)
- ax.set_yticks(ticks)
- labelsize = 15
- ax.tick_params(labelsize=labelsize)
- plot_name = 'qqplot/{}'.format(massif_name_corrected)
- handles, labels = ax.get_legend_handles_labels()
- ax.legend(handles[::-1], labels[::-1], prop={'size': labelsize})
- self.studies.show_or_save_to_file(plot_name=plot_name, show=self.show, no_title=True)
- plt.close()
+ def plot_independent(self):
+ for ensemble_fit in self.ensemble_class_to_ensemble_fit.values():
+ ensemble_fit.plot()
diff --git a/projects/projected_snowfall/trends/main_projected_trends_separately.py b/projects/projected_snowfall/trends/main_projected_trends_separately.py
deleted file mode 100644
index 8ed7dcc680526d0b8ea5f7de0b4418ad1225088b..0000000000000000000000000000000000000000
--- a/projects/projected_snowfall/trends/main_projected_trends_separately.py
+++ /dev/null
@@ -1,95 +0,0 @@
-import time
-from typing import List
-
-import matplotlib as mpl
-import numpy as np
-
-from extreme_data.meteo_france_data.adamont_data.abstract_adamont_study import WrongYearMinOrYearMax
-from extreme_data.meteo_france_data.adamont_data.adamont.adamont_snowfall import AdamontSnowfall
-from extreme_data.meteo_france_data.adamont_data.adamont_scenario import AdamontScenario
-from extreme_data.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal, CrocusSnowLoad3Days
-from projects.altitude_spatial_model.altitudes_fit.plots.plot_coherence_curves import plot_coherence_curves
-from projects.altitude_spatial_model.altitudes_fit.plots.plot_histogram_altitude_studies import \
- plot_histogram_all_models_against_altitudes, plot_histogram_all_trends_against_altitudes, \
- plot_shoe_plot_changes_against_altitude
-
-mpl.rcParams['text.usetex'] = True
-mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
-
-from projects.altitude_spatial_model.altitudes_fit.utils_altitude_studies_visualizer import load_visualizer_list
-
-from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitude_group import altitudes_for_groups
-from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.plot_total_aic import plot_individual_aic
-
-from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day, SafranSnowfall3Days, \
- SafranSnowfall5Days, SafranSnowfall7Days, SafranDateFirstSnowfall, SafranPrecipitation1Day, SafranPrecipitation3Days
-from extreme_data.meteo_france_data.scm_models_data.utils import Season
-
-
-def main():
- study_class = AdamontSnowfall
- scenario = AdamontScenario.rcp85_extended
- gcm_rcm_couples = [('CNRM-CM5', 'CCLM4-8-17'), ('CNRM-CM5', 'RCA4')][1:]
-
- fast = True
- if fast is None:
- massif_names = None
- altitudes_list = altitudes_for_groups[1:2]
- elif fast:
- massif_names = ['Vanoise', 'Haute-Maurienne', 'Vercors'][:1]
- altitudes_list = altitudes_for_groups[1:3]
- else:
- massif_names = None
- altitudes_list = altitudes_for_groups[:]
-
- # One plot for each massif name
- for massif_name in massif_names:
- print(massif_name)
- main_loop(altitudes_list, [massif_name], gcm_rcm_couples, study_class, scenario)
-
-
-def main_loop(altitudes_list, massif_names, gcm_rcm_couples, study_class, scenario):
- assert isinstance(altitudes_list, List)
- assert isinstance(altitudes_list[0], List)
- altitudes_for_return_levels = [altitudes[-2] for altitudes in altitudes_list]
- print(altitudes_for_return_levels)
- season = Season.annual
- first_year_min, first_year_max = 1951, 2010
- nb_years = first_year_max - first_year_min + 1
- temporal_windows = [(first_year_min + i, first_year_max + i) for i in [30 * j for j in range(4)]]
- all_changes_in_return_levels = []
- for gcm_rcm_couple in gcm_rcm_couples:
- print('Inner', gcm_rcm_couple)
- changes_in_return_levels = []
- for temporal_window in temporal_windows:
- year_min, year_max = temporal_window
- try:
- visualizer_list = load_visualizer_list(season, study_class, altitudes_list, massif_names,
- scenario=scenario, gcm_rcm_couple=gcm_rcm_couple,
- year_min=year_min, year_max=year_max)
-
- for visualizer in visualizer_list:
- visualizer.studies.plot_maxima_time_series(massif_names)
- massif_name = massif_names[0]
- changes_for_temporal_window = [
- v.massif_name_to_one_fold_fit[massif_name].changes_of_moment(altitudes=[a],
- year=year_max,
- nb_years=nb_years,
- order=None
- )[0]
- for a, v in zip(altitudes_for_return_levels, visualizer_list)]
-
- except WrongYearMinOrYearMax:
- changes_for_temporal_window = [np.nan for _ in altitudes_for_return_levels]
-
- print(changes_for_temporal_window)
- changes_in_return_levels.append(changes_for_temporal_window)
- changes_in_return_levels = np.array(changes_in_return_levels)
- all_changes_in_return_levels.append(changes_in_return_levels)
- all_changes_in_return_levels = np.array(all_changes_in_return_levels)
-
- return all_changes_in_return_levels, temporal_windows, altitudes_for_return_levels, nb_years
-
-
-if __name__ == '__main__':
- main()