[contrasting project] add normalized precipitation. add folder spatial trends

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):

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):

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'

projects/contrasting_trends_in_snow_loads/main_result.py → 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)
 

projects/contrasting_trends_in_snow_loads/main_spatial_relative_change_in_maxima_at_fixed_altitude.py → projects/contrasting_trends_in_snow_loads/spatial trends/main_spatial_relative_change_in_maxima_at_fixed_altitude.py

+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
+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 density_wrt_altitude():
+def relative_change_in_maxima_wrt_altitude():
     save_to_file = True
-    study_class = [SafranPrecipitation, SafranRainfall, SafranSnowfall, CrocusSnowLoad3Days, CrocusSnowLoadTotal][-2]
-    altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][::-1]
+    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)
+        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:
-            s = study.observations_annual_maxima.df_maxima_gev.loc[massif_name]
-            year_limit = 1987
+            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()
-            relative_change_value = 100 * (df_after - df_before) / df_before
-            massif_name_to_value[massif_name] = relative_change_value
+            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()]) + 5
+        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='Relative changes in mean annual maxima\n'
-                                    'of {}\n between 1958-1987 and 1988-2017\n  at {}m (%)\n'
-                                    ''.format(variable_name, study.altitude)
+                              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 = 'relative_changes_in_maxima'
+        study_visualizer.plot_name = '{}_changes_in_maxima'.format(prefix)
         study_visualizer.show_or_save_to_file()
         ax.clear()
         plt.close()
 
 
 if __name__ == '__main__':
-    density_wrt_altitude()
+    relative_change_in_maxima_wrt_altitude()
     # test()