[projections] add relative bias, and some improvement on the name

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',

projects/projected_snowfall/comparison_with_scm/main_comparison_on_quantile_period.py

 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__':