[contrasting project] add daily snowfall fraction plot. add test_all_data

extreme_data/meteo_france_data/scm_models_data/abstract_study.py

@@ -198,8 +198,10 @@ class AbstractStudy(object):
         return list(chain(*list(self.year_to_days.values())))
 
     @property
-    def all_daily_series(self):
-        all_daily_series = np.concatenate(list(self.year_to_daily_time_serie_array.values()))
+    def all_daily_series(self) -> np.ndarray:
+        """Return an array of approximate shape (total_number_of_days, 23) x """
+        all_daily_series = np.concatenate([ time_serie_array
+                                            for time_serie_array in self.year_to_daily_time_serie_array.values()])
         assert len(all_daily_series) == len(self.all_days)
         return all_daily_series
 

projects/contrasting_trends_in_snow_loads/snowfall fraction of total precipitation/daily snowfall fraction.py

+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.ndimage import gaussian_filter, gaussian_filter1d
+
+from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranRainfall1Day, SafranTemperature, \
+    SafranSnowfall1Day
+
+
+def plot_snowfall_fraction(ax, altitude, temperature_array, snowfall_array, rainfall_array):
+    if ax is None:
+        ax = plt.gca()
+    lim = 10
+    bins = np.linspace(-lim, lim, num=(4 * 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=0.5)
+    x = bins[:-1] + 0.125
+
+    ax.plot(x, new_snowfall_fractions, label=altitude)
+    ax.set_ylabel('Snowfall fraction (%)')
+    ax.set_xlabel('Daily surface air temperature (T)')
+    ax.legend()
+
+
+def daily_snowfall_fraction(ax, altitude, year_min=1959, year_max=2019):
+    temperature_study = SafranTemperature(altitude=altitude)
+    study_rainfall = SafranRainfall1Day(altitude=altitude)
+    study_snowfall = SafranSnowfall1Day(altitude=altitude)
+    all_time_series = [temperature_study.all_daily_series, study_snowfall.all_daily_series,
+                       study_rainfall.all_daily_series]
+    plot_snowfall_fraction(ax, *[np.concatenate(t) for t in all_time_series])
+
+
+def daily_snowfall_fraction_wrt_altitude():
+    ax = plt.gca()
+    for altitude in [900, 1800, 2700]:
+        daily_snowfall_fraction(ax, altitude)
+    plt.show()
+
+
+if __name__ == '__main__':
+    daily_snowfall_fraction_wrt_altitude()
+    # daily_snowfall_fraction(altitude=900)

test/test_extreme_data/test_meteo_france_data/test_SCM_study.py

@@ -82,6 +82,12 @@ class TestSCMSafranSnowfall(TestSCMStudy):
         # Assert that the massif names are the same between SAFRAN and the coordinate file
         assert not set(self.study.study_massif_names).symmetric_difference(set(df_centroid['NOM']))
 
+    def test_all_data(self):
+        all_daily_series = self.study.all_daily_series
+        self.assertEqual(all_daily_series.ndim, 2)
+        self.assertEqual(all_daily_series.shape[1], 23)
+        self.assertEqual(all_daily_series.shape[0], 22280)
+
 
 class TestSCMPrecipitation(TestSCMStudy):