From 5b9f98f2b85012564ab54b8918c4bf5aa9d56d12 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Fri, 8 Jan 2021 09:52:29 +0100
Subject: [PATCH] [contrasting] harmonize coherence plot & preliminary plot.
add verification script for return levels.
---
.../plots/plot_coherence_curves.py | 7 +++-
.../preliminary_analysis.py | 4 ++-
.../verification_return_level.py | 34 +++++++++++++++++++
3 files changed, 43 insertions(+), 2 deletions(-)
create mode 100644 projects/altitude_spatial_model/verification_return_level.py
diff --git a/projects/altitude_spatial_model/altitudes_fit/plots/plot_coherence_curves.py b/projects/altitude_spatial_model/altitudes_fit/plots/plot_coherence_curves.py
index 21692283..7e47ac6e 100644
--- a/projects/altitude_spatial_model/altitudes_fit/plots/plot_coherence_curves.py
+++ b/projects/altitude_spatial_model/altitudes_fit/plots/plot_coherence_curves.py
@@ -40,14 +40,19 @@ def plot_coherence_curves(massif_names, visualizer_list: List[AltitudesStudiesVi
is_altitudinal,
year)
label = labels[i]
+ label = label.capitalize()
# Set labels
fontsize_label = 15
+
+ for a in [ax, ax2, ax3]:
+ a.tick_params(labelsize=fontsize_label)
+
if elevation_as_xaxis:
ax3.set_xlabel('Elevation (m)', fontsize=fontsize_label)
ax2.set_ylabel(label, fontsize=fontsize_label)
else:
- ax2.set_ylabel('Elevation(m)', fontsize=fontsize_label)
+ ax2.set_ylabel('Elevation (m)', fontsize=fontsize_label)
if i == 3:
ax.set_xlabel(label, fontsize=fontsize_label)
else:
diff --git a/projects/altitude_spatial_model/preliminary_analysis.py b/projects/altitude_spatial_model/preliminary_analysis.py
index 679f99a8..44260de4 100644
--- a/projects/altitude_spatial_model/preliminary_analysis.py
+++ b/projects/altitude_spatial_model/preliminary_analysis.py
@@ -45,11 +45,13 @@ class PointwiseGevStudyVisualizer(AltitudesStudies):
print(param_name, np.mean([c for c in massif_name_to_linear_coef.values()]))
# Display x label
- elevation_ticks = [1000 * i for i in range(1, 4)]
+ 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)
diff --git a/projects/altitude_spatial_model/verification_return_level.py b/projects/altitude_spatial_model/verification_return_level.py
new file mode 100644
index 00000000..6d5bc92a
--- /dev/null
+++ b/projects/altitude_spatial_model/verification_return_level.py
@@ -0,0 +1,34 @@
+import numpy as np
+
+from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day
+
+
+def compute_return_level():
+ count_all = 0
+ count_exceedance = 0
+ diff_list = []
+ for altitude in [600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600][:]:
+ study = SafranSnowfall1Day(altitude=altitude, year_max=2008)
+ for massif_name in study.study_massif_names:
+ if massif_name in study.massif_name_to_stationary_gev_params:
+
+ gev_params = study.massif_name_to_stationary_gev_params[massif_name]
+ return_level100 = gev_params.return_level(return_period=100)
+ annual_maxima = study.massif_name_to_annual_maxima[massif_name]
+ annual_maxima = sorted(annual_maxima)
+ max_annual_maxima = max(annual_maxima)
+
+ count_all += 1
+ relative_diff =100 * (return_level100 - max_annual_maxima) / max_annual_maxima
+ diff_list.append(relative_diff)
+
+ if return_level100 < max_annual_maxima:
+ count_exceedance += 1
+ print(altitude, massif_name, annual_maxima[-2:], return_level100, gev_params.shape)
+ percent_of_exceedance = 100 * count_exceedance / count_all
+ print(percent_of_exceedance)
+ mean_relatif_diff = np.mean(diff_list)
+ print(mean_relatif_diff)
+
+if __name__ == '__main__':
+ compute_return_level()
\ No newline at end of file
--
GitLab