From 00a0bd4fae72e97c65873d9e4a27545d57892a2c Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Tue, 18 Feb 2020 11:47:58 +0100
Subject: [PATCH] [paper 2] add snow load in 1 day. add plot for contrasting
trend curve for each climatic regions
---
.../scm_models_data/crocus/crocus.py | 6 +-
.../crocus/crocus_variables.py | 10 +++
papers/contrasting_snow_loads/main_result.py | 77 ++++++++++++++++++
...tive_change_in_maxima_at_fixed_altitude.py | 4 +-
.../plot_contrasting_trend_curves.py | 78 +++++++++++++++++++
...dy_visualizer_for_non_stationary_trends.py | 16 ++++
6 files changed, 188 insertions(+), 3 deletions(-)
create mode 100644 papers/contrasting_snow_loads/main_result.py
create mode 100644 papers/contrasting_snow_loads/plot_contrasting_trend_curves.py
diff --git a/experiment/meteo_france_data/scm_models_data/crocus/crocus.py b/experiment/meteo_france_data/scm_models_data/crocus/crocus.py
index 35acef69..aafab19e 100644
--- a/experiment/meteo_france_data/scm_models_data/crocus/crocus.py
+++ b/experiment/meteo_france_data/scm_models_data/crocus/crocus.py
@@ -7,7 +7,7 @@ from experiment.meteo_france_data.scm_models_data.abstract_study import Abstract
from experiment.meteo_france_data.scm_models_data.crocus.crocus_variables import CrocusTotalSweVariable, \
CrocusDepthVariable, CrocusRecentSweVariableThreeDays, TotalSnowLoadVariable, RecentSnowLoadVariableThreeDays, \
CrocusSnowLoadEurocodeVariable, CrocusDensityVariable, RecentSnowLoadVariableFiveDays, \
- RecentSnowLoadVariableSevenDays
+ RecentSnowLoadVariableSevenDays, RecentSnowLoadVariableOneDay
class Crocus(AbstractStudy):
@@ -20,6 +20,7 @@ class Crocus(AbstractStudy):
RecentSnowLoadVariableThreeDays, TotalSnowLoadVariable,
CrocusSnowLoadEurocodeVariable,
CrocusDensityVariable,
+ RecentSnowLoadVariableOneDay,
RecentSnowLoadVariableFiveDays,
RecentSnowLoadVariableSevenDays]
super().__init__(variable_class, *args, **kwargs)
@@ -60,6 +61,9 @@ class CrocusSnowLoadTotal(Crocus):
def __init__(self, *args, **kwargs):
Crocus.__init__(self, TotalSnowLoadVariable, *args, **kwargs)
+class CrocusSnowLoad1Day(CrocusSweTotal):
+ def __init__(self, *args, **kwargs):
+ Crocus.__init__(self, RecentSnowLoadVariableOneDay, *args, **kwargs)
class CrocusSnowLoad3Days(CrocusSweTotal):
def __init__(self, *args, **kwargs):
diff --git a/experiment/meteo_france_data/scm_models_data/crocus/crocus_variables.py b/experiment/meteo_france_data/scm_models_data/crocus/crocus_variables.py
index 02c9d3b1..ec226fae 100644
--- a/experiment/meteo_france_data/scm_models_data/crocus/crocus_variables.py
+++ b/experiment/meteo_france_data/scm_models_data/crocus/crocus_variables.py
@@ -20,6 +20,14 @@ class CrocusTotalSweVariable(CrocusVariable):
return 'WSN_T_ISBA'
+class CrocusRecentSweVariableOneDay(CrocusTotalSweVariable):
+ NAME = 'Snow Water Equivalent last 1 day'
+
+ @classmethod
+ def keyword(cls):
+ return 'SWE_1DY_ISBA'
+
+
class CrocusRecentSweVariableThreeDays(CrocusTotalSweVariable):
NAME = 'Snow Water Equivalent last 3 days'
@@ -52,6 +60,8 @@ class AbstractSnowLoadVariable(CrocusVariable):
snow_pressure = self.snow_load_multiplication_factor * super().daily_time_serie_array
return snow_pressure
+class RecentSnowLoadVariableOneDay(AbstractSnowLoadVariable, CrocusRecentSweVariableOneDay):
+ NAME = 'Snow load last 1 day'
class RecentSnowLoadVariableThreeDays(AbstractSnowLoadVariable, CrocusRecentSweVariableThreeDays):
NAME = 'Snow load last 3 days'
diff --git a/papers/contrasting_snow_loads/main_result.py b/papers/contrasting_snow_loads/main_result.py
new file mode 100644
index 00000000..dba2c438
--- /dev/null
+++ b/papers/contrasting_snow_loads/main_result.py
@@ -0,0 +1,77 @@
+from multiprocessing.pool import Pool
+
+import matplotlib as mpl
+
+from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal, CrocusSnowLoad3Days, \
+ CrocusSnowLoad5Days, CrocusSnowLoad7Days, CrocusSnowLoad1Day
+from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import \
+ ConfidenceIntervalMethodFromExtremes
+from papers.contrasting_snow_loads.plot_contrasting_trend_curves import plot_contrasting_trend_curves
+from papers.exceeding_snow_loads.paper_main_utils import load_altitude_to_visualizer
+from papers.exceeding_snow_loads.paper_utils import paper_study_classes, paper_altitudes
+from papers.exceeding_snow_loads.result_trends_and_return_levels.main_result_trends_and_return_levels import \
+ compute_minimized_aic
+from papers.exceeding_snow_loads.result_trends_and_return_levels.plot_selection_curves import plot_selection_curves
+from papers.exceeding_snow_loads.result_trends_and_return_levels.plot_trend_curves import plot_trend_curves, \
+ plot_trend_map
+from papers.exceeding_snow_loads.result_trends_and_return_levels.plot_uncertainty_curves import plot_uncertainty_massifs
+from papers.exceeding_snow_loads.result_trends_and_return_levels.plot_uncertainty_histogram import \
+ plot_uncertainty_histogram
+from root_utils import NB_CORES
+
+mpl.rcParams['text.usetex'] = True
+mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
+
+
+def intermediate_result(altitudes, massif_names=None,
+ model_subsets_for_uncertainty=None, uncertainty_methods=None,
+ study_class=CrocusSnowLoad3Days,
+ multiprocessing=False):
+ """
+ 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)
+ # 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(NB_CORES) 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)
+
+
+def major_result():
+ uncertainty_methods = [ConfidenceIntervalMethodFromExtremes.my_bayes,
+ ConfidenceIntervalMethodFromExtremes.ci_mle][1:]
+ massif_names = None
+ model_subsets_for_uncertainty = None
+ altitudes = paper_altitudes
+ altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700]
+ study_classes = [CrocusSnowLoad1Day, CrocusSnowLoad3Days, CrocusSnowLoad5Days, CrocusSnowLoad7Days][:]
+ for study_class in study_classes:
+ intermediate_result(altitudes, massif_names, model_subsets_for_uncertainty,
+ uncertainty_methods, study_class)
+
+
+if __name__ == '__main__':
+ major_result()
+ # intermediate_result(altitudes=[1500, 1800][:], massif_names=None,
+ # uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
+ # ConfidenceIntervalMethodFromExtremes.ci_mle][1:],
+ # multiprocessing=True)
\ No newline at end of file
diff --git a/papers/contrasting_snow_loads/main_spatial_relative_change_in_maxima_at_fixed_altitude.py b/papers/contrasting_snow_loads/main_spatial_relative_change_in_maxima_at_fixed_altitude.py
index 4451ab32..fe763ba0 100644
--- a/papers/contrasting_snow_loads/main_spatial_relative_change_in_maxima_at_fixed_altitude.py
+++ b/papers/contrasting_snow_loads/main_spatial_relative_change_in_maxima_at_fixed_altitude.py
@@ -9,11 +9,11 @@ from experiment.meteo_france_data.scm_models_data.visualization.study_visualizat
StudyVisualizer
import matplotlib.pyplot as plt
-from experiment.exceeding_snow_loads.discussion_data_comparison_with_eurocode.crocus_study_comparison_with_eurocode import \
+from papers.exceeding_snow_loads.discussion_data_comparison_with_eurocode.crocus_study_comparison_with_eurocode import \
CrocusDifferenceSnowLoad, \
CrocusSnowDensityAtMaxofSwe, CrocusDifferenceSnowLoadRescaledAndEurocodeToSeeSynchronization, \
CrocusSnowDepthAtMaxofSwe, CrocusSnowDepthDifference
-from experiment.exceeding_snow_loads.paper_utils import dpi_paper1_figure
+from papers.exceeding_snow_loads.paper_utils import dpi_paper1_figure
def test():
diff --git a/papers/contrasting_snow_loads/plot_contrasting_trend_curves.py b/papers/contrasting_snow_loads/plot_contrasting_trend_curves.py
new file mode 100644
index 00000000..616bb728
--- /dev/null
+++ b/papers/contrasting_snow_loads/plot_contrasting_trend_curves.py
@@ -0,0 +1,78 @@
+from typing import Dict
+import matplotlib.pyplot as plt
+
+from experiment.meteo_france_data.scm_models_data.abstract_extended_study import AbstractExtendedStudy
+from experiment.meteo_france_data.scm_models_data.visualization.utils import create_adjusted_axes
+from papers.exceeding_snow_loads.paper_utils import dpi_paper1_figure
+from papers.exceeding_snow_loads.study_visualizer_for_non_stationary_trends import \
+ StudyVisualizerForNonStationaryTrends
+
+
+def plot_contrasting_trend_curves(altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends]):
+ """
+ Plot a single trend curves
+ :return:
+ """
+ visualizer = list(altitude_to_visualizer.values())[0]
+
+ ax = create_adjusted_axes(1, 1)
+ ax_twinx = ax.twinx()
+ ax_twiny = ax.twiny()
+
+ trend_summary_values = list(zip(*[v.trend_summary_contrasting_values() 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.yaxis.grid()
+
+ label_curve = (visualizer.label).replace('change', 'decrease')
+ ax_twinx.set_ylabel(label_curve.replace('', ''), fontsize=legend_fontsize)
+ for region_name, mean_change in zip(AbstractExtendedStudy.region_names, mean_changes):
+ if len(mean_changes) > 1:
+ label = region_name
+ 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})
+
+ # Save plot
+ visualizer.plot_name = 'Trend curves'
+ visualizer.show_or_save_to_file(no_title=True, dpi=dpi_paper1_figure)
+ plt.close()
diff --git a/papers/exceeding_snow_loads/study_visualizer_for_non_stationary_trends.py b/papers/exceeding_snow_loads/study_visualizer_for_non_stationary_trends.py
index ef742b53..a7f3e056 100644
--- a/papers/exceeding_snow_loads/study_visualizer_for_non_stationary_trends.py
+++ b/papers/exceeding_snow_loads/study_visualizer_for_non_stationary_trends.py
@@ -429,3 +429,19 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
# if massif_name_to_region_name[m] == region_name]
# mean_decreases.append(compute_mean_decrease(change_values))
return (self.altitude, percentage_decrease, percentage_decrease_significative, *mean_decreases)
+
+ def trend_summary_contrasting_values(self):
+ # trend_tests = list(self.massif_name_to_trend_test_that_minimized_aic.values())
+ # decreasing_trend_tests = [t for t in trend_tests if t.time_derivative_of_return_level < 0]
+ # percentage_decrease = 100 * len(decreasing_trend_tests) / len(trend_tests)
+ # significative_decrease_trend_tests = [t for t in decreasing_trend_tests if t.is_significant]
+ # percentage_decrease_significative = 100 * len(significative_decrease_trend_tests) / len(trend_tests)
+ compute_mean_change = lambda l: np.mean(np.array(list(l)))
+ mean_changes = [compute_mean_change(self.massif_name_to_relative_change_value.values())]
+ # Compute mean relatives per regions (for the moment i don't add the region means)
+ massif_name_to_region_name = AbstractExtendedStudy.massif_name_to_region_name
+ for region_name in AbstractExtendedStudy.real_region_names:
+ change_values = [v for m, v in self.massif_name_to_relative_change_value.items()
+ if massif_name_to_region_name[m] == region_name]
+ mean_changes.append(compute_mean_change(change_values))
+ return (self.altitude, *mean_changes)
--
GitLab