From 1188c60803c34dc04e9b76083cb97949fbcbce1f Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Tue, 23 Mar 2021 17:26:46 +0100
Subject: [PATCH] [projection swe] add sensitivity for return levels.
---
.../adamont_data/cmip5/temperature_to_year.py | 4 +-
.../visualizer_for_sensitivity.py | 46 +++++++++++++++++++
...es_visualizer_for_non_stationary_models.py | 18 ++++++--
extreme_trend/one_fold_fit/one_fold_fit.py | 11 +++++
.../discussion/main_sensitivity.py | 14 +++---
5 files changed, 80 insertions(+), 13 deletions(-)
diff --git a/extreme_data/meteo_france_data/adamont_data/cmip5/temperature_to_year.py b/extreme_data/meteo_france_data/adamont_data/cmip5/temperature_to_year.py
index 86c4b975..0a476c32 100644
--- a/extreme_data/meteo_france_data/adamont_data/cmip5/temperature_to_year.py
+++ b/extreme_data/meteo_france_data/adamont_data/cmip5/temperature_to_year.py
@@ -69,7 +69,7 @@ def plot_nb_data(is_temperature_interval, is_shift_interval):
ax = plt.gca()
for gcm in get_gcm_list(adamont_version=2)[:]:
- for i, scenario in enumerate(rcp_scenarios[2:]):
+ for i, scenario in enumerate(rcp_scenarios[1:]):
first_scenario = i == 0
plot_nb_data_one_line(ax, gcm, scenario, left_limit, right_limit,
first_scenario, is_temperature_interval)
@@ -93,7 +93,7 @@ def plot_nb_data(is_temperature_interval, is_shift_interval):
def get_interval_limits(is_temperature_interval, is_shift_interval):
if is_temperature_interval:
temp_min = np.arange(0, 2, 0.5)
- temp_max = temp_min + 1.5
+ temp_max = temp_min + 2
left_limit, right_limit = temp_min, temp_max
else:
shift = 25
diff --git a/extreme_trend/ensemble_fit/visualizer_for_sensitivity.py b/extreme_trend/ensemble_fit/visualizer_for_sensitivity.py
index 453e8ece..70f12b95 100644
--- a/extreme_trend/ensemble_fit/visualizer_for_sensitivity.py
+++ b/extreme_trend/ensemble_fit/visualizer_for_sensitivity.py
@@ -59,6 +59,7 @@ class VisualizerForSensivity(object):
if year_min_and_year_max[0] is not None:
gcm_to_year_min_and_year_max[gcm] = year_min_and_year_max
+ print(gcm_to_year_min_and_year_max)
visualizer = VisualizerForProjectionEnsemble(
altitudes_list, gcm_rcm_couples, study_class, season, scenario,
model_classes=model_classes,
@@ -84,9 +85,54 @@ class VisualizerForSensivity(object):
merge_visualizer_str_list.append(AbstractEnsembleFit.Together_merge)
for merge_visualizer_str in merge_visualizer_str_list:
self.sensitivity_plot_percentages(merge_visualizer_str)
+ self.sensitivity_plot_return_levels(merge_visualizer_str)
for relative in [True, False]:
self.sensitivity_plot_changes(merge_visualizer_str, relative)
+ def sensitivity_plot_return_levels(self, merge_visualizer_str):
+ ax = plt.gca()
+ for altitudes in self.altitudes_list:
+ altitude_class = get_altitude_class_from_altitudes(altitudes)
+ self.interval_plot_return_levels(ax, altitude_class, merge_visualizer_str)
+
+ ticks_labels = get_ticks_labels_for_interval(self.is_temperature_interval, self.is_shift_interval)
+ name = 'Return levels at the end of the interval'
+ ax.set_ylabel(name)
+ ax.set_xlabel('Interval used to compute the trends ')
+ ax.set_xticks(self.right_limits)
+ ax.set_xticklabels(ticks_labels)
+ lim_left, lim_right = ax.get_xlim()
+ ax.hlines(0, xmin=lim_left, xmax=lim_right, linestyles='dashed')
+ ax.legend(prop={'size': 7}, loc='upper center', ncol=2)
+ # ax.set_ylim((0, 122))
+ # ax.set_yticks([i * 10 for i in range(11)])
+ self.save_plot(merge_visualizer_str, name)
+
+ def interval_plot_return_levels(self, ax, altitude_class, merge_visualizer_str):
+ linestyle = get_linestyle_for_altitude_class(altitude_class)
+
+ mean_return_levels = []
+ for v in self.right_limit_to_visualizer.values():
+ merge_visualizer = self.get_merge_visualizer(altitude_class, v, merge_visualizer_str)
+ mean_return_level = merge_visualizer.mean_return_level(self.massif_names)
+ mean_return_levels.append(mean_return_level)
+
+ label = altitude_class().formula
+ ax.plot(self.right_limits, mean_return_levels, linestyle=linestyle, label=label, color='orange')
+
+ def interval_plot_changes(self, ax, altitude_class, merge_visualizer_str, relative):
+ linestyle = get_linestyle_for_altitude_class(altitude_class)
+
+ mean_changes = []
+ for v in self.right_limit_to_visualizer.values():
+ merge_visualizer = self.get_merge_visualizer(altitude_class, v, merge_visualizer_str)
+ changes, non_stationary_changes = merge_visualizer.all_changes(self.massif_names, relative=relative,
+ with_significance=False)
+ mean_changes.append(np.mean(changes))
+
+ label = altitude_class().formula
+ ax.plot(self.right_limits, mean_changes, linestyle=linestyle, label=label, color='darkgreen')
+
def sensitivity_plot_changes(self, merge_visualizer_str, relative):
ax = plt.gca()
for altitudes in self.altitudes_list:
diff --git a/extreme_trend/one_fold_fit/altitudes_studies_visualizer_for_non_stationary_models.py b/extreme_trend/one_fold_fit/altitudes_studies_visualizer_for_non_stationary_models.py
index c555a33c..73b674f5 100644
--- a/extreme_trend/one_fold_fit/altitudes_studies_visualizer_for_non_stationary_models.py
+++ b/extreme_trend/one_fold_fit/altitudes_studies_visualizer_for_non_stationary_models.py
@@ -266,7 +266,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
moment = moment.replace('moment',
'{} in {}'.format(OneFoldFit.get_moment_str(order=order), OneFoldFit.last_year))
plot_name = 'Model {} annual maxima of {}'.format(moment,
- SCM_STUDY_CLASS_TO_ABBREVIATION[self.studies.study_class])
+ 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)
@@ -323,9 +323,9 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
raise NotImplementedError
self.plot_map(massif_name_to_value=massif_name_to_best_coef,
label='Coef/{} plot for {} {}'.format(coef_name,
- SCM_STUDY_CLASS_TO_ABBREVIATION[
- type(self.study)],
- self.study.variable_unit),
+ 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)
@@ -388,7 +388,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
massif_name.replace('_', ' '), self.study.variable_unit))
peak_year_folder = 'Peak year ' + ylabel
plot_name = '{}/Peak year for {}'.format(peak_year_folder,
- massif_name.replace('_', ''))
+ 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()
@@ -462,6 +462,14 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
self.studies.show_or_save_to_file(plot_name=plot_name, show=self.show)
plt.close()
+ def mean_return_level(self, massif_names):
+ valid_massif_names = self.get_valid_names(massif_names)
+ return_levels = []
+ 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]:
+ return_levels.append(one_fold.return_level_last_temporal_coordinate)
+ return np.mean(return_levels)
+
def all_trends(self, massif_names, with_significance=True, with_relative_change=False):
"""return percents which contain decrease, significant decrease, increase, significant increase percentages"""
valid_massif_names = self.get_valid_names(massif_names)
diff --git a/extreme_trend/one_fold_fit/one_fold_fit.py b/extreme_trend/one_fold_fit/one_fold_fit.py
index b62a83d3..5146cb04 100644
--- a/extreme_trend/one_fold_fit/one_fold_fit.py
+++ b/extreme_trend/one_fold_fit/one_fold_fit.py
@@ -29,6 +29,7 @@ from extreme_fit.model.result_from_model_fit.result_from_extremes.eurocode_retur
EurocodeConfidenceIntervalFromExtremes
from extreme_trend.one_fold_fit.altitude_group import DefaultAltitudeGroup, altitudes_for_groups
from root_utils import NB_CORES, batch
+from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_covariate_for_fit import TimeTemporalCovariate
from spatio_temporal_dataset.coordinates.temporal_coordinates.temperature_covariate import \
AnomalyTemperatureWithSplineTemporalCovariate
@@ -382,6 +383,16 @@ class OneFoldFit(object):
return is_significant, altitude_and_year_to_return_level_mean_estimate, altitude_and_year_to_return_level_confidence_interval
+ @property
+ def return_level_last_temporal_coordinate(self):
+ df_temporal_covariate = self.dataset.coordinates.df_temporal_coordinates_for_fit(temporal_covariate_for_fit=self.temporal_covariate_for_fit,
+ drop_duplicates=False)
+ last_temporal_coordinate = df_temporal_covariate.loc[:, AbstractCoordinates.COORDINATE_T].max()
+ print('last temporal coordinate', last_temporal_coordinate)
+ altitude = self.altitude_group.reference_altitude
+ coordinate = np.array([altitude, last_temporal_coordinate])
+ return self.get_return_level(self.best_function_from_fit, coordinate)
+
@property
def bootstrap_fitted_functions_from_fit(self):
print('nb of bootstrap for confidence interval=', AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP)
diff --git a/projects/projected_extreme_snowfall/discussion/main_sensitivity.py b/projects/projected_extreme_snowfall/discussion/main_sensitivity.py
index 70c7de45..b125e224 100644
--- a/projects/projected_extreme_snowfall/discussion/main_sensitivity.py
+++ b/projects/projected_extreme_snowfall/discussion/main_sensitivity.py
@@ -25,7 +25,7 @@ from extreme_fit.model.utils import set_seed_for_test
from extreme_data.meteo_france_data.adamont_data.adamont.adamont_safran import AdamontSnowfall
from extreme_data.meteo_france_data.adamont_data.adamont_scenario import AdamontScenario, get_gcm_rcm_couples, \
- rcp_scenarios
+ rcp_scenarios, rcm_scenarios_extended
from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_covariate_for_fit import \
TimeTemporalCovariate
@@ -44,24 +44,26 @@ def main():
set_seed_for_test()
AbstractExtractEurocodeReturnLevel.ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY = 0.2
- fast = None
+ fast = True
scenarios = [AdamontScenario.rcp85]
+ scenarios = rcm_scenarios_extended[1:]
scenarios = rcp_scenarios[1:]
+ if fast in [None, True]:
+ scenarios = scenarios[:1]
+
for scenario in scenarios:
gcm_rcm_couples = get_gcm_rcm_couples(scenario)
if fast is None:
- scenarios = scenarios[:1]
massif_names = None
gcm_rcm_couples = gcm_rcm_couples[4:6]
AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP = 10
altitudes_list = altitudes_for_groups[1:3]
elif fast:
- scenarios = scenarios[:1]
massif_names = ['Vanoise', 'Haute-Maurienne']
gcm_rcm_couples = gcm_rcm_couples[4:6]
AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP = 10
- altitudes_list = altitudes_for_groups[1:3]
+ altitudes_list = altitudes_for_groups[2:3]
else:
massif_names = None
altitudes_list = altitudes_for_groups[:]
@@ -73,7 +75,7 @@ def main():
print('Scenario is', scenario)
model_classes = ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS
- assert scenario in rcp_scenarios
+ assert scenario is not AdamontScenario.histo
remove_physically_implausible_models = True
temp_cov = True
temporal_covariate_for_fit = AnomalyTemperatureWithSplineTemporalCovariate if temp_cov else TimeTemporalCovariate
--
GitLab