diff --git a/spatio_temporal_dataset/coordinates/axis_coordinates/__init__.py b/experiment/__init__.py
similarity index 100%
rename from spatio_temporal_dataset/coordinates/axis_coordinates/__init__.py
rename to experiment/__init__.py
diff --git a/spatio_temporal_dataset/marginals/__init__.py b/experiment/robustness_plot/__init__.py
similarity index 100%
rename from spatio_temporal_dataset/marginals/__init__.py
rename to experiment/robustness_plot/__init__.py
diff --git a/extreme_estimator/robustness_plot/display_item.py b/experiment/robustness_plot/display_item.py
similarity index 100%
rename from extreme_estimator/robustness_plot/display_item.py
rename to experiment/robustness_plot/display_item.py
diff --git a/spatio_temporal_dataset/stations/__init__.py b/experiment/robustness_plot/estimation_robustness/__init__.py
similarity index 100%
rename from spatio_temporal_dataset/stations/__init__.py
rename to experiment/robustness_plot/estimation_robustness/__init__.py
diff --git a/experiment/robustness_plot/estimation_robustness/max_stable_process_plot.py b/experiment/robustness_plot/estimation_robustness/max_stable_process_plot.py
new file mode 100644
index 0000000000000000000000000000000000000000..23b35b5b236c8fae7847dd37af78062d0dcc37ee
--- /dev/null
+++ b/experiment/robustness_plot/estimation_robustness/max_stable_process_plot.py
@@ -0,0 +1,55 @@
+from extreme_estimator.estimator.max_stable_estimator import MaxStableEstimator
+from extreme_estimator.extreme_models.max_stable_model.abstract_max_stable_model import AbstractMaxStableModel
+from extreme_estimator.extreme_models.max_stable_model.max_stable_models import Smith
+from experiment.robustness_plot.display_item import DisplayItem
+from experiment.robustness_plot.multiple_plot import MultiplePlot
+from experiment.robustness_plot.single_plot import SinglePlot
+from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
+from spatio_temporal_dataset.coordinates.spatial_coordinates.generated_coordinates import CircleCoordinatesRadius1
+from spatio_temporal_dataset.dataset.simulation_dataset import MaxStableDataset
+
+
+class MaxStableDisplayItem(DisplayItem):
+
+    def display_name_from_value(self, value: AbstractMaxStableModel):
+        return value.cov_mod
+
+
+class CoordinateDisplayItem(DisplayItem):
+
+    def display_name_from_value(self, value: AbstractCoordinates):
+        return str(value).split('.')[-1].split("'")[0]
+
+
+class MaxStableProcessPlot(object):
+    MaxStableModelItem = MaxStableDisplayItem('max_stable_model', Smith)
+    CoordinateClassItem = CoordinateDisplayItem('coordinate_class', CircleCoordinatesRadius1)
+    NbStationItem = DisplayItem('Number of stations', 50)
+    NbObservationItem = DisplayItem('nb_obs', 60)
+
+    def msp_spatial_ordinates(self, **kwargs_single_point) -> dict:
+        # Get the argument from kwargs
+        max_stable_model = self.MaxStableModelItem.value_from_kwargs(
+            **kwargs_single_point)  # type: AbstractMaxStableModel
+        coordinate_class = self.CoordinateClassItem.value_from_kwargs(**kwargs_single_point)
+        nb_station = self.NbStationItem.value_from_kwargs(**kwargs_single_point)
+        nb_obs = self.NbObservationItem.value_from_kwargs(**kwargs_single_point)
+        # Run the estimation
+        spatial_coordinates = coordinate_class.from_nb_points(nb_points=nb_station)
+        dataset = MaxStableDataset.from_sampling(nb_obs=nb_obs, max_stable_model=max_stable_model,
+                                                 coordinates=spatial_coordinates)
+        estimator = MaxStableEstimator(dataset, max_stable_model)
+        estimator.fit()
+        return estimator.scalars(max_stable_model.params_sample)
+
+
+class SingleMaxStableProcessPlot(SinglePlot, MaxStableProcessPlot):
+
+    def compute_value_from_kwargs_single_point(self, **kwargs_single_point):
+        return self.msp_spatial_ordinates(**kwargs_single_point)
+
+
+class MultipleMaxStableProcessPlot(MultiplePlot, MaxStableProcessPlot):
+
+    def compute_value_from_kwargs_single_point(self, **kwargs_single_point):
+        return self.msp_spatial_ordinates(**kwargs_single_point)
\ No newline at end of file
diff --git a/spatio_temporal_dataset/marginals/spatio_temporal_marginals.py b/experiment/robustness_plot/estimation_robustness/spatial_robustness/__init__.py
similarity index 100%
rename from spatio_temporal_dataset/marginals/spatio_temporal_marginals.py
rename to experiment/robustness_plot/estimation_robustness/spatial_robustness/__init__.py
diff --git a/experiment/robustness_plot/estimation_robustness/spatial_robustness/alps_msp_robustness.py b/experiment/robustness_plot/estimation_robustness/spatial_robustness/alps_msp_robustness.py
new file mode 100644
index 0000000000000000000000000000000000000000..69b4d61921955387ce5d2fbad6f481dec63cbd2e
--- /dev/null
+++ b/experiment/robustness_plot/estimation_robustness/spatial_robustness/alps_msp_robustness.py
@@ -0,0 +1,61 @@
+from extreme_estimator.estimator.abstract_estimator import AbstractEstimator
+from extreme_estimator.extreme_models.max_stable_model.max_stable_models import Smith, BrownResnick
+from experiment.robustness_plot.estimation_robustness.max_stable_process_plot import MultipleMaxStableProcessPlot, MaxStableProcessPlot
+from experiment.robustness_plot.single_plot import SinglePlot
+from spatio_temporal_dataset.coordinates.spatial_coordinates.alps_station_2D_coordinates import \
+    AlpsStation2DCoordinatesBetweenZeroAndOne, AlpsStationCoordinatesBetweenZeroAndTwo
+from spatio_temporal_dataset.coordinates.spatial_coordinates.generated_coordinates import CircleCoordinatesRadius1, \
+    CircleCoordinatesRadius2
+
+
+# def single_spatial_robustness_alps():
+#     spatial_robustness = SingleMspSpatial(grid_row_item=SingleMspSpatial.NbObservationItem,
+#                                           grid_column_item=SingleMspSpatial.SpatialCoordinateClassItem,
+#                                           plot_row_item=SingleMspSpatial.NbStationItem,
+#                                           plot_label_item=SingleMspSpatial.MaxStableModelItem)
+#     # Put only the parameter that will vary
+#     spatial_robustness.robustness_grid_plot(**{
+#         SingleMspSpatial.NbStationItem.name: list(range(43, 87, 15)),
+#         SingleMspSpatial.NbObservationItem.name: [10],
+#         SingleMspSpatial.MaxStableModelItem.name: [Smith(), BrownResnick()][:],
+#         SingleMspSpatial.SpatialCoordinateClassItem.name: [CircleCoordinatesRadius1,
+#                                                            AlpsStationCoordinatesBetweenZeroAndOne][:],
+#     })
+
+
+def multiple_spatial_robustness_alps():
+    nb_observation = 20
+    nb_sample = 1
+    plot_name = 'fast_result'
+    nb_stations = list(range(43, 87, 15))
+    # nb_stations = [10, 20, 30]
+
+    spatial_robustness = MultipleMaxStableProcessPlot(
+        grid_column_item=MaxStableProcessPlot.CoordinateClassItem,
+        plot_row_item=MaxStableProcessPlot.NbStationItem,
+        plot_label_item=MaxStableProcessPlot.MaxStableModelItem,
+        nb_samples=nb_sample,
+        main_title="Max stable analysis with {} years of temporal_observations".format(nb_observation),
+        plot_png_filename=plot_name
+    )
+    # Load all the models
+    msp_models = [Smith(), BrownResnick()]
+    # for covariance_function in CovarianceFunction:
+    #     msp_models.extend([ExtremalT(covariance_function=covariance_function)])
+
+    # Put only the parameter that will vary
+    spatial_robustness.robustness_grid_plot(**{
+        SinglePlot.OrdinateItem.name: [AbstractEstimator.MAE_ERROR, AbstractEstimator.DURATION],
+        MaxStableProcessPlot.NbStationItem.name: nb_stations,
+        MaxStableProcessPlot.NbObservationItem.name: nb_observation,
+        MaxStableProcessPlot.MaxStableModelItem.name: msp_models,
+        MaxStableProcessPlot.CoordinateClassItem.name: [CircleCoordinatesRadius1,
+                                                        CircleCoordinatesRadius2,
+                                                        AlpsStation2DCoordinatesBetweenZeroAndOne,
+                                                        AlpsStationCoordinatesBetweenZeroAndTwo][:],
+    })
+
+
+if __name__ == '__main__':
+    # single_spatial_robustness_alps()
+    multiple_spatial_robustness_alps()
diff --git a/experiment/robustness_plot/estimation_robustness/unidimensional_robustness/__init__.py b/experiment/robustness_plot/estimation_robustness/unidimensional_robustness/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/experiment/robustness_plot/estimation_robustness/unidimensional_robustness/unidimensional_robustness.py b/experiment/robustness_plot/estimation_robustness/unidimensional_robustness/unidimensional_robustness.py
new file mode 100644
index 0000000000000000000000000000000000000000..846127f1612a2d325ea506e8eb6b30548228ac30
--- /dev/null
+++ b/experiment/robustness_plot/estimation_robustness/unidimensional_robustness/unidimensional_robustness.py
@@ -0,0 +1,61 @@
+from extreme_estimator.estimator.abstract_estimator import AbstractEstimator
+from extreme_estimator.extreme_models.max_stable_model.max_stable_models import Smith, BrownResnick
+from experiment.robustness_plot.estimation_robustness.max_stable_process_plot import MultipleMaxStableProcessPlot, MaxStableProcessPlot
+from experiment.robustness_plot.single_plot import SinglePlot
+from spatio_temporal_dataset.coordinates.spatial_coordinates.alps_station_2D_coordinates import \
+    AlpsStation2DCoordinatesBetweenZeroAndOne, AlpsStationCoordinatesBetweenZeroAndTwo
+from spatio_temporal_dataset.coordinates.spatial_coordinates.generated_coordinates import CircleCoordinatesRadius1, \
+    CircleCoordinatesRadius2
+
+
+# def single_spatial_robustness_alps():
+#     spatial_robustness = SingleMspSpatial(grid_row_item=SingleMspSpatial.NbObservationItem,
+#                                           grid_column_item=SingleMspSpatial.SpatialCoordinateClassItem,
+#                                           plot_row_item=SingleMspSpatial.NbStationItem,
+#                                           plot_label_item=SingleMspSpatial.MaxStableModelItem)
+#     # Put only the parameter that will vary
+#     spatial_robustness.robustness_grid_plot(**{
+#         SingleMspSpatial.NbStationItem.name: list(range(43, 87, 15)),
+#         SingleMspSpatial.NbObservationItem.name: [10],
+#         SingleMspSpatial.MaxStableModelItem.name: [Smith(), BrownResnick()][:],
+#         SingleMspSpatial.SpatialCoordinateClassItem.name: [CircleCoordinatesRadius1,
+#                                                            AlpsStationCoordinatesBetweenZeroAndOne][:],
+#     })
+
+
+def multiple_unidimensional_robustness():
+    nb_observation = 20
+    nb_sample = 1
+    plot_name = 'fast_result'
+    nb_stations = list(range(43, 87, 15))
+    # nb_stations = [10, 20, 30]
+
+    spatial_robustness = MultipleMaxStableProcessPlot(
+        grid_column_item=MaxStableProcessPlot.CoordinateClassItem,
+        plot_row_item=MaxStableProcessPlot.NbStationItem,
+        plot_label_item=MaxStableProcessPlot.MaxStableModelItem,
+        nb_samples=nb_sample,
+        main_title="Max stable analysis with {} years of temporal_observations".format(nb_observation),
+        plot_png_filename=plot_name
+    )
+    # Load all the models
+    msp_models = [Smith(), BrownResnick()]
+    # for covariance_function in CovarianceFunction:
+    #     msp_models.extend([ExtremalT(covariance_function=covariance_function)])
+
+    # Put only the parameter that will vary
+    spatial_robustness.robustness_grid_plot(**{
+        SinglePlot.OrdinateItem.name: [AbstractEstimator.MAE_ERROR, AbstractEstimator.DURATION],
+        MaxStableProcessPlot.NbStationItem.name: nb_stations,
+        MaxStableProcessPlot.NbObservationItem.name: nb_observation,
+        MaxStableProcessPlot.MaxStableModelItem.name: msp_models,
+        MaxStableProcessPlot.CoordinateClassItem.name: [CircleCoordinatesRadius1,
+                                                        CircleCoordinatesRadius2,
+                                                        AlpsStation2DCoordinatesBetweenZeroAndOne,
+                                                        AlpsStationCoordinatesBetweenZeroAndTwo][:],
+    })
+
+
+if __name__ == '__main__':
+    # single_spatial_robustness_alps()
+    multiple_unidimensional_robustness()
diff --git a/extreme_estimator/robustness_plot/multiple_plot.py b/experiment/robustness_plot/multiple_plot.py
similarity index 96%
rename from extreme_estimator/robustness_plot/multiple_plot.py
rename to experiment/robustness_plot/multiple_plot.py
index ae98c0c09875d7808c63124dcdefcabc24b2d53a..e0b733f2597d64cc576af9ab5dc2ed2aec0535b7 100644
--- a/extreme_estimator/robustness_plot/multiple_plot.py
+++ b/experiment/robustness_plot/multiple_plot.py
@@ -1,4 +1,4 @@
-from extreme_estimator.robustness_plot.single_plot import SinglePlot
+from experiment.robustness_plot.single_plot import SinglePlot
 
 
 class MultiplePlot(SinglePlot):
