diff --git a/extreme_data/meteo_france_data/scm_models_data/abstract_study.py b/extreme_data/meteo_france_data/scm_models_data/abstract_study.py
index 4d95d941c08143caa3e0a54b35ee4617f6fe87a7..424bfb1441fcb47e41b32e30b4ba86ca6e7fb1b7 100644
--- a/extreme_data/meteo_france_data/scm_models_data/abstract_study.py
+++ b/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()
diff --git a/projects/contrasting_trends_in_snow_loads/gorman_figures/figure 3/figure3_daily snowfall fraction.py b/projects/contrasting_trends_in_snow_loads/gorman_figures/figure 3/figure3_daily snowfall fraction.py
index 672b84086c380a951f1825f70b3e87f5825960d0..1243fc450ada9b103f94d9858ce7325066ac72ca 100644
--- a/projects/contrasting_trends_in_snow_loads/gorman_figures/figure 3/figure3_daily snowfall fraction.py
+++ b/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()
diff --git a/projects/contrasting_trends_in_snow_loads/temperature_of_snowfall_maxima.py b/projects/contrasting_trends_in_snow_loads/temperature_of_snowfall_maxima.py
new file mode 100644
index 0000000000000000000000000000000000000000..6de908ab65930a8ed64fb0fa927d12ab5a061926
--- /dev/null
+++ b/projects/contrasting_trends_in_snow_loads/temperature_of_snowfall_maxima.py
@@ -0,0 +1,35 @@
+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])