diff --git a/extreme_estimator/R_fit/gev_fit/abstract_margin_model.py b/extreme_estimator/R_fit/gev_fit/abstract_margin_model.py
new file mode 100644
index 0000000000000000000000000000000000000000..dc3e41a7ef16a8d90576d8cd54671a41b6b87765
--- /dev/null
+++ b/extreme_estimator/R_fit/gev_fit/abstract_margin_model.py
@@ -0,0 +1,77 @@
+from extreme_estimator.R_fit.gev_fit.gev_mle_fit import GevMleFit, mle_gev
+from extreme_estimator.R_fit.utils import get_loaded_r
+
+
+def frechet_unitary_transformation(data, location, scale, shape):
+    """
+    Compute the unitary Frechet transformed data
+    (1 + \zeta \frac{z - \mu}{\sigma}) ^ {\frac{1}{\zeta}}
+    """
+    assert False
+    # todo: there is already a function doing that in R
+    return (1 + shape * (data - location) / scale) ** (1 / shape)
+
+
+class GevParameters(object):
+
+    def __init__(self, location, scale, shape):
+        self.location = location
+        self.scale = scale
+        self.shape = shape
+
+
+def frechet_unitary_transformation_from_gev_parameters(data, gev_parameters: GevParameters):
+    return frechet_unitary_transformation(data, gev_parameters.location)
+
+
+class AbstractMarginModel(object):
+    GEV_SCALE = GevMleFit.GEV_SCALE
+    GEV_LOCATION = GevMleFit.GEV_LOCATION
+    GEV_SHAPE = GevMleFit.GEV_SHAPE
+    GEV_PARAMETERS = [GEV_LOCATION, GEV_SCALE, GEV_SHAPE]
+
+    def __init__(self):
+        """
+        Class to fit a GEV a list of data. Compute also the corresponding unitary data
+
+        :param coordinate: Represents the spatio-temporal spatial_coordinates of the marginals
+        :param data:  array of data corresponding to this position (and potentially its neighborhood)
+        """
+        self.default_params = {gev_param: 1.0 for gev_param in self.GEV_PARAMETERS}
+        self.r = get_loaded_r()
+
+    # Define the method to sample/fit a single gev
+
+    def rgev(self, nb_obs, loc, scale, shape):
+        gev_params = {
+            self.GEV_LOCATION: loc,
+            self.GEV_SCALE: scale,
+            self.GEV_SHAPE: shape,
+        }
+        return self.r.rgev(nb_obs, **gev_params)
+
+    def fitgev(self, x_gev, estimator=GevMleFit):
+        mle_params = mle_gev(x_gev=x_gev)
+
+    # Define the method to sample/fit all marginals globally in the child classes
+
+    def fitmargin(self, maxima, coord):
+        df_gev_params = None
+        return df_gev_params
+
+    def rmargin(self, nb_obs, coord):
+        pass
+
+    def get_maxima(self, maxima_normalized, coord):
+        pass
+
+    def get_maxima_normalized(self, maxima, df_gev_params):
+        pass
+
+
+class SmoothMarginModel(AbstractMarginModel):
+    pass
+
+
+class ConstantMarginModel(SmoothMarginModel):
+    pass
diff --git a/extreme_estimator/R_fit/gev_fit/gev_marginal.py b/extreme_estimator/R_fit/gev_fit/gev_marginal.py
deleted file mode 100644
index 8b15be5850ea08103b6c5f20d437baf9f5b2c7c4..0000000000000000000000000000000000000000
--- a/extreme_estimator/R_fit/gev_fit/gev_marginal.py
+++ /dev/null
@@ -1,51 +0,0 @@
-from R.gev_fit.gev_mle_fit import GevMleFit
-
-
-def frechet_unitary_transformation(data, location, scale, shape):
-    """
-    Compute the unitary Frechet transformed data
-    (1 + \zeta \frac{z - \mu}{\sigma}) ^ {\frac{1}{\zeta}}
-    """
-    # todo: there is already a function doing that in R
-    return (1 + shape * (data - location) / scale) ** (1 / shape)
-
-
-class GevParameters(object):
-
-    def __init__(self, location, scale, shape):
-        self.location = location
-        self.scale = scale
-        self.shape = shape
-
-
-def frechet_unitary_transformation_from_gev_parameters(data, gev_parameters: GevParameters):
-    return frechet_unitary_transformation(data, gev_parameters.location)
-
-
-class GevMarginal(object):
-
-    def __init__(self, coordinate, data, estimator=GevMleFit):
-        """
-        Class to fit a GEV a list of data. Compute also the corresponding unitary data
-
-        :param coordinate: Represents the spatio-temporal spatial_coordinates of the marginals
-        :param data:  array of data corresponding to this position (and potentially its neighborhood)
-        """
-        self.coordinate = coordinate
-        self.data = data
-        self.gev_estimator = estimator(x_gev=data)
-        self.gev_parameters_estimated = [self.location, self.scale, self.shape]
-        self.frechet_unitary_data = frechet_unitary_transformation(data=data, location=self.location,
-                                                                   scale=self.scale, shape=self.shape)
-
-    @property
-    def location(self):
-        return self.gev_estimator.location
-
-    @property
-    def scale(self):
-        return self.gev_estimator.scale
-
-    @property
-    def shape(self):
-        return self.gev_estimator.shape
diff --git a/extreme_estimator/R_fit/gev_fit/gev_mle_fit.py b/extreme_estimator/R_fit/gev_fit/gev_mle_fit.py
index e56280b3bdd1703be99fb7ac4156cef9e065c2ee..bfd5413ad37744fa67dad3392eb7519c615ec164 100644
--- a/extreme_estimator/R_fit/gev_fit/gev_mle_fit.py
+++ b/extreme_estimator/R_fit/gev_fit/gev_mle_fit.py
@@ -6,10 +6,6 @@ import os.path as op
 # Defining some constants
 from extreme_estimator.R_fit.utils import get_associated_r_file
 