diff --git a/extreme_estimator/robustness_plot/single_plot.py b/experiment/robustness_plot/single_plot.py
similarity index 98%
rename from extreme_estimator/robustness_plot/single_plot.py
rename to experiment/robustness_plot/single_plot.py
index 06e6f2e25597128a90fef861ef3ecf37b8e4e253..9dee3b785ee25def4145e0988a5d0b7144432505 100644
--- a/extreme_estimator/robustness_plot/single_plot.py
+++ b/experiment/robustness_plot/single_plot.py
@@ -1,13 +1,12 @@
 import os
 import os.path as op
-import random
 
 import matplotlib.pyplot as plt
 import numpy as np
 from itertools import product
 
 from extreme_estimator.estimator.abstract_estimator import AbstractEstimator
-from extreme_estimator.robustness_plot.display_item import DisplayItem
+from experiment.robustness_plot.display_item import DisplayItem
 from utils import get_full_path
 
 plt.style.use('seaborn-white')
diff --git a/extreme_estimator/estimator/abstract_estimator.py b/extreme_estimator/estimator/abstract_estimator.py
index 7f705474ae76fc3d5989302d5f60f51f83eec0a8..3e22c0e27b5eeb7bd75d418861e8adb3265ae0f9 100644
--- a/extreme_estimator/estimator/abstract_estimator.py
+++ b/extreme_estimator/estimator/abstract_estimator.py
@@ -12,7 +12,7 @@ class AbstractEstimator(object):
     # - The optimization method for each part of the process
 
     def __init__(self, dataset: AbstractDataset):
-        self.dataset = dataset
+        self.dataset = dataset  # type: AbstractDataset
         self.additional_information = dict()
 
     def fit(self):
diff --git a/extreme_estimator/estimator/full_estimator.py b/extreme_estimator/estimator/full_estimator.py
index 8c6fea1d509d273261ccdd9d97b1ea705193cb60..04fe3ad9b13d5618e9ab7a285eb5857db278b333 100644
--- a/extreme_estimator/estimator/full_estimator.py
+++ b/extreme_estimator/estimator/full_estimator.py
@@ -25,7 +25,7 @@ class SmoothMarginalsThenUnitaryMsp(AbstractFullEstimator):
         self.margin_estimator.fit()
         # Compute the maxima_frech
         maxima_frech = AbstractMarginModel.gev2frech(maxima_gev=self.dataset.maxima_gev,
-                                                     coordinates=self.dataset.coordinates,
+                                                     coordinates=self.dataset.coordinates_values,
                                                      margin_function=self.margin_estimator.margin_function_fitted)
         # Update maxima frech field through the dataset object
         self.dataset.maxima_frech = maxima_frech
