From 5d45d5c3edcbab8e6185b3fe2861098f629fd3c8 Mon Sep 17 00:00:00 2001 From: Le Roux Erwan <erwan.le-roux@irstea.fr> Date: Tue, 10 Nov 2020 16:22:33 +0100 Subject: [PATCH] [snow projection] improve drastically the loading time for the abstract_adamont_study.py based on a simple conversion to numpy array. remove fit_method from plot_map argument. --- .../adamont_data/abstract_adamont_study.py | 18 +++-- .../adamont_data/adamont_studies.py | 22 ++++-- .../scm_models_data/abstract_study.py | 9 ++- .../visualization/plot_utils.py | 2 +- .../visualization/study_visualizer.py | 6 +- ...es_visualizer_for_non_stationary_models.py | 2 +- .../preliminary_analysis.py | 2 +- .../comparison_historical_visualizer.py | 73 ++++++++++++++++--- .../comparison_with_scm/comparison_plot.py | 17 +++-- ...main_comparison_on_beginning_projection.py | 43 ----------- ... => main_comparison_on_rcp_85_extended.py} | 38 +++++----- 11 files changed, 131 insertions(+), 101 deletions(-) delete mode 100644 projects/projected_snowfall/comparison_with_scm/main_comparison_on_beginning_projection.py rename projects/projected_snowfall/comparison_with_scm/{main_comparison_on_historical_period.py => main_comparison_on_rcp_85_extended.py} (72%) diff --git a/extreme_data/meteo_france_data/adamont_data/abstract_adamont_study.py b/extreme_data/meteo_france_data/adamont_data/abstract_adamont_study.py index a292f9cc..43eae3c0 100644 --- a/extreme_data/meteo_france_data/adamont_data/abstract_adamont_study.py +++ b/extreme_data/meteo_france_data/adamont_data/abstract_adamont_study.py @@ -76,7 +76,9 @@ class AbstractAdamontStudy(AbstractStudy): # Load data & year list data_list, data_year_list = [], [] for dataset, real_scenario in zip(self.datasets, self.adamont_real_scenarios): - data_list.extend(dataset.variables[self.variable_class.keyword()]) + data = dataset.variables[self.variable_class.keyword()] + data = np.array(data) + data_list.extend(data) data_year_list.extend(self.winter_years_for_each_time_step(real_scenario, dataset)) # Remark. The last winter year for the HISTO scenario correspond to 2006. # Thus, the last value is not taken into account @@ -87,13 +89,17 @@ class AbstractAdamontStudy(AbstractStudy): for year_data, data in zip(data_year_list, data_list): if self.year_min <= year_data <= self.year_max: year_to_data_list[year_data].append(data) - # Load the variable object - year_to_variable_object = OrderedDict() - for year in self.ordered_years: - variable_array = np.array(year_to_data_list[year]) - year_to_variable_object[year] = self.variable_class(variable_array) + # Load efficiently the variable object + # Multiprocessing is set to False, because this is not a time consuming part + data_list_list = [year_to_data_list[year] for year in self.ordered_years] + year_to_variable_object = self.efficient_variable_loading(self.ordered_years, data_list_list, multiprocessing=False) return year_to_variable_object + def load_variable_object(self, data_list): + variable_array = np.array(data_list) + variable_object = self.variable_class(variable_array) + return variable_object + @cached_property def flat_mask(self): zs_list = [int(e) for e in np.array(self.dataset.variables['ZS'])] diff --git a/extreme_data/meteo_france_data/adamont_data/adamont_studies.py b/extreme_data/meteo_france_data/adamont_data/adamont_studies.py index e8a7d84e..55eea8e9 100644 --- a/extreme_data/meteo_france_data/adamont_data/adamont_studies.py +++ b/extreme_data/meteo_france_data/adamont_data/adamont_studies.py @@ -34,6 +34,10 @@ class AdamontStudies(object): def study(self) -> AbstractStudy: return self.study_list[0] + @property + def nb_ensemble_members(self): + return len(self.gcm_rcm_couples) + # Some plots def show_or_save_to_file(self, plot_name, show=False, no_title=False, tight_layout=None): @@ -42,7 +46,7 @@ class AdamontStudies(object): study_visualizer.show_or_save_to_file(add_classic_title=False, dpi=500, no_title=no_title, tight_layout=tight_layout) - def plot_maxima_time_series(self, massif_names=None, scm_study=None): + def plot_maxima_time_series_adamont(self, massif_names=None, scm_study=None): massif_names = massif_names if massif_names is not None else self.study.all_massif_names() for massif_names in massif_names: self._plot_maxima_time_series(massif_names, scm_study) @@ -75,15 +79,19 @@ class AdamontStudies(object): except KeyError: pass - ax.xaxis.set_ticks([year for year in y if year % 10 == 0]) - ax.grid() + ticks = [year for year in x if year % 10 == 0] + ax.xaxis.set_ticks(ticks) + ax.yaxis.grid() + ax.set_xlim((min(x), max(x))) ax.tick_params(axis='both', which='major', labelsize=13) handles, labels = ax.get_legend_handles_labels() ax.legend(handles[::-1], labels[::-1], ncol=2) - plot_name = 'Annual maxima of {} in {}'.format(ADAMONT_STUDY_CLASS_TO_ABBREVIATION[self.study_class], - massif_name.replace('_', ' ')) - ax.set_ylabel('{} ({})'.format(plot_name, self.study.variable_unit), fontsize=15) - ax.set_xlabel('years', fontsize=15) + plot_name = 'Annual maxima of {} in {} at {} m'.format(ADAMONT_STUDY_CLASS_TO_ABBREVIATION[self.study_class], + massif_name.replace('_', ' '), + self.study.altitude) + fontsize = 13 + ax.set_ylabel('{} ({})'.format(plot_name, self.study.variable_unit), fontsize=fontsize) + ax.set_xlabel('years', fontsize=fontsize) plot_name = 'time series/' + plot_name self.show_or_save_to_file(plot_name=plot_name, show=False, no_title=True, tight_layout=True) ax.clear() 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 1436a908..27a0ff49 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 @@ -369,11 +369,14 @@ class AbstractStudy(object): def year_to_variable_object(self) -> OrderedDict: # Map each year to the variable array path_files, ordered_years = self.ordered_years_and_path_files - if self.multiprocessing: + return self.efficient_variable_loading(ordered_years, path_files, multiprocessing=self.multiprocessing) + + def efficient_variable_loading(self, ordered_years, arguments, multiprocessing): + if multiprocessing: with Pool(NB_CORES) as p: - variables = p.map(self.load_variable_object, path_files) + variables = p.map(self.load_variable_object, arguments) else: - variables = [self.load_variable_object(path_file) for path_file in path_files] + variables = [self.load_variable_object(argument) for argument in arguments] return OrderedDict(zip(ordered_years, variables)) def instantiate_variable_object(self, variable_array) -> AbstractVariable: diff --git a/extreme_data/meteo_france_data/scm_models_data/visualization/plot_utils.py b/extreme_data/meteo_france_data/scm_models_data/visualization/plot_utils.py index 2a1bc186..ecab9092 100644 --- a/extreme_data/meteo_france_data/scm_models_data/visualization/plot_utils.py +++ b/extreme_data/meteo_france_data/scm_models_data/visualization/plot_utils.py @@ -34,7 +34,7 @@ def plot_against_altitude(x_ticks, ax, massif_id, massif_name, values, altitude= ax.fill_between(x_ticks, lower_bound, upper_bound, color=color, alpha=0.2) -def load_plot(cmap, graduation, label, massif_name_to_value, altitude, fit_method, add_x_label=True, +def load_plot(cmap, graduation, label, massif_name_to_value, altitude, add_x_label=True, negative_and_positive_values=True, massif_name_to_text=None, add_colorbar=True, max_abs_change=None, xlabel=None, fontsize_label=10): if max_abs_change is None: diff --git a/extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py b/extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py index 877c78c0..8ce26dec 100644 --- a/extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py +++ b/extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py @@ -728,13 +728,13 @@ class StudyVisualizer(VisualizationParameters): # PLot functions that should be common - def plot_map(self, cmap, fit_method, graduation, label, massif_name_to_value, plot_name, add_x_label=True, + def plot_map(self, cmap, graduation, label, massif_name_to_value, plot_name, add_x_label=True, negative_and_positive_values=True, massif_name_to_text=None, altitude=None, add_colorbar=True, max_abs_change=None, xlabel=None, fontsize_label=10): if altitude is None: altitude = self.study.altitude if len(massif_name_to_value) > 0: - load_plot(cmap, graduation, label, massif_name_to_value, altitude, fitmethod_to_str(fit_method), + load_plot(cmap, graduation, label, massif_name_to_value, altitude, add_x_label=add_x_label, negative_and_positive_values=negative_and_positive_values, massif_name_to_text=massif_name_to_text, add_colorbar=add_colorbar, max_abs_change=max_abs_change, xlabel=xlabel, @@ -751,6 +751,6 @@ class StudyVisualizer(VisualizationParameters): # Regroup the plot by type of plot also plot_name2 = '{}/{}'.format(plot_name.split()[0], plot_name) for plot_name in [plot_name1, plot_name2]: - self.plot_map(cmap, fit_method, graduation, label, massif_name_to_value, plot_name, add_x_label, negative_and_positive_values, + self.plot_map(cmap, graduation, label, massif_name_to_value, plot_name, add_x_label, negative_and_positive_values, massif_name_to_text, ) diff --git a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py index de092623..e174d78c 100644 --- a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py +++ b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py @@ -170,7 +170,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer): negative_and_positive_values = self.moment_names.index(method_name) > 0 # Plot the map - self.plot_map(cmap=plt.cm.coolwarm, fit_method=self.fit_method, graduation=graduation, + self.plot_map(cmap=plt.cm.coolwarm, graduation=graduation, label=ylabel, massif_name_to_value=massif_name_to_value, plot_name=plot_name, add_x_label=True, negative_and_positive_values=negative_and_positive_values, diff --git a/projects/altitude_spatial_model/preliminary_analysis.py b/projects/altitude_spatial_model/preliminary_analysis.py index 9b33c8d5..935e2728 100644 --- a/projects/altitude_spatial_model/preliminary_analysis.py +++ b/projects/altitude_spatial_model/preliminary_analysis.py @@ -107,7 +107,7 @@ class PointwiseGevStudyVisualizer(AltitudesStudies): } if param_name == GevParams.SHAPE: print(massif_name_to_linear_coef) - visualizer.plot_map(cmap=plt.cm.coolwarm, fit_method=self.study.fit_method, + visualizer.plot_map(cmap=plt.cm.coolwarm, graduation=gev_param_name_to_graduation[param_name], label=label, massif_name_to_value=massif_name_to_linear_coef, plot_name=label.replace('/', ' every '), add_x_label=False, diff --git a/projects/projected_snowfall/comparison_with_scm/comparison_historical_visualizer.py b/projects/projected_snowfall/comparison_with_scm/comparison_historical_visualizer.py index ee8c31e3..72beeb55 100644 --- a/projects/projected_snowfall/comparison_with_scm/comparison_historical_visualizer.py +++ b/projects/projected_snowfall/comparison_with_scm/comparison_historical_visualizer.py @@ -30,7 +30,8 @@ class ComparisonHistoricalVisualizer(StudyVisualizer): def get_values(self, study_method, massif_name): """ - Return an array of size (nb_ensembles + 1) x nb_observations + Return an array "values" of size (nb_ensembles + 1) x nb_observations + Return gcm_rcm_couples of size nb_ensembles :param study_method: :param massif_name: :return: @@ -45,18 +46,49 @@ class ComparisonHistoricalVisualizer(StudyVisualizer): pass return np.array(values), gcm_rcm_couples + def compute_bias_list_in_the_mean(self, massif_name, relative, study_method): + values, gcm_rcm_couples = self.get_values(study_method, massif_name) + mean_values = np.mean(values, axis=1) + bias_in_the_mean = (mean_values - mean_values[0])[1:] + if relative: + bias_in_the_mean *= 100 / mean_values[0] + return bias_in_the_mean, gcm_rcm_couples + + # Map massif name to bias list (ordered by the gcm_rcm_couples) + + def massif_name_to_bias_list_in_the_mean(self, plot_maxima=True, relative=False): + study_method = self.get_study_method(plot_maxima) + massif_name_to_bias_list = {} + for massif_name in self.massif_names: + bias_list, gcm_rcm_couples = self.compute_bias_list_in_the_mean(massif_name, relative, study_method) + massif_name_to_bias_list[massif_name] = bias_list + return massif_name_to_bias_list + + @property + def massif_name_to_rank(self): + massif_name_to_rank = {} + for massif_name, bias_list in self.massif_name_to_bias_list_in_the_mean(plot_maxima=True).items(): + # Count the number of bias negative + nb_of_negative = sum([b < 0 for b in bias_list]) + # Rank starts to 1 + massif_name_to_rank[massif_name] = float(1 + nb_of_negative) + return massif_name_to_rank + + # Map gcm_rcm_couple to bias list (ordered by the massif_names) + + @cached_property + def gcm_rcm_couple_to_bias_list_in_the_mean_maxima(self): + return self.gcm_rcm_couple_to_bias_list_for_the_mean(plot_maxima=True, relative=False) + @cached_property - def mean_bias_maxima(self): - return self.get_mean_bias(plot_maxima=True) + def gcm_rcm_couple_to_relative_bias_list_in_the_mean_maxima(self): + return self.gcm_rcm_couple_to_bias_list_for_the_mean(plot_maxima=True, relative=True) - def get_mean_bias(self, plot_maxima=True): - # Return an array: nb_ensemble x nb_massif + def gcm_rcm_couple_to_bias_list_for_the_mean(self, plot_maxima=True, relative=False): study_method = self.get_study_method(plot_maxima) gcm_rcm_couple_to_bias_list = {couple: [] for couple in self.adamont_studies.gcm_rcm_couples} for massif_name in self.massif_names: - values, gcm_rcm_couples = self.get_values(study_method, massif_name) - mean_values = np.mean(values, axis=1) - bias = (mean_values - mean_values[0])[1:] + bias, gcm_rcm_couples = self.compute_bias_list_in_the_mean(massif_name, relative, study_method) for b, couple in zip(bias, gcm_rcm_couples): gcm_rcm_couple_to_bias_list[couple].append(b) return gcm_rcm_couple_to_bias_list @@ -107,8 +139,8 @@ class ComparisonHistoricalVisualizer(StudyVisualizer): plt.close() def shoe_plot_bias_maxima_comparison(self): - couples = list(self.mean_bias_maxima.keys()) - values = list(self.mean_bias_maxima.values()) + couples = list(self.gcm_rcm_couple_to_bias_list_in_the_mean_maxima.keys()) + values = list(self.gcm_rcm_couple_to_bias_list_in_the_mean_maxima.values()) colors = [gcm_rcm_couple_to_color[couple] for couple in couples] labels = [gcm_rcm_couple_to_str(couple) for couple in couples] @@ -128,3 +160,24 @@ class ComparisonHistoricalVisualizer(StudyVisualizer): self.show_or_save_to_file(add_classic_title=False, no_title=True, tight_layout=True) ax.clear() plt.close() + + def plot_map_with_the_rank(self): + massif_name_to_value = self.massif_name_to_rank + max_abs_change = self.adamont_studies.nb_ensemble_members + 1 + ylabel = 'Rank of the mean maxima\n,' \ + 'which is between 1 (lowest) and {} (largest)'.format(max_abs_change) + plot_name = ylabel + self.plot_map(cmap=plt.cm.coolwarm, graduation=1, + label=ylabel, + massif_name_to_value=massif_name_to_value, + plot_name=plot_name, add_x_label=True, + negative_and_positive_values=False, + altitude=self.altitude, + add_colorbar=True, + max_abs_change=max_abs_change, + massif_name_to_text={m: str(v) for m, v in massif_name_to_value.items()}, + # xlabel=self.altitude_group.xlabel, + ) + + def plot_map_with_the_bias_in_the_mean(self, relative=True): + pass diff --git a/projects/projected_snowfall/comparison_with_scm/comparison_plot.py b/projects/projected_snowfall/comparison_with_scm/comparison_plot.py index 28707f1f..5ebafb6b 100644 --- a/projects/projected_snowfall/comparison_with_scm/comparison_plot.py +++ b/projects/projected_snowfall/comparison_with_scm/comparison_plot.py @@ -8,22 +8,27 @@ from projects.projected_snowfall.comparison_with_scm.comparison_historical_visua ComparisonHistoricalVisualizer -def individual_plot(v): - # v.adamont_studies.plot_maxima_time_series(v.massif_names, v.scm_study) - v.shoe_plot_bias_maxima_comparison() +def individual_plot(visualizer: ComparisonHistoricalVisualizer): + # visualizer.adamont_studies.plot_maxima_time_series_adamont(visualizer.massif_names, visualizer.scm_study) + # visualizer.shoe_plot_bias_maxima_comparison() + # for relative in [True, False]: + # visualizer.plot_map_with_the_bias_in_the_mean(relative) + + visualizer.plot_map_with_the_rank() # for plot_maxima in [True, False][:1]: - # v.plot_comparison(plot_maxima) + # visualizer.plot_comparison(plot_maxima) def collective_plot(altitude_to_visualizer): - bias_of_the_mean_with_the_altitude(altitude_to_visualizer) + # bias_of_the_mean_with_the_altitude(altitude_to_visualizer) + pass def bias_of_the_mean_with_the_altitude(altitude_to_visualizer: Dict[int, ComparisonHistoricalVisualizer]): visualizer = list(altitude_to_visualizer.values())[0] altitudes = list(altitude_to_visualizer.keys()) for couple in visualizer.adamont_studies.gcm_rcm_couples: - values = [v.mean_bias_maxima[couple] for v in altitude_to_visualizer.values()] + values = [v.gcm_rcm_couple_to_bias_list_in_the_mean_maxima[couple] for v in altitude_to_visualizer.values()] ax = plt.gca() width = 10 diff --git a/projects/projected_snowfall/comparison_with_scm/main_comparison_on_beginning_projection.py b/projects/projected_snowfall/comparison_with_scm/main_comparison_on_beginning_projection.py deleted file mode 100644 index c5d682f3..00000000 --- a/projects/projected_snowfall/comparison_with_scm/main_comparison_on_beginning_projection.py +++ /dev/null @@ -1,43 +0,0 @@ - - -import matplotlib as mpl - -from projects.projected_snowfall.comparison_with_scm.main_comparison_on_historical_period import load_vizualizer_histo - -mpl.use('Agg') -mpl.rcParams['text.usetex'] = True -mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}'] -from collections import OrderedDict - -from extreme_data.meteo_france_data.adamont_data.adamont.adamont_snowfall import AdamontSnowfall -from extreme_data.meteo_france_data.adamont_data.adamont_scenario import gcm_rcm_couple_to_full_name -from extreme_data.meteo_france_data.adamont_data.adamont_studies import AdamontStudies -from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day -from extreme_data.meteo_france_data.scm_models_data.utils import Season -from projects.projected_snowfall.comparison_with_scm.comparison_historical_visualizer import \ - ComparisonHistoricalVisualizer - - -def main(): - fast = True - # Set the year_min and year_max for the comparison - if fast: - year_min = [2006][0] - year_max = [2019][0] - massif_names = ['Vanoise'] - altitudes = [1800] - else: - year_min = [2006][0] - year_max = [2019][0] - massif_names = None - altitudes = [900, 1800, 2700] - - # Load visualizers - for altitude in altitudes: - study_class_couple = [(SafranSnowfall1Day, AdamontSnowfall)][0] - v = load_vizualizer_histo(altitude, massif_names, study_class_couple, year_max, year_min) - v.adamont_studies.plot_maxima_time_series(v.massif_names, v.scm_study) - - -if __name__ == '__main__': - main() diff --git a/projects/projected_snowfall/comparison_with_scm/main_comparison_on_historical_period.py b/projects/projected_snowfall/comparison_with_scm/main_comparison_on_rcp_85_extended.py similarity index 72% rename from projects/projected_snowfall/comparison_with_scm/main_comparison_on_historical_period.py rename to projects/projected_snowfall/comparison_with_scm/main_comparison_on_rcp_85_extended.py index f6ef4e35..894613de 100644 --- a/projects/projected_snowfall/comparison_with_scm/main_comparison_on_historical_period.py +++ b/projects/projected_snowfall/comparison_with_scm/main_comparison_on_rcp_85_extended.py @@ -1,6 +1,7 @@ from collections import OrderedDict import matplotlib as mpl + mpl.use('Agg') mpl.rcParams['text.usetex'] = True mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}'] @@ -18,14 +19,14 @@ from projects.projected_snowfall.comparison_with_scm.comparison_historical_visua def main(): - fast = False + fast = True # Set the year_min and year_max for the comparison if fast is None: year_min = [1982, 1950][1] - massif_names = ['Vanoise', 'Vercors', 'Mont-Blanc'] - altitudes = [900, 1800] + massif_names = ['Vanoise'] + altitudes = [1200, 1500, 1800, 2100, 2400] elif fast: - year_min = [1982, 1950][1] + year_min = [1982, 1950][0] massif_names = ['Vanoise'] altitudes = [1800] else: @@ -34,34 +35,31 @@ def main(): altitudes = [600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600][:] # Load visualizers - - # Individual plot - # for altitude in altitudes: - # v = load_visualizer_histo(altitude, massif_names, year_min) - # individual_plot(v) - - # # Collective plot altitude_to_visualizer = OrderedDict() for altitude in altitudes: - altitude_to_visualizer[altitude] = load_visualizer_histo(altitude, massif_names, year_min) + visualizer = load_visualizer(altitude, massif_names, year_min) + altitude_to_visualizer[altitude] = visualizer + # Individual plot + individual_plot(visualizer) + # Collective plot collective_plot(altitude_to_visualizer) -def load_visualizer_histo(altitude, massif_names, year_min): +def load_visualizer(altitude, massif_names, year_min) -> ComparisonHistoricalVisualizer: year_min = max(1959, year_min) - year_max = 2005 + year_max = 2019 study_class_couple = [(SafranSnowfall1Day, AdamontSnowfall)][0] - return load_vizualizer_histo(altitude, massif_names, study_class_couple, year_max, year_min) - - -def load_vizualizer_histo(altitude, massif_names, study_class_couple, year_max, year_min): scm_study_class, adamont_study_class = study_class_couple season = Season.annual + adamont_scenario = AdamontScenario.rcp85_extended + + # Loading part scm_study = scm_study_class(altitude=altitude, year_min=year_min, year_max=year_max, season=season) - gcm_rcm_couples = load_gcm_rcm_couples_for_year_min_and_year_max(year_min, year_max, adamont_scenario=AdamontScenario.histo) + gcm_rcm_couples = load_gcm_rcm_couples_for_year_min_and_year_max(year_min, year_max, + adamont_scenario=adamont_scenario) adamont_studies = AdamontStudies(adamont_study_class, gcm_rcm_couples, altitude=altitude, year_min=year_min, year_max=year_max, season=season, - scenario=AdamontScenario.histo) + scenario=adamont_scenario) visualizer = ComparisonHistoricalVisualizer(scm_study, adamont_studies, massif_names=massif_names) return visualizer -- GitLab