[paper 3] add main difference between reanalysis and simulations. move adamont data.

adamont/single_simulation.py

-from cached_property import cached_property
-from netCDF4._netCDF4 import Dataset
-
-
-class SingleSimulation(object):
-
-    def __init__(self, nc_path):
-        self.nc_path = nc_path
-
-    @cached_property
-    def dataset(self):
-        return Dataset(self.nc_path)
-
-
-
-    def massif_name_and_altitude_to_return_level(self):
-        return {}
-
-

adamont/ensemble_simulation.py → experiment/meteo_france_data/adamont_data/ensemble_simulation.py

@@ -2,9 +2,10 @@ import os
 from typing import List
 import os.path as op
 
+import numpy as np
 from cached_property import cached_property
 
-from adamont.single_simulation import SingleSimulation
+from experiment.meteo_france_data.adamont_data.single_simulation import SingleSimulation
 
 ADAMONT_PATH = r"/home/erwan/Documents/projects/spatiotemporalextremes/local/spatio_temporal_datasets/ADAMONT"
 
@@ -25,14 +26,19 @@ class EnsembleSimulation(object):
         assert first_winter_required_for_histo <= 2004
 
         # Load simulations
-        self.simulations = [SingleSimulation(op.join(ADAMONT_PATH, nc_file)) for nc_file in self.nc_files]
-
-    """Ce problème affecte toutes lessimulations HISTORIQUE CORDEX réalisées en utilisant le forçage CNRM-CM5: CCLM4-8-17: ALADIN53 et RCA4"""
+        # todo: so far i am using one ensemble member
+        self.simulations = [SingleSimulation(nc_path, self.parameter,
+                                             self.first_winter_required_for_histo,
+                                             self.last_year_for_histo) for nc_path in self.nc_paths][:1]
 
     @cached_property
     def simulations_path(self):
         return op.join(ADAMONT_PATH, self.parameter, self.scenario)
 
+    @cached_property
+    def nc_paths(self):
+        return [op.join(ADAMONT_PATH, self.parameter, self.scenario, nc_file) for nc_file in self.nc_files]
+
     @cached_property
     def nc_files(self) -> List[str]:
         nc_files = []
@@ -57,9 +63,31 @@ class EnsembleSimulation(object):
     def massif_name_and_altitude_to_mean_return_level(self):
         return {}
 
+    @property
+    def first_simulation(self):
+        return self.simulations[0]
+
+    @property
+    def massif_name_and_altitude(self):
+        pass
+
+    @cached_property
+    def massif_name_and_altitude_to_mean_average_annual_maxima(self):
+        d = {}
+        for m, a in self.first_simulation.massif_name_and_altitude_to_average_maxima.keys():
+            d[(m, a)] = np.mean([s.massif_name_and_altitude_to_average_maxima[(m, a)] for s in self.simulations])
+        return d
 
 
 if __name__ == '__main__':
-    ensemble = EnsembleSimulation(first_winter_required_for_histo=1985)
+    # np.array(d.variables['SNOWSWE'])
+    ensemble = EnsembleSimulation(first_winter_required_for_histo=1958)
     print(len(ensemble.simulations))
     print(ensemble.simulations_names)
+    s = ensemble.first_simulation
+    d = s.dataset
+    # print(s.massif_name_and_altitude_to_annual_maxima_time_series)
+    # print(s.massif_name_and_altitude_to_average_maxima)
+    print(ensemble.massif_name_and_altitude_to_mean_average_annual_maxima)
+    print(s.years)
+    # print(d)

experiment/meteo_france_data/adamont_data/single_simulation.py