diff --git a/extreme_estimator/estimator/margin_estimator.py b/extreme_estimator/estimator/margin_estimator.py
index 82fdc39fbdf040d24a41d22cf992e4d294acf3f0..ce4d9a3ea54130e6be4078f4de210b2cc4a7fcc7 100644
--- a/extreme_estimator/estimator/margin_estimator.py
+++ b/extreme_estimator/estimator/margin_estimator.py
@@ -33,4 +33,4 @@ class SmoothMarginEstimator(AbstractMarginEstimator):
 
     def _fit(self):
         self._margin_function_fitted = self.margin_model.fitmargin_from_maxima_gev(maxima_gev=self.dataset.maxima_gev,
-                                                                                   coordinates=self.dataset.coordinates)
+                                                                                   coordinates=self.dataset.coordinates_values)
diff --git a/extreme_estimator/estimator/max_stable_estimator.py b/extreme_estimator/estimator/max_stable_estimator.py
index 4ca92d18ee48b44a5bee8d581326b78d19a231ea..d83809a3cd2a3c8107732710f234707688829cda 100644
--- a/extreme_estimator/estimator/max_stable_estimator.py
+++ b/extreme_estimator/estimator/max_stable_estimator.py
@@ -19,8 +19,7 @@ class MaxStableEstimator(AbstractMaxStableEstimator):
         assert self.dataset.maxima_frech is not None
         self.max_stable_params_fitted = self.max_stable_model.fitmaxstab(
             maxima_frech=self.dataset.maxima_frech,
-            df_coordinates=self.dataset.spatial_coordinates.df_coordinates)
-            # coord=self.dataset.coordinates)
+            df_coordinates=self.dataset.coordinates.df_coordinates)
 
     def _error(self, true_max_stable_params: dict):
         absolute_errors = {param_name: np.abs(param_true_value - self.max_stable_params_fitted[param_name])
diff --git a/extreme_estimator/extreme_models/margin_model/abstract_margin_model.py b/extreme_estimator/extreme_models/margin_model/abstract_margin_model.py
index 46d04cecb7ee211bb349cfc2bb282ed18e9610f6..5ca3d45bb6af8b1a1c698ae476aa5e51bb787cb2 100644
--- a/extreme_estimator/extreme_models/margin_model/abstract_margin_model.py
+++ b/extreme_estimator/extreme_models/margin_model/abstract_margin_model.py
@@ -3,23 +3,23 @@ import numpy as np
 from extreme_estimator.extreme_models.abstract_model import AbstractModel
 from extreme_estimator.extreme_models.margin_model.margin_function.abstract_margin_function import AbstractMarginFunction
 from extreme_estimator.gev_params import GevParams
-from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractSpatialCoordinates
+from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
 
 
 class AbstractMarginModel(AbstractModel):
 
-    def __init__(self, spatial_coordinates: AbstractSpatialCoordinates, params_start_fit=None, params_sample=None):
+    def __init__(self, coordinates: AbstractCoordinates, params_start_fit=None, params_sample=None):
         super().__init__(params_start_fit, params_sample)
-        self.spatial_coordinates = spatial_coordinates
+        self.coordinates = coordinates
         self.margin_function_sample = None  # type: AbstractMarginFunction
         self.margin_function_start_fit = None  # type: AbstractMarginFunction
         self.load_margin_functions()
 
     def load_margin_functions(self, margin_function_class: type = None):
         assert margin_function_class is not None
-        self.margin_function_sample = margin_function_class(spatial_coordinates=self.spatial_coordinates,
+        self.margin_function_sample = margin_function_class(coordinates=self.coordinates,
                                                             default_params=GevParams.from_dict(self.params_sample))
-        self.margin_function_start_fit = margin_function_class(spatial_coordinates=self.spatial_coordinates,
+        self.margin_function_start_fit = margin_function_class(coordinates=self.coordinates,
                                                                default_params=GevParams.from_dict(
                                                                    self.params_start_fit))
 
diff --git a/extreme_estimator/extreme_models/margin_model/margin_function/abstract_margin_function.py b/extreme_estimator/extreme_models/margin_model/margin_function/abstract_margin_function.py
index 865aee3042147b897168c58d5bf5ebd062fced91..8b768b35e197f40bbf1ea08bf51d2b78960b8ef5 100644
--- a/extreme_estimator/extreme_models/margin_model/margin_function/abstract_margin_function.py
+++ b/extreme_estimator/extreme_models/margin_model/margin_function/abstract_margin_function.py
@@ -3,14 +3,14 @@ import matplotlib.pyplot as plt
 import numpy as np
 
 from extreme_estimator.gev_params import GevParams
-from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractSpatialCoordinates
+from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
 
 
 class AbstractMarginFunction(object):
     """ Class of function mapping points from a space S (could be 1D, 2D,...) to R^3 (the 3 parameters of the GEV)"""
 
-    def __init__(self, spatial_coordinates: AbstractSpatialCoordinates, default_params: GevParams):
-        self.spatial_coordinates = spatial_coordinates
+    def __init__(self, coordinates: AbstractCoordinates, default_params: GevParams):
+        self.coordinates = coordinates
         self.default_params = default_params.to_dict()
 
     def get_gev_params(self, coordinate: np.ndarray) -> GevParams:
@@ -20,8 +20,8 @@ class AbstractMarginFunction(object):
     # Visualization function
 
     def visualize_2D(self, gev_param_name=GevParams.GEV_LOC, ax=None, show=False):
-        x = self.spatial_coordinates.x_coordinates
-        y = self.spatial_coordinates.y_coordinates
+        x = self.coordinates.x_coordinates
+        y = self.coordinates.y_coordinates
         grid = self.get_grid_2D(x, y)
         gev_param_idx = GevParams.GEV_PARAM_NAMES.index(gev_param_name)
         if ax is None:
diff --git a/extreme_estimator/extreme_models/margin_model/margin_function/independent_margin_function.py b/extreme_estimator/extreme_models/margin_model/margin_function/independent_margin_function.py
index dc7d27c53c7e6745127ee55fac776af7aa2535b5..65adfcf3f2984998269b9934bba434a5f741e1e3 100644
--- a/extreme_estimator/extreme_models/margin_model/margin_function/independent_margin_function.py
+++ b/extreme_estimator/extreme_models/margin_model/margin_function/independent_margin_function.py
@@ -7,15 +7,15 @@ from extreme_estimator.extreme_models.margin_model.margin_function.param_functio
 from extreme_estimator.gev_params import GevParams
 from extreme_estimator.extreme_models.margin_model.margin_function.abstract_margin_function import \
     AbstractMarginFunction
-from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractSpatialCoordinates
+from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
 
 
 class IndependentMarginFunction(AbstractMarginFunction):
     """Margin Function where each parameter of the GEV are modeled independently"""
 
-    def __init__(self, spatial_coordinates: AbstractSpatialCoordinates, default_params: GevParams):
+    def __init__(self, coordinates: AbstractCoordinates, default_params: GevParams):
         """Attribute 'gev_param_name_to_param_function' maps each GEV parameter to its corresponding function"""
-        super().__init__(spatial_coordinates, default_params)
+        super().__init__(coordinates, default_params)
         self.gev_param_name_to_param_function = None  # type: Dict[str, ParamFunction]
 
     def get_gev_params(self, coordinate: np.ndarray) -> GevParams:
@@ -34,7 +34,7 @@ class LinearMarginFunction(IndependentMarginFunction):
         On the extremal point along all the dimension, the GEV parameters will equal the default_param value
         Then, it will augment linearly along a single linear axis"""
 
-    def __init__(self, spatial_coordinates: AbstractSpatialCoordinates,
+    def __init__(self, coordinates: AbstractCoordinates,
                  default_params: GevParams,
                  gev_param_name_to_linear_axes: Dict[str, List[int]],
                  gev_param_name_and_axis_to_start_fit: Dict[Tuple[str, int], float] = None):
@@ -42,13 +42,13 @@ class LinearMarginFunction(IndependentMarginFunction):
         -Attribute 'gev_param_name_to_linear_axis'        maps each GEV parameter to its corresponding function
         -Attribute 'gev_param_name_and_axis_to_start_fit' maps each tuple (GEV parameter, axis) to its start value for
             fitting (by default equal to 1). Also start value for the intercept is equal to 0 by default."""
-        super().__init__(spatial_coordinates, default_params)
+        super().__init__(coordinates, default_params)
         self.gev_param_name_and_axis_to_start_fit = gev_param_name_and_axis_to_start_fit
         self.gev_param_name_to_linear_axes = gev_param_name_to_linear_axes
 
         # Check the axes are well-defined with respect to the coordinates
         for axes in self.gev_param_name_to_linear_axes.values():
-            assert all([axis < np.ndim(spatial_coordinates.coordinates) for axis in axes])
+            assert all([axis < np.ndim(coordinates.coordinates_values) for axis in axes])
 
         # Build gev_parameter_to_param_function dictionary
         self.gev_param_name_to_param_function = {}  # type: Dict[str, ParamFunction]
@@ -60,7 +60,7 @@ class LinearMarginFunction(IndependentMarginFunction):
             # Otherwise, we fit a LinearParamFunction
             else:
                 param_function = LinearParamFunction(linear_axes=self.gev_param_name_to_linear_axes[gev_param_name],
-                                                     coordinates=self.spatial_coordinates.coordinates,
+                                                     coordinates=self.coordinates.coordinates_values,
                                                      start=self.default_params[gev_param_name])
                 # Some check on the Linear param function
                 if gev_param_name == GevParams.GEV_SCALE:
@@ -79,7 +79,7 @@ class LinearMarginFunction(IndependentMarginFunction):
         :return:
         """
         fit_marge_form_dict = {}
-        axis_to_name = {i: name for i, name in enumerate(AbstractSpatialCoordinates.COORDINATE_NAMES)}
+        axis_to_name = {i: name for i, name in enumerate(AbstractCoordinates.COORDINATE_NAMES)}
         for gev_param_name in GevParams.GEV_PARAM_NAMES:
             axes = self.gev_param_name_to_linear_axes.get(gev_param_name, None)
             formula_str = '1' if axes is None else '+'.join([axis_to_name[axis] for axis in axes])
diff --git a/extreme_estimator/extreme_models/margin_model/smooth_margin_model.py b/extreme_estimator/extreme_models/margin_model/smooth_margin_model.py
index aeec0112747f237c0dd81054479fc06354090d71..15f16c562075e129a262eec5a9e89dd9504d2275 100644
--- a/extreme_estimator/extreme_models/margin_model/smooth_margin_model.py
+++ b/extreme_estimator/extreme_models/margin_model/smooth_margin_model.py
@@ -13,10 +13,10 @@ class LinearMarginModel(AbstractMarginModel):
     def load_margin_functions(self, gev_param_name_to_linear_axis=None):
         self.default_params_sample = GevParams(1.0, 1.0, 1.0).to_dict()
         self.default_params_start_fit = GevParams(1.0, 1.0, 1.0).to_dict()
-        self.margin_function_sample = LinearMarginFunction(spatial_coordinates=self.spatial_coordinates,
+        self.margin_function_sample = LinearMarginFunction(coordinates=self.coordinates,
                                                            default_params=GevParams.from_dict(self.params_sample),
                                                            gev_param_name_to_linear_axes=gev_param_name_to_linear_axis)
-        self.margin_function_start_fit = LinearMarginFunction(spatial_coordinates=self.spatial_coordinates,
+        self.margin_function_start_fit = LinearMarginFunction(coordinates=self.coordinates,
                                                               default_params=GevParams.from_dict(self.params_start_fit),
                                                               gev_param_name_to_linear_axes=gev_param_name_to_linear_axis)
 
diff --git a/extreme_estimator/robustness_plot/msp_robustness.py b/extreme_estimator/robustness_plot/msp_robustness.py
deleted file mode 100644
index dc03ce853d4777337176b7a74ca9db661172cee5..0000000000000000000000000000000000000000
--- a/extreme_estimator/robustness_plot/msp_robustness.py
+++ /dev/null
@@ -1,112 +0,0 @@
-from extreme_estimator.extreme_models.max_stable_model.abstract_max_stable_model import AbstractMaxStableModel
-from extreme_estimator.extreme_models.max_stable_model.max_stable_models import Smith, BrownResnick
-from extreme_estimator.estimator.abstract_estimator import AbstractEstimator
-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 MaxStableDataset
-from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractSpatialCoordinates
-from spatio_temporal_dataset.coordinates.spatial_coordinates.alps_station_2D_coordinates import \
-    AlpsStation2DCoordinatesBetweenZeroAndOne, AlpsStationCoordinatesBetweenZeroAndTwo
-from spatio_temporal_dataset.coordinates.spatial_coordinates.generated_coordinates import CircleCoordinatesRadius1, \
-    CircleCoordinatesRadius2
-from extreme_estimator.robustness_plot.display_item import DisplayItem
-
-
-class MaxStableDisplayItem(DisplayItem):
-
-    def display_name_from_value(self, value: AbstractMaxStableModel):
-        return value.cov_mod
-
-
-class SpatialCoordinateDisplayItem(DisplayItem):
-
-    def display_name_from_value(self, value: AbstractSpatialCoordinates):
-        return str(value).split('.')[-1].split("'")[0]
-
-
-class MspSpatial(object):
-    MaxStableModelItem = MaxStableDisplayItem('max_stable_model', Smith)
-    SpatialCoordinateClassItem = SpatialCoordinateDisplayItem('spatial_coordinate_class', CircleCoordinatesRadius1)
-    NbStationItem = DisplayItem('Number of stations', 50)
-    NbObservationItem = DisplayItem('nb_obs', 60)
-
-    def msp_spatial_ordinates(self, **kwargs_single_point) -> dict:
-        # Get the argument from kwargs
-        max_stable_model = self.MaxStableModelItem.value_from_kwargs(
-            **kwargs_single_point)  # type: AbstractMaxStableModel
-        spatial_coordinate_class = self.SpatialCoordinateClassItem.value_from_kwargs(**kwargs_single_point)
-        nb_station = self.NbStationItem.value_from_kwargs(**kwargs_single_point)
-        nb_obs = self.NbObservationItem.value_from_kwargs(**kwargs_single_point)
-        # Run the estimation
-        spatial_coordinates = spatial_coordinate_class.from_nb_points(nb_points=nb_station)
-        dataset = MaxStableDataset.from_sampling(nb_obs=nb_obs, max_stable_model=max_stable_model,
-                                                 spatial_coordinates=spatial_coordinates)
-        estimator = MaxStableEstimator(dataset, max_stable_model)
-        estimator.fit()
-        return estimator.scalars(max_stable_model.params_sample)
-
-
-class SingleMspSpatial(SinglePlot, MspSpatial):
-
-    def compute_value_from_kwargs_single_point(self, **kwargs_single_point):
-        return self.msp_spatial_ordinates(**kwargs_single_point)
-
-
-class MultipleMspSpatial(MultiplePlot, MspSpatial):
-
-    def compute_value_from_kwargs_single_point(self, **kwargs_single_point):
-        return self.msp_spatial_ordinates(**kwargs_single_point)
-
-
-# def single_spatial_robustness_alps():
-#     spatial_robustness = SingleMspSpatial(grid_row_item=SingleMspSpatial.NbObservationItem,
-#                                           grid_column_item=SingleMspSpatial.SpatialCoordinateClassItem,
-#                                           plot_row_item=SingleMspSpatial.NbStationItem,
-#                                           plot_label_item=SingleMspSpatial.MaxStableModelItem)
-#     # Put only the parameter that will vary
-#     spatial_robustness.robustness_grid_plot(**{
-#         SingleMspSpatial.NbStationItem.name: list(range(43, 87, 15)),
-#         SingleMspSpatial.NbObservationItem.name: [10],
-#         SingleMspSpatial.MaxStableModelItem.name: [Smith(), BrownResnick()][:],
-#         SingleMspSpatial.SpatialCoordinateClassItem.name: [CircleCoordinatesRadius1,
-#                                                            AlpsStationCoordinatesBetweenZeroAndOne][:],
-#     })
-
-
-def multiple_spatial_robustness_alps():
-    nb_observation = 60
-    nb_sample = 10
-    plot_name = 'fast_result'
-    nb_stations = list(range(43, 87, 15))
-    # nb_stations = [10, 20, 30]
-
-    spatial_robustness = MultipleMspSpatial(
-        grid_column_item=MspSpatial.SpatialCoordinateClassItem,
-        plot_row_item=MspSpatial.NbStationItem,
-        plot_label_item=MspSpatial.MaxStableModelItem,
-        nb_samples=nb_sample,
-        main_title="Max stable analysis with {} years of observations".format(nb_observation),
-        plot_png_filename=plot_name
-    )
-    # Load all the models
-    msp_models = [Smith(), BrownResnick()]
-    # for covariance_function in CovarianceFunction:
-    #     msp_models.extend([ExtremalT(covariance_function=covariance_function)])
-
-    # Put only the parameter that will vary
-    spatial_robustness.robustness_grid_plot(**{
-        SinglePlot.OrdinateItem.name: [AbstractEstimator.MAE_ERROR, AbstractEstimator.DURATION],
-        MspSpatial.NbStationItem.name: nb_stations,
-        MspSpatial.NbObservationItem.name: nb_observation,
-        MspSpatial.MaxStableModelItem.name: msp_models,
-        MspSpatial.SpatialCoordinateClassItem.name: [CircleCoordinatesRadius1,
-                                                     CircleCoordinatesRadius2,
-                                                     AlpsStation2DCoordinatesBetweenZeroAndOne,
-                                                     AlpsStationCoordinatesBetweenZeroAndTwo][:],
-    })
-
-
-if __name__ == '__main__':
-    # single_spatial_robustness_alps()
-    multiple_spatial_robustness_alps()
diff --git a/spatio_temporal_dataset/coordinates/abstract_coordinates.py b/spatio_temporal_dataset/coordinates/abstract_coordinates.py
index 3b6d31beafa52645f15b67f1a51e85732629019b..a51e8ccb615941cce3086bfce26ed0793e73cf49 100644
--- a/spatio_temporal_dataset/coordinates/abstract_coordinates.py
+++ b/spatio_temporal_dataset/coordinates/abstract_coordinates.py
@@ -7,7 +7,7 @@ import matplotlib.pyplot as plt
 from mpl_toolkits.mplot3d import Axes3D
 
 
-class AbstractSpatialCoordinates(object):
+class AbstractCoordinates(object):
     # Columns
     COORDINATE_X = 'coord_x'
     COORDINATE_Y = 'coord_y'
@@ -58,7 +58,7 @@ class AbstractSpatialCoordinates(object):
     @classmethod
     def from_nb_points(cls, nb_points: int, **kwargs):
         # Call the default class method from csv
-        coordinates = cls.from_csv()  # type: AbstractSpatialCoordinates
+        coordinates = cls.from_csv()  # type: AbstractCoordinates
         # Sample randomly nb_points coordinates
         nb_coordinates = len(coordinates)
         if nb_points > nb_coordinates:
@@ -72,12 +72,12 @@ class AbstractSpatialCoordinates(object):
         ind = self.s_split == split_str
         return self.df_coordinates.loc[ind]
 
-    def coordinates_values(self, df_coordinates: pd.DataFrame) -> np.ndarray:
+    def _coordinates_values(self, df_coordinates: pd.DataFrame) -> np.ndarray:
         return df_coordinates.loc[:, self.coordinates_columns(df_coordinates)].values
 
     @property
-    def coordinates(self) -> np.ndarray:
-        return self.coordinates_values(df_coordinates=self.df_coordinates)
+    def coordinates_values(self) -> np.ndarray:
+        return self._coordinates_values(df_coordinates=self.df_coordinates)
 
     @property
     def x_coordinates(self) -> np.ndarray:
@@ -89,11 +89,11 @@ class AbstractSpatialCoordinates(object):
 
     @property
     def coordinates_train(self) -> np.ndarray:
-        return self.coordinates_values(df_coordinates=self.df_coordinates_split(self.TRAIN_SPLIT_STR))
+        return self._coordinates_values(df_coordinates=self.df_coordinates_split(self.TRAIN_SPLIT_STR))
 
     @property
     def coordinates_test(self) -> np.ndarray:
-        return self.coordinates_values(df_coordinates=self.df_coordinates_split(self.TEST_SPLIT_STR))
+        return self._coordinates_values(df_coordinates=self.df_coordinates_split(self.TEST_SPLIT_STR))
 
     @property
     def index(self):
@@ -103,20 +103,20 @@ class AbstractSpatialCoordinates(object):
 
     def visualization_1D(self):
         assert len(self.coordinates_columns(self.df_coordinates)) >= 1
-        x = self.coordinates[:]
+        x = self.coordinates_values[:]
         y = np.zeros(len(x))
         plt.scatter(x, y)
         plt.show()
 
     def visualization_2D(self):
         assert len(self.coordinates_columns(self.df_coordinates)) >= 2
-        x, y = self.coordinates[:, 0], self.coordinates[:, 1]
+        x, y = self.coordinates_values[:, 0], self.coordinates_values[:, 1]
         plt.scatter(x, y)
         plt.show()
 
     def visualization_3D(self):
         assert len(self.coordinates_columns(self.df_coordinates)) == 3
-        x, y, z = self.coordinates[:, 0], self.coordinates[:, 1], self.coordinates[:, 2]
+        x, y, z = self.coordinates_values[:, 0], self.coordinates_values[:, 1], self.coordinates_values[:, 2]
         fig = plt.figure()
         ax = fig.add_subplot(111, projection='3d')  # type: Axes3D
         ax.scatter(x, y, z, marker='^')
diff --git a/spatio_temporal_dataset/coordinates/spatial_coordinates/alps_station_2D_coordinates.py b/spatio_temporal_dataset/coordinates/spatial_coordinates/alps_station_2D_coordinates.py
index db3d3156d730b3aac0655b58090364b70a660d3b..629e24f94824abbef266000689f7da7ea55bf922 100644
--- a/spatio_temporal_dataset/coordinates/spatial_coordinates/alps_station_2D_coordinates.py
+++ b/spatio_temporal_dataset/coordinates/spatial_coordinates/alps_station_2D_coordinates.py
@@ -19,7 +19,7 @@ class AlpsStation2DCoordinatesBetweenZeroAndOne(AlpsStation2DCoordinates):
     @classmethod
     def from_csv(cls, csv_file='coord-lambert2'):
         coord = super().from_csv(csv_file)
-        return TransformedCoordinates.from_coordinates(spatial_coordinates=coord,
+        return TransformedCoordinates.from_coordinates(coordinates=coord,
                                                        transformation_function=BetweenZeroAndOne2DNormalization())
 
 
diff --git a/spatio_temporal_dataset/coordinates/spatial_coordinates/alps_station_3D_coordinates.py b/spatio_temporal_dataset/coordinates/spatial_coordinates/alps_station_3D_coordinates.py
index 02a8062ca42a66b69c2a1684ae3a0346d0e961cc..d5ef2ad10d125cda6db89b5e5469032c02cd6cb7 100644
--- a/spatio_temporal_dataset/coordinates/spatial_coordinates/alps_station_3D_coordinates.py
+++ b/spatio_temporal_dataset/coordinates/spatial_coordinates/alps_station_3D_coordinates.py
@@ -2,13 +2,13 @@ import os.path as op
 
 import pandas as pd
 
-from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractSpatialCoordinates
+from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
 from spatio_temporal_dataset.coordinates.transformed_coordinates.tranformation_3D import AnisotropyTransformation
 from spatio_temporal_dataset.coordinates.transformed_coordinates.transformed_coordinates import TransformedCoordinates
 from utils import get_full_path
 
 
-class AlpsStation3DCoordinates(AbstractSpatialCoordinates):
+class AlpsStation3DCoordinates(AbstractCoordinates):
     """
     For the Alps Stations, X, Y coordinates are in Lambert 2. Altitude is in meters
     """
@@ -43,5 +43,5 @@ class AlpsStation3DCoordinatesWithAnisotropy(AlpsStation3DCoordinates):
     @classmethod
     def from_csv(cls, csv_file='coord-lambert2'):
         coord = super().from_csv(csv_file)
-        return TransformedCoordinates.from_coordinates(spatial_coordinates=coord,
+        return TransformedCoordinates.from_coordinates(coordinates=coord,
                                                        transformation_function=AnisotropyTransformation())
diff --git a/spatio_temporal_dataset/coordinates/spatial_coordinates/generated_coordinates.py b/spatio_temporal_dataset/coordinates/spatial_coordinates/generated_coordinates.py
index 883504562f99934783b9e08ece41d7711a48eded..e9d939afd2f4dc503a3448aa552610ef3bf88e81 100644
--- a/spatio_temporal_dataset/coordinates/spatial_coordinates/generated_coordinates.py
+++ b/spatio_temporal_dataset/coordinates/spatial_coordinates/generated_coordinates.py
@@ -3,11 +3,11 @@ import numpy as np
 import pandas as pd
 
 from extreme_estimator.extreme_models.utils import get_loaded_r
-from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractSpatialCoordinates
+from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
 import matplotlib.pyplot as plt
 
 
-class CircleCoordinatesRadius1(AbstractSpatialCoordinates):
+class CircleCoordinatesRadius1(AbstractCoordinates):
 
     @classmethod
     def from_nb_points(cls, nb_points, max_radius=1.0):
diff --git a/spatio_temporal_dataset/coordinates/transformed_coordinates/transformed_coordinates.py b/spatio_temporal_dataset/coordinates/transformed_coordinates/transformed_coordinates.py
index b1e4c7a6b1441008f70a45f2b065d25a7fbfd21f..4cc302d98d3854f889274d14f626927c3e10996e 100644
--- a/spatio_temporal_dataset/coordinates/transformed_coordinates/transformed_coordinates.py
+++ b/spatio_temporal_dataset/coordinates/transformed_coordinates/transformed_coordinates.py
@@ -1,14 +1,14 @@
-from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractSpatialCoordinates
+from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
 from spatio_temporal_dataset.coordinates.transformed_coordinates.abstract_transformation import AbstractTransformation
 
 
-class TransformedCoordinates(AbstractSpatialCoordinates):
+class TransformedCoordinates(AbstractCoordinates):
 
     @classmethod
-    def from_coordinates(cls, spatial_coordinates: AbstractSpatialCoordinates,
+    def from_coordinates(cls, coordinates: AbstractCoordinates,
                          transformation_function: AbstractTransformation):
-        df_coordinates_transformed = spatial_coordinates.df_coordinates.copy()
+        df_coordinates_transformed = coordinates.df_coordinates.copy()
         df_coordinates_transformed = transformation_function.transform(df_coord=df_coordinates_transformed)
-        return cls(df_coordinates=df_coordinates_transformed, s_split=spatial_coordinates.s_split)
+        return cls(df_coordinates=df_coordinates_transformed, s_split=coordinates.s_split)
 
 
diff --git a/spatio_temporal_dataset/coordinates/unidimensional_coordinates/__init__.py b/spatio_temporal_dataset/coordinates/unidimensional_coordinates/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/spatio_temporal_dataset/coordinates/axis_coordinates/axis_coordinates.py b/spatio_temporal_dataset/coordinates/unidimensional_coordinates/axis_coordinates.py
similarity index 79%
rename from spatio_temporal_dataset/coordinates/axis_coordinates/axis_coordinates.py
rename to spatio_temporal_dataset/coordinates/unidimensional_coordinates/axis_coordinates.py
index 21153b1e66ad8c9dad592b8d07cf3be21d4e236c..d0541f86fb53e9edba43972c6279c02d929e700a 100644
--- a/spatio_temporal_dataset/coordinates/axis_coordinates/axis_coordinates.py
+++ b/spatio_temporal_dataset/coordinates/unidimensional_coordinates/axis_coordinates.py
@@ -3,14 +3,14 @@ import pandas as pd
 import numpy as np
 
 from extreme_estimator.extreme_models.utils import get_loaded_r
-from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractSpatialCoordinates
+from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
 
 
-class AxisCoordinates(AbstractSpatialCoordinates):
+class UniDimensionalCoordinates(AbstractCoordinates):
     pass
 
 
-class UniformAxisCoordinates(AxisCoordinates):
+class UniformAxisCoordinates(UniDimensionalCoordinates):
 
     @classmethod
     def from_nb_points(cls, nb_points, start=0.0, end=1.0):
diff --git a/spatio_temporal_dataset/dataset/abstract_dataset.py b/spatio_temporal_dataset/dataset/abstract_dataset.py
index aa46eb827430e65acf6531a1b3ab1baa33856dff..f42453efab75cb81dd6022e294feb88f7bcaa876 100644
--- a/spatio_temporal_dataset/dataset/abstract_dataset.py
+++ b/spatio_temporal_dataset/dataset/abstract_dataset.py
@@ -3,24 +3,24 @@ import numpy as np
 import os.path as op
 import pandas as pd
 from spatio_temporal_dataset.temporal_observations.abstract_temporal_observations import AbstractTemporalObservations
-from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractSpatialCoordinates
+from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
 
 
 class AbstractDataset(object):
 
-    def __init__(self, temporal_observations: AbstractTemporalObservations, spatial_coordinates: AbstractSpatialCoordinates):
+    def __init__(self, temporal_observations: AbstractTemporalObservations, coordinates: AbstractCoordinates):
         # is_same_index = temporal_observations.index == coordinates.index  # type: pd.Series
         # assert is_same_index.all()
         self.temporal_observations = temporal_observations
-        self.spatial_coordinates = spatial_coordinates
+        self.coordinates = coordinates
 
     @classmethod
     def from_csv(cls, csv_path: str):
         assert op.exists(csv_path)
         df = pd.read_csv(csv_path)
         temporal_maxima = AbstractTemporalObservations.from_df(df)
-        spatial_coordinates = AbstractSpatialCoordinates.from_df(df)
-        return cls(temporal_maxima, spatial_coordinates)
+        coordinates = AbstractCoordinates.from_df(df)
+        return cls(temporal_maxima, coordinates)
 
     def to_csv(self, csv_path: str):
         dirname = op.dirname(csv_path)
@@ -31,15 +31,15 @@ class AbstractDataset(object):
     @property
     def df_dataset(self) -> pd.DataFrame:
         # Merge dataframes with the maxima and with the coordinates
-        return self.temporal_observations.df_maxima_gev.join(self.spatial_coordinates.df_coordinates)
+        return self.temporal_observations.df_maxima_gev.join(self.coordinates.df_coordinates)
 
     @property
     def df_coordinates(self):
-        return self.spatial_coordinates.df_coordinates
+        return self.coordinates.df_coordinates
 
     @property
-    def coordinates(self):
-        return self.spatial_coordinates.coordinates
+    def coordinates_values(self):
+        return self.coordinates.coordinates_values
 
     @property
     def maxima_gev(self) -> np.ndarray:
diff --git a/spatio_temporal_dataset/dataset/simulation_dataset.py b/spatio_temporal_dataset/dataset/simulation_dataset.py
index e40d77b1ea98e2ed683306397a2d70cead7e7794..0b489d3a9f527775aa1441984ae2b65b567a07d5 100644
--- a/spatio_temporal_dataset/dataset/simulation_dataset.py
+++ b/spatio_temporal_dataset/dataset/simulation_dataset.py
@@ -1,7 +1,7 @@
 from extreme_estimator.extreme_models.margin_model.abstract_margin_model import AbstractMarginModel
 from extreme_estimator.extreme_models.max_stable_model.abstract_max_stable_model import AbstractMaxStableModel
 from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
-from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractSpatialCoordinates
+from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
 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
@@ -14,10 +14,10 @@ class SimulatedDataset(AbstractDataset):
     """
 
     def __init__(self, temporal_observations: AbstractTemporalObservations,
-                 spatial_coordinates: AbstractSpatialCoordinates,
+                 coordinates: AbstractCoordinates,
                  max_stable_model: AbstractMaxStableModel = None,
                  margin_model: AbstractMarginModel = None):
-        super().__init__(temporal_observations, spatial_coordinates)
+        super().__init__(temporal_observations, 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
@@ -26,31 +26,25 @@ class SimulatedDataset(AbstractDataset):
 class MaxStableDataset(SimulatedDataset):
 
     @classmethod
-    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)
+    def from_sampling(cls, nb_obs: int, max_stable_model: AbstractMaxStableModel, coordinates: AbstractCoordinates):
+        temporal_obs = MaxStableAnnualMaxima.from_sampling(nb_obs, max_stable_model, coordinates)
+        return cls(temporal_observations=temporal_obs, coordinates=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)
+    def from_sampling(cls, nb_obs: int, margin_model: AbstractMarginModel,coordinates: AbstractCoordinates):
+        temporal_obs = MarginAnnualMaxima.from_sampling(nb_obs, coordinates, margin_model)
+        return cls(temporal_observations=temporal_obs, coordinates=coordinates, margin_model=margin_model)
 
 
 class FullSimulatedDataset(SimulatedDataset):
 
     @classmethod
     def from_double_sampling(cls, nb_obs: int, max_stable_model: AbstractMaxStableModel,
-                             spatial_coordinates: AbstractSpatialCoordinates,
+                             coordinates: AbstractCoordinates,
                              margin_model: AbstractMarginModel):
         temporal_obs = FullAnnualMaxima.from_double_sampling(nb_obs, max_stable_model,
-                                                             spatial_coordinates, margin_model)
-        return cls(temporal_obs, spatial_coordinates, max_stable_model)
+                                                             coordinates, margin_model)
+        return cls(temporal_obs, coordinates, max_stable_model)
diff --git a/spatio_temporal_dataset/marginals/abstract_marginals.py b/spatio_temporal_dataset/marginals/abstract_marginals.py
deleted file mode 100644
index 63a9b49d475e18c34ac041d99bc19d489422f3f0..0000000000000000000000000000000000000000
--- a/spatio_temporal_dataset/marginals/abstract_marginals.py
+++ /dev/null
@@ -1,48 +0,0 @@
-from typing import List
-
-from R.gev_fit.gev_marginal import GevMarginal, frechet_unitary_transformation
-from spatio_temporal_dataset.stations.station import Station
-
-
-class AbstractMarginals(object):
-
-    def __init__(self, stations: List[Station]):
-        self.stations = stations
-        self.gev_marginals = []  # type: List[GevMarginal]
-
-        #  Compute some temporal arguments
-        self.years_of_study = self.stations[0].year_of_study
-        self.temporal_window_size = 20
-
-
-        self.smoothing_function = None
-
-    def sample_marginals(self, smoothing_function):
-        pass
-
-    @property
-    def gev_parameters(self):
-        return [gev_marginal.gev_parameters_estimated for gev_marginal in self.gev_marginals]
-
-
-class FrechetUnitaryTransformationFunction(object):
-    pass
-
-class NoSmoothing(object):
-
-    def gev_parameters(self, coordinate):
-        pass
-
-class ContinuousSmoothing(object):
-
-    def frechet_unitary_transformation(self, data, coordinate):
-        gev_parameters = self.gev_parameters(coordinate)
-        transformed_data = frechet_unitary_transformation(data, gev_parameters)
-
-
-    def gev_parameters(self, coordinate):
-        return
-        return ()
-
-if __name__ == '__main__':
-    pass
\ No newline at end of file
diff --git a/spatio_temporal_dataset/marginals/spatial_marginals.py b/spatio_temporal_dataset/marginals/spatial_marginals.py
deleted file mode 100644
index 6530e0cdfd9ecf89a201c3f1454df54dcead4964..0000000000000000000000000000000000000000
--- a/spatio_temporal_dataset/marginals/spatial_marginals.py
+++ /dev/null
@@ -1,14 +0,0 @@
-from spatio_temporal_dataset.marginals.abstract_marginals import AbstractMarginals
-
-from spatio_temporal_dataset.coordinates.abstract_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
-
-
-class SimulatedSpatialMarginal(SpatialMarginal):
-    pass
diff --git a/spatio_temporal_dataset/spatio_temporal_data_handler.py b/spatio_temporal_dataset/spatio_temporal_data_handler.py
deleted file mode 100644
index b94df33b796673b636ec98dbe9833b99f07bff15..0000000000000000000000000000000000000000
--- a/spatio_temporal_dataset/spatio_temporal_data_handler.py
+++ /dev/null
@@ -1,59 +0,0 @@
-from typing import List
-import pandas as pd
-from spatio_temporal_dataset.marginals.abstract_marginals import AbstractMarginals
-from spatio_temporal_dataset.marginals.spatial_marginals import SpatialMarginal
-from spatio_temporal_dataset.stations.station import Station, load_stations_from_dataframe
-from spatio_temporal_dataset.stations.station_distance import EuclideanDistance2D, StationDistance
-import pickle
-import os.path as op
-from itertools import combinations
-
-
-class SpatioTemporalDataHandler(object):
-
-    def __init__(self, marginals_class: type, stations: List[Station], station_distance: StationDistance):
-        self.stations = stations
-
-        #  Compute once the distances between stations
-        for station1, station2 in combinations(self.stations, 2):
-            distance = station_distance.compute_distance(station1, station2)
-            station1.distance[station2] = distance
-            station2.distance[station1] = distance
-
-        # Compute the marginals
-        self.marginals = marginals_class(self.stations)  # type: AbstractMarginals
-
-        # Define the max stable
-        # self.max_stable =
-
-        print(self.marginals.gev_parameters)
-
-    @classmethod
-    def from_dataframe(cls, df):
-        return cls.from_spatial_dataframe(df)
-
-    @classmethod
-    def from_spatial_dataframe(cls, df):
-        stations = load_stations_from_dataframe(df)
-        marginal_class = SpatialMarginal
-        station_distance = EuclideanDistance2D()
-        return cls(marginals_class=marginal_class, stations=stations, station_distance=station_distance)
-
-
-def get_spatio_temporal_data_handler(pickle_path: str, load_pickle: bool = True, dump_pickle: bool = False, *args) \
-        -> SpatioTemporalDataHandler:
-    # Either load or dump pickle of a SpatioTemporalDataHandler object
-    assert load_pickle or dump_pickle
-    if load_pickle:
-        assert op.exists(pickle_path) and not dump_pickle
-        spatio_temporal_experiment = pickle.load(pickle_path)
-    else:
-        assert not op.exists(pickle_path)
-        spatio_temporal_experiment = SpatioTemporalDataHandler(*args)
-        pickle.dump(spatio_temporal_experiment, file=pickle_path)
-    return spatio_temporal_experiment
-
-
-if __name__ == '__main__':
-    df = pd.DataFrame(1, index=['station1', 'station2'], columns=['200' + str(i) for i in range(18)])
-    xp = SpatioTemporalDataHandler.from_dataframe(df)
diff --git a/spatio_temporal_dataset/stations/station.py b/spatio_temporal_dataset/stations/station.py
deleted file mode 100644
index 9e3f6aa1123350edb5548916c35c7b8a695174ca..0000000000000000000000000000000000000000
--- a/spatio_temporal_dataset/stations/station.py
+++ /dev/null
@@ -1,29 +0,0 @@
-import pandas as pd
-import numpy as np
-
-
-class Station(object):
-    def __init__(self, name: str, annual_maxima: pd.Series, longitude=np.nan, latitude=np.nan, altitude=np.nan):
-        self.annual_maxima = annual_maxima
-        self.year_of_study = list(annual_maxima.index)
-        self.name = name
-        self.altitude = altitude
-        self.latitude = latitude
-        self.longitude = longitude
-        self.distance = {}
-
-
-def load_stations_from_dataframe(df):
-    return [Station(name=i, annual_maxima=row) for i, row in df.iterrows()]
-
-def load_station_from_two_dataframe(df, location_df):
-    pass
-
-
-if __name__ == '__main__':
-    df = pd.DataFrame(1, index=['station1', 'station2'], columns=['200' + str(i) for i in range(9)])
-    stations = load_stations_from_dataframe(df)
-    station = stations[0]
-    print(station.name)
-    print(station.annual_maxima)
-    print(station.year_of_study)
diff --git a/spatio_temporal_dataset/stations/station_distance.py b/spatio_temporal_dataset/stations/station_distance.py
deleted file mode 100644
index d4a40eb5ef918d278e4eb5483b11a2fec93c2019..0000000000000000000000000000000000000000
--- a/spatio_temporal_dataset/stations/station_distance.py
+++ /dev/null
@@ -1,22 +0,0 @@
-from spatio_temporal_dataset.stations.station import Station
-import numpy as np
-
-
-class StationDistance(object):
-
-    @classmethod
-    def compute_distance(self, station1: Station, station2: Station) -> float:
-        return np.nan
-
-
-def euclidean_distance(arr1: np.ndarray, arr2: np.ndarray) -> float:
-    return np.linalg.norm(arr1 - arr2)
-
-
-class EuclideanDistance2D(StationDistance):
-
-    @classmethod
-    def compute_distance(self, station1: Station, station2: Station) -> float:
-        print(station1.latitude)
-        stations_coordinates = [np.array([station.latitude, station.longitude]) for station in [station1, station2]]
-        return euclidean_distance(*stations_coordinates)
diff --git a/spatio_temporal_dataset/temporal_observations/annual_maxima_observations.py b/spatio_temporal_dataset/temporal_observations/annual_maxima_observations.py
index 5ee31b26b0c26d24daba48bb11aaf8d73c928d6d..6eb582b71e2d27b9555e8e4e483c27891a07626c 100644
--- a/spatio_temporal_dataset/temporal_observations/annual_maxima_observations.py
+++ b/spatio_temporal_dataset/temporal_observations/annual_maxima_observations.py
@@ -2,7 +2,7 @@ import pandas as pd
 
 from extreme_estimator.extreme_models.margin_model.abstract_margin_model import AbstractMarginModel
 from extreme_estimator.extreme_models.max_stable_model.abstract_max_stable_model import AbstractMaxStableModel
-from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractSpatialCoordinates
+from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
 from spatio_temporal_dataset.temporal_observations.abstract_temporal_observations import AbstractTemporalObservations
 
 
@@ -17,20 +17,19 @@ class AnnualMaxima(AbstractTemporalObservations):
 class MarginAnnualMaxima(AnnualMaxima):
 
     @classmethod
-    def from_sampling(cls, nb_obs: int, spatial_coordinates: AbstractSpatialCoordinates,
+    def from_sampling(cls, nb_obs: int, coordinates: AbstractCoordinates,
                       margin_model: AbstractMarginModel):
-        maxima_gev = margin_model.rmargin_from_nb_obs(nb_obs=nb_obs, coordinates=spatial_coordinates.coordinates)
-        df_maxima_gev = pd.DataFrame(data=maxima_gev, index=spatial_coordinates.index)
+        maxima_gev = margin_model.rmargin_from_nb_obs(nb_obs=nb_obs, coordinates=coordinates.coordinates_values)
+        df_maxima_gev = pd.DataFrame(data=maxima_gev, index=coordinates.index)
         return cls(df_maxima_gev=df_maxima_gev)
 
 
 class MaxStableAnnualMaxima(AbstractTemporalObservations):
 
     @classmethod
-    def from_sampling(cls, nb_obs: int, max_stable_model: AbstractMaxStableModel,
-                      spatial_coordinates: AbstractSpatialCoordinates):
-        maxima_frech = max_stable_model.rmaxstab(nb_obs=nb_obs, coordinates=spatial_coordinates.coordinates)
-        df_maxima_frech = pd.DataFrame(data=maxima_frech, index=spatial_coordinates.index)
+    def from_sampling(cls, nb_obs: int, max_stable_model: AbstractMaxStableModel, coordinates: AbstractCoordinates):
+        maxima_frech = max_stable_model.rmaxstab(nb_obs=nb_obs, coordinates=coordinates.coordinates_values)
+        df_maxima_frech = pd.DataFrame(data=maxima_frech, index=coordinates.index)
         return cls(df_maxima_frech=df_maxima_frech)
 
 
@@ -38,11 +37,10 @@ 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)
+                             coordinates: AbstractCoordinates, margin_model: AbstractMarginModel):
+        max_stable_annual_maxima = super().from_sampling(nb_obs, max_stable_model, coordinates)
         #  Compute df_maxima_gev from df_maxima_frech
         maxima_gev = margin_model.rmargin_from_maxima_frech(maxima_frech=max_stable_annual_maxima.maxima_frech,
-                                                            coordinates=spatial_coordinates.coordinates)
-        max_stable_annual_maxima.df_maxima_gev = pd.DataFrame(data=maxima_gev, index=spatial_coordinates.index)
+                                                            coordinates=coordinates.coordinates_values)
+        max_stable_annual_maxima.df_maxima_gev = pd.DataFrame(data=maxima_gev, index=coordinates.index)
         return max_stable_annual_maxima
