From 1ed4faf429d44645e6ece55d5b2866b53fe1f04c Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Wed, 15 Apr 2020 09:59:22 +0200
Subject: [PATCH] [contrasting project] add normalized precipitation. add
 folder spatial trends

---
 .../scm_models_data/abstract_study.py         | 21 ++++--
 .../scm_models_data/safran/safran.py          | 20 +++++-
 .../scm_models_data/safran/safran_variable.py | 17 ++++-
 ...tive_change_in_maxima_at_fixed_altitude.py | 52 --------------
 .../spatial trends/__init__.py                |  0
 .../{ => spatial trends}/main_result.py       | 20 +++---
 ...tive_change_in_maxima_at_fixed_altitude.py | 70 +++++++++++++++++++
 .../plot_contrasting_trend_curves.py          |  0
 8 files changed, 129 insertions(+), 71 deletions(-)
 delete mode 100644 projects/contrasting_trends_in_snow_loads/main_spatial_relative_change_in_maxima_at_fixed_altitude.py
 create mode 100644 projects/contrasting_trends_in_snow_loads/spatial trends/__init__.py
 rename projects/contrasting_trends_in_snow_loads/{ => spatial trends}/main_result.py (88%)
 create mode 100644 projects/contrasting_trends_in_snow_loads/spatial trends/main_spatial_relative_change_in_maxima_at_fixed_altitude.py
 rename projects/contrasting_trends_in_snow_loads/{ => spatial trends}/plot_contrasting_trend_curves.py (100%)

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 805de9df..234c12dd 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
@@ -200,8 +200,8 @@ class AbstractStudy(object):
     @property
     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()])
+        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
 
@@ -211,6 +211,10 @@ class AbstractStudy(object):
     def observations_annual_maxima(self) -> AnnualMaxima:
         return AnnualMaxima(df_maxima_gev=pd.DataFrame(self.year_to_annual_maxima, index=self.study_massif_names))
 
+    @cached_property
+    def observations_annual_mean(self) -> pd.DataFrame:
+        return pd.DataFrame(self.year_to_annual_mean, index=self.study_massif_names)
+
     def annual_maxima_and_years(self, massif_name) -> Tuple[np.ndarray, np.ndarray]:
         df = self.observations_annual_maxima.df_maxima_gev
         return df.loc[massif_name].values, np.array(df.columns)
@@ -223,6 +227,14 @@ class AbstractStudy(object):
             year_to_annual_maxima[year] = time_serie.max(axis=0)
         return year_to_annual_maxima
 
+    @cached_property
+    def year_to_annual_mean(self) -> OrderedDict:
+        # Map each year to an array of size nb_massif
+        year_to_annual_mean = OrderedDict()
+        for year, time_serie in self._year_to_max_daily_time_serie.items():
+            year_to_annual_mean[year] = time_serie.mean(axis=0)
+        return year_to_annual_mean
+
     @cached_property
     def year_to_annual_maxima_index(self) -> OrderedDict:
         # Map each year to an array of size nb_massif
@@ -238,7 +250,7 @@ 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()
@@ -390,7 +402,8 @@ class AbstractStudy(object):
 
     @property
     def missing_massif_name(self):
-        return set(self.all_massif_names(self.reanalysis_path, self.dbf_filename)) - set(self.altitude_to_massif_names[self.altitude])
+        return set(self.all_massif_names(self.reanalysis_path, self.dbf_filename)) - set(
+            self.altitude_to_massif_names[self.altitude])
 
     @property
     def column_mask(self):
diff --git a/extreme_data/meteo_france_data/scm_models_data/safran/safran.py b/extreme_data/meteo_france_data/scm_models_data/safran/safran.py
index 7512559d..e0fad5e0 100644
--- a/extreme_data/meteo_france_data/scm_models_data/safran/safran.py
+++ b/extreme_data/meteo_france_data/scm_models_data/safran/safran.py
@@ -8,14 +8,16 @@ from extreme_data.meteo_france_data.scm_models_data.abstract_variable import Abs
 from extreme_data.meteo_france_data.scm_models_data.safran.cumulated_study import CumulatedStudy
 from extreme_data.meteo_france_data.scm_models_data.safran.safran_variable import SafranSnowfallVariable, \
     SafranRainfallVariable, SafranTemperatureVariable, SafranTotalPrecipVariable, \
