From e9e777bc65977fbee11c5daad8dcb8c90a97e56b Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Wed, 25 Mar 2020 19:24:43 +0100
Subject: [PATCH] [contrasting project] add main_distribution_wps.py to study
 the (spatial and temporal) distribution of weather patterns that are causing
 maxima in the different region of the French Alps

---
 .../scm_models_data/abstract_study.py         | 33 ++++++++
 .../main_distribution_wps.py                  | 75 +++++++++++++++++++
 2 files changed, 108 insertions(+)
 create mode 100644 projects/contrasting_trends_in_snow_loads/weather_types_analysis/main_distribution_wps.py

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 ddafc5d3..68873376 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
@@ -77,6 +77,10 @@ class AbstractStudy(object):
 
     """ Time """
 
+    @property
+    def nb_years(self):
+        return self.year_max - self.year_min + 1
+
     @cached_property
     def year_to_days(self) -> OrderedDict:
         # Map each year to the 'days since year-08-01 06:00:00'
@@ -101,6 +105,35 @@ class AbstractStudy(object):
             year_to_wps[year] = self.df_weather_types.loc[days].iloc[:, 0].values
         return year_to_wps
 
+    def wps_for_top_annual_maxima(self, nb_top, massif_ids):
+        d = self.massif_id_to_df_ordered_by_maxima
+        wps = pd.Series(np.concatenate([d[massif_id]['WP'].values[:nb_top]
+                                        for massif_id in massif_ids]))
+        s_normalized = wps.value_counts(normalize=True) * 100
+        s_normalized = s_normalized.round()
+        s_not_normalized = wps.value_counts()
+        df = pd.concat([s_normalized, s_not_normalized], axis=1)
+        df.columns = ['Percentage', 'Nb massifs concerned']
+        df.index.name = 'Number Top={}'.format(nb_top)
+        return df
+
+    @cached_property
+    def massif_id_to_df_ordered_by_maxima(self):
+        df_annual_maxima = pd.DataFrame(self.year_to_annual_maxima)
+        df_wps = pd.DataFrame(self.year_to_wp_for_annual_maxima)
+        massif_id_to_df_ordered_by_maxima = {}
+        for massif_id, s_annual_maxima in df_annual_maxima.iterrows():
+            s_annual_maxima.sort_values(inplace=True, ascending=False)
+            d = {
+                'Year': s_annual_maxima.index,
+                'Maxima': s_annual_maxima.values,
+                'WP': df_wps.loc[massif_id, s_annual_maxima.index],
+            }
+            df = pd.DataFrame(d)
+            df.set_index('Year', inplace=True)
+            massif_id_to_df_ordered_by_maxima[massif_id] = df
+        return massif_id_to_df_ordered_by_maxima
+
     @cached_property
     def year_to_wp_for_annual_maxima(self):
         year_to_wp_for_annual_maxima = OrderedDict()
diff --git a/projects/contrasting_trends_in_snow_loads/weather_types_analysis/main_distribution_wps.py b/projects/contrasting_trends_in_snow_loads/weather_types_analysis/main_distribution_wps.py
new file mode 100644
index 00000000..8191e44a
--- /dev/null
+++ b/projects/contrasting_trends_in_snow_loads/weather_types_analysis/main_distribution_wps.py
@@ -0,0 +1,75 @@
+import pandas as pd
+import numpy as np
+
+from extreme_data.meteo_france_data.scm_models_data.abstract_extended_study import AbstractExtendedStudy
+from extreme_data.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoad3Days, CrocusSnowLoad1Day
+from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranPrecipitation1Day
+
+
+def main_spatial_distribution_wps(study_class, year_min=1954, year_max=2008):
+    study = study_class(altitude=1800, year_min=year_min, year_max=year_max)
+    for region_name in AbstractExtendedStudy.region_names:
+        massifs_ids = AbstractExtendedStudy.region_name_to_massif_ids[region_name]
+        print('\n \n', region_name, '\n')
+        for nb_top in [study.nb_years, 5, 1][1:2]:
+            print(study.wps_for_top_annual_maxima(nb_top=nb_top, massif_ids=massifs_ids), '\n')
+
+
+"""
+Some analysis:
+
+At altitude 1800m, 
+for the precipitation in 1 day
+between year_min=1954, year_max=2008
+
+If we look at the 5 biggest maxima for each massif, 
+and look to which Weather pattern they correspond
+then we do some percentage for each climatic region
+
+                    Percentage  Nb massifs concerned
+ Northern Alps 
+Steady Oceanic            77.0                    27
+Atlantic Wave             11.0                     4
+
+ Central Alps 
+Steady Oceanic            57.0                    20
+South Circulation         17.0                     6
+East Return               14.0                     5
+
+ Southern Alps 
+South Circulation         43.0                    13
+Central Depression        37.0                    11
+East Return               17.0                     5
+
+ Extreme South Alps 
+Central Depression        53.0                     8
+South Circulation         47.0                     7 
+
+"""
+
+
+def main_temporal_distribution_wps(study_class, year_min=1954, year_max=2008):
+    altitude = 1800
+    study_before = study_class(altitude=altitude, year_min=year_min, year_max=1981)
+    study_after = study_class(altitude=altitude, year_min=1981, year_max=2008)
+    for region_name in AbstractExtendedStudy.region_names:
+        massifs_ids = AbstractExtendedStudy.region_name_to_massif_ids[region_name]
+        print('\n \n', region_name, '\n')
+        for nb_top in [study_before.nb_years, 10, 5, 1][1:2]:
+            print(study_before.wps_for_top_annual_maxima(nb_top=nb_top, massif_ids=massifs_ids), '\n')
+            print(study_after.wps_for_top_annual_maxima(nb_top=nb_top, massif_ids=massifs_ids), '\n')
+
+"""
+There is no real stationarity in the percentage of the kind of storms that are causing extreme.
+
+the reparittion of storm before and after (no matter the nb_top consider 10, or all) are the same.
+even for the local region it is the same.
+
+-> so what is really changing is probably the strength associated to each kind of storm.
+"""
+
+
+if __name__ == '__main__':
+    study_class = [CrocusSnowLoad1Day, SafranPrecipitation1Day][-1]
+    # main_spatial_distribution_wps(study_class)
+    main_temporal_distribution_wps(study_class)
-- 
GitLab