diff --git a/test/test_extreme_estimator/test_R_model/test_margin_function.py b/test/test_extreme_estimator/test_R_model/test_margin_function.py
index 12b099e34bffb71f20ea8ff69eed0572f92d4f1b..e789085e46f6745e3c11b4680584ffa98e3857e9 100644
--- a/test/test_extreme_estimator/test_R_model/test_margin_function.py
+++ b/test/test_extreme_estimator/test_R_model/test_margin_function.py
@@ -10,7 +10,7 @@ class TestLinearMarginModel(unittest.TestCase):
 
     def test_visualization_2D(self):
         spatial_coordinates = CircleCoordinatesRadius1.from_nb_points(nb_points=50)
-        margin_model = LinearShapeAxis0MarginModel(spatial_coordinates=spatial_coordinates)
+        margin_model = LinearShapeAxis0MarginModel(coordinates=spatial_coordinates)
         for gev_param_name in GevParams.GEV_PARAM_NAMES:
             margin_model.margin_function_sample.visualize_2D(gev_param_name=gev_param_name, show=self.DISPLAY)
         # maxima_gev = margin_model.rmargin_from_nb_obs(nb_obs=10, coordinates=coordinates)
diff --git a/test/test_extreme_estimator/test_estimator/test_full_estimators.py b/test/test_extreme_estimator/test_estimator/test_full_estimators.py
index d112bb6172b68219e9b11b08117e7999ce0b9f2f..a5a51dfd74b6c5d5d7e7ea058c020788e06ba31f 100644
--- a/test/test_extreme_estimator/test_estimator/test_full_estimators.py
+++ b/test/test_extreme_estimator/test_estimator/test_full_estimators.py
@@ -22,7 +22,7 @@ class TestFullEstimators(unittest.TestCase):
     def test_full_estimators(self):
         for margin_model, max_stable_model in product(self.smooth_margin_models, self.max_stable_models):
             dataset = FullSimulatedDataset.from_double_sampling(nb_obs=10, margin_model=margin_model,
-                                                                spatial_coordinates=self.spatial_coordinates,
+                                                                coordinates=self.spatial_coordinates,
                                                                 max_stable_model=max_stable_model)
 
             for estimator_class in self.FULL_ESTIMATORS:
