An error occurred while loading the file. Please try again.
-
Le Roux Erwan authored35620a94
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
import matplotlib.pyplot as plt
import numpy as np
from cached_property import cached_property
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 projects.altitude_spatial_model.altitudes_fit.altitudes_studies import AltitudesStudies
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):
def __init__(self, study_class, altitudes, spatial_transformation_class=None, temporal_transformation_class=None,
**kwargs_study):
super().__init__(study_class, altitudes, spatial_transformation_class, temporal_transformation_class,
**kwargs_study)
# self.altitudes_for_temporal_hypothesis = [min(self.altitudes), 2100, max(self.altitudes)]
self.altitudes_for_temporal_hypothesis = [600, 1500, 2400, 3300]
def plot_gev_params_against_altitude(self):
for j, param_name in enumerate(GevParams.PARAM_NAMES + [100]):
ax = plt.gca()
massif_name_to_linear_coef = {}
massif_name_to_r2_score = {}
massif_names = self.study.all_massif_names()[:]
for i in range(8):
for massif_name in massif_names[i::8]:
linear_coef, _, r2 = self._plot_gev_params_against_altitude_one_massif(ax, massif_name, param_name)
massif_name_to_linear_coef[massif_name] = 100 * linear_coef[0]
massif_name_to_r2_score[massif_name] = str(round(r2, 2))
print(param_name, np.mean([c for c in massif_name_to_linear_coef.values()]))
# Display x label
xticks = [1000 * i for i in range(1, 5)]
ax.set_xticks(xticks)
fontsize_label = 15
ax.tick_params(labelsize=fontsize_label)
# ax.set_xlabel('Altitude')
# Compute the y label
if param_name in GevParams.PARAM_NAMES:
ylabel = GevParams.full_name_from_param_name(param_name) + ' parameter'
else:
ylabel = '{}-year return levels'.format(param_name)
# Add units
if param_name == GevParams.SHAPE:
unit = 'no unit'
else:
unit = self.study.variable_unit
ylabel += ' ({})'.format(unit)
# Display the y label on the twin axis
if param_name in [100, GevParams.SCALE]:
ax2 = ax.twinx()
ax2.set_yticks(ax.get_yticks())
ax2.set_ylim(ax.get_ylim())
ax2.set_ylabel(ylabel, fontsize=fontsize_label)
ax2.tick_params(labelsize=fontsize_label)
ax.set_yticks([])
tight_layout = False
else:
ax.tick_params(labelsize=fontsize_label)
tight_layout = True
ax.set_ylabel(ylabel, fontsize=fontsize_label)
# Make room for the ylabel
if param_name == 100:
plt.gcf().subplots_adjust(right=0.85)
plot_name = '{} change with altitude'.format(param_name)
# # Display the legend
# ax.legend(labelspacing=2.5, ncol=8, handlelength=12, markerscale=0.7, bbox_to_anchor=(1.05, 1), loc='upper left',
# prop={'size': 2}, fontsize='x-large')
# plt.gcf().subplots_adjust(right=0.15)
# ax.set_yticks([])
# ax.set_ylabel('')
# plt.show()
self.show_or_save_to_file(plot_name, no_title=True, tight_layout=tight_layout, show=False)
ax.clear()
plt.close()
# Plot map of slope for each massif
visualizer = StudyVisualizer(study=self.study, show=False, save_to_file=True)
idx = 8 if param_name == GevParams.SHAPE else 1
label = 'Altitude gradient for the {}'.format(ylabel[:-idx] + '/100m)')
gev_param_name_to_graduation = {
GevParams.LOC: 0.5,
GevParams.SCALE: 0.1,
GevParams.SHAPE: 0.01,
100: 1,
}
if param_name == GevParams.SHAPE:
print(massif_name_to_linear_coef)
visualizer.plot_map(cmap=plt.cm.coolwarm, fit_method=self.study.fit_method,
graduation=gev_param_name_to_graduation[param_name],
label=label, massif_name_to_value=massif_name_to_linear_coef,
plot_name=label.replace('/', ' every '), add_x_label=False,
negative_and_positive_values=param_name == GevParams.SHAPE,
add_colorbar=True,
massif_name_to_text=massif_name_to_r2_score,
)
plt.close()
def _plot_gev_params_against_altitude_one_massif(self, ax, massif_name, param_name):
altitudes = []
params = []
# confidence_intervals = []
for altitude, study in self.altitude_to_study.items():
if massif_name in study.massif_name_to_stationary_gev_params:
gev_params = study.massif_name_to_stationary_gev_params[massif_name]
altitudes.append(altitude)
if param_name in GevParams.PARAM_NAMES:
param = gev_params.to_dict()[param_name]
else:
assert isinstance(param_name, int)
param = gev_params.return_level(return_period=param_name)
params.append(param)
# confidence_intervals.append(gev_params.param_name_to_confidence_interval[param_name])
massif_id = self.study.all_massif_names().index(massif_name)
plot_against_altitude(altitudes, ax, massif_id, massif_name, params, fill=False)
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))
if __name__ == '__main__':
altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300]
altitudes = [600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900]
# altitudes = paper_altitudes
# altitudes = [1800, 2100]
visualizer = PointwiseGevStudyVisualizer(SafranSnowfall1Day, altitudes=altitudes)
visualizer.plot_gev_params_against_altitude()
# visualizer.plot_gev_params_against_time_for_all_altitudes()
# visualizer.plot_gev_params_against_time_for_all_massifs()
# visualizer.plot_time_derivative_against_time()