[contrasting] add plot for temperature when maxima happen

extreme_data/meteo_france_data/scm_models_data/abstract_study.py

@@ -251,6 +251,14 @@ class AbstractStudy(object):
             year_to_annual_maxima[year] = time_serie.argmax(axis=0)
         return year_to_annual_maxima
 
+    @cached_property
+    def year_to_annual_maxima_tuple_indices_for_daily_time_series(self):
+        year_to_annual_maxima_indices_for_daily_time_series = OrderedDict()
+        for year in self.ordered_years:
+            l = [(idx, i) for i, idx in enumerate(self.year_to_annual_maxima_index[year])]
+            year_to_annual_maxima_indices_for_daily_time_series[year] = l
+        return year_to_annual_maxima_indices_for_daily_time_series
+
     @cached_property
     def massif_name_to_annual_maxima_ordered_index(self):
         massif_name_to_annual_maxima_ordered_index = OrderedDict()

projects/contrasting_trends_in_snow_loads/gorman_figures/figure 3/figure3_daily snowfall fraction.py

@@ -9,10 +9,11 @@ Altitudes_groups = [[300, 600, 900], [1200, 1500, 1800], [2100, 2400, 2700], [30
 
 
 def compute_smoother_snowfall_fraction(altitudes, year_min, year_max, lim, massif_names=None):
-    all_time_series = np.concatenate([get_time_series(altitude, year_min, year_max, massif_names) for altitude in altitudes], axis=1)
+    all_time_series = np.concatenate(
+        [get_time_series(altitude, year_min, year_max, massif_names) for altitude in altitudes], axis=1)
     temperature_array, snowfall_array, rainfall_array = all_time_series
     # nb_bins_for_one_celsius = 4
-    nb_bins_for_one_celsius = 2
+    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 = []
@@ -20,20 +21,31 @@ def compute_smoother_snowfall_fraction(altitudes, year_min, year_max, lim, massi
         mask = inds == j
         fraction = np.mean(snowfall_array[mask]) / np.mean(snowfall_array[mask] + rainfall_array[mask])
         snowfall_fractions.append(fraction)
-    x = bins[:-1] + 0.125
-    return np.array(snowfall_fractions), x
+    center_bins = bins[:-1] + 0.5 * (1 / nb_bins_for_one_celsius)
+    return np.array(snowfall_fractions), center_bins
 
 
-def daily_snowfall_fraction(ax, altitudes, year_min=1959, year_max=2019, massif_names=None):
+def daily_snowfall_fraction(ax, altitudes, year_min=1959, year_max=2019, massif_names=None, sigma=1.0):
     lim = 6
-    snowfall_fractions, x = compute_smoother_snowfall_fraction(altitudes, year_min, year_max, lim=lim,
-                                                               massif_names=massif_names)
-    sigma = 1.0
+    snowfall_fractions, center_bins = compute_smoother_snowfall_fraction(altitudes, year_min, year_max, lim=lim,
+                                                                         massif_names=massif_names)
+
+    # Smooth the data afterwards
     snowfall_fractions = 100 * gaussian_filter1d(snowfall_fractions, sigma=sigma)
+    # Extend the data by interpolation
+    center_bins_interpolation = np.arange(-lim, lim, step=0.1)
+    snowfall_fractions_interpolated = np.interp(center_bins_interpolation, center_bins, snowfall_fractions)
+    # Compute the optimal temperature
+    h = np.exp(0.06 * center_bins_interpolation) * snowfall_fractions_interpolated
+    temperature_optimal = center_bins_interpolation[np.argmax(h)]
+
+
     # Plot results
     label = '{}-{} at {}'.format(year_min, year_max, ' & '.join(['{} m'.format(a) for a in altitudes]))
     ax.set_xlim([-lim, lim])
-    ax.plot(x, snowfall_fractions, label=label, linewidth=5)
+    p = ax.plot(center_bins, snowfall_fractions, label=label, linewidth=5)
+    color = p[0].get_color()
+    ax.plot([temperature_optimal], [0], color=color, marker='x', markersize=10)
     ax.set_ylabel('Snowfall fraction (%)')
     end_plot(ax, sigma)
 
@@ -46,14 +58,16 @@ def snowfall_fraction_difference_between_periods(ax, altitudes, year_min=1959, y
     last_percentage_of_interest = 0.5
     mask = snowfall_fractions_past > last_percentage_of_interest
     x = x[mask]
-    snowfall_percentage = 100 * (snowfall_fractions_recent[mask] - snowfall_fractions_past[mask]) / snowfall_fractions_past[mask]
+    snowfall_percentage = 100 * (snowfall_fractions_recent[mask] - snowfall_fractions_past[mask]) / \
+                          snowfall_fractions_past[mask]
     sigma = 1.0
     snowfall_percentage = gaussian_filter1d(snowfall_percentage, sigma=sigma)
     label = '{}'.format(' & '.join(['{} m'.format(a) for a in altitudes]))
     ax.plot(x, snowfall_percentage, label=label, linewidth=4)
     ax.set_ylabel('Relative change of Snowfall fraction {}-{}\n'
                   'compared to Snowfall fraction {}-{}\n'
-                  'until the latter reach a fraction of 0.5 (%)'.format(year_middle + 1, year_max, year_min, year_middle))
+                  'until the latter reach a fraction of 0.5 (%)'.format(year_middle + 1, year_max, year_min,
+                                                                        year_middle))
     end_plot(ax, sigma)
 
 
@@ -68,9 +82,6 @@ def get_time_series(altitude, year_min, year_max, massif_names=None):
     temperature_study = SafranTemperature(altitude=altitude, year_min=year_min, year_max=year_max)
     study_rainfall = SafranRainfall1Day(altitude=altitude, year_min=year_min, year_max=year_max)
     study_snowfall = SafranSnowfall1Day(altitude=altitude, year_min=year_min, year_max=year_max)
-    studies = [temperature_study, study_rainfall, study_snowfall]
-    if massif_names is None:
-        all
     all_time_series = [temperature_study.all_daily_series,
                        study_snowfall.all_daily_series,
                        study_rainfall.all_daily_series]
@@ -85,12 +96,12 @@ def get_time_series(altitude, year_min, year_max, massif_names=None):
 
 def daily_snowfall_fraction_wrt_altitude(fast=True):
     ax = plt.gca()
+    groups = [[a] for a in [900, 1200, 1500][:]]
     groups = [[1800]] if fast else Altitudes_groups
-    groups = [[a] for a in [900, 1200, 1500]]
     massif_names = None
-    massif_names = ['Vercors']
+    # massif_names = ['Vercors']
     for altitudes in groups:
-        daily_snowfall_fraction(ax, altitudes, massif_names=massif_names, year_max=1999)
+        daily_snowfall_fraction(ax, altitudes, massif_names=massif_names)
     plt.show()
 
 

projects/contrasting_trends_in_snow_loads/temperature_of_snowfall_maxima.py

+from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day, SafranTemperature
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+def distribution_temperature_of_maxima_of_snowfall(altitudes):
+    # Load the temperature corresponding to the maxima of snowfall
+    altitude_to_maxima_temperature = {}
+    for altitude in altitudes:
+        print(altitude)
+        snowfall_study = SafranSnowfall1Day(altitude=altitude)
+        temperature_study = SafranTemperature(altitude=altitude)
+        # temperature_study = SafranTemperature(altitude=altitude)
+        maxima_temperatures = []
+        for year in snowfall_study.ordered_years[:2]:
+            annual_maxima_index = snowfall_study.year_to_annual_maxima_tuple_indices_for_daily_time_series[year]
+            a = temperature_study.year_to_daily_time_serie_array[year]
+            temp = [a[tuple_idx] for tuple_idx in annual_maxima_index]
+            maxima_temperatures.append(temp)
+        altitude_to_maxima_temperature[altitude] = np.concatenate(maxima_temperatures)
+
+    # Plot the box plot
+    ax = plt.gca()
+    width = 150
+    ax.boxplot(altitude_to_maxima_temperature.values(), positions=altitudes, widths=width)
+    ax.set_xlim([min(altitudes) - width, max(altitudes) + width])
+    ax.set_ylabel('Temperature for maxima of snowfall ({})'.format(temperature_study.variable_class.UNIT))
+    ax.set_xlabel('Altitude (m)')
+    ax.grid()
+
+    plt.show()
+
+
+if __name__ == '__main__':
+    distribution_temperature_of_maxima_of_snowfall(altitudes=[900, 1200, 1500, 1800, 2100, 2400, 2700, 3000])