[contrasting] improve snowfall plot

extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py

@@ -710,13 +710,14 @@ class StudyVisualizer(VisualizationParameters):
 
     def plot_map(self, cmap, fit_method, graduation, label, massif_name_to_value, plot_name, add_x_label=True,
                  negative_and_positive_values=True, massif_name_to_text=None):
-        load_plot(cmap, graduation, label, massif_name_to_value, self.study.altitude, fitmethod_to_str(fit_method),
-                  add_x_label=add_x_label, negative_and_positive_values=negative_and_positive_values,
-                  massif_name_to_text=massif_name_to_text)
-        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()
+        if len(massif_name_to_value) > 0:
+            load_plot(cmap, graduation, label, massif_name_to_value, self.study.altitude, fitmethod_to_str(fit_method),
+                      add_x_label=add_x_label, negative_and_positive_values=negative_and_positive_values,
+                      massif_name_to_text=massif_name_to_text)
+            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,
                       add_x_label=True, negative_and_positive_values=True, massif_name_to_text=None):

extreme_trend/abstract_gev_trend_test.py

@@ -492,4 +492,19 @@ class AbstractGevTrendTest(object):
         old_mean = self.get_unconstrained_gev_params(year=self.years[-nb_years]).mean
         return last_mean - old_mean
 
+    def relative_change_in_mean_for_the_last_x_years(self, nb_years):
+        last_mean = self.unconstrained_estimator_gev_params_last_year.mean
+        old_mean = self.get_unconstrained_gev_params(year=self.years[-nb_years]).mean
+        return 100 * (last_mean - old_mean) / old_mean
+
+    def change_in_50_year_return_level_for_the_last_x_years(self, nb_years, return_period=50):
+        last_mean = self.unconstrained_estimator_gev_params_last_year.return_level(return_period=return_period)
+        old_mean = self.get_unconstrained_gev_params(year=self.years[-nb_years]).return_level(return_period=return_period)
+        return last_mean - old_mean
+
+    def relative_change_in_50_year_return_level_for_the_last_x_years(self, nb_years, return_period=50):
+        last_mean = self.unconstrained_estimator_gev_params_last_year.return_level(return_period=return_period)
+        old_mean = self.get_unconstrained_gev_params(year=self.years[-nb_years]).return_level(return_period=return_period)
+        return 100 * (last_mean - old_mean) / old_mean
+
 

extreme_trend/visualizers/study_visualizer_for_non_stationary_trends.py

@@ -177,8 +177,12 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
             trend_test_that_minimized_aic = sorted_trend_test[0]
             massif_name_to_trend_test_that_minimized_aic[massif_name] = trend_test_that_minimized_aic
             # Extract the stationary model that minimized AIC
-            stationary_trend_test_that_minimized_aic = [t for t in sorted_trend_test if type(t) in
-                                                        [GumbelVersusGumbel, GevStationaryVersusGumbel]][0]
+            stationary_trend_tests_that_minimized_aic = [t for t in sorted_trend_test if type(t) in
+                                                        [GumbelVersusGumbel, GevStationaryVersusGumbel]]
+            if len(stationary_trend_tests_that_minimized_aic) == 0:
+                stationary_trend_test_that_minimized_aic = None
+            else:
+                stationary_trend_test_that_minimized_aic = stationary_trend_tests_that_minimized_aic[0]
             massif_name_to_stationary_trend_test_that_minimized_aic[
                 massif_name] = stationary_trend_test_that_minimized_aic
             # Extract the Gumbel model that minimized AIC

extreme_trend/visualizers/utils.py

 from collections import OrderedDict
 
-from extreme_data.meteo_france_data.scm_models_data.utils import Season
+from extreme_data.meteo_france_data.scm_models_data.utils import Season, FrenchRegion
 from extreme_fit.model.margin_model.utils import \
     MarginFitMethod
 from extreme_trend.visualizers.study_visualizer_for_non_stationary_trends import \
