[contrasted] add massif_name_to_annual_total to abstract_study. add plot that...

[contrasted] add massif_name_to_annual_total to abstract_study. add plot that compare side by side the total snowfall & maxima of snowfall.

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)
 

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)
 

projects/altitude_spatial_model/altitudes_fit/plots/compute_histogram_change_in_total_snowfall.py

+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)
+

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()