diff --git a/extreme_data/meteo_france_data/scm_models_data/abstract_study.py b/extreme_data/meteo_france_data/scm_models_data/abstract_study.py
index 083f3c9eb6b7e8f7816b177c9473894d6395156c..8390c92d6ffa012cc7fe99f91faa0242a4496917 100644
--- a/extreme_data/meteo_france_data/scm_models_data/abstract_study.py
+++ b/extreme_data/meteo_france_data/scm_models_data/abstract_study.py
@@ -335,6 +335,15 @@ class AbstractStudy(object):
year_to_annual_mean[year] = self.apply_annual_aggregation(time_serie)
return year_to_annual_mean
+ @cached_property
+ def massif_name_to_annual_total(self):
+ # Map each massif to an array of size nb_years
+ massif_name_to_annual_total = OrderedDict()
+ for i, massif_name in enumerate(self.study_massif_names):
+ maxima = np.array([self.year_to_annual_total[year][i] for year in self.ordered_years])
+ massif_name_to_annual_total[massif_name] = maxima
+ return massif_name_to_annual_total
+
def apply_annual_aggregation(self, time_serie):
return self.annual_aggregation_function(time_serie, axis=0)
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 5bc947e3580c62bbf2efdf66fd64c44f7314826f..92bcc6d858c62ad1781d572398812244af2cdf7d 100644
--- a/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
+++ b/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
@@ -7,7 +7,7 @@ from extreme_data.meteo_france_data.scm_models_data.crocus.crocus import CrocusS
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_shoe_plot_changes_against_altitude
+ plot_shoe_plot_changes_against_altitude, plot_shoe_plot_changes_against_altitude_for_maxima_and_total
mpl.rcParams['text.usetex'] = True
mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
@@ -32,7 +32,7 @@ def main():
altitudes_list = altitudes_for_groups[1:2]
elif fast:
massif_names = ['Vanoise', 'Haute-Maurienne', 'Vercors'][:]
- altitudes_list = altitudes_for_groups[3:]
+ altitudes_list = altitudes_for_groups[2:]
else:
massif_names = None
altitudes_list = altitudes_for_groups[:]
@@ -55,21 +55,19 @@ 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_trends_against_altitudes(massif_names, visualizer_list)
for relative in [True, False]:
- plot_shoe_plot_changes_against_altitude(massif_names, visualizer_list, relative=relative)
+ # plot_shoe_plot_changes_against_altitude(massif_names, visualizer_list, relative=relative)
+ plot_shoe_plot_changes_against_altitude_for_maxima_and_total(massif_names, visualizer_list, relative=relative)
# plot_coherence_curves(massif_names, visualizer_list)
- plot_coherence_curves(['Vanoise'], visualizer_list)
+ # plot_coherence_curves(['Vanoise'], visualizer_list)
def plot_visualizer(massif_names, visualizer):
# Plot time series
# visualizer.studies.plot_maxima_time_series(massif_names)
visualizer.studies.plot_maxima_time_series(['Vanoise'])
- # Plot moments against altitude
- # for std in [True, False][:]:
- # 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)
diff --git a/projects/altitude_spatial_model/altitudes_fit/plots/compute_histogram_change_in_total_snowfall.py b/projects/altitude_spatial_model/altitudes_fit/plots/compute_histogram_change_in_total_snowfall.py
new file mode 100644
index 0000000000000000000000000000000000000000..65e3e902c7a3bd22528c70dd38b101ffa541cadf
--- /dev/null
+++ b/projects/altitude_spatial_model/altitudes_fit/plots/compute_histogram_change_in_total_snowfall.py
@@ -0,0 +1,66 @@
+from typing import List
+import matplotlib.pyplot as plt
+import numpy as np
+
+from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall
+from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitudes_studies_visualizer_for_non_stationary_models import \
+ AltitudesStudiesVisualizerForNonStationaryModels
+from projects.contrasting_trends_in_snow_loads.article2_snowfall_versus_time_and_altitude.snowfall_plot import \
+ fit_linear_regression
+
+
+def compute_changes_in_total_snowfall(visualizer_list: List[
+ AltitudesStudiesVisualizerForNonStationaryModels], relative=True):
+ changes_in_total_snowfall = []
+ for visualizer in visualizer_list:
+ changes_for_the_visualizer = []
+ for massif_name in visualizer.massif_name_to_one_fold_fit.keys():
+ changes_massif = []
+ for altitude, study in visualizer.studies.altitude_to_study.items():
+ print(altitude)
+ if massif_name in study.study_massif_names:
+ change = compute_change_in_total(study, massif_name, relative, plot=False)
+ changes_massif.append(change)
+ print(len(changes_massif), 'length')
+ mean_change = np.mean(changes_massif)
+ changes_for_the_visualizer.append(mean_change)
+
+
+ changes_in_total_snowfall.append(changes_for_the_visualizer)
+ return changes_in_total_snowfall
+
+def compute_change_in_total(study, massif_name, relative, plot=False):
+ annual_total = study.massif_name_to_annual_total[massif_name]
+ a, b, r2score = fit_linear_regression(study.ordered_years, annual_total)
+ years_for_change = [1969, 2019]
+ values = [a * y + b for y in years_for_change]
+ change = values[1] - values[0]
+ if relative:
+ change *= 100 / values[0]
+ if plot:
+ ax = plt.gca()
+ ax.plot(study.ordered_years, annual_total)
+ linear_plot = [a * y + b for y in study.ordered_years]
+ ax.plot(study.ordered_years, linear_plot)
+ ax.plot(years_for_change, values, linewidth=0, marker='o')
+ ax.set_xlim((study.year_min, study.year_max))
+ ax.set_ylabel('Total of snowfall for the {} at {} ({})'.format(massif_name, study.altitude, study.variable_unit))
+ plt.show()
+ ax.clear()
+ plt.close()
+ return change
+
+
+if __name__ == '__main__':
+ altitude = 3000
+ year_min = 1959
+ year_max = 2019
+ study = SafranSnowfall(altitude=altitude, year_min=year_min, year_max=year_max)
+ print(study.study_massif_names)
+ # print(study.massif_name_to_annual_maxima)
+ # print(study.year_to_daily_time_serie_array[1959].shape)
+ # print(study.massif_name_to_daily_time_series['Vanoise'].shape)
+ # study._save_excel_with_longitutde_and_latitude()
+ print(study.massif_name_to_annual_total['Vanoise'])
+ compute_change_in_total(study, 'Vanoise', relative=True, plot=True)
+
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 473ebb3441ee76c751d0de3a0ea34345837439d8..1e839ae8f23be5486d17f6b53dba2dc06755f8b4 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
@@ -9,6 +9,8 @@ from matplotlib.lines import Line2D
from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitudes_studies_visualizer_for_non_stationary_models import \
AltitudesStudiesVisualizerForNonStationaryModels
from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.one_fold_fit import OneFoldFit
+from projects.altitude_spatial_model.altitudes_fit.plots.compute_histogram_change_in_total_snowfall import \
+ compute_changes_in_total_snowfall
def plot_histogram_all_models_against_altitudes(massif_names, visualizer_list: List[
@@ -171,6 +173,66 @@ def plot_shoe_plot_changes_against_altitude(massif_names, visualizer_list: List[
plt.close()
+def plot_shoe_plot_changes_against_altitude_for_maxima_and_total(massif_names, visualizer_list: List[
+ AltitudesStudiesVisualizerForNonStationaryModels],
+ relative=False):
+ visualizer = visualizer_list[0]
+
+ all_changes_total= compute_changes_in_total_snowfall(visualizer_list,
+ relative)
+
+ all_changes_extreme = [v.all_changes(massif_names, relative=relative) for v in visualizer_list]
+ all_changes_extreme = list(zip(*all_changes_extreme))[:1][0]
+
+ all_changes = [all_changes_extreme, all_changes_total]
+ labels = ['{}-year return levels'.format(OneFoldFit.return_period), 'Total snowfall']
+ colors = ['darkgreen', 'royalblue']
+ 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 = 15
+ labelsize = 10
+ linewidth = 3
+
+ x = np.array([3 * width * (i + 1) for i in range(len(nb_massifs))])
+ for j, (changes, label, color) in enumerate(list(zip(all_changes, labels, colors))):
+ print(len(changes), changes, label)
+ positions = x + (2 * j - 1) * 0.5 * width
+ bplot = ax.boxplot(changes, positions=positions, widths=width, patch_artist=True, showmeans=True)
+ for patch in bplot['boxes']:
+ patch.set_facecolor(color)
+
+ custom_lines = [Line2D([0], [0], color=color, lw=4) for color in colors]
+ loc = 'lower right' if relative else 'upper left'
+ ax.legend(custom_lines, labels, loc=loc)
+
+ start = 'Relative changes' if relative else 'Changes'
+ unit = '\%' if relative else visualizer.study.variable_unit
+ ax.set_ylabel('{} between 1969 and 2019 ({})'.format(start, unit),
+ fontsize=legend_fontsize)
+ ax.set_xlabel('Elevation', 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 = 2 * width
+ ax.set_xlim((min(x) - shift, max(x) + shift))
+
+ # I could display the number of massif used to build each box plot.
+ # plot_nb_massif_on_upper_axis(ax, labelsize, legend_fontsize, nb_massifs, x)
+
+ visualizer.plot_name = 'All ' + start + ' and total '
+ visualizer.show_or_save_to_file(add_classic_title=False, no_title=True)
+
+ plt.close()
+
+
def plot_nb_massif_on_upper_axis(ax, labelsize, legend_fontsize, nb_massifs, x):
# Plot number of massifs on the upper axis
ax_twiny = ax.twiny()