From a5ff99e17d3e602b3c16f493bc5c6b277068d6f2 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Tue, 7 Apr 2020 18:15:49 +0200
Subject: [PATCH] [contrasting project] add
 SafranNormalizedPreciptationRateOnWetDays.

---
 .../scm_models_data/safran/safran.py          | 34 +++++++++++++++++--
 .../scm_models_data/safran/safran_variable.py | 21 ++++++++++++
 .../figure3_daily snowfall fraction.py        |  4 +--
 .../test_meteo_france_data/test_SCM_study.py  | 11 ++++--
 4 files changed, 63 insertions(+), 7 deletions(-)

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 8eecfc81..7512559d 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
@@ -1,3 +1,5 @@
+from collections import OrderedDict
+
 import numpy as np
 
 from extreme_data.meteo_france_data.scm_models_data.abstract_extended_study import AbstractExtendedStudy
@@ -5,14 +7,15 @@ from extreme_data.meteo_france_data.scm_models_data.abstract_study import Abstra
 from extreme_data.meteo_france_data.scm_models_data.abstract_variable import AbstractVariable
 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
+    SafranRainfallVariable, SafranTemperatureVariable, SafranTotalPrecipVariable, \
+    SafranNormalizedPrecipitationRateOnWetDaysVariable
 
 
 class Safran(AbstractStudy):
 
     def __init__(self, variable_class: type, *args, **kwargs):
         assert variable_class in [SafranSnowfallVariable, SafranRainfallVariable, SafranTemperatureVariable,
-                                  SafranTotalPrecipVariable]
+                                  SafranTotalPrecipVariable, SafranNormalizedPrecipitationRateOnWetDaysVariable]
         super().__init__(variable_class, *args, **kwargs)
         self.model_name = 'Safran'
 
@@ -84,6 +87,31 @@ class SafranRainfall7Days(SafranRainfall):
         super().__init__(nb_consecutive_days=7, **kwargs)
 
 
+class SafranNormalizedPreciptationRateOnWetDays(CumulatedStudy, Safran):
+
+    def __init__(self, **kwargs):
+        super().__init__(SafranNormalizedPrecipitationRateOnWetDaysVariable, **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)
+
+    @property
+    def _year_to_daily_time_serie_array(self) -> OrderedDict:
+        # Filter and keep only values different than np.nan
+        year_to_time_series = super()._year_to_daily_time_serie_array
+        updated_year_to_time_series = OrderedDict()
+        for year, time_serie in year_to_time_series.items():
+            time_serie_without_nan = time_serie[~np.isnan(time_serie)]
+            assert not np.isnan(time_serie_without_nan).any()
+            updated_year_to_time_series[year] = time_serie_without_nan
+        return updated_year_to_time_series
+
+
 class SafranPrecipitation(CumulatedStudy, Safran):
 
     def __init__(self, **kwargs):
@@ -138,7 +166,7 @@ if __name__ == '__main__':
     altitude = 900
     year_min = 1959
     year_max = 2000
-    study = SafranRainfall1Day(altitude, year_min=year_min, year_max=year_max)
+    study = SafranNormalizedPreciptationRateOnWetDays(altitude=altitude, year_min=year_min, year_max=year_max)
     d = study.year_to_dataset_ordered_dict[1959]
     print(d.keywords)
     print(d.variables.keys())
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 aa7f78d3..6c279cd3 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
@@ -94,6 +94,27 @@ class SafranTotalPrecipVariable(AbstractVariable):
         return self._daily_time_serie_array
 
 
+class SafranNormalizedPrecipitationRateOnWetDaysVariable(AbstractVariable):
+
+    def __init__(self, temperature_variable_array, snow_variable_array, rain_variable_array, nb_consecutive_days):
+        super().__init__(None)
+        temperature = SafranTemperatureVariable(temperature_variable_array)
+        total_precipitation = SafranTotalPrecipVariable(snow_variable_array, rain_variable_array, nb_consecutive_days)
+        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):
+        return [SafranTemperatureVariable.keyword(), SafranSnowfallVariable.keyword(), SafranRainfallVariable.keyword()]
+
+    @property
+    def daily_time_serie_array(self) -> np.ndarray:
+        return self._daily_time_serie_array
+
+
 class SafranTemperatureVariable(AbstractVariable):
     NAME = 'Temperature'
     UNIT = 'Celsius Degrees'
diff --git a/projects/contrasting_trends_in_snow_loads/gorman_figures/figure3_daily snowfall fraction.py b/projects/contrasting_trends_in_snow_loads/gorman_figures/figure3_daily snowfall fraction.py
index 6a1b5cf1..0b06ee9a 100644
--- a/projects/contrasting_trends_in_snow_loads/gorman_figures/figure3_daily snowfall fraction.py	
+++ b/projects/contrasting_trends_in_snow_loads/gorman_figures/figure3_daily snowfall fraction.py	
@@ -93,6 +93,6 @@ def daily_snowfall_fraction_wrt_time():
 
 
 if __name__ == '__main__':
-    daily_snowfall_fraction_wrt_altitude(fast=False)
+    # daily_snowfall_fraction_wrt_altitude(fast=False)
     # daily_snowfall_fraction_wrt_time()
-    # ratio_of_past_and_recent_daily_snowfall_fraction_wrt_altitude()
+    ratio_of_past_and_recent_daily_snowfall_fraction_wrt_altitude()
diff --git a/test/test_extreme_data/test_meteo_france_data/test_SCM_study.py b/test/test_extreme_data/test_meteo_france_data/test_SCM_study.py
index 16cfbe9c..b9d999d3 100644
--- a/test/test_extreme_data/test_meteo_france_data/test_SCM_study.py
+++ b/test/test_extreme_data/test_meteo_france_data/test_SCM_study.py
@@ -1,4 +1,5 @@
 import os.path as op
+import numpy as np
 import unittest
 from random import sample
 
@@ -7,7 +8,7 @@ import pandas as pd
 from extreme_data.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoad3Days
 from extreme_data.meteo_france_data.scm_models_data.safran.cumulated_study import NB_DAYS
 from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall, SafranTemperature, \
-    SafranPrecipitation, SafranSnowfall3Days, SafranRainfall3Days
+    SafranPrecipitation, SafranSnowfall3Days, SafranRainfall3Days, SafranNormalizedPreciptationRateOnWetDays
 from extreme_data.meteo_france_data.scm_models_data.utils import SeasonForTheMaxima
 from extreme_data.meteo_france_data.scm_models_data.visualization.main_study_visualizer import \
     study_iterator_global, SCM_STUDIES, ALL_ALTITUDES
@@ -55,6 +56,13 @@ class TestSCMAllStudy(unittest.TestCase):
             study_class(altitude=altitude, year_min=year_min, year_max=year_max)
 
 
+class TestSCMSafranNormalizedPrecipitationRateOnWetDays(unittest.TestCase):
+
+    def test_annual_maxima(self):
+        study = SafranNormalizedPreciptationRateOnWetDays(year_max=1960)
+        self.assertFalse(np.isnan(study.year_to_annual_maxima[1959]).any())
+
+
 class TestSCMStudy(unittest.TestCase):
 
     def setUp(self) -> None:
@@ -88,7 +96,6 @@ class TestSCMSafranSnowfall(TestSCMStudy):
         self.assertEqual(all_daily_series.shape[1], 23)
         self.assertEqual(all_daily_series.shape[0], 22280)
 
-
 class TestSCMPrecipitation(TestSCMStudy):
 
     def setUp(self) -> None:
-- 
GitLab