[contrasting project] add figure 1

extreme_data/meteo_france_data/scm_models_data/abstract_study.py

@@ -238,12 +238,20 @@ class AbstractStudy(object):
             ordered_index = [self.year_to_annual_maxima_index[year][i] for year in years]
             massif_name_to_annual_maxima_ordered_index[massif_name] = ordered_index
         return massif_name_to_annual_maxima_ordered_index
+    
+    @cached_property
+    def massif_name_to_annual_maxima(self):
+        massif_name_to_annual_maxima = OrderedDict()
+        for i, massif_name in enumerate(self.study_massif_names):
+            maxima = np.array([self.year_to_annual_maxima[year][i] for year in self.ordered_years])
+            massif_name_to_annual_maxima[massif_name] = maxima
+        return massif_name_to_annual_maxima
 
     @cached_property
     def massif_name_to_annual_maxima_ordered_years(self):
         massif_name_to_annual_maxima_ordered_years = OrderedDict()
-        for i, massif_name in enumerate(self.study_massif_names):
-            maxima = np.array([self.year_to_annual_maxima[year][i] for year in self.ordered_years])
+        for massif_name in self.study_massif_names:
+            maxima = self.massif_name_to_annual_maxima[massif_name]
             annual_maxima_ordered_index = np.argsort(maxima)
             annual_maxima_ordered_years = [self.ordered_years[idx] for idx in annual_maxima_ordered_index]
             massif_name_to_annual_maxima_ordered_years[massif_name] = annual_maxima_ordered_years

extreme_fit/distribution/gev/gev_params.py

@@ -37,6 +37,9 @@ class GevParams(AbstractExtremeParams):
         else:
             return r.qgev(p, self.location, self.scale, self.shape)[0]
 
+    def return_level(self, return_period):
+        return self.quantile(1 - 1 / return_period)
+
     def density(self, x):
         if self.has_undefined_parameters:
             return np.nan

projects/contrasting_trends_in_snow_loads/gorman_figures/figure1_mean_ratio_return_level_ratio.py

+from collections import OrderedDict
+
+import numpy as np
+
+from extreme_data.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
+from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day
+from extreme_data.meteo_france_data.scm_models_data.visualization.create_shifted_cmap import get_shifted_map, \
+    get_colors, shiftedColorMap
+from extreme_fit.estimator.margin_estimator.utils import fitted_stationary_gev
+from extreme_fit.model.margin_model.linear_margin_model.abstract_temporal_linear_margin_model import \
+    TemporalMarginFitMethod
+
+import matplotlib.pyplot as plt
+
+from extreme_trend.visualizers.study_visualizer_for_non_stationary_trends import StudyVisualizerForNonStationaryTrends
+
+
+def mean_ratio(altitude):
+    pass
+
+
+def massif_name_to_return_level(altitude, year_min, year_max, return_period):
+    study = SafranSnowfall1Day(altitude=altitude, year_min=year_min, year_max=year_max)
+    massif_name_to_return_levels = OrderedDict()
+    for massif_name, maxima in study.massif_name_to_annual_maxima.items():
+        gev_param = fitted_stationary_gev(maxima, fit_method=TemporalMarginFitMethod.extremes_fevd_mle)
+        return_level = gev_param.return_level(return_period=return_period)
+        massif_name_to_return_levels[massif_name] = return_level
+    return massif_name_to_return_levels
+
+
+def plot_return_level_ratio(altitude, return_period=30, year_min=1959, year_middle=1989, year_max=2019):
+    massif_name_to_return_level_past = massif_name_to_return_level(altitude, year_min, year_middle, return_period)
+    massif_name_to_return_level_recent = massif_name_to_return_level(altitude, year_middle + 1, year_max,
+                                                                     return_period)
+    massif_name_to_ratio = {m: v / massif_name_to_return_level_past[m] for m, v in
+                            massif_name_to_return_level_recent.items()}
+    max_ratio = max([e for e in massif_name_to_ratio.values()])
+    min_ratio = min([e for e in massif_name_to_ratio.values()])
+    midpoint = (max_ratio - 1.0) / (max_ratio - 0)
+    num = 3
+    ticks = np.concatenate([np.linspace(0, midpoint, num), np.linspace(midpoint, 1, num)])
+    labels = np.concatenate([np.linspace(min_ratio, 1.0, num), np.linspace(1.0, max_ratio, num)])
+    cmap = shiftedColorMap(plt.cm.bwr, midpoint=midpoint, name='shifted')
+    ax = plt.gca()
+    print(massif_name_to_ratio)
+    massif_name_to_color = {m: get_colors([v], cmap, min_ratio, max_ratio)[0]
+                for m, v in massif_name_to_ratio.items()}
+    ticks_values_and_labels = ticks, labels
+    AbstractStudy.visualize_study(ax=ax,
+                                  massif_name_to_value=massif_name_to_ratio,
+                                  massif_name_to_color=massif_name_to_color,
+                                  replace_blue_by_white=True,
+                                  axis_off=False,
+                                  cmap=cmap,
+                                  show_label=False,
+                                  add_colorbar=True,
+                                  show=False,
+                                  vmin=min_ratio,
+                                  vmax=max_ratio,
+                                  ticks_values_and_labels=ticks_values_and_labels
+                                  )
+    plt.show()
+
+
+if __name__ == '__main__':
+    plot_return_level_ratio(altitude=1800)

