An error occurred while loading the file. Please try again.
-
Le Roux Erwan authored0f030611
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
import matplotlib.pyplot as plt
import os.path as op
from itertools import chain
import numpy as np
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 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
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)
def plot_gev_params_against_altitude(self):
legend = False
elevation_as_xaxis = False
param_names = GevParams.PARAM_NAMES + [100]
if legend:
param_names = param_names[:1]
for j, param_name in enumerate(param_names):
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,
elevation_as_xaxis,
legend=legend)
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
elevation_ticks = [500 * i for i in range(1, 8)]
if elevation_as_xaxis:
ax.set_xticks(elevation_ticks)
else:
ax.set_yticks(elevation_ticks)
if j == 2:
ax.set_xlim(right=0.45)
fontsize_label = 15
ax.tick_params(labelsize=fontsize_label)
# Compute the y label
if param_name in GevParams.PARAM_NAMES:
value_label = GevParams.full_name_from_param_name(param_name) + ' parameter'
else:
value_label = '{}-year return levels'.format(param_name)
# Add units
if param_name == GevParams.SHAPE:
unit = 'no unit'
else:
unit = self.study.variable_unit
value_label += ' ({})'.format(unit)
value_label = value_label.capitalize()
if elevation_as_xaxis:
# 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(value_label, 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(value_label, fontsize=fontsize_label)
# Make room for the ylabel
if param_name == 100:
plt.gcf().subplots_adjust(right=0.85)
else:
if param_name in [GevParams.LOC, GevParams.SCALE]:
ax2 = ax.twiny()
ax2.set_xticks(ax.get_xticks())
ax.set_xticks([])
else:
ax2 = ax
ax2.set_xlim(ax.get_xlim())
ax2.set_xlabel(value_label, fontsize=fontsize_label)
ax2.tick_params(labelsize=fontsize_label)
if param_name in [100, GevParams.SCALE]:
ax3 = ax2.twinx()
ax3.set_yticks(ax2.get_yticks())
ax3.set_ylim(ax2.get_ylim())
ax3.tick_params(labelsize=fontsize_label)
ax2.set_yticks([])
ax3.set_ylabel('Elevation (m)', fontsize=fontsize_label)
else:
ax.set_ylabel('Elevation (m)', fontsize=fontsize_label)
tight_layout = False
plot_name = '{} change with altitude'.format(param_name)
if isinstance(self.study, AbstractAdamontStudy):
plot_name = op.join(plot_name, _gcm_rcm_couple_adamont_v2_to_full_name[self.study.gcm_rcm_couple])
# # Display the legend
if 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')
# ax.legend(labelspacing=1, ncol=8, handlelength=5, bbox_to_anchor=(1.05, 1), loc='upper left',
# prop={'size': 4}, fontsize='xx-large', columnspacing=0.5)
ax.legend(ncol=8, bbox_to_anchor=(1.05, 1), loc='upper left',
prop={'size': 3.5}, handlelength=5, fontsize='xx-large', columnspacing=0.5,
handletextpad=0.5)
# handles, labels = ax.get_legend_handles_labels()
# print(type(handles))
# handles = np.array(handles).reshape((3, 8)).transpose().flatten()
# labels = np.array(handles).reshape((3, 8)).transpose().flatten()
# ax.legend(handles, labels)
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
the = ' the' if param_name in GevParams.PARAM_NAMES else ''
label = 'Elevation gradient for\n{} {}'.format(the, value_label[:-idx] + '/100m)')
plot_name = label.replace('/', ' every ')
if isinstance(self.study, AbstractAdamontStudy):
plot_name = op.join(plot_name, _gcm_rcm_couple_adamont_v2_to_full_name[self.study.gcm_rcm_couple])
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,
graduation=gev_param_name_to_graduation[param_name],
label=label, massif_name_to_value=massif_name_to_linear_coef,
plot_name=plot_name, add_x_label=False,
negative_and_positive_values=param_name == GevParams.SHAPE,
add_colorbar=True,
massif_name_to_text=massif_name_to_r2_score,
fontsize_label=13,
)
plt.close()
def _plot_gev_params_against_altitude_one_massif(self, ax, massif_name, param_name, elevation_as_xaxis,
legend=False):
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)
massif_id = self.study.all_massif_names().index(massif_name)
plot_against_altitude(altitudes, ax, massif_id, massif_name, params, fill=False,
elevation_as_xaxis=elevation_as_xaxis,
legend=legend)
return fit_linear_regression(altitudes, params)
def main_paper2():
altitudes = list(chain.from_iterable(altitudes_for_groups))
# 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()
def main_paper3():
altitudes = list(chain.from_iterable(altitudes_for_groups))
# altitudes = [1200, 1500, 1800]
for scenario in rcp_scenarios[2:]:
gcm_rcm_couples = get_gcm_rcm_couples(scenario)
gcm_rcm_couples =[('CNRM-CM5', 'CCLM4-8-17')]
for gcm_rcm_couple in gcm_rcm_couples:
visualizer = PointwiseGevStudyVisualizer(AdamontSnowfall, altitudes=altitudes, scenario=scenario,
gcm_rcm_couple=gcm_rcm_couple)
visualizer.plot_gev_params_against_altitude()
if __name__ == '__main__':
main_paper3()