-    SafranNormalizedPrecipitationRateOnWetDaysVariable
+    SafranNormalizedPrecipitationRateOnWetDaysVariable, SafranNormalizedPrecipitationRateVariable
 
 
 class Safran(AbstractStudy):
 
     def __init__(self, variable_class: type, *args, **kwargs):
         assert variable_class in [SafranSnowfallVariable, SafranRainfallVariable, SafranTemperatureVariable,
-                                  SafranTotalPrecipVariable, SafranNormalizedPrecipitationRateOnWetDaysVariable]
+                                  SafranTotalPrecipVariable,
+                                  SafranNormalizedPrecipitationRateVariable,
+                                  SafranNormalizedPrecipitationRateOnWetDaysVariable]
         super().__init__(variable_class, *args, **kwargs)
         self.model_name = 'Safran'
 
@@ -87,6 +89,20 @@ class SafranRainfall7Days(SafranRainfall):
         super().__init__(nb_consecutive_days=7, **kwargs)
 
 
+
+class SafranNormalizedPreciptationRate(CumulatedStudy, Safran):
+
+    def __init__(self, **kwargs):
+        super().__init__(SafranNormalizedPrecipitationRateVariable, **kwargs)
+
+    def load_variable_array(self, dataset):
+        return [np.array(dataset.variables[k]) for k in self.load_keyword()]
+
+    def instantiate_variable_object(self, variable_array) -> AbstractVariable:
+        variable_array_temperature, variable_array_snowfall, variable_array_rainfall = variable_array
+        return self.variable_class(variable_array_temperature,
+                                   variable_array_snowfall, variable_array_rainfall, self.nb_consecutive_days)
+
 class SafranNormalizedPreciptationRateOnWetDays(CumulatedStudy, Safran):
 
     def __init__(self, **kwargs):
diff --git a/extreme_data/meteo_france_data/scm_models_data/safran/safran_variable.py b/extreme_data/meteo_france_data/scm_models_data/safran/safran_variable.py
index 6c279cd3..66dcc8f5 100644
--- a/extreme_data/meteo_france_data/scm_models_data/safran/safran_variable.py
+++ b/extreme_data/meteo_france_data/scm_models_data/safran/safran_variable.py
@@ -77,6 +77,7 @@ class SafranRainfallVariable(SafranSnowfallVariable):
 
 
 class SafranTotalPrecipVariable(AbstractVariable):
+    NAME = 'Precipitation'
 
     def __init__(self, snow_variable_array, rain_variable_array, nb_consecutive_days):
         super().__init__(None)
@@ -94,7 +95,10 @@ class SafranTotalPrecipVariable(AbstractVariable):
         return self._daily_time_serie_array
 
 
-class SafranNormalizedPrecipitationRateOnWetDaysVariable(AbstractVariable):
+class SafranNormalizedPrecipitationRateVariable(AbstractVariable):
+    NAME = 'Normalized Precip'
+
+
 
     def __init__(self, temperature_variable_array, snow_variable_array, rain_variable_array, nb_consecutive_days):
         super().__init__(None)
@@ -103,8 +107,6 @@ class SafranNormalizedPrecipitationRateOnWetDaysVariable(AbstractVariable):
         beta = 0.06
         self._daily_time_serie_array = np.exp(-beta * temperature.daily_time_serie_array) \
                                        * total_precipitation.daily_time_serie_array
-        mask_for_nan_values = total_precipitation.daily_time_serie_array < 0.01
-        self._daily_time_serie_array[mask_for_nan_values] = np.nan
 
     @classmethod
     def keyword(cls):
@@ -115,6 +117,15 @@ class SafranNormalizedPrecipitationRateOnWetDaysVariable(AbstractVariable):
         return self._daily_time_serie_array
 
 
