From fff8059acccc38b832f2a6aac9e749409f240010 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Mon, 6 Apr 2020 21:36:58 +0200
Subject: [PATCH] [contrasting project] add figure 1
---
.../scm_models_data/abstract_study.py | 12 +++-
extreme_fit/distribution/gev/gev_params.py | 3 +
.../gorman_figures/__init__.py | 0
.../figure1_mean_ratio_return_level_ratio.py | 67 +++++++++++++++++++
.../figure3_daily snowfall fraction.py} | 47 +++++++------
5 files changed, 103 insertions(+), 26 deletions(-)
create mode 100644 projects/contrasting_trends_in_snow_loads/gorman_figures/__init__.py
create mode 100644 projects/contrasting_trends_in_snow_loads/gorman_figures/figure1_mean_ratio_return_level_ratio.py
rename projects/contrasting_trends_in_snow_loads/{snowfall fraction of total precipitation/daily snowfall fraction.py => gorman_figures/figure3_daily snowfall fraction.py} (69%)
diff --git a/extreme_data/meteo_france_data/scm_models_data/abstract_study.py b/extreme_data/meteo_france_data/scm_models_data/abstract_study.py
index 54e7da27..adbebbdf 100644
--- a/extreme_data/meteo_france_data/scm_models_data/abstract_study.py
+++ b/extreme_data/meteo_france_data/scm_models_data/abstract_study.py
@@ -238,12 +238,20 @@ class AbstractStudy(object):
ordered_index = [self.year_to_annual_maxima_index[year][i] for year in years]
massif_name_to_annual_maxima_ordered_index[massif_name] = ordered_index
return massif_name_to_annual_maxima_ordered_index
+
+ @cached_property
+ def massif_name_to_annual_maxima(self):
+ massif_name_to_annual_maxima = OrderedDict()
+ for i, massif_name in enumerate(self.study_massif_names):
+ maxima = np.array([self.year_to_annual_maxima[year][i] for year in self.ordered_years])
+ massif_name_to_annual_maxima[massif_name] = maxima
+ return massif_name_to_annual_maxima
@cached_property
def massif_name_to_annual_maxima_ordered_years(self):
massif_name_to_annual_maxima_ordered_years = OrderedDict()
- for i, massif_name in enumerate(self.study_massif_names):
- maxima = np.array([self.year_to_annual_maxima[year][i] for year in self.ordered_years])
+ for massif_name in self.study_massif_names:
+ maxima = self.massif_name_to_annual_maxima[massif_name]
annual_maxima_ordered_index = np.argsort(maxima)
annual_maxima_ordered_years = [self.ordered_years[idx] for idx in annual_maxima_ordered_index]
massif_name_to_annual_maxima_ordered_years[massif_name] = annual_maxima_ordered_years
diff --git a/extreme_fit/distribution/gev/gev_params.py b/extreme_fit/distribution/gev/gev_params.py
index 79dbe892..ce34b36b 100644
--- a/extreme_fit/distribution/gev/gev_params.py
+++ b/extreme_fit/distribution/gev/gev_params.py
@@ -37,6 +37,9 @@ class GevParams(AbstractExtremeParams):
else:
return r.qgev(p, self.location, self.scale, self.shape)[0]
+ def return_level(self, return_period):
+ return self.quantile(1 - 1 / return_period)
+
def density(self, x):
if self.has_undefined_parameters:
return np.nan
diff --git a/projects/contrasting_trends_in_snow_loads/gorman_figures/__init__.py b/projects/contrasting_trends_in_snow_loads/gorman_figures/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/projects/contrasting_trends_in_snow_loads/gorman_figures/figure1_mean_ratio_return_level_ratio.py b/projects/contrasting_trends_in_snow_loads/gorman_figures/figure1_mean_ratio_return_level_ratio.py
new file mode 100644
index 00000000..fd5608f8
--- /dev/null
+++ b/projects/contrasting_trends_in_snow_loads/gorman_figures/figure1_mean_ratio_return_level_ratio.py
@@ -0,0 +1,67 @@
+from collections import OrderedDict
+
+import numpy as np
+
+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
+from extreme_data.meteo_france_data.scm_models_data.visualization.create_shifted_cmap import get_shifted_map, \
+ get_colors, shiftedColorMap
+from extreme_fit.estimator.margin_estimator.utils import fitted_stationary_gev
+from extreme_fit.model.margin_model.linear_margin_model.abstract_temporal_linear_margin_model import \
+ TemporalMarginFitMethod
+
+import matplotlib.pyplot as plt
+
+from extreme_trend.visualizers.study_visualizer_for_non_stationary_trends import StudyVisualizerForNonStationaryTrends
+
+
+def mean_ratio(altitude):
+ pass
+
+
+def massif_name_to_return_level(altitude, year_min, year_max, return_period):
+ study = SafranSnowfall1Day(altitude=altitude, year_min=year_min, year_max=year_max)
+ massif_name_to_return_levels = OrderedDict()
+ for massif_name, maxima in study.massif_name_to_annual_maxima.items():
+ gev_param = fitted_stationary_gev(maxima, fit_method=TemporalMarginFitMethod.extremes_fevd_mle)
+ return_level = gev_param.return_level(return_period=return_period)
+ massif_name_to_return_levels[massif_name] = return_level
+ return massif_name_to_return_levels
+
+
+def plot_return_level_ratio(altitude, return_period=30, year_min=1959, year_middle=1989, year_max=2019):
+ massif_name_to_return_level_past = massif_name_to_return_level(altitude, year_min, year_middle, return_period)
+ massif_name_to_return_level_recent = massif_name_to_return_level(altitude, year_middle + 1, year_max,
+ return_period)
+ massif_name_to_ratio = {m: v / massif_name_to_return_level_past[m] for m, v in
+ massif_name_to_return_level_recent.items()}
+ max_ratio = max([e for e in massif_name_to_ratio.values()])
+ min_ratio = min([e for e in massif_name_to_ratio.values()])
+ midpoint = (max_ratio - 1.0) / (max_ratio - 0)
+ num = 3
+ ticks = np.concatenate([np.linspace(0, midpoint, num), np.linspace(midpoint, 1, num)])
+ labels = np.concatenate([np.linspace(min_ratio, 1.0, num), np.linspace(1.0, max_ratio, num)])
+ cmap = shiftedColorMap(plt.cm.bwr, midpoint=midpoint, name='shifted')
+ ax = plt.gca()
+ print(massif_name_to_ratio)
+ massif_name_to_color = {m: get_colors([v], cmap, min_ratio, max_ratio)[0]
+ for m, v in massif_name_to_ratio.items()}
+ ticks_values_and_labels = ticks, labels
+ AbstractStudy.visualize_study(ax=ax,
+ massif_name_to_value=massif_name_to_ratio,
+ massif_name_to_color=massif_name_to_color,
+ replace_blue_by_white=True,
+ axis_off=False,
+ cmap=cmap,
+ show_label=False,
+ add_colorbar=True,
+ show=False,
+ vmin=min_ratio,
+ vmax=max_ratio,
+ ticks_values_and_labels=ticks_values_and_labels
+ )
+ plt.show()
+
+
+if __name__ == '__main__':
+ plot_return_level_ratio(altitude=1800)
diff --git a/projects/contrasting_trends_in_snow_loads/snowfall fraction of total precipitation/daily snowfall fraction.py b/projects/contrasting_trends_in_snow_loads/gorman_figures/figure3_daily snowfall fraction.py
similarity index 69%
rename from projects/contrasting_trends_in_snow_loads/snowfall fraction of total precipitation/daily snowfall fraction.py
rename to projects/contrasting_trends_in_snow_loads/gorman_figures/figure3_daily snowfall fraction.py
index 8403554b..6a1b5cf1 100644
--- a/projects/contrasting_trends_in_snow_loads/snowfall fraction of total precipitation/daily snowfall fraction.py
+++ b/projects/contrasting_trends_in_snow_loads/gorman_figures/figure3_daily snowfall fraction.py
@@ -8,51 +8,50 @@ from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranR
Altitudes_groups = [[300, 600, 900], [1200, 1500, 1800], [2100, 2400, 2700], [3000, 3300, 3600]]
-def compute_smoother_snowfall_fraction(altitudes, year_min, year_max):
+def compute_smoother_snowfall_fraction(altitudes, year_min, year_max, lim):
all_time_series = np.concatenate([get_time_series(altitude, year_min, year_max) for altitude in altitudes], axis=1)
temperature_array, snowfall_array, rainfall_array = all_time_series
- paper_setting = False
- if paper_setting:
- lim = 10
- sigma = 0.5
- else:
- lim = 7.5
- sigma = 1
- bins = np.linspace(-lim, lim, num=(4 * 2 * lim) + 1)
+ nb_bins_for_one_celsius = 4
+ bins = np.linspace(-lim, lim, num=int((nb_bins_for_one_celsius * 2 * lim) + 1))
inds = np.digitize(temperature_array, bins)
snowfall_fractions = []
for j in range(1, len(bins)):
mask = inds == j
fraction = np.mean(snowfall_array[mask]) / np.mean(snowfall_array[mask] + rainfall_array[mask])
snowfall_fractions.append(fraction)
- new_snowfall_fractions = 100 * gaussian_filter1d(snowfall_fractions, sigma=sigma)
x = bins[:-1] + 0.125
- return new_snowfall_fractions, sigma, x
+ return np.array(snowfall_fractions), x
def daily_snowfall_fraction(ax, altitudes, year_min=1959, year_max=2019):
- snowfall_fractions, sigma, x = compute_smoother_snowfall_fraction(altitudes, year_min, year_max)
+ snowfall_fractions, x = compute_smoother_snowfall_fraction(altitudes, year_min, year_max, lim=7)
+ sigma = 1.0
+ snowfall_fractions = 100 * gaussian_filter1d(snowfall_fractions, sigma=sigma)
# Plot results
label = '{}-{} at {}'.format(year_min, year_max, ' & '.join(['{} m'.format(a) for a in altitudes]))
ax.plot(x, snowfall_fractions, label=label, linewidth=4)
ax.set_ylabel('Snowfall fraction (%)')
- ax.set_xlabel('Daily surface air temperature (Celsius)')
- ax.set_title('Snowfall fraction is smoothed with a gaussian filter with standard deviation={}'.format(sigma))
- ax.grid(b=True)
- ax.legend()
+ end_plot(ax, sigma)
def snowfall_fraction_difference_between_periods(ax, altitudes, year_min=1959, year_middle=1989, year_max=2019):
- snowfall_fractions_past, sigma, x = compute_smoother_snowfall_fraction(altitudes, year_min, year_middle)
- snowfall_fractions_recent, sigma, x = compute_smoother_snowfall_fraction(altitudes, year_middle + 1, year_max)
- mask = snowfall_fractions_past == 0
+ lim = 4
+ snowfall_fractions_past, x = compute_smoother_snowfall_fraction(altitudes, year_min, year_middle, lim=lim)
+ snowfall_fractions_recent, x = compute_smoother_snowfall_fraction(altitudes, year_middle + 1, year_max, lim=lim)
+ # Ensure that we do not divide by zero
+ mask = snowfall_fractions_past > 0
x = x[mask]
snowfall_ratio = snowfall_fractions_recent[mask] / snowfall_fractions_past[mask]
+ sigma = 2.0
+ snowfall_ratio = gaussian_filter1d(snowfall_ratio, sigma=sigma)
label = '{}'.format(' & '.join(['{} m'.format(a) for a in altitudes]))
ax.plot(x, snowfall_ratio, label=label, linewidth=4)
- ax.set_ylabel('Ratio of Snowfall fraction {}-{} divided by Snowfall fraction {}-{}'.format(year_min, year_middle,
- year_middle + 1,
- year_max))
+ ax.set_ylabel('Ratio of Snowfall fraction {}-{}\n'
+ 'divided by Snowfall fraction {}-{}'.format(year_middle + 1, year_max, year_min, year_middle))
+ end_plot(ax, sigma)
+
+
+def end_plot(ax, sigma):
ax.set_xlabel('Daily surface air temperature (Celsius)')
ax.set_title('Snowfall fraction is smoothed with a gaussian filter with standard deviation={}'.format(sigma))
ax.grid(b=True)
@@ -94,6 +93,6 @@ def daily_snowfall_fraction_wrt_time():
if __name__ == '__main__':
- # daily_snowfall_fraction_wrt_altitude(fast=True)
+ daily_snowfall_fraction_wrt_altitude(fast=False)
# daily_snowfall_fraction_wrt_time()
- ratio_of_past_and_recent_daily_snowfall_fraction_wrt_altitude()
+ # ratio_of_past_and_recent_daily_snowfall_fraction_wrt_altitude()
--
GitLab