@@ -12,11 +12,12 @@ def load_altitude_to_visualizer(altitudes, massif_names, model_subsets_for_uncer
                                 study_visualizer_class=StudyVisualizerForNonStationaryTrends,
                                 save_to_file=True,
                                 multiprocessing=True,
-                                season=Season.annual):
+                                season=Season.annual,
+                                french_region=FrenchRegion.alps):
     fit_method = MarginFitMethod.extremes_fevd_mle
     altitude_to_visualizer = OrderedDict()
     for altitude in altitudes:
-        study = study_class(altitude=altitude, multiprocessing=multiprocessing, season=season)
+        study = study_class(altitude=altitude, multiprocessing=multiprocessing, season=season, french_region=french_region)
         study_visualizer = study_visualizer_class(study=study, multiprocessing=multiprocessing,
                                                   save_to_file=save_to_file, uncertainty_massif_names=massif_names,
                                                   uncertainty_methods=uncertainty_methods,

projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/main_snowfall_article.py

@@ -3,7 +3,7 @@ from multiprocessing.pool import Pool
 
 import matplotlib as mpl
 
-from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day
+from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day, SafranPrecipitation1Day
 from projects.contrasting_trends_in_snow_loads.article2_snowfall_versus_time_and_altitude.shape_plot import shape_plot
 from projects.contrasting_trends_in_snow_loads.article2_snowfall_versus_time_and_altitude.snowfall_plot import \
     plot_snowfall_mean, plot_snowfall_change_mean
@@ -74,25 +74,25 @@ def intermediate_result(altitudes, massif_names=None,
             _ = compute_minimized_aic(visualizer)
 
     # Plots
-    validation_plot(altitude_to_visualizer, order_derivative=0)
-    validation_plot(altitude_to_visualizer, order_derivative=1)
+    # validation_plot(altitude_to_visualizer, order_derivative=0)
+    # validation_plot(altitude_to_visualizer, order_derivative=1)
     plot_snowfall_mean(altitude_to_visualizer)
     plot_selection_curves(altitude_to_visualizer, paper1=False)
     plot_snowfall_change_mean(altitude_to_visualizer)
-    shape_plot(altitude_to_visualizer)
+    # shape_plot(altitude_to_visualizer)
 
 
 def major_result():
     uncertainty_methods = [ConfidenceIntervalMethodFromExtremes.ci_mle][:]
     # massif_names = ['Beaufortain', 'Vercors']
     massif_names = None
-    study_classes = [SafranSnowfall1Day]
+    study_classes = [SafranSnowfall1Day, SafranPrecipitation1Day][::-1]
     model_subsets_for_uncertainty = None
     altitudes = paper_altitudes
-    # altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900]
+    altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900]
     # altitudes = [900, 1200, 1500, 1800][:2]
     # altitudes = [1800, 2100, 2400, 2700][:2]
-    altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000]
+    # altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000]
     # altitudes = draft_altitudes
 
     for study_class in study_classes:

projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/snowfall_plot.py

@@ -23,74 +23,87 @@ def fit_linear_regression(x, y):
 def plot_snowfall_change_mean(altitude_to_visualizer: Dict[int, StudyVisualizerForMeanValues]):
     visualizer = list(altitude_to_visualizer.values())[0]
     study = visualizer.study
