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 d706fcf9a9b6eda9c7622f16d409302ea5b59019..c3bf33a42519e7c1385e7178bf526d3dff73f3b0 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 81a95f7083a1a8e924e0b63665a8af2cde305b01..fab3a2d26d9f3fd5df2f6ddce2778780a0d7daeb 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 1b4df798dee5fdfc620565fee6490927dbf5b69c..477f210111e75641dfa815302ebe193ae43ae6ff 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 276d670676f387333afdfb7329c5e69352780c20..b67ca98fc9692eb89e6cc57a9bb74ccdd652428d 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 93f7a7f4fae2dfc97e8ff7f7eb28975b001c0c0c..f4791e07d45dd8e23303b6241191fe88ce2aeefd 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 e82328e34e99a991840e60d5523dd5478165e484..537398a3de1620dc415b623430446ddbae4d41cf 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 eb59e3fb4ae1179548960a4bb9ee9c0c1bbf5a48..4f72286be42c8d1f92c015aefb786106fc7d8fc9 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)