[paper 2] add snow load in 1 day. add plot for contrasting trend curve for each climatic regions

experiment/meteo_france_data/scm_models_data/crocus/crocus.py

@@ -7,7 +7,7 @@ from experiment.meteo_france_data.scm_models_data.abstract_study import Abstract
 from experiment.meteo_france_data.scm_models_data.crocus.crocus_variables import CrocusTotalSweVariable, \
     CrocusDepthVariable, CrocusRecentSweVariableThreeDays, TotalSnowLoadVariable, RecentSnowLoadVariableThreeDays, \
     CrocusSnowLoadEurocodeVariable, CrocusDensityVariable, RecentSnowLoadVariableFiveDays, \
-    RecentSnowLoadVariableSevenDays
+    RecentSnowLoadVariableSevenDays, RecentSnowLoadVariableOneDay
 
 
 class Crocus(AbstractStudy):
@@ -20,6 +20,7 @@ class Crocus(AbstractStudy):
                                   RecentSnowLoadVariableThreeDays, TotalSnowLoadVariable,
                                   CrocusSnowLoadEurocodeVariable,
                                   CrocusDensityVariable,
+                                  RecentSnowLoadVariableOneDay,
                                   RecentSnowLoadVariableFiveDays,
                                   RecentSnowLoadVariableSevenDays]
         super().__init__(variable_class, *args, **kwargs)
@@ -60,6 +61,9 @@ class CrocusSnowLoadTotal(Crocus):
     def __init__(self, *args, **kwargs):
         Crocus.__init__(self, TotalSnowLoadVariable, *args, **kwargs)
 
+class CrocusSnowLoad1Day(CrocusSweTotal):
+    def __init__(self, *args, **kwargs):
+        Crocus.__init__(self, RecentSnowLoadVariableOneDay, *args, **kwargs)
 
 class CrocusSnowLoad3Days(CrocusSweTotal):
     def __init__(self, *args, **kwargs):

experiment/meteo_france_data/scm_models_data/crocus/crocus_variables.py

@@ -20,6 +20,14 @@ class CrocusTotalSweVariable(CrocusVariable):
         return 'WSN_T_ISBA'
 
 
+class CrocusRecentSweVariableOneDay(CrocusTotalSweVariable):
+    NAME = 'Snow Water Equivalent last 1 day'
+
+    @classmethod
+    def keyword(cls):
+        return 'SWE_1DY_ISBA'
+
+
 class CrocusRecentSweVariableThreeDays(CrocusTotalSweVariable):
     NAME = 'Snow Water Equivalent last 3 days'
 
@@ -52,6 +60,8 @@ class AbstractSnowLoadVariable(CrocusVariable):
         snow_pressure = self.snow_load_multiplication_factor * super().daily_time_serie_array
         return snow_pressure
 
+class RecentSnowLoadVariableOneDay(AbstractSnowLoadVariable, CrocusRecentSweVariableOneDay):
+    NAME = 'Snow load last 1 day'
 
 class RecentSnowLoadVariableThreeDays(AbstractSnowLoadVariable, CrocusRecentSweVariableThreeDays):
     NAME = 'Snow load last 3 days'

papers/contrasting_snow_loads/main_result.py

