diff --git a/extreme_data/meteo_france_data/scm_models_data/visualization/main_study_visualizer.py b/extreme_data/meteo_france_data/scm_models_data/visualization/main_study_visualizer.py index 9e346c0c97042fa70054c1a72607257a68adb5c2..9130c65cc940f8c8f8dc0586288e38cc8f2123d6 100644 --- a/extreme_data/meteo_france_data/scm_models_data/visualization/main_study_visualizer.py +++ b/extreme_data/meteo_france_data/scm_models_data/visualization/main_study_visualizer.py @@ -18,9 +18,11 @@ from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranS SafranRainfall, \ SafranTemperature, SafranPrecipitation, SafranSnowfall1Day, SafranSnowfall3Days, SafranSnowfall5Days, \ SafranSnowfall7Days, SafranPrecipitation1Day, SafranPrecipitation3Days, SafranPrecipitation5Days, \ - SafranPrecipitation7Days, SafranDateFirstSnowfall + SafranPrecipitation7Days, SafranDateFirstSnowfall, SafranSnowfallCenterOnDay1dayMeanRate, \ + SafranSnowfallCenterOnDay1day from extreme_data.meteo_france_data.scm_models_data.safran.safran_max_snowf import SafranSnowfall2020, \ SafranSnowfall2019 +from extreme_data.meteo_france_data.scm_models_data.safran.safran_variable import SafranSnowfallVariableCenterOnDay from extreme_data.meteo_france_data.scm_models_data.visualization.study_visualizer import \ StudyVisualizer from projects.exceeding_snow_loads.section_discussion.crocus_study_comparison_with_eurocode import \ @@ -39,13 +41,14 @@ SCM_STUDIES = [SafranSnowfall, CrocusSweTotal, CrocusDepth, CrocusSwe3Days] SCM_STUDIES_NAMES = [get_display_name_from_object_type(k) for k in SCM_STUDIES] SCM_STUDY_NAME_TO_SCM_STUDY = dict(zip(SCM_STUDIES_NAMES, SCM_STUDIES)) - # I keep the scm study separated from the adamont study (for the tests) SCM_STUDY_CLASS_TO_ABBREVIATION = { SafranSnowfall: 'SF3', SafranSnowfall1Day: 'daily snowfall', SafranSnowfall2020: 'daily snowfall', SafranSnowfall2019: 'daily snowfall', + SafranSnowfallCenterOnDay1dayMeanRate: 'daily snowfall', + SafranSnowfallCenterOnDay1day: 'daily snowfall', GapBetweenSafranSnowfall2019And2020: 'daily snowfall\n bias = SAFRAN 2020 minus SAFRAN 2019', GapBetweenSafranSnowfall2019AndMySafranSnowfall2019Recentered: 'daily snowfall\n my SAFRAN 2019 recentered minus SAFRAN 2019', GapBetweenSafranSnowfall2019AndMySafranSnowfall2019NotRecentered: 'daily snowfall\n my SAFRAN 2019 notrecentered minus SAFRAN 2019', diff --git a/projects/projected_snowfall/comparison_with_scm/main_comparison_on_quantile_period.py b/projects/projected_snowfall/comparison_with_scm/main_comparison_on_quantile_period.py index a4289f85b2d2adc4ecea3954fe8f8fd18eb62979..8317e6e9f54391a8b8db6c5db309f8dca3f6076d 100644 --- a/projects/projected_snowfall/comparison_with_scm/main_comparison_on_quantile_period.py +++ b/projects/projected_snowfall/comparison_with_scm/main_comparison_on_quantile_period.py @@ -1,4 +1,5 @@ import numpy as np +import os.path as op import matplotlib from extreme_data.meteo_france_data.adamont_data.adamont.adamont_snowfall import AdamontSnowfall @@ -6,8 +7,11 @@ from extreme_data.meteo_france_data.adamont_data.adamont_gcm_rcm_couples import get_year_min_and_year_max_used_to_compute_quantile from extreme_data.meteo_france_data.adamont_data.adamont_scenario import AdamontScenario, gcm_rcm_couple_to_color, \ gcm_rcm_couple_to_str, load_gcm_rcm_couples -from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day -from extreme_data.meteo_france_data.scm_models_data.safran.safran_max_snowf import SafranSnowfall2020 +from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day, \ + SafranSnowfallCenterOnDay1dayMeanRate, SafranSnowfallCenterOnDay1day +from extreme_data.meteo_france_data.scm_models_data.safran.safran_max_snowf import SafranSnowfall2020, \ + SafranSnowfall2019 +from extreme_data.meteo_france_data.scm_models_data.safran.safran_variable import SafranSnowfallVariableCenterOnDay matplotlib.use('Agg') @@ -22,6 +26,7 @@ def compute_bias_and_display_it(ax, adamont_altitude_studies: AltitudesStudies, gcm_rcm_couple, massif_names=None, + relative_bias=False, ): bias_in_the_mean_maxima = [] altitudes = [] @@ -44,6 +49,8 @@ def compute_bias_and_display_it(ax, if m in adamont_study.massif_name_to_annual_maxima])) bias = mean_maxima_adamont - mean_maxima_reanalysis + if relative_bias: + bias *= 100 / mean_maxima_reanalysis bias_in_the_mean_maxima.append(bias) color = gcm_rcm_couple_to_color[gcm_rcm_couple] @@ -55,12 +62,14 @@ def compute_bias_and_display_it(ax, def main_comparaison_plot(): altitudes = [600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600][:] - for adamont_version in [1, 2][1:]: + + for adamont_version in [1, 2][:]: print('version:', adamont_version) gcm_rcm_couples = load_gcm_rcm_couples(adamont_scenario=AdamontScenario.histo, adamont_version=adamont_version) - study_class = SafranSnowfall2020 if adamont_version == 2 else SafranSnowfall1Day + study_class_for_adamont_v1 = SafranSnowfall1Day + study_class = SafranSnowfall2020 if adamont_version == 2 else study_class_for_adamont_v1 comparaison_study_class = 'SAFRAN 2020' if adamont_version == 2 else 'SAFRAN 2019' # Faster to load once the two cases @@ -78,61 +87,67 @@ def main_comparaison_plot(): else: list_of_massis_names = [None] + [[m] for m in reanalysis_altitude_studies_1981.study.all_massif_names()] - for massif_names in list_of_massis_names[:]: - ax = plt.gca() - bias_in_the_mean = [] - list_altitudes_for_bias = [] - for gcm_rcm_couple in gcm_rcm_couples: - print(massif_names, gcm_rcm_couple) - gcm, rcm = gcm_rcm_couple - years_reanalysis, years_model = get_year_min_and_year_max_used_to_compute_quantile(gcm) - assert years_reanalysis[0] in [1981, 1988] - if years_reanalysis[0] == 1981: - reanalysis_altitude_studies = reanalysis_altitude_studies_1981 - else: - reanalysis_altitude_studies = reanalysis_altitude_studies_1988 - adamont_altitude_studies = AltitudesStudies(study_class=AdamontSnowfall, - altitudes=altitudes, - year_min=years_model[0], - year_max=years_model[1], - scenario=AdamontScenario.histo, - gcm_rcm_couple=gcm_rcm_couple, - adamont_version=adamont_version) - bias, altitudes_for_bias = compute_bias_and_display_it(ax, reanalysis_altitude_studies, - adamont_altitude_studies, gcm_rcm_couple, massif_names) - bias_in_the_mean.append(bias) - list_altitudes_for_bias.append(altitudes_for_bias) - - # Assert the all the bias have been computed for the same altitudes - altitudes_for_bias = list_altitudes_for_bias[0] - for alti in list_altitudes_for_bias: - assert alti == altitudes_for_bias - - bias_in_the_mean = np.array(bias_in_the_mean) - min_bias, median_bias, max_bias = [f(bias_in_the_mean, axis=0) for f in [np.min, np.median, np.max]] - - # Plot the range for the bias, and the median - ax.yaxis.set_ticks(altitudes) - color = 'k' - ax.plot(median_bias, altitudes_for_bias, label='Median bias', color=color, linewidth=4) - # ax.fill_betweenx(altitudes, min_bias, max_bias, label='Range for the bias', alpha=0.2, color='whitesmoke') - ax.vlines(0, ymin=altitudes[0], ymax=altitudes[-1], color='k', linestyles='dashed') - study_str = STUDY_CLASS_TO_ABBREVIATION[type(reanalysis_altitude_studies.study)] - ax.set_ylim(top=altitudes[-1] + 1300) - study = adamont_altitude_studies.study - ax.legend(ncol=3, prop={'size': 7}) - ax.set_ylabel('Altitude (m)', fontsize=10) - massif_str = 'all massifs' if massif_names is None else 'the {} massif'.format(massif_names[0]) - title = 'Bias in the mean annual maxima of {} of {}\n' \ - 'for ADAMONT v{}' \ - ' against {} on the quantile mapping period ({})'.format(study_str, massif_str, - adamont_version, - comparaison_study_class, - study.variable_unit) - plot_name = title - ax.set_xlabel(title, fontsize=10) - reanalysis_altitude_studies.show_or_save_to_file(plot_name=plot_name, no_title=True) - plt.close() + for relative_bias in [True, False][:1]: + for massif_names in list_of_massis_names[:]: + ax = plt.gca() + bias_in_the_mean = [] + list_altitudes_for_bias = [] + for gcm_rcm_couple in gcm_rcm_couples: + print(massif_names, gcm_rcm_couple) + gcm, rcm = gcm_rcm_couple + years_reanalysis, years_model = get_year_min_and_year_max_used_to_compute_quantile(gcm) + assert years_reanalysis[0] in [1981, 1988] + if years_reanalysis[0] == 1981: + reanalysis_altitude_studies = reanalysis_altitude_studies_1981 + else: + reanalysis_altitude_studies = reanalysis_altitude_studies_1988 + adamont_altitude_studies = AltitudesStudies(study_class=AdamontSnowfall, + altitudes=altitudes, + year_min=years_model[0], + year_max=years_model[1], + scenario=AdamontScenario.histo, + gcm_rcm_couple=gcm_rcm_couple, + adamont_version=adamont_version) + bias, altitudes_for_bias = compute_bias_and_display_it(ax, reanalysis_altitude_studies, + adamont_altitude_studies, gcm_rcm_couple, massif_names, + relative_bias) + bias_in_the_mean.append(bias) + list_altitudes_for_bias.append(altitudes_for_bias) + + # Assert the all the bias have been computed for the same altitudes + altitudes_for_bias = list_altitudes_for_bias[0] + for alti in list_altitudes_for_bias: + assert alti == altitudes_for_bias + + bias_in_the_mean = np.array(bias_in_the_mean) + min_bias, median_bias, max_bias = [f(bias_in_the_mean, axis=0) for f in [np.min, np.median, np.max]] + + # Plot the range for the bias, and the median + ax.yaxis.set_ticks(altitudes) + color = 'k' + ax.plot(median_bias, altitudes_for_bias, label='Median bias', color=color, linewidth=4) + # ax.fill_betweenx(altitudes, min_bias, max_bias, label='Range for the bias', alpha=0.2, color='whitesmoke') + ax.vlines(0, ymin=altitudes[0], ymax=altitudes[-1], color='k', linestyles='dashed') + study_str = STUDY_CLASS_TO_ABBREVIATION[type(reanalysis_altitude_studies.study)] + ax.set_ylim(top=altitudes[-1] + 1300) + study = adamont_altitude_studies.study + ax.legend(ncol=3, prop={'size': 7}) + ax.set_ylabel('Altitude (m)', fontsize=10) + massif_str = 'all massifs' if massif_names is None else 'the {} massif'.format(massif_names[0]) + unit = '%' if relative_bias else study.variable_unit + bias_name = 'Relative bias' if relative_bias else 'Bias' + title = '{} in the mean annual maxima of {} of {}\n' \ + 'for ADAMONT v{}' \ + ' against {} on the quantile mapping period ({})'.format(bias_name, + study_str, massif_str, + adamont_version, + comparaison_study_class, + unit) + folder = 'relative bias' if relative_bias else 'absolute bias' + plot_name = op.join(folder, title) + ax.set_xlabel(title, fontsize=10) + reanalysis_altitude_studies.show_or_save_to_file(plot_name=plot_name, no_title=True) + plt.close() if __name__ == '__main__':