diff --git a/adamont/single_simulation.py b/adamont/single_simulation.py deleted file mode 100644 index ab588f32a78f536aafbb0ed205fb522e37577497..0000000000000000000000000000000000000000 --- a/adamont/single_simulation.py +++ /dev/null @@ -1,19 +0,0 @@ -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 {} - - diff --git a/experiment/meteo_france_data/adamont_data/__init__.py b/experiment/meteo_france_data/adamont_data/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/adamont/download_adamont.py b/experiment/meteo_france_data/adamont_data/download_adamont.py similarity index 100% rename from adamont/download_adamont.py rename to experiment/meteo_france_data/adamont_data/download_adamont.py diff --git a/adamont/ensemble_simulation.py b/experiment/meteo_france_data/adamont_data/ensemble_simulation.py similarity index 61% rename from adamont/ensemble_simulation.py rename to experiment/meteo_france_data/adamont_data/ensemble_simulation.py index ea80242c9eb8b12d0209cdd324195da6d7c0bfd9..bdd1b692d5b4a1b669d59cfcd292017525e6745e 100644 --- a/adamont/ensemble_simulation.py +++ b/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) diff --git a/experiment/meteo_france_data/adamont_data/single_simulation.py b/experiment/meteo_france_data/adamont_data/single_simulation.py new file mode 100644 index 0000000000000000000000000000000000000000..0292ed6cd5eb86f33b5a4da22b091fc28eeb9906 --- /dev/null +++ b/experiment/meteo_france_data/adamont_data/single_simulation.py @@ -0,0 +1,86 @@ +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()} diff --git a/experiment/paper 2 - contrasting/main_spatial_relative_change_in_maxima_at_fixed_altitude.py b/experiment/paper 2 - contrasting/main_spatial_relative_change_in_maxima_at_fixed_altitude.py index d98427cc1cda9ab703df1721ee83ab966f49677c..296cbfdb4635784609798bab91e1678526955953 100644 --- a/experiment/paper 2 - contrasting/main_spatial_relative_change_in_maxima_at_fixed_altitude.py +++ b/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) diff --git a/experiment/paper 3/main_difference_between_reanalysis_and_simulations.py b/experiment/paper 3/main_difference_between_reanalysis_and_simulations.py new file mode 100644 index 0000000000000000000000000000000000000000..220143476743f42d53f3659ed7625d00db779a2f --- /dev/null +++ b/experiment/paper 3/main_difference_between_reanalysis_and_simulations.py @@ -0,0 +1,76 @@ +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()