projects/contrasting_trends_in_snow_loads/snowfall fraction of total precipitation/daily snowfall fraction.py → projects/contrasting_trends_in_snow_loads/gorman_figures/figure3_daily snowfall fraction.py

@@ -8,51 +8,50 @@ from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranR
 Altitudes_groups = [[300, 600, 900], [1200, 1500, 1800], [2100, 2400, 2700], [3000, 3300, 3600]]
 
 
-def compute_smoother_snowfall_fraction(altitudes, year_min, year_max):
+def compute_smoother_snowfall_fraction(altitudes, year_min, year_max, lim):
     all_time_series = np.concatenate([get_time_series(altitude, year_min, year_max) for altitude in altitudes], axis=1)
     temperature_array, snowfall_array, rainfall_array = all_time_series
-    paper_setting = False
-    if paper_setting:
-        lim = 10
-        sigma = 0.5
-    else:
-        lim = 7.5
-        sigma = 1
-    bins = np.linspace(-lim, lim, num=(4 * 2 * lim) + 1)
+    nb_bins_for_one_celsius = 4
+    bins = np.linspace(-lim, lim, num=int((nb_bins_for_one_celsius * 2 * lim) + 1))
     inds = np.digitize(temperature_array, bins)
     snowfall_fractions = []
     for j in range(1, len(bins)):
         mask = inds == j
         fraction = np.mean(snowfall_array[mask]) / np.mean(snowfall_array[mask] + rainfall_array[mask])
         snowfall_fractions.append(fraction)
-    new_snowfall_fractions = 100 * gaussian_filter1d(snowfall_fractions, sigma=sigma)
     x = bins[:-1] + 0.125
-    return new_snowfall_fractions, sigma, x
+    return np.array(snowfall_fractions), x
 
 
 def daily_snowfall_fraction(ax, altitudes, year_min=1959, year_max=2019):
-    snowfall_fractions, sigma, x = compute_smoother_snowfall_fraction(altitudes, year_min, year_max)
+    snowfall_fractions, x = compute_smoother_snowfall_fraction(altitudes, year_min, year_max, lim=7)
+    sigma = 1.0
+    snowfall_fractions = 100 * gaussian_filter1d(snowfall_fractions, sigma=sigma)
     # Plot results
     label = '{}-{} at {}'.format(year_min, year_max, ' & '.join(['{} m'.format(a) for a in altitudes]))
     ax.plot(x, snowfall_fractions, label=label, linewidth=4)
     ax.set_ylabel('Snowfall fraction (%)')
-    ax.set_xlabel('Daily surface air temperature (Celsius)')
-    ax.set_title('Snowfall fraction is smoothed with a gaussian filter with standard deviation={}'.format(sigma))
-    ax.grid(b=True)
-    ax.legend()
+    end_plot(ax, sigma)
 
 
 def snowfall_fraction_difference_between_periods(ax, altitudes, year_min=1959, year_middle=1989, year_max=2019):
-    snowfall_fractions_past, sigma, x = compute_smoother_snowfall_fraction(altitudes, year_min, year_middle)
-    snowfall_fractions_recent, sigma, x = compute_smoother_snowfall_fraction(altitudes, year_middle + 1, year_max)
-    mask = snowfall_fractions_past == 0
+    lim = 4
+    snowfall_fractions_past, x = compute_smoother_snowfall_fraction(altitudes, year_min, year_middle, lim=lim)
+    snowfall_fractions_recent, x = compute_smoother_snowfall_fraction(altitudes, year_middle + 1, year_max, lim=lim)
+    # Ensure that we do not divide by zero
+    mask = snowfall_fractions_past > 0
     x = x[mask]
     snowfall_ratio = snowfall_fractions_recent[mask] / snowfall_fractions_past[mask]
+    sigma = 2.0
+    snowfall_ratio = gaussian_filter1d(snowfall_ratio, sigma=sigma)
     label = '{}'.format(' & '.join(['{} m'.format(a) for a in altitudes]))
     ax.plot(x, snowfall_ratio, label=label, linewidth=4)
-    ax.set_ylabel('Ratio of Snowfall fraction {}-{} divided by Snowfall fraction {}-{}'.format(year_min, year_middle,
-                                                                                               year_middle + 1,
-                                                                                               year_max))
+    ax.set_ylabel('Ratio of Snowfall fraction {}-{}\n'
+                  'divided by Snowfall fraction {}-{}'.format(year_middle + 1, year_max, year_min, year_middle))
+    end_plot(ax, sigma)
+
+
+def end_plot(ax, sigma):
     ax.set_xlabel('Daily surface air temperature (Celsius)')
     ax.set_title('Snowfall fraction is smoothed with a gaussian filter with standard deviation={}'.format(sigma))
     ax.grid(b=True)
@@ -94,6 +93,6 @@ def daily_snowfall_fraction_wrt_time():
 
 
 if __name__ == '__main__':
-    # daily_snowfall_fraction_wrt_altitude(fast=True)
+    daily_snowfall_fraction_wrt_altitude(fast=False)
     # daily_snowfall_fraction_wrt_time()
-    ratio_of_past_and_recent_daily_snowfall_fraction_wrt_altitude()
+    # ratio_of_past_and_recent_daily_snowfall_fraction_wrt_altitude()