Commit d37cee4b authored by Le Roux Erwan's avatar Le Roux Erwan
Browse files

[SCM] add function to compare with durand paper. refactor objects

parent c7416d72
No related merge requests found
Showing with 38 additions and 36 deletions
+38 -36
......@@ -41,7 +41,7 @@ class AbstractExtendedStudy(AbstractStudy):
""" Properties """
@property
def _year_to_daily_time_serie(self) -> OrderedDict:
def _year_to_daily_time_serie_array(self) -> OrderedDict:
return self._year_to_extended_time_serie(aggregation_function=np.mean)
@property
......@@ -50,7 +50,7 @@ class AbstractExtendedStudy(AbstractStudy):
def _year_to_extended_time_serie(self, aggregation_function) -> OrderedDict:
year_to_extended_time_serie = OrderedDict()
for year, old_time_serie in super()._year_to_daily_time_serie.items():
for year, old_time_serie in super()._year_to_daily_time_serie_array.items():
new_time_serie = np.zeros([len(old_time_serie), len(self.safran_massif_names)])
new_time_serie[:, self.nb_region_names:] = old_time_serie
for i, region_name in enumerate(self.region_names):
......
......@@ -38,7 +38,7 @@ class AbstractStudy(object):
@property
def df_all_daily_time_series_concatenated(self) -> pd.DataFrame:
df_list = [pd.DataFrame(time_serie, columns=self.safran_massif_names) for time_serie in
self.year_to_daily_time_serie.values()]
self.year_to_daily_time_serie_array.values()]
df_concatenated = pd.concat(df_list)
return df_concatenated
......@@ -63,12 +63,11 @@ class AbstractStudy(object):
for year, nc_file in sorted(nc_files, key=lambda t: t[0]):
if self.year_min <= year < self.year_max:
year_to_dataset[year] = Dataset(op.join(self.safran_full_path, nc_file))
print(year_to_dataset.keys())
return year_to_dataset
@cached_property
def year_to_daily_time_serie(self) -> OrderedDict:
return self._year_to_daily_time_serie
def year_to_daily_time_serie_array(self) -> OrderedDict:
return self._year_to_daily_time_serie_array
@cached_property
def year_to_annual_maxima(self) -> OrderedDict:
......@@ -82,7 +81,8 @@ class AbstractStudy(object):
def year_to_annual_total(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_daily_time_serie.items():
for year, time_serie in self._year_to_daily_time_serie_array.items():
print(time_serie.shape)
year_to_annual_mean[year] = self.annual_aggregation_function(time_serie, axis=0)
return year_to_annual_mean
......@@ -92,18 +92,18 @@ class AbstractStudy(object):
""" Private methods to be overwritten """
@property
def _year_to_daily_time_serie(self) -> OrderedDict:
def _year_to_daily_time_serie_array(self) -> OrderedDict:
# Map each year to a matrix of size 365-nb_days_consecutive+1 x nb_massifs
year_to_variable = {year: self.instantiate_variable_object(dataset) for year, dataset in
self.year_to_dataset_ordered_dict.items()}
year_to_daily_time_serie = OrderedDict()
year_to_daily_time_serie_array = OrderedDict()
for year in self.year_to_dataset_ordered_dict.keys():
year_to_daily_time_serie[year] = year_to_variable[year].daily_time_serie
return year_to_daily_time_serie
year_to_daily_time_serie_array[year] = year_to_variable[year].daily_time_serie_array
return year_to_daily_time_serie_array
@property
def _year_to_max_daily_time_serie(self) -> OrderedDict:
return self._year_to_daily_time_serie
return self._year_to_daily_time_serie_array
......
......@@ -9,6 +9,6 @@ class AbstractVariable(object):
self.altitude = altitude
@property
def daily_time_serie(self):
def daily_time_serie_array(self):
# Return an array of size length of time series x nb_massif
raise NotImplementedError
\ No newline at end of file
......@@ -10,9 +10,9 @@ class Crocus(AbstractStudy):
In the Crocus data, there is no 'massifsList' variable, thus we assume massifs are ordered just like Safran data
"""
def __init__(self, variable_class, altitude=1800):
def __init__(self, variable_class, *args, **kwargs):
assert variable_class in [CrocusSweVariable, CrocusDepthVariable]
super().__init__(variable_class, altitude)
super().__init__(variable_class, *args, **kwargs)
self.model_name = 'Crocus'
@property
......@@ -26,8 +26,8 @@ class Crocus(AbstractStudy):
class CrocusSwe(Crocus):
def __init__(self, altitude=1800):
super().__init__(CrocusSweVariable, altitude)
def __init__(self, *args, **kwargs):
super().__init__(CrocusSweVariable, *args, **kwargs)
class ExtendedCrocusSwe(AbstractExtendedStudy, CrocusSwe):
......@@ -36,8 +36,8 @@ class ExtendedCrocusSwe(AbstractExtendedStudy, CrocusSwe):
class CrocusDepth(Crocus):
def __init__(self, altitude=1800):
super().__init__(CrocusDepthVariable, altitude)
def __init__(self, *args, **kwargs):
super().__init__(CrocusDepthVariable, *args, **kwargs)
class ExtendedCrocusDepth(AbstractExtendedStudy, CrocusDepth):
......@@ -51,5 +51,5 @@ if __name__ == '__main__':
time_arr = np.array(d.variables['time'])
print(time_arr)
# print(d)
a = study.year_to_daily_time_serie[1960]
a = study.year_to_daily_time_serie_array[1960]
print(a.shape)
......@@ -10,7 +10,7 @@ class CrocusVariable(AbstractVariable):
self.variable_name = variable_name
@property
def daily_time_serie(self):
def daily_time_serie_array(self):
time_serie_every_6_hours = np.array(self.dataset.variables[self.variable_name])[:, 0, :]
if self.altitude == 2400:
time_serie_daily = time_serie_every_6_hours
......
......@@ -13,7 +13,7 @@ def fit_mle_gev_for_all_safran_and_different_days():
# safran = Safran(safran_alti, nb_day)
safran = ExtendedSafranSnowfall(safran_alti, nb_day)
df = safran.df_gev_mle_each_massif
df.index += ' Safran{} with {} days'.format(safran.altitude, safran.nb_days_of_snowfall)
df.index += ' Safran{} with {} days'.format(safran.altitude, safran.nb_consecutive_days)
dfs.append(df)
df = pd.concat(dfs)
path = r'/home/erwan/Documents/projects/spatiotemporalextremes/local/spatio_temporal_datasets/results/fit_mle_massif/fit_mle_gev_{}.csv'
......
......@@ -10,22 +10,23 @@ from experiment.meteo_france_SCM_study.safran.safran_variable import SafranSnowf
class Safran(AbstractStudy):
def __init__(self, variable_class: type, *args, **kwargs):
assert variable_class in [SafranSnowfallVariable, SafranPrecipitationVariable, SafranTemperatureVariable]
super().__init__(variable_class, *args, **kwargs)
self.model_name = 'Safran'
class SafranFrequency(Safran):
def __init__(self, variable_class: type, nb_days_of_snowfall=1, *args, **kwargs):
def __init__(self, variable_class: type, nb_consecutive_days=1, *args, **kwargs):
super().__init__(variable_class, *args, **kwargs)
self.nb_days_of_snowfall = nb_days_of_snowfall
self.nb_consecutive_days = nb_consecutive_days
def instantiate_variable_object(self, dataset) -> AbstractVariable:
return self.variable_class(dataset, self.nb_days_of_snowfall)
return self.variable_class(dataset, self.nb_consecutive_days)
@property
def variable_name(self):
return super().variable_name + ' cumulated over {} days'.format(self.nb_days_of_snowfall)
return super().variable_name + ' cumulated over {} days'.format(self.nb_consecutive_days)
def annual_aggregation_function(self, *args, **kwargs):
return np.sum(*args, **kwargs)
......
......@@ -35,7 +35,7 @@ class SafranSnowfallVariable(AbstractVariable):
self.daily_snowfall = [sum(hourly_snowfall[24 * i:24 * (i + 1)]) for i in range(nb_days)]
@property
def daily_time_serie(self):
def daily_time_serie_array(self):
# Aggregate the daily snowfall by the number of consecutive days
shifted_list = [self.daily_snowfall[i:] for i in range(self.nb_consecutive_days_of_snowfall)]
# First element of shifted_list is of length n, Second element of length n-1, Third element n-2....
......@@ -62,6 +62,6 @@ class SafranTemperatureVariable(AbstractVariable):
@property
def daily_time_serie(self):
def daily_time_serie_array(self):
return self.daily_temperature
......@@ -22,7 +22,7 @@ def study_iterator(study_class, only_first_one=False, both_altitude=False, verbo
for alti in AbstractStudy.ALTITUDES[::1]:
if verbose:
print('alti: {}, nb_day: {}'.format(alti, nb_day))
study = study_class(alti, nb_day) if is_safran_study else study_class(alti)
study = study_class(altitude=alti, nb_consecutive_days=nb_day) if is_safran_study else study_class(alti)
yield study
if only_first_one and not both_altitude:
break
......
......@@ -129,7 +129,7 @@ class StudyVisualizer(object):
# Plot some additional quantiles from the correspond Annual Maxima law
if self.plot_block_maxima_quantiles:
# This formula can only be applied if we have a daily time serie
assert len(self.study.year_to_daily_time_serie[1958]) in [365, 366]
assert len(self.study.year_to_daily_time_serie_array[1958]) in [365, 366]
p = p ** (1 / 365)
x_level = all_massif_data[int(p * len(all_massif_data))]
name_to_xlevel_and_color[BLOCK_MAXIMA_DISPLAY_NAME + name] = (x_level, color)
......@@ -172,10 +172,10 @@ class StudyVisualizer(object):
def get_all_massif_data(self, massif_id):
if self.year_for_kde_plot is not None:
all_massif_data = self.study.year_to_daily_time_serie[self.year_for_kde_plot][:, massif_id]
all_massif_data = self.study.year_to_daily_time_serie_array[self.year_for_kde_plot][:, massif_id]
else:
all_massif_data = np.concatenate(
[data[:, massif_id] for data in self.study.year_to_daily_time_serie.values()])
[data[:, massif_id] for data in self.study.year_to_daily_time_serie_array.values()])
all_massif_data = np.sort(all_massif_data)
return all_massif_data
......@@ -197,7 +197,7 @@ class StudyVisualizer(object):
# Counting the sum of 3-consecutive days of snowfall does not have any physical meaning,
# as we are counting twice some days
color_mean = 'g'
tuples_x_y = [(year, np.mean(data[:, massif_id])) for year, data in self.study.year_to_daily_time_serie.items()]
tuples_x_y = [(year, np.mean(data[:, massif_id])) for year, data in self.study.year_to_daily_time_serie_array.items()]
x, y = list(zip(*tuples_x_y))
x, y = average_smoothing_with_sliding_window(x, y, window_size_for_smoothing=self.window_size_for_smoothing)
ax.plot(x, y, color=color_mean)
......
......@@ -30,7 +30,7 @@ class TestSCMStudy(unittest.TestCase):
def test_scm_daily_data(self):
for study in load_scm_studies():
time_serie = study.year_to_daily_time_serie[1958]
time_serie = study.year_to_daily_time_serie_array[1958]
self.assertTrue(time_serie.ndim == 2, msg='for {} ndim={}'.format(study.__repr__(), time_serie.ndim))
self.assertTrue(len(time_serie) in [365, 366],
msg="current time serie length for {} is {}".format(study.__repr__(), len(time_serie)))
......
......@@ -96,12 +96,13 @@ def load_test_spatiotemporal_coordinates(nb_points, nb_steps, train_split_ratio=
def load_safran_studies(altitudes) -> List[Safran]:
nb_days_list = [1]
return [SafranSnowfall(safran_altitude, nb_days) for safran_altitude in altitudes for nb_days in nb_days_list]
return [SafranSnowfall(altitude=safran_altitude, nb_consecutive_days=nb_days)
for safran_altitude in altitudes for nb_days in nb_days_list]
def load_crocus_studies(altitudes) -> List[Crocus]:
crocus_classes = [CrocusSwe, CrocusDepth][:]
return [crocus_class(altitude) for crocus_class, altitude in product(crocus_classes, altitudes)]
return [crocus_class(altitude=altitude) for crocus_class, altitude in product(crocus_classes, altitudes)]
def load_scm_studies() -> List[AbstractStudy]:
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment