From 86fe53c82591854033a1117d6be5303e5513bf5c Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Tue, 7 Apr 2020 16:22:32 +0200
Subject: [PATCH] [contrasting project] create
StudyVisualizerForReturnLevelChange.
---
.../scm_models_data/abstract_study.py | 4 +-
.../visualization/create_shifted_cmap.py | 43 ++++-
...dy_visualizer_for_non_stationary_trends.py | 12 +-
.../figure1_mean_ratio_return_level_ratio.py | 157 ++++++++++--------
4 files changed, 139 insertions(+), 77 deletions(-)
diff --git a/extreme_data/meteo_france_data/scm_models_data/abstract_study.py b/extreme_data/meteo_france_data/scm_models_data/abstract_study.py
index adbebbdf..805de9df 100644
--- a/extreme_data/meteo_france_data/scm_models_data/abstract_study.py
+++ b/extreme_data/meteo_france_data/scm_models_data/abstract_study.py
@@ -461,6 +461,7 @@ class AbstractStudy(object):
massif_name_to_marker_style=None,
marker_style_to_label_name=None,
ticks_values_and_labels=None,
+ fontsize_label=15,
):
if ax is None:
ax = plt.gca()
@@ -541,7 +542,8 @@ 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, ticks_values_and_labels=ticks_values_and_labels)
+ create_colorbase_axis(ax, label, cmap, norm, ticks_values_and_labels=ticks_values_and_labels,
+ fontsize=fontsize_label)
if axis_off:
plt.axis('off')
diff --git a/extreme_data/meteo_france_data/scm_models_data/visualization/create_shifted_cmap.py b/extreme_data/meteo_france_data/scm_models_data/visualization/create_shifted_cmap.py
index 9534f0fc..ed734262 100644
--- a/extreme_data/meteo_france_data/scm_models_data/visualization/create_shifted_cmap.py
+++ b/extreme_data/meteo_france_data/scm_models_data/visualization/create_shifted_cmap.py
@@ -24,7 +24,7 @@ def get_shifted_map(vmin, vmax, cmap=plt.cm.bwr):
return shifted_cmap
-def create_colorbase_axis(ax, label, cmap, norm, ticks_values_and_labels=None):
+def create_colorbase_axis(ax, label, cmap, norm, ticks_values_and_labels=None, fontsize=15):
divider = make_axes_locatable(ax)
cax = divider.append_axes('right', size='5%', pad=0.0)
ticks = ticks_values_and_labels[0] if ticks_values_and_labels is not None else None
@@ -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, fontsize=15)
+ cb.set_label(label, fontsize=fontsize)
return norm
@@ -68,6 +68,45 @@ def imshow_shifted(ax, gev_param_name, values, visualization_extend, mask_2D=Non
ax.imshow(masked_array, extent=visualization_extend, cmap=shifted_cmap, origin='lower')
+def ticks_values_and_labels_for_percentages(graduation, max_abs_change):
+ positive_ticks = []
+ tick = graduation
+ while tick < 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 / max_abs_change) + 1) / 2 for t in all_ticks_labels]
+ return ticks_values, all_ticks_labels
+
+
+def ticks_and_labels_centered_on_one(max_ratio, min_ratio):
+ """When we compute some ratio of two values.
+ Then if we want to make a plot, when the color change from blue to red is at 1,
+ then you should use this function"""
+ # Option to have a number of graduation constant
+ m = max(max_ratio / 1.0, 1.0 / min_ratio)
+ max_ratio = 1.0 * m
+ min_ratio = 1.0 / m
+ # Build the middle point
+ midpoint = (max_ratio - 1.0) / (max_ratio - 0)
+ graduation = 0.1
+ # Build lower graduation
+ n = int(np.math.floor((1.0 - min_ratio) / graduation)) + 1
+ a1 = midpoint / (1.0 - min_ratio)
+ b1 = midpoint - 1.0 * a1
+ xlist1 = [1.0 - i * graduation for i in range(n)]
+ y_list1 = [a1 * x + b1 for x in xlist1]
+ # Build upper graduation
+ n = int(np.math.floor((max_ratio - 1.0) / graduation)) + 1
+ xlist2 = [1.0 + i * graduation for i in range(n)]
+ a2 = (1 - midpoint) / (max_ratio - 1.0)
+ b2 = 1.0 - a2 * max_ratio
+ y_list2 = [a2 * x + b2 for x in xlist2]
+ labels = xlist1 + xlist2
+ ticks = y_list1 + y_list2
+ labels = [np.round(l, 1) for l in labels]
+ return labels, max_ratio, midpoint, min_ratio, ticks
+
# from: https://stackoverflow.com/questions/7404116/defining-the-midpoint-of-a-colormap-in-matplotlib/20528097
def shiftedColorMap(cmap, start=0, midpoint=0.5, stop=1.0, name='shiftedcmap'):
diff --git a/extreme_trend/visualizers/study_visualizer_for_non_stationary_trends.py b/extreme_trend/visualizers/study_visualizer_for_non_stationary_trends.py
index f9a2c83d..e55b047c 100644
--- a/extreme_trend/visualizers/study_visualizer_for_non_stationary_trends.py
+++ b/extreme_trend/visualizers/study_visualizer_for_non_stationary_trends.py
@@ -11,7 +11,8 @@ from extreme_data.eurocode_data.eurocode_region import C2, C1, E
from extreme_data.eurocode_data.massif_name_to_departement import massif_name_to_eurocode_region
from extreme_data.eurocode_data.utils import EUROCODE_QUANTILE, EUROCODE_RETURN_LEVEL_STR, \
YEAR_OF_INTEREST_FOR_RETURN_LEVEL
-from extreme_data.meteo_france_data.scm_models_data.visualization.create_shifted_cmap import get_shifted_map, get_colors
+from extreme_data.meteo_france_data.scm_models_data.visualization.create_shifted_cmap import get_shifted_map, \
+ get_colors, ticks_values_and_labels_for_percentages
from extreme_data.meteo_france_data.scm_models_data.abstract_extended_study import AbstractExtendedStudy
from extreme_data.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
from extreme_data.meteo_france_data.scm_models_data.visualization.study_visualizer import \
@@ -254,14 +255,7 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
@property
def ticks_values_and_labels(self):
- positive_ticks = []
- tick = self.graduation
- while tick < self._max_abs_change:
- positive_ticks.append(round(tick, 1))
- tick += self.graduation
- all_ticks_labels = [-t for t in positive_ticks] + [0] + positive_ticks
- ticks_values = [((t / self._max_abs_change) + 1) / 2 for t in all_ticks_labels]
- return ticks_values, all_ticks_labels
+ return ticks_values_and_labels_for_percentages(self.graduation, self._max_abs_change)
@property
def graduation(self):
diff --git a/projects/contrasting_trends_in_snow_loads/gorman_figures/figure1_mean_ratio_return_level_ratio.py b/projects/contrasting_trends_in_snow_loads/gorman_figures/figure1_mean_ratio_return_level_ratio.py
index 431316ba..996c4e72 100644
--- a/projects/contrasting_trends_in_snow_loads/gorman_figures/figure1_mean_ratio_return_level_ratio.py
+++ b/projects/contrasting_trends_in_snow_loads/gorman_figures/figure1_mean_ratio_return_level_ratio.py
@@ -5,7 +5,8 @@ import numpy as np
from extreme_data.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day
from extreme_data.meteo_france_data.scm_models_data.visualization.create_shifted_cmap import get_shifted_map, \
- get_colors, shiftedColorMap
+ get_colors, shiftedColorMap, ticks_values_and_labels_for_percentages
+from extreme_data.meteo_france_data.scm_models_data.visualization.study_visualizer import StudyVisualizer
from extreme_fit.estimator.margin_estimator.utils import fitted_stationary_gev
from extreme_fit.model.margin_model.linear_margin_model.abstract_temporal_linear_margin_model import \
TemporalMarginFitMethod
@@ -15,70 +16,96 @@ import matplotlib.pyplot as plt
from extreme_trend.visualizers.study_visualizer_for_non_stationary_trends import StudyVisualizerForNonStationaryTrends
-def mean_ratio(altitude):
- pass
-
-
-def massif_name_to_return_level(altitude, year_min, year_max, return_period):
- study = SafranSnowfall1Day(altitude=altitude, year_min=year_min, year_max=year_max)
- massif_name_to_return_levels = OrderedDict()
- for massif_name, maxima in study.massif_name_to_annual_maxima.items():
- gev_param = fitted_stationary_gev(maxima, fit_method=TemporalMarginFitMethod.extremes_fevd_mle)
- return_level = gev_param.return_level(return_period=return_period)
- massif_name_to_return_levels[massif_name] = return_level
- return massif_name_to_return_levels
-
-
-def plot_return_level_ratio(altitude, return_period=30, year_min=1959, year_middle=1989, year_max=2019):
- massif_name_to_return_level_past = massif_name_to_return_level(altitude, year_min, year_middle, return_period)
- massif_name_to_return_level_recent = massif_name_to_return_level(altitude, year_middle + 1, year_max,
- return_period)
- massif_name_to_ratio = {m: v / massif_name_to_return_level_past[m] for m, v in
- massif_name_to_return_level_recent.items()}
- max_ratio = max([e for e in massif_name_to_ratio.values()])
- min_ratio = min([e for e in massif_name_to_ratio.values()])
- # Option to have a number of graduation constant
- m = max(max_ratio / 1.0, 1.0 / min_ratio)
- max_ratio = 1.0 * m
- min_ratio = 1.0 / m
- # Build the middle point
- midpoint = (max_ratio - 1.0) / (max_ratio - 0)
- graduation = 0.1
- # Build lower graduation
- n = int(np.math.floor((1.0 - min_ratio) / graduation)) + 1
- a1 = midpoint / (1.0 - min_ratio)
- b1 = midpoint - 1.0 * a1
- xlist1 = [1.0 - i * graduation for i in range(n)]
- y_list1 = [a1 * x + b1 for x in xlist1]
- # Build upper graduation
- n = int(np.math.floor((max_ratio - 1.0) / graduation)) + 1
- xlist2 = [1.0 + i * graduation for i in range(n)]
- a2 = (1 - midpoint) / (max_ratio - 1.0)
- b2 = 1.0 - a2 * max_ratio
- y_list2 = [a2 * x + b2 for x in xlist2]
- labels = xlist1 + xlist2
- ticks = y_list1 + y_list2
- labels = [np.round(l, 1) for l in labels]
- cmap = shiftedColorMap(plt.cm.bwr, midpoint=midpoint, name='shifted')
- ax = plt.gca()
- massif_name_to_color = {m: get_colors([v], cmap, min_ratio, max_ratio)[0]
- for m, v in massif_name_to_ratio.items()}
- ticks_values_and_labels = ticks, labels
- AbstractStudy.visualize_study(ax=ax,
- massif_name_to_value=massif_name_to_ratio,
- massif_name_to_color=massif_name_to_color,
- replace_blue_by_white=True,
- axis_off=False,
- cmap=cmap,
- show_label=False,
- add_colorbar=True,
- show=False,
- vmin=min_ratio,
- vmax=max_ratio,
- ticks_values_and_labels=ticks_values_and_labels
- )
- plt.show()
+class StudyVisualizerForReturnLevelChange(StudyVisualizer):
+
+ def __init__(self, study_class, altitude, return_period=30, year_min=1959, year_middle=1989, year_max=2019, save_to_file=False,
+ fit_method=TemporalMarginFitMethod.extremes_fevd_mle):
+ self.fit_method = fit_method
+ self._year_min = year_min
+ self._year_middle = year_middle
+ self._year_max = year_max
+ self.study_before = study_class(altitude=altitude, year_min=self.year_min_before, year_max=self.year_max_before)
+ # Study will always refer to study before
+ super().__init__(self.study_before, save_to_file=save_to_file, show=not save_to_file)
+ self.study_after = study_class(altitude=altitude, year_min=self.year_min_after, year_max=self.year_max_after)
+ self.return_period = return_period
+
+ @property
+ def year_min_before(self):
+ return self._year_min
+
+ @property
+ def year_max_before(self):
+ return self._year_middle
+
+ @property
+ def year_min_after(self):
+ return self._year_middle + 1
+
+ @property
+ def year_max_after(self):
+ return self._year_max
+
+ def massif_name_to_return_level(self, study):
+ massif_name_to_return_levels = OrderedDict()
+ for massif_name, maxima in study.massif_name_to_annual_maxima.items():
+ gev_param = fitted_stationary_gev(maxima, fit_method=self.fit_method)
+ return_level = gev_param.return_level(return_period=self.return_period)
+ massif_name_to_return_levels[massif_name] = return_level
+ return massif_name_to_return_levels
+
+ def plot_return_level_percentage(self):
+ massif_name_to_return_level_past = self.massif_name_to_return_level(self.study_before)
+ massif_name_to_return_level_recent = self.massif_name_to_return_level(self.study_after)
+ massif_name_to_percentage = {
+ m: 100 * (v - massif_name_to_return_level_past[m]) / massif_name_to_return_level_past[m]
+ for m, v in
+ massif_name_to_return_level_recent.items()}
+
+ max_abs_change = max([abs(e) for e in massif_name_to_percentage.values()])
+ ticks, labels = ticks_values_and_labels_for_percentages(graduation=10, max_abs_change=max_abs_change)
+ min_ratio = -max_abs_change
+ max_ratio = max_abs_change
+ cmap = get_shifted_map(min_ratio, max_ratio)
+ ax = plt.gca()
+ massif_name_to_color = {m: get_colors([v], cmap, min_ratio, max_ratio)[0]
+ for m, v in massif_name_to_percentage.items()}
+ ticks_values_and_labels = ticks, labels
+ label = 'Relative change in {} return level between \n' \
+ 'a GEV fitted on {}-{} and a GEV fitted on {}-{} (%)'.format(self.return_period,
+ self.year_min_before,
+ self.year_max_before,
+ self.year_min_after,
+ self.year_max_after)
+ AbstractStudy.visualize_study(ax=ax,
+ massif_name_to_value=massif_name_to_percentage,
+ massif_name_to_color=massif_name_to_color,
+ replace_blue_by_white=True,
+ axis_off=False,
+ cmap=cmap,
+ show_label=False,
+ add_colorbar=True,
+ show=False,
+ vmin=min_ratio,
+ vmax=max_ratio,
+ ticks_values_and_labels=ticks_values_and_labels,
+ label=label,
+ fontsize_label=10,
+ )
+ ax.get_xaxis().set_visible(True)
+ ax.set_xticks([])
+ ax.set_xlabel('Altitude = {}m'.format(self.study.altitude), fontsize=15)
+ self.plot_name = 'change in return levels'
+ self.show_or_save_to_file(add_classic_title=False, tight_layout=True, no_title=True,
+ dpi=500)
+ plt.close()
if __name__ == '__main__':
- plot_return_level_ratio(altitude=1800)
+ for altitude in [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][:]:
+ study_visualizer = StudyVisualizerForReturnLevelChange(study_class=SafranSnowfall1Day, altitude=altitude,
+ return_period=30,
+ save_to_file=True,
+ fit_method=TemporalMarginFitMethod.extremes_fevd_l_moments)
+ study_visualizer.plot_return_level_percentage()
+
--
GitLab