From b0131941c7bd1ec5c53d26d6adc3a8ff06af51de Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Wed, 29 Apr 2020 19:30:58 +0200
Subject: [PATCH] [contrasting] add time derivative for the linear models. add
one_axis_param_function.py
---
.../visualization/plot_utils.py | 3 +-
.../param_function/one_axis_param_function.py | 39 +++++++++++++++++++
.../function/param_function/param_function.py | 25 +++++-------
.../main_snowfall_article.py | 8 ++--
.../snowfall_plot.py | 28 ++++++++++++-
.../validation_plot.py | 1 +
6 files changed, 84 insertions(+), 20 deletions(-)
create mode 100644 extreme_fit/function/param_function/one_axis_param_function.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 b029c576..70e8fa68 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
@@ -17,7 +17,8 @@ def plot_against_altitude(altitudes, ax, massif_id, massif_name, values):
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)
+ massif_name_str = ' '.join(massif_name.split('_'))
+ ax.plot(altitudes, values, color=color, linewidth=2, label=massif_name_str, linestyle=linestyle)
def load_plot(cmap, graduation, label, massif_name_to_value, altitude, fit_method, add_x_label=True):
diff --git a/extreme_fit/function/param_function/one_axis_param_function.py b/extreme_fit/function/param_function/one_axis_param_function.py
new file mode 100644
index 00000000..c7292fee
--- /dev/null
+++ b/extreme_fit/function/param_function/one_axis_param_function.py
@@ -0,0 +1,39 @@
+import numpy as np
+
+
+class AbstractOneAxisParamFunction(object):
+
+ def __init__(self, dim: int, coordinates: np.ndarray):
+ self.dim = dim
+ self.t_min = coordinates[:, dim].min()
+ self.t_max = coordinates[:, dim].max()
+
+ def get_param_value(self, coordinate: np.ndarray, assert_range=True) -> float:
+ t = coordinate[self.dim]
+ if assert_range:
+ assert self.t_min <= t <= self.t_max, '{} is out of bounds ({}, {})'.format(t, self.t_min, self.t_max)
+ return self._get_param_value(t)
+
+ def _get_param_value(self, t) -> float:
+ raise NotImplementedError
+
+
+class LinearOneAxisParamFunction(AbstractOneAxisParamFunction):
+
+ def __init__(self, dim: int, coordinates: np.ndarray, coef: float = 0.01):
+ super().__init__(dim, coordinates)
+ self.coef = coef
+
+ def _get_param_value(self, t) -> float:
+ return t * self.coef
+
+
+class QuadraticOneAxisParamFunction(AbstractOneAxisParamFunction):
+
+ def __init__(self, dim: int, coordinates: np.ndarray, coef1: float = 0.01, coef2: float = 0.01):
+ super().__init__(dim, coordinates)
+ self.coef1 = coef1
+ self.coef2 = coef2
+
+ def _get_param_value(self, t) -> float:
+ return np.power(t, 2) * self.coef2 + t * self.coef1
diff --git a/extreme_fit/function/param_function/param_function.py b/extreme_fit/function/param_function/param_function.py
index 575ff747..1211dd9e 100644
--- a/extreme_fit/function/param_function/param_function.py
+++ b/extreme_fit/function/param_function/param_function.py
@@ -1,6 +1,7 @@
from typing import List
import numpy as np
from extreme_fit.function.param_function.linear_coef import LinearCoef
+from extreme_fit.function.param_function.one_axis_param_function import LinearOneAxisParamFunction
from extreme_fit.function.param_function.spline_coef import SplineCoef
@@ -22,20 +23,8 @@ class ConstantParamFunction(AbstractParamFunction):
def get_param_value(self, coordinate: np.ndarray) -> float:
return self.constant
-
-class LinearOneAxisParamFunction(AbstractParamFunction):
-
- def __init__(self, dim: int, coordinates: np.ndarray, coef: float = 0.01):
- self.dim = dim
- self.t_min = coordinates[:, dim].min()
- self.t_max = coordinates[:, dim].max()
- self.coef = coef
-
- def get_param_value(self, coordinate: np.ndarray) -> float:
- t = coordinate[self.dim]
- if self.OUT_OF_BOUNDS_ASSERT:
- assert self.t_min <= t <= self.t_max, '{} is out of bounds ({}, {})'.format(t, self.t_min, self.t_max)
- return t * self.coef
+ def get_first_derivative_param_value(self, coordinate: np.ndarray, dim: int) -> float:
+ return 0
class LinearParamFunction(AbstractParamFunction):
@@ -48,14 +37,20 @@ class LinearParamFunction(AbstractParamFunction):
param_function = LinearOneAxisParamFunction(dim=dim, coordinates=coordinates,
coef=self.linear_coef.get_coef(idx=dim))
self.linear_one_axis_param_functions.append(param_function)
+ self.dim_to_linear_one_axis_param_function = dict(zip(dims, self.linear_one_axis_param_functions))
def get_param_value(self, coordinate: np.ndarray) -> float:
# Add the intercept and the value with respect to each axis
gev_param_value = self.linear_coef.intercept
for linear_one_axis_param_function in self.linear_one_axis_param_functions:
- gev_param_value += linear_one_axis_param_function.get_param_value(coordinate)
+ gev_param_value += linear_one_axis_param_function.get_param_value(coordinate, self.OUT_OF_BOUNDS_ASSERT)
return gev_param_value
+ def get_first_derivative_param_value(self, coordinate: np.ndarray, dim: int) -> float:
+ return self.dim_to_linear_one_axis_param_function[dim].coef
+
+
+
class SplineParamFunction(AbstractParamFunction):
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 33c73e08..3a6a9fec 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,7 +4,8 @@ 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.snowfall_plot import \
+ plot_snowfall_mean, plot_snowfall_time_derivative_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
@@ -71,8 +72,9 @@ def intermediate_result(altitudes, massif_names=None,
_ = compute_minimized_aic(visualizer)
# Plots
- validation_plot(altitude_to_visualizer)
- plot_snowfall_mean(altitude_to_visualizer)
+ # validation_plot(altitude_to_visualizer)
+ # plot_snowfall_mean(altitude_to_visualizer)
+ plot_snowfall_time_derivative_mean(altitude_to_visualizer)
def major_result():
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
index dfbbf81e..a51a4701 100644
--- 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
@@ -20,6 +20,30 @@ def fit_linear_regression(x, y):
return a, b, r2_score
+def plot_snowfall_time_derivative_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=True)
+ # Augmentation every km
+ massif_name_to_augmentation_every_km = {m: a * 1000 for m, a in massif_name_to_a.items()}
+ visualizer.plot_abstract_fast(massif_name_to_augmentation_every_km,
+ label='Augmentation of time derivative of mean annual maxima for every km of elevation ()',
+ add_x_label=False)
+ # Value at 2000 m
+ massif_name_to_mean_at_2000 = {m: a * 2000 + massif_name_to_b[m] for m, a in massif_name_to_a.items()}
+ visualizer.plot_abstract_fast(massif_name_to_mean_at_2000, label='Time derivative of mean annual maxima of at 2000 m ()',
+ add_x_label=False)
+ # Altitude for the change of dynamic
+ massif_name_to_altitude_change_dynamic = {m: - massif_name_to_b[m] / a for m, a in massif_name_to_a.items()}
+ # Keep only those that are in a reasonable range
+ massif_name_to_altitude_change_dynamic = {m: d for m, d in massif_name_to_altitude_change_dynamic.items()
+ if 0 < d < 3000}
+ visualizer.plot_abstract_fast(massif_name_to_altitude_change_dynamic, label='Altitude for the change of dynamic ()',
+ add_x_label=False, graduation=500)
+ # R2 score
+ visualizer.plot_abstract_fast(massif_name_to_r2_score, label='r2', graduation=0.1, add_x_label=False)
+
+
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
@@ -60,6 +84,8 @@ def plot_mean(altitude_to_visualizer: Dict[int, StudyVisualizerForMeanValues], d
values = [t.mean_value(year=year) for t in trend_tests]
massif_name_to_linear_regression_result[massif_name] = fit_linear_regression(altitudes_massif, values)
plot_values_against_altitudes(ax, altitudes_massif, massif_id, massif_name, moment, study, values, visualizer)
+ ax.legend(prop={'size': 7}, ncol=3)
+ visualizer.show_or_save_to_file(dpi=500, add_classic_title=False)
plt.close()
return [{m: t[i][0] if i == 0 else t[i]
@@ -76,4 +102,4 @@ def plot_values_against_altitudes(ax, altitudes, massif_id, massif_name, moment,
# 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/validation_plot.py b/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/validation_plot.py
index 88a5a5b5..5fa8ee8b 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
@@ -20,6 +20,7 @@ def validation_plot(altitude_to_visualizer: Dict[int, StudyVisualizerForMeanValu
# Shoe plot with respect to the altitude.
plot_shoe_relative_differences_distribution(altitude_to_relative_differences, altitudes)
study_visualizer.show_or_save_to_file(add_classic_title=False, dpi=500)
+ plt.close()
def plot_shoe_relative_differences_distribution(altitude_to_relative_differences, altitudes):
--
GitLab