From dadce483813746a89f03702110975cf2fcc64952 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Mon, 20 Jul 2020 11:30:48 +0200
Subject: [PATCH] [exceeding] fix code to create a qqplot. (sort only once the
empirical quantiles have been computed). Create ALTITUDE_TO_GREY_MASSIF in
utils to ensure that the exceedance plots are always done with the same
massifs everywhere (and these massifs are only the one that are not grey for
the non-stationary gev & gumbel)
---
extreme_trend/abstract_gev_trend_test.py | 15 ++++---
...dy_visualizer_for_non_stationary_trends.py | 43 ++++++++++++-------
.../main_result_trends_and_return_levels.py | 14 +++---
.../plot_uncertainty_curves.py | 2 +-
.../plot_uncertainty_histogram.py | 3 +-
projects/exceeding_snow_loads/utils.py | 9 ++++
.../test_exceeding_snow_loads/test_results.py | 2 +-
7 files changed, 59 insertions(+), 29 deletions(-)
diff --git a/extreme_trend/abstract_gev_trend_test.py b/extreme_trend/abstract_gev_trend_test.py
index 85df6b20..c18c003c 100644
--- a/extreme_trend/abstract_gev_trend_test.py
+++ b/extreme_trend/abstract_gev_trend_test.py
@@ -369,7 +369,6 @@ class AbstractGevTrendTest(object):
def qqplot_wrt_standard_gumbel(self, massif_name, altitude):
ax = plt.gca()
- size = 15
standard_gumbel_quantiles = self.get_standard_gumbel_quantiles()
unconstrained_empirical_quantiles = self.compute_empirical_quantiles(self.unconstrained_estimator)
# constrained_empirical_quantiles = self.compute_empirical_quantiles(self.constrained_estimator)
@@ -385,9 +384,11 @@ class AbstractGevTrendTest(object):
ax.plot(standard_gumbel_quantiles, unconstrained_empirical_quantiles, linestyle='None',
label=label_generic + '(selected model is ${}$)'.format(self.label), marker='o')
- ax.set_xlabel("Standard Gumbel quantile", fontsize=size)
- ax.set_ylabel("Standard Empirical quantile", fontsize=size)
- ax.legend(loc='lower right', prop={'size': 10})
+ size_label = 20
+ ax.set_xlabel("Theoretical quantile", fontsize=size_label)
+ ax.set_ylabel("Empirical quantile", fontsize=size_label)
+ # size_legend = 14
+ # ax.legend(loc='lower right', prop={'size': size_legend})
ax.set_xlim(ax_lim)
ax.set_ylim(ax_lim)
@@ -396,7 +397,7 @@ class AbstractGevTrendTest(object):
ax.set_xticks(ticks)
ax.set_yticks(ticks)
# ax.grid()
- ax.tick_params(labelsize=size)
+ ax.tick_params(labelsize=15)
def get_standard_gumbel_quantiles(self):
# Standard Gumbel quantiles
@@ -413,12 +414,14 @@ class AbstractGevTrendTest(object):
def compute_empirical_quantiles(self, estimator):
empirical_quantiles = []
- for year, maximum in sorted(zip(self.years, self.maxima), key=lambda t: t[1]):
+ # for year, maximum in sorted(zip(self.years, self.maxima), key=lambda t: t[1]):
+ for year, maximum in zip(self.years, self.maxima):
gev_param = estimator.function_from_fit.get_params(
coordinate=np.array([year]),
is_transformed=False)
maximum_standardized = gev_param.gumbel_standardization(maximum)
empirical_quantiles.append(maximum_standardized)
+ empirical_quantiles = sorted(empirical_quantiles)
return empirical_quantiles
# For some visualizations
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 7f99652d..1058570d 100644
--- a/extreme_trend/visualizers/study_visualizer_for_non_stationary_trends.py
+++ b/extreme_trend/visualizers/study_visualizer_for_non_stationary_trends.py
@@ -17,7 +17,7 @@ from extreme_data.meteo_france_data.scm_models_data.abstract_extended_study impo
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 \
StudyVisualizer
-from projects.exceeding_snow_loads.utils import NON_STATIONARY_TREND_TEST_PAPER_1
+from projects.exceeding_snow_loads.utils import NON_STATIONARY_TREND_TEST_PAPER_1, ALTITUDE_TO_GREY_MASSIF
from extreme_trend.abstract_gev_trend_test import AbstractGevTrendTest
from extreme_trend.trend_test_one_parameter.gumbel_trend_test_one_parameter import \
GumbelLocationTrendTest, GevStationaryVersusGumbel, GumbelScaleTrendTest, GumbelVersusGumbel
@@ -82,6 +82,7 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
# or only the stationary gumbel model which correspond to the building standards
self.model_subsets_for_uncertainty = [ModelSubsetForUncertainty.stationary_gumbel,
ModelSubsetForUncertainty.non_stationary_gumbel_and_gev][:]
+ assert isinstance(self.model_subsets_for_uncertainty, list)
if self.uncertainty_methods is None:
self.uncertainty_methods = [ConfidenceIntervalMethodFromExtremes.my_bayes,
ConfidenceIntervalMethodFromExtremes.ci_mle][1:]
@@ -127,8 +128,19 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
return {m: EUROCODE_QUANTILE for m in self.massif_name_to_psnow.keys()}
@property
- def massif_names_fitted(self):
- return list(self.massif_name_to_years_and_maxima_for_model_fitting.keys())
+ def intersection_of_massif_names_fitted(self) -> List[str]:
+ all_massif = set(self.study.all_massif_names())
+ missing_massifs_given = set(ALTITUDE_TO_GREY_MASSIF[self.altitude])
+ if ModelSubsetForUncertainty.non_stationary_gumbel_and_gev in self.model_subsets_for_uncertainty:
+ # Check the assertion for the missing massifs
+ # For the paper 1, I should enter in this loop only for the ground snow load
+ # (because i should not use this ModelSubsetForUncertainty for the Eurocode)
+ massifs_found = set(self.massif_name_to_trend_test_that_minimized_aic.keys())
+ missing_massifs_found = all_massif - massifs_found
+ assert missing_massifs_found == missing_massifs_given
+ # Return by default this result
+ intersection_massif_names = all_massif - missing_massifs_given
+ return list(intersection_massif_names)
@cached_property
def massif_name_to_years_and_maxima_for_model_fitting(self):
@@ -160,14 +172,20 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
return massif_name_to_trend_test_that_minimized_aic
- def get_sorted_trend_test(self, massif_name):
+
+
+ def get_trend_trend_test(self, massif_name, trend_test_classes):
x, y = self.massif_name_to_years_and_maxima_for_model_fitting[massif_name]
- starting_year = None
quantile_level = self.massif_name_to_eurocode_quantile_level_in_practice[massif_name]
all_trend_test = [
- t(years=x, maxima=y, starting_year=starting_year, quantile_level=quantile_level,
+ t(years=x, maxima=y, starting_year=None, quantile_level=quantile_level,
fit_method=self.fit_method)
- for t in self.non_stationary_trend_test] # type: List[AbstractGevTrendTest]
+ for t in trend_test_classes] # type: List[AbstractGevTrendTest]
+ return all_trend_test
+
+ def get_sorted_trend_test(self, massif_name):
+ # Compute trend test
+ all_trend_test = self.get_trend_trend_test(massif_name, self.non_stationary_trend_test)
# Exclude GEV models whose shape parameter is not in the support of the prior distribution for GMLE
if self.select_only_acceptable_shape_parameter:
acceptable_shape_parameter = lambda s: -0.5 <= s <= 0.5 # physically acceptable prior
@@ -332,12 +350,7 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
model_subset_for_uncertainty=ModelSubsetForUncertainty.non_stationary_gumbel_and_gev) \
-> Dict[str, EurocodeConfidenceIntervalFromExtremes]:
# Compute for the uncertainty massif names
- massifs_names = set(self.massif_name_to_years_and_maxima_for_model_fitting.keys()). \
- intersection(self.uncertainty_massif_names)
- # Update the name of the massif (because some massifs might have been removed by anderson test)
- if model_subset_for_uncertainty is ModelSubsetForUncertainty.non_stationary_gumbel_and_gev\
- and self.select_only_model_that_pass_anderson_test:
- massifs_names = massifs_names.intersection(set(self.massif_name_to_trend_test_that_minimized_aic.keys()))
+ massifs_names = self.intersection_of_massif_names_fitted
arguments = [
[self.massif_name_to_years_and_maxima_for_model_fitting[m],
self.massif_name_and_model_subset_to_model_class(m, model_subset_for_uncertainty),
@@ -384,7 +397,7 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
triplet = [(massif_name_to_eurocode_region[massif_name],
self.massif_name_to_eurocode_values[massif_name],
self.triplet_to_eurocode_uncertainty[(ci_method, model_subset_for_uncertainty, massif_name)])
- for massif_name in self.massif_names_fitted]
+ for massif_name in self.intersection_of_massif_names_fitted]
# First array for histogram
a = 100 * np.array([(uncertainty.confidence_interval[0] > eurocode,
uncertainty.mean_estimate > eurocode,
@@ -422,7 +435,7 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
trend_test.qqplot_wrt_standard_gumbel(massif_name, self.altitude)
self.plot_name = 'qpplot_plot_{}_{}_{}'.format(self.altitude, massif_name,
trend_test.unconstrained_estimator_gev_params_last_year.shape)
- self.show_or_save_to_file(add_classic_title=False, no_title=True)
+ self.show_or_save_to_file(add_classic_title=False, no_title=True, tight_layout=True)
plt.close()
def return_level_plot(self, ax, ax2, massif_name, color=None):
diff --git a/projects/exceeding_snow_loads/section_results/main_result_trends_and_return_levels.py b/projects/exceeding_snow_loads/section_results/main_result_trends_and_return_levels.py
index efc2f2d4..6aba0ed6 100644
--- a/projects/exceeding_snow_loads/section_results/main_result_trends_and_return_levels.py
+++ b/projects/exceeding_snow_loads/section_results/main_result_trends_and_return_levels.py
@@ -11,11 +11,11 @@ from projects.exceeding_snow_loads.section_results.plot_uncertainty_histogram im
mpl.rcParams['text.usetex'] = True
mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
-from extreme_data.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal
+from extreme_data.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal, CrocusSnowLoadEurocode
from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import \
ConfidenceIntervalMethodFromExtremes
from extreme_trend.visualizers.study_visualizer_for_non_stationary_trends import \
- StudyVisualizerForNonStationaryTrends
+ StudyVisualizerForNonStationaryTrends, ModelSubsetForUncertainty
from extreme_trend.visualizers.utils import load_altitude_to_visualizer
from projects.exceeding_snow_loads.section_results.plot_uncertainty_curves import plot_uncertainty_massifs
from projects.exceeding_snow_loads.utils import paper_study_classes, paper_altitudes
@@ -87,12 +87,16 @@ def major_result():
# massif_names = ['Beaufortain', 'Vercors']
massif_names = None
study_classes = paper_study_classes[:1]
- model_subsets_for_uncertainty = None
# study_classes = [CrocusSnowLoad3Days, CrocusSnowLoad5Days, CrocusSnowLoad7Days][::-1]
- # altitudes = [300, 600, 900, 1800, 2700][:2]
+ altitudes = [300, 600, 900, 1800, 2700][:2]
+ altitudes = [300, 600, 900, 1200, 1500, 1800]
altitudes = paper_altitudes
- # altitudes = [900, 1800, 2700][:1]
+ # altitudes = [900, 1800, 270{{0][:1]
for study_class in study_classes:
+ if study_class == CrocusSnowLoadEurocode:
+ model_subsets_for_uncertainty = [ModelSubsetForUncertainty.stationary_gumbel]
+ else:
+ model_subsets_for_uncertainty = None
intermediate_result(altitudes, massif_names, model_subsets_for_uncertainty,
uncertainty_methods, study_class,
multiprocessing=True)
diff --git a/projects/exceeding_snow_loads/section_results/plot_uncertainty_curves.py b/projects/exceeding_snow_loads/section_results/plot_uncertainty_curves.py
index 61ea9ccc..78a3e7ba 100644
--- a/projects/exceeding_snow_loads/section_results/plot_uncertainty_curves.py
+++ b/projects/exceeding_snow_loads/section_results/plot_uncertainty_curves.py
@@ -154,7 +154,7 @@ def add_title(ax, eurocode_region, massif_name, non_stationary_context):
def plot_tdrl_bars(altitude_to_visualizer, ax, massif_name, valid_altitudes, legend_size, fontsize):
visualizers = [v for a, v in altitude_to_visualizer.items()
- if a in valid_altitudes and massif_name in v.massif_names_fitted]
+ if a in valid_altitudes and massif_name in v.intersection_of_massif_names_fitted]
if len(visualizers) > 0:
tdrl_values = [v.massif_name_to_tdrl_value[massif_name] for v in visualizers]
# Plot bars
diff --git a/projects/exceeding_snow_loads/section_results/plot_uncertainty_histogram.py b/projects/exceeding_snow_loads/section_results/plot_uncertainty_histogram.py
index 3b53c98b..d6d31e60 100644
--- a/projects/exceeding_snow_loads/section_results/plot_uncertainty_histogram.py
+++ b/projects/exceeding_snow_loads/section_results/plot_uncertainty_histogram.py
@@ -63,7 +63,8 @@ def plot_histogram(altitude_to_visualizer, model_subset_for_uncertainty):
ax_twiny.tick_params(labelsize=fontsize_label)
ax_twiny.set_xlim(ax.get_xlim())
ax_twiny.set_xticks(altitudes)
- nb_massif_names = [len(v.massif_names_fitted) for v in altitude_to_visualizer.values()]
+
+ nb_massif_names = [len(v.intersection_of_massif_names_fitted) for v in altitude_to_visualizer.values()]
ax_twiny.set_xticklabels(nb_massif_names)
ax_twiny.set_xlabel('Total number of massifs at each altitude (for the percentage)', fontsize=fontsize_label)
diff --git a/projects/exceeding_snow_loads/utils.py b/projects/exceeding_snow_loads/utils.py
index c28ee205..dd14d667 100644
--- a/projects/exceeding_snow_loads/utils.py
+++ b/projects/exceeding_snow_loads/utils.py
@@ -29,6 +29,15 @@ NON_STATIONARY_TREND_TEST_PAPER_1 = [GumbelVersusGumbel,
GevLocationAgainstGumbel, GevScaleAgainstGumbel,
GevLocationAndScaleTrendTestAgainstGumbel]
+ALTITUDE_TO_GREY_MASSIF = {
+ 300: ['Devoluy', 'Queyras', 'Champsaur', 'Grandes-Rousses', 'Ubaye', 'Pelvoux', 'Haute-Maurienne', 'Maurienne', 'Vanoise', 'Thabor', 'Mercantour', 'Haut_Var-Haut_Verdon', 'Haute-Tarentaise', 'Parpaillon', 'Belledonne', 'Oisans'],
+ 600: ['Queyras', 'Pelvoux', 'Haute-Maurienne', 'Mercantour', 'Haut_Var-Haut_Verdon', 'Thabor'],
+ 900: [],
+ 1200: [],
+ 1500: [],
+ 1800: [],
+}
+
NON_STATIONARY_TREND_TEST_PAPER_2 = [
# Gumbel models
GumbelVersusGumbel,
diff --git a/test/test_projects/test_exceeding_snow_loads/test_results.py b/test/test_projects/test_exceeding_snow_loads/test_results.py
index c636f28b..b7434aa2 100644
--- a/test/test_projects/test_exceeding_snow_loads/test_results.py
+++ b/test/test_projects/test_exceeding_snow_loads/test_results.py
@@ -14,7 +14,7 @@ class TestResults(unittest.TestCase):
def test_run_intermediate_results(self):
# Load data
- altitude_to_visualizer = load_altitude_to_visualizer(altitudes=[300, 600], massif_names=['Chartreuse'],
+ altitude_to_visualizer = load_altitude_to_visualizer(altitudes=[300, 600], massif_names=None,
study_class=CrocusSnowLoadTotal,
model_subsets_for_uncertainty=None,
uncertainty_methods=[
--
GitLab