Commit 9f042157 authored by Le Roux Erwan's avatar Le Roux Erwan
Browse files

[snow projection] improve rank & bias plots. remove dataset attribute for...

[snow projection] improve rank & bias plots. remove dataset attribute for abstract_adamont_study.py.
parent 5d45d5c3
No related merge requests found
Showing with 127 additions and 36 deletions
+127 -36
...@@ -102,23 +102,24 @@ class AbstractAdamontStudy(AbstractStudy): ...@@ -102,23 +102,24 @@ class AbstractAdamontStudy(AbstractStudy):
@cached_property @cached_property
def flat_mask(self): 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 return np.array(zs_list) == self.altitude
@cached_property @cached_property
def study_massif_names(self) -> List[str]: 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] return [massif_number_to_massif_name[massif_id] for massif_id in massif_ids]
@cached_property @cached_property
def datasets(self): def datasets(self):
return [Dataset(file_path) for file_path in self.nc_file_paths] return [Dataset(file_path) for file_path in self.nc_file_paths]
@property
def dataset(self):
# todo: improve that
return self.datasets[0]
# PATHS # PATHS
@property @property
...@@ -150,4 +151,5 @@ class AbstractAdamontStudy(AbstractStudy): ...@@ -150,4 +151,5 @@ class AbstractAdamontStudy(AbstractStudy):
scenario_name, scenario_name,
self.region_name, suffix_nc_file) self.region_name, suffix_nc_file)
file_paths.append(op.join(files_path, 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 return file_paths
...@@ -67,7 +67,8 @@ def get_suffix_for_the_nc_file(adamont_scenario, gcm_rcm_couple): ...@@ -67,7 +67,8 @@ def get_suffix_for_the_nc_file(adamont_scenario, gcm_rcm_couple):
def scenario_to_str(adamont_scenario): 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): def scenario_to_real_scenarios(adamont_scenario):
......
...@@ -2,9 +2,11 @@ import matplotlib.pyplot as plt ...@@ -2,9 +2,11 @@ import matplotlib.pyplot as plt
import numpy as np import numpy as np
from cached_property import cached_property from cached_property import cached_property
from matplotlib.lines import Line2D 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, \ 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.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.abstract_study import AbstractStudy
from extreme_data.meteo_france_data.scm_models_data.visualization.main_study_visualizer import \ 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 ...@@ -14,6 +16,7 @@ from projects.altitude_spatial_model.altitudes_fit.altitudes_studies import Alti
class ComparisonHistoricalVisualizer(StudyVisualizer): class ComparisonHistoricalVisualizer(StudyVisualizer):
study: AbstractAdamontStudy
def __init__(self, scm_study: AbstractStudy, def __init__(self, scm_study: AbstractStudy,
adamont_studies: AdamontStudies, adamont_studies: AdamontStudies,
...@@ -61,19 +64,24 @@ class ComparisonHistoricalVisualizer(StudyVisualizer): ...@@ -61,19 +64,24 @@ class ComparisonHistoricalVisualizer(StudyVisualizer):
massif_name_to_bias_list = {} massif_name_to_bias_list = {}
for massif_name in self.massif_names: for massif_name in self.massif_names:
bias_list, gcm_rcm_couples = self.compute_bias_list_in_the_mean(massif_name, relative, study_method) 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 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 @property
def massif_name_to_rank(self): def massif_name_to_rank(self):
massif_name_to_rank = {} massif_name_to_rank = {}
for massif_name, bias_list in self.massif_name_to_bias_list_in_the_mean(plot_maxima=True).items(): 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 # Count the number of bias negative
nb_of_negative = sum([b < 0 for b in bias_list]) nb_of_negative = sum([b < 0 for b in bias_list])
# Rank starts to 1 # Rank starts to 0
massif_name_to_rank[massif_name] = float(1 + nb_of_negative) massif_name_to_rank[massif_name] = float(nb_of_negative)
return massif_name_to_rank return massif_name_to_rank
# Map gcm_rcm_couple to bias list (ordered by the massif_names) # Map gcm_rcm_couple to bias list (ordered by the massif_names)
@cached_property @cached_property
...@@ -163,10 +171,12 @@ class ComparisonHistoricalVisualizer(StudyVisualizer): ...@@ -163,10 +171,12 @@ class ComparisonHistoricalVisualizer(StudyVisualizer):
def plot_map_with_the_rank(self): def plot_map_with_the_rank(self):
massif_name_to_value = self.massif_name_to_rank massif_name_to_value = self.massif_name_to_rank
max_abs_change = self.adamont_studies.nb_ensemble_members + 1 max_abs_change = self.adamont_studies.nb_ensemble_members
ylabel = 'Rank of the mean maxima\n,' \ ylabel = 'Rank of the mean maxima\n' \
'which is between 1 (lowest) and {} (largest)'.format(max_abs_change) 'for the period {}-{} and the scenario {}\n' \
plot_name = ylabel '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, self.plot_map(cmap=plt.cm.coolwarm, graduation=1,
label=ylabel, label=ylabel,
massif_name_to_value=massif_name_to_value, massif_name_to_value=massif_name_to_value,
...@@ -175,9 +185,32 @@ class ComparisonHistoricalVisualizer(StudyVisualizer): ...@@ -175,9 +185,32 @@ class ComparisonHistoricalVisualizer(StudyVisualizer):
altitude=self.altitude, altitude=self.altitude,
add_colorbar=True, add_colorbar=True,
max_abs_change=max_abs_change, 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, # xlabel=self.altitude_group.xlabel,
) )
def plot_map_with_the_bias_in_the_mean(self, relative=True): def plot_map_with_the_mean_bias_in_the_mean(self, relative=True):
pass 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,
)
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from typing import Dict from typing import Dict
import numpy as np
from matplotlib.lines import Line2D 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 \ from projects.projected_snowfall.comparison_with_scm.comparison_historical_visualizer import \
ComparisonHistoricalVisualizer ComparisonHistoricalVisualizer
...@@ -11,17 +13,38 @@ from projects.projected_snowfall.comparison_with_scm.comparison_historical_visua ...@@ -11,17 +13,38 @@ from projects.projected_snowfall.comparison_with_scm.comparison_historical_visua
def individual_plot(visualizer: ComparisonHistoricalVisualizer): def individual_plot(visualizer: ComparisonHistoricalVisualizer):
# visualizer.adamont_studies.plot_maxima_time_series_adamont(visualizer.massif_names, visualizer.scm_study) # visualizer.adamont_studies.plot_maxima_time_series_adamont(visualizer.massif_names, visualizer.scm_study)
# visualizer.shoe_plot_bias_maxima_comparison() # visualizer.shoe_plot_bias_maxima_comparison()
# for relative in [True, False]: for relative in [True, False]:
# visualizer.plot_map_with_the_bias_in_the_mean(relative) visualizer.plot_map_with_the_mean_bias_in_the_mean(relative)
visualizer.plot_map_with_the_rank() visualizer.plot_map_with_the_rank()
# for relative
# for plot_maxima in [True, False][:1]: # for plot_maxima in [True, False][:1]:
# visualizer.plot_comparison(plot_maxima) # visualizer.plot_comparison(plot_maxima)
def collective_plot(altitude_to_visualizer): def collective_plot(altitude_to_visualizer: Dict[int, ComparisonHistoricalVisualizer]):
# bias_of_the_mean_with_the_altitude(altitude_to_visualizer) visualizer = list(altitude_to_visualizer.values())[0]
pass 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]): def bias_of_the_mean_with_the_altitude(altitude_to_visualizer: Dict[int, ComparisonHistoricalVisualizer]):
......
...@@ -19,25 +19,51 @@ from projects.projected_snowfall.comparison_with_scm.comparison_historical_visua ...@@ -19,25 +19,51 @@ from projects.projected_snowfall.comparison_with_scm.comparison_historical_visua
def main(): def main():
fast = True fast = 7
year_max = 2019
# Set the year_min and year_max for the comparison # Set the year_min and year_max for the comparison
if fast is None: if fast is 1:
year_min = [1982, 1950][1]
massif_names = ['Vanoise']
altitudes = [1200, 1500, 1800, 2100, 2400]
elif fast:
year_min = [1982, 1950][0] year_min = [1982, 1950][0]
massif_names = ['Vanoise'] massif_names = ['Vanoise']
altitudes = [1800] 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: else:
year_max = 2005
massif_names = None massif_names = None
year_min = [1982, 1950][1] year_min = 1982
altitudes = [600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600][:] altitudes = [600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600][:]
# Load visualizers # Load visualizers
altitude_to_visualizer = OrderedDict() altitude_to_visualizer = OrderedDict()
for altitude in altitudes: 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 altitude_to_visualizer[altitude] = visualizer
# Individual plot # Individual plot
individual_plot(visualizer) individual_plot(visualizer)
...@@ -45,13 +71,19 @@ def main(): ...@@ -45,13 +71,19 @@ def main():
collective_plot(altitude_to_visualizer) 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_min = max(1959, year_min)
year_max = 2019 year_max = min(2019, year_max)
study_class_couple = [(SafranSnowfall1Day, AdamontSnowfall)][0] study_class_couple = [(SafranSnowfall1Day, AdamontSnowfall)][0]
scm_study_class, adamont_study_class = study_class_couple scm_study_class, adamont_study_class = study_class_couple
season = Season.annual 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 # Loading part
scm_study = scm_study_class(altitude=altitude, year_min=year_min, year_max=year_max, season=season) scm_study = scm_study_class(altitude=altitude, year_min=year_min, year_max=year_max, season=season)
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment