diff --git a/extreme_fit/model/margin_model/polynomial_margin_model/utils.py b/extreme_fit/model/margin_model/polynomial_margin_model/utils.py
index 7a61879eac82d1c93f243d4265164d24faca06f8..75671432d34f946e60cd169e96d73d1e117f8e2b 100644
--- a/extreme_fit/model/margin_model/polynomial_margin_model/utils.py
+++ b/extreme_fit/model/margin_model/polynomial_margin_model/utils.py
@@ -81,13 +81,13 @@ ALTITUDINAL_GEV_MODELS_LOCATION_QUADRATIC_MINIMUM = [
# NonStationaryAltitudinalLocationLinear,
NonStationaryAltitudinalLocationQuadratic,
NonStationaryAltitudinalLocationCubic,
- NonStationaryAltitudinalLocationOrder4,
+ # NonStationaryAltitudinalLocationOrder4,
# # First order cross term
# NonStationaryCrossTermForLocation,
# NonStationaryAltitudinalLocationLinearCrossTermForLocation,
NonStationaryAltitudinalLocationQuadraticCrossTermForLocation,
NonStationaryAltitudinalLocationCubicCrossTermForLocation,
- NonStationaryAltitudinalLocationOrder4CrossTermForLocation,
+ # NonStationaryAltitudinalLocationOrder4CrossTermForLocation,
]
diff --git a/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py b/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
index f34480115283fb87afca6cf768cd46b4d3cfe65c..f22c1703dd40299744944672a84e58e7957f9305 100644
--- a/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
+++ b/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
@@ -37,8 +37,8 @@ def plot_moments(studies, massif_names=None):
def plot_altitudinal_fit(studies, massif_names=None):
# model_classes = ALTITUDINAL_GEV_MODELS_LOCATION_ONLY_SCALE_ALTITUDES
- # model_classes = ALTITUDINAL_GEV_MODELS_LOCATION_QUADRATIC_MINIMUM
- model_classes = ALTITUDINAL_GEV_MODELS_LOCATION
+ model_classes = ALTITUDINAL_GEV_MODELS_LOCATION_QUADRATIC_MINIMUM
+ # model_classes = ALTITUDINAL_GEV_MODELS_LOCATION
# model_classes = ALTITUDINAL_GEV_MODELS_LOCATION_CUBIC_MINIMUM
# model_classes = ALTITUDINAL_GEV_MODELS_QUADRATIC
visualizer = AltitudesStudiesVisualizerForNonStationaryModels(studies=studies,
diff --git a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
index 38afd43d980db6ed995b384676cac7962146e3dc..1bce4578aae56629bb8f06edba7d433603b81b78 100644
--- a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
+++ b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
@@ -2,6 +2,7 @@ from typing import List, Dict
import matplotlib.pyplot as plt
import numpy as np
+from cached_property import cached_property
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 \
@@ -9,6 +10,7 @@ from extreme_data.meteo_france_data.scm_models_data.visualization.main_study_vis
from extreme_data.meteo_france_data.scm_models_data.visualization.plot_utils import plot_against_altitude
from extreme_data.meteo_france_data.scm_models_data.visualization.study_visualizer import StudyVisualizer
from extreme_fit.distribution.gev.gev_params import GevParams
+from extreme_fit.function.margin_function.abstract_margin_function import AbstractMarginFunction
from extreme_fit.function.param_function.linear_coef import LinearCoef
from extreme_fit.model.margin_model.polynomial_margin_model.spatio_temporal_polynomial_model import \
AbstractSpatioTemporalPolynomialModel
@@ -40,7 +42,8 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
self._massif_name_to_one_fold_fit = {}
for massif_name in self.massif_names:
dataset = studies.spatio_temporal_dataset(massif_name=massif_name)
- old_fold_fit = OneFoldFit(massif_name, dataset, model_classes, self.fit_method, self.temporal_covariate_for_fit)
+ old_fold_fit = OneFoldFit(massif_name, dataset, model_classes, self.fit_method,
+ self.temporal_covariate_for_fit)
self._massif_name_to_one_fold_fit[massif_name] = old_fold_fit
@property
@@ -144,9 +147,10 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
def plot_altitude_for_the_peak(self):
pass
- def plot_year_for_the_peak(self):
+ def plot_year_for_the_peak(self, plot_mean=True):
t_list = 1959 + np.arange(141)
# t_list = 1800 + np.arange(400)
+ return_period = 50
for massif_name, one_fold_fit in self.massif_name_to_one_fold_fit.items():
ax = plt.gca()
# One plot for each altitude
@@ -158,21 +162,101 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
nearest_altitude = self.studies.altitudes[i]
nearest_study = self.studies.altitude_to_study[nearest_altitude]
if massif_name in nearest_study.study_massif_names:
- mean_list = []
+ y_list = []
for t in t_list:
coordinate = np.array([altitude, t])
gev_params = one_fold_fit.best_function_from_fit.get_params(coordinate, is_transformed=False)
- mean_list.append(gev_params.mean)
+ if plot_mean:
+ y = gev_params.mean
+ else:
+ y = gev_params.return_level(return_period=return_period)
+ y_list.append(y)
label = '{} m'.format(altitude)
- ax.plot(t_list, mean_list, label=label)
+ ax.plot(t_list, y_list, label=label)
ax.legend()
# Modify the limits of the y axis
lim_down, lim_up = ax.get_ylim()
ax_lim = (0, lim_up)
ax.set_ylim(ax_lim)
ax.set_xlabel('Year')
- ax.set_ylabel('Mean annual maxima of {}'.format(SCM_STUDY_CLASS_TO_ABBREVIATION[type(self.study)]))
- plot_name = '{}/{}/Peak year for {}'.format(OneFoldFit.folder_for_plots, 'Peak year', massif_name.replace('_', ''))
+ if plot_mean:
+ ylabel = 'Mean annual maxima'
+ else:
+ ylabel = '{}-year return level'.format(return_period)
+ ax.set_ylabel('{} of {}'.format(ylabel, SCM_STUDY_CLASS_TO_ABBREVIATION[type(self.study)]))
+ peak_year_folder = 'Peak year ' + ylabel
+ plot_name = '{}/{}/Peak year for {}'.format(OneFoldFit.folder_for_plots, peak_year_folder,
+ massif_name.replace('_', ''))
self.studies.show_or_save_to_file(plot_name=plot_name, show=self.show)
plt.close()
+ # Plots "altitude switch" and "peak year"
+
+ @property
+ def massif_name_to_is_decreasing_parabol(self):
+ # For the test we only activate the Mont-Blanc massif
+ d = {massif_name: False for massif_name in self.massif_name_to_one_fold_fit.keys()}
+ for massif_name in ['Mont-Blanc', 'Vanoise', 'Aravis', 'Beaufortain', 'Chablais']:
+ d[massif_name] = True
+ return d
+
+ @property
+ def massif_name_to_altitudes_switch_and_peak_years(self):
+ return {massif_name: self.compute_couple_peak_year_and_altitude_switch(massif_name)
+ for massif_name, is_decreasing_parabol in self.massif_name_to_is_decreasing_parabol.items()
+ if is_decreasing_parabol}
+
+ def compute_couple_peak_year_and_altitude_switch(self, massif_name):
+ # Get the altitude limits
+ altitudes = self.study.massif_name_to_altitudes[massif_name]
+ # use a step of 100m for instance
+ step = 10
+ altitudes = list(np.arange(min(altitudes), max(altitudes) + step, step))
+ # Get all the correspond peak years
+ margin_function = self.massif_name_to_one_fold_fit[massif_name].best_function_from_fit
+ peak_years = []
+ year_left = 1900
+ switch_altitudes = []
+ for altitude in altitudes:
+ year_left = self.compute_peak_year(margin_function, altitude, year_left)
+ if year_left > 2020:
+ break
+ peak_years.append(year_left)
+ switch_altitudes.append(altitude)
+ print(switch_altitudes)
+ print(peak_years)
+ return switch_altitudes, peak_years
+
+ def compute_peak_year(self, margin_function: AbstractMarginFunction, altitude, year_left):
+ year_right = year_left + 0.1
+ mean_left = margin_function.get_params(np.array([altitude, year_left])).mean
+ mean_right = margin_function.get_params(np.array([altitude, year_right])).mean
+ print(year_left, year_right, mean_left, mean_right)
+ if mean_right < mean_left:
+ return year_left
+ else:
+ return self.compute_peak_year(margin_function, altitude, year_right)
+
+ def plot_peak_year_against_altitude(self):
+ ax = plt.gca()
+ for massif_name, (altitudes, peak_years) in self.massif_name_to_altitudes_switch_and_peak_years.items():
+ ax.plot(altitudes, peak_years, label=massif_name)
+ ax.legend()
+ ax.set_xlabel('Altitude')
+ ax.set_ylabel('Peak years')
+ plot_name = 'Peak Years'
+ self.studies.show_or_save_to_file(plot_name=plot_name, show=self.show)
+ plt.close()
+
+ def plot_altitude_switch_against_peak_year(self):
+ ax = plt.gca()
+ for massif_name, (altitudes, peak_years) in self.massif_name_to_altitudes_switch_and_peak_years.items():
+ ax.plot(peak_years, altitudes, label=massif_name)
+ ax.legend()
+ ax.set_xlabel('Peak years')
+ ax.set_ylabel('Altitude')
+ plot_name = 'Switch altitude'
+ self.studies.show_or_save_to_file(plot_name=plot_name, show=self.show)
+ plt.close()
+
+
diff --git a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/plot_total_aic.py b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/plot_total_aic.py
index d3f93142d4af9b4eb041bc41b9734c8fa6f8faf4..c097413797a231084e55d660b446af736d17423c 100644
--- a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/plot_total_aic.py
+++ b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/plot_total_aic.py
@@ -12,9 +12,12 @@ from projects.exceeding_snow_loads.utils import dpi_paper1_figure
def plots(visualizer):
visualizer.plot_shape_map()
- visualizer.plot_year_for_the_peak()
+ for plot_mean in [True, False]:
+ visualizer.plot_year_for_the_peak(plot_mean=plot_mean)
visualizer.plot_moments()
visualizer.plot_best_coef_maps()
+ visualizer.plot_peak_year_against_altitude()
+ visualizer.plot_altitude_switch_against_peak_year()
def plot_individual_aic(visualizer):
diff --git a/test/test_extreme_data/test_nasa_data/test_mean_global_temperature.py b/test/test_extreme_data/test_nasa_data/test_mean_global_temperature.py
index e30e4a3d05cb3c8e710f469952926d669bc3113e..c91196b4fdb0c09c1e6dba4335ba8f37db7bed5f 100644
--- a/test/test_extreme_data/test_nasa_data/test_mean_global_temperature.py
+++ b/test/test_extreme_data/test_nasa_data/test_mean_global_temperature.py
@@ -14,5 +14,6 @@ class TestMeanGlobalTemperatures(unittest.TestCase):
value = list(d.values())[0]
self.assertIsInstance(value, float)
+
if __name__ == '__main__':
unittest.main()