From 1ec6c1a3974db3147e6275f9ae24d0326eb0ab31 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Wed, 29 Apr 2020 17:45:42 +0200
Subject: [PATCH] [contrasting] add snowfall_plot.py
---
.../visualization/plot_utils.py | 17 ++++-
.../visualization/study_visualizer.py | 1 +
.../independent_margin_function.py | 8 +++
.../parametric_margin_function.py | 6 +-
.../function/param_function/param_function.py | 5 +-
extreme_trend/abstract_gev_trend_test.py | 20 ++++++
.../altitudes_fit/altitudes_studies.py | 18 +----
.../main_snowfall_article.py | 7 +-
.../snowfall_plot.py | 72 +++++++++++++++++++
.../study_visualizer_for_mean_values.py | 2 -
.../validation_plot.py | 5 +-
11 files changed, 134 insertions(+), 27 deletions(-)
create mode 100644 projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/snowfall_plot.py
diff --git a/extreme_data/meteo_france_data/scm_models_data/visualization/plot_utils.py b/extreme_data/meteo_france_data/scm_models_data/visualization/plot_utils.py
index 341b53ce..c16465d2 100644
--- a/extreme_data/meteo_france_data/scm_models_data/visualization/plot_utils.py
+++ b/extreme_data/meteo_france_data/scm_models_data/visualization/plot_utils.py
@@ -1,10 +1,25 @@
import matplotlib.pyplot as plt
+import matplotlib.colors as mcolors
from extreme_data.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
from extreme_data.meteo_france_data.scm_models_data.visualization.create_shifted_cmap import \
ticks_values_and_labels_for_percentages, get_shifted_map, get_colors
+def plot_against_altitude(altitudes, ax, massif_id, massif_name, values):
+ di = massif_id // 8
+ if di == 0:
+ linestyle = '-'
+ elif di == 1:
+ linestyle = 'dotted'
+ else:
+ linestyle = '--'
+ colors = list(mcolors.TABLEAU_COLORS)
+ colors[-3:-1] = [] # remove gray and olive
+ color = colors[massif_id % 8]
+ ax.plot(altitudes, values, color=color, linewidth=2, label=massif_name, linestyle=linestyle)
+
+
def load_plot(cmap, graduation, label, massif_name_to_value, altitude, fit_method):
max_abs_change = max([abs(e) for e in massif_name_to_value.values()])
ticks, labels = ticks_values_and_labels_for_percentages(graduation=graduation, max_abs_change=max_abs_change)
@@ -33,4 +48,4 @@ def load_plot(cmap, graduation, label, massif_name_to_value, altitude, fit_metho
ax.get_xaxis().set_visible(True)
ax.set_xticks([])
ax.set_xlabel('Altitude = {}m'.format(altitude), fontsize=15)
- ax.set_title('Fit method is {}'.format(fit_method))
\ No newline at end of file
+ ax.set_title('Fit method is {}'.format(fit_method))
diff --git a/extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py b/extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py
index 38481079..ed322e21 100644
--- a/extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py
+++ b/extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py
@@ -1,4 +1,5 @@
import os
+import matplotlib.colors as mcolors
import os.path as op
from collections import OrderedDict
from multiprocessing.pool import Pool
diff --git a/extreme_fit/function/margin_function/independent_margin_function.py b/extreme_fit/function/margin_function/independent_margin_function.py
index b1d6b54a..79db8d38 100644
--- a/extreme_fit/function/margin_function/independent_margin_function.py
+++ b/extreme_fit/function/margin_function/independent_margin_function.py
@@ -30,6 +30,14 @@ class IndependentMarginFunction(AbstractMarginFunction):
for param_name, param_function in self.param_name_to_param_function.items()}
return self.params_class.from_dict(params)
+ def get_first_derivative_param(self, coordinate: np.ndarray, is_transformed: bool, dim: int = 0):
+ transformed_coordinate = coordinate if is_transformed else self.transform(coordinate)
+ return {
+ param_name: param_function.get_first_derivative_param_value(transformed_coordinate, dim)
+ for param_name, param_function in self.param_name_to_param_function.items()
+ }
+
+
diff --git a/extreme_fit/function/margin_function/parametric_margin_function.py b/extreme_fit/function/margin_function/parametric_margin_function.py
index ccf2a83d..17e5dcfd 100644
--- a/extreme_fit/function/margin_function/parametric_margin_function.py
+++ b/extreme_fit/function/margin_function/parametric_margin_function.py
@@ -67,6 +67,10 @@ class ParametricMarginFunction(IndependentMarginFunction):
return self.coordinates.temporal_coordinates.transformation.transform_array(np.array([self.starting_point]))
def get_params(self, coordinate: np.ndarray, is_transformed: bool = True) -> GevParams:
+ coordinate = self.shift_coordinates_if_needed(coordinate, is_transformed)
+ return super().get_params(coordinate, is_transformed=is_transformed)
+
+ def shift_coordinates_if_needed(self, coordinate, is_transformed):
if self.starting_point is not None:
starting_point = self.transformed_starting_point if is_transformed else self.starting_point
# Shift temporal coordinate to enable to model temporal trend with starting point
@@ -74,7 +78,7 @@ class ParametricMarginFunction(IndependentMarginFunction):
assert 0 <= self.coordinates.idx_temporal_coordinates < len(coordinate)
if coordinate[self.coordinates.idx_temporal_coordinates] < starting_point:
coordinate[self.coordinates.idx_temporal_coordinates] = starting_point
- return super().get_params(coordinate, is_transformed=is_transformed)
+ return coordinate
@classmethod
def from_coef_dict(cls, coordinates: AbstractCoordinates, param_name_to_dims: Dict[str, List[int]],
diff --git a/extreme_fit/function/param_function/param_function.py b/extreme_fit/function/param_function/param_function.py
index 470b1c69..575ff747 100644
--- a/extreme_fit/function/param_function/param_function.py
+++ b/extreme_fit/function/param_function/param_function.py
@@ -8,7 +8,10 @@ class AbstractParamFunction(object):
OUT_OF_BOUNDS_ASSERT = True
def get_param_value(self, coordinate: np.ndarray) -> float:
- pass
+ raise NotImplementedError
+
+ def get_first_derivative_param_value(self, coordinate: np.ndarray, dim: int) -> float:
+ raise NotImplementedError
class ConstantParamFunction(AbstractParamFunction):
diff --git a/extreme_trend/abstract_gev_trend_test.py b/extreme_trend/abstract_gev_trend_test.py
index 3ba90563..eac300d0 100644
--- a/extreme_trend/abstract_gev_trend_test.py
+++ b/extreme_trend/abstract_gev_trend_test.py
@@ -238,6 +238,13 @@ class AbstractGevTrendTest(object):
is_transformed=False)
return gev_params_year
+ def get_unconstrained_time_derivative_gev_params(self, year):
+ gev_params_year = self.unconstrained_estimator.function_from_fit.get_first_derivative_param(
+ coordinate=np.array([year]),
+ is_transformed=False,
+ dim=0)
+ return gev_params_year
+
def linear_extension(self, ax, q, quantiles, sorted_maxima):
# Extend the curve linear a bit if the return period 50 is not in the quantiles
def compute_slope_intercept(x, y):
@@ -387,3 +394,16 @@ class AbstractGevTrendTest(object):
average_mean_value = np.mean(mean_values)
assert isinstance(average_mean_value, float)
return average_mean_value
+
+ # Derivative mean value in some year
+
+ def mean_value(self, year):
+ return self.get_unconstrained_gev_params(year).mean
+
+ def first_derivative_mean_value(self, year):
+ param_name_to_value = self.get_unconstrained_time_derivative_gev_params(year)
+ assert param_name_to_value[GevParams.SHAPE] == 0, 'The formula does not work if is not zero'
+ gev_params = self.get_unconstrained_gev_params(year)
+ gev_params.location = param_name_to_value[GevParams.LOC]
+ gev_params.scale = param_name_to_value[GevParams.SCALE]
+ return gev_params.mean
diff --git a/projects/altitude_spatial_model/altitudes_fit/altitudes_studies.py b/projects/altitude_spatial_model/altitudes_fit/altitudes_studies.py
index 284a7622..34a8b7ff 100644
--- a/projects/altitude_spatial_model/altitudes_fit/altitudes_studies.py
+++ b/projects/altitude_spatial_model/altitudes_fit/altitudes_studies.py
@@ -7,6 +7,7 @@ from cached_property import cached_property
from extreme_data.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
from extreme_data.meteo_france_data.scm_models_data.visualization.main_study_visualizer import \
SCM_STUDY_CLASS_TO_ABBREVIATION
+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 spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
from spatio_temporal_dataset.coordinates.spatial_coordinates.abstract_spatial_coordinates import \
@@ -18,7 +19,6 @@ from spatio_temporal_dataset.coordinates.temporal_coordinates.generated_temporal
from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
from spatio_temporal_dataset.spatio_temporal_observations.annual_maxima_observations import AnnualMaxima
import matplotlib.pyplot as plt
-import matplotlib.colors as mcolors
class AltitudesStudies(object):
@@ -157,18 +157,4 @@ class AltitudesStudies(object):
moment = function(annual_maxima)
mean_moment.append(moment)
altitudes.append(altitude)
- self.plot_against_altitude(altitudes, ax, massif_id, massif_name, mean_moment)
-
- @staticmethod
- def plot_against_altitude(altitudes, ax, massif_id, massif_name, mean_moment):
- di = massif_id // 8
- if di == 0:
- linestyle = '-'
- elif di == 1:
- linestyle = 'dotted'
- else:
- linestyle = '--'
- colors = list(mcolors.TABLEAU_COLORS)
- colors[-3:-1] = [] # remove gray and olive
- color = colors[massif_id % 8]
- ax.plot(altitudes, mean_moment, color=color, linewidth=2, label=massif_name, linestyle=linestyle)
+ plot_against_altitude(altitudes, ax, massif_id, massif_name, mean_moment)
diff --git a/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/main_snowfall_article.py b/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/main_snowfall_article.py
index a2b9366d..cbc4c2af 100644
--- a/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/main_snowfall_article.py
+++ b/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/main_snowfall_article.py
@@ -4,6 +4,7 @@ from multiprocessing.pool import Pool
import matplotlib as mpl
from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day
+from projects.contrasting_trends_in_snow_loads.snowfall_versus_time_and_altitude.snowfall_plot import plot_snowfall_mean
from projects.contrasting_trends_in_snow_loads.snowfall_versus_time_and_altitude.study_visualizer_for_mean_values import \
StudyVisualizerForMeanValues
from projects.contrasting_trends_in_snow_loads.snowfall_versus_time_and_altitude.validation_plot import validation_plot
@@ -70,16 +71,18 @@ def intermediate_result(altitudes, massif_names=None,
_ = compute_minimized_aic(visualizer)
# Plots
- validation_plot(altitude_to_visualizer)
+ # validation_plot(altitude_to_visualizer)
+ plot_snowfall_mean(altitude_to_visualizer)
def major_result():
uncertainty_methods = [ConfidenceIntervalMethodFromExtremes.ci_mle][:]
- massif_names = ['Beaufortain', 'Vercors']
+ # massif_names = ['Beaufortain', 'Vercors']
massif_names = None
study_classes = [SafranSnowfall1Day]
model_subsets_for_uncertainty = None
altitudes = paper_altitudes
+ altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000]
# altitudes = [900, 1200]
for study_class in study_classes:
diff --git a/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/snowfall_plot.py b/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/snowfall_plot.py
new file mode 100644
index 00000000..5d17a958
--- /dev/null
+++ b/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/snowfall_plot.py
@@ -0,0 +1,72 @@
+from typing import Dict
+import matplotlib.pyplot as plt
+
+import numpy as np
+from sklearn.linear_model import LinearRegression
+
+from extreme_data.meteo_france_data.scm_models_data.visualization.main_study_visualizer import \
+ SCM_STUDY_CLASS_TO_ABBREVIATION
+from extreme_data.meteo_france_data.scm_models_data.visualization.plot_utils import plot_against_altitude
+from projects.contrasting_trends_in_snow_loads.snowfall_versus_time_and_altitude.study_visualizer_for_mean_values import \
+ StudyVisualizerForMeanValues
+
+
+def fit_linear_regression(x, y):
+ X = np.array(x).reshape(-1, 1)
+ reg = LinearRegression().fit(X, y)
+ r2_score = reg.score(X, y)
+ a = reg.coef_
+ b = reg.intercept_
+ return a, b, r2_score
+
+
+def plot_snowfall_mean(altitude_to_visualizer: Dict[int, StudyVisualizerForMeanValues]):
+ visualizer = list(altitude_to_visualizer.values())[0]
+ # Plot the curve for the evolution of the mean
+ massif_name_to_a, massif_name_to_b, massif_name_to_r2_score = plot_mean(altitude_to_visualizer, derivative=False)
+ # Plot map with the coefficient a
+ visualizer.plot_abstract_fast(massif_name_to_a, label='a')
+ visualizer.plot_abstract_fast(massif_name_to_b, label='b')
+ visualizer.plot_abstract_fast(massif_name_to_r2_score, label='r2')
+
+
+def plot_mean(altitude_to_visualizer: Dict[int, StudyVisualizerForMeanValues], derivative=False):
+ massif_name_to_linear_regression_result = {}
+
+ altitudes = list(altitude_to_visualizer.keys())
+ visualizers = list(altitude_to_visualizer.values())
+ visualizer = visualizers[0]
+ study = visualizer.study
+ year = study.year_max
+
+ for massif_id, massif_name in enumerate(visualizer.study.all_massif_names()):
+ altitudes_massif = [a for a, v in altitude_to_visualizer.items()
+ if massif_name in v.massif_name_to_trend_test_that_minimized_aic]
+ if len(altitudes_massif) >= 2:
+ trend_tests = [altitude_to_visualizer[a].massif_name_to_trend_test_that_minimized_aic[massif_name]
+ for a in altitudes_massif]
+ if derivative:
+ moment = 'time derivative of the mean'
+ values = [t.first_derivative_mean_value(year=year) for t in trend_tests]
+ else:
+ moment = 'mean'
+ values = [t.mean_value(year=year) for t in trend_tests]
+ massif_name_to_linear_regression_result[massif_name] = fit_linear_regression(altitudes, values)
+ plot_values_against_altitudes(altitudes_massif, massif_id, massif_name, moment, study, values, visualizer)
+
+ return [{m: t[i][0] if i == 0 else t[i]
+ for m, t in massif_name_to_linear_regression_result.items()} for i in range(3)]
+
+
+def plot_values_against_altitudes(altitudes, massif_id, massif_name, moment, study, values, visualizer):
+ ax = plt.gca()
+ plot_against_altitude(altitudes=altitudes, ax=ax, massif_id=massif_id, massif_name=massif_name, values=values)
+ plot_name = '{} annual maxima of {}'.format(moment, SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)])
+ ax.set_ylabel('{} ({})'.format(plot_name, study.variable_unit), fontsize=15)
+ ax.set_xlabel('altitudes', fontsize=15)
+ lim_down, lim_up = ax.get_ylim()
+ lim_up += (lim_up - lim_down) / 3
+ ax.set_ylim([lim_down, lim_up])
+ ax.tick_params(axis='both', which='major', labelsize=13)
+ visualizer.plot_name = plot_name
+ visualizer.show_or_save_to_file(dpi=500, add_classic_title=False)
diff --git a/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/study_visualizer_for_mean_values.py b/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/study_visualizer_for_mean_values.py
index 84dd12be..334ec2ca 100644
--- a/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/study_visualizer_for_mean_values.py
+++ b/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/study_visualizer_for_mean_values.py
@@ -39,7 +39,6 @@ class StudyVisualizerForMeanValues(StudyVisualizerForNonStationaryTrends):
massif_name_to_parameter_value = {}
for massif_name, trend_test in self.massif_name_to_trend_test_that_minimized_aic.items():
parameter_value = trend_test.unconstrained_average_mean_value
- print(type(parameter_value), parameter_value)
massif_name_to_parameter_value[massif_name] = parameter_value
return massif_name_to_parameter_value
@@ -56,4 +55,3 @@ class StudyVisualizerForMeanValues(StudyVisualizerForNonStationaryTrends):
def plot_abstract_fast(self, massif_name_to_value, label, graduation=10.0, cmap=plt.cm.coolwarm):
super().plot_abstract(massif_name_to_value, label, label, self.fit_method, graduation, cmap)
-
diff --git a/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/validation_plot.py b/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/validation_plot.py
index 58d27a21..6b630252 100644
--- a/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/validation_plot.py
+++ b/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/validation_plot.py
@@ -32,7 +32,6 @@ def plot_shoe_relative_differences_distribution(altitude_to_relative_differences
ax.set_xlabel('Altitude (m)')
ax.legend()
ax.grid()
- ax.set_ylim([-8, 5])
def plot_relative_difference_map(visualizer: StudyVisualizerForMeanValues):
@@ -40,10 +39,8 @@ def plot_relative_difference_map(visualizer: StudyVisualizerForMeanValues):
visualizer.plot_abstract_fast(massif_name_to_value=visualizer.massif_name_to_empirical_mean,
label='Empirical' + label)
- print(visualizer.massif_name_to_parametric_mean)
visualizer.plot_abstract_fast(massif_name_to_value=visualizer.massif_name_to_parametric_mean,
label='Model' + label)
- print(visualizer.massif_name_to_relative_difference_for_mean)
visualizer.plot_abstract_fast(massif_name_to_value=visualizer.massif_name_to_relative_difference_for_mean,
- label='Relative difference of' + label)
+ label='Relative difference of' + label, graduation=1)
return list(visualizer.massif_name_to_relative_difference_for_mean.values())
--
GitLab