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 3bb1802926afbe8f88a84aeabbe4b475c5bbcf4c..82b712a59427463f9d3a824e220149bb75f56c24 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
@@ -35,7 +35,8 @@ def plot_against_altitude(x_ticks, ax, massif_id, massif_name, values, altitude=
def load_plot(cmap, graduation, label, massif_name_to_value, altitude, fit_method, add_x_label=True,
- negative_and_positive_values=True, massif_name_to_text=None, add_colorbar=True, max_abs_change=None):
+ negative_and_positive_values=True, massif_name_to_text=None, add_colorbar=True, max_abs_change=None,
+ xlabel=None):
if max_abs_change is None:
max_abs_change = max([abs(e) for e in massif_name_to_value.values()])
if negative_and_positive_values:
@@ -82,5 +83,8 @@ def load_plot(cmap, graduation, label, massif_name_to_value, altitude, fit_metho
ax.get_xaxis().set_visible(True)
ax.set_xticks([])
if add_x_label:
- ax.set_xlabel('Altitude = {}m'.format(altitude), fontsize=15)
+ if xlabel is None:
+ ax.set_xlabel('Altitude = {}m'.format(altitude), fontsize=15)
+ else:
+ ax.set_xlabel(xlabel, fontsize=10)
# 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 2ca68c259b72296f721aec607d05869b16a1dddc..f8d3172c5b263256c673dd72ee4ccab27b7fc1da 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
@@ -717,14 +717,14 @@ class StudyVisualizer(VisualizationParameters):
def plot_map(self, cmap, fit_method, graduation, label, massif_name_to_value, plot_name, add_x_label=True,
negative_and_positive_values=True, massif_name_to_text=None, altitude=None, add_colorbar=True,
- max_abs_change=None):
+ max_abs_change=None, xlabel=None):
if altitude is None:
altitude = self.study.altitude
if len(massif_name_to_value) > 0:
load_plot(cmap, graduation, label, massif_name_to_value, altitude, fitmethod_to_str(fit_method),
add_x_label=add_x_label, negative_and_positive_values=negative_and_positive_values,
massif_name_to_text=massif_name_to_text,
- add_colorbar=add_colorbar, max_abs_change=max_abs_change)
+ add_colorbar=add_colorbar, max_abs_change=max_abs_change, xlabel=xlabel)
self.plot_name = plot_name
# self.show_or_save_to_file(add_classic_title=False, tight_layout=True, no_title=True, dpi=500)
self.show_or_save_to_file(add_classic_title=False, no_title=True, dpi=500, tight_layout=True)
diff --git a/extreme_fit/distribution/gev/ci_fevd_fixed.R b/extreme_fit/distribution/gev/ci_fevd_fixed.R
index 785b3ca3d8a7c347581a99aa853ca4eb2f1dd959..2dfb7b60d4b91d71af9a385217789f3268901c1f 100644
--- a/extreme_fit/distribution/gev/ci_fevd_fixed.R
+++ b/extreme_fit/distribution/gev/ci_fevd_fixed.R
@@ -424,6 +424,7 @@ ci.rl.ns.fevd.mle_fixed <- function (x, alpha = 0.05, return.period = 100, metho
qcov.base = qcov.base)
z.alpha <- qnorm(alpha/2, lower.tail = FALSE)
cov.theta <- parcov.fevd(x)
+ # cov.theta <- parcov.fevd_fixed(x)
if (is.null(cov.theta))
stop("ci: Sorry, unable to calculate the parameter covariance matrix. Maybe try a different method.")
var.theta <- diag(cov.theta)
@@ -454,6 +455,60 @@ ci.rl.ns.fevd.mle_fixed <- function (x, alpha = 0.05, return.period = 100, metho
return(out)
}
+parcov.fevd_fixed <- function (x)
+{
+ cov.theta <- NULL
+ theta.hat <- x$results$par
+ theta.names <- names(theta.hat)
+ if (is.element("log.scale", theta.names)) {
+ theta.hat[theta.names == "log.scale"] <- exp(theta.hat[theta.names ==
+ "log.scale"])
+ theta.names[theta.names == "log.scale"] <- "scale"
+ names(theta.hat) <- theta.names
+ phiU <- FALSE
+ }
+ else phiU <- x$par.models$log.scale
+ y <- datagrabber(x)
+ if (x$data.name[2] != "") {
+ xdat <- y[, 1]
+ data <- y[, -1]
+ }
+ else {
+ xdat <- y[, 1]
+ data <- NULL
+ }
+ designs <- setup.design(x)
+ if (x$method != "GMLE") {
+ hold <- try(suppressWarnings(optimHess(theta.hat, oevd,
+ gr = grlevd, o = x, des = designs, x = xdat, data = data,
+ u = x$threshold, npy = x$npy, phi = phiU, blocks = x$blocks)),
+ silent = TRUE)
+ if ((class(hold) != "try-error") && all(!is.na(hold))) {
+ cov.theta <- try(suppressWarnings(solve(hold)), silent = TRUE)
+ print(cov.theta)
+ # if (any(diag(cov.theta) <= 0))
+ # re.do <- TRUE
+ # else re.do <- FALSE
+ re.do = TRUE
+ }
+ else re.do <- TRUE
+ }
+ else re.do <- TRUE
+ if (re.do) {
+ hold <- try(optimHess(theta.hat, oevd, o = x, des = designs,
+ x = xdat, data = data, u = x$threshold, npy = x$npy,
+ phi = phiU, blocks = x$blocks), silent = TRUE)
+ if (class(hold) != "try-error" && all(!is.na(hold))) {
+ cov.theta <- try(solve(hold), silent = TRUE)
+ if (class(cov.theta) == "try-error" || any(diag(cov.theta) <=
+ 0))
+ cov.theta <- NULL
+ }
+ }
+ return(cov.theta)
+}
+
+
return.level.ns.fevd.mle_fixed <- function (x, return.period = c(2, 20, 100), ..., alpha = 0.05,
method = c("normal"), do.ci = FALSE, verbose = FALSE, qcov = NULL,
qcov.base = NULL)
diff --git a/extreme_fit/distribution/gev/main_fevd_mle_two_covariates.R b/extreme_fit/distribution/gev/main_fevd_mle_two_covariates.R
index 49fcd28b201945faff7fe46f6b2751f99ffb0525..02360381c4aaae915e7e1589e7697f0133f53879 100644
--- a/extreme_fit/distribution/gev/main_fevd_mle_two_covariates.R
+++ b/extreme_fit/distribution/gev/main_fevd_mle_two_covariates.R
@@ -11,8 +11,8 @@ source('ci_fevd_fixed.R')
source('summary_fevd_fixed.R')
# Sample from a GEV
set.seed(42)
-N <- 50
-loc = 0; scale = 1; shape <- 1
+N <- 100
+loc = 0; scale = 1; shape <- 0.1
x_gev <- rgev(N, loc = loc, scale = scale, shape = shape)
# start_loc = 0; start_scale = 1; start_shape = 1
# N <- 50
@@ -22,14 +22,14 @@ print(N)
coord <- matrix(ncol=2, nrow = N)
coord[,1]=seq(0,N-1,1)
coord[,2]=seq(0,N-1,1)
-print(coord)
+
colnames(coord) = c("X", "T")
coord = data.frame(coord, stringsAsFactors = TRUE)
# res = fevd_fixed(x_gev, data=coord, method='MLE', verbose=TRUE, use.phi=FALSE)
# res = fevd_fixed(x_gev, data=coord, location.fun= ~T, scale.fun= ~T, method='MLE', type="GEV", verbose=FALSE, use.phi=FALSE)
# res = fevd_fixed(x_gev, data=coord, location.fun= ~sin(X) + cos(T), method='MLE', type="GEV", verbose=FALSE, use.phi=FALSE)
# res = fevd_fixed(x_gev, data=coord, location.fun= ~poly(X * T, 1, raw = TRUE), method='MLE', type="Gumbel", verbose=FALSE, use.phi=FALSE)
-res = fevd_fixed(x_gev, data=coord, location.fun= ~poly(X, 1, raw = TRUE) + poly(T, 2, raw = TRUE) , method='MLE', type="Gumbel", verbose=FALSE, use.phi=FALSE)
+res = fevd_fixed(x_gev, data=coord, location.fun= ~poly(X, 1, raw = TRUE) + poly(T, 1, raw = TRUE) , method='MLE', type="Gumbel", verbose=FALSE, use.phi=FALSE)
# print(res)
print(summary.fevd.mle_fixed(res))
# print(summary(res)$AIC)
@@ -44,6 +44,18 @@ print(summary.fevd.mle_fixed(res))
# print(m[1])
+# v = make.qcov(res, vals = list(mu0 = c(0.0), mu1 = c(0.0), mu2 = c(0.0), sigma = c(0.0)))
+v = make.qcov(res, vals = list(mu1 = c(0.0), mu2 = c(0.0)))
+
+res_ci = ci.fevd.mle_fixed(res, alpha = 0.05, type = c("return.level"),
+ return.period = 50, method = "normal", xrange = NULL, nint = 20, verbose = FALSE,
+ tscale = FALSE, return.samples = FALSE, qcov=v)
+print(res_ci)
+
+
+
+
+
# Confidence interval staionary
# method = "proflik"
# res_ci = ci.fevd.mle(res, alpha = 0.05, type = c("return.level"),
diff --git a/projects/altitude_spatial_model/altitudes_fit/altitudes_studies.py b/projects/altitude_spatial_model/altitudes_fit/altitudes_studies.py
index 6221655315816d1c053a51adf8f16a65839d174c..b72a90d26253604733128c4d075c3c5315c88704 100644
--- a/projects/altitude_spatial_model/altitudes_fit/altitudes_studies.py
+++ b/projects/altitude_spatial_model/altitudes_fit/altitudes_studies.py
@@ -43,9 +43,11 @@ class AltitudesStudies(object):
# Dataset Loader
def spatio_temporal_dataset(self, massif_name, s_split_spatial: pd.Series = None,
- s_split_temporal: pd.Series = None):
+ s_split_temporal: pd.Series = None,
+ massif_altitudes=None):
coordinate_values_to_maxima = {}
- massif_altitudes = self.massif_name_to_altitudes[massif_name]
+ if massif_altitudes is None:
+ massif_altitudes = self.massif_name_to_altitudes[massif_name]
if len(massif_altitudes) == 0:
print('{} has no data for these altitudes: {}'.format(massif_name, self.altitudes))
for altitude in massif_altitudes:
diff --git a/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py b/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
index c99cd54855bc127b22326b2e8e8a1b5f2569efcc..650ed6d86b0960fcb43a9a1998986ea82a818e1b 100644
--- a/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
+++ b/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
@@ -9,7 +9,6 @@ from projects.altitude_spatial_model.altitudes_fit.plot_histogram_altitude_studi
mpl.rcParams['text.usetex'] = True
mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
-
from projects.altitude_spatial_model.altitudes_fit.utils_altitude_studies_visualizer import load_visualizer_list
from extreme_fit.model.margin_model.polynomial_margin_model.utils import \
@@ -29,13 +28,20 @@ from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitudes_s
def main():
study_classes = [SafranSnowfall1Day, SafranSnowfall3Days, SafranSnowfall5Days, SafranSnowfall7Days][:1]
# study_classes = [SafranPrecipitation1Day][:1]
-
- massif_names = None
- # massif_names = ['Mercantour', 'Vercors', 'Ubaye']
-
seasons = [Season.annual, Season.winter, Season.spring, Season.automn][:1]
- main_loop(altitudes_for_groups[:], massif_names, seasons, study_classes)
+ fast = None
+ if fast is None:
+ massif_names = None
+ altitudes_list = altitudes_for_groups[:2]
+ elif fast:
+ massif_names = ['Mercantour', 'Vercors', 'Ubaye']
+ altitudes_list = altitudes_for_groups[:2]
+ else:
+ massif_names = None
+ altitudes_list = altitudes_for_groups
+
+ main_loop(altitudes_list, massif_names, seasons, study_classes)
def main_loop(altitudes_list, massif_names, seasons, study_classes):
@@ -59,7 +65,7 @@ def plot_visualizers(massif_names, visualizer_list):
def plot_visualizer(massif_names, visualizer):
# Plot time series
- # visualizer.studies.plot_maxima_time_series(massif_names=massif_names)
+ visualizer.studies.plot_maxima_time_series(massif_names=massif_names)
# Plot moments against altitude
# for std in [True, False][:]:
# for change in [True, False, None]:
diff --git a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitude_group.py b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitude_group.py
index 786a8684223d0bf7ec0c78c69edc405dff4508e5..5bf1bc5cfa9c51c6aa766199a149380438663576 100644
--- a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitude_group.py
+++ b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitude_group.py
@@ -2,26 +2,13 @@ from enum import Enum
# The order is important
altitudes_for_groups = [
- [300, 600, 900],
+ [300, 600, 900][1:],
[1200, 1500, 1800],
[2100, 2400, 2700],
[3000, 3300, 3600, 3900]
]
-# altitudes_for_groups = [
-# [900, 1200, 1500],
-# [1800, 2100, 2400],
-# [2700, 3000, 3300]
-# ]
-
-# altitudes_for_groups = [
-# [600, 900, 1200, 1500, 1800],
-# [1500, 1800, 2100, 2400, 2700],
-# [2400, 2700, 3000, 3300, 3600]
-# ]
-
-
class AbstractAltitudeGroup(object):
@property
@@ -32,6 +19,18 @@ class AbstractAltitudeGroup(object):
def reference_altitude(self):
raise NotImplementedError
+ @property
+ def xlabel(self):
+ # warning: the label could not correspond to all massifs, some might have been fitted with less data
+ # idx = get_index_group_from_reference_altitude(reference_altitude)
+ # min_altitude, *_, max_altitude = altitudes_for_groups[idx]
+ i = self.reference_altitude // 1000
+ min_altitude, max_altitude = 1000 * i, 1000 * (i + 1)
+ return 'Altitude = {} m\n' \
+ 'Estimated with maxima between {} m and {} m'.format(self.reference_altitude,
+ min_altitude,
+ max_altitude)
+
class LowAltitudeGroup(AbstractAltitudeGroup):
@@ -87,7 +86,6 @@ class DefaultAltitudeGroup(AbstractAltitudeGroup):
def reference_altitude(self):
return 500
-
def get_altitude_group_from_altitudes(altitudes):
s = set(altitudes)
if s == set(altitudes_for_groups[0]):
diff --git a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
index ae27f5c5cb09d924143ed5429531229270ed6b94..74ad0ebbe4163ebc33a9821925f7b3e3f38a9e46 100644
--- a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
+++ b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
@@ -51,12 +51,12 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
# Load one fold fit
self._massif_name_to_one_fold_fit = {}
for massif_name in self.massif_names:
- if self.load_condition(massif_name):
+ # Load valid massif altitudes
+ massif_altitudes = self.get_massif_altitudes(massif_name)
+ if self.load_condition(massif_altitudes):
# Load dataset
- dataset = studies.spatio_temporal_dataset(massif_name=massif_name)
- dataset_without_zero = AbstractDataset.remove_zeros(dataset.observations,
- dataset.coordinates)
- old_fold_fit = OneFoldFit(massif_name, dataset_without_zero, model_classes, self.fit_method,
+ dataset = studies.spatio_temporal_dataset(massif_name=massif_name, massif_altitudes=massif_altitudes)
+ old_fold_fit = OneFoldFit(massif_name, dataset, model_classes, self.fit_method,
self.temporal_covariate_for_fit,
type(self.altitude_group),
self.display_only_model_that_pass_anderson_test)
@@ -69,10 +69,25 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
moment_names = ['moment', 'changes_for_moment', 'relative_changes_for_moment'][:]
orders = [1, 2, None][2:]
- def load_condition(self, massif_name):
- altitudes_for_massif = [altitude for altitude, study in self.studies.altitude_to_study.items()
+ def get_massif_altitudes(self, massif_name):
+ valid_altitudes = [altitude for altitude, study in self.studies.altitude_to_study.items()
if massif_name in study.study_massif_names]
- return (self.altitude_group.reference_altitude in altitudes_for_massif) and (len(altitudes_for_massif) >= 2)
+ massif_altitudes = []
+ for altitude in valid_altitudes:
+ study = self.studies.altitude_to_study[altitude]
+ annual_maxima = study.massif_name_to_annual_maxima[massif_name]
+ percentage_of_non_zeros = 100 * np.count_nonzero(annual_maxima) / len(annual_maxima)
+ if percentage_of_non_zeros > 90:
+ massif_altitudes.append(altitude)
+ else:
+ print(massif_name, altitude, percentage_of_non_zeros)
+ return massif_altitudes
+
+ def load_condition(self, massif_altitudes):
+ # At least two altitudes for the estimated, including the reference altitude
+ reference_altitude_is_in_altitudes = (self.altitude_group.reference_altitude in massif_altitudes)
+ at_least_two_altitudes = (len(massif_altitudes) >= 2)
+ return reference_altitude_is_in_altitudes and at_least_two_altitudes
@property
def massif_name_to_one_fold_fit(self) -> Dict[str, OneFoldFit]:
@@ -150,6 +165,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
add_colorbar=add_colorbar,
max_abs_change=max_abs_change,
massif_name_to_text=massif_name_to_text,
+ xlabel=self.altitude_group.xlabel,
)
@property
@@ -250,6 +266,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
add_colorbar=self.add_colorbar,
max_abs_change=self._max_abs_for_shape,
fit_method=self.fit_method,
+ xlabel=self.altitude_group.xlabel,
)
def plot_altitude_for_the_peak(self):
@@ -374,24 +391,25 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
nb_valid_massif_names = len(valid_massif_names)
nbs = np.zeros(4)
- relative_changes = [[], []]
+ relative_changes = []
for one_fold in [one_fold for m, one_fold in self.massif_name_to_one_fold_fit.items()
if m in valid_massif_names]:
- if one_fold.change_for_reference_altitude == 0:
- continue
# Compute relative changes
- relative_changes[0].append(one_fold.relative_change_for_reference_altitude)
- if one_fold.is_significant:
- relative_changes[1].append(one_fold.relative_change_for_reference_altitude)
+ relative_changes.append(one_fold.relative_change_in_return_level_for_reference_altitude)
+ # Compute nb of non stationary models
+ if one_fold.change_in_return_level_for_reference_altitude == 0:
+ continue
# Compute nbs
- idx = 0 if one_fold.change_for_reference_altitude < 0 else 2
+ idx = 0 if one_fold.change_in_return_level_for_reference_altitude < 0 else 2
nbs[idx] += 1
if one_fold.is_significant:
nbs[idx + 1] += 1
percents = 100 * nbs / nb_valid_massif_names
- mean_relative_changes = [np.mean(l) for l in relative_changes]
- return [nb_valid_massif_names] + mean_relative_changes + list(percents)
+ mean_relative_change = np.mean(relative_changes)
+ median_relative_change = np.median(relative_changes)
+ print('mean', mean_relative_change, relative_changes)
+ return [nb_valid_massif_names, mean_relative_change, median_relative_change] + list(percents)
def get_valid_names(self, massif_names):
valid_massif_names = set(self.massif_name_to_one_fold_fit.keys())
diff --git a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py
index badae11f99d6f44c2668d1142ab4065c5c65c897..f0e5cbf9c8978ccae210fc832a94bd40d9f67848 100644
--- a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py
+++ b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py
@@ -28,7 +28,7 @@ from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
class OneFoldFit(object):
SIGNIFICANCE_LEVEL = 0.05
best_estimator_minimizes_total_aic = False
- return_period = 50
+ return_period = 100
def __init__(self, massif_name: str, dataset: AbstractDataset, models_classes,
fit_method=MarginFitMethod.extremes_fevd_mle, temporal_covariate_for_fit=None,
@@ -89,12 +89,12 @@ class OneFoldFit(object):
return [self.get_moment(altitude, year, order) for altitude in altitudes]
@property
- def change_for_reference_altitude(self) -> float:
- return self.changes_for_moment(altitudes=[self.altitude_plot])[0]
+ def change_in_return_level_for_reference_altitude(self) -> float:
+ return self.changes_for_moment(altitudes=[self.altitude_plot], order=None)[0]
@property
- def relative_change_for_reference_altitude(self) -> float:
- return self.relative_changes_for_moment(altitudes=[self.altitude_plot])[0]
+ def relative_change_in_return_level_for_reference_altitude(self) -> float:
+ return self.relative_changes_for_moment(altitudes=[self.altitude_plot], order=None)[0]
def changes_for_moment(self, altitudes, year=2019, nb_years=50, order=1):
changes = []
diff --git a/projects/altitude_spatial_model/altitudes_fit/plot_histogram_altitude_studies.py b/projects/altitude_spatial_model/altitudes_fit/plot_histogram_altitude_studies.py
index 0fe005e6fde945d1bcd7f87d68f28d78c6a9f5df..9f44468897fe47cceaf3fe6a2a59558ed877c912 100644
--- a/projects/altitude_spatial_model/altitudes_fit/plot_histogram_altitude_studies.py
+++ b/projects/altitude_spatial_model/altitudes_fit/plot_histogram_altitude_studies.py
@@ -6,6 +6,7 @@ import matplotlib.pyplot as plt
from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitudes_studies_visualizer_for_non_stationary_models import \
AltitudesStudiesVisualizerForNonStationaryModels
+from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.one_fold_fit import OneFoldFit
def plot_histogram_all_models_against_altitudes(massif_names, visualizer_list: List[
@@ -32,8 +33,8 @@ def plot_histogram_all_models_against_altitudes(massif_names, visualizer_list: L
x_shifted = [tick_middle + width * shift / 2 for shift in range(-3, 5, 2)]
percentages = model_name_to_percentages[model_name]
colors = ['white', 'yellow', 'orange', 'red']
- labels = ['{} m (\% out of {} massifs)'.format(v.altitude_group.reference_altitude,
- len(v.get_valid_names(massif_names))) for v in visualizer_list]
+ labels = ['{} m - {} m (\% out of {} massifs)'.format(1000 * i, 1000 * (i+1),
+ len(v.get_valid_names(massif_names))) for i, v in enumerate(visualizer_list)]
for x, color, percentage, label in zip(x_shifted, colors, percentages, labels):
ax.bar([x], [percentage], width=width, label=label,
linewidth=linewidth, edgecolor='black', color=color)
@@ -53,9 +54,12 @@ def plot_histogram_all_models_against_altitudes(massif_names, visualizer_list: L
plt.close()
-def plot_histogram_all_trends_against_altitudes(massif_names, visualizer_list):
+def plot_histogram_all_trends_against_altitudes(massif_names, visualizer_list: List[
+ AltitudesStudiesVisualizerForNonStationaryModels]):
+ visualizer = visualizer_list[0]
+
all_trends = [v.all_trends(massif_names) for v in visualizer_list]
- nb_massifs, mean_relative_change, mean_significant_relative_change, *all_l = zip(*all_trends)
+ nb_massifs, mean_relative_changes, median_relative_changes, *all_l = zip(*all_trends)
plt.close()
ax = plt.gca()
@@ -87,17 +91,14 @@ def plot_histogram_all_trends_against_altitudes(massif_names, visualizer_list):
# PLot mean relative change against altitude
ax_twinx = ax.twinx()
- colors = ['grey', 'black']
- relative_changes = [mean_relative_change, mean_significant_relative_change]
- labels = ['All massifs with a trend', 'Only massifs with a significant trend']
- for color, y, label in zip(colors, relative_changes, labels):
- print(x, y)
- ax_twinx.plot(x, y, label=label, color=color)
+ ax_twinx.plot(x, mean_relative_changes, label='Mean relative change', color='black')
+ ax_twinx.plot(x, median_relative_changes, label='Median relative change', color='grey')
ax_twinx.legend(loc='upper right', prop={'size': size})
+ ylabel = 'Relative change of {}-year return levels ({})'.format(OneFoldFit.return_period,
+ visualizer.study.variable_unit)
+ ax_twinx.set_ylabel(ylabel, fontsize=legend_fontsize)
- ax_twinx.set_ylabel('Mean relative change average on massif with trends', fontsize=legend_fontsize)
- visualizer = visualizer_list[0]
visualizer.plot_name = 'All trends'
visualizer.show_or_save_to_file(add_classic_title=False, no_title=True)