+from multiprocessing.pool import Pool
+
+import matplotlib as mpl
+
+from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal, CrocusSnowLoad3Days, \
+    CrocusSnowLoad5Days, CrocusSnowLoad7Days, CrocusSnowLoad1Day
+from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import \
+    ConfidenceIntervalMethodFromExtremes
+from papers.contrasting_snow_loads.plot_contrasting_trend_curves import plot_contrasting_trend_curves
+from papers.exceeding_snow_loads.paper_main_utils import load_altitude_to_visualizer
+from papers.exceeding_snow_loads.paper_utils import paper_study_classes, paper_altitudes
+from papers.exceeding_snow_loads.result_trends_and_return_levels.main_result_trends_and_return_levels import \
+    compute_minimized_aic
+from papers.exceeding_snow_loads.result_trends_and_return_levels.plot_selection_curves import plot_selection_curves
+from papers.exceeding_snow_loads.result_trends_and_return_levels.plot_trend_curves import plot_trend_curves, \
+    plot_trend_map
+from papers.exceeding_snow_loads.result_trends_and_return_levels.plot_uncertainty_curves import plot_uncertainty_massifs
+from papers.exceeding_snow_loads.result_trends_and_return_levels.plot_uncertainty_histogram import \
+    plot_uncertainty_histogram
+from root_utils import NB_CORES
+
+mpl.rcParams['text.usetex'] = True
+mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
+
+
+def intermediate_result(altitudes, massif_names=None,
+                        model_subsets_for_uncertainty=None, uncertainty_methods=None,
+                        study_class=CrocusSnowLoad3Days,
+                        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)
+    # 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
+    plot_contrasting_trend_curves(altitude_to_visualizer)
+
+
+def major_result():
+    uncertainty_methods = [ConfidenceIntervalMethodFromExtremes.my_bayes,
+                           ConfidenceIntervalMethodFromExtremes.ci_mle][1:]
+    massif_names = None
+    model_subsets_for_uncertainty = None
+    altitudes = paper_altitudes
+    altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700]
+    study_classes = [CrocusSnowLoad1Day, CrocusSnowLoad3Days, CrocusSnowLoad5Days, CrocusSnowLoad7Days][:]
+    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=[1500, 1800][:], massif_names=None,
+    #                     uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
+    #                                          ConfidenceIntervalMethodFromExtremes.ci_mle][1:],
+    #                     multiprocessing=True)
\ No newline at end of file

papers/contrasting_snow_loads/main_spatial_relative_change_in_maxima_at_fixed_altitude.py

@@ -9,11 +9,11 @@ from experiment.meteo_france_data.scm_models_data.visualization.study_visualizat
     StudyVisualizer
 import matplotlib.pyplot as plt
 
-from experiment.exceeding_snow_loads.discussion_data_comparison_with_eurocode.crocus_study_comparison_with_eurocode import \
+from papers.exceeding_snow_loads.discussion_data_comparison_with_eurocode.crocus_study_comparison_with_eurocode import \
     CrocusDifferenceSnowLoad, \
     CrocusSnowDensityAtMaxofSwe, CrocusDifferenceSnowLoadRescaledAndEurocodeToSeeSynchronization, \
     CrocusSnowDepthAtMaxofSwe, CrocusSnowDepthDifference
-from experiment.exceeding_snow_loads.paper_utils import dpi_paper1_figure
+from papers.exceeding_snow_loads.paper_utils import dpi_paper1_figure
 
 
 def test():

papers/contrasting_snow_loads/plot_contrasting_trend_curves.py

+from typing import Dict
+import matplotlib.pyplot as plt
+
+from experiment.meteo_france_data.scm_models_data.abstract_extended_study import AbstractExtendedStudy
+from experiment.meteo_france_data.scm_models_data.visualization.utils import create_adjusted_axes
+from papers.exceeding_snow_loads.paper_utils import dpi_paper1_figure
+from papers.exceeding_snow_loads.study_visualizer_for_non_stationary_trends import \
+    StudyVisualizerForNonStationaryTrends
+
+
+def plot_contrasting_trend_curves(altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends]):
+    """
+    Plot a single trend curves
+    :return:
+    """
+    visualizer = list(altitude_to_visualizer.values())[0]
+
+    ax = create_adjusted_axes(1, 1)
+    ax_twinx = ax.twinx()
+    ax_twiny = ax.twiny()
+
+    trend_summary_values = list(zip(*[v.trend_summary_contrasting_values() for v in altitude_to_visualizer.values()]))
+    altitudes, *mean_changes = trend_summary_values
+
+    # parameters
+    width = 150
+    size = 20
+    legend_fontsize = 20
+    color = 'white'
+    labelsize = 15
+    linewidth = 3
+    # ax.bar(altitudes, percent_decrease, width=width, color=color, edgecolor='blue', label='decreasing trend',
+    #        linewidth=linewidth)
+    # ax.bar(altitudes, percent_decrease_signi, width=width, color=color, edgecolor='black',
+    #        label='significative decreasing trend',
+    #        linewidth=linewidth)
+    # ax.legend(loc='upper left', prop={'size': size})
+
+    for ax_horizontal in [ax, ax_twiny]:
+        if ax_horizontal == ax_twiny:
+            ax_horizontal.plot(altitudes, [0 for _ in altitudes], linewidth=0)
+        else:
+            ax_horizontal.set_xlabel('Altitude', fontsize=legend_fontsize)
+        ax_horizontal.set_xticks(altitudes)
+        # ax_horizontal.set_xlim([700, 5000])
+        ax_horizontal.tick_params(labelsize=labelsize)
+
+    # Set the number of massifs on the upper axis
+    ax_twiny.set_xticklabels([v.study.nb_study_massif_names for v in altitude_to_visualizer.values()])
+    ax_twiny.set_xlabel('Total number of massifs at each altitude (for the percentage)', fontsize=legend_fontsize)
+
+    # ax.set_ylabel('Massifs with decreasing trend (\%)', fontsize=legend_fontsize)
+    # max_percent = int(max(percent_decrease))
+    # n = 2 + (max_percent // 10)
+    # ax_ticks = [10 * i for i in range(n)]
+    # # upper_lim = max_percent + 3
+    # upper_lim = n + 5
+    # ax_lim = [0, upper_lim]
+    # for axis in [ax, ax_twinx]:
+    #     axis.set_ylim(ax_lim)
+    #     axis.set_yticks(ax_ticks)
+    #     axis.tick_params(labelsize=labelsize)
+    ax.yaxis.grid()
+
+    label_curve = (visualizer.label).replace('change', 'decrease')
+    ax_twinx.set_ylabel(label_curve.replace('', ''), fontsize=legend_fontsize)
+    for region_name, mean_change in zip(AbstractExtendedStudy.region_names, mean_changes):
+        if len(mean_changes) > 1:
+            label = region_name
+        else:
+            label = 'Mean relative change'
+        ax_twinx.plot(altitudes, mean_change, label=label, linewidth=linewidth, marker='o')
+        ax_twinx.legend(loc='upper right', prop={'size': size})
+
+    # Save plot
+    visualizer.plot_name = 'Trend curves'
+    visualizer.show_or_save_to_file(no_title=True, dpi=dpi_paper1_figure)
+    plt.close()

papers/exceeding_snow_loads/study_visualizer_for_non_stationary_trends.py

@@ -429,3 +429,19 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
         #                      if massif_name_to_region_name[m] == region_name]
         #     mean_decreases.append(compute_mean_decrease(change_values))
         return (self.altitude, percentage_decrease, percentage_decrease_significative, *mean_decreases)
+
+    def trend_summary_contrasting_values(self):
+        # trend_tests = list(self.massif_name_to_trend_test_that_minimized_aic.values())
+        # decreasing_trend_tests = [t for t in trend_tests if t.time_derivative_of_return_level < 0]
+        # percentage_decrease = 100 * len(decreasing_trend_tests) / len(trend_tests)
+        # significative_decrease_trend_tests = [t for t in decreasing_trend_tests if t.is_significant]
+        # percentage_decrease_significative = 100 * len(significative_decrease_trend_tests) / len(trend_tests)
+        compute_mean_change = lambda l: np.mean(np.array(list(l)))
+        mean_changes = [compute_mean_change(self.massif_name_to_relative_change_value.values())]
+        # Compute mean relatives per regions (for the moment i don't add the region means)
+        massif_name_to_region_name = AbstractExtendedStudy.massif_name_to_region_name
+        for region_name in AbstractExtendedStudy.real_region_names:
+            change_values = [v for m, v in self.massif_name_to_relative_change_value.items()
+                             if massif_name_to_region_name[m] == region_name]
+            mean_changes.append(compute_mean_change(change_values))
+        return (self.altitude, *mean_changes)