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 2c4068f30cf5ecebe31974e40243c7f456951edd..c0fbb7fb988841d37c149b1c50c717c77daaff3a 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
@@ -81,7 +81,8 @@ 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))
+ new_tick = round(tick, 2)
+ positive_ticks.append(new_tick)
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]
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
index 82b712a59427463f9d3a824e220149bb75f56c24..8d67f062928f00030cf546dd91c5cc6d14fc0462 100644
--- 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
@@ -11,7 +11,7 @@ def plot_against_altitude(x_ticks, ax, massif_id, massif_name, values, altitude=
if massif_name_as_labels:
di = massif_id // 8
if di == 0:
- linestyle = '-'
+ linestyle = '-.'
elif di == 1:
linestyle = 'dotted'
else:
@@ -27,7 +27,7 @@ def plot_against_altitude(x_ticks, ax, massif_id, massif_name, values, altitude=
linestyle = None
label = '{} m'.format(altitude)
if not fill:
- ax.plot(x_ticks, values, color=color, linewidth=2, label=label, linestyle=linestyle)
+ ax.plot(x_ticks, values, color=color, linewidth=2, label=label, linestyle=linestyle, marker='o')
else:
lower_bound, upper_bound = zip(*values)
# ax.fill_between(altitudes, lower_bound, upper_bound, color=color, alpha=0.2, label=label + '95\% confidence interval')
diff --git a/extreme_fit/distribution/gev/main_fevd_mle_two_covariates.R b/extreme_fit/distribution/gev/main_fevd_mle_two_covariates.R
index d3579a69962ea918342791791f0cb3e680e4d655..18dcac06efb3dd2fbbb97d133daa93d6f687e1a2 100644
--- a/extreme_fit/distribution/gev/main_fevd_mle_two_covariates.R
+++ b/extreme_fit/distribution/gev/main_fevd_mle_two_covariates.R
@@ -56,7 +56,9 @@ spatial_cov = c(100.0, 100.0)
temporal_cov = c(500.0, 500.0)
#
spatial_cov = c(100.0, 100.0, 100.0, 100.0)
-temporal_cov = c(500.0, 500.0, 500.0, 500.0)
+spatial_cov = c(1000.0, 1000.0, 1000.0, 1000)
+# temporal_cov = c(500.0, 500.0, 500.0, 500.0)
+temporal_cov = c(300.0, 300.0, 300.0, 300.0)
# spatial_cov = c(200.0)
# temporal_cov = c(500.0)
@@ -64,11 +66,18 @@ temporal_cov = c(500.0, 500.0, 500.0, 500.0)
v = make.qcov(res, vals = list( mu1 = spatial_cov, mu2 = temporal_cov, sigma2 =spatial_cov, sigma1 = temporal_cov ))
# v = make.qcov(res, vals = list(mu1 = c(0.0)))
-res_ci = ci.fevd.mle_fixed(res, alpha = 0.05, type = c("return.level"),
+
+
+res_ci = ci.fevd.mle_fixed(res, alpha = 0.05, type = c("parameter"),
return.period = 50, method = "normal", xrange = NULL, nint = 20, verbose = FALSE,
tscale = FALSE, return.samples = FALSE, qcov=v)
print(res_ci)
+# res_ci = ci.fevd.mle_fixed(res, alpha = 0.05, type = c("return.level"),
+# return.period = 50, method = "normal", xrange = NULL, nint = 20, verbose = FALSE,
+# tscale = FALSE, return.samples = FALSE, qcov=v)
+# print(res_ci)
+
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 4f5cf39bed754d32e861d6b34963b448a6c9a5ef..df1542dc6de41e8d9aeb7e13b86e859494bc316b 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
@@ -52,13 +52,14 @@ class AbstractResultFromExtremes(AbstractResultFromModelFit):
r_matrix = self.name_to_value[name]
return pd.DataFrame(np.array(r_matrix), columns=r.colnames(r_matrix))
- def confidence_interval_method(self, quantile_level, alpha_interval, transformed_temporal_covariate, ci_method):
+ def confidence_interval_method(self, quantile_level, alpha_interval, transformed_temporal_covariate, ci_method,
+ return_level_interval=True):
return_period = round(1 / (1 - quantile_level))
common_kwargs = {
'return.period': return_period,
'alpha': alpha_interval,
'tscale': False,
- 'type': r.c("return.level")
+ 'type': r.c("return.level") if return_level_interval else r.c("parameter")
}
if self.param_name_to_dim:
if isinstance(transformed_temporal_covariate, (int, float, np.int, np.float, np.int64)):
diff --git a/projects/altitude_spatial_model/altitudes_fit/altitudes_studies.py b/projects/altitude_spatial_model/altitudes_fit/altitudes_studies.py
index b72a90d26253604733128c4d075c3c5315c88704..73dda0cc9fcec446f8c911b1d38a015c0d36480e 100644
--- a/projects/altitude_spatial_model/altitudes_fit/altitudes_studies.py
+++ b/projects/altitude_spatial_model/altitudes_fit/altitudes_studies.py
@@ -140,10 +140,10 @@ class AltitudesStudies(object):
ax.legend(handles[::-1], labels[::-1])
plot_name = 'Annual maxima of {} in {}'.format(SCM_STUDY_CLASS_TO_ABBREVIATION[self.study_class],
massif_name.replace('_', ' '))
- ax.set_ylabel('{} ({})'.format(plot_name, self.study.variable_unit), fontsize=15)
+ # ax.set_ylabel('{} ({})'.format(plot_name, self.study.variable_unit), fontsize=15)
+ # ax.set_xlabel('years', fontsize=15)
plot_name = 'time series/' + plot_name
- ax.set_xlabel('years', fontsize=15)
- self.show_or_save_to_file(plot_name=plot_name, show=show, no_title=True)
+ self.show_or_save_to_file(plot_name=plot_name, show=show, no_title=True, tight_layout=True)
ax.clear()
plt.close()
diff --git a/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py b/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
index d4b20888120afe71359a9ace1749de4908ac963c..4b0f672b1892cbeec3a638c13ceb8de0045d3fc5 100644
--- a/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
+++ b/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
@@ -5,7 +5,8 @@ import matplotlib as mpl
from projects.altitude_spatial_model.altitudes_fit.plots.plot_coherence_curves import plot_coherence_curves
from projects.altitude_spatial_model.altitudes_fit.plots.plot_histogram_altitude_studies import \
- plot_histogram_all_models_against_altitudes, plot_histogram_all_trends_against_altitudes
+ plot_histogram_all_models_against_altitudes, plot_histogram_all_trends_against_altitudes, \
+ plot_shoe_plot_changes_against_altitude
mpl.rcParams['text.usetex'] = True
mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
@@ -55,8 +56,9 @@ def main_loop(altitudes_list, massif_names, seasons, study_classes):
def plot_visualizers(massif_names, visualizer_list):
plot_histogram_all_trends_against_altitudes(massif_names, visualizer_list)
- # plot_histogram_all_models_against_altitudes(massif_names, visualizer_list)
- plot_coherence_curves(massif_names, visualizer_list)
+ plot_histogram_all_models_against_altitudes(massif_names, visualizer_list)
+ plot_shoe_plot_changes_against_altitude(massif_names, visualizer_list)
+ # plot_coherence_curves(massif_names, visualizer_list)
pass
@@ -68,7 +70,7 @@ def plot_visualizer(massif_names, visualizer):
# for change in [True, False, None]:
# studies.plot_mean_maxima_against_altitude(massif_names=massif_names, std=std, change=change)
# Plot the results for the model that minimizes the individual aic
- # plot_individual_aic(visualizer)
+ plot_individual_aic(visualizer)
# Plot the results for the model that minimizes the total aic
# plot_total_aic(model_classes, visualizer)
pass
diff --git a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitude_group.py b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitude_group.py
index 5bf1bc5cfa9c51c6aa766199a149380438663576..4d9a644f961613faa46238f28d9279e47f619292 100644
--- a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitude_group.py
+++ b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitude_group.py
@@ -19,21 +19,33 @@ class AbstractAltitudeGroup(object):
def reference_altitude(self):
raise NotImplementedError
+ @property
+ def group_id(self):
+ raise NotImplementedError
+
@property
def xlabel(self):
# warning: the label could not correspond to all massifs, some might have been fitted with less data
# idx = get_index_group_from_reference_altitude(reference_altitude)
# min_altitude, *_, max_altitude = altitudes_for_groups[idx]
i = self.reference_altitude // 1000
- min_altitude, max_altitude = 1000 * i, 1000 * (i + 1)
- return 'Altitude = {} m\n' \
- 'Estimated with maxima between {} m and {} m'.format(self.reference_altitude,
- min_altitude,
- max_altitude)
+ if self.group_id == 1:
+ formula = 'below 1000 m'
+ elif self.group_id == 4:
+ formula = 'above 3000 m'
+ else:
+ min_altitude, max_altitude = 1000 * i, 1000 * (i + 1)
+ formula = 'between {} m and {} m'.format(min_altitude, max_altitude)
+ return 'Altitude = {} m. Selected models were estimated with\n' \
+ 'altitude group {}, i.e. with maxima {}'.format(self.reference_altitude, self.group_id, formula)
class LowAltitudeGroup(AbstractAltitudeGroup):
+ @property
+ def group_id(self):
+ return 1
+
@property
def name(self):
return 'low'
@@ -45,6 +57,10 @@ class LowAltitudeGroup(AbstractAltitudeGroup):
class MidAltitudeGroup(AbstractAltitudeGroup):
+ @property
+ def group_id(self):
+ return 2
+
@property
def name(self):
return 'mid'
@@ -56,6 +72,10 @@ class MidAltitudeGroup(AbstractAltitudeGroup):
class HighAltitudeGroup(AbstractAltitudeGroup):
+ @property
+ def group_id(self):
+ return 3
+
@property
def name(self):
return 'high'
@@ -67,6 +87,10 @@ class HighAltitudeGroup(AbstractAltitudeGroup):
class VeyHighAltitudeGroup(AbstractAltitudeGroup):
+ @property
+ def group_id(self):
+ return 4
+
@property
def name(self):
return 'very high'
@@ -86,6 +110,7 @@ class DefaultAltitudeGroup(AbstractAltitudeGroup):
def reference_altitude(self):
return 500
+
def get_altitude_group_from_altitudes(altitudes):
s = set(altitudes)
if s == set(altitudes_for_groups[0]):
diff --git a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
index 5fdb9258e85c18f5c0641d398edde73e4ca29213..6ee246b24fae743a2b362f8a5fb1cd2f9526faa2 100644
--- a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
+++ b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
@@ -70,12 +70,12 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
self._method_name_and_order_to_max_abs = {}
self._max_abs_for_shape = None
- moment_names = ['moment', 'changes_for_moment', 'relative_changes_for_moment'][:]
+ moment_names = ['moment', 'changes_of_moment', 'relative_changes_of_moment'][:]
orders = [1, 2, None][2:]
def get_massif_altitudes(self, massif_name):
valid_altitudes = [altitude for altitude, study in self.studies.altitude_to_study.items()
- if massif_name in study.study_massif_names]
+ if massif_name in study.study_massif_names]
massif_altitudes = []
for altitude in valid_altitudes:
study = self.studies.altitude_to_study[altitude]
@@ -135,13 +135,11 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
# Plot settings
moment = ' '.join(method_name.split('_'))
moment = moment.replace('moment', '{} in 2019'.format(OneFoldFit.get_moment_str(order=order)))
- plot_name = '{}{} for annual maxima of {}'.format(OneFoldFit.folder_for_plots, moment,
- SCM_STUDY_CLASS_TO_ABBREVIATION[
- self.studies.study_class])
+ plot_name = '{}{} '.format(OneFoldFit.folder_for_plots, moment)
massif_name_to_text = self.massif_name_to_best_name
if 'change' in method_name:
- plot_name += '\nbetween {} and {}'.format(2019 - 50, 2019)
+ plot_name += ' between {} and {}'.format(2019 - 50, 2019)
if 'relative' not in method_name:
# Put the relative score as text on the plot for the change.
massif_name_to_text = {m: ('+' if v > 0 else '') + str(int(v)) + '\%' for m, v in
@@ -174,6 +172,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
@property
def add_colorbar(self):
+ # return isinstance(self.altitude_group, (VeyHighAltitudeGroup))
return isinstance(self.altitude_group, (VeyHighAltitudeGroup, MidAltitudeGroup))
def plot_against_years(self, method_name, order):
@@ -258,9 +257,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
def plot_shape_map(self):
- label = 'Shape parameter for {} maxima of {} in 2019'.format(self.study.season_name,
- SCM_STUDY_CLASS_TO_ABBREVIATION[
- type(self.study)])
+ label = 'Shape parameter in 2019 (no unit)'
self.plot_map(massif_name_to_value=self.massif_name_to_shape,
label=label,
plot_name=label,
@@ -395,11 +392,8 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
nb_valid_massif_names = len(valid_massif_names)
nbs = np.zeros(4)
- changes = []
for one_fold in [one_fold for m, one_fold in self.massif_name_to_one_fold_fit.items()
- if m in valid_massif_names]:
- # Compute changes
- changes.append(one_fold.change_in_return_level_for_reference_altitude)
+ if m in valid_massif_names]:
# Compute nb of non stationary models
if one_fold.change_in_return_level_for_reference_altitude == 0:
continue
@@ -410,9 +404,18 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
nbs[idx + 1] += 1
percents = 100 * nbs / nb_valid_massif_names
- mean_change = np.mean(changes)
- median_change = np.median(changes)
- return [nb_valid_massif_names, mean_change, median_change] + list(percents)
+ return [nb_valid_massif_names] + list(percents)
+
+ def all_changes(self, massif_names):
+ """return percents which contain decrease, significant decrease, increase, significant increase percentages"""
+ valid_massif_names = self.get_valid_names(massif_names)
+ changes = []
+ for one_fold in [one_fold for m, one_fold in self.massif_name_to_one_fold_fit.items()
+ if m in valid_massif_names]:
+ # Compute changes
+ changes.append(one_fold.change_in_return_level_for_reference_altitude)
+
+ return changes
def get_valid_names(self, massif_names):
valid_massif_names = set(self.massif_name_to_one_fold_fit.keys())
diff --git a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py
index a0db58c0959f55ef88405629194c288948e34ee6..0f4a644ae0730cae9c57a6a871cf3c980a85ac24 100644
--- a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py
+++ b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py
@@ -94,13 +94,13 @@ class OneFoldFit(object):
@property
def change_in_return_level_for_reference_altitude(self) -> float:
- return self.changes_for_moment(altitudes=[self.altitude_plot], order=None)[0]
+ return self.changes_of_moment(altitudes=[self.altitude_plot], order=None)[0]
@property
def relative_change_in_return_level_for_reference_altitude(self) -> float:
- return self.relative_changes_for_moment(altitudes=[self.altitude_plot], order=None)[0]
+ return self.relative_changes_of_moment(altitudes=[self.altitude_plot], order=None)[0]
- def changes_for_moment(self, altitudes, year=2019, nb_years=50, order=1):
+ def changes_of_moment(self, altitudes, year=2019, nb_years=50, order=1):
changes = []
for altitude in altitudes:
mean_after = self.get_moment(altitude, year, order)
@@ -109,7 +109,7 @@ class OneFoldFit(object):
changes.append(change)
return changes
- def relative_changes_for_moment(self, altitudes, year=2019, nb_years=50, order=1):
+ def relative_changes_of_moment(self, altitudes, year=2019, nb_years=50, order=1):
relative_changes = []
for altitude in altitudes:
mean_after = self.get_moment(altitude, year, order)
diff --git a/projects/altitude_spatial_model/altitudes_fit/plots/plot_coherence_curves.py b/projects/altitude_spatial_model/altitudes_fit/plots/plot_coherence_curves.py
index 71ad064db4f7fa3d5192f6a58f74a3bf31a5b88a..858edaeb47ac94888f596743642a04bcd6f3a7fb 100644
--- a/projects/altitude_spatial_model/altitudes_fit/plots/plot_coherence_curves.py
+++ b/projects/altitude_spatial_model/altitudes_fit/plots/plot_coherence_curves.py
@@ -16,6 +16,10 @@ def plot_coherence_curves(massif_names, visualizer_list: List[
for massif_name in all_valid_names:
_, axes = plt.subplots(2, 2)
axes = axes.flatten()
+ for i, ax in enumerate(axes):
+ if i % 2 == 1:
+ ax.set_yticks([])
+ axes = [ax if i % 2 == 0 else ax.twinx() for i, ax in enumerate(axes)]
colors = ['blue', 'yellow', 'green']
labels = ['Altitudinal-temporal model in 2019', 'Altitudinal-temporal model in 1969', 'Stationary model']
altitudinal_model = [True, True, False]
@@ -29,7 +33,8 @@ def plot_coherence_curves(massif_names, visualizer_list: List[
def plot_coherence_curve(axes, massif_name, visualizer_list: List[AltitudesStudiesVisualizerForNonStationaryModels],
is_altitudinal, color, global_label, year):
- x_all_list, values_all_list, labels, all_bound_list = load_all_list(massif_name, visualizer_list, is_altitudinal, year)
+ x_all_list, values_all_list, labels, all_bound_list = load_all_list(massif_name, visualizer_list, is_altitudinal,
+ year)
for i, label in enumerate(labels):
ax = axes[i]
# Plot with complete line
@@ -37,9 +42,11 @@ def plot_coherence_curve(axes, massif_name, visualizer_list: List[AltitudesStudi
value_list_i = value_list[i]
label_plot = global_label if j == 0 else None
if is_altitudinal:
- ax.plot(x_list, value_list_i, linestyle='solid', color=color, label=label_plot)
+ ax.plot(x_list, value_list_i, linestyle='solid', color=color)
+ # ax.plot(x_list, value_list_i, linestyle='solid', color=color, label=label_plot)
else:
- ax.plot(x_list, value_list_i, linestyle='None', color=color, label=label_plot, marker='x')
+ # ax.plot(x_list, value_list_i, linestyle='None', color=color, label=label_plot, marker='o')
+ ax.plot(x_list, value_list_i, linestyle='None', color=color, marker='o')
ax.plot(x_list, value_list_i, linestyle='dotted', color=color)
# Plot with dotted line
@@ -52,16 +59,21 @@ def plot_coherence_curve(axes, massif_name, visualizer_list: List[AltitudesStudi
# Plot confidence interval
if i == 3:
- for x_list, bounds in zip(x_all_list, all_bound_list):
+ for j, (x_list, bounds) in enumerate(list(zip(x_all_list, all_bound_list))):
if len(bounds) > 0:
lower_bound, upper_bound = bounds
+ # model_name = 'altitudinal-temporal' if is_altitudinal else 'stationary'
+ # fill_label = "95\% confidence interval for the {} model".format(model_name) if j == 0 else None
+ # ax.fill_between(x_list, lower_bound, upper_bound, color=color, alpha=0.2, label=fill_label)
ax.fill_between(x_list, lower_bound, upper_bound, color=color, alpha=0.2)
+ # if is_altitudinal:
+ # min, max = ax.get_ylim()
+ # ax.set_ylim([min, 1.5 * max])
+ # ax.legend(prop={'size': 5})
ax.set_ylabel(label)
-
- ax.legend(prop={'size': 5})
- if i >= 2:
- ax.set_xlabel('Altitude')
+ # if i >= 2:
+ # ax.set_xlabel('Altitude')
def load_all_list(massif_name, visualizer_list, altitudinal_model=True, year=2019):
@@ -76,15 +88,20 @@ def load_all_list(massif_name, visualizer_list, altitudinal_model=True, year=201
one_fold_fit = visualizer.massif_name_to_one_fold_fit[massif_name]
altitudes_list = list(range(min_altitude, max_altitude, 10))
gev_params_list = [one_fold_fit.get_gev_params(altitude, year) for altitude in altitudes_list]
- confidence_interval_values = [one_fold_fit.best_confidence_interval(altitude, year) for altitude in altitudes_list]
+ confidence_interval_values = [one_fold_fit.best_confidence_interval(altitude, year) for altitude in
+ altitudes_list]
else:
assert OneFoldFit.return_period == 100, 'change the call below'
altitudes_list, study_list_valid = zip(*[(a, s) for a, s in visualizer.studies.altitude_to_study.items()
- if massif_name in s.massif_name_to_stationary_gev_params_and_confidence_for_return_level_100[0]])
- gev_params_list = [study.massif_name_to_stationary_gev_params_and_confidence_for_return_level_100[0][massif_name]
- for study in study_list_valid]
- confidence_interval_values = [study.massif_name_to_stationary_gev_params_and_confidence_for_return_level_100[1][massif_name]
- for study in study_list_valid]
+ if massif_name in
+ s.massif_name_to_stationary_gev_params_and_confidence_for_return_level_100[
+ 0]])
+ gev_params_list = [
+ study.massif_name_to_stationary_gev_params_and_confidence_for_return_level_100[0][massif_name]
+ for study in study_list_valid]
+ confidence_interval_values = [
+ study.massif_name_to_stationary_gev_params_and_confidence_for_return_level_100[1][massif_name]
+ for study in study_list_valid]
# Checks
values = [(gev_params.location, gev_params.scale, gev_params.shape,
@@ -98,7 +115,13 @@ def load_all_list(massif_name, visualizer_list, altitudinal_model=True, year=201
all_values_list.append(values_list)
all_altitudes_list.append(altitudes_list)
all_bound_list.append(list(zip(*bound_list)))
- labels = ['location parameter', 'scale parameter', 'shape pameter', '{}-year return level'.format(OneFoldFit.return_period)]
- return all_altitudes_list, all_values_list, labels, all_bound_list
+ labels = ['location', 'scale', 'shape', '{}-year return levels'.format(OneFoldFit.return_period)]
+ for i, label in enumerate(labels):
+ if i < 3:
+ label += ' parameter'
+ unit = 'no unit' if i == 2 else visualizer_list[0].study.variable_unit
+ label += ' ({})'.format(unit)
+ labels[i] = label
+ return all_altitudes_list, all_values_list, labels, all_bound_list
diff --git a/projects/altitude_spatial_model/altitudes_fit/plots/plot_histogram_altitude_studies.py b/projects/altitude_spatial_model/altitudes_fit/plots/plot_histogram_altitude_studies.py
index b524e119e3599a9e24be614b773e966048e6cba4..58a56bb36fd6b6e90fc8f3b58bf91aefcec87599 100644
--- a/projects/altitude_spatial_model/altitudes_fit/plots/plot_histogram_altitude_studies.py
+++ b/projects/altitude_spatial_model/altitudes_fit/plots/plot_histogram_altitude_studies.py
@@ -21,10 +21,16 @@ def plot_histogram_all_models_against_altitudes(massif_names, visualizer_list: L
model_name_to_significant_percentages_for_v = v.model_name_to_percentages(massif_names, only_significant=True)
for model_name in model_names:
model_name_to_percentages[model_name].append(model_name_to_percentages_for_v[model_name])
- model_name_to_percentages_significant[model_name].append(model_name_to_significant_percentages_for_v[model_name])
+ model_name_to_percentages_significant[model_name].append(
+ model_name_to_significant_percentages_for_v[model_name])
# Sort model based on their mean percentage.
model_name_to_mean_percentage = {m: np.mean(a) for m, a in model_name_to_percentages.items()}
+ model_name_to_mean_percentage_significant = {m: np.mean(a) for m, a in
+ model_name_to_percentages_significant.items()}
sorted_model_names = sorted(model_names, key=lambda m: model_name_to_mean_percentage[m], reverse=True)
+ for model_name in sorted_model_names:
+ print(model_name_to_mean_percentage[model_name], model_name_to_mean_percentage_significant[model_name],
+ model_name)
# Plot part
ax = plt.gca()
@@ -40,9 +46,11 @@ def plot_histogram_all_models_against_altitudes(massif_names, visualizer_list: L
percentages = model_name_to_percentages[model_name]
percentages_significant = model_name_to_percentages_significant[model_name]
colors = ['white', 'yellow', 'orange', 'red']
- labels = ['{} m - {} m (\% out of {} massifs)'.format(1000 * i, 1000 * (i+1),
- len(v.get_valid_names(massif_names))) for i, v in enumerate(visualizer_list)]
- for x, color, percentage, label,percentage_significant in zip(x_shifted, colors, percentages, labels, percentages_significant):
+ labels = ['{} m - {} m (\% out of {} massifs)'.format(1000 * i, 1000 * (i + 1),
+ len(v.get_valid_names(massif_names))) for i, v in
+ enumerate(visualizer_list)]
+ for x, color, percentage, label, percentage_significant in zip(x_shifted, colors, percentages, labels,
+ percentages_significant):
ax.bar([x], [percentage], width=width, label=label,
linewidth=2 * linewidth, edgecolor='black', color=color)
heights = list(range(0, math.ceil(percentage_significant), 1))[::-1]
@@ -69,7 +77,7 @@ def plot_histogram_all_trends_against_altitudes(massif_names, visualizer_list: L
visualizer = visualizer_list[0]
all_trends = [v.all_trends(massif_names) for v in visualizer_list]
- nb_massifs, mean_changes, median_changes, *all_l = zip(*all_trends)
+ nb_massifs, *all_l = zip(*all_trends)
plt.close()
ax = plt.gca()
@@ -99,22 +107,48 @@ def plot_histogram_all_trends_against_altitudes(massif_names, visualizer_list: L
plot_nb_massif_on_upper_axis(ax, labelsize, legend_fontsize, nb_massifs, x)
- # PLot mean relative change against altitude
- ax_twinx = ax.twinx()
- ax_twinx.plot(x, mean_changes, label='Mean change', color='black')
- # ax_twinx.plot(x, median_changes, label='Median change', color='grey')
- ax_twinx.legend(loc='upper right', prop={'size': size})
- ylabel = 'Change of {}-year return levels ({})'.format(OneFoldFit.return_period, visualizer.study.variable_unit)
- ax_twinx.set_ylabel(ylabel, fontsize=legend_fontsize)
+ visualizer.plot_name = 'All trends'
+ visualizer.show_or_save_to_file(add_classic_title=False, no_title=True)
+ plt.close()
- low, up = ax_twinx.get_ylim()
- low2, up2 = ax.get_ylim()
- new_lim = min(low, low2), max(up, up2)
- ax.set_ylim(new_lim)
- ax_twinx.set_ylim(new_lim)
- visualizer.plot_name = 'All trends'
+def plot_shoe_plot_changes_against_altitude(massif_names, visualizer_list: List[
+ AltitudesStudiesVisualizerForNonStationaryModels]):
+ visualizer = visualizer_list[0]
+
+ all_changes = [v.all_changes(massif_names) for v in visualizer_list]
+ nb_massifs = [len(v.get_valid_names(massif_names)) for v in visualizer_list]
+
+ plt.close()
+ ax = plt.gca()
+ width = 5
+ size = 8
+ legend_fontsize = 10
+ labelsize = 10
+ linewidth = 3
+
+ x = np.array([3 * width * (i + 1) for i in range(len(nb_massifs))])
+
+ ax.boxplot(all_changes, positions=x, widths=width)
+
+ ax.legend(loc='upper left', prop={'size': size})
+ ax.set_ylabel('Changes between 1969 and 2019 in {}-year return levels ({})'.format(OneFoldFit.return_period,
+ visualizer.study.variable_unit),
+ fontsize=legend_fontsize)
+ ax.set_xlabel('Altitudes', fontsize=legend_fontsize)
+ ax.tick_params(axis='both', which='major', labelsize=labelsize)
+ ax.set_xticks(x)
+ ax.yaxis.grid()
+ altitudes = [v.altitude_group.reference_altitude for v in visualizer_list]
+ ax.set_xticklabels([str(a) for a in altitudes])
+
+ shift = 300
+ ax.set_xlim((min(altitudes) - shift, max(altitudes) + shift))
+
+ plot_nb_massif_on_upper_axis(ax, labelsize, legend_fontsize, nb_massifs, x)
+
+ visualizer.plot_name = 'All changes'
visualizer.show_or_save_to_file(add_classic_title=False, no_title=True)
plt.close()
diff --git a/projects/altitude_spatial_model/preliminary_analysis.py b/projects/altitude_spatial_model/preliminary_analysis.py
index bef0bf888382ba46a51e1652b9ee13184c5716c4..581cabdb9d31f08dd75728d9288e9e44abe51638 100644
--- a/projects/altitude_spatial_model/preliminary_analysis.py
+++ b/projects/altitude_spatial_model/preliminary_analysis.py
@@ -22,37 +22,89 @@ class PointwiseGevStudyVisualizer(AltitudesStudies):
self.altitudes_for_temporal_hypothesis = [600, 1500, 2400, 3300]
def plot_gev_params_against_altitude(self):
- for param_name in GevParams.PARAM_NAMES[:]:
+ for j, param_name in enumerate(GevParams.PARAM_NAMES + [100]):
ax = plt.gca()
+
massif_name_to_linear_coef = {}
massif_name_to_r2_score = {}
- for massif_name in self.study.all_massif_names()[:]:
- linear_coef, _, r2 = self._plot_gev_params_against_altitude_one_massif(ax, massif_name, param_name)
- massif_name_to_linear_coef[massif_name] = 1000 * linear_coef[0]
- massif_name_to_r2_score[massif_name] = str(round(r2, 2))
- print(massif_name_to_linear_coef, massif_name_to_r2_score)
- # Plot change against altitude
- # ax.legend(prop={'size': 7}, ncol=3)
- ax.set_xlabel('Altitude')
- ax.set_ylabel(GevParams.full_name_from_param_name(param_name) + ' parameter for a stationary GEV distribution')
+ massif_names = self.study.all_massif_names()[:]
+ for i in range(8):
+ for massif_name in massif_names[i::8]:
+ linear_coef, _, r2 = self._plot_gev_params_against_altitude_one_massif(ax, massif_name, param_name)
+ massif_name_to_linear_coef[massif_name] = 100 * linear_coef[0]
+ massif_name_to_r2_score[massif_name] = str(round(r2, 2))
+
+ # Display x label
+ xticks = [1000 * i for i in range(1, 5)]
+ ax.set_xticks(xticks)
+ fontsize_label = 15
+ ax.tick_params(labelsize=fontsize_label)
+
+ # ax.set_xlabel('Altitude')
+
+ # Compute the y label
+ if param_name in GevParams.PARAM_NAMES:
+ ylabel = GevParams.full_name_from_param_name(param_name) + ' parameter'
+ else:
+ ylabel = '{}-year return levels'.format(param_name)
+ # Add units
+ if param_name == GevParams.SHAPE:
+ unit = 'no unit'
+ else:
+ unit = self.study.variable_unit
+ ylabel += ' ({})'.format(unit)
+
+ # Display the y label on the twin axis
+ if param_name in [100, GevParams.SCALE]:
+ ax2 = ax.twinx()
+ ax2.set_yticks(ax.get_yticks())
+ ax2.set_ylim(ax.get_ylim())
+ ax2.set_ylabel(ylabel, fontsize=fontsize_label)
+ ax2.tick_params(labelsize=fontsize_label)
+ ax.set_yticks([])
+ tight_layout = False
+ else:
+ ax.tick_params(labelsize=fontsize_label)
+ tight_layout = True
+ ax.set_ylabel(ylabel, fontsize=fontsize_label)
+ # Make room for the ylabel
+ if param_name == 100:
+ plt.gcf().subplots_adjust(right=0.85)
+
plot_name = '{} change with altitude'.format(param_name)
- self.show_or_save_to_file(plot_name, no_title=True, tight_layout=True, show=False)
+
+ # # Display the legend
+ # ax.legend(labelspacing=2.5, ncol=8, handlelength=12, markerscale=0.7, bbox_to_anchor=(1.05, 1), loc='upper left',
+ # prop={'size': 2}, fontsize='x-large')
+ # plt.gcf().subplots_adjust(right=0.15)
+ # ax.set_yticks([])
+ # ax.set_ylabel('')
+
+ # plt.show()
+ self.show_or_save_to_file(plot_name, no_title=True, tight_layout=tight_layout, show=False)
ax.clear()
+ plt.close()
+
# Plot map of slope for each massif
visualizer = StudyVisualizer(study=self.study, show=False, save_to_file=True)
- ylabel = 'Linear slope for the {} parameter (change every 1000 m of altitude)'.format(param_name)
+ idx = 8 if param_name == GevParams.SHAPE else 1
+ label = 'Altitude gradient for the {}'.format(ylabel[:-idx] + '/100m)')
gev_param_name_to_graduation = {
- GevParams.LOC: 5,
- GevParams.SCALE: 1,
- GevParams.SHAPE: 0.1,
+ GevParams.LOC: 0.5,
+ GevParams.SCALE: 0.1,
+ GevParams.SHAPE: 0.01,
+ 100: 1,
}
+ if param_name == GevParams.SHAPE:
+ print(massif_name_to_linear_coef)
visualizer.plot_map(cmap=plt.cm.coolwarm, fit_method=self.study.fit_method,
graduation=gev_param_name_to_graduation[param_name],
- label=ylabel, massif_name_to_value=massif_name_to_linear_coef,
- plot_name=ylabel, add_x_label=False,
+ label=label, massif_name_to_value=massif_name_to_linear_coef,
+ plot_name=label.replace('/', ' every '), add_x_label=False,
negative_and_positive_values=param_name == GevParams.SHAPE,
add_colorbar=True,
massif_name_to_text=massif_name_to_r2_score,
+
)
plt.close()
@@ -64,7 +116,12 @@ class PointwiseGevStudyVisualizer(AltitudesStudies):
if massif_name in study.massif_name_to_stationary_gev_params:
gev_params = study.massif_name_to_stationary_gev_params[massif_name]
altitudes.append(altitude)
- params.append(gev_params.to_dict()[param_name])
+ if param_name in GevParams.PARAM_NAMES:
+ param = gev_params.to_dict()[param_name]
+ else:
+ assert isinstance(param_name, int)
+ param = gev_params.return_level(return_period=param_name)
+ params.append(param)
# confidence_intervals.append(gev_params.param_name_to_confidence_interval[param_name])
massif_id = self.study.all_massif_names().index(massif_name)
plot_against_altitude(altitudes, ax, massif_id, massif_name, params, fill=False)
@@ -195,7 +252,7 @@ if __name__ == '__main__':
altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300]
altitudes = [600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900]
# altitudes = paper_altitudes
- # altitudes = [1500, 1800]
+ altitudes = [1800, 2100]
visualizer = PointwiseGevStudyVisualizer(SafranSnowfall1Day, altitudes=altitudes)
visualizer.plot_gev_params_against_altitude()
# visualizer.plot_gev_params_against_time_for_all_altitudes()