-GEV_SCALE = 'scale'
-GEV_LOCATION = 'loc'
-GEV_SHAPE = 'shape'
-
 
 def mle_gev(x_gev: np.ndarray, start_loc=0, start_scale=1, start_shape=0):
     assert np.ndim(x_gev) == 1
@@ -25,10 +21,13 @@ def mle_gev(x_gev: np.ndarray, start_loc=0, start_scale=1, start_shape=0):
 
 
 class GevMleFit(object):
+    GEV_SCALE = 'scale'
+    GEV_LOCATION = 'loc'
+    GEV_SHAPE = 'shape'
 
     def __init__(self, x_gev: np.ndarray, start_loc=0, start_scale=1, start_shape=0):
         self.x_gev = x_gev
         self.mle_params = mle_gev(x_gev, start_loc, start_scale, start_shape)
-        self.shape = self.mle_params[GEV_SHAPE]
-        self.scale = self.mle_params[GEV_SCALE]
-        self.location = self.mle_params[GEV_LOCATION]
\ No newline at end of file
+        self.shape = self.mle_params[self.GEV_SHAPE]
+        self.scale = self.mle_params[self.GEV_SCALE]
+        self.location = self.mle_params[self.GEV_LOCATION]
diff --git a/extreme_estimator/R_fit/max_stable_fit/abstract_max_stable_model.py b/extreme_estimator/R_fit/max_stable_fit/abstract_max_stable_model.py
index bd2c3723d7275024f4039062afb38ca15d73941d..8009570bb0a018a5ea11f11a2c5c56d397f891bc 100644
--- a/extreme_estimator/R_fit/max_stable_fit/abstract_max_stable_model.py
+++ b/extreme_estimator/R_fit/max_stable_fit/abstract_max_stable_model.py
@@ -16,7 +16,8 @@ class AbstractMaxStableModel(object):
         self.user_params_sample = params_sample
         self.r = get_loaded_r()
 
-    def fitmaxstab(self, maxima: np.ndarray, coord: np.ndarray, fit_marge=False):
+    def fitmaxstab(self, maxima_normalized: np.ndarray, coord: np.ndarray, fit_marge=False):
+        assert all([isinstance(arr, np.ndarray) for arr in [maxima_normalized, coord]])
         #  Specify the fit params
         fit_params = {
             'fit.marge': fit_marge,
@@ -25,7 +26,7 @@ class AbstractMaxStableModel(object):
         # Run the fitmaxstab in R
         # todo: find how to specify the optim function to use
         try:
-            res = self.r.fitmaxstab(np.transpose(maxima), coord, **self.cov_mod_param, **fit_params)  # type: ListVector
+            res = self.r.fitmaxstab(np.transpose(maxima_normalized), coord, **self.cov_mod_param, **fit_params)  # type: ListVector
         except rpy2.rinterface.RRuntimeError as error:
             raise Exception('Some R exception have been launched at RunTime: {}'.format(error.__repr__()))
         # todo: maybe if the convergence was not successful I could try other starting point several times
@@ -81,10 +82,4 @@ class AbstractMaxStableModelWithCovarianceFunction(AbstractMaxStableModel):
             'range': 3,
             'smooth': 0.5,
             'nugget': 0.5
-        }
-
-
-if __name__ == '__main__':
-     print([c.name for c in CovarianceFunction])
-     for covariance_function in CovarianceFunction:
-         print(covariance_function)
\ No newline at end of file
+        }
\ No newline at end of file
diff --git a/extreme_estimator/robustness_plot/msp_robustness.py b/extreme_estimator/robustness_plot/msp_robustness.py
index b73fb0f43ee8a5c378166ed33002835e3ed5f98f..5940a62e37f8c0a9bd28b6168739d86e4ad33cc6 100644
--- a/extreme_estimator/robustness_plot/msp_robustness.py
+++ b/extreme_estimator/robustness_plot/msp_robustness.py
@@ -1,13 +1,13 @@
 from extreme_estimator.R_fit.max_stable_fit.abstract_max_stable_model import AbstractMaxStableModel, CovarianceFunction
 from extreme_estimator.R_fit.max_stable_fit.max_stable_models import Smith, BrownResnick, Schlather, ExtremalT
 from extreme_estimator.estimator.abstract_estimator import AbstractEstimator
-from extreme_estimator.estimator.unitary_msp_estimator import MaxStableEstimator
+from extreme_estimator.estimator.max_stable_estimator import MaxStableEstimator
 from extreme_estimator.robustness_plot.multiple_plot import MultiplePlot
 from extreme_estimator.robustness_plot.single_plot import SinglePlot
 from spatio_temporal_dataset.dataset.simulation_dataset import SimulatedDataset
-from spatio_temporal_dataset.spatial_coordinates.abstract_coordinates import AbstractSpatialCoordinates
-from spatio_temporal_dataset.spatial_coordinates.alps_station_coordinates import \
-    AlpsStationCoordinatesBetweenZeroAndOne, AlpsStationCoordinatesBetweenZeroAndTwo
+from spatio_temporal_dataset.spatial_coordinates.abstract_spatial_coordinates import AbstractSpatialCoordinates
+from spatio_temporal_dataset.spatial_coordinates.alps_station_2D_coordinates import \
+    AlpsStation2DCoordinatesBetweenZeroAndOne, AlpsStationCoordinatesBetweenZeroAndTwo
 from spatio_temporal_dataset.spatial_coordinates.generated_coordinates import CircleCoordinatesRadius1, \
     CircleCoordinatesRadius2
 from extreme_estimator.robustness_plot.display_item import DisplayItem
@@ -102,7 +102,7 @@ def multiple_spatial_robustness_alps():
         MspSpatial.MaxStableModelItem.name: msp_models,
         MspSpatial.SpatialCoordinateClassItem.name: [CircleCoordinatesRadius1,
                                                      CircleCoordinatesRadius2,
-                                                     AlpsStationCoordinatesBetweenZeroAndOne,
+                                                     AlpsStation2DCoordinatesBetweenZeroAndOne,
                                                      AlpsStationCoordinatesBetweenZeroAndTwo][:],
     })
 