-    # Plot the curve for the evolution of the mean
-    massif_name_to_a, massif_name_to_b, massif_name_to_r2_score = plot_mean(altitude_to_visualizer, derivative=True)
-    # Augmentation every km
-    massif_name_to_augmentation_every_km = {m: a * 1000 for m, a in massif_name_to_a.items()}
-    visualizer.plot_abstract_fast(massif_name_to_augmentation_every_km,
-                                  label='Augmentation of time derivative of mean annual maxima of {}\n for every km of elevation ({})'.format(
-                                      SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)], study.variable_unit),
-                                  add_x_label=False)
-    # Value at 2000 m
-    massif_name_to_mean_at_2000 = {m: a * 2000 + massif_name_to_b[m] for m, a in massif_name_to_a.items()}
-    visualizer.plot_abstract_fast(massif_name_to_mean_at_2000,
-                                  label='Time derivative of mean annual maxima \nof {} at 2000 m ({})'.format(
-                                      SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)], study.variable_unit),
-                                  add_x_label=False,
-                                  add_text=True, massif_name_to_text={k: str(v) for k,v in massif_name_to_r2_score.items()})
-    # Altitude for the change of dynamic
-    massif_name_to_altitude_change_dynamic = {m: - massif_name_to_b[m] / a for m, a in massif_name_to_a.items()}
-    # Keep only those that are in a reasonable range
-    massif_name_to_altitude_change_dynamic = {m: d for m, d in massif_name_to_altitude_change_dynamic.items()
-                                              if 0 < d < 3000}
-    visualizer.plot_abstract_fast(massif_name_to_altitude_change_dynamic,
-                                  label='Altitude for the change of dynamic (m)',
-                                  add_x_label=False, graduation=500,
-                                  add_text=True, massif_name_to_text={k: str(v) for k,v in massif_name_to_r2_score.items()})
-    # R2 score
-    # visualizer.plot_abstract_fast(massif_name_to_r2_score, label='r2 time derivative of the mean', graduation=0.1,
-    #                               add_x_label=False,
-    #                               negative_and_positive_values=False,
-    #                               )
+
+    variables = ['mean annual maxima', '50 year return level']
+    mean_indicators = [True, False]
+    for variable, mean_indicator in zip(variables, mean_indicators):
+        for relative in [False, True]:
+            if relative:
+                variable += str(' (relative)')
+            # Plot the curve for the evolution of the mean
+            massif_name_to_a, massif_name_to_b, massif_name_to_r2_score = plot_mean(altitude_to_visualizer,
+                                                                                    derivative=True, mean_indicator=mean_indicator,
+                                                                                    relative=relative)
+            # Augmentation every km
+            massif_name_to_augmentation_every_km = {m: a * 1000 for m, a in massif_name_to_a.items()}
+            massif_to_r2_score_text = {k: str(round(v, 2)) for k, v in massif_name_to_r2_score.items()}
+
+            visualizer.plot_abstract_fast(massif_name_to_augmentation_every_km,
+                                          label='Slope for changes of {} of {}\n for every km of elevation ({})'.format(
+                                              variable, SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)], study.variable_unit),
+                                          add_x_label=False)
+            # Value at 2000 m
+            massif_name_to_mean_at_2000 = {m: a * 2000 + massif_name_to_b[m] for m, a in massif_name_to_a.items()}
+            visualizer.plot_abstract_fast(massif_name_to_mean_at_2000,
+                                          label='Changes of {} \nof {} at 2000 m ({})'.format(
+                                              variable, SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)], study.variable_unit),
+                                          add_x_label=False,
+                                          add_text=True, massif_name_to_text=massif_to_r2_score_text)
+            # Altitude for the change of dynamic
+            massif_name_to_altitude_change_dynamic = {m: - massif_name_to_b[m] / a for m, a in massif_name_to_a.items()}
+            # Keep only those that are in a reasonable range
+            massif_name_to_altitude_change_dynamic = {m: d for m, d in massif_name_to_altitude_change_dynamic.items()
+                                                      if 0 < d < 3000}
+            visualizer.plot_abstract_fast(massif_name_to_altitude_change_dynamic,
+                                          label=('Altitude when the changes equal zero for %s (m)' % variable),
+                                          add_x_label=False, graduation=500,
+                                          add_text=True, massif_name_to_text=massif_to_r2_score_text)
 
 
 def plot_snowfall_mean(altitude_to_visualizer: Dict[int, StudyVisualizerForMeanValues]):
     visualizer = list(altitude_to_visualizer.values())[0]
     study = visualizer.study