+from datetime import datetime
+
+import numpy as np
+import os.path as op
+from cached_property import cached_property
+from netCDF4._netCDF4 import Dataset
+from datetime import timedelta
+
+from experiment.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
+
+
+class SingleSimulation(object):
+
+    def __init__(self, nc_path, parameter, fist_year, last_year):
+        self.fist_year = fist_year
+        self.last_year = last_year
+        self.parameter = parameter
+        self.nc_path = nc_path
+
+    @cached_property
+    def dataset(self):
+        return Dataset(self.nc_path)
+
+    @cached_property
+    def winter_year(self):
+        start = datetime(year=1900, month=1, day=1, hour=0, minute=0, second=0)
+        seconds_after_start = np.array(self.dataset.variables['TIME'])
+        dates = [start + timedelta(seconds=s) for s in seconds_after_start]
+        winter_year = [date.year - 1 if date.month < 8 else date.year for date in dates]
+        return np.array(winter_year)
+
+    @cached_property
+    def years(self):
+        return sorted([year for year in set(self.winter_year) if self.fist_year <= year <= self.last_year])
+
+    def massif_name_and_altitude_to_return_level(self):
+        return {}
+
+    @property
+    def massif_number_to_massif_name(self):
+        # from adamont_data metadata
+        s = """1	Chablais
+        2	Aravis
+        3	Mont-Blanc
+        4	Bauges
+        5	Beaufortain
+        6	Haute-Tarentaise
+        7	Chartreuse
+        8	Belledonne
+        9	Maurienne
+        10	Vanoise
+        11	Haute-Maurienne
+        12	Grandes-Rousses
+        13	Thabor
+        14	Vercors
+        15	Oisans
+        16	Pelvoux
+        17	Queyras
+        18	Devoluy
+        19	Champsaur
+        20	Parpaillon
+        21	Ubaye
+        22	Haut_Var-Haut_Verdon
+        23	Mercantour"""
+        l = s.split('\n')
+        return dict([e.split() for e in l])
+
+    @cached_property
+    def massif_name_and_altitude_to_annual_maxima_time_series(self):
+        all_values = np.array(self.dataset.variables[self.parameter])
+        zs_list = [int(e) for e in np.array(self.dataset.variables['ZS'])]
+        massif_number_list = np.array(self.dataset.variables['MASSIF_NUMBER'])
+        massif_name_list = [self.massif_number_to_massif_name[str(n)] for n in massif_number_list]
+        d = {}
+        for year in self.years:
+            indexes = np.where(self.winter_year == year)[0]
+            winter_values = all_values[indexes, 0, :]
+            assert len(winter_values) in [365, 366]
+            for time_serie, zs, massif_name in zip(winter_values.transpose(), zs_list, massif_name_list):
+                # print(zs, massif_name, len(time_serie))
+                d[(massif_name, zs)] = time_serie
+        return d
+
+    @cached_property
+    def massif_name_and_altitude_to_average_maxima(self):
+        return {t: np.mean(s) for t, s in self.massif_name_and_altitude_to_annual_maxima_time_series.items()}

experiment/paper 2 - contrasting/main_spatial_relative_change_in_maxima_at_fixed_altitude.py

@@ -38,8 +38,8 @@ def density_wrt_altitude():
             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()
+            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)

experiment/paper 3/main_difference_between_reanalysis_and_simulations.py

+from experiment.meteo_france_data.adamont_data.ensemble_simulation import EnsembleSimulation
+from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoad3Days, \
+    CrocusSweTotal
+from experiment.meteo_france_data.scm_models_data.visualization.study_visualization.study_visualizer import \
+    StudyVisualizer
+import matplotlib.pyplot as plt
+
+
+def test():
+    study = CrocusSnowLoad3Days(altitude=1200)
+    study_visualizer = StudyVisualizer(study)
+    study_visualizer.visualize_max_graphs_poster('Queyras', altitude='noope', snow_abbreviation="ok", color='red')
+    plt.show()
+
+
+def density_wrt_altitude():
+    save_to_file = True
+    study_class = CrocusSweTotal
+    altitudes = [600, 900, 1200, 1500, 1800, 2100, 2400, 2700][::-1]
+
+    ensemble = EnsembleSimulation(first_winter_required_for_histo=1958)
+
+    # None means relative change, true means reanalysis values, false means simulations values
+    # mode = False
+
+    for mode in [None, True, False]:
+        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 = {}
+            df = study.observations_annual_maxima.df_maxima_gev.iloc[:2004]
+            for massif_name in study.study_massif_names:
+                s_reanalysis = df.loc[massif_name].mean()
+                if (massif_name, altitude) in ensemble.massif_name_and_altitude_to_mean_average_annual_maxima:
+                    s_simulation = ensemble.massif_name_and_altitude_to_mean_average_annual_maxima[(massif_name, altitude)]
+                    relative_change_value = 100 * (s_simulation - s_reanalysis) / s_reanalysis
+                    if mode == None:
+                        value = relative_change_value
+                    else:
+                        value = s_reanalysis if mode else s_simulation
+                    massif_name_to_value[massif_name] = 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
+            label_relative_change = 'Relative changes in mean annual maxima\n' \
+                        'of {}\n between reanalysis and historical simulation\n' \
+                        'for 1958-2004  at {}m (%)\n' \
+                        ''.format(variable_name, study.altitude)
+            if mode is None:
+                label = label_relative_change
+            else:
+                label = 'Mean max SWE for ' + ('reanalysis' if mode else 'simulation')
+
+            vmin = -max_values if mode is None else 0.1
+            study.visualize_study(ax=ax, massif_name_to_value=massif_name_to_value,
+                                  vmin=vmin, vmax=max_values,
+                                  add_colorbar=True,
+                                  replace_blue_by_white=False,
+                                  show=False,
+                                  label=label
+                                  )
+            study_visualizer = StudyVisualizer(study, save_to_file=save_to_file)
+            study_visualizer.plot_name = 'relative_changes_in_maxima' if mode is None else label
+            study_visualizer.show_or_save_to_file()
+            ax.clear()
+            plt.close()
+
+
+if __name__ == '__main__':
+    density_wrt_altitude()
+    # test()