From c820df29fdbae0cb0375136c79a23c167abd3486 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Mon, 19 Oct 2020 11:55:51 +0200
Subject: [PATCH] Add Study For DateFirstSnowfall, and try to apply my latest
code on it
---
.../scm_models_data/abstract_study.py | 15 ++++---
.../scm_models_data/safran/safran.py | 43 ++++++++++++++-----
.../scm_models_data/safran/safran_variable.py | 21 +++++++++
.../visualization/main_study_visualizer.py | 3 +-
.../abstract_margin_estimator.py | 2 +-
.../altitudes_fit/main_altitudes_studies.py | 22 ++++++----
.../test_adamont_study.py | 1 -
7 files changed, 79 insertions(+), 28 deletions(-)
diff --git a/extreme_data/meteo_france_data/scm_models_data/abstract_study.py b/extreme_data/meteo_france_data/scm_models_data/abstract_study.py
index d706fcf9..c3bf33a4 100644
--- a/extreme_data/meteo_france_data/scm_models_data/abstract_study.py
+++ b/extreme_data/meteo_france_data/scm_models_data/abstract_study.py
@@ -339,19 +339,22 @@ class AbstractStudy(object):
year_to_daily_time_serie_array = OrderedDict()
for year in self.ordered_years:
# Check daily data
- daily_time_serie = self.year_to_variable_object[year].daily_time_serie_array
- print(daily_time_serie.shape)
- assert daily_time_serie.shape[0] in [365, 366]
- assert daily_time_serie.shape[1] == len(self.column_mask)
+ daily_time_serie = self.daily_time_series(year)
# Filter only the data corresponding:
# 1: to treturnhe start_index and last_index of the season
# 2: to the massifs for the altitude of interest
- first_index, last_index = self.year_to_first_index_and_last_index[year]
- daily_time_serie = daily_time_serie[first_index:last_index + 1, self.column_mask]
assert daily_time_serie.shape == (len(self.year_to_days[year]), len(self.study_massif_names))
year_to_daily_time_serie_array[year] = daily_time_serie
return year_to_daily_time_serie_array
+ def daily_time_series(self, year):
+ daily_time_serie = self.year_to_variable_object[year].daily_time_serie_array
+ assert daily_time_serie.shape[0] in [365, 366]
+ assert daily_time_serie.shape[1] == len(self.column_mask)
+ first_index, last_index = self.year_to_first_index_and_last_index[year]
+ daily_time_serie = daily_time_serie[first_index:last_index + 1, self.column_mask]
+ return daily_time_serie
+
""" Load Variables and Datasets """
@cached_property
diff --git a/extreme_data/meteo_france_data/scm_models_data/safran/safran.py b/extreme_data/meteo_france_data/scm_models_data/safran/safran.py
index 81a95f70..fab3a2d2 100644
--- a/extreme_data/meteo_france_data/scm_models_data/safran/safran.py
+++ b/extreme_data/meteo_france_data/scm_models_data/safran/safran.py
@@ -1,4 +1,5 @@
from collections import OrderedDict
+from typing import List, Union, Dict
import numpy as np
@@ -8,7 +9,8 @@ from extreme_data.meteo_france_data.scm_models_data.abstract_variable import Abs
from extreme_data.meteo_france_data.scm_models_data.safran.cumulated_study import CumulatedStudy
from extreme_data.meteo_france_data.scm_models_data.safran.safran_variable import SafranSnowfallVariable, \
SafranRainfallVariable, SafranTemperatureVariable, SafranTotalPrecipVariable, \
- SafranNormalizedPrecipitationRateOnWetDaysVariable, SafranNormalizedPrecipitationRateVariable
+ SafranNormalizedPrecipitationRateOnWetDaysVariable, SafranNormalizedPrecipitationRateVariable, \
+ SafranDateFirstSnowfallVariable
class Safran(AbstractStudy):
@@ -17,7 +19,8 @@ class Safran(AbstractStudy):
SafranTemperatureVariable,
SafranTotalPrecipVariable,
SafranNormalizedPrecipitationRateVariable,
- SafranNormalizedPrecipitationRateOnWetDaysVariable]
+ SafranNormalizedPrecipitationRateOnWetDaysVariable,
+ SafranDateFirstSnowfallVariable]
def __init__(self, variable_class: type, *args, **kwargs):
assert variable_class in self.SAFRAN_VARIABLES
@@ -40,6 +43,28 @@ class SafranSnowfall1Day(SafranSnowfall):
super().__init__(nb_consecutive_days=1, **kwargs)
+class SafranDateFirstSnowfall(Safran, CumulatedStudy):
+
+ def __init__(self, **kwargs):
+ super().__init__(SafranDateFirstSnowfallVariable, nb_consecutive_days=1, **kwargs)
+ self.massif_id_to_remove = set()
+ for year in self.ordered_years:
+ s = super().daily_time_series(year)
+ column_has_nan = np.isnan(s).any(axis=0)
+ index_with_nan = list(np.nonzero(column_has_nan)[0])
+ if len(index_with_nan) > 0:
+ self.massif_id_to_remove.update(set(index_with_nan))
+ self.massif_id_to_keep = tuple([i for i in range(s.shape[1])
+ if i not in self.massif_id_to_remove])
+
+ def daily_time_series(self, year):
+ return super().daily_time_series(year)[:, self.massif_id_to_keep]
+
+ @property
+ def study_massif_names(self) -> List[str]:
+ return [m for i, m in enumerate(super().study_massif_names) if i not in self.massif_id_to_remove]
+
+
class SafranSnowfall3Days(SafranSnowfall):
def __init__(self, **kwargs):
@@ -182,12 +207,10 @@ class SafranTemperature(Safran):
if __name__ == '__main__':
- altitude = 900
+ altitude = 600
year_min = 1959
- year_max = 2000
- study = SafranNormalizedPreciptationRateOnWetDays(altitude=altitude, year_min=year_min, year_max=year_max)
- d = study.year_to_dataset_ordered_dict[1959]
- print(d.keywords)
- print(d.variables.keys())
- print(study.year_to_annual_maxima[1959])
- print(study.ordered_years)
+ year_max = 2019
+ study = SafranDateFirstSnowfall(altitude=altitude, year_min=year_min, year_max=year_max)
+ print(study.study_massif_names)
+ # print(study.year_to_annual_maxima[1959])
+ # print(study.ordered_years)
diff --git a/extreme_data/meteo_france_data/scm_models_data/safran/safran_variable.py b/extreme_data/meteo_france_data/scm_models_data/safran/safran_variable.py
index 1b4df798..477f2101 100644
--- a/extreme_data/meteo_france_data/scm_models_data/safran/safran_variable.py
+++ b/extreme_data/meteo_france_data/scm_models_data/safran/safran_variable.py
@@ -67,6 +67,27 @@ class SafranSnowfallVariable(AbstractVariable):
return snowfall_in_consecutive_days
+class SafranDateFirstSnowfallVariable(SafranSnowfallVariable):
+ NAME = 'Date First Snow'
+ UNIT = 'days'
+
+ @property
+ def daily_time_serie_array(self) -> np.ndarray:
+ daily_time_series = super().daily_time_serie_array
+ new_daily_time_series = []
+ for i, s in enumerate(daily_time_series.transpose()):
+ dates_with_snow = np.nonzero(s)[0]
+ if len(dates_with_snow) > 0:
+ min = np.min(dates_with_snow)
+ else:
+ min = np.nan
+ # first_date = 1 - min / 366
+ first_date = 1 - min / 366
+ first_date_repeated = np.ones(len(s)) * first_date
+ new_daily_time_series.append(first_date_repeated)
+ new_daily_time_series_array = np.array(new_daily_time_series).transpose()
+ return new_daily_time_series_array
+
class SafranRainfallVariable(SafranSnowfallVariable):
"""Warning: this corresponds to water falling. Total precipitaiton equals Rainfall + Snowfall"""
NAME = 'Rainfall'
diff --git a/extreme_data/meteo_france_data/scm_models_data/visualization/main_study_visualizer.py b/extreme_data/meteo_france_data/scm_models_data/visualization/main_study_visualizer.py
index 276d6706..b67ca98f 100644
--- a/extreme_data/meteo_france_data/scm_models_data/visualization/main_study_visualizer.py
+++ b/extreme_data/meteo_france_data/scm_models_data/visualization/main_study_visualizer.py
@@ -14,7 +14,7 @@ from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranS
SafranRainfall, \
SafranTemperature, SafranPrecipitation, SafranSnowfall1Day, SafranSnowfall3Days, SafranSnowfall5Days, \
SafranSnowfall7Days, SafranPrecipitation1Day, SafranPrecipitation3Days, SafranPrecipitation5Days, \
- SafranPrecipitation7Days
+ SafranPrecipitation7Days, SafranDateFirstSnowfall
from extreme_data.meteo_france_data.scm_models_data.visualization.study_visualizer import \
StudyVisualizer
from projects.exceeding_snow_loads.section_discussion.crocus_study_comparison_with_eurocode import \
@@ -58,6 +58,7 @@ SCM_STUDY_CLASS_TO_ABBREVIATION = {
CrocusSnowDepthAtMaxofSwe: 'HS at max of GSL',
CrocusSnowDensity: 'Density',
SafranSnowfallSimulationRCP85: 'SF1 RCP85 projections',
+ SafranDateFirstSnowfall: 'SF1 first date'
}
altitude_massif_name_and_study_class_for_poster = [
diff --git a/extreme_fit/estimator/margin_estimator/abstract_margin_estimator.py b/extreme_fit/estimator/margin_estimator/abstract_margin_estimator.py
index 93f7a7f4..f4791e07 100644
--- a/extreme_fit/estimator/margin_estimator/abstract_margin_estimator.py
+++ b/extreme_fit/estimator/margin_estimator/abstract_margin_estimator.py
@@ -93,7 +93,7 @@ class LinearMarginEstimator(AbstractMarginEstimator):
def aic(self, split=Split.all):
aic = 2 * self.margin_model.nb_params + 2 * self.nllh(split=split)
- npt.assert_almost_equal(self.result_from_model_fit.aic, aic, decimal=1)
+ npt.assert_almost_equal(self.result_from_model_fit.aic, aic, decimal=0)
return aic
def n(self, split=Split.all):
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 e82328e3..537398a3 100644
--- a/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
+++ b/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
@@ -3,6 +3,7 @@ from typing import List
import matplotlib as mpl
+from extreme_data.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal, CrocusSnowLoad3Days
from projects.altitude_spatial_model.altitudes_fit.plots.plot_coherence_curves import plot_coherence_curves
from projects.altitude_spatial_model.altitudes_fit.plots.plot_histogram_altitude_studies import \
plot_histogram_all_models_against_altitudes, plot_histogram_all_trends_against_altitudes, \
@@ -17,25 +18,28 @@ from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitude_gr
from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.plot_total_aic import plot_individual_aic
from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day, SafranSnowfall3Days, \
- SafranSnowfall5Days, SafranSnowfall7Days
+ SafranSnowfall5Days, SafranSnowfall7Days, SafranDateFirstSnowfall, SafranPrecipitation1Day, SafranPrecipitation3Days
from extreme_data.meteo_france_data.scm_models_data.utils import Season
def main():
study_classes = [SafranSnowfall1Day, SafranSnowfall3Days, SafranSnowfall5Days, SafranSnowfall7Days][:1]
- # study_classes = [SafranPrecipitation1Day][:1]
+ # study_classes = [SafranDateFirstSnowfall]
+ # study_classes = [CrocusSnowLoadTotal]
+ study_classes = [SafranPrecipitation1Day, CrocusSnowLoadTotal, SafranDateFirstSnowfall][2:]
+ study_classes = [CrocusSnowLoad3Days, SafranSnowfall3Days, SafranSnowfall5Days]
seasons = [Season.annual, Season.winter, Season.spring, Season.automn][:1]
fast = False
if fast is None:
massif_names = None
- altitudes_list = altitudes_for_groups[:2]
+ altitudes_list = altitudes_for_groups[2:3]
elif fast:
- massif_names = ['Vanoise', 'Haute-Maurienne', 'Vercors'][:1]
- altitudes_list = altitudes_for_groups[2:]
+ massif_names = ['Vanoise', 'Haute-Maurienne', 'Vercors'][:]
+ altitudes_list = altitudes_for_groups[2:3]
else:
massif_names = None
- altitudes_list = altitudes_for_groups
+ altitudes_list = altitudes_for_groups[:]
main_loop(altitudes_list, massif_names, seasons, study_classes)
@@ -55,7 +59,7 @@ def main_loop(altitudes_list, massif_names, seasons, study_classes):
def plot_visualizers(massif_names, visualizer_list):
- # plot_histogram_all_trends_against_altitudes(massif_names, visualizer_list)
+ plot_histogram_all_trends_against_altitudes(massif_names, visualizer_list)
for relative in [True, False]:
plot_shoe_plot_changes_against_altitude(massif_names, visualizer_list, relative=relative)
# plot_coherence_curves(massif_names, visualizer_list)
@@ -64,13 +68,13 @@ 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]:
# studies.plot_mean_maxima_against_altitude(massif_names=massif_names, std=std, change=change)
# Plot the results for the model that minimizes the individual aic
- # plot_individual_aic(visualizer)
+ plot_individual_aic(visualizer)
# Plot the results for the model that minimizes the total aic
# plot_total_aic(model_classes, visualizer)
pass
diff --git a/test/test_extreme_data/test_meteo_france_data/test_adamont_study.py b/test/test_extreme_data/test_meteo_france_data/test_adamont_study.py
index eb59e3fb..4f72286b 100644
--- a/test/test_extreme_data/test_meteo_france_data/test_adamont_study.py
+++ b/test/test_extreme_data/test_meteo_france_data/test_adamont_study.py
@@ -6,7 +6,6 @@ from extreme_data.meteo_france_data.adamont_data.snowfall_simulation import Safr
class TestAdamontStudy(unittest.TestCase):
def test_year_to_date(self):
-
study = SafranSnowfallSimulationRCP85(altitude=900)
self.assertTrue(True)
--
GitLab