From 4faa115a861d5500f4ce9c827cefcb0d403fc0e3 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Tue, 10 Dec 2019 14:57:56 +0100
Subject: [PATCH] [PAPER 1] create relative change plot
---
.../plot/create_shifted_cmap.py | 2 +-
.../main_result_trends_and_return_levels.py | 27 ++++---
.../plot_uncertainty_curves.py | 4 +-
...dy_visualizer_for_non_stationary_trends.py | 77 +++++++++++++------
.../abstract_gev_trend_test.py | 26 +++++--
5 files changed, 97 insertions(+), 39 deletions(-)
diff --git a/experiment/meteo_france_data/plot/create_shifted_cmap.py b/experiment/meteo_france_data/plot/create_shifted_cmap.py
index 73c4a88e..dd49798c 100644
--- a/experiment/meteo_france_data/plot/create_shifted_cmap.py
+++ b/experiment/meteo_france_data/plot/create_shifted_cmap.py
@@ -32,7 +32,7 @@ def create_colorbase_axis(ax, label, cmap, norm, ticks_values_and_labels=None):
if ticks_values_and_labels is not None:
cb.ax.set_yticklabels([str(t) for t in ticks_values_and_labels[1]])
if isinstance(label, str):
- cb.set_label(label)
+ cb.set_label(label, fontsize=15)
return norm
diff --git a/experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py b/experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py
index 796edeae..65332ae3 100644
--- a/experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py
+++ b/experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py
@@ -17,14 +17,14 @@ from root_utils import NB_CORES
mpl.rcParams['text.usetex'] = True
mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
-
+import matplotlib.pyplot as plt
def minor_result(altitude):
"""Plot trends for a single altitude to be fast"""
visualizer = StudyVisualizerForNonStationaryTrends(CrocusSnowLoadTotal(altitude=altitude), multiprocessing=True,
)
visualizer.plot_trends()
- # plt.show()
+ plt.show()
def compute_minimized_aic(visualizer):
_ = visualizer.massif_name_to_minimized_aic_non_stationary_trend_test
@@ -50,7 +50,7 @@ def intermediate_result(altitudes, massif_names=None,
# Load variable object efficiently
for v in altitude_to_visualizer.values():
_ = v.study.year_to_variable_object
- # Plot trends
+ # Compute minimized value efficiently
visualizers = list(altitude_to_visualizer.values())
if multiprocessing:
with Pool(NB_CORES) as p:
@@ -59,9 +59,13 @@ def intermediate_result(altitudes, massif_names=None,
for visualizer in visualizers:
_ = compute_minimized_aic(visualizer)
# Compute common max value for the colorbar
- max_abs_tdrl = max([visualizer.max_abs_tdrl for visualizer in visualizers])
- for visualizer in altitude_to_visualizer.values():
- visualizer.plot_trends(max_abs_tdrl)
+ altitudes_for_plot_trend = [900, 1800, 2700]
+ visualizers_for_altitudes = [visualizer
+ for altitude, visualizer in altitude_to_visualizer.items()
+ if altitude in altitudes_for_plot_trend]
+ max_abs_tdrl = max([visualizer.max_abs_change for visualizer in visualizers_for_altitudes])
+ for visualizer in visualizers_for_altitudes:
+ visualizer.plot_trends(max_abs_tdrl, add_colorbar=visualizer.study.altitude==2700)
# Plot graph
plot_uncertainty_massifs(altitude_to_visualizer)
@@ -73,7 +77,7 @@ def major_result():
uncertainty_methods = [ConfidenceIntervalMethodFromExtremes.my_bayes,
ConfidenceIntervalMethodFromExtremes.ci_mle][:]
massif_names = None
- for study_class in paper_study_classes[1:2]:
+ for study_class in paper_study_classes[:1]:
if study_class == CrocusSnowLoadEurocode:
non_stationary_uncertainty = [False]
else:
@@ -82,16 +86,21 @@ def major_result():
if __name__ == '__main__':
- major_result()
+ # major_result()
# intermediate_result(altitudes=paper_altitudes, massif_names=['Maurienne'],
# uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
# ConfidenceIntervalMethodFromExtremes.ci_mle][1:],
# non_stationary_uncertainty=[False, True][1:],
# multiprocessing=True)
+ intermediate_result(altitudes=paper_altitudes, massif_names=['Maurienne'],
+ uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
+ ConfidenceIntervalMethodFromExtremes.ci_mle][:],
+ non_stationary_uncertainty=[False, True][:],
+ multiprocessing=True)
# intermediate_result(altitudes=[900, 1200], massif_names=None)
# intermediate_result(ALL_ALTITUDES_WITHOUT_NAN)
# intermediate_result(paper_altitudes)
- # minor_result(altitude=600)
+ # minor_result(altitude=900)
# intermediate_result(altitudes=[1500, 1800], massif_names=['Chartreuse'],
# uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.ci_mle,
# ConfidenceIntervalMethodFromExtremes.ci_bayes],
diff --git a/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py b/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py
index 5ad18222..d1da622f 100644
--- a/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py
+++ b/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py
@@ -138,7 +138,9 @@ def plot_tdrl_bars(altitude_to_visualizer, ax, massif_name, valid_altitudes, leg
if len(visualizers) > 0:
tdrl_values = [v.massif_name_to_tdrl_value[massif_name] for v in visualizers]
# Plot bars
- colors = [v.massif_name_to_tdrl_color[massif_name] for v in visualizers]
+ # colors = [v.massif_name_to_color[massif_name] for v in visualizers]
+ # From snow on, we set a black color for the bars
+ colors = ['white' for v in visualizers]
ax.bar(valid_altitudes, tdrl_values, width=150, color=colors, label=visualizers[0].label_tdrl_bar,
edgecolor='black', hatch='//')
# Plot markers
diff --git a/experiment/paper_past_snow_loads/study_visualizer_for_non_stationary_trends.py b/experiment/paper_past_snow_loads/study_visualizer_for_non_stationary_trends.py
index da62927e..4fa5978d 100644
--- a/experiment/paper_past_snow_loads/study_visualizer_for_non_stationary_trends.py
+++ b/experiment/paper_past_snow_loads/study_visualizer_for_non_stationary_trends.py
@@ -38,12 +38,14 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
uncertainty_methods=None,
non_stationary_contexts=None,
uncertainty_massif_names=None,
- effective_temporal_covariate=2017):
+ effective_temporal_covariate=2017,
+ relative_change_trend_plot=True):
super().__init__(study, show, save_to_file, only_one_graph, only_first_row, vertical_kde_plot,
year_for_kde_plot, plot_block_maxima_quantiles, temporal_non_stationarity,
transformation_class, verbose, multiprocessing, complete_non_stationary_trend_analysis,
normalization_under_one_observations, score_class)
# Add some attributes
+ self.relative_change_trend_plot = relative_change_trend_plot
self.effective_temporal_covariate = effective_temporal_covariate
self.non_stationary_contexts = non_stationary_contexts
self.uncertainty_methods = uncertainty_methods
@@ -60,7 +62,7 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
self.non_stationary_trend_test = [GevLocationTrendTest, GevScaleTrendTest, GevLocationAndScaleTrendTest]
self.non_stationary_trend_test_to_marker = dict(zip(self.non_stationary_trend_test, ["s", "^", "D"]))
- self.global_max_abs_tdrl = None
+ self.global_max_abs_change = None
# Utils
@@ -105,53 +107,64 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
massif_name_to_trend_test_that_minimized_aic[massif_name] = trend_test_that_minimized_aic
return massif_name_to_trend_test_that_minimized_aic
-
# Part 1 - Trends
@property
- def max_abs_tdrl(self):
- return max(abs(min(self.massif_name_to_tdrl_value.values())), max(self.massif_name_to_tdrl_value.values()))
+ def max_abs_change(self):
+ return max(abs(min(self.massif_name_to_change_value.values())), max(self.massif_name_to_change_value.values()))
@cached_property
- def _max_abs_tdrl(self):
- return self.global_max_abs_tdrl if self.global_max_abs_tdrl is not None else self.max_abs_tdrl
+ def _max_abs_change(self):
+ return self.global_max_abs_change if self.global_max_abs_change is not None else self.max_abs_change
- def plot_trends(self, max_abs_tdrl=None):
+ def plot_trends(self, max_abs_tdrl=None, add_colorbar=True):
if max_abs_tdrl is not None:
- self.global_max_abs_tdrl = max_abs_tdrl
- ax = self.study.visualize_study(massif_name_to_value=self.massif_name_to_tdrl_value,
+ self.global_max_abs_change = max_abs_tdrl
+ ax = self.study.visualize_study(massif_name_to_value=self.massif_name_to_change_value,
replace_blue_by_white=False,
axis_off=False, show_label=False,
- add_colorbar=True,
+ add_colorbar=add_colorbar,
massif_name_to_marker_style=self.massif_name_to_marker_style,
- massif_name_to_color=self.massif_name_to_tdrl_color,
+ massif_name_to_color=self.massif_name_to_color,
cmap=self.cmap,
show=False,
ticks_values_and_labels=self.ticks_values_and_labels,
- label=self.label_tdrl_bar + ' for {}'.format(EUROCODE_RETURN_LEVEL_STR))
+ label=self.label)
ax.get_xaxis().set_visible(True)
ax.set_xticks([])
- ax.set_xlabel('Altitude = {}m'.format(self.study.altitude), fontsize=12)
+ ax.set_xlabel('Altitude = {}m'.format(self.study.altitude), fontsize=15)
self.plot_name = 'tdlr_trends'
self.show_or_save_to_file(add_classic_title=False, tight_layout=True, no_title=True,
- dpi=dpi_paper1_figure)
+ dpi=500)
plt.close()
+ @property
+ def label(self):
+ if self.relative_change_trend_plot:
+ label_tdlr_bar = 'Relative change between {} and {}'.format(self.initial_year, self.final_year)
+ else:
+ label_tdlr_bar = self.label_tdrl_bar
+ label = label_tdlr_bar + '\nfor {}'.format(EUROCODE_RETURN_LEVEL_STR)
+ if self.relative_change_trend_plot:
+ # change units
+ label = label.split('(')[0] + '(\%)'
+ return label
+
@property
def label_tdrl_bar(self):
return 'Change in {} years'.format(AbstractGevTrendTest.nb_years_for_quantile_evolution)
@property
def ticks_values_and_labels(self):
- graduation = 0.2
+ graduation = 10 if self.relative_change_trend_plot else 0.2
positive_ticks = []
tick = graduation
- while tick < self._max_abs_tdrl:
+ while tick < self._max_abs_change:
positive_ticks.append(round(tick, 1))
tick += graduation
all_ticks_labels = [-t for t in positive_ticks] + [0] + positive_ticks
- ticks_values = [((t / self._max_abs_tdrl) + 1) / 2 for t in all_ticks_labels]
+ ticks_values = [((t / self._max_abs_change) + 1) / 2 for t in all_ticks_labels]
return ticks_values, all_ticks_labels
@cached_property
@@ -159,14 +172,34 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
return {m: t.time_derivative_of_return_level for m, t in
self.massif_name_to_minimized_aic_non_stationary_trend_test.items()}
+ @cached_property
+ def massif_name_to_relative_change_value(self):
+ return {m: t.relative_change_in_return_level(initial_year=self.initial_year, final_year=self.final_year)
+ for m, t in self.massif_name_to_minimized_aic_non_stationary_trend_test.items()}
+
+ @property
+ def initial_year(self):
+ return self.final_year - 50
+
+ @property
+ def final_year(self):
+ return 2010
+
+ @cached_property
+ def massif_name_to_change_value(self):
+ if self.relative_change_trend_plot:
+ return self.massif_name_to_relative_change_value
+ else:
+ return self.massif_name_to_tdrl_value
+
@cached_property
def cmap(self):
- return get_shifted_map(-self._max_abs_tdrl, self._max_abs_tdrl)
+ return get_shifted_map(-self._max_abs_change, self._max_abs_change)
@cached_property
- def massif_name_to_tdrl_color(self):
- return {m: get_colors([v], self.cmap, -self._max_abs_tdrl, self._max_abs_tdrl)[0]
- for m, v in self.massif_name_to_tdrl_value.items()}
+ def massif_name_to_color(self):
+ return {m: get_colors([v], self.cmap, -self._max_abs_change, self._max_abs_change)[0]
+ for m, v in self.massif_name_to_change_value.items()}
@cached_property
def massif_name_to_marker_style(self):
diff --git a/experiment/trend_analysis/univariate_test/abstract_gev_trend_test.py b/experiment/trend_analysis/univariate_test/abstract_gev_trend_test.py
index 59184411..71ff1725 100644
--- a/experiment/trend_analysis/univariate_test/abstract_gev_trend_test.py
+++ b/experiment/trend_analysis/univariate_test/abstract_gev_trend_test.py
@@ -143,17 +143,31 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
return self.constrained_estimator.margin_function_from_fit.get_gev_params(coordinate=np.array([1958]),
is_transformed=False)
+ def time_derivative_times_years(self, nb_years):
+ # Compute the slope strength
+ slope = self._slope_strength()
+ # Delta T must in the same unit as were the parameter of slope mu1 and sigma1
+ slope *= nb_years * self.coordinates.transformed_distance_between_two_successive_years[0]
+ return slope
+
@property
def time_derivative_of_return_level(self):
if self.crashed:
return 0.0
else:
- # Compute the slope strength
- slope = self._slope_strength()
- # Delta T must in the same unit as were the parameter of slope mu1 and sigma1
- slope *= self.nb_years_for_quantile_evolution * \
- self.coordinates.transformed_distance_between_two_successive_years[0]
- return slope
+ return self.time_derivative_times_years(self.nb_years_for_quantile_evolution)
+
+ def relative_change_in_return_level(self, initial_year, final_year):
+ return_level_values = []
+ for year in [initial_year, final_year]:
+ gev_params = self.unconstrained_estimator.margin_function_from_fit.get_gev_params(coordinate=np.array([year]),
+ is_transformed=False)
+ return_level_values.append(gev_params.quantile(self.quantile_level))
+ change_until_final_year = self.time_derivative_times_years(nb_years=final_year - initial_year)
+ change_in_between = (return_level_values[1] - return_level_values[0])
+ np.testing.assert_almost_equal(change_until_final_year, change_in_between, decimal=5)
+ initial_return_level = return_level_values[0]
+ return 100 * change_until_final_year / initial_return_level
def _slope_strength(self):
raise NotImplementedError
--
GitLab