diff --git a/extreme_estimator/extreme_models/max_stable_model/abstract_max_stable_model.py b/extreme_estimator/extreme_models/max_stable_model/abstract_max_stable_model.py
index 029eebc99e50cbdf6f92c52ed1c620a23c165cbf..57cd069e2974791b8bc26b9eb5ed2dabb29fbb22 100644
--- a/extreme_estimator/extreme_models/max_stable_model/abstract_max_stable_model.py
+++ b/extreme_estimator/extreme_models/max_stable_model/abstract_max_stable_model.py
@@ -2,11 +2,13 @@ from enum import Enum
 
 import numpy as np
 import pandas as pd
+from rpy2.rinterface import RRuntimeWarning
+from rpy2.rinterface._rinterface import RRuntimeError
 
 from extreme_estimator.extreme_models.abstract_model import AbstractModel
 from extreme_estimator.extreme_models.result_from_fit import ResultFromFit, ResultFromSpatialExtreme
 from extreme_estimator.extreme_models.utils import r, safe_run_r_estimator, get_coord, \
-    get_margin_formula
+    get_margin_formula, SafeRunException
 from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
 
 
@@ -71,7 +73,8 @@ class AbstractMaxStableModel(AbstractModel):
                                    **fit_params)
         return ResultFromSpatialExtreme(res)
 
-    def rmaxstab(self, nb_obs: int, coordinates_values: np.ndarray, use_rmaxstab_with_2_coordinates=False) -> np.ndarray:
+    def rmaxstab(self, nb_obs: int, coordinates_values: np.ndarray,
+                 use_rmaxstab_with_2_coordinates=False) -> np.ndarray:
         """
         Return an numpy of maxima. With rows being the stations and columns being the years of maxima
         """
@@ -79,8 +82,13 @@ class AbstractMaxStableModel(AbstractModel):
             # When we have more than 2 coordinates, then sample based on the 2 first coordinates only
             coordinates_values = coordinates_values[:, :2]
 
-        maxima_frech = np.array(
-            r.rmaxstab(nb_obs, coordinates_values, *list(self.cov_mod_param.values()), **self.params_sample))
+        try:
+            maxima_frech = np.array(
+                r.rmaxstab(nb_obs, coordinates_values, *list(self.cov_mod_param.values()), **self.params_sample))
+        except (RRuntimeError, RRuntimeWarning) as e:
+            raise SafeRunException('\n\nSome R exception have been launched at RunTime: \n{} \n{}'.format(e.__repr__(),
+                                   'To sample from 3D data we advise to set the function argument '
+                                   '"use_rmaxstab_with_2_coordinates" to True'))
 
         return np.transpose(maxima_frech)
 
diff --git a/test/test_spatio_temporal_dataset/test_dataset.py b/test/test_spatio_temporal_dataset/test_dataset.py
index 340c360f6f77408326b6e7dce4e2700f96124683..538f9f0c7ea94dde9398fc02e72086f797e76e9b 100644
--- a/test/test_spatio_temporal_dataset/test_dataset.py
+++ b/test/test_spatio_temporal_dataset/test_dataset.py
@@ -5,7 +5,7 @@ import numpy as np
 from rpy2.rinterface import RRuntimeError
 
 from extreme_estimator.extreme_models.margin_model.linear_margin_model import LinearNonStationaryLocationMarginModel
-from extreme_estimator.extreme_models.utils import set_seed_for_test
+from extreme_estimator.extreme_models.utils import set_seed_for_test, SafeRunException
 from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
 from spatio_temporal_dataset.coordinates.transformed_coordinates.transformation.uniform_normalization import \
     BetweenZeroAndOneNormalization
@@ -30,7 +30,7 @@ class TestDataset(unittest.TestCase):
 
     def test_max_stable_dataset_crash_R3(self):
         """Test to warn me when spatialExtremes handles R3"""
-        with self.assertRaises(RRuntimeError):
+        with self.assertRaises(SafeRunException):
             smith_process = load_test_max_stable_models()[0]
             coordinates = load_test_3D_spatial_coordinates(nb_points=self.nb_points)[0]
             MaxStableDataset.from_sampling(nb_obs=self.nb_obs,
@@ -39,7 +39,9 @@ class TestDataset(unittest.TestCase):
 
     def test_max_stable_dataset_R3_with_R2_spatial_structure(self):
         """Test to create a dataset in R3, but where the spatial structure is only with respect to the 2 fisrt coordinates"""
-        coordinates = load_test_3D_spatial_coordinates(nb_points=self.nb_points, transformation_class=BetweenZeroAndOneNormalization)[0]
+        coordinates = \
+        load_test_3D_spatial_coordinates(nb_points=self.nb_points, transformation_class=BetweenZeroAndOneNormalization)[
+            0]
         smith_process = load_test_max_stable_models()[0]
         MaxStableDataset.from_sampling(nb_obs=self.nb_obs,
                                        max_stable_model=smith_process,
@@ -47,7 +49,6 @@ class TestDataset(unittest.TestCase):
                                        use_rmaxstab_with_2_coordinates=True)
 
 
-
 class TestSpatioTemporalDataset(unittest.TestCase):
     nb_obs = 2
     nb_points = 3