diff --git a/extreme_data/meteo_france_data/scm_models_data/visualization/plot_utils.py b/extreme_data/meteo_france_data/scm_models_data/visualization/plot_utils.py
new file mode 100644
index 0000000000000000000000000000000000000000..341b53cec5697a81fe6051ee1907b145d415ed4e
--- /dev/null
+++ b/extreme_data/meteo_france_data/scm_models_data/visualization/plot_utils.py
@@ -0,0 +1,36 @@
+import matplotlib.pyplot as plt
+
+from extreme_data.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
+from extreme_data.meteo_france_data.scm_models_data.visualization.create_shifted_cmap import \
+ ticks_values_and_labels_for_percentages, get_shifted_map, get_colors
+
+
+def load_plot(cmap, graduation, label, massif_name_to_value, altitude, fit_method):
+ max_abs_change = max([abs(e) for e in massif_name_to_value.values()])
+ ticks, labels = ticks_values_and_labels_for_percentages(graduation=graduation, max_abs_change=max_abs_change)
+ min_ratio = -max_abs_change
+ max_ratio = max_abs_change
+ cmap = get_shifted_map(min_ratio, max_ratio, cmap)
+ massif_name_to_color = {m: get_colors([v], cmap, min_ratio, max_ratio)[0]
+ for m, v in massif_name_to_value.items()}
+ ticks_values_and_labels = ticks, labels
+ ax = plt.gca()
+ AbstractStudy.visualize_study(ax=ax,
+ massif_name_to_value=massif_name_to_value,
+ 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(altitude), fontsize=15)
+ ax.set_title('Fit method is {}'.format(fit_method))
\ No newline at end of file
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 603b40200e999176155211a97e5c09291c179be4..38481079b2ff9e48f14a87d38ca12b35452deb9d 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
@@ -12,7 +12,9 @@ import numpy as np
import pandas as pd
import seaborn as sns
+from extreme_data.meteo_france_data.scm_models_data.visualization.plot_utils import load_plot
from extreme_fit.estimator.margin_estimator.utils import fitted_stationary_gev
+from extreme_fit.model.margin_model.utils import fitmethod_to_str
from extreme_fit.model.result_from_model_fit.result_from_extremes.eurocode_return_level_uncertainties import \
EurocodeConfidenceIntervalFromExtremes, compute_eurocode_confidence_interval
from extreme_data.meteo_france_data.scm_models_data.abstract_extended_study import AbstractExtendedStudy
@@ -320,7 +322,6 @@ class StudyVisualizer(VisualizationParameters):
zip([massif_name for _, massif_name in massif_ids_and_names], altitudes_and_res))
return massif_name_to_eurocode_return_level_uncertainty
-
def visualize_all_mean_and_max_graphs(self):
specified_massif_ids = [self.study.study_massif_names.index(massif_name)
for massif_name in
@@ -409,7 +410,8 @@ class StudyVisualizer(VisualizationParameters):
tuples_x_y = [(year, np.mean(data[:, massif_id])) for year, data in
self.study.year_to_daily_time_serie_array.items()]
x, y = list(zip(*tuples_x_y))
- x, y = self.average_smoothing_with_sliding_window(x, y, window_size_for_smoothing=self.window_size_for_smoothing)
+ x, y = self.average_smoothing_with_sliding_window(x, y,
+ window_size_for_smoothing=self.window_size_for_smoothing)
ax.plot(x, y, color=color_mean)
ax.set_ylabel('mean'.format(self.window_size_for_smoothing), color=color_mean)
massif_name = self.study.study_massif_names[massif_id]
@@ -423,7 +425,8 @@ class StudyVisualizer(VisualizationParameters):
tuples_x_y = [(year, annual_maxima[massif_id]) for year, annual_maxima in
self.study.year_to_annual_maxima.items()]
x, y = list(zip(*tuples_x_y))
- x, y = self.average_smoothing_with_sliding_window(x, y, window_size_for_smoothing=self.window_size_for_smoothing)
+ x, y = self.average_smoothing_with_sliding_window(x, y,
+ window_size_for_smoothing=self.window_size_for_smoothing)
self.massif_id_to_smooth_maxima[massif_id] = (x, y)
return self.massif_id_to_smooth_maxima[massif_id]
@@ -701,3 +704,21 @@ class StudyVisualizer(VisualizationParameters):
x = np.array(x)[nb_to_delete:-nb_to_delete]
assert len(x) == len(y), "{} vs {}".format(len(x), len(y))
return x, y
+
+ # PLot functions that should be common
+
+ def plot_map(self, cmap, fit_method, graduation, label, massif_name_to_value, plot_name):
+ load_plot(cmap, graduation, label, massif_name_to_value, self.study.altitude, fitmethod_to_str(fit_method))
+ 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):
+ # Regroup the plot by altitudes
+ plot_name1 = '{}/{}'.format(self.study.altitude, plot_name)
+ # Regroup the plot by type of plot also
+ plot_name2 = '{}/{}'.format(plot_name.split()[0], plot_name)
+ for plot_name in [plot_name1, plot_name2]:
+ self.plot_map(cmap, fit_method, graduation, label, massif_name_to_value, plot_name)
+
diff --git a/extreme_fit/distribution/gev/gev_params.py b/extreme_fit/distribution/gev/gev_params.py
index 2638464108741d18a5723e31d09c8787c10bffa0..fce8d5444abee0a1cda2b972512ff028c5fbfce2 100644
--- a/extreme_fit/distribution/gev/gev_params.py
+++ b/extreme_fit/distribution/gev/gev_params.py
@@ -3,7 +3,7 @@ import matplotlib.pyplot as plt
from typing import List
from cached_property import cached_property
-from mpmath import euler, pi
+from mpmath import euler
from extreme_fit.distribution.abstract_extreme_params import AbstractExtremeParams
from extreme_fit.distribution.abstract_params import AbstractParams
@@ -86,13 +86,15 @@ class GevParams(AbstractExtremeParams):
@property
def mean(self) -> float:
if self.has_undefined_parameters:
- return np.nan
+ mean = np.nan
elif self.shape >= 1:
- return np.inf
+ mean = np.inf
elif self.shape == 0:
- return self.location + self.scale * euler
+ mean = self.location + self.scale * float(euler)
else:
- return self.location + self.scale * (self.g(k=1) - 1) / self.shape
+ mean = self.location + self.scale * (self.g(k=1) - 1) / self.shape
+ assert isinstance(mean, float)
+ return mean
@property
def variance(self) -> float:
@@ -101,7 +103,7 @@ class GevParams(AbstractExtremeParams):
elif self.shape >= 0.5:
return np.inf
elif self.shape == 0.0:
- return (self.scale * pi) ** 2 / 6
+ return (self.scale * np.pi) ** 2 / 6
else:
return ((self.scale / self.shape) ** 2) * (self.g(k=2) - self.g(k=1) ** 2)
diff --git a/extreme_fit/model/margin_model/utils.py b/extreme_fit/model/margin_model/utils.py
index 0985928a7c51251dc05001cd6ed863109df8eda5..e713e18c5bbc84d23784c23e261763796822f4a0 100644
--- a/extreme_fit/model/margin_model/utils.py
+++ b/extreme_fit/model/margin_model/utils.py
@@ -18,5 +18,6 @@ FEVD_MARGIN_FIT_METHOD_TO_ARGUMENT_STR = {
}
FEVD_MARGIN_FIT_METHODS = set(FEVD_MARGIN_FIT_METHOD_TO_ARGUMENT_STR.keys())
+
def fitmethod_to_str(fit_method):
- return str(fit_method).split('.')[-1]
+ return ' '.join(str(fit_method).split('.')[-1].split('_'))
diff --git a/extreme_trend/abstract_gev_trend_test.py b/extreme_trend/abstract_gev_trend_test.py
index bcee5a5ec132c80fd28e6c3085997097ddfdf8be..3ba905635a49b0ca019011c23c58efb6547d336b 100644
--- a/extreme_trend/abstract_gev_trend_test.py
+++ b/extreme_trend/abstract_gev_trend_test.py
@@ -46,7 +46,7 @@ class AbstractGevTrendTest(object):
@cached_property
def unconstrained_estimator(self):
return fitted_linear_margin_estimator(self.unconstrained_model_class, self.coordinates, self.dataset,
- self.starting_year, self.fit_method)
+ self.starting_year, self.fit_method)
# Likelihood ratio test
@@ -108,9 +108,7 @@ class AbstractGevTrendTest(object):
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.function_from_fit.get_params(
- coordinate=np.array([year]),
- is_transformed=False)
+ gev_params = self.get_unconstrained_gev_params(year)
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])
@@ -163,8 +161,6 @@ class AbstractGevTrendTest(object):
ax.plot(standard_gumbel_quantiles, sorted_maxima, linestyle='None',
label='Empirical model', marker='o', color='black')
-
-
ax_twiny = ax.twiny()
return_periods = [10, 25, 50]
quantiles = self.get_standard_quantiles_for_return_periods(return_periods, psnow)
@@ -182,11 +178,14 @@ class AbstractGevTrendTest(object):
end_real_proba = 1 - (0.02 / psnow)
stationary = True
if stationary:
- self.plot_model(ax, None, end_proba=end_real_proba, label='Selected model\nwhich is ${}$'.format(self.label),
+ self.plot_model(ax, None, end_proba=end_real_proba,
+ label='Selected model\nwhich is ${}$'.format(self.label),
color='grey')
else:
- self.plot_model(ax, 1959, end_proba=end_real_proba, label='Selected model, which is ${}$'.format(self.label))
- self.plot_model(ax, 2019, end_proba=end_real_proba, label='Selected model, which is ${}$'.format(self.label))
+ self.plot_model(ax, 1959, end_proba=end_real_proba,
+ label='Selected model, which is ${}$'.format(self.label))
+ self.plot_model(ax, 2019, end_proba=end_real_proba,
+ label='Selected model, which is ${}$'.format(self.label))
# Plot for the discarded model
# if 'Verdon' in massif_name and altitude == 300:
@@ -204,9 +203,6 @@ class AbstractGevTrendTest(object):
# + 'with $\zeta_0=0.84$',
# color=color)
-
-
-
ax_lim = [-1.5, 4]
ax.set_xlim(ax_lim)
ax_twiny.set_xlim(ax_lim)
@@ -231,13 +227,17 @@ class AbstractGevTrendTest(object):
label = 'Y({})'.format(year) if year is not None else label
if year is None:
year = 2019
- gev_params_year = self.unconstrained_estimator.function_from_fit.get_params(
- coordinate=np.array([year]),
- is_transformed=False)
+ gev_params_year = self.get_unconstrained_gev_params(year)
extended_maxima = [gev_params_year.gumbel_inverse_standardization(q) for q in extended_quantiles]
ax.plot(extended_quantiles, extended_maxima, linestyle='-', label=label, color=color, linewidth=5)
+ def get_unconstrained_gev_params(self, year):
+ gev_params_year = self.unconstrained_estimator.function_from_fit.get_params(
+ coordinate=np.array([year]),
+ is_transformed=False)
+ return gev_params_year
+
def linear_extension(self, ax, q, quantiles, sorted_maxima):
# Extend the curve linear a bit if the return period 50 is not in the quantiles
def compute_slope_intercept(x, y):
@@ -368,9 +368,22 @@ class AbstractGevTrendTest(object):
plt.gca().set_ylim(bottom=0)
def get_gev_params_with_big_shape_and_correct_shape(self):
- gev_params = self.unconstrained_estimator.function_from_fit.get_params(coordinate=np.array([YEAR_OF_INTEREST_FOR_RETURN_LEVEL]),
- is_transformed=False) # type: GevParams
+ gev_params = self.unconstrained_estimator.function_from_fit.get_params(
+ coordinate=np.array([YEAR_OF_INTEREST_FOR_RETURN_LEVEL]),
+ is_transformed=False) # type: GevParams
gev_params_with_corrected_shape = GevParams(loc=gev_params.location,
scale=gev_params.scale,
shape=0.5)
return gev_params, gev_params_with_corrected_shape
+
+ # Mean value
+
+ @property
+ def unconstrained_average_mean_value(self) -> float:
+ mean_values = []
+ for year in self.years:
+ mean = self.get_unconstrained_gev_params(year).mean
+ mean_values.append(mean)
+ average_mean_value = np.mean(mean_values)
+ assert isinstance(average_mean_value, float)
+ return average_mean_value
diff --git a/projects/altitude_spatial_model/altitudes_fit/altitudes_studies.py b/projects/altitude_spatial_model/altitudes_fit/altitudes_studies.py
index e0772e2eccbf3d86493308ecd81eacc2dab4c695..284a762252bde6c79d27e68012fd2c347e18ecaf 100644
--- a/projects/altitude_spatial_model/altitudes_fit/altitudes_studies.py
+++ b/projects/altitude_spatial_model/altitudes_fit/altitudes_studies.py
@@ -159,7 +159,8 @@ class AltitudesStudies(object):
altitudes.append(altitude)
self.plot_against_altitude(altitudes, ax, massif_id, massif_name, mean_moment)
- def plot_against_altitude(self, altitudes, ax, massif_id, massif_name, mean_moment):
+ @staticmethod
+ def plot_against_altitude(altitudes, ax, massif_id, massif_name, mean_moment):
di = massif_id // 8
if di == 0:
linestyle = '-'
diff --git a/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/main_snowfall_article.py b/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/main_snowfall_article.py
new file mode 100644
index 0000000000000000000000000000000000000000..a2b9366d2dcc9b80643de8b6a1640c11a40fe00d
--- /dev/null
+++ b/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/main_snowfall_article.py
@@ -0,0 +1,95 @@
+
+from multiprocessing.pool import Pool
+
+import matplotlib as mpl
+
+from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day
+from projects.contrasting_trends_in_snow_loads.snowfall_versus_time_and_altitude.study_visualizer_for_mean_values import \
+ StudyVisualizerForMeanValues
+from projects.contrasting_trends_in_snow_loads.snowfall_versus_time_and_altitude.validation_plot import validation_plot
+from projects.exceeding_snow_loads.section_results.plot_trend_curves import plot_trend_map
+
+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_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
+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
+from root_utils import NB_CORES
+
+
+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()
+
+
+def compute_minimized_aic(visualizer):
+ _ = visualizer.massif_name_to_trend_test_that_minimized_aic
+ return True
+
+
+def intermediate_result(altitudes, massif_names=None,
+ model_subsets_for_uncertainty=None, uncertainty_methods=None,
+ study_class=SafranSnowfall1Day,
+ multiprocessing=False):
+ """
+ Plot all the trends for all altitudes
+ And enable to plot uncertainty plot for some specific massif_names, uncertainty methods to be fast
+ :param altitudes:
+ :param massif_names:
+ :param non_stationary_uncertainty:
+ :param uncertainty_methods:
+ :param study_class:
+ :return:
+ """
+ # Load altitude to visualizer
+ altitude_to_visualizer = load_altitude_to_visualizer(altitudes, massif_names, model_subsets_for_uncertainty,
+ study_class, uncertainty_methods,
+ study_visualizer_class=StudyVisualizerForMeanValues)
+ # Load variable object efficiently
+ for v in altitude_to_visualizer.values():
+ _ = v.study.year_to_variable_object
+ # Compute minimized value efficiently
+ visualizers = list(altitude_to_visualizer.values())
+ if multiprocessing:
+ with Pool(NB_CORES) as p:
+ _ = p.map(compute_minimized_aic, visualizers)
+ else:
+ for visualizer in visualizers:
+ _ = compute_minimized_aic(visualizer)
+
+ # Plots
+ validation_plot(altitude_to_visualizer)
+
+
+def major_result():
+ uncertainty_methods = [ConfidenceIntervalMethodFromExtremes.ci_mle][:]
+ massif_names = ['Beaufortain', 'Vercors']
+ massif_names = None
+ study_classes = [SafranSnowfall1Day]
+ model_subsets_for_uncertainty = None
+ altitudes = paper_altitudes
+ # altitudes = [900, 1200]
+
+ for study_class in study_classes:
+ intermediate_result(altitudes, massif_names, model_subsets_for_uncertainty,
+ uncertainty_methods, study_class)
+
+
+if __name__ == '__main__':
+ major_result()
+ # intermediate_result(altitudes=[300], massif_names=None,
+ # uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
+ # ConfidenceIntervalMethodFromExtremes.ci_mle][1:],
+ # multiprocessing=True)
diff --git a/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/study_visualizer_for_mean_values.py b/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/study_visualizer_for_mean_values.py
new file mode 100644
index 0000000000000000000000000000000000000000..84dd12bedd34de8d307567a67ca2ddd0d83171d7
--- /dev/null
+++ b/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/study_visualizer_for_mean_values.py
@@ -0,0 +1,59 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from cached_property import cached_property
+
+from extreme_data.eurocode_data.utils import YEAR_OF_INTEREST_FOR_RETURN_LEVEL
+from extreme_data.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
+from extreme_fit.model.margin_model.utils import MarginFitMethod
+from extreme_trend.visualizers.study_visualizer_for_non_stationary_trends import StudyVisualizerForNonStationaryTrends
+
+
+class StudyVisualizerForMeanValues(StudyVisualizerForNonStationaryTrends):
+
+ def __init__(self, study: AbstractStudy, show=True, save_to_file=False, only_one_graph=False, only_first_row=False,
+ vertical_kde_plot=False, year_for_kde_plot=None, plot_block_maxima_quantiles=False,
+ temporal_non_stationarity=False, transformation_class=None, verbose=False, multiprocessing=False,
+ complete_non_stationary_trend_analysis=False, normalization_under_one_observations=True,
+ uncertainty_methods=None, model_subsets_for_uncertainty=None, uncertainty_massif_names=None,
+ effective_temporal_covariate=YEAR_OF_INTEREST_FOR_RETURN_LEVEL, relative_change_trend_plot=True,
+ non_stationary_trend_test_to_marker=None, fit_method=MarginFitMethod.extremes_fevd_mle,
+ select_only_acceptable_shape_parameter=True, fit_gev_only_on_non_null_maxima=False,
+ fit_only_time_series_with_ninety_percent_of_non_null_values=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, uncertainty_methods, model_subsets_for_uncertainty,
+ uncertainty_massif_names, effective_temporal_covariate, relative_change_trend_plot,
+ non_stationary_trend_test_to_marker, fit_method, select_only_acceptable_shape_parameter,
+ fit_gev_only_on_non_null_maxima, fit_only_time_series_with_ninety_percent_of_non_null_values)
+
+ @cached_property
+ def massif_name_to_empirical_mean(self):
+ massif_name_to_empirical_value = {}
+ for massif_name, maxima in self.study.massif_name_to_annual_maxima.items():
+ massif_name_to_empirical_value[massif_name] = np.mean(maxima)
+ return massif_name_to_empirical_value
+
+ @cached_property
+ def massif_name_to_parametric_mean(self):
+ massif_name_to_parameter_value = {}
+ for massif_name, trend_test in self.massif_name_to_trend_test_that_minimized_aic.items():
+ parameter_value = trend_test.unconstrained_average_mean_value
+ print(type(parameter_value), parameter_value)
+ massif_name_to_parameter_value[massif_name] = parameter_value
+ return massif_name_to_parameter_value
+
+ @cached_property
+ def massif_name_to_relative_difference_for_mean(self):
+ massif_name_to_relative_difference = {}
+ for massif_name in self.massif_name_to_trend_test_that_minimized_aic.keys():
+ e = self.massif_name_to_empirical_mean[massif_name]
+ p = self.massif_name_to_parametric_mean[massif_name]
+ relative_diference = 100 * (p-e) / e
+ massif_name_to_relative_difference[massif_name] = relative_diference
+ return massif_name_to_relative_difference
+
+ def plot_abstract_fast(self, massif_name_to_value, label, graduation=10.0, cmap=plt.cm.coolwarm):
+ super().plot_abstract(massif_name_to_value, label, label, self.fit_method, graduation, cmap)
+
+
diff --git a/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/validation_plot.py b/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/validation_plot.py
new file mode 100644
index 0000000000000000000000000000000000000000..58d27a2154c803ea699b0253db7b74e45b31d770
--- /dev/null
+++ b/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/validation_plot.py
@@ -0,0 +1,49 @@
+from typing import Dict
+
+import matplotlib.pyplot as plt
+
+from projects.contrasting_trends_in_snow_loads.snowfall_versus_time_and_altitude.study_visualizer_for_mean_values import \
+ StudyVisualizerForMeanValues
+from projects.ogorman.gorman_figures.figure1.study_visualizer_for_double_stationary_fit import \
+ StudyVisualizerForReturnLevelChange
+
+
+def validation_plot(altitude_to_visualizer: Dict[int, StudyVisualizerForMeanValues]):
+ # Plot the mean empirical, the mean parametric and the relative difference between the two
+ altitudes = list(altitude_to_visualizer.keys())
+ study_visualizer = list(altitude_to_visualizer.values())[0]
+ altitude_to_relative_differences = {}
+ # Plot map for the repartition of the difference
+ for altitude, visualizer in altitude_to_visualizer.items():
+ altitude_to_relative_differences[altitude] = plot_relative_difference_map(visualizer)
+ study_visualizer.show_or_save_to_file(add_classic_title=False, dpi=500)
+ # Shoe plot with respect to the altitude.
+ plot_shoe_relative_differences_distribution(altitude_to_relative_differences, altitudes)
+ study_visualizer.show_or_save_to_file(add_classic_title=False, dpi=500)
+
+
+def plot_shoe_relative_differences_distribution(altitude_to_relative_differences, altitudes):
+ ax = plt.gca()
+ width = 150
+ ax.boxplot([altitude_to_relative_differences[a] for a in altitudes], positions=altitudes, widths=width)
+ ax.set_xlim([min(altitudes) - width, max(altitudes) + width])
+ ylabel = 'Relative difference between empirical mean and parametric mean (%)'
+ ax.set_ylabel(ylabel)
+ ax.set_xlabel('Altitude (m)')
+ ax.legend()
+ ax.grid()
+ ax.set_ylim([-8, 5])
+
+
+def plot_relative_difference_map(visualizer: StudyVisualizerForMeanValues):
+ label = ' mean annual maxima of {} ({})'.format('', '')
+ visualizer.plot_abstract_fast(massif_name_to_value=visualizer.massif_name_to_empirical_mean,
+ label='Empirical' + label)
+
+ print(visualizer.massif_name_to_parametric_mean)
+ visualizer.plot_abstract_fast(massif_name_to_value=visualizer.massif_name_to_parametric_mean,
+ label='Model' + label)
+ print(visualizer.massif_name_to_relative_difference_for_mean)
+ visualizer.plot_abstract_fast(massif_name_to_value=visualizer.massif_name_to_relative_difference_for_mean,
+ label='Relative difference of' + label)
+ return list(visualizer.massif_name_to_relative_difference_for_mean.values())
diff --git a/projects/ogorman/gorman_figures/figure1/study_visualizer_for_double_stationary_fit.py b/projects/ogorman/gorman_figures/figure1/study_visualizer_for_double_stationary_fit.py
index cd52f2fa7f1264a9e3f38c27663984a1e0c31f75..3fd5d2bd5cda20b4087f1f0d9dbe3d1aa55226d0 100644
--- a/projects/ogorman/gorman_figures/figure1/study_visualizer_for_double_stationary_fit.py
+++ b/projects/ogorman/gorman_figures/figure1/study_visualizer_for_double_stationary_fit.py
@@ -9,8 +9,7 @@ from extreme_fit.model.margin_model.utils import \
import matplotlib.pyplot as plt
-from projects.ogorman.gorman_figures.figure1 import \
- ResultFromDoubleStationaryFit
+from projects.ogorman.gorman_figures.figure1.result_from_stationary_fit import ResultFromDoubleStationaryFit
class StudyVisualizerForReturnLevelChange(StudyVisualizer):
@@ -79,6 +78,7 @@ class StudyVisualizerForReturnLevelChange(StudyVisualizer):
self.plot_abstract(
massif_name_to_value=result.massif_name_to_difference_return_level_and_maxima,
label=label, plot_name=plot_name,
+ fit_method=self.fit_method,
cmap=plt.cm.coolwarm)
def plot_shape_values(self):
@@ -110,44 +110,5 @@ class StudyVisualizerForReturnLevelChange(StudyVisualizer):
plot_name = 'Change {} in return levels'.format(b)
self.plot_abstract(
massif_name_to_value=massif_name_to_value,
- label=label, plot_name=plot_name)
-
- def plot_abstract(self, massif_name_to_value, label, plot_name, graduation=10.0, cmap=plt.cm.bwr):
- plot_name1 = '{}/{}'.format(self.study.altitude, plot_name)
- plot_name2 = '{}/{}'.format(plot_name.split()[0], plot_name)
- for plot_name in [plot_name1, plot_name2]:
- self.load_plot(cmap, graduation, label, massif_name_to_value)
- self.plot_name = plot_name
- self.show_or_save_to_file(add_classic_title=False, tight_layout=True, no_title=True,
- dpi=500)
- plt.close()
-
- def load_plot(self, cmap, graduation, label, massif_name_to_value):
- max_abs_change = max([abs(e) for e in massif_name_to_value.values()])
- ticks, labels = ticks_values_and_labels_for_percentages(graduation=graduation, max_abs_change=max_abs_change)
- min_ratio = -max_abs_change
- max_ratio = max_abs_change
- cmap = get_shifted_map(min_ratio, max_ratio, cmap)
- massif_name_to_color = {m: get_colors([v], cmap, min_ratio, max_ratio)[0]
- for m, v in massif_name_to_value.items()}
- ticks_values_and_labels = ticks, labels
- ax = plt.gca()
- AbstractStudy.visualize_study(ax=ax,
- massif_name_to_value=massif_name_to_value,
- 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)
- ax.set_title('Fit method is {}'.format(fitmethod_to_str(self.fit_method)))
+ label=label, plot_name=plot_name,
+ fit_method=self.fit_method)