diff --git a/experiment/meteo_france_data/scm_models_data/abstract_study.py b/experiment/meteo_france_data/scm_models_data/abstract_study.py
index 9506fb209f6e9dad3bb0e80219430b7068aa7b05..abe1fbef39ba6c3d943374b6b7db22b18da4169b 100644
--- a/experiment/meteo_france_data/scm_models_data/abstract_study.py
+++ b/experiment/meteo_france_data/scm_models_data/abstract_study.py
@@ -39,6 +39,7 @@ f = io.StringIO()
with redirect_stdout(f):
from simpledbf import Dbf5
+filled_marker_legend_list = ['Filled marker =', 'Selected model is significant', 'w.r.t $\mathcal{M}_0$']
class AbstractStudy(object):
"""
@@ -117,6 +118,33 @@ class AbstractStudy(object):
year_to_annual_maxima[year] = time_serie.max(axis=0)
return year_to_annual_maxima
+ @cached_property
+ def year_to_winter_annual_maxima(self) -> OrderedDict:
+ # Map each year to an array of size nb_massif
+ year_to_annual_maxima = OrderedDict()
+ for year, time_serie in self._year_to_max_daily_time_serie.items():
+ year_to_annual_maxima[year] = time_serie[91:-61].max(axis=0)
+ return year_to_annual_maxima
+
+ @property
+ def observations_winter_annual_maxima(self) -> AnnualMaxima:
+ return AnnualMaxima(df_maxima_gev=pd.DataFrame(self.year_to_winter_annual_maxima, index=self.study_massif_names))
+
+ @cached_property
+ def year_to_summer_annual_maxima(self) -> OrderedDict:
+ # Map each year to an array of size nb_massif
+ year_to_annual_maxima = OrderedDict()
+ for year, time_serie in self._year_to_max_daily_time_serie.items():
+ annual_maxima = np.concatenate((time_serie[:91], time_serie[-61:]), axis=0).max(axis=0)
+ print(annual_maxima)
+ year_to_annual_maxima[year] = annual_maxima
+ return year_to_annual_maxima
+
+
+ @property
+ def observations_summer_annual_maxima(self) -> AnnualMaxima:
+ return AnnualMaxima(df_maxima_gev=pd.DataFrame(self.year_to_summer_annual_maxima, index=self.study_massif_names))
+
@cached_property
def year_to_annual_maxima_index(self) -> OrderedDict:
# Map each year to an array of size nb_massif
@@ -298,7 +326,7 @@ class AbstractStudy(object):
@classmethod
def visualize_study(cls, ax=None, massif_name_to_value: Union[None, Dict[str, float]] = None,
show=True, fill=True,
- replace_blue_by_white=True,
+ replace_blue_by_white=False,
label=None, add_text=False, cmap=None, add_colorbar=False, vmax=100, vmin=0,
default_color_for_missing_massif='gainsboro',
default_color_for_nan_values='w',
@@ -406,7 +434,7 @@ class AbstractStudy(object):
legend_elements = cls.get_legend_for_model_symbol(marker_style_to_label_name, markersize=8)
ax.legend(handles=legend_elements[:], loc='upper right')
# ax.legend(handles=legend_elements[4:], bbox_to_anchor=(0.01, 0.03), loc='lower left')
- ax.annotate("Filled symbol = significant trend w.r.t $\mathcal{M}_0$",
+ ax.annotate(' '.join(filled_marker_legend_list),
xy=(0.05, 0.015), xycoords='axes fraction', fontsize=7)
if show:
diff --git a/experiment/paper 2 - contrasting/__init__.py b/experiment/paper 2 - contrasting/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/experiment/paper 2 - contrasting/main_spatial_relative_change_in_maxima_at_fixed_altitude.py b/experiment/paper 2 - contrasting/main_spatial_relative_change_in_maxima_at_fixed_altitude.py
new file mode 100644
index 0000000000000000000000000000000000000000..b2c1c40640f2d9a7c04cd88eb5b6c0232d4194be
--- /dev/null
+++ b/experiment/paper 2 - contrasting/main_spatial_relative_change_in_maxima_at_fixed_altitude.py
@@ -0,0 +1,65 @@
+from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusDepth, CrocusSnowLoad3Days, \
+ CrocusSnowLoadTotal
+from experiment.meteo_france_data.scm_models_data.crocus.crocus_variables import CrocusDepthVariable
+from experiment.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall, SafranRainfall, SafranTotalPrecip
+from experiment.meteo_france_data.scm_models_data.safran.safran_variable import SafranTotalPrecipVariable
+from experiment.meteo_france_data.scm_models_data.visualization.study_visualization.main_study_visualizer import \
+ study_iterator_global, SCM_STUDY_CLASS_TO_ABBREVIATION, snow_density_str, ALL_ALTITUDES_WITHOUT_NAN
+from experiment.meteo_france_data.scm_models_data.visualization.study_visualization.study_visualizer import \
+ StudyVisualizer
+import matplotlib.pyplot as plt
+
+from experiment.paper_past_snow_loads.discussion_data_comparison_with_eurocode.crocus_study_comparison_with_eurocode import \
+ CrocusDifferenceSnowLoad, \
+ CrocusSnowDensityAtMaxofSwe, CrocusDifferenceSnowLoadRescaledAndEurocodeToSeeSynchronization, \
+ CrocusSnowDepthAtMaxofSwe, CrocusSnowDepthDifference
+from experiment.paper_past_snow_loads.paper_utils import dpi_paper1_figure
+
+
+def density_wrt_altitude():
+ """
+ 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 = [SafranTotalPrecip, SafranRainfall, SafranSnowfall, CrocusSnowLoad3Days, CrocusSnowLoadTotal][0]
+ # altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][::-1]
+ altitudes = ALL_ALTITUDES_WITHOUT_NAN[1:]
+ altitudes = [900]
+
+
+
+ for altitude in altitudes:
+
+ ax = plt.gca()
+ study = study_class(altitude=altitude, nb_consecutive_days=3)
+ massif_name_to_value = {}
+ for massif_name in study.study_massif_names:
+ s = study.observations_summer_annual_maxima.df_maxima_gev.loc[massif_name]
+ year_limit = 1987
+ print(s)
+ df_before, df_after = s.loc[:year_limit], s.loc[year_limit+1:]
+ print(df_before, df_after)
+
+ df_before, df_after = df_before.median(), df_after.median()
+ relative_change_value = 100 * (df_after - df_before) / df_before
+ massif_name_to_value[massif_name] = relative_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()]) + 5
+ print(max_values)
+ 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,
+ label='Relative changes for \n{}\n at {}m (%)'.format(study.variable_name, study.altitude)
+ )
+ plt.show()
+ ax.clear()
+
+
+if __name__ == '__main__':
+ density_wrt_altitude()
diff --git a/experiment/paper_past_snow_loads/data/main_example_swe_total_plot.py b/experiment/paper_past_snow_loads/data/main_example_swe_total_plot.py
index d912e4c8b765d6ed549681d27d022da2e89008bb..eb32d1cca9f91f2395040e492610c06263d17945 100644
--- a/experiment/paper_past_snow_loads/data/main_example_swe_total_plot.py
+++ b/experiment/paper_past_snow_loads/data/main_example_swe_total_plot.py
@@ -8,21 +8,11 @@ from experiment.meteo_france_data.scm_models_data.visualization.study_visualizat
from experiment.paper_past_snow_loads.paper_utils import dpi_paper1_figure
-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 = True
- study_class = CrocusSnowLoadTotal
- examples_for_the_paper = True
- ax = plt.gca()
+def tuples_for_examples_paper1(examples_for_the_paper=True):
if examples_for_the_paper:
- ax.set_ylim([0, 20])
- ax.set_yticks(list(range(0, 21, 2)))
+
marker_altitude_massif_name_for_paper1 = [
('magenta', 900, 'Ubaye'),
('darkmagenta', 1800, 'Vercors'),
@@ -33,6 +23,26 @@ def max_graph_annual_maxima_poster():
('yellow', 600, 'Ubaye'),
('purple', 600, 'Parpaillon'),
]
+ return marker_altitude_massif_name_for_paper1
+
+tuples_for_examples_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 = True
+ study_class = CrocusSnowLoadTotal
+
+ examples_for_the_paper = True
+
+ ax = plt.gca()
+ ax.set_ylim([0, 20])
+ ax.set_yticks(list(range(0, 21, 2)))
+
+ marker_altitude_massif_name_for_paper1 = tuples_for_examples_paper1(examples_for_the_paper)
for color, altitude, massif_name in marker_altitude_massif_name_for_paper1[::-1]:
for study in study_iterator_global([study_class], altitudes=[altitude]):
diff --git a/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py b/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py
index 1948b79daf61f9e1476886fdaac900649e95e02e..622a282b80e3203419f2225bb4c16b89b58f4978 100644
--- a/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py
+++ b/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py
@@ -4,7 +4,7 @@ import matplotlib.pyplot as plt
import numpy as np
from experiment.eurocode_data.utils import EUROCODE_RETURN_LEVEL_STR, EUROCODE_ALTITUDES
-from experiment.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
+from experiment.meteo_france_data.scm_models_data.abstract_study import AbstractStudy, filled_marker_legend_list
from experiment.meteo_france_data.scm_models_data.visualization.study_visualization.main_study_visualizer import \
SCM_STUDY_CLASS_TO_ABBREVIATION
from experiment.paper_past_snow_loads.paper_utils import dpi_paper1_figure, ModelSubsetForUncertainty
@@ -41,12 +41,8 @@ def plot_subgroup_uncertainty_massifs(altitude_to_visualizer: Dict[int, StudyVis
"""Create a plot with a maximum of 4 massif names
We will save the plot at this level
"""
- visualizer = list(altitude_to_visualizer.values())[0]
nb_massif_names = len(massif_names)
assert nb_massif_names <= 5
- # axes = create_adjusted_axes(nb_massif_names, visualizer.nb_contexts)
- # if nb_massif_names == 1:
- # axes = [axes]
for massif_name in massif_names:
plot_single_uncertainty_massif(altitude_to_visualizer, massif_name)
@@ -120,8 +116,15 @@ def plot_single_uncertainty_massif_and_non_stationary_context(ax, massif_name, m
ax.legend(loc=2, prop={'size': legend_size})
# ax.set_ylim([-1, 16])
ax.set_xlim([200, 1900])
- if massif_name == 'Maurienne':
- ax.set_ylim([-1.5, 13])
+ if massif_name in ['Maurienne', 'Chartreuse', 'Beaufortain']:
+ ax.set_ylim([-1.5, 13.5])
+ ax.set_yticks([2 * i for i in range(7)])
+ if massif_name in ['Vercors']:
+ ax.set_ylim([-1, 10])
+ ax.set_yticks([2 * i for i in range(6)])
+
+
+
# add_title(ax, eurocode_region, massif_name, non_stationary_context)
ax.set_xticks(altitudes)
ax.tick_params(labelsize=fontsize_label)
@@ -173,8 +176,9 @@ def plot_tdrl_bars(altitude_to_visualizer, ax, massif_name, valid_altitudes, leg
ax2 = ax.twinx()
# ax2.legend(handles=legend_elements, bbox_to_anchor=(0.93, 0.7), loc='upper right')
# ax2.annotate("Filled symbol = significant trend ", xy=(0.85, 0.5), xycoords='axes fraction', fontsize=7)
- ax2.legend(handles=legend_elements, loc='upper right', prop={'size': legend_size})
- ax2.annotate("Filled symbol =\nsignificant trend \nw.r.t $\mathcal{M}_0$", xy=(0.6, 0.85), xycoords='axes fraction', fontsize=fontsize)
+ ax2.legend(handles=legend_elements, loc='center left', prop={'size': legend_size})
+ # ax2.annotate("Filled symbol =\nsignificant trend \nw.r.t $\mathcal{M}_0$", xy=(0.6, 0.85), xycoords='axes fraction', fontsize=fontsize)
+ ax2.annotate('\n'.join(filled_marker_legend_list), xy=(0.23, 0.43), xycoords='axes fraction', fontsize=fontsize)
ax2.set_yticks([])
diff --git a/experiment/paper_past_snow_loads/study_visualizer_for_non_stationary_trends.py b/experiment/paper_past_snow_loads/study_visualizer_for_non_stationary_trends.py
index fb37f57b56dd0dac12684329ad2dadd2ef56fdb1..efd0742ef5d8b823730ef6faaece0e56c499fe7f 100644
--- a/experiment/paper_past_snow_loads/study_visualizer_for_non_stationary_trends.py
+++ b/experiment/paper_past_snow_loads/study_visualizer_for_non_stationary_trends.py
@@ -9,6 +9,7 @@ from cached_property import cached_property
from experiment.eurocode_data.massif_name_to_departement import massif_name_to_eurocode_region
from experiment.eurocode_data.utils import EUROCODE_QUANTILE, EUROCODE_RETURN_LEVEL_STR
from experiment.meteo_france_data.plot.create_shifted_cmap import get_shifted_map, get_colors
+from experiment.meteo_france_data.scm_models_data.abstract_extended_study import AbstractExtendedStudy
from experiment.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
from experiment.meteo_france_data.scm_models_data.visualization.study_visualization.study_visualizer import \
StudyVisualizer
@@ -190,7 +191,8 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
ax.get_xaxis().set_visible(True)
ax.set_xticks([])
ax.set_xlabel('Altitude = {}m'.format(self.study.altitude), fontsize=15)
- self.plot_name = 'tdlr_trends_w' + 'o' if not add_colorbar else '' + '_colorbar'
+ middle_word = 'o' if (not add_colorbar and self.study.altitude == 2700) else ''
+ self.plot_name = 'tdlr_trends_w' + middle_word + '_colorbar'
self.show_or_save_to_file(add_classic_title=False, tight_layout=True, no_title=True,
dpi=500)
plt.close()
@@ -366,11 +368,29 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
for eurocode, uncertainty in eurocode_and_uncertainties])
return 100 * np.mean(a, axis=0)
- # Part 3 - Zeros
-
- # Part 4 - QQPLOT
+ # Part 3 - QQPLOT
def qqplot(self, massif_name, color=None):
trend_test = self.massif_name_to_trend_test_that_minimized_aic[massif_name]
marker = self.massif_name_to_marker_style[massif_name]
- trend_test.qqplot_wrt_standard_gumbel(marker, color)
\ No newline at end of file
+ trend_test.qqplot_wrt_standard_gumbel(marker, color)
+
+ # Part 4 - Trend plot
+
+ @property
+ def altitude(self):
+ return self.study.altitude
+
+ def trend_summary_labels(self):
+ return 'altitude', ''
+
+ def trend_summary_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)
+ massif_name_to_region_name = AbstractExtendedStudy.massif_name_to_region_name
+ mean_relative = np.mean(np.array(list(self.massif_name_to_relative_change_value.values())))
+
+ return self.altitude, percentage_decrease, percentage_decrease_significative
\ No newline at end of file
diff --git a/experiment/trend_analysis/univariate_test/extreme_trend_test/abstract_gev_trend_test.py b/experiment/trend_analysis/univariate_test/extreme_trend_test/abstract_gev_trend_test.py
index 6bc14d7a2d395206b91641a1549ec8bac4a84da4..9a07a18e1c5fe518fe1cf04b80ac03ae18b447d3 100644
--- a/experiment/trend_analysis/univariate_test/extreme_trend_test/abstract_gev_trend_test.py
+++ b/experiment/trend_analysis/univariate_test/extreme_trend_test/abstract_gev_trend_test.py
@@ -210,6 +210,7 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
# Some display function
def qqplot_wrt_standard_gumbel(self, marker, color=None):
+ size = 20
# Standard Gumbel quantiles
standard_gumbel_distribution = GevParams(loc=0, scale=1, shape=0)
n = len(self.years)
@@ -220,9 +221,9 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
plt.plot(standard_gumbel_quantiles, constrained_empirical_quantiles, 'x', label='Gumbel model')
plt.plot(standard_gumbel_quantiles, unconstrained_empirical_quantiles, linestyle='None',
label='Selected model', **marker)
- plt.xlabel("Standard Gumbel quantiles")
- plt.ylabel("Empirical quantiles")
- plt.legend()
+ plt.xlabel("Standard Gumbel quantiles", fontsize=15)
+ plt.ylabel("Empirical quantiles", fontsize=15)
+ plt.legend(loc='upper left', prop={'size': size})
plt.show()
def compute_empirical_quantiles(self, estimator):