+class SafranNormalizedPrecipitationRateOnWetDaysVariable(SafranNormalizedPrecipitationRateVariable):
+
+    def __init__(self, temperature_variable_array, snow_variable_array, rain_variable_array, nb_consecutive_days):
+        super().__init__(temperature_variable_array, snow_variable_array, rain_variable_array, nb_consecutive_days)
+        total_precipitation = SafranTotalPrecipVariable(snow_variable_array, rain_variable_array, nb_consecutive_days)
+        mask_for_nan_values = total_precipitation.daily_time_serie_array < 0.01
+        self._daily_time_serie_array[mask_for_nan_values] = np.nan
+
+
 class SafranTemperatureVariable(AbstractVariable):
     NAME = 'Temperature'
     UNIT = 'Celsius Degrees'
diff --git a/projects/contrasting_trends_in_snow_loads/main_spatial_relative_change_in_maxima_at_fixed_altitude.py b/projects/contrasting_trends_in_snow_loads/main_spatial_relative_change_in_maxima_at_fixed_altitude.py
deleted file mode 100644
index 91260c0b..00000000
--- a/projects/contrasting_trends_in_snow_loads/main_spatial_relative_change_in_maxima_at_fixed_altitude.py
+++ /dev/null
@@ -1,52 +0,0 @@
-from extreme_data.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoad3Days, \
-    CrocusSnowLoadTotal
-from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall, SafranRainfall, SafranPrecipitation
-from extreme_data.meteo_france_data.scm_models_data.visualization.study_visualizer import \
-    StudyVisualizer
-import matplotlib.pyplot as plt
-
-
-def density_wrt_altitude():
-    save_to_file = True
-    study_class = [SafranPrecipitation, SafranRainfall, SafranSnowfall, CrocusSnowLoad3Days, CrocusSnowLoadTotal][-2]
-    altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][::-1]
-
-    for altitude in altitudes:
-
-        ax = plt.gca()
-        study = study_class(altitude=altitude)
-        # study = study_class(altitude=altitude, nb_consecutive_days=3)
-        massif_name_to_value = {}
-        for massif_name in study.study_massif_names:
-            s = study.observations_annual_maxima.df_maxima_gev.loc[massif_name]
-            year_limit = 1987
-            df_before, df_after = s.loc[:year_limit], s.loc[year_limit + 1:]
-            df_before, df_after = df_before.mean(), df_after.mean()
-            # df_before, df_after = df_before.median(), df_after.median()
-            relative_change_value = 100 * (df_after - df_before) / df_before
-            massif_name_to_value[massif_name] = relative_change_value
-        print(massif_name_to_value)
-        # Plot
-        # massif_name_to_value = {m: i for i, m in enumerate(study.study_massif_names)}
-        max_values = max([abs(e) for e in massif_name_to_value.values()]) + 5
-        print(max_values)
-        variable_name = study.variable_name
-        study.visualize_study(ax=ax, massif_name_to_value=massif_name_to_value,
-                              vmin=-max_values, vmax=max_values,
-                              add_colorbar=True,
-                              replace_blue_by_white=False,
-                              show=False,
-                              label='Relative changes in mean annual maxima\n'
-                                    'of {}\n between 1958-1987 and 1988-2017\n  at {}m (%)\n'
-                                    ''.format(variable_name, study.altitude)
-                              )
-        study_visualizer = StudyVisualizer(study, save_to_file=save_to_file)
-        study_visualizer.plot_name = 'relative_changes_in_maxima'
-        study_visualizer.show_or_save_to_file()
-        ax.clear()
-        plt.close()
-
-
-if __name__ == '__main__':
-    density_wrt_altitude()
-    # test()
diff --git a/projects/contrasting_trends_in_snow_loads/spatial trends/__init__.py b/projects/contrasting_trends_in_snow_loads/spatial trends/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/projects/contrasting_trends_in_snow_loads/main_result.py b/projects/contrasting_trends_in_snow_loads/spatial trends/main_result.py
similarity index 88%
rename from projects/contrasting_trends_in_snow_loads/main_result.py
rename to projects/contrasting_trends_in_snow_loads/spatial trends/main_result.py
index ec17b4e1..9df6f230 100644
--- a/projects/contrasting_trends_in_snow_loads/main_result.py
+++ b/projects/contrasting_trends_in_snow_loads/spatial trends/main_result.py	
@@ -2,13 +2,13 @@ from multiprocessing.pool import Pool
 
 import matplotlib as mpl
 