diff --git a/test/test_extreme_estimator/test_estimator/test_margin_estimators.py b/test/test_extreme_estimator/test_estimator/test_margin_estimators.py
index c3a9b54304aaf9bbbf0addbe65b692acd9c9f1b8..89cfc17f7d8b41c0b32d2cef28c68def7d04be4f 100644
--- a/test/test_extreme_estimator/test_estimator/test_margin_estimators.py
+++ b/test/test_extreme_estimator/test_estimator/test_margin_estimators.py
@@ -23,12 +23,12 @@ class TestSmoothMarginEstimator(unittest.TestCase):
 
     @classmethod
     def load_smooth_margin_models(cls, spatial_coordinates):
-        return [margin_class(spatial_coordinates=spatial_coordinates) for margin_class in cls.MARGIN_TYPES]
+        return [margin_class(coordinates=spatial_coordinates) for margin_class in cls.MARGIN_TYPES]
 
     def test_dependency_estimators(self):
         for margin_model in self.smooth_margin_models:
             dataset = MarginDataset.from_sampling(nb_obs=10, margin_model=margin_model,
-                                                  spatial_coordinates=self.spatial_coordinates)
+                                                  coordinates=self.spatial_coordinates)
             # Fit estimator
             estimator = SmoothMarginEstimator(dataset=dataset, margin_model=margin_model)
             estimator.fit()
