From 883067a19011f43b9af8ac1f5fbcf0f49c0bbb8a Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Mon, 25 Nov 2019 16:02:39 +0100
Subject: [PATCH] [PAPER 1] add ticks to the colorbar for the trend. add
histogram & marker to the uncertainty plot
---
.../plot/create_shifted_cmap.py | 7 ++-
.../scm_models_data/abstract_study.py | 3 +-
.../eurocode_visualizer.py | 20 ++++++-
.../main_result_trends_and_return_levels.py | 17 +++---
...dy_visualizer_for_non_stationary_trends.py | 58 +++++++++++++++----
.../abstract_result_from_extremes.py | 6 +-
.../result_from_mle_extremes.py | 5 +-
7 files changed, 86 insertions(+), 30 deletions(-)
diff --git a/experiment/meteo_france_data/plot/create_shifted_cmap.py b/experiment/meteo_france_data/plot/create_shifted_cmap.py
index 8edc0375..73c4a88e 100644
--- a/experiment/meteo_france_data/plot/create_shifted_cmap.py
+++ b/experiment/meteo_france_data/plot/create_shifted_cmap.py
@@ -24,10 +24,13 @@ def get_shifted_map(vmin, vmax):
return shifted_cmap
-def create_colorbase_axis(ax, label, cmap, norm):
+def create_colorbase_axis(ax, label, cmap, norm, ticks_values_and_labels=None):
divider = make_axes_locatable(ax)
cax = divider.append_axes('right', size='5%', pad=0.0)
- cb = cbar.ColorbarBase(cax, cmap=cmap, norm=norm)
+ ticks = ticks_values_and_labels[0] if ticks_values_and_labels is not None else None
+ cb = cbar.ColorbarBase(cax, cmap=cmap, norm=norm, ticks=ticks)
+ 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)
return norm
diff --git a/experiment/meteo_france_data/scm_models_data/abstract_study.py b/experiment/meteo_france_data/scm_models_data/abstract_study.py
index e194878b..bb1429a9 100644
--- a/experiment/meteo_france_data/scm_models_data/abstract_study.py
+++ b/experiment/meteo_france_data/scm_models_data/abstract_study.py
@@ -306,6 +306,7 @@ class AbstractStudy(object):
massif_name_to_hatch_boolean_list=None,
norm=None,
massif_name_to_marker_style=None,
+ ticks_values_and_labels=None,
):
if ax is None:
ax = plt.gca()
@@ -391,7 +392,7 @@ class AbstractStudy(object):
# create the colorbar only at the end
if add_colorbar:
if len(set(values)) > 1:
- create_colorbase_axis(ax, label, cmap, norm)
+ create_colorbase_axis(ax, label, cmap, norm, ticks_values_and_labels=ticks_values_and_labels)
if axis_off:
plt.axis('off')
diff --git a/experiment/paper_past_snow_loads/result_trends_and_return_levels/eurocode_visualizer.py b/experiment/paper_past_snow_loads/result_trends_and_return_levels/eurocode_visualizer.py
index fcf66c6f..1c917c6e 100644
--- a/experiment/paper_past_snow_loads/result_trends_and_return_levels/eurocode_visualizer.py
+++ b/experiment/paper_past_snow_loads/result_trends_and_return_levels/eurocode_visualizer.py
@@ -38,9 +38,10 @@ def plot_uncertainty_massifs(altitude_to_visualizer: Dict[int, StudyVisualizerFo
"""
visualizer = list(altitude_to_visualizer.values())[0]
# Subdivide massif names in group of 3
- m = 3
+ m = 1
uncertainty_massif_names = visualizer.uncertainty_massif_names
n = (len(uncertainty_massif_names) // m) + 1
+ print('total nb of massif', n)
for i in list(range(n))[:]:
massif_names = uncertainty_massif_names[m * i: m * (i + 1)]
print(massif_names)
@@ -97,7 +98,7 @@ def plot_single_uncertainty_massif_and_non_stationary_context(ax, massif_name, n
""" Generic function that might be used by many other more global functions"""
altitudes = list(altitude_to_visualizer.keys())
visualizer = list(altitude_to_visualizer.values())[0]
- colors = ['tab:green', 'tab:olive']
+ colors = ['tab:green', 'tab:brown'][::-1]
alpha = 0.2
# Display the EUROCODE return level
eurocode_region = massif_name_to_eurocode_region[massif_name]()
@@ -106,14 +107,27 @@ def plot_single_uncertainty_massif_and_non_stationary_context(ax, massif_name, n
for j, (color, uncertainty_method) in enumerate(zip(colors, visualizer.uncertainty_methods)):
# Plot eurocode standards only for the first loop
if j == 0:
+ # Plot eurocode norm
eurocode_region.plot_eurocode_snow_load_on_ground_characteristic_value_variable_action(ax,
altitudes=altitudes)
+ # Plot bars of TDRL only in the non stationary case
+ if non_stationary_context:
+ tdrl_values = [v.massif_name_to_tdrl_value[massif_name] for v in altitude_to_visualizer.values()]
+ # Plot bars
+ colors = [v.massif_name_to_tdrl_color[massif_name] for v in altitude_to_visualizer.values()]
+ ax.bar(altitudes, tdrl_values, width=150, color=colors)
+ # Plot markers
+ markers_kwargs = [v.massif_name_to_marker_style(markersize=4)[massif_name]
+ for v, tdrl_value in zip(altitude_to_visualizer.values(), tdrl_values)]
+ for altitude, marker_kwargs, value in zip(altitudes, markers_kwargs, tdrl_values):
+ ax.plot([altitude], [value / 2], **marker_kwargs)
+
# Plot uncertainties
plot_valid_return_level_uncertainties(alpha, altitude_to_visualizer, altitudes, ax, color, massif_name,
non_stationary_context, uncertainty_method)
ax.legend(loc=2)
- ax.set_ylim([0.0, 16])
+ ax.set_ylim([-1, 16])
massif_name_str = massif_name.replace('_', ' ')
eurocode_region_str = get_display_name_from_object_type(type(eurocode_region))
is_non_stationary_model = non_stationary_context if isinstance(non_stationary_context,
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 6c867dfb..70ad497e 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
@@ -4,6 +4,7 @@ import os.path as op
import matplotlib as mpl
import matplotlib.pyplot as plt
+from experiment.meteo_france_data.plot.create_shifted_cmap import get_shifted_map
from experiment.meteo_france_data.scm_models_data.visualization.study_visualization.study_visualizer import \
StudyVisualizer
from experiment.paper_past_snow_loads.result_trends_and_return_levels.eurocode_visualizer import \
@@ -23,6 +24,7 @@ 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()
def intermediate_result(altitudes, massif_names=None,
@@ -67,16 +69,17 @@ def major_result():
if __name__ == '__main__':
# minor_result(altitude=1800)
- intermediate_result(altitudes=[1500, 1800], massif_names=['Chartreuse'],
- uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.ci_mle],
- non_stationary_uncertainty=[False])
+ # intermediate_result(altitudes=[1500, 1800], massif_names=['Chartreuse'],
+ # uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.ci_mle,
+ # ConfidenceIntervalMethodFromExtremes.ci_bayes],
+ # non_stationary_uncertainty=[True])
# intermediate_result(altitudes=[1500, 1800], massif_names=None,
# uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.ci_mle],
# non_stationary_uncertainty=[False])
# intermediate_result(altitudes=[300, 600, 900, 1200, 1500, 1800], massif_names=None,
# uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.ci_mle],
# non_stationary_uncertainty=[False])
- # intermediate_result(altitudes=[300, 600, 900, 1200, 1500, 1800], massif_names=['Vercors', 'Oisans', 'Devoluy'],
- # uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.ci_mle,
- # ConfidenceIntervalMethodFromExtremes.ci_bayes],
- # non_stationary_uncertainty=[False, True])
+ intermediate_result(altitudes=[300, 600, 900, 1200, 1500, 1800], massif_names=None,
+ uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.ci_mle,
+ ConfidenceIntervalMethodFromExtremes.ci_bayes],
+ non_stationary_uncertainty=[False, True])
diff --git a/experiment/paper_past_snow_loads/result_trends_and_return_levels/study_visualizer_for_non_stationary_trends.py b/experiment/paper_past_snow_loads/result_trends_and_return_levels/study_visualizer_for_non_stationary_trends.py
index 85272623..0bcacfd0 100644
--- a/experiment/paper_past_snow_loads/result_trends_and_return_levels/study_visualizer_for_non_stationary_trends.py
+++ b/experiment/paper_past_snow_loads/result_trends_and_return_levels/study_visualizer_for_non_stationary_trends.py
@@ -7,6 +7,7 @@ import numpy as np
from cached_property import cached_property
from experiment.eurocode_data.utils import EUROCODE_QUANTILE
+from experiment.meteo_france_data.plot.create_shifted_cmap import get_shifted_map, get_colors
from experiment.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
from experiment.meteo_france_data.scm_models_data.visualization.study_visualization.study_visualizer import \
StudyVisualizer
@@ -55,6 +56,8 @@ 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
+
# Utils
@cached_property
@@ -91,8 +94,9 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
massif_name_to_trend_test_that_minimized_aic = {}
for massif_name, (x, y) in self.massif_name_to_non_null_years_and_maxima.items():
quantile_level = self.massif_name_to_eurocode_quantile_level_in_practice[massif_name]
- non_stationary_trend_test = [t(years=x, maxima=y, starting_year=starting_year, quantile_level=quantile_level)
- for t in self.non_stationary_trend_test]
+ non_stationary_trend_test = [
+ t(years=x, maxima=y, starting_year=starting_year, quantile_level=quantile_level)
+ for t in self.non_stationary_trend_test]
trend_test_that_minimized_aic = sorted(non_stationary_trend_test, key=lambda t: t.aic)[0]
massif_name_to_trend_test_that_minimized_aic[massif_name] = trend_test_that_minimized_aic
return massif_name_to_trend_test_that_minimized_aic
@@ -103,32 +107,62 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
def max_abs_tdrl(self):
return max(abs(min(self.massif_name_to_tdrl_value.values())), max(self.massif_name_to_tdrl_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 plot_trends(self, max_abs_tdrl=None):
- if max_abs_tdrl is None:
- max_abs_tdrl = self.max_abs_tdrl
- assert max_abs_tdrl > 0
+ if max_abs_tdrl is not None:
+ self.global_max_abs_tdrl = max_abs_tdrl
self.study.visualize_study(massif_name_to_value=self.massif_name_to_tdrl_value,
- vmin=-max_abs_tdrl, vmax=max_abs_tdrl,
- replace_blue_by_white=False, axis_off=True, show_label=False,
+ replace_blue_by_white=False, axis_off=False, show_label=False,
add_colorbar=True,
- massif_name_to_marker_style=self.massif_name_to_marker_style,
- show=self.show)
+ massif_name_to_marker_style=self.massif_name_to_marker_style(),
+ massif_name_to_color=self.massif_name_to_tdrl_color,
+ cmap=self.cmap,
+ show=self.show,
+ ticks_values_and_labels=self.ticks_values_and_labels,
+ label=self.label_trend_bar)
self.plot_name = 'tdlr_trends'
self.show_or_save_to_file(add_classic_title=False, tight_layout=False, no_title=True)
plt.close()
+ @property
+ def label_trend_bar(self):
+ return 'Change in {} years for Eurocode quantile of SL (kN $m^-2$)'.format(AbstractGevTrendTest.nb_years_for_quantile_evolution)
+
+ @property
+ def ticks_values_and_labels(self):
+ graduation = 0.2
+ positive_ticks = []
+ tick = graduation
+ while tick < self._max_abs_tdrl:
+ positive_ticks.append(round(tick, 1))
+ tick += graduation
+ all_ticks_labels = [-t for t in positive_ticks] + [0] + positive_ticks
+ ticks_values = [(t / 2 * self._max_abs_tdrl) + 0.5 for t in all_ticks_labels]
+ return ticks_values, all_ticks_labels
+
@cached_property
def massif_name_to_tdrl_value(self):
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 cmap(self):
+ return get_shifted_map(-self._max_abs_tdrl, self._max_abs_tdrl)
- @property
- def massif_name_to_marker_style(self):
+ @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_marker_style(self, markersize=5):
d = {}
for m, t in self.massif_name_to_minimized_aic_non_stationary_trend_test.items():
d[m] = {'marker': self.non_stationary_trend_test_to_marker[type(t)],
'color': 'k',
- 'markersize': 5,
+ 'markersize': markersize,
'fillstyle': 'full' if t.is_significant else 'none'}
return d
diff --git a/extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_result_from_extremes.py b/extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_result_from_extremes.py
index 240d6535..10c67292 100644
--- a/extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_result_from_extremes.py
+++ b/extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_result_from_extremes.py
@@ -1,5 +1,6 @@
import numpy as np
import pandas as pd
+import rpy2
from rpy2 import robjects
from extreme_fit.distribution.gev.gev_params import GevParams
@@ -48,7 +49,10 @@ class AbstractResultFromExtremes(AbstractResultFromModelFit):
qcov = r("make.qcov")(self.result_from_fit,
**kwargs)
common_kwargs['qcov'] = qcov
- mean_estimate, confidence_interval = self._confidence_interval_method(common_kwargs, ci_method, return_period)
+ try:
+ mean_estimate, confidence_interval = self._confidence_interval_method(common_kwargs, ci_method, return_period)
+ except rpy2.rinterface.RRuntimeError:
+ mean_estimate, confidence_interval = np.nan, (np.nan, np.nan)
return mean_estimate, confidence_interval
def _confidence_interval_method(self, common_kwargs, ci_method, return_period):
diff --git a/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_mle_extremes.py b/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_mle_extremes.py
index 243eb2e9..56a167bf 100644
--- a/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_mle_extremes.py
+++ b/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_mle_extremes.py
@@ -24,10 +24,7 @@ class ResultFromMleExtremes(AbstractResultFromExtremes):
'method': method_name,
# xrange = NULL, nint = 20
}
- try:
- res = r.ci(self.result_from_fit, **mle_ci_parameters, **common_kwargs)
- except rpy2.rinterface.RRuntimeError:
- return np.nan, (np.nan, np.nan)
+ res = r.ci(self.result_from_fit, **mle_ci_parameters, **common_kwargs)
if self.is_non_stationary:
a = np.array(res)[0]
lower, mean_estimate, upper, _ = a
--
GitLab