-from extreme_data.meteo_france_data.scm_models_data.utils import SeasonForTheMaxima
-from extreme_trend.visualizers.utils import load_altitude_to_visualizer
-
 mpl.use('Agg')
 mpl.rcParams['text.usetex'] = True
 mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
 
+from extreme_trend.visualizers.utils import load_altitude_to_visualizer
+from projects.exceeding_snow_loads.section_results.plot_trend_curves import plot_trend_map
+
 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 SafranPrecipitation3Days, \
     SafranPrecipitation1Day, SafranPrecipitation5Days, SafranPrecipitation7Days, SafranSnowfall1Day, \
@@ -19,11 +19,8 @@ from extreme_data.meteo_france_data.scm_models_data.crocus.crocus import CrocusS
     CrocusSnowLoad5Days, CrocusSnowLoad7Days, CrocusSnowLoad1Day
 from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import \
     ConfidenceIntervalMethodFromExtremes
-from projects.contrasting_trends_in_snow_loads.plot_contrasting_trend_curves import plot_contrasting_trend_curves, \
-    plot_contrasting_trend_curves_massif
 from projects.exceeding_snow_loads.section_results.main_result_trends_and_return_levels import \
     compute_minimized_aic
-from root_utils import NB_CORES
 
 
 def intermediate_result(altitudes, massif_names=None,
@@ -58,7 +55,9 @@ def intermediate_result(altitudes, massif_names=None,
 
     # Plots
     # plot_contrasting_trend_curves(altitude_to_visualizer, all_regions=True)
-    plot_contrasting_trend_curves_massif(altitude_to_visualizer, all_regions=True)
+    # plot_contrasting_trend_curves_massif(altitude_to_visualizer, all_regions=True)
+    # plot_trend_curves(altitude_to_visualizer)
+    plot_trend_map(altitude_to_visualizer)
 
 
 def major_result():
@@ -68,15 +67,16 @@ def major_result():
     model_subsets_for_uncertainty = None
     # altitudes = paper_altitudes
     # altitudes = paper_altitudes
-    altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][:]
+    # altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][:]
+    altitudes = [1200, 1500, 1800, 2100, 2400, 2700][:]
     snow_load_classes = [CrocusSnowLoad1Day, CrocusSnowLoad3Days, CrocusSnowLoad5Days, CrocusSnowLoad7Days][:]
     precipitation_classes = [SafranPrecipitation1Day, SafranPrecipitation3Days, SafranPrecipitation5Days,
                              SafranPrecipitation7Days][:]
     snowfall_classes = [SafranSnowfall1Day, SafranSnowfall3Days, SafranSnowfall5Days, SafranSnowfall7Days]
     rainfall_classes = [SafranRainfall1Day, SafranRainfall3Days, SafranRainfall5Days, SafranRainfall7Days]
     study_classes = precipitation_classes + snow_load_classes
-    # study_classes = snowfall_classes + rainfall_classes
-    for study_class in snowfall_classes[:1]:
+    study_classes = snowfall_classes + rainfall_classes
+    for study_class in precipitation_classes:
         intermediate_result(altitudes, massif_names, model_subsets_for_uncertainty,
                             uncertainty_methods, study_class, multiprocessing=True)
 
diff --git a/projects/contrasting_trends_in_snow_loads/spatial trends/main_spatial_relative_change_in_maxima_at_fixed_altitude.py b/projects/contrasting_trends_in_snow_loads/spatial trends/main_spatial_relative_change_in_maxima_at_fixed_altitude.py
new file mode 100644
index 00000000..60a7a867
--- /dev/null
+++ b/projects/contrasting_trends_in_snow_loads/spatial trends/main_spatial_relative_change_in_maxima_at_fixed_altitude.py	
@@ -0,0 +1,70 @@
+import pandas as pd
+from extreme_data.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoad3Days, \
+    CrocusSnowLoadTotal
+from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall, SafranRainfall, \
+    SafranPrecipitation, SafranPrecipitation1Day, SafranSnowfall1Day, SafranTemperature, \
+    SafranNormalizedPreciptationRateOnWetDays, SafranNormalizedPreciptationRate
+from extreme_data.meteo_france_data.scm_models_data.safran.safran_variable import \
+    SafranNormalizedPrecipitationRateOnWetDaysVariable
+from extreme_data.meteo_france_data.scm_models_data.utils import SeasonForTheMaxima
+from extreme_data.meteo_france_data.scm_models_data.visualization.study_visualizer import \
+    StudyVisualizer
+import matplotlib.pyplot as plt
+
+
+def relative_change_in_maxima_wrt_altitude():
+    save_to_file = True
+    study_class = [SafranSnowfall1Day, SafranPrecipitation1Day, SafranPrecipitation, SafranRainfall,
+                   SafranSnowfall, CrocusSnowLoad3Days, CrocusSnowLoadTotal,
+                   SafranTemperature, SafranNormalizedPreciptationRateOnWetDays,
+                   SafranNormalizedPreciptationRate][-1]
+    altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][::-1]
+    relative = False
+
+    for altitude in altitudes:
+
+        ax = plt.gca()
+        study = study_class(altitude=altitude, season=SeasonForTheMaxima.winter_extended)
+        # study = study_class(altitude=altitude, nb_consecutive_days=3)
+        massif_name_to_value = {}
+        for massif_name in study.study_massif_names:
+            if study_class is SafranTemperature:
+                s = study.observations_annual_mean.loc[massif_name]
+            else:
+                s = study.observations_annual_maxima.df_maxima_gev.loc[massif_name]
+            year_limit = 1989
+            df_before, df_after = s.loc[:year_limit], s.loc[year_limit + 1:]
+            df_before, df_after = df_before.mean(), df_after.mean()
+            # df_before, df_after = df_before.median(), df_after.median()
+            if relative:
+                change_value = 100 * (df_after - df_before) / df_before
+            else:
+                change_value = (df_after - df_before)
+            massif_name_to_value[massif_name] = change_value
+        print(massif_name_to_value)
+        # Plot
+        # massif_name_to_value = {m: i for i, m in enumerate(study.study_massif_names)}
+        max_values = max([abs(e) for e in massif_name_to_value.values()]) * 1.05
+        print(max_values)
+        variable_name = study.variable_name
+        prefix = 'Relative' if relative else ''
+        str_season = str(study.season).split('.')[-1]
+        study.visualize_study(ax=ax, massif_name_to_value=massif_name_to_value,
+                              vmin=-max_values, vmax=max_values,
+                              add_colorbar=True,
+                              replace_blue_by_white=False,
+                              show=False,
+                              label='{} changes in mean {} maxima\n'
+                                    'of {}\n between 1959-1989 and 1990-2019\n  at {}m (%)\n'
+                                    ''.format(prefix, str_season, variable_name, study.altitude)
+                              )
+        study_visualizer = StudyVisualizer(study, save_to_file=save_to_file)
+        study_visualizer.plot_name = '{}_changes_in_maxima'.format(prefix)
+        study_visualizer.show_or_save_to_file()
+        ax.clear()
+        plt.close()
+
+
+if __name__ == '__main__':
+    relative_change_in_maxima_wrt_altitude()
+    # test()
diff --git a/projects/contrasting_trends_in_snow_loads/plot_contrasting_trend_curves.py b/projects/contrasting_trends_in_snow_loads/spatial trends/plot_contrasting_trend_curves.py
similarity index 100%
rename from projects/contrasting_trends_in_snow_loads/plot_contrasting_trend_curves.py
rename to projects/contrasting_trends_in_snow_loads/spatial trends/plot_contrasting_trend_curves.py
-- 
GitLab