From dd68ea52c7c96e82b642eac964720bd555231ff8 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Wed, 10 Mar 2021 16:51:11 +0100
Subject: [PATCH] [projection snowfall] add preliminary analysis for the work
on model with effect related to the GCM or RCM type.
---
.../adamont_data/cmip5/temperature_to_year.py | 4 +-
.../visualizer_for_projection_ensemble.py | 42 ++++++++++
extreme_trend/one_fold_fit/one_fold_fit.py | 7 +-
.../main_sensitivity.py | 38 ++++-----
.../preliminary_works/__init__.py | 0
.../main_preliminary_works.py | 80 +++++++++++++++++++
6 files changed, 148 insertions(+), 23 deletions(-)
create mode 100644 projects/projected_extreme_snowfall/preliminary_works/__init__.py
create mode 100644 projects/projected_extreme_snowfall/preliminary_works/main_preliminary_works.py
diff --git a/extreme_data/meteo_france_data/adamont_data/cmip5/temperature_to_year.py b/extreme_data/meteo_france_data/adamont_data/cmip5/temperature_to_year.py
index 564877ee..6453d9f0 100644
--- a/extreme_data/meteo_france_data/adamont_data/cmip5/temperature_to_year.py
+++ b/extreme_data/meteo_france_data/adamont_data/cmip5/temperature_to_year.py
@@ -69,7 +69,7 @@ def plot_nb_data(is_temperature_interval, is_shift_interval):
ax = plt.gca()
for gcm in get_gcm_list(adamont_version=2)[:]:
- for i, scenario in enumerate(rcp_scenarios[:2]):
+ for i, scenario in enumerate(rcp_scenarios[2:]):
plot_nb_data_one_line(ax, gcm, scenario, left_limit, right_limit,
i == 0, is_temperature_interval)
@@ -118,6 +118,6 @@ def get_ticks_labels_for_interval(is_temperature_interval, is_shift_interval):
if __name__ == '__main__':
for shift_interval in [False, True]:
- for temp_interval in [False, True]:
+ for temp_interval in [False, True][1:]:
print("shift = {}, temp_inteval = {}".format(shift_interval, temp_interval))
plot_nb_data(is_temperature_interval=temp_interval, is_shift_interval=shift_interval)
diff --git a/extreme_trend/ensemble_fit/visualizer_for_projection_ensemble.py b/extreme_trend/ensemble_fit/visualizer_for_projection_ensemble.py
index 42b38f05..bf53041d 100644
--- a/extreme_trend/ensemble_fit/visualizer_for_projection_ensemble.py
+++ b/extreme_trend/ensemble_fit/visualizer_for_projection_ensemble.py
@@ -1,6 +1,11 @@
+import matplotlib.pyplot as plt
from collections import OrderedDict
from typing import List
+from extreme_data.meteo_france_data.adamont_data.adamont_gcm_rcm_couples import gcm_rcm_couple_to_color
+from extreme_data.meteo_france_data.adamont_data.adamont_scenario import gcm_rcm_couple_to_str
+from extreme_data.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
+from extreme_fit.distribution.gev.gev_params import GevParams
from extreme_fit.model.margin_model.polynomial_margin_model.spatio_temporal_polynomial_model import \
AbstractSpatioTemporalPolynomialModel
from extreme_fit.model.margin_model.utils import MarginFitMethod
@@ -129,3 +134,40 @@ class VisualizerForProjectionEnsemble(object):
return [ensemble_class_to_ensemble_fit[ensemble_class]
for ensemble_class_to_ensemble_fit
in self.altitude_class_to_ensemble_class_to_ensemble_fit.values()]
+
+ def plot_preliminary_first_part(self):
+ if self.massif_names is None:
+ massif_names = AbstractStudy.all_massif_names()
+ else:
+ massif_names = self.massif_names
+ assert isinstance(massif_names, list)
+ # Plot for all parameters
+ for param_name in GevParams.PARAM_NAMES:
+ for degree in [0, 1]:
+ for massif_name in massif_names:
+ self.plot_preliminary_first_part_for_one_massif(massif_name, param_name, degree)
+
+ def plot_preliminary_first_part_for_one_massif(self, massif_name, param_name, degree):
+ # Retrieve the data
+ ensemble_fit: IndependentEnsembleFit
+ gcm_rcm_couple_to_data = {
+ c: [] for c in self.gcm_rcm_couples
+ }
+ for ensemble_fit in self.ensemble_fits(IndependentEnsembleFit):
+ for gcm_rcm_couple in self.gcm_rcm_couples:
+ visualizer = ensemble_fit.gcm_rcm_couple_to_visualizer[gcm_rcm_couple]
+ one_fold_fit = visualizer.massif_name_to_one_fold_fit[massif_name]
+ coef = one_fold_fit.best_coef(param_name, 0, degree)
+ altitude = visualizer.altitude_group.reference_altitude
+ gcm_rcm_couple_to_data[gcm_rcm_couple].append((altitude, coef))
+ # Plot
+ ax = plt.gca()
+ for gcm_rcm_couple, data in gcm_rcm_couple_to_data.items():
+ altitudes, coefs = list(zip(*data))
+ color = gcm_rcm_couple_to_color[gcm_rcm_couple]
+ label = gcm_rcm_couple_to_str(gcm_rcm_couple)
+ ax.plot(coefs, altitudes, color=color, label=label, marker='o')
+ ax.legend()
+ visualizer.plot_name = '{}/{}_{}'.format(param_name, degree, massif_name)
+ visualizer.show_or_save_to_file(no_title=True)
+ plt.close()
\ No newline at end of file
diff --git a/extreme_trend/one_fold_fit/one_fold_fit.py b/extreme_trend/one_fold_fit/one_fold_fit.py
index 82a3a39e..b62a83d3 100644
--- a/extreme_trend/one_fold_fit/one_fold_fit.py
+++ b/extreme_trend/one_fold_fit/one_fold_fit.py
@@ -246,8 +246,11 @@ class OneFoldFit(object):
def best_coef(self, param_name, dim, degree):
try:
coef = self.best_function_from_fit.param_name_to_coef[param_name] # type: PolynomialAllCoef
- coef = coef.dim_to_polynomial_coef[dim] # type: PolynomialCoef
- coef = coef.idx_to_coef[degree]
+ if coef.dim_to_polynomial_coef is None:
+ return coef.intercept
+ else:
+ coef = coef.dim_to_polynomial_coef[dim] # type: PolynomialCoef
+ coef = coef.idx_to_coef[degree]
return coef
except (TypeError, KeyError):
return None
diff --git a/projects/projected_extreme_snowfall/main_sensitivity.py b/projects/projected_extreme_snowfall/main_sensitivity.py
index 39b2d243..a2e5e60c 100644
--- a/projects/projected_extreme_snowfall/main_sensitivity.py
+++ b/projects/projected_extreme_snowfall/main_sensitivity.py
@@ -8,17 +8,16 @@ from extreme_trend.ensemble_fit.together_ensemble_fit.together_ensemble_fit impo
matplotlib.use('Agg')
import matplotlib as mpl
+
mpl.rcParams['text.usetex'] = True
mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
-
from extreme_trend.ensemble_fit.independent_ensemble_fit.independent_ensemble_fit import IndependentEnsembleFit
from extreme_trend.ensemble_fit.visualizer_for_projection_ensemble import VisualizerForProjectionEnsemble
from extreme_trend.ensemble_fit.visualizer_for_sensitivity import VisualizerForSensivity
from spatio_temporal_dataset.coordinates.temporal_coordinates.temperature_covariate import \
AnomalyTemperatureWithSplineTemporalCovariate
-
from extreme_fit.model.margin_model.polynomial_margin_model.utils import \
ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS
@@ -27,8 +26,8 @@ from extreme_fit.model.utils import set_seed_for_test
from extreme_data.meteo_france_data.adamont_data.adamont.adamont_safran import AdamontSnowfall
from extreme_data.meteo_france_data.adamont_data.adamont_scenario import AdamontScenario, get_gcm_rcm_couples, \
rcp_scenarios
-from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_covariate_for_fit import TimeTemporalCovariate
-
+from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_covariate_for_fit import \
+ TimeTemporalCovariate
from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_extract_eurocode_return_level import \
AbstractExtractEurocodeReturnLevel
@@ -42,8 +41,6 @@ def main():
start = time.time()
study_class = AdamontSnowfall
ensemble_fit_classes = [IndependentEnsembleFit, TogetherEnsembleFit][1:]
- temporal_covariate_for_fit = [TimeTemporalCovariate,
- AnomalyTemperatureWithSplineTemporalCovariate][0]
set_seed_for_test()
AbstractExtractEurocodeReturnLevel.ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY = 0.2
@@ -71,29 +68,32 @@ def main():
assert isinstance(altitudes_list, List)
assert isinstance(altitudes_list[0], List)
print('Scenario is', scenario)
- print('Covariate is {}'.format(temporal_covariate_for_fit))
model_classes = ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS
assert scenario in rcp_scenarios
remove_physically_implausible_models = True
+ temp_cov = False
+ temporal_covariate_for_fit = AnomalyTemperatureWithSplineTemporalCovariate if temp_cov else TimeTemporalCovariate
+ print('Covariate is {}'.format(temporal_covariate_for_fit))
- visualizer = VisualizerForSensivity(
- altitudes_list, gcm_rcm_couples, study_class, Season.annual, scenario,
- model_classes=model_classes,
- ensemble_fit_classes=ensemble_fit_classes,
- massif_names=massif_names,
- temporal_covariate_for_fit=temporal_covariate_for_fit,
- remove_physically_implausible_models=remove_physically_implausible_models,
- is_temperature_interval=False,
- is_shift_interval=False,
- )
- visualizer.plot()
+ for is_temperature_interval in [True, False][1:]:
+ for is_shift_interval in [True, False][1:]:
+ visualizer = VisualizerForSensivity(
+ altitudes_list, gcm_rcm_couples, study_class, Season.annual, scenario,
+ model_classes=model_classes,
+ ensemble_fit_classes=ensemble_fit_classes,
+ massif_names=massif_names,
+ temporal_covariate_for_fit=temporal_covariate_for_fit,
+ remove_physically_implausible_models=remove_physically_implausible_models,
+ is_temperature_interval=is_temperature_interval,
+ is_shift_interval=is_shift_interval,
+ )
+ visualizer.plot()
end = time.time()
duration = str(datetime.timedelta(seconds=end - start))
print('Total duration', duration)
-
if __name__ == '__main__':
main()
diff --git a/projects/projected_extreme_snowfall/preliminary_works/__init__.py b/projects/projected_extreme_snowfall/preliminary_works/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/projects/projected_extreme_snowfall/preliminary_works/main_preliminary_works.py b/projects/projected_extreme_snowfall/preliminary_works/main_preliminary_works.py
new file mode 100644
index 00000000..4d9d61a2
--- /dev/null
+++ b/projects/projected_extreme_snowfall/preliminary_works/main_preliminary_works.py
@@ -0,0 +1,80 @@
+import datetime
+import time
+from typing import List
+import matplotlib
+
+from extreme_fit.distribution.gev.gev_params import GevParams
+from extreme_fit.model.margin_model.polynomial_margin_model.gev_altitudinal_models import StationaryAltitudinal
+from extreme_trend.ensemble_fit.together_ensemble_fit.together_ensemble_fit import TogetherEnsembleFit
+
+matplotlib.use('Agg')
+import matplotlib as mpl
+
+mpl.rcParams['text.usetex'] = True
+mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
+
+from extreme_trend.ensemble_fit.independent_ensemble_fit.independent_ensemble_fit import IndependentEnsembleFit
+from extreme_trend.ensemble_fit.visualizer_for_projection_ensemble import VisualizerForProjectionEnsemble
+from extreme_trend.ensemble_fit.visualizer_for_sensitivity import VisualizerForSensivity
+from spatio_temporal_dataset.coordinates.temporal_coordinates.temperature_covariate import \
+ AnomalyTemperatureWithSplineTemporalCovariate
+
+from extreme_fit.model.margin_model.polynomial_margin_model.utils import \
+ ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS
+
+from extreme_fit.model.utils import set_seed_for_test
+
+from extreme_data.meteo_france_data.adamont_data.adamont.adamont_safran import AdamontSnowfall
+from extreme_data.meteo_france_data.adamont_data.adamont_scenario import AdamontScenario, get_gcm_rcm_couples, \
+ rcp_scenarios
+from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_covariate_for_fit import \
+ TimeTemporalCovariate
+
+from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_extract_eurocode_return_level import \
+ AbstractExtractEurocodeReturnLevel
+
+from extreme_trend.one_fold_fit.altitude_group import altitudes_for_groups
+
+from extreme_data.meteo_france_data.scm_models_data.utils import Season
+
+
+def main(scenario=AdamontScenario.rcp85):
+ start = time.time()
+ study_class = AdamontSnowfall
+ set_seed_for_test()
+ gcm_rcm_couples = get_gcm_rcm_couples(scenario)
+ model_classes = [StationaryAltitudinal]
+
+ fast = False
+
+ if fast is None:
+ massif_names = ['Vanoise']
+ altitudes_list = altitudes_for_groups[:]
+ elif fast:
+ massif_names = ['Vanoise']
+ gcm_rcm_couples = gcm_rcm_couples[:1]
+ altitudes_list = altitudes_for_groups[2:]
+
+ else:
+ massif_names = None
+ altitudes_list = altitudes_for_groups[:]
+
+ visualizer = VisualizerForProjectionEnsemble(
+ altitudes_list, gcm_rcm_couples, study_class, Season.annual, scenario,
+ model_classes=model_classes,
+ ensemble_fit_classes=[IndependentEnsembleFit],
+ massif_names=massif_names,
+ temporal_covariate_for_fit=None,
+ remove_physically_implausible_models=True,
+ gcm_to_year_min_and_year_max=None,
+ )
+ visualizer.plot_preliminary_first_part()
+
+
+ end = time.time()
+ duration = str(datetime.timedelta(seconds=end - start))
+ print('Total duration', duration)
+
+
+if __name__ == '__main__':
+ main()
--
GitLab