-    # Plot the curve for the evolution of the mean
-    massif_name_to_a, massif_name_to_b, massif_name_to_r2_score = plot_mean(altitude_to_visualizer, derivative=False)
-    # Compute values of interest
-    massif_name_to_mean_at_2000 = {m: a * 2000 + massif_name_to_b[m] for m, a in massif_name_to_a.items()}
-    massif_name_to_mean_at_1000 = {m: a * 1000 + massif_name_to_b[m] for m, a in massif_name_to_a.items()}
-    massif_name_to_augmentation_every_km = {m: a * 1000 for m, a in massif_name_to_a.items()}
-    massif_name_to_relative_augmentation_between_2000_and_3000_every_km = {m: 100 * (a * 1000) / (a * 2000 + massif_name_to_b[m]) for m, a in massif_name_to_a.items()}
-    massif_name_to_relative_augmentation_between_1000_and_2000_every_km = {m: 100 * (a * 1000) / (a * 1000 + massif_name_to_b[m]) for m, a in massif_name_to_a.items()}
-
-    # Augmentation every km
-    ds = [massif_name_to_augmentation_every_km, massif_name_to_relative_augmentation_between_1000_and_2000_every_km,
-          massif_name_to_relative_augmentation_between_2000_and_3000_every_km]
-    prefixs = ['Augmentation', 'Relative augmentation between 1000m and 2000m', 'Relative augmentation between 2000m and 3000m']
-    graduations = [1.0, 10.0, 10.0]
-    for d, prefix in zip(ds, prefixs):
-        visualizer.plot_abstract_fast(d,
-                                      graduation=1.0,
-                                      label=prefix + ' of mean annual maxima of {} \nfor every km of elevation ({})'.format(
-                                          SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)], study.variable_unit),
+
+    variables = ['mean annual maxima', '50 year return level']
+    mean_indicators = [True, False]
+    for variable, mean_indicator in zip(variables, mean_indicators):
+        # Plot the curve for the evolution of the mean
+        massif_name_to_a, massif_name_to_b, massif_name_to_r2_score = plot_mean(altitude_to_visualizer,
+                                                                                derivative=False,
+                                                                                mean_indicator=mean_indicator)
+        # Compute values of interest
+        massif_name_to_mean_at_2000 = {m: a * 2000 + massif_name_to_b[m] for m, a in massif_name_to_a.items()}
+        massif_name_to_augmentation_every_km = {m: a * 1000 for m, a in massif_name_to_a.items()}
+        massif_name_to_relative_augmentation_between_2000_and_3000_every_km = {
+            m: 100 * (a * 1000) / (a * 2000 + massif_name_to_b[m]) for m, a in massif_name_to_a.items()}
+        massif_name_to_relative_augmentation_between_1000_and_2000_every_km = {
+            m: 100 * (a * 1000) / (a * 1000 + massif_name_to_b[m]) for m, a in massif_name_to_a.items()}
+        massif_to_r2_score_text = {k: str(round(v, 2)) for k, v in massif_name_to_r2_score.items()}
+        # Augmentation every km
+        ds = [massif_name_to_augmentation_every_km, massif_name_to_relative_augmentation_between_1000_and_2000_every_km,
+              massif_name_to_relative_augmentation_between_2000_and_3000_every_km]
+        prefixs = ['Augmentation', 'Relative augmentation between 1000m and 2000m',
+                   'Relative augmentation between 2000m and 3000m']
+        graduations = [1.0, 10.0, 10.0]
+        for d, prefix, graduation in zip(ds, prefixs, graduations):
+            visualizer.plot_abstract_fast(d,
+                                          graduation=10.0,
+                                          label=prefix + ' of {} of {} \nfor every km of elevation ({})'.format(
+                                              variable, SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)],
+                                              study.variable_unit),
+                                          add_x_label=False, negative_and_positive_values=False,
+                                          add_text=True,
+                                          massif_name_to_text=massif_to_r2_score_text
+                                          )
+        # Value at 2000 m
+        visualizer.plot_abstract_fast(massif_name_to_mean_at_2000, label='{} of {} at 2000 m ()'.format(
+            variable, SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)], study.variable_unit),
                                       add_x_label=False, negative_and_positive_values=False,
