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 b029c576a3653017a0e35e7c8c77d2b5c2f851cd..70e8fa68e568b46714a3cac2d774a8376b2256e6 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 0000000000000000000000000000000000000000..c7292feef4de711c4e7892344b06363766fbb807 --- /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 575ff7477925dd7272c3ad2789e70d985483c05a..1211dd9edf91dd4eb86aa8b871b9853590445ed2 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 33c73e08ecd5391ef601eb1d3d51482280ac8ebc..3a6a9feccf6cc6f2999aa0f1eea41b4b2924c21e 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 dfbbf81ed3123f5a5b5581050027cd1bfc61a4ef..a51a470135b80e19ccc934a41dd72f384908c3ce 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 88a5a5b5f178b4acc503226f8eedc0eff0c46c4f..5fa8ee8bb790eb7846073a95ed0af46c72e40458 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):