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 805de9df240561c04336b11be4e311e2c99b8d59..234c12ddd8ae3abfaa905b64aba3018dc311128b 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 @@ -200,8 +200,8 @@ class AbstractStudy(object): @property def all_daily_series(self) -> np.ndarray: """Return an array of approximate shape (total_number_of_days, 23) x """ - all_daily_series = np.concatenate([ time_serie_array - for time_serie_array in self.year_to_daily_time_serie_array.values()]) + all_daily_series = np.concatenate([time_serie_array + for time_serie_array in self.year_to_daily_time_serie_array.values()]) assert len(all_daily_series) == len(self.all_days) return all_daily_series @@ -211,6 +211,10 @@ class AbstractStudy(object): def observations_annual_maxima(self) -> AnnualMaxima: return AnnualMaxima(df_maxima_gev=pd.DataFrame(self.year_to_annual_maxima, index=self.study_massif_names)) + @cached_property + def observations_annual_mean(self) -> pd.DataFrame: + return pd.DataFrame(self.year_to_annual_mean, index=self.study_massif_names) + def annual_maxima_and_years(self, massif_name) -> Tuple[np.ndarray, np.ndarray]: df = self.observations_annual_maxima.df_maxima_gev return df.loc[massif_name].values, np.array(df.columns) @@ -223,6 +227,14 @@ class AbstractStudy(object): year_to_annual_maxima[year] = time_serie.max(axis=0) return year_to_annual_maxima + @cached_property + def year_to_annual_mean(self) -> OrderedDict: + # Map each year to an array of size nb_massif + year_to_annual_mean = OrderedDict() + for year, time_serie in self._year_to_max_daily_time_serie.items(): + year_to_annual_mean[year] = time_serie.mean(axis=0) + return year_to_annual_mean + @cached_property def year_to_annual_maxima_index(self) -> OrderedDict: # Map each year to an array of size nb_massif @@ -238,7 +250,7 @@ class AbstractStudy(object): ordered_index = [self.year_to_annual_maxima_index[year][i] for year in years] massif_name_to_annual_maxima_ordered_index[massif_name] = ordered_index return massif_name_to_annual_maxima_ordered_index - + @cached_property def massif_name_to_annual_maxima(self): massif_name_to_annual_maxima = OrderedDict() @@ -390,7 +402,8 @@ class AbstractStudy(object): @property def missing_massif_name(self): - return set(self.all_massif_names(self.reanalysis_path, self.dbf_filename)) - set(self.altitude_to_massif_names[self.altitude]) + return set(self.all_massif_names(self.reanalysis_path, self.dbf_filename)) - set( + self.altitude_to_massif_names[self.altitude]) @property def column_mask(self): 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 7512559d2d0450b0ed74417cbd4af30848c6c209..e0fad5e0d7a9e650ef0f64c64d2052ed09074b56 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 @@ -8,14 +8,16 @@ 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 + SafranNormalizedPrecipitationRateOnWetDaysVariable, SafranNormalizedPrecipitationRateVariable class Safran(AbstractStudy): def __init__(self, variable_class: type, *args, **kwargs): assert variable_class in [SafranSnowfallVariable, SafranRainfallVariable, SafranTemperatureVariable, - SafranTotalPrecipVariable, SafranNormalizedPrecipitationRateOnWetDaysVariable] + SafranTotalPrecipVariable, + SafranNormalizedPrecipitationRateVariable, + SafranNormalizedPrecipitationRateOnWetDaysVariable] super().__init__(variable_class, *args, **kwargs) self.model_name = 'Safran' @@ -87,6 +89,20 @@ class SafranRainfall7Days(SafranRainfall): super().__init__(nb_consecutive_days=7, **kwargs) + +class SafranNormalizedPreciptationRate(CumulatedStudy, Safran): + + def __init__(self, **kwargs): + super().__init__(SafranNormalizedPrecipitationRateVariable, **kwargs) + + def load_variable_array(self, dataset): + return [np.array(dataset.variables[k]) for k in self.load_keyword()] + + def instantiate_variable_object(self, variable_array) -> AbstractVariable: + variable_array_temperature, variable_array_snowfall, variable_array_rainfall = variable_array + return self.variable_class(variable_array_temperature, + variable_array_snowfall, variable_array_rainfall, self.nb_consecutive_days) + class SafranNormalizedPreciptationRateOnWetDays(CumulatedStudy, Safran): def __init__(self, **kwargs): 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 6c279cd3c67fe287a3237a95967a9700180235cc..66dcc8f541f51ac4415e715a9419fb86273278d5 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 @@ -77,6 +77,7 @@ class SafranRainfallVariable(SafranSnowfallVariable): class SafranTotalPrecipVariable(AbstractVariable): + NAME = 'Precipitation' def __init__(self, snow_variable_array, rain_variable_array, nb_consecutive_days): super().__init__(None) @@ -94,7 +95,10 @@ class SafranTotalPrecipVariable(AbstractVariable): return self._daily_time_serie_array -class SafranNormalizedPrecipitationRateOnWetDaysVariable(AbstractVariable): +class SafranNormalizedPrecipitationRateVariable(AbstractVariable): + NAME = 'Normalized Precip' + + def __init__(self, temperature_variable_array, snow_variable_array, rain_variable_array, nb_consecutive_days): super().__init__(None) @@ -103,8 +107,6 @@ class SafranNormalizedPrecipitationRateOnWetDaysVariable(AbstractVariable): beta = 0.06 self._daily_time_serie_array = np.exp(-beta * temperature.daily_time_serie_array) \ * total_precipitation.daily_time_serie_array - mask_for_nan_values = total_precipitation.daily_time_serie_array < 0.01 - self._daily_time_serie_array[mask_for_nan_values] = np.nan @classmethod def keyword(cls): @@ -115,6 +117,15 @@ class SafranNormalizedPrecipitationRateOnWetDaysVariable(AbstractVariable): return self._daily_time_serie_array +class SafranNormalizedPrecipitationRateOnWetDaysVariable(SafranNormalizedPrecipitationRateVariable): + + def __init__(self, temperature_variable_array, snow_variable_array, rain_variable_array, nb_consecutive_days): + super().__init__(temperature_variable_array, snow_variable_array, rain_variable_array, nb_consecutive_days) + total_precipitation = SafranTotalPrecipVariable(snow_variable_array, rain_variable_array, nb_consecutive_days) + mask_for_nan_values = total_precipitation.daily_time_serie_array < 0.01 + self._daily_time_serie_array[mask_for_nan_values] = np.nan + + class SafranTemperatureVariable(AbstractVariable): NAME = 'Temperature' UNIT = 'Celsius Degrees' diff --git a/projects/contrasting_trends_in_snow_loads/main_spatial_relative_change_in_maxima_at_fixed_altitude.py b/projects/contrasting_trends_in_snow_loads/main_spatial_relative_change_in_maxima_at_fixed_altitude.py deleted file mode 100644 index 91260c0b93074d2fc48be64ba2cd5db611d4e2e0..0000000000000000000000000000000000000000 --- a/projects/contrasting_trends_in_snow_loads/main_spatial_relative_change_in_maxima_at_fixed_altitude.py +++ /dev/null @@ -1,52 +0,0 @@ -from extreme_data.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoad3Days, \ - CrocusSnowLoadTotal -from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall, SafranRainfall, SafranPrecipitation -from extreme_data.meteo_france_data.scm_models_data.visualization.study_visualizer import \ - StudyVisualizer -import matplotlib.pyplot as plt - - -def density_wrt_altitude(): - save_to_file = True - study_class = [SafranPrecipitation, SafranRainfall, SafranSnowfall, CrocusSnowLoad3Days, CrocusSnowLoadTotal][-2] - altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][::-1] - - for altitude in altitudes: - - ax = plt.gca() - study = study_class(altitude=altitude) - # study = study_class(altitude=altitude, nb_consecutive_days=3) - massif_name_to_value = {} - for massif_name in study.study_massif_names: - s = study.observations_annual_maxima.df_maxima_gev.loc[massif_name] - year_limit = 1987 - df_before, df_after = s.loc[:year_limit], s.loc[year_limit + 1:] - df_before, df_after = df_before.mean(), df_after.mean() - # df_before, df_after = df_before.median(), df_after.median() - relative_change_value = 100 * (df_after - df_before) / df_before - massif_name_to_value[massif_name] = relative_change_value - print(massif_name_to_value) - # Plot - # massif_name_to_value = {m: i for i, m in enumerate(study.study_massif_names)} - max_values = max([abs(e) for e in massif_name_to_value.values()]) + 5 - print(max_values) - variable_name = study.variable_name - study.visualize_study(ax=ax, massif_name_to_value=massif_name_to_value, - vmin=-max_values, vmax=max_values, - add_colorbar=True, - replace_blue_by_white=False, - show=False, - label='Relative changes in mean annual maxima\n' - 'of {}\n between 1958-1987 and 1988-2017\n at {}m (%)\n' - ''.format(variable_name, study.altitude) - ) - study_visualizer = StudyVisualizer(study, save_to_file=save_to_file) - study_visualizer.plot_name = 'relative_changes_in_maxima' - study_visualizer.show_or_save_to_file() - ax.clear() - plt.close() - - -if __name__ == '__main__': - density_wrt_altitude() - # test() diff --git a/projects/contrasting_trends_in_snow_loads/spatial trends/__init__.py b/projects/contrasting_trends_in_snow_loads/spatial trends/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/projects/contrasting_trends_in_snow_loads/main_result.py b/projects/contrasting_trends_in_snow_loads/spatial trends/main_result.py similarity index 88% rename from projects/contrasting_trends_in_snow_loads/main_result.py rename to projects/contrasting_trends_in_snow_loads/spatial trends/main_result.py index ec17b4e187a170db73d9acc650a09f4448d94fe2..9df6f230e024e82169b400fd5ec5f3860b2a6579 100644 --- a/projects/contrasting_trends_in_snow_loads/main_result.py +++ b/projects/contrasting_trends_in_snow_loads/spatial trends/main_result.py @@ -2,13 +2,13 @@ from multiprocessing.pool import Pool import matplotlib as mpl -from extreme_data.meteo_france_data.scm_models_data.utils import SeasonForTheMaxima -from extreme_trend.visualizers.utils import load_altitude_to_visualizer - mpl.use('Agg') mpl.rcParams['text.usetex'] = True mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}'] +from extreme_trend.visualizers.utils import load_altitude_to_visualizer +from projects.exceeding_snow_loads.section_results.plot_trend_curves import plot_trend_map + from extreme_data.meteo_france_data.scm_models_data.abstract_study import AbstractStudy from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranPrecipitation3Days, \ SafranPrecipitation1Day, SafranPrecipitation5Days, SafranPrecipitation7Days, SafranSnowfall1Day, \ @@ -19,11 +19,8 @@ from extreme_data.meteo_france_data.scm_models_data.crocus.crocus import CrocusS CrocusSnowLoad5Days, CrocusSnowLoad7Days, CrocusSnowLoad1Day from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import \ ConfidenceIntervalMethodFromExtremes -from projects.contrasting_trends_in_snow_loads.plot_contrasting_trend_curves import plot_contrasting_trend_curves, \ - plot_contrasting_trend_curves_massif from projects.exceeding_snow_loads.section_results.main_result_trends_and_return_levels import \ compute_minimized_aic -from root_utils import NB_CORES def intermediate_result(altitudes, massif_names=None, @@ -58,7 +55,9 @@ def intermediate_result(altitudes, massif_names=None, # Plots # plot_contrasting_trend_curves(altitude_to_visualizer, all_regions=True) - plot_contrasting_trend_curves_massif(altitude_to_visualizer, all_regions=True) + # plot_contrasting_trend_curves_massif(altitude_to_visualizer, all_regions=True) + # plot_trend_curves(altitude_to_visualizer) + plot_trend_map(altitude_to_visualizer) def major_result(): @@ -68,15 +67,16 @@ def major_result(): model_subsets_for_uncertainty = None # altitudes = paper_altitudes # altitudes = paper_altitudes - altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][:] + # altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][:] + altitudes = [1200, 1500, 1800, 2100, 2400, 2700][:] snow_load_classes = [CrocusSnowLoad1Day, CrocusSnowLoad3Days, CrocusSnowLoad5Days, CrocusSnowLoad7Days][:] precipitation_classes = [SafranPrecipitation1Day, SafranPrecipitation3Days, SafranPrecipitation5Days, SafranPrecipitation7Days][:] snowfall_classes = [SafranSnowfall1Day, SafranSnowfall3Days, SafranSnowfall5Days, SafranSnowfall7Days] rainfall_classes = [SafranRainfall1Day, SafranRainfall3Days, SafranRainfall5Days, SafranRainfall7Days] study_classes = precipitation_classes + snow_load_classes - # study_classes = snowfall_classes + rainfall_classes - for study_class in snowfall_classes[:1]: + study_classes = snowfall_classes + rainfall_classes + for study_class in precipitation_classes: intermediate_result(altitudes, massif_names, model_subsets_for_uncertainty, uncertainty_methods, study_class, multiprocessing=True) diff --git a/projects/contrasting_trends_in_snow_loads/spatial trends/main_spatial_relative_change_in_maxima_at_fixed_altitude.py b/projects/contrasting_trends_in_snow_loads/spatial trends/main_spatial_relative_change_in_maxima_at_fixed_altitude.py new file mode 100644 index 0000000000000000000000000000000000000000..60a7a867f1be65c3273aacd34d4c11ccdaeed771 --- /dev/null +++ b/projects/contrasting_trends_in_snow_loads/spatial trends/main_spatial_relative_change_in_maxima_at_fixed_altitude.py @@ -0,0 +1,70 @@ +import pandas as pd +from extreme_data.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoad3Days, \ + CrocusSnowLoadTotal +from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall, SafranRainfall, \ + SafranPrecipitation, SafranPrecipitation1Day, SafranSnowfall1Day, SafranTemperature, \ + SafranNormalizedPreciptationRateOnWetDays, SafranNormalizedPreciptationRate +from extreme_data.meteo_france_data.scm_models_data.safran.safran_variable import \ + SafranNormalizedPrecipitationRateOnWetDaysVariable +from extreme_data.meteo_france_data.scm_models_data.utils import SeasonForTheMaxima +from extreme_data.meteo_france_data.scm_models_data.visualization.study_visualizer import \ + StudyVisualizer +import matplotlib.pyplot as plt + + +def relative_change_in_maxima_wrt_altitude(): + save_to_file = True + study_class = [SafranSnowfall1Day, SafranPrecipitation1Day, SafranPrecipitation, SafranRainfall, + SafranSnowfall, CrocusSnowLoad3Days, CrocusSnowLoadTotal, + SafranTemperature, SafranNormalizedPreciptationRateOnWetDays, + SafranNormalizedPreciptationRate][-1] + altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][::-1] + relative = False + + for altitude in altitudes: + + ax = plt.gca() + study = study_class(altitude=altitude, season=SeasonForTheMaxima.winter_extended) + # study = study_class(altitude=altitude, nb_consecutive_days=3) + massif_name_to_value = {} + for massif_name in study.study_massif_names: + if study_class is SafranTemperature: + s = study.observations_annual_mean.loc[massif_name] + else: + s = study.observations_annual_maxima.df_maxima_gev.loc[massif_name] + year_limit = 1989 + df_before, df_after = s.loc[:year_limit], s.loc[year_limit + 1:] + df_before, df_after = df_before.mean(), df_after.mean() + # df_before, df_after = df_before.median(), df_after.median() + if relative: + change_value = 100 * (df_after - df_before) / df_before + else: + change_value = (df_after - df_before) + massif_name_to_value[massif_name] = change_value + print(massif_name_to_value) + # Plot + # massif_name_to_value = {m: i for i, m in enumerate(study.study_massif_names)} + max_values = max([abs(e) for e in massif_name_to_value.values()]) * 1.05 + print(max_values) + variable_name = study.variable_name + prefix = 'Relative' if relative else '' + str_season = str(study.season).split('.')[-1] + study.visualize_study(ax=ax, massif_name_to_value=massif_name_to_value, + vmin=-max_values, vmax=max_values, + add_colorbar=True, + replace_blue_by_white=False, + show=False, + label='{} changes in mean {} maxima\n' + 'of {}\n between 1959-1989 and 1990-2019\n at {}m (%)\n' + ''.format(prefix, str_season, variable_name, study.altitude) + ) + study_visualizer = StudyVisualizer(study, save_to_file=save_to_file) + study_visualizer.plot_name = '{}_changes_in_maxima'.format(prefix) + study_visualizer.show_or_save_to_file() + ax.clear() + plt.close() + + +if __name__ == '__main__': + relative_change_in_maxima_wrt_altitude() + # test() diff --git a/projects/contrasting_trends_in_snow_loads/plot_contrasting_trend_curves.py b/projects/contrasting_trends_in_snow_loads/spatial trends/plot_contrasting_trend_curves.py similarity index 100% rename from projects/contrasting_trends_in_snow_loads/plot_contrasting_trend_curves.py rename to projects/contrasting_trends_in_snow_loads/spatial trends/plot_contrasting_trend_curves.py