diff --git a/spatio_temporal_dataset/dataset/abstract_dataset.py b/spatio_temporal_dataset/dataset/abstract_dataset.py
index b270fc133b29e9f9662eb7393f28e837ed9422b1..855389c6a886dfe6fb341dbdef57681ea8d409cc 100644
--- a/spatio_temporal_dataset/dataset/abstract_dataset.py
+++ b/spatio_temporal_dataset/dataset/abstract_dataset.py
@@ -1,24 +1,27 @@
 import os
+import numpy as np
 import os.path as op
 import pandas as pd
-from spatio_temporal_dataset.temporal_maxima.temporal_maxima import TemporalMaxima
-from spatio_temporal_dataset.spatial_coordinates.abstract_coordinates import AbstractSpatialCoordinates
+from spatio_temporal_dataset.temporal_observations.abstract_temporal_observations import AbstractTemporalObservations
+from spatio_temporal_dataset.spatial_coordinates.abstract_spatial_coordinates import AbstractSpatialCoordinates
 
 
 class AbstractDataset(object):
 
-    def __init__(self, temporal_maxima: TemporalMaxima, spatial_coordinates: AbstractSpatialCoordinates):
-        assert (temporal_maxima.index == spatial_coordinates.index).all()
-        self.temporal_maxima = temporal_maxima
+    def __init__(self, temporal_observations: AbstractTemporalObservations, spatial_coordinates: AbstractSpatialCoordinates):
+        # assert
+        # is_same_index = temporal_observations.index == spatial_coordinates.index  # type: pd.Series
+        # assert is_same_index.all()
+        self.temporal_observations = temporal_observations
         self.spatial_coordinates = spatial_coordinates
 
     @classmethod
     def from_csv(cls, csv_path: str):
         assert op.exists(csv_path)
         df = pd.read_csv(csv_path)
-        temporal_maxima = TemporalMaxima.from_df(df)
+        temporal_maxima = AbstractTemporalObservations.from_df(df)
         spatial_coordinates = AbstractSpatialCoordinates.from_df(df)
-        return cls(temporal_maxima=temporal_maxima, spatial_coordinates=spatial_coordinates)
+        return cls(temporal_maxima, spatial_coordinates)
 
     def to_csv(self, csv_path: str):
         dirname = op.dirname(csv_path)
@@ -29,16 +32,21 @@ class AbstractDataset(object):
     @property
     def df_dataset(self) -> pd.DataFrame:
         # Merge dataframes with the maxima and with the coordinates
-        return self.temporal_maxima.df_maxima.join(self.spatial_coordinates.df_coord)
+        return self.temporal_observations.df_maxima.join(self.spatial_coordinates.df_coord)
 
     @property
     def coord(self):
         return self.spatial_coordinates.coord
 
     @property
-    def maxima(self):
-        return self.temporal_maxima.maxima
+    def maxima(self) -> np.ndarray:
+        return self.temporal_observations.maxima
 
+    @property
+    def maxima_normalized(self):
+        return self.temporal_observations.maxima_normalized
+
+    @maxima_normalized.setter
+    def maxima_normalized(self, maxima_normalized_to_set):
+        self.temporal_observations.maxima_normalized = maxima_normalized_to_set
 
-class RealDataset(AbstractDataset):
-    pass
diff --git a/spatio_temporal_dataset/dataset/simulation_dataset.py b/spatio_temporal_dataset/dataset/simulation_dataset.py
index 6f8a1a6f5a0b8fda64103a7466fff1ce65a01032..e76668b282d94ad2a775a73d310ced359bde5f3f 100644
--- a/spatio_temporal_dataset/dataset/simulation_dataset.py
+++ b/spatio_temporal_dataset/dataset/simulation_dataset.py
@@ -1,26 +1,55 @@
+from extreme_estimator.R_fit.gev_fit.abstract_margin_model import AbstractMarginModel
 from extreme_estimator.R_fit.max_stable_fit.abstract_max_stable_model import AbstractMaxStableModel
-import pandas as pd
 from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