-                                      add_text=True,
-                                      massif_name_to_text={k: str(v) for k, v in massif_name_to_r2_score.items()}
-                                      )
-    # Value at 2000 m
-    visualizer.plot_abstract_fast(massif_name_to_mean_at_2000, label='Mean annual maxima of {} at 2000 m ()'.format(
-        SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)], study.variable_unit),
-                                  add_x_label=False, negative_and_positive_values=False,
-                                  add_text=True, massif_name_to_text={k: str(v) for k,v in massif_name_to_r2_score.items()})
-    # # R2 score
-    # visualizer.plot_abstract_fast(massif_name_to_r2_score, label='r2 mean', graduation=0.1,
-    #                               add_x_label=False, negative_and_positive_values=False)
-
-
-def plot_mean(altitude_to_visualizer: Dict[int, StudyVisualizerForMeanValues], derivative=False):
+                                      add_text=True, massif_name_to_text=massif_to_r2_score_text)
+
+
+def plot_mean(altitude_to_visualizer: Dict[int, StudyVisualizerForMeanValues], derivative=False, mean_indicator=True,
+              relative=False):
     ax = plt.gca()
     massif_name_to_linear_regression_result = {}
 
@@ -104,29 +117,44 @@ def plot_mean(altitude_to_visualizer: Dict[int, StudyVisualizerForMeanValues], d
         if len(altitudes_massif) >= 2:
             trend_tests = [altitude_to_visualizer[a].massif_name_to_trend_test_that_minimized_aic[massif_name]
                            for a in altitudes_massif]
+            nb_years = 50
+            return_period = 50
+            if mean_indicator:
+                indicator_str = 'mean'
+            else:
+                indicator_str = '50 year return level'
+            res = [(a, t)
+                   for i, (a, t) in enumerate(zip(altitudes_massif, trend_tests))
+                   if not t.unconstrained_model_is_stationary]
             if derivative:
-                nb_years = 10
-                res = [(a, t.change_in_mean_for_the_last_x_years(nb_years=nb_years))
-                       for i, (a, t) in enumerate(zip(altitudes_massif, trend_tests))
-                       if t.is_significant]
+                if mean_indicator:
+                    if relative:
+                        res = [(a, t.relative_change_in_mean_for_the_last_x_years(nb_years=nb_years)) for a, t in res]
+                    else:
+                        res = [(a, t.change_in_mean_for_the_last_x_years(nb_years=nb_years)) for a, t in res]
+                else:
+                    if relative:
+                        res = [(a, t.relative_change_in_50_year_return_level_for_the_last_x_years(nb_years=nb_years, return_period=return_period)) for a, t in res]
+                    else:
+                        res = [(a, t.change_in_50_year_return_level_for_the_last_x_years(nb_years=nb_years, return_period=return_period)) for a, t in res]
                 if len(res) > 0:
                     altitudes_values, values = zip(*res)
                 else:
                     altitudes_values, values = [], []
-                moment = 'Change in the last {} years  \nfor significative models'.format(nb_years)
-                # res = [(a, t.change_in_mean_for_the_last_x_years(nb_years=nb_years))
-                #        for i, (a, t) in enumerate(zip(altitudes_massif, trend_tests))
-                #        if not t.unconstrained_model_is_stationary]
-                # altitudes_values, values = zip(*res)
-                # moment = 'Change in the last {} years  \nfor non-stationary models'.format(nb_years)
+                indicator_str = 'Change of the {} for the last {} years  \nfor non-stationary models'.format(indicator_str, nb_years)
+                if relative:
+                    indicator_str += ' (relative)'
             else:
-                moment = 'mean'
-                values = [t.unconstrained_estimator_gev_params_last_year.mean for t in trend_tests]
+                if mean_indicator:
+                    values = [t.unconstrained_estimator_gev_params_last_year.mean for t in trend_tests]
+                else:
+                    values = [t.unconstrained_estimator_gev_params_last_year.return_level(return_period=return_period)
+                              for t in trend_tests]
                 altitudes_values = altitudes_massif
             # Plot
             if len(altitudes_values) >= 2:
                 massif_name_to_linear_regression_result[massif_name] = fit_linear_regression(altitudes_values, values)
-                plot_values_against_altitudes(ax, altitudes_values, massif_id, massif_name, moment, study, values,
+                plot_values_against_altitudes(ax, altitudes_values, massif_id, massif_name, indicator_str, study, values,
                                               visualizer)
     ax.legend(prop={'size': 7}, ncol=3)
     visualizer.show_or_save_to_file(dpi=500, add_classic_title=False)