[projections] add test_massif_names_adamont_v2 to check that all massifs are...

[projections] add test_massif_names_adamont_v2 to check that all massifs are the same than for SAFRAN 2020.
add assertion to check to the year_min and year_max for each scenario.
compute bias in the mean maxima at the masssif level

extreme_data/meteo_france_data/adamont_data/abstract_adamont_study.py

@@ -13,7 +13,7 @@ from extreme_data.meteo_france_data.adamont_data.adamont.adamont_variables impor
 from extreme_data.meteo_france_data.adamont_data.adamont_gcm_rcm_couples import get_gcm_rcm_couple_adamont_to_full_name
 from extreme_data.meteo_france_data.adamont_data.adamont_scenario import scenario_to_str, AdamontScenario, \
     get_year_min_and_year_max_from_scenario, get_suffix_for_the_nc_file, \
-    scenario_to_real_scenarios, get_year_max
+    scenario_to_real_scenarios, get_year_max, adamont_scenarios_real
 from extreme_data.meteo_france_data.adamont_data.utils.utils import massif_number_to_massif_name
 
 from extreme_data.utils import DATA_PATH
@@ -89,6 +89,12 @@ class AbstractAdamontStudy(AbstractStudy):
             year_to_annual_maxima = OrderedDict()
             year_min, year_max = get_year_min_and_year_max_from_scenario(self.scenario, self.gcm_rcm_couple)
             years = list(range(year_min, year_max + 1))
+            if self.scenario in adamont_scenarios_real:
+                time = dataset.variables['time']
+                msg = 'len_years={} while len_time={},' \
+                      'check year_min and year_max, ' \
+                      'check in debug mode the time field of the daatset to see the starting date'.format(years, time)
+                assert len(years) == len(time), msg
             for year, maxima in zip(years, annual_maxima):
                 if self.year_min <= year <= self.year_max:
                     year_to_annual_maxima[year] = maxima
@@ -110,11 +116,10 @@ class AbstractAdamontStudy(AbstractStudy):
         scenario_name = self._scenario_to_str_adamont_v2(scenario)
         directory = self.gcm_rcm_full_name + '_' + scenario_name
         filename = self.nc_filename_adamont_v2(scenario)
-        print(directory)
-        print(filename)
         full_path = op.join(ADAMONT_v2_WEBPATH, directory, filename)
         # Download file
         request = 'wget {} -P {}'.format(full_path, path_folder)
+        print(request)
         subprocess.run(request, shell=True)
 
     def nc_filename_adamont_v2(self, scenario):

extreme_data/meteo_france_data/adamont_data/adamont_scenario.py

@@ -43,7 +43,8 @@ def get_year_min(adamont_scenario, gcm_rcm_couple):
             year_min = 1982
         elif rcm == 'RCA4':
             year_min = 1971
-        elif gcm_rcm_couple in [('NorESM1-M', 'HIRHAM5'), ('IPSL-CM5A-MR', 'WRF331F')]:
+        elif gcm_rcm_couple in [('NorESM1-M', 'HIRHAM5'), ('IPSL-CM5A-MR', 'WRF331F'), ('CNRM-CM5', 'ALADIN63'),
+                                ('IPSL-CM5A-MR', 'WRF381P')]:
             year_min = 1952
         else:
             year_min = 1951

projects/projected_snowfall/comparison_with_scm/main_comparison_on_quantile_period.py

@@ -20,29 +20,44 @@ from projects.altitude_spatial_model.altitudes_fit.altitudes_studies import Alti
 def compute_bias_and_display_it(ax,
                                 altitude_studies_reanalysis: AltitudesStudies,
                                 adamont_altitude_studies: AltitudesStudies,
-                                gcm_rcm_couple
+                                gcm_rcm_couple,
+                                massif_names=None,
                                 ):
     bias_in_the_mean_maxima = []
     altitudes = []
     for altitude, study_reanalysis in altitude_studies_reanalysis.altitude_to_study.items():
-        altitudes.append(altitude)
         adamont_study = adamont_altitude_studies.altitude_to_study[altitude]