-from spatio_temporal_dataset.temporal_maxima.temporal_maxima import TemporalMaxima
-from spatio_temporal_dataset.spatial_coordinates.abstract_coordinates import AbstractSpatialCoordinates
-from spatio_temporal_dataset.spatial_coordinates.generated_coordinates import CircleCoordinatesRadius1
+from spatio_temporal_dataset.spatial_coordinates.abstract_spatial_coordinates import AbstractSpatialCoordinates
+from spatio_temporal_dataset.temporal_observations.abstract_temporal_observations import AbstractTemporalObservations
+from spatio_temporal_dataset.temporal_observations.annual_maxima_observations import \
+    MaxStableAnnualMaxima, AnnualMaxima, MarginAnnualMaxima, FullAnnualMaxima
 
 
 class SimulatedDataset(AbstractDataset):
-
-    def __init__(self, temporal_maxima: TemporalMaxima, spatial_coordinates: AbstractSpatialCoordinates,
-                 max_stable_model: AbstractMaxStableModel):
-        super().__init__(temporal_maxima, spatial_coordinates)
+    """
+    Class SimulatedDataset gives access to:
+        -the max_stable_model AND/OR marginal_model that was used for sampling
+    """
+
+    def __init__(self, temporal_observations: AbstractTemporalObservations,
+                 spatial_coordinates: AbstractSpatialCoordinates,
+                 max_stable_model: AbstractMaxStableModel = None,
+                 margin_model: AbstractMarginModel = None):
+        super().__init__(temporal_observations, spatial_coordinates)
+        assert margin_model is not None or max_stable_model is not None
+        self.margin_model = margin_model
         self.max_stable_model = max_stable_model
 
+
+class MaxStableDataset(SimulatedDataset):
+
     @classmethod
-    def from_max_stable_sampling(cls, nb_obs: int, max_stable_model:AbstractMaxStableModel, spatial_coordinates: AbstractSpatialCoordinates):
-        maxima = max_stable_model.rmaxstab(nb_obs=nb_obs, coord=spatial_coordinates.coord)
-        df_maxima = pd.DataFrame(data=maxima, index=spatial_coordinates.index)
-        temporal_maxima = TemporalMaxima(df_maxima=df_maxima)
-        return cls(temporal_maxima=temporal_maxima, spatial_coordinates=spatial_coordinates, max_stable_model=max_stable_model)
+    def from_sampling(cls, nb_obs: int, max_stable_model: AbstractMaxStableModel,
+                      spatial_coordinates: AbstractSpatialCoordinates):
+        temporal_obs = MaxStableAnnualMaxima.from_sampling(nb_obs, max_stable_model, spatial_coordinates)
+        return cls(temporal_observations=temporal_obs,
+                   spatial_coordinates=spatial_coordinates,
+                   max_stable_model=max_stable_model)
 
 
+class MarginDataset(SimulatedDataset):
+
+    @classmethod
+    def from_sampling(cls, nb_obs: int, margin_model: AbstractMarginModel,
+                      spatial_coordinates: AbstractSpatialCoordinates):
+        temporal_obs = MarginAnnualMaxima.from_sampling(nb_obs, spatial_coordinates, margin_model)
+        return cls(temporal_observations=temporal_obs,
+                   spatial_coordinates=spatial_coordinates,
+                   margin_model=margin_model)
 
 
+class FullSimulatedDataset(SimulatedDataset):
 
+    @classmethod
+    def from_sampling(cls, nb_obs: int, max_stable_model: AbstractMaxStableModel,
+                      spatial_coordinates: AbstractSpatialCoordinates):
+        temporal_obs = FullAnnualMaxima.from_double_sampling(nb_obs, max_stable_model,
+                                                             spatial_coordinates)
+        return cls(temporal_obs, spatial_coordinates, max_stable_model)
diff --git a/spatio_temporal_dataset/marginals/spatial_marginals.py b/spatio_temporal_dataset/marginals/spatial_marginals.py
index ef411a53da08a8cc530e8ab9f0a653fc3b66e7dd..91b2bcf053a23c68e9c6176b18551e4cd45802c8 100644
--- a/spatio_temporal_dataset/marginals/spatial_marginals.py
+++ b/spatio_temporal_dataset/marginals/spatial_marginals.py
@@ -1,13 +1,14 @@
-from typing import List
-
 from spatio_temporal_dataset.marginals.abstract_marginals import AbstractMarginals
