From 3b257b005e863b35c5e4fc6dc4d6fb83b985da14 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Tue, 23 Mar 2021 15:36:23 +0100
Subject: [PATCH] [projection swe] modify interval (and data constraint) for
temperature sensitivity. add invidual plots in sensitivity. add changes plot
in sensitivity.
---
.../adamont_data/cmip5/plot_temperatures.py | 2 +-
.../adamont_data/cmip5/temperature_to_year.py | 19 ++--
.../ensemble_fit/abstract_ensemble_fit.py | 2 +-
.../visualizer_for_projection_ensemble.py | 6 +-
.../visualizer_for_sensitivity.py | 89 ++++++++++++++-----
.../discussion/main_sensitivity.py | 15 ++--
.../results/main_projections_ensemble.py | 2 +-
7 files changed, 93 insertions(+), 42 deletions(-)
diff --git a/extreme_data/meteo_france_data/adamont_data/cmip5/plot_temperatures.py b/extreme_data/meteo_france_data/adamont_data/cmip5/plot_temperatures.py
index e8007d1b..2c68e0c8 100644
--- a/extreme_data/meteo_france_data/adamont_data/cmip5/plot_temperatures.py
+++ b/extreme_data/meteo_france_data/adamont_data/cmip5/plot_temperatures.py
@@ -88,6 +88,6 @@ def plot_temperature_for_rcp_gcm(ax, gcm, scenario, year_min, year_max, linestyl
if __name__ == '__main__':
- for anomaly in [True, False][:]:
+ for anomaly in [True, False][:1]:
main_plot_temperature_with_spline_on_top(anomaly=anomaly)
# main_plot_temperature_all(anomaly=True, spline=True)
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 6453d9f0..86c4b975 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
@@ -18,8 +18,8 @@ def get_year_min_and_year_max(gcm, scenario, left_limit, right_limit, is_tempera
else:
year_min, year_max = left_limit, right_limit
- # A minimum of 30 years of data is needed to find a trend
- if year_max - year_min + 1 >= 30:
+ # A minimum of 10 years of data is needed for each GCM/RCM
+ if year_max - year_min + 1 >= 10:
return year_min, year_max
else:
return None, None
@@ -70,10 +70,11 @@ 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:]):
+ first_scenario = i == 0
plot_nb_data_one_line(ax, gcm, scenario, left_limit, right_limit,
- i == 0, is_temperature_interval)
+ first_scenario, is_temperature_interval)
- ax.legend()
+ ax.legend(loc='upper left')
ticks_labels = get_ticks_labels_for_interval(is_temperature_interval, is_shift_interval)
ax.set_xticks(right_limit)
ax.set_xticklabels(ticks_labels)
@@ -84,15 +85,15 @@ def plot_nb_data(is_temperature_interval, is_shift_interval):
Line2D([0], [0], color='k', lw=1, label=scenario_to_str(s),
linestyle=get_linestyle_from_scenario(s)) for s in adamont_scenarios_real
]
- ax2.legend(handles=legend_elements, loc='upper center')
+ ax2.legend(handles=legend_elements, loc='lower right')
ax2.set_yticks([])
plt.show()
def get_interval_limits(is_temperature_interval, is_shift_interval):
if is_temperature_interval:
- temp_min = np.arange(0, 3, 1)
- temp_max = temp_min + 2
+ temp_min = np.arange(0, 2, 0.5)
+ temp_max = temp_min + 1.5
left_limit, right_limit = temp_min, temp_max
else:
shift = 25
@@ -103,7 +104,7 @@ def get_interval_limits(is_temperature_interval, is_shift_interval):
if not is_shift_interval:
min_interval_left = min(left_limit)
left_limit = [min_interval_left for _ in right_limit]
- return left_limit, right_limit
+ return left_limit[:2], right_limit[:2]
def get_ticks_labels_for_interval(is_temperature_interval, is_shift_interval):
@@ -117,7 +118,7 @@ def get_ticks_labels_for_interval(is_temperature_interval, is_shift_interval):
if __name__ == '__main__':
- for shift_interval in [False, True]:
+ for shift_interval in [False, True][:1]:
for temp_interval in [False, True][1:]:
print("shift = {}, temp_inteval = {}".format(shift_interval, temp_interval))
plot_nb_data(is_temperature_interval=temp_interval, is_shift_interval=shift_interval)
diff --git a/extreme_trend/ensemble_fit/abstract_ensemble_fit.py b/extreme_trend/ensemble_fit/abstract_ensemble_fit.py
index 115b0700..b9348041 100644
--- a/extreme_trend/ensemble_fit/abstract_ensemble_fit.py
+++ b/extreme_trend/ensemble_fit/abstract_ensemble_fit.py
@@ -29,7 +29,7 @@ class AbstractEnsembleFit(object):
# Set appropriate setting
OneFoldFit.last_year = 2100
- OneFoldFit.nb_years = 95
+ OneFoldFit.nb_years = OneFoldFit.last_year - 2005
@property
def altitudes(self):
diff --git a/extreme_trend/ensemble_fit/visualizer_for_projection_ensemble.py b/extreme_trend/ensemble_fit/visualizer_for_projection_ensemble.py
index a75e81b0..70ba3371 100644
--- a/extreme_trend/ensemble_fit/visualizer_for_projection_ensemble.py
+++ b/extreme_trend/ensemble_fit/visualizer_for_projection_ensemble.py
@@ -31,7 +31,9 @@ class VisualizerForProjectionEnsemble(object):
confidence_interval_based_on_delta_method=False,
remove_physically_implausible_models=False,
gcm_to_year_min_and_year_max=None,
+ interval_str_prefix='',
):
+ self.interval_str_prefix = interval_str_prefix
self.altitudes_list = altitudes_list
self.temporal_covariate_for_fit = temporal_covariate_for_fit
self.scenario = scenario
@@ -119,14 +121,14 @@ class VisualizerForProjectionEnsemble(object):
]
if key in merge_keys:
for v in visualizer_list:
- v.studies.study.gcm_rcm_couple = (key, "merge")
+ v.studies.study.gcm_rcm_couple = ("{} {}".format(key, "merge"), self.interval_str_prefix)
self.plot_for_visualizer_list(visualizer_list)
def plot_together(self):
visualizer_list = [together_ensemble_fit.visualizer
for together_ensemble_fit in self.ensemble_fits(TogetherEnsembleFit)]
for v in visualizer_list:
- v.studies.study.gcm_rcm_couple = ("together", "merge")
+ v.studies.study.gcm_rcm_couple = ("together merge", self.interval_str_prefix)
self.plot_for_visualizer_list(visualizer_list)
def ensemble_fits(self, ensemble_class):
diff --git a/extreme_trend/ensemble_fit/visualizer_for_sensitivity.py b/extreme_trend/ensemble_fit/visualizer_for_sensitivity.py
index bd696188..453e8ece 100644
--- a/extreme_trend/ensemble_fit/visualizer_for_sensitivity.py
+++ b/extreme_trend/ensemble_fit/visualizer_for_sensitivity.py
@@ -2,6 +2,8 @@ from collections import OrderedDict
import matplotlib.pyplot as plt
from typing import List, Dict
+import numpy as np
+
from extreme_data.meteo_france_data.adamont_data.cmip5.temperature_to_year import get_interval_limits, \
get_year_min_and_year_max, get_ticks_labels_for_interval
from extreme_data.meteo_france_data.scm_models_data.utils import Season
@@ -40,11 +42,14 @@ class VisualizerForSensivity(object):
self.massif_names = massif_names
self.left_limits, self.right_limits = get_interval_limits(self.is_temperature_interval,
self.is_shift_interval)
+
self.left_limit_to_right_limit = OrderedDict(zip(self.left_limits, self.right_limits))
self.right_limit_to_visualizer = {} # type: Dict[float, VisualizerForProjectionEnsemble]
for left_limit, right_limit in zip(self.left_limits, self.right_limits):
- print("Interval is", left_limit, right_limit)
+ interval_str_prefix = "{}-{}".format(left_limit, right_limit)
+ interval_str_prefix = interval_str_prefix.replace('.', ',')
+ print("Interval is", interval_str_prefix)
# Build gcm_to_year_min_and_year_max
gcm_to_year_min_and_year_max = {}
gcm_list = list(set([g for g, r in gcm_rcm_couples]))
@@ -64,29 +69,61 @@ class VisualizerForSensivity(object):
massif_names=massif_names,
temporal_covariate_for_fit=temporal_covariate_for_fit,
remove_physically_implausible_models=remove_physically_implausible_models,
- gcm_to_year_min_and_year_max=gcm_to_year_min_and_year_max
+ gcm_to_year_min_and_year_max=gcm_to_year_min_and_year_max,
+ interval_str_prefix=interval_str_prefix,
)
self.right_limit_to_visualizer[right_limit] = visualizer
def plot(self):
- # todo: before reactivating the subplot, i should ensure that we can modify the prefix
- # todo: so that we can have all the subplot for the merge visualizer
- # , and not just the plots for the last t_min
- # for visualizer in self.temp_min_to_visualizer.values():
- # visualizer.plot()
+ for visualizer in self.right_limit_to_visualizer.values():
+ visualizer.plot()
merge_visualizer_str_list = []
if IndependentEnsembleFit in self.ensemble_fit_classes:
merge_visualizer_str_list.extend([AbstractEnsembleFit.Median_merge, AbstractEnsembleFit.Mean_merge])
if TogetherEnsembleFit in self.ensemble_fit_classes:
merge_visualizer_str_list.append(AbstractEnsembleFit.Together_merge)
for merge_visualizer_str in merge_visualizer_str_list:
- self.sensitivity_plot(merge_visualizer_str)
+ self.sensitivity_plot_percentages(merge_visualizer_str)
+ for relative in [True, False]:
+ self.sensitivity_plot_changes(merge_visualizer_str, relative)
- def sensitivity_plot(self, merge_visualizer_str):
+ def sensitivity_plot_changes(self, merge_visualizer_str, relative):
ax = plt.gca()
for altitudes in self.altitudes_list:
altitude_class = get_altitude_class_from_altitudes(altitudes)
- self.interval_plot(ax, altitude_class, merge_visualizer_str)
+ self.interval_plot_changes(ax, altitude_class, merge_visualizer_str, relative)
+
+ ticks_labels = get_ticks_labels_for_interval(self.is_temperature_interval, self.is_shift_interval)
+ name = 'relative changes' if relative else "changes"
+ 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_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_percentages(self, merge_visualizer_str):
+ ax = plt.gca()
+ for altitudes in self.altitudes_list:
+ altitude_class = get_altitude_class_from_altitudes(altitudes)
+ self.interval_plot_percentages(ax, altitude_class, merge_visualizer_str)
ticks_labels = get_ticks_labels_for_interval(self.is_temperature_interval, self.is_shift_interval)
ax.set_ylabel('Percentages of massifs (\%)')
@@ -96,16 +133,9 @@ class VisualizerForSensivity(object):
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)])
- merge_visualizer = self.first_merge_visualizer(merge_visualizer_str)
- temp_cov = self.temporal_covariate_for_fit is AnomalyTemperatureWithSplineTemporalCovariate
- merge_visualizer.plot_name = 'Sensitivity plot with ' \
- 'shift={} temp_interval={}, temp_cov={}'.format(self.is_shift_interval,
- self.is_temperature_interval,
- temp_cov)
- merge_visualizer.show_or_save_to_file(no_title=True)
- plt.close()
+ self.save_plot(merge_visualizer_str, "percentages")
- def interval_plot(self, ax, altitude_class, merge_visualizer_str):
+ def interval_plot_percentages(self, ax, altitude_class, merge_visualizer_str):
linestyle = get_linestyle_for_altitude_class(altitude_class)
increasing_key = 'increasing'
decreasing_key = 'decreasing'
@@ -130,16 +160,19 @@ class VisualizerForSensivity(object):
color = label_to_color[label]
ax.plot(self.right_limits, l, label=label_improved, color=color, linestyle=linestyle)
+ # Merge visualizer
+
def get_merge_visualizer(self, altitude_class, visualizer_projection: VisualizerForProjectionEnsemble,
merge_visualizer_str):
if merge_visualizer_str in [AbstractEnsembleFit.Median_merge, AbstractEnsembleFit.Mean_merge]:
independent_ensemble_fit = \
- visualizer_projection.altitude_class_to_ensemble_class_to_ensemble_fit[altitude_class][
- IndependentEnsembleFit]
+ visualizer_projection.altitude_class_to_ensemble_class_to_ensemble_fit[altitude_class][
+ IndependentEnsembleFit]
merge_visualizer = independent_ensemble_fit.merge_function_name_to_visualizer[merge_visualizer_str]
else:
together_ensemble_fit = \
- visualizer_projection.altitude_class_to_ensemble_class_to_ensemble_fit[altitude_class][TogetherEnsembleFit]
+ visualizer_projection.altitude_class_to_ensemble_class_to_ensemble_fit[altitude_class][
+ TogetherEnsembleFit]
merge_visualizer = together_ensemble_fit.visualizer
merge_visualizer.studies.study.gcm_rcm_couple = (merge_visualizer_str, "merge")
return merge_visualizer
@@ -148,3 +181,15 @@ class VisualizerForSensivity(object):
altitude_class = get_altitude_class_from_altitudes(self.altitudes_list[0])
visualizer_projection = list(self.right_limit_to_visualizer.values())[0]
return self.get_merge_visualizer(altitude_class, visualizer_projection, merge_visualizer_str)
+
+ # Utils
+
+ def save_plot(self, merge_visualizer_str, name):
+ merge_visualizer = self.first_merge_visualizer(merge_visualizer_str)
+ temp_cov = self.temporal_covariate_for_fit is AnomalyTemperatureWithSplineTemporalCovariate
+ merge_visualizer.plot_name = 'Sensitivity plot for {} with ' \
+ 'shift={} temp_interval={}, temp_cov={}'.format(name, self.is_shift_interval,
+ self.is_temperature_interval,
+ temp_cov)
+ merge_visualizer.show_or_save_to_file(no_title=True)
+ plt.close()
diff --git a/projects/projected_extreme_snowfall/discussion/main_sensitivity.py b/projects/projected_extreme_snowfall/discussion/main_sensitivity.py
index a2e5e60c..70c7de45 100644
--- a/projects/projected_extreme_snowfall/discussion/main_sensitivity.py
+++ b/projects/projected_extreme_snowfall/discussion/main_sensitivity.py
@@ -44,21 +44,24 @@ def main():
set_seed_for_test()
AbstractExtractEurocodeReturnLevel.ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY = 0.2
- fast = False
- scenarios = rcp_scenarios[-1:] if fast is False else [AdamontScenario.rcp85]
+ fast = None
+ scenarios = [AdamontScenario.rcp85]
+ scenarios = rcp_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[3:]
+ 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]
+ altitudes_list = altitudes_for_groups[1:3]
else:
massif_names = None
altitudes_list = altitudes_for_groups[:]
@@ -72,11 +75,11 @@ def main():
model_classes = ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS
assert scenario in rcp_scenarios
remove_physically_implausible_models = True
- temp_cov = False
+ temp_cov = True
temporal_covariate_for_fit = AnomalyTemperatureWithSplineTemporalCovariate if temp_cov else TimeTemporalCovariate
print('Covariate is {}'.format(temporal_covariate_for_fit))
- for is_temperature_interval in [True, False][1:]:
+ for is_temperature_interval in [True, False][:1]:
for is_shift_interval in [True, False][1:]:
visualizer = VisualizerForSensivity(
altitudes_list, gcm_rcm_couples, study_class, Season.annual, scenario,
diff --git a/projects/projected_extreme_snowfall/results/main_projections_ensemble.py b/projects/projected_extreme_snowfall/results/main_projections_ensemble.py
index 01b20911..edc30d24 100644
--- a/projects/projected_extreme_snowfall/results/main_projections_ensemble.py
+++ b/projects/projected_extreme_snowfall/results/main_projections_ensemble.py
@@ -45,7 +45,7 @@ def main():
set_seed_for_test()
AbstractExtractEurocodeReturnLevel.ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY = 0.2
- fast = False
+ fast = True
scenarios = rcp_scenarios[::-1] if fast is False else [AdamontScenario.rcp85]
for scenario in scenarios:
--
GitLab