diff --git a/extreme_fit/utils.py b/extreme_fit/utils.py new file mode 100644 index 0000000000000000000000000000000000000000..21d3db969eebc2aa144ce7dd4494760c0c48ed88 --- /dev/null +++ b/extreme_fit/utils.py @@ -0,0 +1,11 @@ +import numpy as np +from sklearn.linear_model import LinearRegression + + +def fit_linear_regression(x, y): + X = np.array(x).reshape(-1, 1) + reg = LinearRegression().fit(X, y) + r2_score = reg.score(X, y) + a = reg.coef_ + b = reg.intercept_ + return a, b, r2_score diff --git a/projects/altitude_spatial_model/altitudes_fit/plots/compute_histogram_change_in_total_snowfall.py b/projects/altitude_spatial_model/altitudes_fit/plots/compute_histogram_change_in_total_snowfall.py index aa9786142618415488c8aaf1c8d742fa1dab5b76..c7d69c79e2ed661b891acb1a74dc5b8a3a8722cd 100644 --- a/projects/altitude_spatial_model/altitudes_fit/plots/compute_histogram_change_in_total_snowfall.py +++ b/projects/altitude_spatial_model/altitudes_fit/plots/compute_histogram_change_in_total_snowfall.py @@ -3,10 +3,9 @@ import matplotlib.pyplot as plt import numpy as np from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall +from extreme_fit.utils import fit_linear_regression from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitudes_studies_visualizer_for_non_stationary_models import \ AltitudesStudiesVisualizerForNonStationaryModels -from projects.contrasting_trends_in_snow_loads.article2_snowfall_versus_time_and_altitude.snowfall_plot import \ - fit_linear_regression def compute_changes_in_total_snowfall(visualizer_list: List[ diff --git a/projects/altitude_spatial_model/delta_graph.py b/projects/altitude_spatial_model/delta_graph.py deleted file mode 100644 index 92f57de2610784ca79fb0a68d355bcb47ceeacd6..0000000000000000000000000000000000000000 --- a/projects/altitude_spatial_model/delta_graph.py +++ /dev/null @@ -1,74 +0,0 @@ -import numpy as np -import matplotlib.pyplot as plt -from extreme_data.meteo_france_data.scm_models_data.abstract_study import AbstractStudy -from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day, SafranTemperature, \ - SafranPrecipitation1Day -from extreme_data.meteo_france_data.scm_models_data.visualization.main_study_visualizer import \ - STUDY_CLASS_TO_ABBREVIATION -from extreme_data.meteo_france_data.scm_models_data.visualization.study_visualizer import StudyVisualizer -from projects.contrasting_trends_in_snow_loads.article2_snowfall_versus_time_and_altitude.snowfall_plot import \ - fit_linear_regression - - -def delta_graph(study_class, altitudes, maxima=True): - ax = plt.gca() - all_delta_point_for_regression = [] - # colors = ['orange', 'red', 'blue', 'green', 'yellow'] - for altitude in altitudes: - delta_points = [] - study = study_class(altitude=altitude) # type : AbstractStudy - study_temperature = SafranTemperature(altitude=altitude) - for massif_name in study.study_massif_names: - if maxima: - values = study.massif_name_to_annual_maxima[massif_name] - else: - values = study.massif_name_to_annual_total[massif_name] - values_temperature = study_temperature.massif_name_to_annual_total[massif_name] - delta_point = (compute_delta(values_temperature, relative=False), compute_delta(values, relative=True)) - delta_points.append(delta_point) - all_delta_point_for_regression.extend(delta_points) - # Plot part - x, y = list(zip(*delta_points)) - ax.scatter(x, y, label='{}m'.format(altitude)) - aggregation_str = 'maxima' if maxima else 'total' - ylabel = 'Delta for {} {}'.format(aggregation_str, STUDY_CLASS_TO_ABBREVIATION[study_class]) - ax.set_ylabel(ylabel) - ax.set_xlabel('Delta mean Temperature') - # Plot linear regression - x_all, y_all = list(zip(*all_delta_point_for_regression)) - a, b, r2_score = fit_linear_regression(x_all, y_all) - a = a[0] - x_plot = np.linspace(start=np.min(x_all), stop=np.max(x_all), num=100) - y_plot = a * x_plot + b - rounded_number = [str(np.round(e, 2)) for e in [a, b, r2_score]] - ax.plot(x_plot, y_plot, label='{} x + {} (R2 = {})'.format(*rounded_number)) - visualizer = StudyVisualizer(study=study, show=False, save_to_file=True) - visualizer.plot_name = ylabel - ax.legend() - - # Show / Save plot - visualizer.show_or_save_to_file(no_title=True) - plt.close() - - -def compute_delta(values, relative=True): - index = 30 - before, after = values[:index], values[index:] - mean_before, mean_after = np.mean(before), np.mean(after) - delta = mean_after - mean_before - if relative: - delta *= 100 / mean_before - return delta - - -if __name__ == '__main__': - fast = False - if fast is None: - altitudes = [900, 1800, 2700] - elif fast: - altitudes = [900] - else: - altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600][:] - for study_class in [SafranSnowfall1Day, SafranPrecipitation1Day]: - for maxima in [True, False]: - delta_graph(study_class, altitudes, maxima=maxima) diff --git a/projects/altitude_spatial_model/preliminary_analysis.py b/projects/altitude_spatial_model/preliminary_analysis.py index dcc3559296f891b259e8e927f1033d16cb17e288..f3642f4b1476ad42c8e3d8d5cf07ff20da627da2 100644 --- a/projects/altitude_spatial_model/preliminary_analysis.py +++ b/projects/altitude_spatial_model/preliminary_analysis.py @@ -3,22 +3,18 @@ import os.path as op from itertools import chain import numpy as np -from cached_property import cached_property from extreme_data.meteo_france_data.adamont_data.abstract_adamont_study import AbstractAdamontStudy from extreme_data.meteo_france_data.adamont_data.adamont.adamont_snowfall import AdamontSnowfall from extreme_data.meteo_france_data.adamont_data.adamont_gcm_rcm_couples import _gcm_rcm_couple_adamont_v2_to_full_name -from extreme_data.meteo_france_data.adamont_data.adamont_scenario import adamont_scenarios_real, \ - get_gcm_rcm_couples, rcp_scenarios +from extreme_data.meteo_france_data.adamont_data.adamont_scenario import get_gcm_rcm_couples, rcp_scenarios from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day 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.utils import fit_linear_regression from projects.altitude_spatial_model.altitudes_fit.altitudes_studies import AltitudesStudies from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitude_group import altitudes_for_groups -from projects.contrasting_trends_in_snow_loads.article2_snowfall_versus_time_and_altitude.snowfall_plot import \ - fit_linear_regression -from projects.exceeding_snow_loads.utils import paper_altitudes class PointwiseGevStudyVisualizer(AltitudesStudies): @@ -193,125 +189,6 @@ class PointwiseGevStudyVisualizer(AltitudesStudies): legend=legend) return fit_linear_regression(altitudes, params) - # plot_against_altitude(altitudes, ax, massif_id, massif_name, confidence_intervals, fill=True) - - # Plot against the time - - # @property - # def year_min_and_max_list(self): - # l = [] - # year_min, year_max = 1959, 1989 - # for shift in range(0, 7): - # l.append((year_min + 5 * shift, year_max + 5 * shift)) - # return l - - # @property - # def min_years_for_plot_x_axis(self): - # return [c[0] for c in self.year_min_and_max_list] - # - # def plot_gev_params_against_time_for_all_altitudes(self): - # for altitude in self.altitudes_for_temporal_hypothesis: - # self._plot_gev_params_against_time_for_one_altitude(altitude) - # - # def _plot_gev_params_against_time_for_one_altitude(self, altitude): - # for param_name in GevParams.PARAM_NAMES[:]: - # ax = plt.gca() - # for massif_name in self.study.all_massif_names()[:]: - # self._plot_gev_params_against_time_for_one_altitude_and_one_massif(ax, massif_name, param_name, - # altitude, - # massif_name_as_labels=True) - # ax.legend(prop={'size': 7}, ncol=3) - # ax.set_xlabel('Year') - # ax.set_ylabel(param_name + ' for altitude={}'.format(altitude)) - # xlabels = ['-'.join([str(e) for e in t]) for t in self.year_min_and_max_list] - # ax.set_xticks(self.min_years_for_plot_x_axis) - # ax.set_xticklabels(xlabels) - # # ax.tick_params(labelsize=5) - # plot_name = '{} change /all with years /for altitude={}'.format(param_name, altitude) - # self.show_or_save_to_file(plot_name, no_title=True, tight_layout=True, show=False) - # ax.clear() - # plt.close() - # - # def _plot_gev_params_against_time_for_one_altitude_and_one_massif(self, ax, massif_name, param_name, altitude, - # massif_name_as_labels): - # study = self.altitude_to_study[altitude] - # if massif_name in study.study_massif_names: - # gev_params_list = study.massif_name_to_gev_param_list(self.year_min_and_max_list)[massif_name] - # params = [gev_params.to_dict()[param_name] for gev_params in gev_params_list] - # # params = np.array(params) - # # param_normalized = params / np.sqrt(np.sum(np.power(params, 2))) - # # confidence_intervals = [gev_params.param_name_to_confidence_interval[param_name] for gev_params in - # # gev_params_list] - # massif_id = self.study.all_massif_names().index(massif_name) - # plot_against_altitude(self.min_years_for_plot_x_axis, ax, massif_id, massif_name, params, - # altitude, False, - # massif_name_as_labels) - # # plot_against_altitude(self.years, ax, massif_id, massif_name, confidence_intervals, True) - # - # # plot for each massif against the time - # - # def plot_gev_params_against_time_for_all_massifs(self): - # for massif_name in self.study.all_massif_names(): - # self._plot_gev_params_against_time_for_one_massif(massif_name) - # - # def _plot_gev_params_against_time_for_one_massif(self, massif_name): - # for param_name in GevParams.PARAM_NAMES[:]: - # ax = plt.gca() - # for altitude in self.altitudes_for_temporal_hypothesis: - # self._plot_gev_params_against_time_for_one_altitude_and_one_massif(ax, massif_name, param_name, - # altitude, - # massif_name_as_labels=False) - # ax.legend() - # ax.set_xlabel('Year') - # ax.set_ylabel(param_name + ' for {}'.format(massif_name)) - # xlabels = ['-'.join([str(e) for e in t]) for t in self.year_min_and_max_list] - # ax.set_xticks(self.min_years_for_plot_x_axis) - # ax.set_xticklabels(xlabels) - # plot_name = '{} change /with years /for {}'.format(param_name, massif_name) - # self.show_or_save_to_file(plot_name, no_title=True, tight_layout=True, show=False) - # ax.clear() - # plt.close() - # - # # PLot for each massif the derivative against the time for each altitude - # - # def plot_time_derivative_against_time(self): - # for param_name in GevParams.PARAM_NAMES[:]: - # ax = plt.gca() - # for massif_name in self.study.all_massif_names()[:]: - # self._plot_gev_params_time_derivative_against_altitude_one_massif(ax, massif_name, param_name) - # ax.legend(prop={'size': 7}, ncol=3) - # ax.set_xlabel('Altitude') - # ax.set_ylabel(param_name) - # plot_name = '{} change /time derivative with altitude'.format(param_name) - # self.show_or_save_to_file(plot_name, no_title=True, tight_layout=True, show=False) - # ax.clear() - # plt.close() - # - # def _plot_gev_params_time_derivative_against_altitude_one_massif(self, ax, massif_name, param_name): - # altitudes = [] - # time_derivatives = [] - # for altitude, study in self.altitude_to_study.items(): - # if (massif_name in study.study_massif_names) and ("Mercan" not in massif_name): - # gev_params_list = study.massif_name_to_gev_param_list(self.year_min_and_max_list)[massif_name] - # params = [gev_params.to_dict()[param_name] for gev_params in gev_params_list] - # x = list(range(len(params))) - # y = params - # a = self.get_robust_slope(x, y) - # time_derivatives.append(a) - # altitudes.append(altitude) - # massif_id = self.study.all_massif_names().index(massif_name) - # plot_against_altitude(altitudes, ax, massif_id, massif_name, time_derivatives, fill=False) - # - # def get_robust_slope(self, x, y): - # a, *_ = fit_linear_regression(x=x, y=y) - # a_list = [a] - # for i in range(len(x)): - # x_copy, y_copy = x[:], y[:] - # x_copy.pop(i) - # y_copy.pop(i) - # a, *_ = fit_linear_regression(x=x_copy, y=y_copy) - # a_list.append(a) - # return np.mean(np.array(a_list)) def main_paper2(): diff --git a/projects/contrasting_trends_in_snow_loads/post thesis committee/__init__.py b/projects/archive/__init__.py similarity index 100% rename from projects/contrasting_trends_in_snow_loads/post thesis committee/__init__.py rename to projects/archive/__init__.py diff --git a/projects/contrasting_trends_in_snow_loads/spatial trends/__init__.py b/projects/archive/check_mcmc_convergence_for_return_levels/__init__.py similarity index 100% rename from projects/contrasting_trends_in_snow_loads/spatial trends/__init__.py rename to projects/archive/check_mcmc_convergence_for_return_levels/__init__.py diff --git a/projects/others/check_mcmc_convergence_for_return_levels/gelman_convergence_test.py b/projects/archive/check_mcmc_convergence_for_return_levels/gelman_convergence_test.py similarity index 100% rename from projects/others/check_mcmc_convergence_for_return_levels/gelman_convergence_test.py rename to projects/archive/check_mcmc_convergence_for_return_levels/gelman_convergence_test.py diff --git a/projects/others/check_mcmc_convergence_for_return_levels/main_bayesian_mcmc.py b/projects/archive/check_mcmc_convergence_for_return_levels/main_bayesian_mcmc.py similarity index 100% rename from projects/others/check_mcmc_convergence_for_return_levels/main_bayesian_mcmc.py rename to projects/archive/check_mcmc_convergence_for_return_levels/main_bayesian_mcmc.py diff --git a/projects/others/check_mcmc_convergence_for_return_levels/main_gelman_convergence_test.py b/projects/archive/check_mcmc_convergence_for_return_levels/main_gelman_convergence_test.py similarity index 100% rename from projects/others/check_mcmc_convergence_for_return_levels/main_gelman_convergence_test.py rename to projects/archive/check_mcmc_convergence_for_return_levels/main_gelman_convergence_test.py diff --git a/projects/ogorman/__init__.py b/projects/archive/ogorman/__init__.py similarity index 100% rename from projects/ogorman/__init__.py rename to projects/archive/ogorman/__init__.py diff --git a/projects/ogorman/gorman_figures/__init__.py b/projects/archive/ogorman/gorman_figures/__init__.py similarity index 100% rename from projects/ogorman/gorman_figures/__init__.py rename to projects/archive/ogorman/gorman_figures/__init__.py diff --git a/projects/ogorman/gorman_figures/daily_snowfall_fraction.py b/projects/archive/ogorman/gorman_figures/daily_snowfall_fraction.py similarity index 100% rename from projects/ogorman/gorman_figures/daily_snowfall_fraction.py rename to projects/archive/ogorman/gorman_figures/daily_snowfall_fraction.py diff --git a/projects/ogorman/gorman_figures/figure1/__init__.py b/projects/archive/ogorman/gorman_figures/figure1/__init__.py similarity index 100% rename from projects/ogorman/gorman_figures/figure1/__init__.py rename to projects/archive/ogorman/gorman_figures/figure1/__init__.py diff --git a/projects/ogorman/gorman_figures/figure1/comparative_curve_wrt_altitude.py b/projects/archive/ogorman/gorman_figures/figure1/comparative_curve_wrt_altitude.py similarity index 97% rename from projects/ogorman/gorman_figures/figure1/comparative_curve_wrt_altitude.py rename to projects/archive/ogorman/gorman_figures/figure1/comparative_curve_wrt_altitude.py index 475494668a9c9def12c8ce6525bbb0fe681ba89a..6986f685e7e444208cf26013929d3e84b90d0004 100644 --- a/projects/ogorman/gorman_figures/figure1/comparative_curve_wrt_altitude.py +++ b/projects/archive/ogorman/gorman_figures/figure1/comparative_curve_wrt_altitude.py @@ -6,7 +6,7 @@ import numpy as np from cached_property import cached_property -from projects.ogorman.gorman_figures.figure1.study_visualizer_for_double_stationary_fit import \ +from projects.archive.ogorman.gorman_figures import \ StudyVisualizerForReturnLevelChange diff --git a/projects/ogorman/gorman_figures/figure1/figure1_mean_ratio_return_level_ratio.py b/projects/archive/ogorman/gorman_figures/figure1/figure1_mean_ratio_return_level_ratio.py similarity index 92% rename from projects/ogorman/gorman_figures/figure1/figure1_mean_ratio_return_level_ratio.py rename to projects/archive/ogorman/gorman_figures/figure1/figure1_mean_ratio_return_level_ratio.py index 51991de1f88a0fe7ca5565cabde972481a50aab3..7fcd9aef2169b09b0c23f70bc14c07729fda4f74 100644 --- a/projects/ogorman/gorman_figures/figure1/figure1_mean_ratio_return_level_ratio.py +++ b/projects/archive/ogorman/gorman_figures/figure1/figure1_mean_ratio_return_level_ratio.py @@ -5,10 +5,10 @@ from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranS from extreme_data.meteo_france_data.scm_models_data.visualization.main_study_visualizer import ALL_ALTITUDES_WITHOUT_NAN from extreme_fit.model.margin_model.utils import \ MarginFitMethod -from projects.ogorman.gorman_figures.figure1 import \ +from projects.archive.ogorman.gorman_figures import \ ComparativeCurveWrtAltitude -from projects.ogorman.gorman_figures.figure1.study_visualizer_for_double_stationary_fit import \ +from projects.archive.ogorman.gorman_figures import \ StudyVisualizerForReturnLevelChange diff --git a/projects/ogorman/gorman_figures/figure1/result_from_stationary_fit.py b/projects/archive/ogorman/gorman_figures/figure1/result_from_stationary_fit.py similarity index 100% rename from projects/ogorman/gorman_figures/figure1/result_from_stationary_fit.py rename to projects/archive/ogorman/gorman_figures/figure1/result_from_stationary_fit.py diff --git a/projects/ogorman/gorman_figures/figure1/study_visualizer_for_double_stationary_fit.py b/projects/archive/ogorman/gorman_figures/figure1/study_visualizer_for_double_stationary_fit.py similarity index 92% rename from projects/ogorman/gorman_figures/figure1/study_visualizer_for_double_stationary_fit.py rename to projects/archive/ogorman/gorman_figures/figure1/study_visualizer_for_double_stationary_fit.py index 3fd5d2bd5cda20b4087f1f0d9dbe3d1aa55226d0..cb66916c0f11ca8fb703f065e192d6d204d9b502 100644 --- a/projects/ogorman/gorman_figures/figure1/study_visualizer_for_double_stationary_fit.py +++ b/projects/archive/ogorman/gorman_figures/figure1/study_visualizer_for_double_stationary_fit.py @@ -1,15 +1,12 @@ import matplotlib -from extreme_data.meteo_france_data.scm_models_data.abstract_study import AbstractStudy -from extreme_data.meteo_france_data.scm_models_data.visualization.create_shifted_cmap import get_shifted_map, \ - get_colors, ticks_values_and_labels_for_percentages from extreme_data.meteo_france_data.scm_models_data.visualization.study_visualizer import StudyVisualizer from extreme_fit.model.margin_model.utils import \ - MarginFitMethod, fitmethod_to_str + MarginFitMethod import matplotlib.pyplot as plt -from projects.ogorman.gorman_figures.figure1.result_from_stationary_fit import ResultFromDoubleStationaryFit +from projects.archive.ogorman.gorman_figures.figure1.result_from_stationary_fit import ResultFromDoubleStationaryFit class StudyVisualizerForReturnLevelChange(StudyVisualizer): diff --git a/projects/ogorman/gorman_figures/temperature_of_snowfall_maxima.py b/projects/archive/ogorman/gorman_figures/temperature_of_snowfall_maxima.py similarity index 97% rename from projects/ogorman/gorman_figures/temperature_of_snowfall_maxima.py rename to projects/archive/ogorman/gorman_figures/temperature_of_snowfall_maxima.py index bd518ce869a5a9843d8a08860b2693a5cc87fdd8..b1cfb28f85fbc736fcfa6b300fd18b287cb05bcd 100644 --- a/projects/ogorman/gorman_figures/temperature_of_snowfall_maxima.py +++ b/projects/archive/ogorman/gorman_figures/temperature_of_snowfall_maxima.py @@ -11,7 +11,7 @@ import numpy as np from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day, SafranTemperature, \ SafranPrecipitation1Day -from projects.ogorman.gorman_figures.daily_snowfall_fraction import \ +from projects.archive.ogorman.gorman_figures import \ compute_daily_snowfall_fraction diff --git a/projects/ogorman/snowfall fraction of total precipitation/__init__.py b/projects/archive/ogorman/snowfall fraction of total precipitation/__init__.py similarity index 100% rename from projects/ogorman/snowfall fraction of total precipitation/__init__.py rename to projects/archive/ogorman/snowfall fraction of total precipitation/__init__.py diff --git a/projects/ogorman/snowfall fraction of total precipitation/short_snow_load_partition.py b/projects/archive/ogorman/snowfall fraction of total precipitation/short_snow_load_partition.py similarity index 100% rename from projects/ogorman/snowfall fraction of total precipitation/short_snow_load_partition.py rename to projects/archive/ogorman/snowfall fraction of total precipitation/short_snow_load_partition.py diff --git a/projects/others/__init__.py b/projects/archive/seasonal_analysis/__init__.py similarity index 100% rename from projects/others/__init__.py rename to projects/archive/seasonal_analysis/__init__.py diff --git a/projects/seasonal_analysis/main_season_repartition_of_maxima.py b/projects/archive/seasonal_analysis/main_season_repartition_of_maxima.py similarity index 100% rename from projects/seasonal_analysis/main_season_repartition_of_maxima.py rename to projects/archive/seasonal_analysis/main_season_repartition_of_maxima.py diff --git a/projects/others/check_mcmc_convergence_for_return_levels/__init__.py b/projects/archive/thesis_report/__init__.py similarity index 100% rename from projects/others/check_mcmc_convergence_for_return_levels/__init__.py rename to projects/archive/thesis_report/__init__.py diff --git a/projects/thesis_report/gev_plot.py b/projects/archive/thesis_report/gev_plot.py similarity index 100% rename from projects/thesis_report/gev_plot.py rename to projects/archive/thesis_report/gev_plot.py diff --git a/projects/thesis_report/simulation_for_quantile_gap.py b/projects/archive/thesis_report/simulation_for_quantile_gap.py similarity index 100% rename from projects/thesis_report/simulation_for_quantile_gap.py rename to projects/archive/thesis_report/simulation_for_quantile_gap.py diff --git a/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/data.py b/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/data.py deleted file mode 100644 index aff3b1848ffe9c1918ed2dd546360539a18b72fc..0000000000000000000000000000000000000000 --- a/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/data.py +++ /dev/null @@ -1,76 +0,0 @@ -from collections import OrderedDict - -import matplotlib.pyplot as plt - -from extreme_data.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal -from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day -from extreme_data.meteo_france_data.scm_models_data.visualization.main_study_visualizer import \ - study_iterator_global, SCM_STUDY_CLASS_TO_ABBREVIATION -from extreme_data.meteo_france_data.scm_models_data.visualization.study_visualizer import \ - StudyVisualizer -from projects.exceeding_snow_loads.utils import dpi_paper1_figure - - -def tuples_for_examples_paper1(examples_for_the_paper=True): - if examples_for_the_paper: - - marker_altitude_massif_name_for_paper1 = [ - ('mediumpurple', 900, 'Ubaye'), - ('darkmagenta', 1800, 'Vercors'), - # ('mediumpurple', 2700, 'Beaufortain'), - ] - else: - marker_altitude_massif_name_for_paper1 = [ - ('magenta', 600, 'Parpaillon'), - ('darkmagenta', 300, 'Devoluy'), - ('mediumpurple', 300, 'Aravis'), - ] - return marker_altitude_massif_name_for_paper1 - - -def max_graph_annual_maxima_poster(): - """ - We choose these massif because each represents a different eurocode region - we also choose them because they belong to a different climatic area - :return: - """ - save_to_file = False - study_class = SafranSnowfall1Day - - examples_for_the_paper = True - - ax = plt.gca() - if examples_for_the_paper: - # ax.set_ylim([0, 20]) - # ax.set_yticks(list(range(0, 21, 2))) - linewidth = 5 - else: - linewidth = 3 - - marker_altitude_massif_name_for_paper1 = tuples_for_examples_paper1(examples_for_the_paper) - - altitude_to_linestyle = OrderedDict() - first_altitude = 900 - second_altitude = 2100 - altitude_to_linestyle[first_altitude] = 'dashed' - altitude_to_linestyle[second_altitude] = 'dotted' - for altitude, linestyle in altitude_to_linestyle.items(): - for study in study_iterator_global([study_class], altitudes=[altitude]): - study_visualizer = StudyVisualizer(study, save_to_file=save_to_file, - verbose=True, - multiprocessing=True) - snow_abbreviation = SCM_STUDY_CLASS_TO_ABBREVIATION[study_class] - tight_pad = {'h_pad': 0.2} - snow_abbreviation = 'max ' + snow_abbreviation - for color, _, massif_name in marker_altitude_massif_name_for_paper1[::-1]: - last_plot = massif_name == 'Ubaye' and altitude == second_altitude - label = '{} massif at {}m'.format(massif_name, altitude) - study_visualizer.visualize_max_graphs_poster(massif_name, altitude, snow_abbreviation, color, label, - last_plot, ax, tight_pad=tight_pad, - dpi=dpi_paper1_figure, - linewidth=linewidth, - linestyle=linestyle) - - -if __name__ == '__main__': - max_graph_annual_maxima_poster() diff --git a/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/main_snowfall_article.py b/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/main_snowfall_article.py deleted file mode 100644 index 7d5c0100b544335d80f03ca07d2186fbaefc5261..0000000000000000000000000000000000000000 --- a/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/main_snowfall_article.py +++ /dev/null @@ -1,158 +0,0 @@ -from multiprocessing.pool import Pool - -import matplotlib as mpl -import numpy as np - -from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day, SafranPrecipitation1Day, \ - SafranPrecipitation3Days, SafranSnowfall3Days -from extreme_data.meteo_france_data.scm_models_data.utils import Season -from extreme_trend.abstract_gev_trend_test import AbstractGevTrendTest -from projects.contrasting_trends_in_snow_loads.article2_snowfall_versus_time_and_altitude.plot_selection_curves_paper2 import \ - plot_selection_curves_paper2 -from projects.contrasting_trends_in_snow_loads.article2_snowfall_versus_time_and_altitude.shape_plot import shape_plot -from projects.contrasting_trends_in_snow_loads.article2_snowfall_versus_time_and_altitude.snowfall_plot import \ - plot_snowfall_mean, plot_snowfall_change_mean -from projects.contrasting_trends_in_snow_loads.article2_snowfall_versus_time_and_altitude.study_visualizer_for_mean_values import \ - StudyVisualizerForMeanValues -from projects.contrasting_trends_in_snow_loads.article2_snowfall_versus_time_and_altitude.study_visualizer_for_mean_values_with_mean_aic import \ - StudyVisualizerForMeanValuesWithMeanAic -from projects.contrasting_trends_in_snow_loads.article2_snowfall_versus_time_and_altitude.validation_plot import \ - validation_plot -from projects.exceeding_snow_loads.section_results.plot_selection_curves import plot_selection_curves -from projects.exceeding_snow_loads.section_results.plot_trend_curves import plot_trend_map - -mpl.rcParams['text.usetex'] = True -mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}'] - -from extreme_data.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal -from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import \ - ConfidenceIntervalMethodFromExtremes -from extreme_trend.visualizers.study_visualizer_for_non_stationary_trends import \ - StudyVisualizerForNonStationaryTrends -from extreme_trend.visualizers.utils import load_altitude_to_visualizer -from projects.exceeding_snow_loads.section_results.plot_uncertainty_curves import plot_uncertainty_massifs -from projects.exceeding_snow_loads.utils import paper_study_classes, paper_altitudes -from root_utils import NB_CORES - -import matplotlib.pyplot as plt - - -def minor_result(altitude): - """Plot trends for a single altitude to be fast""" - visualizer = StudyVisualizerForNonStationaryTrends(CrocusSnowLoadTotal(altitude=altitude), multiprocessing=True, - ) - visualizer.plot_trends() - plt.show() - - -def compute_minimized_aic(visualizer): - if isinstance(visualizer, StudyVisualizerForMeanValuesWithMeanAic): - _ = visualizer.massif_name_and_trend_test_class_to_trend_test - else: - _ = visualizer.massif_name_to_trend_test_that_minimized_aic - return True - - -def intermediate_result(altitudes, massif_names=None, - model_subsets_for_uncertainty=None, - uncertainty_methods=None, - study_class=SafranSnowfall1Day, - multiprocessing=False, study_visualizer_class=StudyVisualizerForMeanValues, - season=Season.annual): - """ - Plot all the trends for all altitudes - And enable to plot uncertainty plot for some specific massif_names, uncertainty methods to be fast - :param altitudes: - :param massif_names: - :param non_stationary_uncertainty: - :param uncertainty_methods: - :param study_class: - :return: - """ - # Load altitude to visualizer - altitude_to_visualizer = load_altitude_to_visualizer(altitudes, massif_names, model_subsets_for_uncertainty, - study_class, uncertainty_methods, - study_visualizer_class=study_visualizer_class, - season=season) - # Load variable object efficiently - for v in altitude_to_visualizer.values(): - _ = v.study.year_to_variable_object - # Compute minimized value efficiently - visualizers = list(altitude_to_visualizer.values()) - if multiprocessing: - with Pool(4) as p: - _ = p.imap(compute_minimized_aic, visualizers) - else: - for visualizer in visualizers: - _ = compute_minimized_aic(visualizer) - # Aggregate the choice for the minimizer - aggregate(visualizers) - - # Plots - validation_plot(altitude_to_visualizer, order_derivative=0) - # validation_plot(altitude_to_visualizer, order_derivative=1) - plot_snowfall_mean(altitude_to_visualizer) - # plot_selection_curves_paper2(altitude_to_visualizer) - plot_snowfall_change_mean(altitude_to_visualizer) - shape_plot(altitude_to_visualizer) - - -def major_result(): - uncertainty_methods = [ConfidenceIntervalMethodFromExtremes.ci_mle][:] - # massif_names = ['Beaufortain', 'Vercors'] - massif_names = None - study_classes = [SafranSnowfall1Day, SafranPrecipitation1Day][::-1][:] - # study_classes = [SafranSnowfall3Days, SafranPrecipitation3Days][::-1] - model_subsets_for_uncertainty = None - altitudes = paper_altitudes - altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000] - # altitudes = [900, 1200, 1500, 1800][:2] - # altitudes = [1800, 2100, 2400, 2700][:3] - # altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000] - # altitudes = draft_altitudes - # for significance_level in [0.1, 0.05][]: - AbstractGevTrendTest.SIGNIFICANCE_LEVEL = 0.05 - study_visualizer_class = [StudyVisualizerForMeanValues, StudyVisualizerForMeanValuesWithMeanAic][0] - season = [Season.annual, Season.winter_extended][1] - for study_class in study_classes: - intermediate_result(altitudes, massif_names, model_subsets_for_uncertainty, - uncertainty_methods, study_class, multiprocessing=False, - study_visualizer_class=study_visualizer_class, - season=season) - - -def aggregate(visualizers): - visualizer = visualizers[0] - if not isinstance(visualizer, StudyVisualizerForMeanValuesWithMeanAic): - return - massif_names = set.union(*[set(v.massif_names) for v in visualizers]) - massif_names = list(massif_names) - trend_tests = visualizer.non_stationary_trend_test - - massif_name_to_trend_test_to_aic_list = {m: {t: [] for t in trend_tests} - for m in massif_names} - for v in visualizers: - for (m, t), t2 in v.massif_name_and_trend_test_class_to_trend_test.items(): - massif_name_to_trend_test_to_aic_list[m][t].append(t2.aic) - - massif_name_to_trend_test_with_minimial_mean_aic = {} - for m in massif_names: - trend_test_and_mean_aic = [(t, np.array(massif_name_to_trend_test_to_aic_list[m][t]).mean()) - for t in trend_tests] - sorted_l = sorted(trend_test_and_mean_aic, key=lambda e: e[1]) - trend_test_with_minimial_mean_aic = sorted_l[0][0] - massif_name_to_trend_test_with_minimial_mean_aic[m] = trend_test_with_minimial_mean_aic - print(massif_name_to_trend_test_with_minimial_mean_aic) - for v in visualizers: - v.massif_name_to_trend_test_with_minimial_mean_aic = {} - for m, t in massif_name_to_trend_test_with_minimial_mean_aic.items(): - if (m, t) in v.massif_name_and_trend_test_class_to_trend_test: - v.massif_name_to_trend_test_with_minimial_mean_aic[m] = v.massif_name_and_trend_test_class_to_trend_test[(m, t)] - - -if __name__ == '__main__': - major_result() - # intermediate_result(altitudes=[300], massif_names=None, - # uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes, - # ConfidenceIntervalMethodFromExtremes.ci_mle][1:], - # multiprocessing=True) diff --git a/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/plot_selection_curves_paper2.py b/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/plot_selection_curves_paper2.py deleted file mode 100644 index 8e54d150acea0d12fb7a304d39439ea60fb0bb45..0000000000000000000000000000000000000000 --- a/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/plot_selection_curves_paper2.py +++ /dev/null @@ -1,73 +0,0 @@ -from typing import Dict -import matplotlib.pyplot as plt - -from extreme_data.meteo_france_data.scm_models_data.visualization.utils import create_adjusted_axes -from extreme_trend.abstract_gev_trend_test import AbstractGevTrendTest -from projects.contrasting_trends_in_snow_loads.article2_snowfall_versus_time_and_altitude.study_visualizer_for_mean_values import \ - StudyVisualizerForMeanValues -from projects.exceeding_snow_loads.section_results.plot_selection_curves import merge_counter, \ - get_ordered_list_from_counter, permute -from projects.exceeding_snow_loads.utils import dpi_paper1_figure, get_trend_test_name -from extreme_trend.visualizers.study_visualizer_for_non_stationary_trends import StudyVisualizerForNonStationaryTrends - - -def plot_selection_curves_paper2(altitude_to_visualizer: Dict[int, StudyVisualizerForMeanValues]): - visualizer = list(altitude_to_visualizer.values())[0] - - ax = create_adjusted_axes(1, 1) - - vs = [v for v in altitude_to_visualizer.values()] - - selected_counter = merge_counter([v.selected_trend_test_class_counter for v in vs]) - selected_and_anderson_counter = merge_counter([v.selected_and_anderson_trend_test_class_counter for v in vs]) - selected_and_anderson_and_likelihood_counter = merge_counter( - [v.selected_and_anderson_and_likelihood_ratio_trend_test_class_counter() for v in vs]) - - total_of_selected_models = sum(selected_counter.values()) - l = sorted(enumerate(visualizer.non_stationary_trend_test), key=lambda e: selected_counter[e[1]]) - permutation = [i for i, v in l][::-1] - - select_list = get_ordered_list_from_counter(selected_counter, total_of_selected_models, visualizer, permutation) - selected_and_anderson_list = get_ordered_list_from_counter(selected_and_anderson_counter, total_of_selected_models, - visualizer, permutation) - selected_and_anderson_and_likelihood_list = get_ordered_list_from_counter( - selected_and_anderson_and_likelihood_counter, total_of_selected_models, visualizer, permutation) - - labels = [get_trend_test_name(t) for t in visualizer.non_stationary_trend_test] - labels = permute(labels, permutation) - print(select_list, sum(select_list)) - - nb_selected_models = sum(select_list) - nb_selected_and_anderson_models = sum(selected_and_anderson_list) - nb_selected_and_anderson_and_likelihood_models = sum(selected_and_anderson_and_likelihood_list) - nb_selected_models_not_passing_any_test = nb_selected_models - nb_selected_and_anderson_models - nb_selected_models_just_passing_anderson = nb_selected_and_anderson_models - nb_selected_and_anderson_and_likelihood_models - - # parameters - width = 5 - size = 30 - legend_fontsize = 30 - labelsize = 15 - linewidth = 3 - x = [10 * (i + 1) for i in range(len(select_list))] - for l, color, label, nb in zip([select_list, selected_and_anderson_list, selected_and_anderson_and_likelihood_list], - ['white', 'grey', 'black'], - ['Not passing any test', 'Just passing anderson test ', - 'Passing both anderson and likelihood tests '], - [nb_selected_models_not_passing_any_test, nb_selected_models_just_passing_anderson, nb_selected_and_anderson_and_likelihood_models]): - label += ' ({} \%)'.format(round(nb, 1)) - ax.bar(x, l, width=width, color=color, edgecolor='black', label=label, linewidth=linewidth) - - ax.legend(loc='upper right', prop={'size': size}) - ax.set_ylabel(' Frequency of selected model w.r.t all time series \n ' - 'i.e. for all massifs and altitudes (\%)', fontsize=legend_fontsize) - ax.set_xlabel('Models', fontsize=legend_fontsize) - ax.tick_params(axis='both', which='major', labelsize=labelsize) - ax.set_xticks(x) - ax.yaxis.grid() - ax.set_xticklabels(labels) - - # Save plot - visualizer.plot_name = 'Selection curves with significance level = {} '.format(AbstractGevTrendTest.SIGNIFICANCE_LEVEL) - visualizer.show_or_save_to_file(no_title=True, dpi=dpi_paper1_figure) - plt.close() diff --git a/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/shape_plot.py b/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/shape_plot.py deleted file mode 100644 index 657bf026d53dfd0f9596d5cd95fd38fe6caf6a09..0000000000000000000000000000000000000000 --- a/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/shape_plot.py +++ /dev/null @@ -1,17 +0,0 @@ -from typing import Dict - -import matplotlib -import matplotlib.pyplot as plt - -from projects.contrasting_trends_in_snow_loads.article2_snowfall_versus_time_and_altitude.study_visualizer_for_mean_values import \ - StudyVisualizerForMeanValues - - -def shape_plot(altitude_to_visualizer: Dict[int, StudyVisualizerForMeanValues]): - # Plot map for the repartition of the difference - for altitude, visualizer in altitude_to_visualizer.items(): - label = ' shape parameter' - visualizer.plot_abstract_fast(massif_name_to_value=visualizer.massif_name_to_model_shape_last_year, - label='Model' + label, negative_and_positive_values=True, add_text=True, - cmap=matplotlib.cm.get_cmap('BrBG_r'), graduation=0.1) - plt.close() diff --git a/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/snowfall_plot.py b/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/snowfall_plot.py deleted file mode 100644 index 03b9c51e3a39ffd2d447ccb80bcdad0735cb4998..0000000000000000000000000000000000000000 --- a/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/snowfall_plot.py +++ /dev/null @@ -1,176 +0,0 @@ -from typing import Dict -import matplotlib.pyplot as plt - -import numpy as np -from sklearn.linear_model import LinearRegression - -from extreme_data.meteo_france_data.scm_models_data.visualization.main_study_visualizer import \ - SCM_STUDY_CLASS_TO_ABBREVIATION -from extreme_data.meteo_france_data.scm_models_data.visualization.plot_utils import plot_against_altitude -from projects.contrasting_trends_in_snow_loads.article2_snowfall_versus_time_and_altitude.study_visualizer_for_mean_values import \ - StudyVisualizerForMeanValues - - -def fit_linear_regression(x, y): - X = np.array(x).reshape(-1, 1) - reg = LinearRegression().fit(X, y) - r2_score = reg.score(X, y) - a = reg.coef_ - b = reg.intercept_ - return a, b, r2_score - - -def plot_snowfall_change_mean(altitude_to_visualizer: Dict[int, StudyVisualizerForMeanValues]): - visualizer = list(altitude_to_visualizer.values())[0] - study = visualizer.study - - variables = ['mean annual maxima', '50 year return level'] - mean_indicators = [True, False] - l = list(zip(variables, mean_indicators)) - for variable, mean_indicator in l[:1]: - for relative in [False, True]: - if relative: - variable += str(' (relative)') - # Plot the curve for the evolution of the mean - massif_name_to_a, massif_name_to_b, massif_name_to_r2_score = plot_mean(altitude_to_visualizer, - derivative=True, mean_indicator=mean_indicator, - relative=relative) - # Augmentation every km - massif_name_to_augmentation_every_km = {m: a * 1000 for m, a in massif_name_to_a.items()} - massif_to_r2_score_text = {k: str(round(v, 2)) for k, v in massif_name_to_r2_score.items()} - - visualizer.plot_abstract_fast(massif_name_to_augmentation_every_km, - label='Slope for changes of {} of {}\n for every km of elevation ({})'.format( - variable, SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)], study.variable_unit), - add_x_label=False, - add_text=True, massif_name_to_text=massif_to_r2_score_text) - # Value at 2000 m - massif_name_to_mean_at_2000 = {m: a * 2000 + massif_name_to_b[m] for m, a in massif_name_to_a.items()} - visualizer.plot_abstract_fast(massif_name_to_mean_at_2000, - label='Changes of {} \nof {} at 2000 m ({})'.format( - variable, SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)], study.variable_unit), - add_x_label=False, - add_text=True, massif_name_to_text=massif_to_r2_score_text) - # Altitude for the change of dynamic - massif_name_to_altitude_change_dynamic = {m: - massif_name_to_b[m] / a for m, a in massif_name_to_a.items()} - # Keep only those that are in a reasonable range - massif_name_to_altitude_change_dynamic = {m: d for m, d in massif_name_to_altitude_change_dynamic.items() - if 0 < d < 3000} - visualizer.plot_abstract_fast(massif_name_to_altitude_change_dynamic, - label=('Altitude when the changes equal zero for %s (m)' % variable), - add_x_label=False, graduation=500, - add_text=True, massif_name_to_text=massif_to_r2_score_text) - - -def plot_snowfall_mean(altitude_to_visualizer: Dict[int, StudyVisualizerForMeanValues]): - visualizer = list(altitude_to_visualizer.values())[0] - study = visualizer.study - - variables = ['mean annual maxima', '50 year return level'] - mean_indicators = [True, False] - for variable, mean_indicator in zip(variables, mean_indicators): - # Plot the curve for the evolution of the mean - massif_name_to_a, massif_name_to_b, massif_name_to_r2_score = plot_mean(altitude_to_visualizer, - derivative=False, - mean_indicator=mean_indicator) - # Compute values of interest - massif_name_to_mean_at_2000 = {m: a * 2000 + massif_name_to_b[m] for m, a in massif_name_to_a.items()} - massif_name_to_augmentation_every_km = {m: a * 1000 for m, a in massif_name_to_a.items()} - massif_name_to_relative_augmentation_between_2000_and_3000_every_km = { - m: 100 * (a * 1000) / (a * 2000 + massif_name_to_b[m]) for m, a in massif_name_to_a.items()} - massif_name_to_relative_augmentation_between_1000_and_2000_every_km = { - m: 100 * (a * 1000) / (a * 1000 + massif_name_to_b[m]) for m, a in massif_name_to_a.items()} - massif_to_r2_score_text = {k: str(round(v, 2)) for k, v in massif_name_to_r2_score.items()} - # Augmentation every km - ds = [massif_name_to_augmentation_every_km, massif_name_to_relative_augmentation_between_1000_and_2000_every_km, - massif_name_to_relative_augmentation_between_2000_and_3000_every_km] - prefixs = ['Augmentation', 'Relative augmentation between 1000m and 2000m', - 'Relative augmentation between 2000m and 3000m'] - graduations = [1.0, 10.0, 10.0] - for d, prefix, graduation in zip(ds, prefixs, graduations): - visualizer.plot_abstract_fast(d, - graduation=10.0, - label=prefix + ' of {} of {} \nfor every km of elevation ({})'.format( - variable, SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)], - study.variable_unit), - add_x_label=False, negative_and_positive_values=False, - add_text=True, - massif_name_to_text=massif_to_r2_score_text - ) - # Value at 2000 m - visualizer.plot_abstract_fast(massif_name_to_mean_at_2000, label='{} of {} at 2000 m ()'.format( - variable, SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)], study.variable_unit), - add_x_label=False, negative_and_positive_values=False, - add_text=True, massif_name_to_text=massif_to_r2_score_text) - - -def plot_mean(altitude_to_visualizer: Dict[int, StudyVisualizerForMeanValues], derivative=False, mean_indicator=True, - relative=False): - ax = plt.gca() - massif_name_to_linear_regression_result = {} - - visualizers = list(altitude_to_visualizer.values()) - visualizer = visualizers[0] - study = visualizer.study - - for massif_id, massif_name in enumerate(visualizer.study.all_massif_names()): - altitudes_massif = [a for a, v in altitude_to_visualizer.items() - if massif_name in v.massif_name_to_trend_test_that_minimized_aic] - if len(altitudes_massif) >= 2: - trend_tests = [altitude_to_visualizer[a].massif_name_to_trend_test_that_minimized_aic[massif_name] - for a in altitudes_massif] - nb_years = 50 - return_period = 50 - if mean_indicator: - indicator_str = 'mean' - else: - indicator_str = '50 year return level' - res = [(a, t) for i, (a, t) in enumerate(zip(altitudes_massif, trend_tests))] - if derivative: - if mean_indicator: - if relative: - res = [(a, t.relative_change_in_mean_for_the_last_x_years(nb_years=nb_years)) for a, t in res] - else: - res = [(a, t.change_in_mean_for_the_last_x_years(nb_years=nb_years)) for a, t in res] - else: - if relative: - res = [(a, t.relative_change_in_50_year_return_level_for_the_last_x_years(nb_years=nb_years, return_period=return_period)) for a, t in res] - else: - res = [(a, t.change_in_50_year_return_level_for_the_last_x_years(nb_years=nb_years, return_period=return_period)) for a, t in res] - if len(res) > 0: - altitudes_values, values = zip(*res) - else: - altitudes_values, values = [], [] - indicator_str = 'Change of the {} for the last {} years \nfor non-stationary models'.format(indicator_str, nb_years) - if relative: - indicator_str += ' (relative)' - else: - if mean_indicator: - values = [t.unconstrained_estimator_gev_params_last_year.mean for t in trend_tests] - else: - values = [t.unconstrained_estimator_gev_params_last_year.return_level(return_period=return_period) - for t in trend_tests] - altitudes_values = altitudes_massif - # Plot - if len(altitudes_values) >= 2: - massif_name_to_linear_regression_result[massif_name] = fit_linear_regression(altitudes_values, values) - plot_values_against_altitudes(ax, altitudes_values, massif_id, massif_name, indicator_str, study, values, - visualizer) - ax.legend(prop={'size': 7}, ncol=3) - visualizer.show_or_save_to_file(dpi=500, add_classic_title=False) - plt.close() - - return [{m: t[i][0] if i == 0 else t[i] - for m, t in massif_name_to_linear_regression_result.items()} for i in range(3)] - - -def plot_values_against_altitudes(ax, altitudes, massif_id, massif_name, moment, study, values, visualizer): - plot_against_altitude(x_ticks=altitudes, ax=ax, massif_id=massif_id, massif_name=massif_name, values=values) - plot_name = '{} annual maxima of {}'.format(moment, SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)]) - ax.set_ylabel('{} ({})'.format(plot_name, study.variable_unit), fontsize=15) - ax.set_xlabel('altitudes', fontsize=15) - # lim_down, lim_up = ax.get_ylim() - # lim_up += (lim_up - lim_down) / 3 - # ax.set_ylim([lim_down, lim_up]) - ax.tick_params(axis='both', which='major', labelsize=13) - visualizer.plot_name = plot_name diff --git a/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/study_visualizer_for_mean_values.py b/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/study_visualizer_for_mean_values.py deleted file mode 100644 index d1a30c5c84ab29232989af25125d48e11186c9a4..0000000000000000000000000000000000000000 --- a/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/study_visualizer_for_mean_values.py +++ /dev/null @@ -1,154 +0,0 @@ -from collections import Counter -from typing import Dict - -import matplotlib.pyplot as plt -import numpy as np -from cached_property import cached_property - -from extreme_data.eurocode_data.utils import YEAR_OF_INTEREST_FOR_RETURN_LEVEL -from extreme_data.meteo_france_data.scm_models_data.abstract_study import AbstractStudy -from extreme_fit.model.margin_model.utils import MarginFitMethod -from extreme_trend.abstract_gev_trend_test import AbstractGevTrendTest -from extreme_trend.visualizers.study_visualizer_for_non_stationary_trends import StudyVisualizerForNonStationaryTrends -from projects.exceeding_snow_loads.utils import NON_STATIONARY_TREND_TEST_PAPER_2 - - -class StudyVisualizerForMeanValues(StudyVisualizerForNonStationaryTrends): - - def __init__(self, study: AbstractStudy, show=True, save_to_file=False, only_one_graph=False, only_first_row=False, - vertical_kde_plot=False, year_for_kde_plot=None, plot_block_maxima_quantiles=False, - temporal_non_stationarity=False, transformation_class=None, verbose=False, multiprocessing=False, - complete_non_stationary_trend_analysis=False, normalization_under_one_observations=True, - uncertainty_methods=None, model_subsets_for_uncertainty=None, uncertainty_massif_names=None, - effective_temporal_covariate=YEAR_OF_INTEREST_FOR_RETURN_LEVEL, relative_change_trend_plot=True, - non_stationary_trend_test_to_marker=None, fit_method=MarginFitMethod.extremes_fevd_mle, - select_only_acceptable_shape_parameter=True, fit_gev_only_on_non_null_maxima=False, - fit_only_time_series_with_ninety_percent_of_non_null_values=True): - - # Change the type of models used in the study - non_stationary_trend_test_to_marker = {c: '' for c in NON_STATIONARY_TREND_TEST_PAPER_2} - - super().__init__(study, show, save_to_file, only_one_graph, only_first_row, vertical_kde_plot, - year_for_kde_plot, plot_block_maxima_quantiles, temporal_non_stationarity, - transformation_class, verbose, multiprocessing, complete_non_stationary_trend_analysis, - normalization_under_one_observations, uncertainty_methods, model_subsets_for_uncertainty, - uncertainty_massif_names, effective_temporal_covariate, relative_change_trend_plot, - non_stationary_trend_test_to_marker, fit_method, select_only_acceptable_shape_parameter, - fit_gev_only_on_non_null_maxima, fit_only_time_series_with_ninety_percent_of_non_null_values) - assert len(self.non_stationary_trend_test) == len(NON_STATIONARY_TREND_TEST_PAPER_2) - - def plot_abstract_fast(self, massif_name_to_value, label, graduation=10.0, cmap=plt.cm.coolwarm, add_x_label=True, - negative_and_positive_values=True, add_text=False, massif_name_to_text=None): - if add_text: - if massif_name_to_text is None: - massif_name_to_text = self.massif_name_to_text - else: - massif_name_to_text = None - super().plot_abstract(massif_name_to_value, label, label, self.fit_method, graduation, cmap, add_x_label, - negative_and_positive_values, massif_name_to_text) - - @property - def massif_name_to_text(self): - d = {} - for m, t in self.massif_name_to_trend_test_that_minimized_aic.items(): - latex_command = 'textbf' if t.is_significant else 'textrm' - d[m] = '$\\' + latex_command + '{' + t.name + '}$' - return d - - # Override the main dict massif_name_to_trend_test_that_minimized_aic - - @property - def super_massif_name_to_trend_test_that_minimized_aic(self): - return super().massif_name_to_trend_test_that_minimized_aic - - @property - def massif_name_to_trend_test_that_minimized_aic_after_anderson_test(self): - return {m: t for m, t in super().massif_name_to_trend_test_that_minimized_aic.items() - if t.goodness_of_fit_anderson_test} - - @property - def massif_name_to_trend_test_that_minimized_aic_after_anderson_test_and_likelihood_ratio_test(self): - return {m: t for m, t in self.massif_name_to_trend_test_that_minimized_aic_after_anderson_test.items() - if t.is_significant} - - @property - def massif_name_to_trend_test_that_minimized_aic(self) -> Dict[str, AbstractGevTrendTest]: - return self.massif_name_to_trend_test_that_minimized_aic_after_anderson_test - - # Counter for the selection process - - @cached_property - def selected_trend_test_class_counter(self): - return Counter([type(t) for t in self.super_massif_name_to_trend_test_that_minimized_aic.values()]) - - @cached_property - def selected_and_anderson_trend_test_class_counter(self): - return Counter([type(t) for t in self.massif_name_to_trend_test_that_minimized_aic_after_anderson_test.values()]) - - def selected_and_anderson_and_likelihood_ratio_trend_test_class_counter(self): - return Counter([type(t) for t in self.massif_name_to_trend_test_that_minimized_aic_after_anderson_test_and_likelihood_ratio_test.values()]) - -# Study of the mean - - @cached_property - def massif_name_to_empirical_mean(self): - massif_name_to_empirical_value = {} - for massif_name, maxima in self.study.massif_name_to_annual_maxima.items(): - massif_name_to_empirical_value[massif_name] = np.mean(maxima) - return massif_name_to_empirical_value - - @cached_property - def massif_name_to_model_shape_last_year(self): - massif_name_to_model_value_last_year = {} - for massif_name, trend_test in self.massif_name_to_trend_test_that_minimized_aic.items(): - massif_name_to_model_value_last_year[ - massif_name] = trend_test.unconstrained_estimator_gev_params_last_year.shape - return massif_name_to_model_value_last_year - - @cached_property - def massif_name_to_model_mean_last_year(self): - massif_name_to_model_value_last_year = {} - for massif_name, trend_test in self.massif_name_to_trend_test_that_minimized_aic.items(): - massif_name_to_model_value_last_year[ - massif_name] = trend_test.unconstrained_estimator_gev_params_last_year.mean - return massif_name_to_model_value_last_year - - @cached_property - def massif_name_to_relative_difference_for_mean(self): - massif_name_to_relative_difference = {} - for massif_name in self.massif_name_to_trend_test_that_minimized_aic.keys(): - e = self.massif_name_to_empirical_mean[massif_name] - m = self.massif_name_to_model_mean_last_year[massif_name] - relative_diference = 100 * (m - e) / e - massif_name_to_relative_difference[massif_name] = relative_diference - return massif_name_to_relative_difference - - # Study of the change in the mean - - @cached_property - def massif_name_to_change_ratio_in_empirical_mean(self): - massif_name_to_empirical_value = {} - for massif_name, maxima in self.study.massif_name_to_annual_maxima.items(): - index = self.study.ordered_years.index(1989) - maxima_before, maxima_after = maxima[:index + 1], maxima[index + 1:] - massif_name_to_empirical_value[massif_name] = np.mean(maxima_after) / np.mean(maxima_before) - return massif_name_to_empirical_value - - @cached_property - def massif_name_to_change_ratio_in_model_mean(self): - massif_name_to_parameter_value = {} - for massif_name, trend_test in self.massif_name_to_trend_test_that_minimized_aic.items(): - model_mean_before = trend_test.unconstrained_average_mean_value(year_min=self.study.year_min, year_max=1989) - model_mean_after = trend_test.unconstrained_average_mean_value(year_min=1990, year_max=self.study.year_max) - massif_name_to_parameter_value[massif_name] = model_mean_after / model_mean_before - return massif_name_to_parameter_value - - @cached_property - def massif_name_to_relative_difference_for_change_ratio_in_mean(self): - massif_name_to_relative_difference = {} - for massif_name in self.massif_name_to_trend_test_that_minimized_aic.keys(): - e = self.massif_name_to_change_ratio_in_empirical_mean[massif_name] - m = self.massif_name_to_change_ratio_in_model_mean[massif_name] - relative_diference = 100 * (m - e) / e - massif_name_to_relative_difference[massif_name] = relative_diference - return massif_name_to_relative_difference diff --git a/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/study_visualizer_for_mean_values_with_mean_aic.py b/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/study_visualizer_for_mean_values_with_mean_aic.py deleted file mode 100644 index 7317f7f6f3c65b3de5b0c89d76e1d9fc663c8421..0000000000000000000000000000000000000000 --- a/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/study_visualizer_for_mean_values_with_mean_aic.py +++ /dev/null @@ -1,57 +0,0 @@ -from collections import Counter -from typing import Dict - -import matplotlib.pyplot as plt -import numpy as np -from cached_property import cached_property - -from extreme_data.eurocode_data.utils import YEAR_OF_INTEREST_FOR_RETURN_LEVEL -from extreme_data.meteo_france_data.scm_models_data.abstract_study import AbstractStudy -from extreme_fit.model.margin_model.utils import MarginFitMethod -from extreme_trend.abstract_gev_trend_test import AbstractGevTrendTest -from extreme_trend.visualizers.study_visualizer_for_non_stationary_trends import StudyVisualizerForNonStationaryTrends -from projects.contrasting_trends_in_snow_loads.article2_snowfall_versus_time_and_altitude.study_visualizer_for_mean_values import \ - StudyVisualizerForMeanValues -from projects.exceeding_snow_loads.utils import NON_STATIONARY_TREND_TEST_PAPER_2 - - -class StudyVisualizerForMeanValuesWithMeanAic(StudyVisualizerForMeanValues): - - def __init__(self, study: AbstractStudy, show=True, save_to_file=False, only_one_graph=False, only_first_row=False, - vertical_kde_plot=False, year_for_kde_plot=None, plot_block_maxima_quantiles=False, - temporal_non_stationarity=False, transformation_class=None, verbose=False, multiprocessing=False, - complete_non_stationary_trend_analysis=False, normalization_under_one_observations=True, - uncertainty_methods=None, model_subsets_for_uncertainty=None, uncertainty_massif_names=None, - effective_temporal_covariate=YEAR_OF_INTEREST_FOR_RETURN_LEVEL, relative_change_trend_plot=True, - non_stationary_trend_test_to_marker=None, fit_method=MarginFitMethod.extremes_fevd_mle, - select_only_acceptable_shape_parameter=True, fit_gev_only_on_non_null_maxima=False, - fit_only_time_series_with_ninety_percent_of_non_null_values=True): - super().__init__(study, show, save_to_file, only_one_graph, only_first_row, vertical_kde_plot, - year_for_kde_plot, plot_block_maxima_quantiles, temporal_non_stationarity, - transformation_class, verbose, multiprocessing, complete_non_stationary_trend_analysis, - normalization_under_one_observations, uncertainty_methods, model_subsets_for_uncertainty, - uncertainty_massif_names, effective_temporal_covariate, relative_change_trend_plot, - non_stationary_trend_test_to_marker, fit_method, select_only_acceptable_shape_parameter, - fit_gev_only_on_non_null_maxima, fit_only_time_series_with_ninety_percent_of_non_null_values) - self.massif_name_to_trend_test_with_minimial_mean_aic = None - - @property - def massif_name_to_trend_test_that_minimized_aic(self) -> Dict[str, AbstractGevTrendTest]: - if self.massif_name_to_trend_test_with_minimial_mean_aic is None: - raise NotImplementedError('Aggregation must be run first') - else: - return self.massif_name_to_trend_test_with_minimial_mean_aic - - @property - def massif_names(self): - return [m for m, _ in self.massif_name_and_trend_test_class_to_trend_test.keys()] - - @cached_property - def massif_name_and_trend_test_class_to_trend_test(self): - d = {} - for massif_name in self.massif_name_to_years_and_maxima_for_model_fitting.keys(): - trend_tests = self.get_sorted_trend_test(massif_name) - for trend_test in trend_tests: - d[(massif_name, type(trend_test))] = trend_test - return d - diff --git a/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/validation_plot.py b/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/validation_plot.py deleted file mode 100644 index 7019880de6b644347763d4b6548b8f954c290513..0000000000000000000000000000000000000000 --- a/projects/contrasting_trends_in_snow_loads/article2_snowfall_versus_time_and_altitude/validation_plot.py +++ /dev/null @@ -1,79 +0,0 @@ -from multiprocessing import Pool -from typing import Dict - -import matplotlib.pyplot as plt - -from extreme_data.meteo_france_data.scm_models_data.visualization.main_study_visualizer import \ - SCM_STUDY_CLASS_TO_ABBREVIATION -from projects.contrasting_trends_in_snow_loads.article2_snowfall_versus_time_and_altitude.study_visualizer_for_mean_values import \ - StudyVisualizerForMeanValues -from projects.ogorman.gorman_figures.figure1.study_visualizer_for_double_stationary_fit import \ - StudyVisualizerForReturnLevelChange -from root_utils import NB_CORES - - -def validation_plot(altitude_to_visualizer: Dict[int, StudyVisualizerForMeanValues], order_derivative=0): - # Plot the mean empirical, the mean parametric and the relative difference between the two - altitudes = list(altitude_to_visualizer.keys()) - study_visualizer = list(altitude_to_visualizer.values())[0] - visualizers = list(altitude_to_visualizer.values()) - if order_derivative == 0: - plot_function = plot_relative_difference_map_order_zero - else: - plot_function = plot_relative_difference_map_order_one - # Plot map for the repartition of the difference - altitude_to_relative_differences = {} - for altitude, visualizer in altitude_to_visualizer.items(): - altitude_to_relative_differences[altitude] = plot_function(visualizer) - # study_visualizer.show_or_save_to_file(add_classic_title=False, dpi=500) - # # Shoe plot with respect to the altitude. - # plot_shoe_relative_differences_distribution(altitude_to_relative_differences, altitudes, study_visualizer, - # order_derivative) - # study_visualizer.show_or_save_to_file(add_classic_title=False, dpi=500) - plt.close() - - -def plot_shoe_relative_differences_distribution(altitude_to_relative_differences, altitudes, visualizer, - order_derivative): - study = visualizer.study - ax = plt.gca() - width = 150 - ax.boxplot([altitude_to_relative_differences[a] for a in altitudes], positions=altitudes, widths=width) - ax.set_xlim([min(altitudes) - width, max(altitudes) + width]) - moment = '' if order_derivative == 0 else 'time derivative of ' - ylabel = 'Global relative difference of the {} model mean \n' \ - 'w.r.t. the {}empirical mean of {} (\%)'.format(moment, moment, - SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)]) - visualizer.plot_trends = ylabel - ax.set_ylabel(ylabel) - ax.set_xlabel('Altitude (m)') - ax.legend() - ax.grid() - - -def plot_relative_difference_map_order_zero(visualizer: StudyVisualizerForMeanValues): - study = visualizer.study - label = ' mean annual maxima of {} ({})'.format(SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)], study.variable_unit) - # visualizer.plot_abstract_fast(massif_name_to_value=visualizer.massif_name_to_empirical_mean, - # label='Empirical' + label, negative_and_positive_values=False) - visualizer.plot_abstract_fast(massif_name_to_value=visualizer.massif_name_to_model_mean_last_year, - label='Model' + label, negative_and_positive_values=False, add_text=True) - # visualizer.plot_abstract_fast(massif_name_to_value=visualizer.massif_name_to_relative_difference_for_mean, - # label='Relative difference of the model mean w.r.t. the empirical mean \n' - # 'for the ' + label, graduation=1) - return list(visualizer.massif_name_to_relative_difference_for_mean.values()) - - -def plot_relative_difference_map_order_one(visualizer: StudyVisualizerForMeanValues): - study = visualizer.study - label = ' time derivative of mean annual maxima of {} ({})'.format(SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)], - study.variable_unit) - # visualizer.plot_abstract_fast(massif_name_to_value=visualizer.massif_name_to_change_ratio_in_empirical_mean, - # label='Empirical' + label, negative_and_positive_values=False, graduation=0.5) - visualizer.plot_abstract_fast(massif_name_to_value=visualizer.massif_name_to_change_ratio_in_model_mean, - label='Model' + label, negative_and_positive_values=False, graduation=0.5) - # visualizer.plot_abstract_fast( - # massif_name_to_value=visualizer.massif_name_to_relative_difference_for_change_ratio_in_mean, - # label='Relative difference of the model mean w.r.t. the empirical mean \n' - # 'for the ' + label, graduation=5) - return list(visualizer.massif_name_to_relative_difference_for_change_ratio_in_mean.values()) diff --git a/projects/contrasting_trends_in_snow_loads/post thesis committee/pointwise_gev_analysis_cluster_level.py b/projects/contrasting_trends_in_snow_loads/post thesis committee/pointwise_gev_analysis_cluster_level.py deleted file mode 100644 index 63457fb441b536dd809905057688d35a2bbbd521..0000000000000000000000000000000000000000 --- a/projects/contrasting_trends_in_snow_loads/post thesis committee/pointwise_gev_analysis_cluster_level.py +++ /dev/null @@ -1,152 +0,0 @@ -import numpy as np -from cached_property import cached_property -import matplotlib.pyplot as plt - -from extreme_data.meteo_france_data.scm_models_data.cluster.clustering_total_precip import TotalPrecipCluster -from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day -from extreme_fit.distribution.gev.gev_params import GevParams -from extreme_fit.estimator.margin_estimator.utils import fitted_stationary_gev -from extreme_fit.model.margin_model.utils import MarginFitMethod -from projects.altitude_spatial_model.altitudes_fit.altitudes_studies import AltitudesStudies - - -class PointWIseGevAnalysisForCluster(AltitudesStudies): - - def __init__(self, study_class, altitudes, - spatial_transformation_class=None, temporal_transformation_class=None, - cluster_class=TotalPrecipCluster, - **kwargs_study): - super().__init__(study_class, altitudes, spatial_transformation_class, temporal_transformation_class, - **kwargs_study) - self.cluster_class = cluster_class - - # Plot for the Altitude - - def cluster_id_to_annual_maxima_list(self, first_year=1959, last_year=2019, altitudes=None): - if altitudes is None: - altitudes = self.altitudes - cluster_id_to_annual_maxima_list = {} - for cluster_id in self.cluster_class.cluster_ids: - annual_maxima_list = [] - massif_names = self.cluster_class.cluster_id_to_massif_names[cluster_id] - for study in [s for a, s in self.altitude_to_study.items() if a in altitudes]: - annual_maxima = [] - massif_ids = [study.massif_name_to_massif_id[m] for m in massif_names if - m in study.massif_name_to_massif_id] - for year in range(first_year, last_year + 1): - annual_maxima.extend(study.year_to_annual_maxima[year][massif_ids]) - annual_maxima_list.append(annual_maxima) - cluster_id_to_annual_maxima_list[cluster_id] = annual_maxima_list - return cluster_id_to_annual_maxima_list - - @property - def cluster_id_to_gev_params_list(self): - cluster_id_to_gev_params_list = {} - for cluster_id, annual_maxima_list in self.cluster_id_to_annual_maxima_list().items(): - gev_params_list = [] - for annual_maxima in annual_maxima_list: - if len(annual_maxima) > 0: - gev_params = fitted_stationary_gev(annual_maxima, fit_method=MarginFitMethod.extremes_fevd_mle) - else: - gev_params = GevParams(0, 0, 0) - assert gev_params.has_undefined_parameters - gev_params_list.append(gev_params) - cluster_id_to_gev_params_list[cluster_id] = gev_params_list - return cluster_id_to_gev_params_list - - def plot_gev_parameters_against_altitude(self): - for param_name in GevParams.PARAM_NAMES[:]: - ax = plt.gca() - for cluster_id, gev_param_list in self.cluster_id_to_gev_params_list.items(): - params, altitudes = zip(*[(gev_param.to_dict()[param_name], altitude) for gev_param, altitude - in zip(gev_param_list, self.altitudes) if - not gev_param.has_undefined_parameters]) - # ax.plot(self.altitudes, params, label=str(cluster_id), linestyle='None', marker='x') - label = self.cluster_class.cluster_id_to_cluster_name[cluster_id] - ax.plot(altitudes, params, label=label, marker='x') - ax.legend() - ax.set_xlabel('Altitude') - ax.set_ylabel(param_name) - plot_name = '{} change with altitude'.format(param_name) - self.show_or_save_to_file(plot_name, no_title=True, tight_layout=True, show=False) - ax.clear() - plt.close() - - # Plot for the time - - @property - def year_min_and_max_list(self): - l = [] - year_min, year_max = 1959, 1989 - for shift in range(0, 7): - l.append((year_min + 5 * shift, year_max + 5 * shift)) - return l - - def cluster_id_to_time_annual_maxima_list(self, altitudes): - cluster_id_to_time_annual_maxima_list = {cluster_id: [] for cluster_id in self.cluster_class.cluster_ids} - for year_min, year_max in self.year_min_and_max_list: - d = self.cluster_id_to_annual_maxima_list(first_year=year_min, - last_year=year_max, - altitudes=altitudes) - for cluster_id, annual_maxima_list in d.items(): - a = np.array(annual_maxima_list) - a = a.flatten() - cluster_id_to_time_annual_maxima_list[cluster_id].append(list(a)) - return cluster_id_to_time_annual_maxima_list - - def cluster_id_to_time_gev_params_list(self, altitudes): - cluster_id_to_gev_params_list = {} - for cluster_id, annual_maxima_list in self.cluster_id_to_time_annual_maxima_list(altitudes=altitudes).items(): - print("cluster", cluster_id) - gev_params_list = [] - for annual_maxima in annual_maxima_list: - if len(annual_maxima) > 20: - print(type(annual_maxima)) - annual_maxima = np.array(annual_maxima) - print(annual_maxima.shape) - annual_maxima = annual_maxima[annual_maxima != 0] - gev_params = fitted_stationary_gev(annual_maxima) - else: - print('here all the time') - gev_params = GevParams(0, 0, 0) - assert gev_params.has_undefined_parameters - gev_params_list.append(gev_params) - cluster_id_to_gev_params_list[cluster_id] = gev_params_list - return cluster_id_to_gev_params_list - - def plot_gev_parameters_against_time(self, altitudes=None): - for param_name in GevParams.PARAM_NAMES[:]: - ax = plt.gca() - for cluster_id, gev_param_list in self.cluster_id_to_time_gev_params_list(altitudes).items(): - print(gev_param_list) - params, years = zip(*[(gev_param.to_dict()[param_name], years) for gev_param, years - in zip(gev_param_list, self.year_min_and_max_list) if - not gev_param.has_undefined_parameters]) - # ax.plot(self.altitudes, params, label=str(cluster_id), linestyle='None', marker='x') - label = self.cluster_class.cluster_id_to_cluster_name[cluster_id] - years = [year[0] for year in years] - ax.plot(years, params, label=label, marker='x') - ax.legend() - ax.set_xlabel('Year') - ax.set_ylabel(param_name) - xlabels = ['-'.join([str(e) for e in t]) for t in self.year_min_and_max_list] - ax.set_xticklabels(xlabels) - plot_name = '{} change with year for altitudes {}'.format(param_name, '+'.join([str(a) for a in altitudes])) - self.show_or_save_to_file(plot_name, no_title=True, tight_layout=True, show=False) - ax.clear() - plt.close() - - -if __name__ == '__main__': - altitudes = [600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600] - # altitudes = [1800, 2100] - pointwise_gev_analysis = PointWIseGevAnalysisForCluster(SafranSnowfall1Day, altitudes=altitudes) - # pointwise_gev_analysis.plot_gev_parameters_against_time(altitudes) - - # pointwise_gev_analysis.plot_gev_parameters_against_altitude() - for altitudes in [[600, 900], - [1200, 1500, 1800], - [2100, 2400, 2700], - [3000, 3300, 3600]][3:]: - print(altitudes) - pointwise_gev_analysis.plot_gev_parameters_against_time(altitudes) diff --git a/projects/contrasting_trends_in_snow_loads/spatial trends/main_constrasting_result.py b/projects/contrasting_trends_in_snow_loads/spatial trends/main_constrasting_result.py deleted file mode 100644 index 57d3cbea1a0522ac5fed6c52922a3abf0c4b53c4..0000000000000000000000000000000000000000 --- a/projects/contrasting_trends_in_snow_loads/spatial trends/main_constrasting_result.py +++ /dev/null @@ -1,96 +0,0 @@ -from multiprocessing.pool import Pool - -import matplotlib as mpl - -mpl.use('Agg') -mpl.rcParams['text.usetex'] = True -mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}'] - -from extreme_trend.visualizers.utils import load_altitude_to_visualizer -from projects.exceeding_snow_loads.section_results.plot_trend_curves import plot_trend_map - -from extreme_data.meteo_france_data.scm_models_data.abstract_study import AbstractStudy -from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranPrecipitation3Days, \ - SafranPrecipitation1Day, SafranPrecipitation5Days, SafranPrecipitation7Days, SafranSnowfall1Day, \ - SafranSnowfall5Days, SafranSnowfall3Days, SafranSnowfall7Days, SafranRainfall1Day, SafranRainfall3Days, \ - SafranRainfall5Days, SafranRainfall7Days - -from extreme_data.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoad3Days, \ - CrocusSnowLoad5Days, CrocusSnowLoad7Days, CrocusSnowLoad1Day -from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import \ - ConfidenceIntervalMethodFromExtremes -from projects.exceeding_snow_loads.section_results.main_result_trends_and_return_levels import \ - compute_minimized_aic - - -def intermediate_result(altitudes, massif_names=None, - model_subsets_for_uncertainty=None, uncertainty_methods=None, - study_class=AbstractStudy, - multiprocessing=False, - save_to_file=True): - """ - Plot all the trends for all altitudes - And enable to plot uncertainty plot for some specific massif_names, uncertainty methods to be fast - :param altitudes: - :param massif_names: - :param non_stationary_uncertainty: - :param uncertainty_methods: - :param study_class: - :return: - """ - # Load altitude to visualizer - altitude_to_visualizer = load_altitude_to_visualizer(altitudes, massif_names, model_subsets_for_uncertainty, - study_class, uncertainty_methods, save_to_file=save_to_file) - # Load variable object efficiently - for v in altitude_to_visualizer.values(): - _ = v.study.year_to_variable_object - # Compute minimized value efficiently - visualizers = list(altitude_to_visualizer.values()) - if multiprocessing: - with Pool(4) as p: - _ = p.map(compute_minimized_aic, visualizers) - else: - for visualizer in visualizers: - _ = compute_minimized_aic(visualizer) - - # Plots - # plot_contrasting_trend_curves(altitude_to_visualizer, all_regions=True) - # plot_contrasting_trend_curves_massif(altitude_to_visualizer, all_regions=True) - # plot_trend_curves(altitude_to_visualizer) - plot_trend_map(altitude_to_visualizer, ) - - -def major_result(): - uncertainty_methods = [ConfidenceIntervalMethodFromExtremes.my_bayes, - ConfidenceIntervalMethodFromExtremes.ci_mle][1:] - massif_names = None - model_subsets_for_uncertainty = None - # altitudes = paper_altitudes - # altitudes = paper_altitudes - altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][:] - # altitudes = [1200, 1500, 1800, 2100, 2400, 2700][:] - snow_load_classes = [CrocusSnowLoad1Day, CrocusSnowLoad3Days, CrocusSnowLoad5Days, CrocusSnowLoad7Days][:] - precipitation_classes = [SafranPrecipitation1Day, SafranPrecipitation3Days, SafranPrecipitation5Days, SafranPrecipitation7Days][:] - snowfall_classes = [SafranSnowfall1Day, SafranSnowfall3Days, SafranSnowfall5Days, SafranSnowfall7Days] - rainfall_classes = [SafranRainfall1Day, SafranRainfall3Days, SafranRainfall5Days, SafranRainfall7Days] - study_classes = precipitation_classes + snow_load_classes - study_classes = snowfall_classes + rainfall_classes - for study_class in precipitation_classes[:1]: - intermediate_result(altitudes, massif_names, model_subsets_for_uncertainty, - uncertainty_methods, study_class, multiprocessing=True) - - -""" - -est ce qu il y a une croissance signifcative en pluie, - -est ce qu'il y a une decroissance signifcatieve à partir d'une certaine altitude -""" - -if __name__ == '__main__': - major_result() - # intermediate_result(altitudes=[1500, 1800][:], massif_names=None, - # uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes, - # ConfidenceIntervalMethodFromExtremes.ci_mle][1:], - # multiprocessing=True, - # save_to_file=False) diff --git a/projects/contrasting_trends_in_snow_loads/spatial trends/main_spatial_relative_change_in_maxima_at_fixed_altitude.py b/projects/contrasting_trends_in_snow_loads/spatial trends/main_spatial_relative_change_in_maxima_at_fixed_altitude.py deleted file mode 100644 index 99425900b880c081d2bf7b755d690ad5d98fdbea..0000000000000000000000000000000000000000 --- a/projects/contrasting_trends_in_snow_loads/spatial trends/main_spatial_relative_change_in_maxima_at_fixed_altitude.py +++ /dev/null @@ -1,70 +0,0 @@ -import pandas as pd -from extreme_data.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoad3Days, \ - CrocusSnowLoadTotal -from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall, SafranRainfall, \ - SafranPrecipitation, SafranPrecipitation1Day, SafranSnowfall1Day, SafranTemperature, \ - SafranNormalizedPreciptationRateOnWetDays, SafranNormalizedPreciptationRate -from extreme_data.meteo_france_data.scm_models_data.safran.safran_variable import \ - SafranNormalizedPrecipitationRateOnWetDaysVariable -from extreme_data.meteo_france_data.scm_models_data.utils import Season -from extreme_data.meteo_france_data.scm_models_data.visualization.study_visualizer import \ - StudyVisualizer -import matplotlib.pyplot as plt - - -def relative_change_in_maxima_wrt_altitude(): - save_to_file = True - study_class = [SafranSnowfall1Day, SafranPrecipitation1Day, SafranPrecipitation, SafranRainfall, - SafranSnowfall, CrocusSnowLoad3Days, CrocusSnowLoadTotal, - SafranTemperature, SafranNormalizedPreciptationRateOnWetDays, - SafranNormalizedPreciptationRate][1] - altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][::-1] - relative = True - - for altitude in altitudes: - - ax = plt.gca() - study = study_class(altitude=altitude, season=Season.winter_extended) - # study = study_class(altitude=altitude, nb_consecutive_days=3) - massif_name_to_value = {} - for massif_name in study.study_massif_names: - if study_class is SafranTemperature: - s = study.observations_annual_mean.loc[massif_name] - else: - s = study.observations_annual_maxima.df_maxima_gev.loc[massif_name] - year_limit = 1989 - df_before, df_after = s.loc[:year_limit], s.loc[year_limit + 1:] - df_before, df_after = df_before.mean(), df_after.mean() - # df_before, df_after = df_before.median(), df_after.median() - if relative: - change_value = 100 * (df_after - df_before) / df_before - else: - change_value = (df_after - df_before) - massif_name_to_value[massif_name] = change_value - print(massif_name_to_value) - # Plot - # massif_name_to_value = {m: i for i, m in enumerate(study.study_massif_names)} - max_values = max([abs(e) for e in massif_name_to_value.values()]) * 1.05 - print(max_values) - variable_name = study.variable_name - prefix = 'Relative' if relative else '' - str_season = str(study.season).split('.')[-1] - study.visualize_study(ax=ax, massif_name_to_value=massif_name_to_value, - vmin=-max_values, vmax=max_values, - add_colorbar=True, - replace_blue_by_white=False, - show=False, - label='{} changes in mean {} maxima\n' - 'of {}\n between 1959-1989 and 1990-2019\n at {}m (%)\n' - ''.format(prefix, str_season, variable_name, study.altitude) - ) - study_visualizer = StudyVisualizer(study, save_to_file=save_to_file) - study_visualizer.plot_name = '{}_changes_in_maxima'.format(prefix) - study_visualizer.show_or_save_to_file() - ax.clear() - plt.close() - - -if __name__ == '__main__': - relative_change_in_maxima_wrt_altitude() - # test() diff --git a/projects/contrasting_trends_in_snow_loads/spatial trends/plot_contrasting_trend_curves.py b/projects/contrasting_trends_in_snow_loads/spatial trends/plot_contrasting_trend_curves.py deleted file mode 100644 index 1a75bb42eff8c49c742feb11bb2b3873b8c0c14e..0000000000000000000000000000000000000000 --- a/projects/contrasting_trends_in_snow_loads/spatial trends/plot_contrasting_trend_curves.py +++ /dev/null @@ -1,144 +0,0 @@ -from typing import Dict -import matplotlib.pyplot as plt - -from extreme_data.meteo_france_data.scm_models_data.abstract_extended_study import AbstractExtendedStudy -from extreme_data.meteo_france_data.scm_models_data.visualization.utils import create_adjusted_axes -from projects.exceeding_snow_loads.utils import dpi_paper1_figure -from extreme_trend.visualizers.study_visualizer_for_non_stationary_trends import \ - StudyVisualizerForNonStationaryTrends - - -def plot_contrasting_trend_curves_massif(altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends], - all_regions=False): - """ - Plot a single trend curves - :return: - """ - visualizers = list(altitude_to_visualizer.values()) - visualizer = visualizers[0] - altitudes = list(altitude_to_visualizer.keys()) - - ax = create_adjusted_axes(1, 1) - # ax_twinx = ax.twinx() - ax_twinx = ax - ax_twiny = ax.twiny() - - # parameters - width = 150 - size = 20 - legend_fontsize = 20 - color = 'white' - labelsize = 15 - linewidth = 3 - - for ax_horizontal in [ax, ax_twiny]: - if ax_horizontal == ax_twiny: - ax_horizontal.plot(altitudes, [0 for _ in altitudes], linewidth=0) - else: - ax_horizontal.set_xlabel('Altitude', fontsize=legend_fontsize) - ax_horizontal.set_xticks(altitudes) - # ax_horizontal.set_xlim([700, 5000]) - ax_horizontal.tick_params(labelsize=labelsize) - - # Set the number of massifs on the upper axis - ax_twiny.set_xticklabels([v.study.nb_study_massif_names for v in altitude_to_visualizer.values()]) - ax_twiny.set_xlabel('Total number of massifs at each altitude (for the percentage)', fontsize=legend_fontsize) - - ax_twinx.yaxis.grid() - - ax_twinx.set_ylabel(visualizer.label, fontsize=legend_fontsize) - for j, massif_name in enumerate(visualizer.study.study_massif_names): - massif_visualizers = [v for v in visualizers if massif_name in v.massif_name_to_change_value] - changes = [v.massif_name_to_relative_change_value[massif_name] for v in massif_visualizers] - massif_altitudes = [v.study.altitude for v in massif_visualizers] - label = massif_name.replace('-', '').replace('_', '') - if j < 10: - linestyle = 'solid' - elif j < 20: - linestyle = 'dashed' - else: - linestyle = 'dotted' - ax_twinx.plot(massif_altitudes, changes, label=label, linewidth=linewidth, marker='o', linestyle=linestyle) - ax_twinx.legend(loc='upper right', prop={'size': 5}) - - ax.axhline(y=0, color='k') - - # Save plot - visualizer.plot_name = 'Trend curves for' + visualizer.study.variable_name - visualizer.show_or_save_to_file(no_title=True, dpi=dpi_paper1_figure, folder_for_variable=False) - plt.close() - - -def plot_contrasting_trend_curves(altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends], - all_regions=False): - """ - Plot a single trend curves - :return: - """ - visualizer = list(altitude_to_visualizer.values())[0] - - ax = create_adjusted_axes(1, 1) - # ax_twinx = ax.twinx() - ax_twinx = ax - ax_twiny = ax.twiny() - - trend_summary_values = list(zip(*[v.trend_summary_contrasting_values(regions=all_regions) for v in altitude_to_visualizer.values()])) - altitudes, *mean_changes = trend_summary_values - - # parameters - width = 150 - size = 20 - legend_fontsize = 20 - color = 'white' - labelsize = 15 - linewidth = 3 - # ax.bar(altitudes, percent_decrease, width=width, color=color, edgecolor='blue', label='decreasing trend', - # linewidth=linewidth) - # ax.bar(altitudes, percent_decrease_signi, width=width, color=color, edgecolor='black', - # label='significative decreasing trend', - # linewidth=linewidth) - # ax.legend(loc='upper left', prop={'size': size}) - - for ax_horizontal in [ax, ax_twiny]: - if ax_horizontal == ax_twiny: - ax_horizontal.plot(altitudes, [0 for _ in altitudes], linewidth=0) - else: - ax_horizontal.set_xlabel('Altitude', fontsize=legend_fontsize) - ax_horizontal.set_xticks(altitudes) - # ax_horizontal.set_xlim([700, 5000]) - ax_horizontal.tick_params(labelsize=labelsize) - - # Set the number of massifs on the upper axis - ax_twiny.set_xticklabels([v.study.nb_study_massif_names for v in altitude_to_visualizer.values()]) - ax_twiny.set_xlabel('Total number of massifs at each altitude (for the percentage)', fontsize=legend_fontsize) - - # ax.set_ylabel('Massifs with decreasing trend (\%)', fontsize=legend_fontsize) - # max_percent = int(max(percent_decrease)) - # n = 2 + (max_percent // 10) - # ax_ticks = [10 * i for i in range(n)] - # # upper_lim = max_percent + 3 - # upper_lim = n + 5 - # ax_lim = [0, upper_lim] - # for axis in [ax, ax_twinx]: - # axis.set_ylim(ax_lim) - # axis.set_yticks(ax_ticks) - # axis.tick_params(labelsize=labelsize) - ax_twinx.yaxis.grid() - - ax_twinx.set_ylabel(visualizer.label, fontsize=legend_fontsize) - for j, (region_name, mean_change) in enumerate(zip(AbstractExtendedStudy.region_names, mean_changes)): - if len(mean_changes) > 2: - label = region_name - elif len(mean_changes) == 2: - label = 'North' if j == 0 else 'South' - else: - label = 'Mean relative change' - ax_twinx.plot(altitudes, mean_change, label=label, linewidth=linewidth, marker='o') - ax_twinx.legend(loc='upper right', prop={'size': size}) - - ax.axhline(y=0, color='k') - - # Save plot - visualizer.plot_name = 'Trend curves for' + visualizer.study.variable_name - visualizer.show_or_save_to_file(no_title=True, dpi=dpi_paper1_figure, folder_for_variable=False) - plt.close() diff --git a/projects/contrasting_trends_in_snow_loads/weather_types_analysis/main_distribution_wps.py b/projects/contrasting_trends_in_snow_loads/weather_types_analysis/main_distribution_wps.py deleted file mode 100644 index 73863db9bdba02a01455d0bde4f88227466bd4b6..0000000000000000000000000000000000000000 --- a/projects/contrasting_trends_in_snow_loads/weather_types_analysis/main_distribution_wps.py +++ /dev/null @@ -1,111 +0,0 @@ -import pandas as pd -import numpy as np - -from extreme_data.meteo_france_data.scm_models_data.abstract_extended_study import AbstractExtendedStudy -from extreme_data.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoad3Days, CrocusSnowLoad1Day -from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranPrecipitation1Day, \ - SafranPrecipitation3Days -from extreme_data.meteo_france_data.scm_models_data.utils import Season - - -def main_spatial_distribution_wps(study_class, year_min=1959, year_max=2008, limit_for_the_percentage=None): - study = study_class(altitude=1800, year_min=year_min, year_max=year_max) - for region_name in AbstractExtendedStudy.region_names: - massif_names = AbstractExtendedStudy.region_name_to_massif_names[region_name] - print('\n \n', region_name, '\n') - for nb_top in [study.nb_years, 10, 1][1:2]: - print(study.df_for_top_annual_maxima(nb_top=nb_top, massif_names=massif_names, limit_for_the_percentage=limit_for_the_percentage), '\n') - - -""" - Northern Alps - % count mean std min median max -Steady Oceanic 92 65 159 35 121 146 282 - - Central Alps - % count mean std min median max -Steady Oceanic 72 51 136 31 97 130 235 -South Circulation 15 11 128 18 105 126 161 - - Southern Alps - % count mean min median max -South Circulation 61 37 147 79 147 235 -Steady Oceanic 13 8 121 93 114 187 - - Extreme South Alps - % count mean min median max -South Circulation 80 24 174 101 166 306 - - -Process finished with exit code 0 - - - - - - -""" - - -def main_temporal_distribution_wps(study_class, year_min=1959, year_max=2008, limit_for_the_percentage=None): - altitude = 1800 - intermediate_year = year_min + round(float(year_max - year_min) / 2) - study_before = study_class(altitude=altitude, year_min=year_min, year_max=intermediate_year) - study_after = study_class(altitude=altitude, year_min=intermediate_year+1, year_max=year_max) - for region_name in AbstractExtendedStudy.region_names: - massif_names = AbstractExtendedStudy.region_name_to_massif_names[region_name] - print('\n \n', '{} ({} massifs)'.format(region_name, len(massif_names)), '\n') - for nb_top in [study_before.nb_years, 10, 1, 10][3:]: - print(study_before.df_for_top_annual_maxima(nb_top=nb_top, massif_names=massif_names, limit_for_the_percentage=limit_for_the_percentage), '\n') - print(study_after.df_for_top_annual_maxima(nb_top=nb_top, massif_names=massif_names, limit_for_the_percentage=limit_for_the_percentage), '\n') - -""" - Northern Alps (7 massifs) - % count mean min median max -Top 10 maxima (1959 -1983) -Steady Oceanic 92 65 131 106 128 202 - -Top 10 maxima (1984 -2008) -Steady Oceanic 87 61 149 103 133 282 - - Central Alps (7 massifs) - % count mean min median max -Top 10 maxima (1959 -1983) -Steady Oceanic 68 48 111 76 109 190 - -Top 10 maxima (1984 -2008) -Steady Oceanic 77 54 123 81 107 235 -South Circulation 14 10 122 93 120 161 - - Southern Alps (6 massifs) - % count mean std min median max -Top 10 maxima (1959 -1983) -South Circulation 50 30 114 32 70 112 197 -Southwest Circulation 13 8 107 24 79 97 140 -Steady Oceanic 13 8 97 14 82 93 122 - -Top 10 maxima (1984 -2008) -South Circulation 58 35 135 73 125 235 -Steady Oceanic 23 14 107 69 104 187 - - Extreme South Alps (3 massifs) - % count mean min median max -Top 10 maxima (1959 -1983) -South Circulation 66 20 141 84 146 212 -Southwest Circulation 13 4 122 95 123 146 - -Top 10 maxima (1984 -2008) -South Circulation 76 23 158 74 138 306 - - - - - -""" - - -if __name__ == '__main__': - limit_percentage = 1 - study_class = [CrocusSnowLoad1Day, SafranPrecipitation1Day, SafranPrecipitation3Days][-1] - main_spatial_distribution_wps(study_class, limit_for_the_percentage=limit_percentage) - # main_temporal_distribution_wps(study_class, limit_for_the_percentage=limit_percentage) diff --git a/projects/seasonal_analysis/__init__.py b/projects/seasonal_analysis/__init__.py deleted file mode 100644 index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000 diff --git a/projects/thesis_report/__init__.py b/projects/thesis_report/__init__.py deleted file mode 100644 index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000