diff --git a/test/test_extreme_estimator/test_estimator/test_max_stable_estimators.py b/test/test_extreme_estimator/test_estimator/test_max_stable_estimators.py
index 050256d6fa9fa4556125ac1cb9563be9b6e1eedf..6eb78d409d371bf4a81e449af9e77b62dce3e627 100644
--- a/test/test_extreme_estimator/test_estimator/test_max_stable_estimators.py
+++ b/test/test_extreme_estimator/test_estimator/test_max_stable_estimators.py
@@ -35,7 +35,7 @@ class TestMaxStableEstimators(unittest.TestCase):
         for max_stable_model in self.max_stable_models:
             dataset = MaxStableDataset.from_sampling(nb_obs=10,
                                                      max_stable_model=max_stable_model,
-                                                     spatial_coordinates=self.spatial_coord)
+                                                     coordinates=self.spatial_coord)
 
             for estimator_class in self.MAX_STABLE_ESTIMATORS:
                 estimator = estimator_class(dataset=dataset, max_stable_model=max_stable_model)
diff --git a/test/test_spatio_temporal_dataset/test_coordinates.py b/test/test_spatio_temporal_dataset/test_coordinates.py
index b4890ee21ced949f01c68b18ea6aabfa5bb6b7af..b2bf236483d9da8083e1a843ac1b85b19f88ec9d 100644
--- a/test/test_spatio_temporal_dataset/test_coordinates.py
+++ b/test/test_spatio_temporal_dataset/test_coordinates.py
@@ -1,6 +1,6 @@
 import unittest
 
-from spatio_temporal_dataset.coordinates.axis_coordinates.axis_coordinates import UniformAxisCoordinates
+from spatio_temporal_dataset.coordinates.unidimensional_coordinates.axis_coordinates import UniformAxisCoordinates
 from spatio_temporal_dataset.coordinates.spatial_coordinates.alps_station_2D_coordinates import \
     AlpsStation2DCoordinatesBetweenZeroAndOne
 from spatio_temporal_dataset.coordinates.spatial_coordinates.alps_station_3D_coordinates import \
@@ -30,7 +30,7 @@ class TestSpatialCoordinates(unittest.TestCase):
         self.assertTrue(True)
 
 
-class TestAxisCoordinates(unittest.TestCase):
+class TestUniDimensionalCoordinates(unittest.TestCase):
     DISPLAY = False
 
     def test_unif(self):
diff --git a/test/test_spatio_temporal_dataset/test_dataset.py b/test/test_spatio_temporal_dataset/test_dataset.py
new file mode 100644
index 0000000000000000000000000000000000000000..320d68778beec576269c7576185453071cced40b
--- /dev/null
+++ b/test/test_spatio_temporal_dataset/test_dataset.py
@@ -0,0 +1,11 @@
+import unittest
+
+
+class TestDataset(unittest.TestCase):
+    pass
+    # already done by the estimator somehow,
+    # maybe I could add some more precise test that check for 1D, 2D and 3D Datasets
+
+
+if __name__ == '__main__':
+    unittest.main()