-from R.gev_fit.gev_marginal import GevMarginal
-from spatio_temporal_dataset.stations.station import Station
+
+from spatio_temporal_dataset.spatial_coordinates.abstract_spatial_coordinates import AbstractSpatialCoordinates
 
 
 class SpatialMarginal(AbstractMarginals):
+    """The main idea is to have on marginal per station"""
+
+    def __init__(self, spatial_coordinates: AbstractSpatialCoordinates):
+        self.spatial_coordinates = spatial_coordinates
+
 
-    def __init__(self, stations: List[Station]):
-        super().__init__(stations)
-        for station in self.stations:
-            self.gev_marginals.append(GevMarginal(coordinate=station, data=station.annual_maxima.values))
\ No newline at end of file
+class SimulatedSpatialMarginal(SpatialMarginal):
+    pass
diff --git a/spatio_temporal_dataset/max_stable/abstract_max_stable.py b/spatio_temporal_dataset/max_stable/abstract_max_stable.py
deleted file mode 100644
index a2fe3404415be5a662a59348957526422383e779..0000000000000000000000000000000000000000
--- a/spatio_temporal_dataset/max_stable/abstract_max_stable.py
+++ /dev/null
@@ -1,10 +0,0 @@
-
-
-class AbstractMaxStableOptimization(object):
-
-    def __init__(self, marginals):
-        pass
-        """
-        We can conclude on the dimension of the max stable
-        """
-
diff --git a/spatio_temporal_dataset/temporal_maxima/temporal_maxima.py b/spatio_temporal_dataset/temporal_maxima/temporal_maxima.py
deleted file mode 100644
index a7c19cb339e1bfeb311e934531b9f561b8a89754..0000000000000000000000000000000000000000
--- a/spatio_temporal_dataset/temporal_maxima/temporal_maxima.py
+++ /dev/null
@@ -1,24 +0,0 @@
-import pandas as pd
-
-class TemporalMaxima(object):
-
-    def __init__(self, df_maxima: pd.DataFrame):
-        """
-        Main attribute of the class is the DataFrame df_maxima
-        Index are stations index, Columns are the year of the maxima
-        """
-        self.df_maxima = df_maxima
-
-    @classmethod
-    def from_df(cls, df):
-        pass
-
-    @property
-    def index(self):
-        return self.df_maxima.index
-
-    @property
-    def maxima(self):
-        return self.df_maxima.values
-
-    # todo: add the transformed_maxima and the not-trasnformed maxima
diff --git a/spatio_temporal_dataset/temporal_maxima/__init__.py b/spatio_temporal_dataset/temporal_observations/__init__.py
similarity index 100%
rename from spatio_temporal_dataset/temporal_maxima/__init__.py
rename to spatio_temporal_dataset/temporal_observations/__init__.py
diff --git a/spatio_temporal_dataset/temporal_observations/abstract_temporal_observations.py b/spatio_temporal_dataset/temporal_observations/abstract_temporal_observations.py
new file mode 100644
index 0000000000000000000000000000000000000000..c9d7c9e210a032ebf545355dfa43efa6d560b05d
--- /dev/null
+++ b/spatio_temporal_dataset/temporal_observations/abstract_temporal_observations.py
@@ -0,0 +1,51 @@
+import pandas as pd
+
+
+class AbstractTemporalObservations(object):
+
+    def __init__(self, df_maxima_normalized: pd.DataFrame = None, df_maxima: pd.DataFrame = None):
+        """
+        Main attribute of the class is the DataFrame df_maxima
+        Index are stations index
+        Columns are the temporal moment of the maxima
+        """
+        self.df_maxima_normalized = df_maxima_normalized
+        self.df_maxima = df_maxima
+
+    @classmethod
+    def from_df(cls, df):
+        pass
+
+    @property
+    def maxima(self):
+        return self.df_maxima.values
+
+    @property
+    def maxima_normalized(self):
+        return self.df_maxima_normalized.values
+
+    @maxima_normalized.setter
+    def maxima_normalized(self, maxima_normalized_to_set):
+        assert self.df_maxima_normalized is None
+        assert maxima_normalized_to_set is not None
+        assert maxima_normalized_to_set.shape == self.maxima.shape
+        self.df_maxima_normalized = pd.DataFrame(data=maxima_normalized_to_set, index=self.df_maxima.index,
+                                                 columns=self.df_maxima.columns)
+
+    @property
+    def column_to_time_index(self):
+        pass
+
+    @property
+    def index(self):
+        return self.df_maxima.index
+
+
+class RealTemporalObservations(object):
+
+    def __init__(self):
+        pass
+
+
+class NormalizedTemporalObservations(object):
+    pass
diff --git a/spatio_temporal_dataset/temporal_observations/alps_precipitation_observations.py b/spatio_temporal_dataset/temporal_observations/alps_precipitation_observations.py
new file mode 100644
index 0000000000000000000000000000000000000000..80d6a8e09c98d25d5dda1a0c3040bf8fae02746d
--- /dev/null
+++ b/spatio_temporal_dataset/temporal_observations/alps_precipitation_observations.py
@@ -0,0 +1,5 @@
+from spatio_temporal_dataset.temporal_observations.abstract_temporal_observations import AbstractTemporalObservations
+
+
+class AlpsPrecipitationObservations(AbstractTemporalObservations):
+    pass
\ No newline at end of file
diff --git a/spatio_temporal_dataset/temporal_observations/annual_maxima_observations.py b/spatio_temporal_dataset/temporal_observations/annual_maxima_observations.py
new file mode 100644
index 0000000000000000000000000000000000000000..7964e829fa2b9ede4e5cf7071a6c42400563d560
--- /dev/null
+++ b/spatio_temporal_dataset/temporal_observations/annual_maxima_observations.py
@@ -0,0 +1,47 @@
+import pandas as pd
+
+from extreme_estimator.R_fit.gev_fit.abstract_margin_model import AbstractMarginModel
+from extreme_estimator.R_fit.max_stable_fit.abstract_max_stable_model import AbstractMaxStableModel
+from spatio_temporal_dataset.spatial_coordinates.abstract_spatial_coordinates import AbstractSpatialCoordinates
+from spatio_temporal_dataset.temporal_observations.abstract_temporal_observations import AbstractTemporalObservations
+
+
+class AnnualMaxima(AbstractTemporalObservations):
+    """
+    Index are stations index
+    Columns are the annual of the maxima
+    """
+    pass
+
+
+class MarginAnnualMaxima(AnnualMaxima):
+
+    @classmethod
+    def from_sampling(cls, nb_obs: int, spatial_coordinates: AbstractSpatialCoordinates,
+                      margin_model: AbstractMarginModel):
+        maxima = margin_model.rmargin(nb_obs=nb_obs, coord=spatial_coordinates.coord)
+        df_maxima = pd.DataFrame(data=maxima, index=spatial_coordinates.index)
+        return cls(df_maxima=df_maxima)
+
+
+class MaxStableAnnualMaxima(AbstractTemporalObservations):
+
+    @classmethod
+    def from_sampling(cls, nb_obs: int, max_stable_model: AbstractMaxStableModel,
+                      spatial_coordinates: AbstractSpatialCoordinates):
+        maxima_normalized = max_stable_model.rmaxstab(nb_obs=nb_obs, coord=spatial_coordinates.coord)
+        df_maxima_normalized = pd.DataFrame(data=maxima_normalized, index=spatial_coordinates.index)
+        return cls(df_maxima_normalized=df_maxima_normalized)
+
+
+class FullAnnualMaxima(MaxStableAnnualMaxima):
+
+    @classmethod
+    def from_double_sampling(cls, nb_obs: int, max_stable_model: AbstractMaxStableModel,
+                             spatial_coordinates: AbstractSpatialCoordinates,
+                             margin_model: AbstractMarginModel):
+        max_stable_annual_maxima = super().from_sampling(nb_obs, max_stable_model, spatial_coordinates)
+        #  Compute df_maxima from df_maxima_normalized
+        maxima = margin_model.get_maxima(max_stable_annual_maxima.maxima_normalized, spatial_coordinates.coord)
+        max_stable_annual_maxima.df_maxima = pd.DataFrame(data=maxima, index=spatial_coordinates.index)
+        return max_stable_annual_maxima