diff --git a/extreme_data/meteo_france_data/adamont_data/adamont_gcm_rcm_couples.py b/extreme_data/meteo_france_data/adamont_data/adamont_gcm_rcm_couples.py index f1a787b9a08c8e85ef47a33b72c0cd4b900ef55a..ce707f95209cdf54b373d7584b6c9c791eaeecc0 100644 --- a/extreme_data/meteo_france_data/adamont_data/adamont_gcm_rcm_couples.py +++ b/extreme_data/meteo_france_data/adamont_data/adamont_gcm_rcm_couples.py @@ -1,4 +1,12 @@ def get_year_min_and_year_max_used_to_compute_quantile(gcm): + """pour les périodes de calcul de quantiles/quantiles, voici le recap : + 1980-2011 pour la réanalyse + 1974-2005 pour le modèle + + sauf GCM=MOHC-HadGEM2-ES + 1987-2011 pour la réanalyse + 1981-2005 pour le modèle + """ if gcm == 'HadGEM2-ES': reanalysis_years = (1988, 2011) model_year = (1982, 2005) diff --git a/extreme_data/meteo_france_data/adamont_data/cmip5/temperature_to_year.py b/extreme_data/meteo_france_data/adamont_data/cmip5/temperature_to_year.py index 61872585cc2a5189b190128c0b9fa5a825d6ec26..f8018362a66561139e5a7f50bfd185a90ae70fb1 100644 --- a/extreme_data/meteo_france_data/adamont_data/cmip5/temperature_to_year.py +++ b/extreme_data/meteo_france_data/adamont_data/cmip5/temperature_to_year.py @@ -96,9 +96,9 @@ def get_interval_limits(is_temperature_interval, is_shift_interval): left_limit, right_limit = temp_min, temp_max else: shift = 25 - nb = 3 - year_min = [1959 + shift * i for i in range(nb)] - year_max = [2020 + shift * i for i in range(nb)] + nb = 1 + year_min = [1950 + shift * i for i in range(nb)] + year_max = [2100 + shift * i for i in range(nb)] left_limit, right_limit = year_min, year_max if not is_shift_interval: min_interval_left = min(left_limit) @@ -122,8 +122,52 @@ def get_ticks_labels_for_interval(is_temperature_interval, is_shift_interval): return ticks_labels -if __name__ == '__main__': +def plot_nb_of_data(): for shift_interval in [False, True][:1]: for temp_interval in [False, True][1:]: print("shift = {}, temp_inteval = {}".format(shift_interval, temp_interval)) plot_nb_data(is_temperature_interval=temp_interval, is_shift_interval=shift_interval) + + +def plot_year_for_temperature_interval(): + left_limit, right_limit = get_interval_limits(True, False) + + ax = plt.gca() + for gcm in get_gcm_list(adamont_version=2)[:]: + for i, scenario in enumerate(rcp_scenarios[2:]): + first_scenario = i == 0 + plot_year_for_temperature_interval_one_line(ax, gcm, scenario, left_limit, right_limit, first_scenario) + + ax.legend(loc='upper left') + ticks_labels = get_ticks_labels_for_interval(True, False) + ax.set_yticks(right_limit) + ax.set_yticklabels(ticks_labels) + ax.set_xlabel('Last year considered') + ax2 = ax.twinx() + legend_elements = [ + Line2D([0], [0], color='k', lw=1, label=scenario_to_str(s), + linestyle=get_linestyle_from_scenario(s)) for s in adamont_scenarios_real + ] + ax2.legend(handles=legend_elements, loc='lower right') + ax2.set_yticks([]) + plt.show() + + +def plot_year_for_temperature_interval_one_line(ax, gcm, scenario, left_limits, right_limits, first_scenario): + year_max_list = [get_year_min_and_year_max(gcm, scenario, left_limit, right_limit, + True)[1] for left_limit, right_limit in zip(left_limits, right_limits)] + + color = gcm_to_color[gcm] + linestyle = get_linestyle_from_scenario(scenario) + + # For the legend + if (len(year_max_list) > 0) and first_scenario: + ax.plot(year_max_list[0], right_limits[0], linestyle='solid', color=color, label=gcm) + + ax.plot(year_max_list, right_limits, linestyle=linestyle, color=color, marker='o') + + +if __name__ == '__main__': + # plot_nb_of_data() + plot_year_for_temperature_interval() + diff --git a/extreme_trend/ensemble_fit/visualizer_for_sensitivity.py b/extreme_trend/ensemble_fit/visualizer_for_sensitivity.py index 0796858790f638e7d904dda8cd77dfa61ffe5c21..37150b16ed99b34e9ad2a74b9358110c904cc673 100644 --- a/extreme_trend/ensemble_fit/visualizer_for_sensitivity.py +++ b/extreme_trend/ensemble_fit/visualizer_for_sensitivity.py @@ -48,7 +48,7 @@ class VisualizerForSensivity(object): self.left_limit_to_right_limit = OrderedDict(zip(self.left_limits, self.right_limits)) self.right_limit_to_visualizer = {} # type: Dict[float, VisualizerForProjectionEnsemble] - + print(self.right_limits) for left_limit, right_limit in zip(self.left_limits, self.right_limits): interval_str_prefix = "{}-{}".format(left_limit, right_limit) interval_str_prefix = interval_str_prefix.replace('.', ',') @@ -88,18 +88,24 @@ class VisualizerForSensivity(object): merge_visualizer_str_list.append(AbstractEnsembleFit.Together_merge) for merge_visualizer_str in merge_visualizer_str_list: self.sensitivity_plot_percentages(merge_visualizer_str) - self.sensitivity_plot_return_levels(merge_visualizer_str) + for only_mean in [True, False]: + self.sensitivity_plot_return_levels(merge_visualizer_str, only_mean) for relative in [True, False]: self.sensitivity_plot_changes(merge_visualizer_str, relative) - def sensitivity_plot_return_levels(self, merge_visualizer_str): + def sensitivity_plot_return_levels(self, merge_visualizer_str, only_mean): ax = plt.gca() for altitudes in self.altitudes_list: altitude_group = get_altitude_class_from_altitudes(altitudes) - self.interval_plot_return_levels(ax, altitude_group, merge_visualizer_str) + self.interval_plot_return_levels(ax, altitude_group, merge_visualizer_str, only_mean) + # todo: dans la fonction return level de sensitivity visualizer faire en sorte de ne pas tracer que le retour de niveau 100, + # mais en tracer d autres pour voir comment ils évoluent. ticks_labels = get_ticks_labels_for_interval(self.is_temperature_interval, self.is_shift_interval) - name = 'Return levels at the end of the interval' + if only_mean: + name = 'Mean return levels at the end of the interval' + else: + name = 'Return levels at the end of the interval' if len(self.altitudes_list) == 1: altitude_group = get_altitude_group_from_altitudes(self.altitudes_list[0]) name += ' at {} m'.format(altitude_group.reference_altitude) @@ -109,7 +115,7 @@ class VisualizerForSensivity(object): ax.set_xticklabels(ticks_labels) self.save_plot(merge_visualizer_str, name) - def interval_plot_return_levels(self, ax, altitude_class, merge_visualizer_str): + def interval_plot_return_levels(self, ax, altitude_class, merge_visualizer_str, only_mean): mean_return_levels = [] massif_names = AbstractStudy.all_massif_names() if self.massif_names is None else self.massif_names @@ -121,18 +127,21 @@ class VisualizerForSensivity(object): massif_name_to_tuple_list[massif_name].append((r, return_level)) mean_return_levels.append(np.mean(list(massif_name_to_return_level.values()))) - for massif_id, massif_name in enumerate(AbstractStudy.all_massif_names()): - if len(massif_name_to_tuple_list[massif_name]) > 0: - x, y = zip(*massif_name_to_tuple_list[massif_name]) - else: - x, y = [], [] - color, linestyle, label = get_color_and_linestyle_from_massif_id(massif_id, massif_name) - ax.plot(x, y, color=color, linestyle=linestyle, label=label) + if not only_mean: + for massif_id, massif_name in enumerate(AbstractStudy.all_massif_names()): + if (massif_name in massif_name_to_tuple_list) and (len(massif_name_to_tuple_list[massif_name]) > 0): + x, y = zip(*massif_name_to_tuple_list[massif_name]) + else: + x, y = [], [] + color, linestyle, label = get_color_and_linestyle_from_massif_id(massif_id, massif_name) + ax.plot(x, y, color=color, linestyle=linestyle, label=label) ax.plot(self.right_limits, mean_return_levels, label="Mean", color='k', linewidth=4) - down, up = ax.get_ylim() - ax.set_ylim((down, up * 1.3)) - ax.legend(prop={'size': 7}, ncol=3, loc="upper center") + + if not only_mean: + down, up = ax.get_ylim() + ax.set_ylim((down, up * 1.3)) + ax.legend(prop={'size': 7}, ncol=3, loc="upper center") def interval_plot_changes(self, ax, altitude_class, merge_visualizer_str, relative): linestyle = get_linestyle_for_altitude_class(altitude_class) diff --git a/extreme_trend/one_fold_fit/one_fold_fit.py b/extreme_trend/one_fold_fit/one_fold_fit.py index 5146cb0495e9d46ec3aae765333c7fae4e91b09b..e8a593839253018e175874095e71cee536f9914a 100644 --- a/extreme_trend/one_fold_fit/one_fold_fit.py +++ b/extreme_trend/one_fold_fit/one_fold_fit.py @@ -388,6 +388,7 @@ class OneFoldFit(object): df_temporal_covariate = self.dataset.coordinates.df_temporal_coordinates_for_fit(temporal_covariate_for_fit=self.temporal_covariate_for_fit, drop_duplicates=False) last_temporal_coordinate = df_temporal_covariate.loc[:, AbstractCoordinates.COORDINATE_T].max() + # todo: améliorer the last temporal coordinate. on recupère la liste des rights_limits, puis on prend la valeur juste au dessus ou égale.""" print('last temporal coordinate', last_temporal_coordinate) altitude = self.altitude_group.reference_altitude coordinate = np.array([altitude, last_temporal_coordinate]) diff --git a/projects/past_extreme_snowfall/section_data_and_results/main_altitudes_studies.py b/projects/past_extreme_snowfall/section_data_and_results/main_altitudes_studies.py index 1801c8d272e79b785f93a9efc6439cee8c73f36f..7c389bd69f648418ced56aa6fc2450f6525d1829 100644 --- a/projects/past_extreme_snowfall/section_data_and_results/main_altitudes_studies.py +++ b/projects/past_extreme_snowfall/section_data_and_results/main_altitudes_studies.py @@ -41,9 +41,10 @@ def main(): study_classes = [SafranSnowfall1Day , SafranSnowfall3Days, SafranSnowfall5Days, SafranSnowfall7Days][:1] - study_classes = [SafranSnowfall2019, SafranSnowfall2020, SafranSnowfallCenterOnDay1day, + study_classes = [SafranSnowfall2020, SafranSnowfall2019, SafranSnowfallCenterOnDay1day, SafranSnowfallNotCenterOnDay1day, - SafranSnowfallCenterOnDay1dayMeanRate, SafranSnowfall1Day][1:2] + SafranSnowfallCenterOnDay1dayMeanRate, SafranSnowfall1Day][1:3] + study_classes = [SafranSnowfallNotCenterOnDay1day, SafranSnowfall2019] seasons = [Season.annual, Season.winter, Season.spring, Season.automn][:1] set_seed_for_test() @@ -88,10 +89,10 @@ def main_loop(altitudes_list, massif_names, seasons, study_classes, model_must_p def plot_visualizers(massif_names, visualizer_list): # plot_histogram_all_models_against_altitudes(massif_names, visualizer_list) - plot_histogram_all_trends_against_altitudes(massif_names, visualizer_list, with_significance=False) + plot_histogram_all_trends_against_altitudes(massif_names, visualizer_list, with_significance=True) # plot_shoe_plot_ratio_interval_size_against_altitude(massif_names, visualizer_list) for relative in [True, False]: - plot_shoe_plot_changes_against_altitude(massif_names, visualizer_list, relative=relative, with_significance=False) + plot_shoe_plot_changes_against_altitude(massif_names, visualizer_list, relative=relative, with_significance=True) # plot_coherence_curves(['Vanoise'], visualizer_list) pass diff --git a/projects/projected_extreme_snowfall/comparison_window_for_the_max/__init__.py b/projects/projected_extreme_snowfall/comparison_window_for_the_max/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/projects/projected_extreme_snowfall/comparison_window_for_the_max/comparison_window_for_the_max.py b/projects/projected_extreme_snowfall/comparison_window_for_the_max/comparison_window_for_the_max.py new file mode 100644 index 0000000000000000000000000000000000000000..a35fbc30b587818bc617d749107b14b8ee2edc03 --- /dev/null +++ b/projects/projected_extreme_snowfall/comparison_window_for_the_max/comparison_window_for_the_max.py @@ -0,0 +1,47 @@ +import matplotlib.pyplot as plt +from extreme_data.meteo_france_data.scm_models_data.altitudes_studies import AltitudesStudies +from extreme_data.meteo_france_data.scm_models_data.safran.gap_between_study import \ + GapBetweenSafranSnowfall2019AndMySafranSnowfall2019 +from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day, SafranSnowfallCenterOnDay, \ + SafranSnowfallCenterOnDay1day +from extreme_data.meteo_france_data.scm_models_data.visualization.main_study_visualizer import \ + STUDY_CLASS_TO_ABBREVIATION + + +def plot_comparison(altitudes, massif_name, stud_class_list): + altitude_studies_list = [AltitudesStudies(study_class, altitudes) for study_class in stud_class_list] + + ax = plt.gca() + + for altitude_studies in altitude_studies_list: + plot_studies(ax, altitude_studies, massif_name) + + labelsize = 10 + + ax.tick_params(axis='both', which='major', labelsize=labelsize) + handles, labels = ax.get_legend_handles_labels() + ax.legend(handles[::-1], labels[::-1], prop={'size': labelsize}) + + plot_name = 'comparaison for {} at {}'.format(massif_name, altitudes[0]) + altitude_studies_list[0].show_or_save_to_file(plot_name=plot_name, show=False, no_title=True, tight_layout=True) + ax.clear() + plt.close() + +def plot_studies(ax, altitude_studies, massif_name): + for altitude, study in list(altitude_studies.altitude_to_study.items()): + x = study.ordered_years + if massif_name in study.massif_name_to_annual_maxima: + y = study.massif_name_to_annual_maxima[massif_name] + label = '{} m'.format(altitude) + ax.plot(x, y, linewidth=2, label=label) + study_class = type(study) + plot_name = 'Annual maxima of {} in {}'.format(STUDY_CLASS_TO_ABBREVIATION[study_class], + massif_name.replace('_', ' ')) + fontsize = 10 + ax.set_ylabel('{} ({})'.format(plot_name, study.variable_unit), fontsize=fontsize) + + +if __name__ == '__main__': + for study_class in [GapBetweenSafranSnowfall2019AndMySafranSnowfall2019, SafranSnowfall1Day, SafranSnowfallCenterOnDay1day][:1]: + for altitudes, massif_name in [([2100, 2400, 2700], 'Chablais'), ([1200, 1500, 1800], 'Vercors')][:]: + plot_comparison(altitudes, massif_name, [study_class]) \ No newline at end of file diff --git a/projects/projected_extreme_snowfall/discussion/main_sensitivity.py b/projects/projected_extreme_snowfall/discussion/main_sensitivity.py index 6ed567fffbb88571767cab01d1eaf7539e36defd..0b466dd125faa5a7b244195e9aa0c6cb252f234e 100644 --- a/projects/projected_extreme_snowfall/discussion/main_sensitivity.py +++ b/projects/projected_extreme_snowfall/discussion/main_sensitivity.py @@ -45,7 +45,7 @@ def main(): set_seed_for_test() AbstractExtractEurocodeReturnLevel.ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY = 0.2 - fast = None + fast = False scenarios = [AdamontScenario.rcp85] scenarios = rcp_scenarios[1:] scenarios = rcm_scenarios_extended[2:][::-1] @@ -59,7 +59,7 @@ def main(): massif_names = None gcm_rcm_couples = gcm_rcm_couples[4:6] AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP = 10 - altitudes_list = altitudes_for_groups[1:3] + altitudes_list = altitudes_for_groups[2:3] elif fast: massif_names = ['Vanoise', 'Haute-Maurienne'] gcm_rcm_couples = gcm_rcm_couples[4:6] @@ -78,11 +78,11 @@ def main(): model_classes = ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS assert scenario is not AdamontScenario.histo remove_physically_implausible_models = True - temp_cov = True + temp_cov = False temporal_covariate_for_fit = AnomalyTemperatureWithSplineTemporalCovariate if temp_cov else TimeTemporalCovariate print('Covariate is {}'.format(temporal_covariate_for_fit)) - for is_temperature_interval in [True, False][:1]: + for is_temperature_interval in [True, False][1:]: for is_shift_interval in [True, False][1:]: visualizer = VisualizerForSensivity( altitudes_list, gcm_rcm_couples, study_class, Season.annual, scenario, diff --git a/projects/projected_extreme_snowfall/discussion/main_sensitivity_one_altitude_range.py b/projects/projected_extreme_snowfall/discussion/main_sensitivity_one_altitude_range.py index d421a5c744ebfb3922476d61650d28ad807ffe30..dce827a629991941ad58c1650722b08bcd75260e 100644 --- a/projects/projected_extreme_snowfall/discussion/main_sensitivity_one_altitude_range.py +++ b/projects/projected_extreme_snowfall/discussion/main_sensitivity_one_altitude_range.py @@ -48,7 +48,7 @@ def main(): fast = False scenarios = [AdamontScenario.rcp85] scenarios = rcp_scenarios[1:] - scenarios = rcm_scenarios_extended[:][::-1] + scenarios = rcm_scenarios_extended[2:][::-1] if fast in [None, True]: scenarios = scenarios[:1] @@ -64,10 +64,10 @@ def main(): massif_names = ['Vanoise', 'Haute-Maurienne', 'Vercors', 'Ubaye'] gcm_rcm_couples = gcm_rcm_couples[4:6] AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP = 10 - altitudes_list = altitudes_for_groups[1:3] + altitudes_list = altitudes_for_groups[1:2] else: massif_names = None - altitudes_list = altitudes_for_groups[:] + altitudes_list = altitudes_for_groups[2:3] assert isinstance(gcm_rcm_couples, list) diff --git a/projects/projected_extreme_snowfall/results/main_projections_ensemble.py b/projects/projected_extreme_snowfall/results/main_projections_ensemble.py index edc30d24a592284d255ed34ebea453fea9511ac7..c44d3cbf17ef9914f913d54e9d64564184c687c0 100644 --- a/projects/projected_extreme_snowfall/results/main_projections_ensemble.py +++ b/projects/projected_extreme_snowfall/results/main_projections_ensemble.py @@ -24,7 +24,7 @@ from extreme_fit.model.utils import set_seed_for_test from extreme_data.meteo_france_data.adamont_data.adamont.adamont_safran import AdamontSnowfall from extreme_data.meteo_france_data.adamont_data.adamont_scenario import AdamontScenario, get_gcm_rcm_couples, \ - rcp_scenarios + rcp_scenarios, rcm_scenarios_extended from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_covariate_for_fit import \ TimeTemporalCovariate @@ -47,6 +47,7 @@ def main(): fast = True scenarios = rcp_scenarios[::-1] if fast is False else [AdamontScenario.rcp85] + scenarios = rcm_scenarios_extended[::-1] for scenario in scenarios: gcm_rcm_couples = get_gcm_rcm_couples(scenario) @@ -72,7 +73,7 @@ def main(): print('Covariate is {}'.format(temporal_covariate_for_fit)) model_classes = ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS - assert scenario in rcp_scenarios + assert scenario is not AdamontScenario.histo visualizer = VisualizerForProjectionEnsemble( altitudes_list, gcm_rcm_couples, study_class, Season.annual, scenario,