diff --git a/experiment/meteo_france_SCM_study/abstract_study.py b/experiment/meteo_france_SCM_study/abstract_study.py
index 12847bb18d8869bdf8d504b2da74abf83a64ef1e..0de8161269a6a1469886ff8450e04202bf0a2a02 100644
--- a/experiment/meteo_france_SCM_study/abstract_study.py
+++ b/experiment/meteo_france_SCM_study/abstract_study.py
@@ -9,6 +9,7 @@ from netCDF4 import Dataset
 
 from experiment.meteo_france_SCM_study.abstract_variable import AbstractVariable
 from experiment.meteo_france_SCM_study.massif import safran_massif_names_from_datasets
+from extreme_estimator.margin_fits.plot.mask_poly import mask_outside_polygon
 from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
 from spatio_temporal_dataset.coordinates.spatial_coordinates.abstract_spatial_coordinates import \
     AbstractSpatialCoordinates
@@ -130,15 +131,23 @@ class AbstractStudy(object):
                          row_massif[AbstractCoordinates.COORDINATE_Y])
                         for _, row_massif in df_massif.iterrows()]
 
-        for coordinate_id in set([tuple[0] for tuple in coord_tuples]):
+        for j, coordinate_id in enumerate(set([tuple[0] for tuple in coord_tuples])):
+            # Retrieve the list of coords (x,y) that define the contour of the massif of id coordinate_id
             l = [coords for idx, *coords in coord_tuples if idx == coordinate_id]
+            # if j == 0:
+            #     mask_outside_polygon(poly_verts=l, ax=ax)
+            # Plot the contour of the massif
             l = list(zip(*l))
             ax.plot(*l, color='black')
+            # Potentially, fill the inside of the polygon with some color
             if fill:
                 massif_name = self.coordinate_id_to_massif_name[coordinate_id]
                 fill_kwargs = massif_name_to_fill_kwargs[massif_name] if massif_name_to_fill_kwargs is not None else {}
                 ax.fill(*l, **fill_kwargs)
+        # Display the center of the massif
         ax.scatter(self.massifs_coordinates.x_coordinates, self.massifs_coordinates.y_coordinates, s=1)
+
+
         if show:
             plt.show()
 
diff --git a/experiment/meteo_france_SCM_study/main_visualize.py b/experiment/meteo_france_SCM_study/main_visualize.py
index 714fb36f925b07cd9bc2e468d6d1919411c75242..2d87709b9e28a0c8bcfd2631b92a9f65bff9c2a0 100644
--- a/experiment/meteo_france_SCM_study/main_visualize.py
+++ b/experiment/meteo_france_SCM_study/main_visualize.py
@@ -48,13 +48,13 @@ def extended_visualization():
 
 
 def normal_visualization():
-    save_to_file = True
-    only_first_one = False
-    for study_class in SCM_STUDIES[:]:
+    save_to_file = False
+    only_first_one = True
+    for study_class in SCM_STUDIES[-1:]:
         for study in study_iterator(study_class, only_first_one=only_first_one):
             study_visualizer = StudyVisualizer(study, save_to_file=save_to_file)
             # study_visualizer.visualize_independent_margin_fits(threshold=[None, 20, 40, 60][0])
-            study_visualizer.visualize_linear_margin_fit(only_first_max_stable=False)
+            study_visualizer.visualize_linear_margin_fit(only_first_max_stable=True)
 
 
 def complete_analysis(only_first_one=False):
diff --git a/extreme_estimator/margin_fits/plot/mask_poly.py b/extreme_estimator/margin_fits/plot/mask_poly.py
new file mode 100644
index 0000000000000000000000000000000000000000..d7a478680b54176a7444f73bae5455140eb196ae
--- /dev/null
+++ b/extreme_estimator/margin_fits/plot/mask_poly.py
@@ -0,0 +1,44 @@
+import matplotlib.pyplot as plt
+import matplotlib.patches as mpatches
+import matplotlib.path as mpath
+
+
+# from: https://stackoverflow.com/questions/3320311/fill-outside-of-polygon-mask-array-where-indicies-are-beyond-a-circular-bounda
+def mask_outside_polygon(poly_verts, ax=None):
+    """
+    Plots a mask on the specified axis ("ax", defaults to plt.gca()) such that
+    all areas outside of the polygon specified by "poly_verts" are masked.
+
+    "poly_verts" must be a list of tuples of the verticies in the polygon in
+    counter-clockwise order.
+
+    Returns the matplotlib.patches.PathPatch instance plotted on the figure.
+    """
+
+    if ax is None:
+        ax = plt.gca()
+
+    # Get current plot limits
+    xlim = ax.get_xlim()
+    ylim = ax.get_ylim()
+
+    # Verticies of the plot boundaries in clockwise order
+    bound_verts = [(xlim[0], ylim[0]), (xlim[0], ylim[1]),
+                   (xlim[1], ylim[1]), (xlim[1], ylim[0]),
+                   (xlim[0], ylim[0])]
+
+    # A series of codes (1 and 2) to tell matplotlib whether to draw a line or
+    # move the "pen" (So that there's no connecting line)
+    bound_codes = [mpath.Path.MOVETO] + (len(bound_verts) - 1) * [mpath.Path.LINETO]
+    poly_codes = [mpath.Path.MOVETO] + (len(poly_verts) - 1) * [mpath.Path.LINETO]
+
+    # Plot the masking patch
+    path = mpath.Path(bound_verts + poly_verts, bound_codes + poly_codes)
+    patch = mpatches.PathPatch(path, facecolor='white', edgecolor='none')
+    patch = ax.add_patch(patch)
+
+    # Reset the plot limits to their original extents
+    ax.set_xlim(xlim)
+    ax.set_ylim(ylim)
+
+    return patch