diff --git a/extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py b/extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py
index cbccba7b6d0376e5f660d5ae647c7230d0ca5e50..78d84b8f401b9ce8ef93cc94dbaf7a13341a0c47 100644
--- a/extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py
+++ b/extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py
@@ -710,13 +710,14 @@ class StudyVisualizer(VisualizationParameters):
def plot_map(self, cmap, fit_method, graduation, label, massif_name_to_value, plot_name, add_x_label=True,
negative_and_positive_values=True, massif_name_to_text=None):
- load_plot(cmap, graduation, label, massif_name_to_value, self.study.altitude, fitmethod_to_str(fit_method),
- add_x_label=add_x_label, negative_and_positive_values=negative_and_positive_values,
- massif_name_to_text=massif_name_to_text)
- self.plot_name = plot_name
- # self.show_or_save_to_file(add_classic_title=False, tight_layout=True, no_title=True, dpi=500)
- self.show_or_save_to_file(add_classic_title=False, no_title=True, dpi=500)
- plt.close()
+ if len(massif_name_to_value) > 0:
+ load_plot(cmap, graduation, label, massif_name_to_value, self.study.altitude, fitmethod_to_str(fit_method),
+ add_x_label=add_x_label, negative_and_positive_values=negative_and_positive_values,
+ massif_name_to_text=massif_name_to_text)
+ self.plot_name = plot_name
+ # self.show_or_save_to_file(add_classic_title=False, tight_layout=True, no_title=True, dpi=500)
+ self.show_or_save_to_file(add_classic_title=False, no_title=True, dpi=500)
+ plt.close()
def plot_abstract(self, massif_name_to_value, label, plot_name, fit_method='', graduation=10.0, cmap=plt.cm.bwr,
add_x_label=True, negative_and_positive_values=True, massif_name_to_text=None):
diff --git a/extreme_trend/abstract_gev_trend_test.py b/extreme_trend/abstract_gev_trend_test.py
index c8a72be9b924e95784a04ff367a6e91ac29f8fb6..e205f35aa7ac202392f441135fd19d81c50589ca 100644
--- a/extreme_trend/abstract_gev_trend_test.py
+++ b/extreme_trend/abstract_gev_trend_test.py
@@ -492,4 +492,19 @@ class AbstractGevTrendTest(object):
old_mean = self.get_unconstrained_gev_params(year=self.years[-nb_years]).mean
return last_mean - old_mean
+ def relative_change_in_mean_for_the_last_x_years(self, nb_years):
+ last_mean = self.unconstrained_estimator_gev_params_last_year.mean
+ old_mean = self.get_unconstrained_gev_params(year=self.years[-nb_years]).mean
+ return 100 * (last_mean - old_mean) / old_mean
+
+ def change_in_50_year_return_level_for_the_last_x_years(self, nb_years, return_period=50):
+ last_mean = self.unconstrained_estimator_gev_params_last_year.return_level(return_period=return_period)
+ old_mean = self.get_unconstrained_gev_params(year=self.years[-nb_years]).return_level(return_period=return_period)
+ return last_mean - old_mean
+
+ def relative_change_in_50_year_return_level_for_the_last_x_years(self, nb_years, return_period=50):
+ last_mean = self.unconstrained_estimator_gev_params_last_year.return_level(return_period=return_period)
+ old_mean = self.get_unconstrained_gev_params(year=self.years[-nb_years]).return_level(return_period=return_period)
+ return 100 * (last_mean - old_mean) / old_mean
+
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 c426001f444d8e119517f96cc7e10a7e7e38a4f8..742e800f70c79d62bfc16ec93a7435dcad973a01 100644
--- a/extreme_trend/visualizers/study_visualizer_for_non_stationary_trends.py
+++ b/extreme_trend/visualizers/study_visualizer_for_non_stationary_trends.py
@@ -177,8 +177,12 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
trend_test_that_minimized_aic = sorted_trend_test[0]
massif_name_to_trend_test_that_minimized_aic[massif_name] = trend_test_that_minimized_aic
# Extract the stationary model that minimized AIC
- stationary_trend_test_that_minimized_aic = [t for t in sorted_trend_test if type(t) in
- [GumbelVersusGumbel, GevStationaryVersusGumbel]][0]
+ stationary_trend_tests_that_minimized_aic = [t for t in sorted_trend_test if type(t) in
+ [GumbelVersusGumbel, GevStationaryVersusGumbel]]
+ if len(stationary_trend_tests_that_minimized_aic) == 0:
+ stationary_trend_test_that_minimized_aic = None
+ else:
+ stationary_trend_test_that_minimized_aic = stationary_trend_tests_that_minimized_aic[0]
massif_name_to_stationary_trend_test_that_minimized_aic[
massif_name] = stationary_trend_test_that_minimized_aic
# Extract the Gumbel model that minimized AIC
diff --git a/extreme_trend/visualizers/utils.py b/extreme_trend/visualizers/utils.py
index a38b202f9fbbe6e318c3c1b9f195bfcee2612f07..f056e22aa07475e882098ed7204488fd2414595f 100644
--- a/extreme_trend/visualizers/utils.py
+++ b/extreme_trend/visualizers/utils.py
@@ -1,6 +1,6 @@
from collections import OrderedDict
-from extreme_data.meteo_france_data.scm_models_data.utils import Season
+from extreme_data.meteo_france_data.scm_models_data.utils import Season, FrenchRegion
from extreme_fit.model.margin_model.utils import \
MarginFitMethod
from extreme_trend.visualizers.study_visualizer_for_non_stationary_trends import \
@@ -12,11 +12,12 @@ def load_altitude_to_visualizer(altitudes, massif_names, model_subsets_for_uncer
study_visualizer_class=StudyVisualizerForNonStationaryTrends,
save_to_file=True,
multiprocessing=True,
- season=Season.annual):
+ season=Season.annual,
+ french_region=FrenchRegion.alps):
fit_method = MarginFitMethod.extremes_fevd_mle
altitude_to_visualizer = OrderedDict()
for altitude in altitudes:
- study = study_class(altitude=altitude, multiprocessing=multiprocessing, season=season)
+ study = study_class(altitude=altitude, multiprocessing=multiprocessing, season=season, french_region=french_region)
study_visualizer = study_visualizer_class(study=study, multiprocessing=multiprocessing,
save_to_file=save_to_file, uncertainty_massif_names=massif_names,
uncertainty_methods=uncertainty_methods,
diff --git a/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/main_snowfall_article.py b/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/main_snowfall_article.py
index c9908b22e14ebe8e8870f532e76194a41fb0bb1d..d46f3ec3fd85e18c251f06361351ec9f8acb009f 100644
--- a/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/main_snowfall_article.py
+++ b/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/main_snowfall_article.py
@@ -3,7 +3,7 @@ from multiprocessing.pool import Pool
import matplotlib as mpl
-from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day
+from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day, SafranPrecipitation1Day
from projects.contrasting_trends_in_snow_loads.article2_snowfall_versus_time_and_altitude.shape_plot import shape_plot
from projects.contrasting_trends_in_snow_loads.article2_snowfall_versus_time_and_altitude.snowfall_plot import \
plot_snowfall_mean, plot_snowfall_change_mean
@@ -74,25 +74,25 @@ def intermediate_result(altitudes, massif_names=None,
_ = compute_minimized_aic(visualizer)
# Plots
- validation_plot(altitude_to_visualizer, order_derivative=0)
- validation_plot(altitude_to_visualizer, order_derivative=1)
+ # validation_plot(altitude_to_visualizer, order_derivative=0)
+ # validation_plot(altitude_to_visualizer, order_derivative=1)
plot_snowfall_mean(altitude_to_visualizer)
plot_selection_curves(altitude_to_visualizer, paper1=False)
plot_snowfall_change_mean(altitude_to_visualizer)
- shape_plot(altitude_to_visualizer)
+ # shape_plot(altitude_to_visualizer)
def major_result():
uncertainty_methods = [ConfidenceIntervalMethodFromExtremes.ci_mle][:]
# massif_names = ['Beaufortain', 'Vercors']
massif_names = None
- study_classes = [SafranSnowfall1Day]
+ study_classes = [SafranSnowfall1Day, SafranPrecipitation1Day][::-1]
model_subsets_for_uncertainty = None
altitudes = paper_altitudes
- # altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900]
+ altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900]
# altitudes = [900, 1200, 1500, 1800][:2]
# altitudes = [1800, 2100, 2400, 2700][:2]
- altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000]
+ # altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000]
# altitudes = draft_altitudes
for study_class in study_classes:
diff --git a/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/snowfall_plot.py b/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/snowfall_plot.py
index e34f87553c6cac70a8ba356748e92e529e7dde0d..fdc653017405e032a246c4519af3ffc7b33ab85d 100644
--- a/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/snowfall_plot.py
+++ b/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/snowfall_plot.py
@@ -23,74 +23,87 @@ def fit_linear_regression(x, y):
def plot_snowfall_change_mean(altitude_to_visualizer: Dict[int, StudyVisualizerForMeanValues]):
visualizer = list(altitude_to_visualizer.values())[0]
study = visualizer.study
- # Plot the curve for the evolution of the mean
- massif_name_to_a, massif_name_to_b, massif_name_to_r2_score = plot_mean(altitude_to_visualizer, derivative=True)
- # Augmentation every km
- massif_name_to_augmentation_every_km = {m: a * 1000 for m, a in massif_name_to_a.items()}
- visualizer.plot_abstract_fast(massif_name_to_augmentation_every_km,
- label='Augmentation of time derivative of mean annual maxima of {}\n for every km of elevation ({})'.format(
- SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)], study.variable_unit),
- add_x_label=False)
- # Value at 2000 m
- massif_name_to_mean_at_2000 = {m: a * 2000 + massif_name_to_b[m] for m, a in massif_name_to_a.items()}
- visualizer.plot_abstract_fast(massif_name_to_mean_at_2000,
- label='Time derivative of mean annual maxima \nof {} at 2000 m ({})'.format(
- SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)], study.variable_unit),
- add_x_label=False,
- add_text=True, massif_name_to_text={k: str(v) for k,v in massif_name_to_r2_score.items()})
- # Altitude for the change of dynamic
- massif_name_to_altitude_change_dynamic = {m: - massif_name_to_b[m] / a for m, a in massif_name_to_a.items()}
- # Keep only those that are in a reasonable range
- massif_name_to_altitude_change_dynamic = {m: d for m, d in massif_name_to_altitude_change_dynamic.items()
- if 0 < d < 3000}
- visualizer.plot_abstract_fast(massif_name_to_altitude_change_dynamic,
- label='Altitude for the change of dynamic (m)',
- add_x_label=False, graduation=500,
- add_text=True, massif_name_to_text={k: str(v) for k,v in massif_name_to_r2_score.items()})
- # R2 score
- # visualizer.plot_abstract_fast(massif_name_to_r2_score, label='r2 time derivative of the mean', graduation=0.1,
- # add_x_label=False,
- # negative_and_positive_values=False,
- # )
+
+ variables = ['mean annual maxima', '50 year return level']
+ mean_indicators = [True, False]
+ for variable, mean_indicator in zip(variables, mean_indicators):
+ for relative in [False, True]:
+ if relative:
+ variable += str(' (relative)')
+ # Plot the curve for the evolution of the mean
+ massif_name_to_a, massif_name_to_b, massif_name_to_r2_score = plot_mean(altitude_to_visualizer,
+ derivative=True, mean_indicator=mean_indicator,
+ relative=relative)
+ # Augmentation every km
+ massif_name_to_augmentation_every_km = {m: a * 1000 for m, a in massif_name_to_a.items()}
+ massif_to_r2_score_text = {k: str(round(v, 2)) for k, v in massif_name_to_r2_score.items()}
+
+ visualizer.plot_abstract_fast(massif_name_to_augmentation_every_km,
+ label='Slope for changes of {} of {}\n for every km of elevation ({})'.format(
+ variable, SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)], study.variable_unit),
+ add_x_label=False)
+ # Value at 2000 m
+ massif_name_to_mean_at_2000 = {m: a * 2000 + massif_name_to_b[m] for m, a in massif_name_to_a.items()}
+ visualizer.plot_abstract_fast(massif_name_to_mean_at_2000,
+ label='Changes of {} \nof {} at 2000 m ({})'.format(
+ variable, SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)], study.variable_unit),
+ add_x_label=False,
+ add_text=True, massif_name_to_text=massif_to_r2_score_text)
+ # Altitude for the change of dynamic
+ massif_name_to_altitude_change_dynamic = {m: - massif_name_to_b[m] / a for m, a in massif_name_to_a.items()}
+ # Keep only those that are in a reasonable range
+ massif_name_to_altitude_change_dynamic = {m: d for m, d in massif_name_to_altitude_change_dynamic.items()
+ if 0 < d < 3000}
+ visualizer.plot_abstract_fast(massif_name_to_altitude_change_dynamic,
+ label=('Altitude when the changes equal zero for %s (m)' % variable),
+ add_x_label=False, graduation=500,
+ add_text=True, massif_name_to_text=massif_to_r2_score_text)
def plot_snowfall_mean(altitude_to_visualizer: Dict[int, StudyVisualizerForMeanValues]):
visualizer = list(altitude_to_visualizer.values())[0]
study = visualizer.study
- # Plot the curve for the evolution of the mean
- massif_name_to_a, massif_name_to_b, massif_name_to_r2_score = plot_mean(altitude_to_visualizer, derivative=False)
- # Compute values of interest
- massif_name_to_mean_at_2000 = {m: a * 2000 + massif_name_to_b[m] for m, a in massif_name_to_a.items()}
- massif_name_to_mean_at_1000 = {m: a * 1000 + massif_name_to_b[m] for m, a in massif_name_to_a.items()}
- massif_name_to_augmentation_every_km = {m: a * 1000 for m, a in massif_name_to_a.items()}
- massif_name_to_relative_augmentation_between_2000_and_3000_every_km = {m: 100 * (a * 1000) / (a * 2000 + massif_name_to_b[m]) for m, a in massif_name_to_a.items()}
- massif_name_to_relative_augmentation_between_1000_and_2000_every_km = {m: 100 * (a * 1000) / (a * 1000 + massif_name_to_b[m]) for m, a in massif_name_to_a.items()}
-
- # Augmentation every km
- ds = [massif_name_to_augmentation_every_km, massif_name_to_relative_augmentation_between_1000_and_2000_every_km,
- massif_name_to_relative_augmentation_between_2000_and_3000_every_km]
- prefixs = ['Augmentation', 'Relative augmentation between 1000m and 2000m', 'Relative augmentation between 2000m and 3000m']
- graduations = [1.0, 10.0, 10.0]
- for d, prefix in zip(ds, prefixs):
- visualizer.plot_abstract_fast(d,
- graduation=1.0,
- label=prefix + ' of mean annual maxima of {} \nfor every km of elevation ({})'.format(
- SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)], study.variable_unit),
+
+ variables = ['mean annual maxima', '50 year return level']
+ mean_indicators = [True, False]
+ for variable, mean_indicator in zip(variables, mean_indicators):
+ # Plot the curve for the evolution of the mean
+ massif_name_to_a, massif_name_to_b, massif_name_to_r2_score = plot_mean(altitude_to_visualizer,
+ derivative=False,
+ mean_indicator=mean_indicator)
+ # Compute values of interest
+ massif_name_to_mean_at_2000 = {m: a * 2000 + massif_name_to_b[m] for m, a in massif_name_to_a.items()}
+ massif_name_to_augmentation_every_km = {m: a * 1000 for m, a in massif_name_to_a.items()}
+ massif_name_to_relative_augmentation_between_2000_and_3000_every_km = {
+ m: 100 * (a * 1000) / (a * 2000 + massif_name_to_b[m]) for m, a in massif_name_to_a.items()}
+ massif_name_to_relative_augmentation_between_1000_and_2000_every_km = {
+ m: 100 * (a * 1000) / (a * 1000 + massif_name_to_b[m]) for m, a in massif_name_to_a.items()}
+ massif_to_r2_score_text = {k: str(round(v, 2)) for k, v in massif_name_to_r2_score.items()}
+ # Augmentation every km
+ ds = [massif_name_to_augmentation_every_km, massif_name_to_relative_augmentation_between_1000_and_2000_every_km,
+ massif_name_to_relative_augmentation_between_2000_and_3000_every_km]
+ prefixs = ['Augmentation', 'Relative augmentation between 1000m and 2000m',
+ 'Relative augmentation between 2000m and 3000m']
+ graduations = [1.0, 10.0, 10.0]
+ for d, prefix, graduation in zip(ds, prefixs, graduations):
+ visualizer.plot_abstract_fast(d,
+ graduation=10.0,
+ label=prefix + ' of {} of {} \nfor every km of elevation ({})'.format(
+ variable, SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)],
+ study.variable_unit),
+ add_x_label=False, negative_and_positive_values=False,
+ add_text=True,
+ massif_name_to_text=massif_to_r2_score_text
+ )
+ # Value at 2000 m
+ visualizer.plot_abstract_fast(massif_name_to_mean_at_2000, label='{} of {} at 2000 m ()'.format(
+ variable, SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)], study.variable_unit),
add_x_label=False, negative_and_positive_values=False,
- add_text=True,
- massif_name_to_text={k: str(v) for k, v in massif_name_to_r2_score.items()}
- )
- # Value at 2000 m
- visualizer.plot_abstract_fast(massif_name_to_mean_at_2000, label='Mean annual maxima of {} at 2000 m ()'.format(
- SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)], study.variable_unit),
- add_x_label=False, negative_and_positive_values=False,
- add_text=True, massif_name_to_text={k: str(v) for k,v in massif_name_to_r2_score.items()})
- # # R2 score
- # visualizer.plot_abstract_fast(massif_name_to_r2_score, label='r2 mean', graduation=0.1,
- # add_x_label=False, negative_and_positive_values=False)
-
-
-def plot_mean(altitude_to_visualizer: Dict[int, StudyVisualizerForMeanValues], derivative=False):
+ add_text=True, massif_name_to_text=massif_to_r2_score_text)
+
+
+def plot_mean(altitude_to_visualizer: Dict[int, StudyVisualizerForMeanValues], derivative=False, mean_indicator=True,
+ relative=False):
ax = plt.gca()
massif_name_to_linear_regression_result = {}
@@ -104,29 +117,44 @@ def plot_mean(altitude_to_visualizer: Dict[int, StudyVisualizerForMeanValues], d
if len(altitudes_massif) >= 2:
trend_tests = [altitude_to_visualizer[a].massif_name_to_trend_test_that_minimized_aic[massif_name]
for a in altitudes_massif]
+ nb_years = 50
+ return_period = 50
+ if mean_indicator:
+ indicator_str = 'mean'
+ else:
+ indicator_str = '50 year return level'
+ res = [(a, t)
+ for i, (a, t) in enumerate(zip(altitudes_massif, trend_tests))
+ if not t.unconstrained_model_is_stationary]
if derivative:
- nb_years = 10
- res = [(a, t.change_in_mean_for_the_last_x_years(nb_years=nb_years))
- for i, (a, t) in enumerate(zip(altitudes_massif, trend_tests))
- if t.is_significant]
+ if mean_indicator:
+ if relative:
+ res = [(a, t.relative_change_in_mean_for_the_last_x_years(nb_years=nb_years)) for a, t in res]
+ else:
+ res = [(a, t.change_in_mean_for_the_last_x_years(nb_years=nb_years)) for a, t in res]
+ else:
+ if relative:
+ res = [(a, t.relative_change_in_50_year_return_level_for_the_last_x_years(nb_years=nb_years, return_period=return_period)) for a, t in res]
+ else:
+ res = [(a, t.change_in_50_year_return_level_for_the_last_x_years(nb_years=nb_years, return_period=return_period)) for a, t in res]
if len(res) > 0:
altitudes_values, values = zip(*res)
else:
altitudes_values, values = [], []
- moment = 'Change in the last {} years \nfor significative models'.format(nb_years)
- # res = [(a, t.change_in_mean_for_the_last_x_years(nb_years=nb_years))
- # for i, (a, t) in enumerate(zip(altitudes_massif, trend_tests))
- # if not t.unconstrained_model_is_stationary]
- # altitudes_values, values = zip(*res)
- # moment = 'Change in the last {} years \nfor non-stationary models'.format(nb_years)
+ indicator_str = 'Change of the {} for the last {} years \nfor non-stationary models'.format(indicator_str, nb_years)
+ if relative:
+ indicator_str += ' (relative)'
else:
- moment = 'mean'
- values = [t.unconstrained_estimator_gev_params_last_year.mean for t in trend_tests]
+ if mean_indicator:
+ values = [t.unconstrained_estimator_gev_params_last_year.mean for t in trend_tests]
+ else:
+ values = [t.unconstrained_estimator_gev_params_last_year.return_level(return_period=return_period)
+ for t in trend_tests]
altitudes_values = altitudes_massif
# Plot
if len(altitudes_values) >= 2:
massif_name_to_linear_regression_result[massif_name] = fit_linear_regression(altitudes_values, values)
- plot_values_against_altitudes(ax, altitudes_values, massif_id, massif_name, moment, study, values,
+ plot_values_against_altitudes(ax, altitudes_values, massif_id, massif_name, indicator_str, study, values,
visualizer)
ax.legend(prop={'size': 7}, ncol=3)
visualizer.show_or_save_to_file(dpi=500, add_classic_title=False)