From 9f042157f5e66a6a36584d64cd02e9e69ca2340c Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Wed, 11 Nov 2020 18:34:02 +0100
Subject: [PATCH] [snow projection] improve rank & bias plots. remove dataset
attribute for abstract_adamont_study.py.
---
.../adamont_data/abstract_adamont_study.py | 16 +++---
.../adamont_data/adamont_scenario.py | 3 +-
.../comparison_historical_visualizer.py | 55 +++++++++++++++----
.../comparison_with_scm/comparison_plot.py | 35 ++++++++++--
.../main_comparison_on_rcp_85_extended.py | 54 ++++++++++++++----
5 files changed, 127 insertions(+), 36 deletions(-)
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 43eae3c0..2cd805a7 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
@@ -102,23 +102,24 @@ class AbstractAdamontStudy(AbstractStudy):
@cached_property
def flat_mask(self):
- zs_list = [int(e) for e in np.array(self.dataset.variables['ZS'])]
+ zs_list = [int(e) for e in np.array(self.datasets[0].variables['ZS'])]
+ if len(self.datasets) > 1:
+ zs_list_bis = [int(e) for e in np.array(self.datasets[1].variables['ZS'])]
+ assert all([a == b for a, b in zip(zs_list, zs_list_bis)])
return np.array(zs_list) == self.altitude
@cached_property
def study_massif_names(self) -> List[str]:
- massif_ids = np.array(self.dataset.variables['MASSIF_NUMBER'])[self.flat_mask]
+ massif_ids = np.array(self.datasets[0].variables['MASSIF_NUMBER'])[self.flat_mask]
+ if len(self.datasets) > 1:
+ massif_ids_bis = np.array(self.datasets[1].variables['MASSIF_NUMBER'])[self.flat_mask]
+ assert all(massif_ids == massif_ids_bis)
return [massif_number_to_massif_name[massif_id] for massif_id in massif_ids]
@cached_property
def datasets(self):
return [Dataset(file_path) for file_path in self.nc_file_paths]
- @property
- def dataset(self):
- # todo: improve that
- return self.datasets[0]
-
# PATHS
@property
@@ -150,4 +151,5 @@ class AbstractAdamontStudy(AbstractStudy):
scenario_name,
self.region_name, suffix_nc_file)
file_paths.append(op.join(files_path, nc_file))
+ assert len(file_paths) <= 2, 'change my code to handle datasets of length larger than'
return file_paths
diff --git a/extreme_data/meteo_france_data/adamont_data/adamont_scenario.py b/extreme_data/meteo_france_data/adamont_data/adamont_scenario.py
index f95decdb..0791a974 100644
--- a/extreme_data/meteo_france_data/adamont_data/adamont_scenario.py
+++ b/extreme_data/meteo_france_data/adamont_data/adamont_scenario.py
@@ -67,7 +67,8 @@ def get_suffix_for_the_nc_file(adamont_scenario, gcm_rcm_couple):
def scenario_to_str(adamont_scenario):
- return str(adamont_scenario).split('.')[-1].upper()
+ return '+'.join([str(real_adamont_scenario).split('.')[-1].upper()
+ for real_adamont_scenario in scenario_to_real_scenarios(adamont_scenario)])
def scenario_to_real_scenarios(adamont_scenario):
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 72beeb55..652a5581 100644
--- a/projects/projected_snowfall/comparison_with_scm/comparison_historical_visualizer.py
+++ b/projects/projected_snowfall/comparison_with_scm/comparison_historical_visualizer.py
@@ -2,9 +2,11 @@ import matplotlib.pyplot as plt
import numpy as np
from cached_property import cached_property
from matplotlib.lines import Line2D
+import os.path as op
+from extreme_data.meteo_france_data.adamont_data.abstract_adamont_study import AbstractAdamontStudy
from extreme_data.meteo_france_data.adamont_data.adamont_scenario import get_color_from_gcm_rcm_couple, \
- gcm_rcm_couple_to_str, gcm_rcm_couple_to_color
+ gcm_rcm_couple_to_str, gcm_rcm_couple_to_color, scenario_to_str
from extreme_data.meteo_france_data.adamont_data.adamont_studies import AdamontStudies
from extreme_data.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
from extreme_data.meteo_france_data.scm_models_data.visualization.main_study_visualizer import \
@@ -14,6 +16,7 @@ from projects.altitude_spatial_model.altitudes_fit.altitudes_studies import Alti
class ComparisonHistoricalVisualizer(StudyVisualizer):
+ study: AbstractAdamontStudy
def __init__(self, scm_study: AbstractStudy,
adamont_studies: AdamontStudies,
@@ -61,19 +64,24 @@ class ComparisonHistoricalVisualizer(StudyVisualizer):
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
+ if len(bias_list) == len(self.adamont_studies.gcm_rcm_couples):
+ massif_name_to_bias_list[massif_name] = bias_list
return massif_name_to_bias_list
+ def massif_name_to_mean_bias_in_the_mean(self, plot_maxima=True, relative=False):
+ return {m: np.mean(l) for m, l in self.massif_name_to_bias_list_in_the_mean(plot_maxima, relative).items()}
+
@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)
+ # Rank starts to 0
+ massif_name_to_rank[massif_name] = float(nb_of_negative)
return massif_name_to_rank
+
# Map gcm_rcm_couple to bias list (ordered by the massif_names)
@cached_property
@@ -163,10 +171,12 @@ class ComparisonHistoricalVisualizer(StudyVisualizer):
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
+ max_abs_change = self.adamont_studies.nb_ensemble_members
+ ylabel = 'Rank of the mean maxima\n' \
+ 'for the period {}-{} and the scenario {}\n' \
+ 'which is between 0 (lowest) and {} (largest)'.format(self.study.year_min, self.study.year_max,
+ scenario_to_str(self.study.scenario), max_abs_change)
+ plot_name = op.join('rank', ylabel)
self.plot_map(cmap=plt.cm.coolwarm, graduation=1,
label=ylabel,
massif_name_to_value=massif_name_to_value,
@@ -175,9 +185,32 @@ class ComparisonHistoricalVisualizer(StudyVisualizer):
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()},
+ massif_name_to_text={m: str(int(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
+ def plot_map_with_the_mean_bias_in_the_mean(self, relative=True):
+ massif_name_to_value = self.massif_name_to_mean_bias_in_the_mean(plot_maxima=True, relative=relative)
+ name = 'relative difference' if relative else 'difference'
+ abreviation = ADAMONT_STUDY_CLASS_TO_ABBREVIATION[type(self.study)]
+ unit = '\%' if relative else self.study.variable_unit
+ ylabel = 'Mean {} of the ensemble members with the reanalysis\n' \
+ 'for the period {}-{} and the scenario {}\n' \
+ 'in the mean annual maxima of {} ({})'.format(name,
+ self.study.year_min, self.study.year_max,
+ scenario_to_str(self.study.scenario),
+ abreviation, unit)
+ plot_name = op.join(name, ylabel)
+ self.plot_map(cmap=plt.cm.coolwarm, graduation=10,
+ label=ylabel,
+ massif_name_to_value=massif_name_to_value,
+ plot_name=plot_name, add_x_label=True,
+ negative_and_positive_values=True,
+ altitude=self.altitude,
+ add_colorbar=True,
+ massif_name_to_text={m: ('+' if v > 0 else '')
+ + str(int(v))
+ + ('\%' if relative else '')
+ for m, v in massif_name_to_value.items()},
+ # xlabel=self.altitude_group.xlabel,
+ )
diff --git a/projects/projected_snowfall/comparison_with_scm/comparison_plot.py b/projects/projected_snowfall/comparison_with_scm/comparison_plot.py
index 5ebafb6b..7c50dced 100644
--- a/projects/projected_snowfall/comparison_with_scm/comparison_plot.py
+++ b/projects/projected_snowfall/comparison_with_scm/comparison_plot.py
@@ -1,9 +1,11 @@
import matplotlib.pyplot as plt
from typing import Dict
+import numpy as np
from matplotlib.lines import Line2D
-from extreme_data.meteo_france_data.adamont_data.adamont_scenario import gcm_rcm_couple_to_color, gcm_rcm_couple_to_str
+from extreme_data.meteo_france_data.adamont_data.adamont_scenario import gcm_rcm_couple_to_color, gcm_rcm_couple_to_str, \
+ scenario_to_str
from projects.projected_snowfall.comparison_with_scm.comparison_historical_visualizer import \
ComparisonHistoricalVisualizer
@@ -11,17 +13,38 @@ from projects.projected_snowfall.comparison_with_scm.comparison_historical_visua
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)
+ for relative in [True, False]:
+ visualizer.plot_map_with_the_mean_bias_in_the_mean(relative)
visualizer.plot_map_with_the_rank()
+ # for relative
# for plot_maxima in [True, False][:1]:
# visualizer.plot_comparison(plot_maxima)
-def collective_plot(altitude_to_visualizer):
- # bias_of_the_mean_with_the_altitude(altitude_to_visualizer)
- pass
+def collective_plot(altitude_to_visualizer: Dict[int, ComparisonHistoricalVisualizer]):
+ visualizer = list(altitude_to_visualizer.values())[0]
+ count_number_of_total_massifs = 0
+ count_number_of_time_the_reanalysis_is_the_smallest = 0
+ count_number_of_time_the_reanalysis_is_the_biggest = 0
+ altitudes = list(altitude_to_visualizer.keys())
+ all_ranks = []
+ for v in altitude_to_visualizer.values():
+ ranks = np.array(list(v.massif_name_to_rank.values()))
+ count_number_of_total_massifs += len(ranks)
+ count_number_of_time_the_reanalysis_is_the_smallest += sum(ranks == 0.0)
+ all_ranks.extend(ranks)
+ print(scenario_to_str(visualizer.study.scenario), visualizer.study.year_min, visualizer.study.year_max)
+ print('Summary for rank for altitudes:', altitudes)
+ print('Mean ranks:', np.mean(all_ranks))
+ print('percentages of time reanalysis is the biggest:',
+ 100 * count_number_of_time_the_reanalysis_is_the_biggest / count_number_of_total_massifs)
+ print('number of time reanalysis is the biggest:', count_number_of_time_the_reanalysis_is_the_biggest,
+ ' out of ', count_number_of_total_massifs, ' time series')
+ print('percentages of time reanalysis is the smallest:',
+ 100 * count_number_of_time_the_reanalysis_is_the_smallest / count_number_of_total_massifs)
+ print('number of time reanalysis is the smallest:', count_number_of_time_the_reanalysis_is_the_smallest,
+ ' out of ', count_number_of_total_massifs, ' time series')
def bias_of_the_mean_with_the_altitude(altitude_to_visualizer: Dict[int, ComparisonHistoricalVisualizer]):
diff --git a/projects/projected_snowfall/comparison_with_scm/main_comparison_on_rcp_85_extended.py b/projects/projected_snowfall/comparison_with_scm/main_comparison_on_rcp_85_extended.py
index 894613de..cda70659 100644
--- a/projects/projected_snowfall/comparison_with_scm/main_comparison_on_rcp_85_extended.py
+++ b/projects/projected_snowfall/comparison_with_scm/main_comparison_on_rcp_85_extended.py
@@ -19,25 +19,51 @@ from projects.projected_snowfall.comparison_with_scm.comparison_historical_visua
def main():
- fast = True
+ fast = 7
+ year_max = 2019
# Set the year_min and year_max for the comparison
- if fast is None:
- year_min = [1982, 1950][1]
- massif_names = ['Vanoise']
- altitudes = [1200, 1500, 1800, 2100, 2400]
- elif fast:
+ if fast is 1:
year_min = [1982, 1950][0]
massif_names = ['Vanoise']
altitudes = [1800]
+ elif fast is 2:
+ year_min = [1982, 1950][0]
+ massif_names = None
+ altitudes = [1800]
+ elif fast is 3:
+ year_min = [1982, 1950][0]
+ massif_names = ['Vanoise']
+ altitudes = [1200, 1500, 1800, 2100, 2400]
+ elif fast is 4:
+ year_max = 2019
+ massif_names = None
+ year_min = 2006
+ altitudes = [600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600][:]
+ elif fast is 5:
+ year_max = 2019
+ massif_names = None
+ year_min = 2013
+ altitudes = [600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600][:]
+ elif fast is 6:
+ year_max = 2012
+ massif_names = None
+ year_min = 2006
+ altitudes = [600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600][:]
+ elif fast is 7:
+ year_max = 2019
+ massif_names = None
+ year_min = 1982
+ altitudes = [600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600][:]
else:
+ year_max = 2005
massif_names = None
- year_min = [1982, 1950][1]
+ year_min = 1982
altitudes = [600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600][:]
# Load visualizers
altitude_to_visualizer = OrderedDict()
for altitude in altitudes:
- visualizer = load_visualizer(altitude, massif_names, year_min)
+ visualizer = load_visualizer(altitude, massif_names, year_min, year_max)
altitude_to_visualizer[altitude] = visualizer
# Individual plot
individual_plot(visualizer)
@@ -45,13 +71,19 @@ def main():
collective_plot(altitude_to_visualizer)
-def load_visualizer(altitude, massif_names, year_min) -> ComparisonHistoricalVisualizer:
+def load_visualizer(altitude, massif_names, year_min, year_max) -> ComparisonHistoricalVisualizer:
year_min = max(1959, year_min)
- year_max = 2019
+ year_max = min(2019, year_max)
study_class_couple = [(SafranSnowfall1Day, AdamontSnowfall)][0]
scm_study_class, adamont_study_class = study_class_couple
season = Season.annual
- adamont_scenario = AdamontScenario.rcp85_extended
+ if year_min <= 2005:
+ if year_max > 2005:
+ adamont_scenario = AdamontScenario.rcp85_extended
+ else:
+ adamont_scenario = AdamontScenario.histo
+ else:
+ adamont_scenario = AdamontScenario.rcp85
# Loading part
scm_study = scm_study_class(altitude=altitude, year_min=year_min, year_max=year_max, season=season)
--
GitLab