-        mean_maxima_adamont = adamont_study.mean_annual_maxima
-        mean_maxima_reanalysis = study_reanalysis.mean_annual_maxima
-        bias = mean_maxima_adamont - mean_maxima_reanalysis
-        bias_in_the_mean_maxima.append(bias)
+        can_compute_biais = (massif_names is None) or any([m in adamont_study.study_massif_names for m in massif_names])
+        if can_compute_biais:
+
+            altitudes.append(altitude)
+
+            if massif_names is None:
+                mean_maxima_adamont = adamont_study.mean_annual_maxima
+                mean_maxima_reanalysis = study_reanalysis.mean_annual_maxima
+            else:
+                mean_maxima_reanalysis = np.mean(np.concatenate([study_reanalysis.massif_name_to_annual_maxima[m]
+                                                                 for m in massif_names
+                                                                 if m in study_reanalysis.massif_name_to_annual_maxima]))
+                mean_maxima_adamont = np.mean(np.concatenate([adamont_study.massif_name_to_annual_maxima[m]
+                                                              for m in massif_names
+                                                              if m in adamont_study.massif_name_to_annual_maxima]))
+
+            bias = mean_maxima_adamont - mean_maxima_reanalysis
+            bias_in_the_mean_maxima.append(bias)
 
     color = gcm_rcm_couple_to_color[gcm_rcm_couple]
     label = gcm_rcm_couple_to_str(gcm_rcm_couple)
     ax.plot(bias_in_the_mean_maxima, altitudes, label=label, color=color)
 
-    return np.array(bias_in_the_mean_maxima)
+    return np.array(bias_in_the_mean_maxima), altitudes
 
 
 def main_comparaison_plot():
     altitudes = [600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600][:]
-    for adamont_version in [1, 2]:
-        ax = plt.gca()
+    for adamont_version in [1, 2][1:]:
+        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
@@ -57,47 +72,67 @@ def main_comparaison_plot():
                                                             altitudes=altitudes,
                                                             year_min=1988,
                                                             year_max=2011)
-        bias_in_the_mean = []
-        for gcm_rcm_couple in gcm_rcm_couples:
-            print(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_in_the_mean.append(compute_bias_and_display_it(ax, reanalysis_altitude_studies,
-                                                                adamont_altitude_studies, gcm_rcm_couple))
-
-        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, 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)]
-        plot_name = 'Bias for annual maxima of {}'.format(study_str)
-        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)
-        ax.set_xlabel('Bias in the mean annual maxima of {} for ADAMONT v{} members\n'
-                      ' against {} on the quantile mapping period ({})'.format(study_str, adamont_version,
-                                                                               comparaison_study_class,
-                                                                               study.variable_unit), fontsize=10)
-        reanalysis_altitude_studies.show_or_save_to_file(plot_name=plot_name, no_title=True)
-        plt.close()
+
+        if adamont_version == 1:
+            list_of_massis_names = [None]
+        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()
 
 
 if __name__ == '__main__':

test/test_extreme_data/test_meteo_france_data/test_adamont_study.py

@@ -3,6 +3,7 @@ import unittest
 from extreme_data.meteo_france_data.adamont_data.adamont_gcm_rcm_couples import get_gcm_rcm_couple_adamont_to_full_name
 from extreme_data.meteo_france_data.adamont_data.adamont_scenario import AdamontScenario
 from extreme_data.meteo_france_data.adamont_data.adamont.adamont_snowfall import AdamontSnowfall
+from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day
 
 
 class TestAdamontStudy(unittest.TestCase):
@@ -27,6 +28,16 @@ class TestAdamontStudy(unittest.TestCase):
                 # print(len(study.massif_name_to_annual_maxima['Vanoise']))
             self.assertTrue(True)
 
+    def test_massifs_names_adamont_v2(self):
+        year_min = 2004
+        adamont_version = 2 # this test will not pass with adamont version 1
+        for altitude in [600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600]:
+            reanalysis_study = SafranSnowfall1Day(altitude=altitude, year_min=year_min)
+            for gcm_rcm_couple in get_gcm_rcm_couple_adamont_to_full_name(adamont_version).keys():
+                adamont_study = AdamontSnowfall(altitude=altitude, adamont_version=adamont_version,
+                                                year_min=year_min, gcm_rcm_couple=gcm_rcm_couple)
+                assert set(adamont_study.study_massif_names) == set(reanalysis_study.study_massif_names)
+
 
 if __name__ == '__main__':
     unittest.main()