diff --git a/.gitignore b/.gitignore
index 519631c..1eeac1a 100644
--- a/.gitignore
+++ b/.gitignore
@@ -3,6 +3,8 @@ __pycache__/
 *.pyc
 *.pyd
 *.so
+*.cpp
+*.c
 
 # Setuptools/distutils distribution folder.
 dist/
diff --git a/seislib/__init__.py b/seislib/__init__.py
index 2c7cf26..b89ed28 100644
--- a/seislib/__init__.py
+++ b/seislib/__init__.py
@@ -15,6 +15,13 @@
 sys.path.insert(1, os.path.dirname(__file__))
 sys.path.insert(1, os.path.join(os.path.dirname(__file__), 'clib'))
 from .tomography.grid import EqualAreaGrid
+from . import an
+from . import colormaps
+from . import eq
+from . import exceptions
+from . import plotting
+from . import tomography
+from . import utils
 #from .tomography.tomography import SeismicTomography
 
 
diff --git a/seislib/an/an_attenuation.py b/seislib/an/an_attenuation.py
index 1c27c21..50e7de9 100644
--- a/seislib/an/an_attenuation.py
+++ b/seislib/an/an_attenuation.py
@@ -26,7 +26,7 @@
 from obspy.signal.invsim import cosine_taper
 from obspy.geodetics.base import gps2dist_azimuth
 from seislib import EqualAreaGrid
-from seislib.utils import load_pickle, save_pickle, resample
+from seislib.utils import load_pickle, save_pickle, remove_file, resample
 from seislib.plotting import add_earth_features
 from seislib.plotting import colormesh, contour, contourf
 from seislib.plotting import plot_stations, plot_map
@@ -70,7 +70,7 @@ class AmbientNoiseAttenuation:
     get_stations_coords(files)
         Retrieves the geographic coordinates associated with each receiver
         
-    prepare_data()
+    prepare_data(recompute=False)
         Retrieves and saves on disk the geographic coordinates
         
     parameterize(cell_size, overlap=0.5, min_no_stations=6, plotting=True)
@@ -96,8 +96,8 @@ class AmbientNoiseAttenuation:
         inversion
      
     plot_stations(ax=None, show=True, oceans_color='water', lands_color='land', 
-                  edgecolor='k', projection='Mercator', color_by_network=True, 
-                  legend_dict={}, **kwargs):
+                  edgecolor='k', projection='Mercator', resolution='110m', 
+                  color_by_network=True, legend_dict={}, **kwargs):
         Maps the seismic receivers for which data are available        
         
     
@@ -105,7 +105,7 @@ class AmbientNoiseAttenuation:
     -------------
     plot_map(mesh, c, ax=None, projection='Mercator', map_boundaries=None, 
              bound_map=True, colorbar=True, show=True, style='colormesh', 
-             add_features=False, cbar_dict={}, **kwargs)
+             add_features=False, resolution='110m', cbar_dict={}, **kwargs)
         Displays an attenuation map
     
     colormesh(mesh, c, ax, **kwargs):
@@ -336,7 +336,7 @@ def get_stations_coords(self, files):
         return coords
     
                 
-    def prepare_data(self):
+    def prepare_data(self, recompute=False):
         """ 
         Saves to disk the geographic coordinates associated with each receiver. 
         These are saved to $self.savedir/stations.pickle
@@ -350,9 +350,21 @@ def prepare_data(self):
               net1.sta2 : (lat2, lon2),
               net2.sta3 : (lat3, lon3)
               }
+            
+        
+        Parameters
+        ----------
+        recompute : bool
+            If True, the station coordinates and times will be removed from
+            disk and recalculated. Otherwise (default), if they are present,
+            they will be loaded into memory, avoiding any computation. This
+            parameter should be set to True whenever one has added files to
+            the source directory
         """
 
-        savefile = os.path.join(self.savedir, 'stations.pickle') 
+        savefile = os.path.join(self.savedir, 'stations.pickle')
+        if recompute:
+            remove_file(savefile)
         if not os.path.exists(savefile):
             coords = self.get_stations_coords(self.files)
             save_pickle(savefile, coords)
@@ -1357,7 +1369,7 @@ def get_related_pixels(mesh, lon1, lon2, lat1, lat2):
     @classmethod
     def plot_map(cls, mesh, c, ax=None, projection='Mercator', map_boundaries=None, 
                  bound_map=True, colorbar=True, show=True, style='colormesh', 
-                 add_features=False, cbar_dict={}, **kwargs):
+                 add_features=False, resolution='110m', cbar_dict={}, **kwargs):
         """ 
         Displays an attenuation map
         
@@ -1403,6 +1415,11 @@ def plot_map(cls, mesh, c, ax=None, projection='Mercator', map_boundaries=None,
             If True, natural Earth features will be added to the GeoAxesSubplot.
             Default is False. If `ax` is None, it is automatically set to True
             
+        resolution : str
+            Resolution of the Earth features displayed in the figure. Passed to
+            cartopy.feature.NaturalEarthFeature. Valid arguments are '110m',
+            '50m', '10m'. Default is '110m'
+            
         cbar_dict : dict
             Keyword arguments passed to matplotlib.colorbar.ColorbarBase
          
@@ -1430,6 +1447,7 @@ def plot_map(cls, mesh, c, ax=None, projection='Mercator', map_boundaries=None,
                         show=show, 
                         style=style, 
                         add_features=add_features, 
+                        resolution=resolution,
                         norm=norm,
                         cmap=cmap,
                         cbar_dict=cbar_dict,
@@ -1536,7 +1554,8 @@ def contour(cls, mesh, c, ax, smoothing=None, **kwargs):
     
     def plot_stations(self, ax=None, show=True, oceans_color='water', 
                       lands_color='land', edgecolor='k', projection='Mercator',
-                      color_by_network=True, legend_dict={}, **kwargs):
+                      resolution='110m', color_by_network=True, legend_dict={}, 
+                      **kwargs):
         """ Maps the seismic receivers for which data are available
         
         Parameters
@@ -1564,6 +1583,11 @@ def plot_stations(self, ax=None, show=True, oceans_color='water',
             Name of the geographic projection used to create the GeoAxesSubplot.
             (Visit the cartopy website for a list of valid projection names.)
             If ax is not None, `projection` is ignored. Default is 'Mercator'
+         
+        resolution : str
+            Resolution of the Earth features displayed in the figure. Passed to
+            cartopy.feature.NaturalEarthFeature. Valid arguments are '110m',
+            '50m', '10m'. Default is '110m'
             
         color_by_network : bool
             If True, each seismic network will have a different color in the
@@ -1589,6 +1613,7 @@ def plot_stations(self, ax=None, show=True, oceans_color='water',
                              lands_color=lands_color, 
                              edgecolor=edgecolor, 
                              projection=projection,
+                             resolution=resolution,
                              color_by_network=color_by_network, 
                              legend_dict=legend_dict,
                              **kwargs)
diff --git a/seislib/an/an_velocity.py b/seislib/an/an_velocity.py
index e3da1bc..a086f1d 100644
--- a/seislib/an/an_velocity.py
+++ b/seislib/an/an_velocity.py
@@ -13,7 +13,7 @@
 from obspy import read
 from obspy.geodetics.base import gps2dist_azimuth
 import matplotlib.pyplot as plt
-from seislib.utils import rotate_stream, load_pickle, save_pickle
+from seislib.utils import rotate_stream, load_pickle, save_pickle, remove_file
 from seislib.an import noisecorr, velocity_filter, extract_dispcurve
 from seislib.plotting import plot_stations
 
@@ -57,7 +57,7 @@ class AmbientNoiseVelocity:
         Retrieves the geographic coordinates and the starting and ending time
         associated with each continuous seismogram
         
-    prepare_data()
+    prepare_data(recompute=False)
         Saves to disk the geographic coordinates and the starting and ending 
         time associated with each continuous seismogram
         
@@ -76,8 +76,8 @@ class AmbientNoiseVelocity:
         at the specified frequency   
         
     plot_stations(ax=None, show=True, oceans_color='water', lands_color='land', 
-                  edgecolor='k', projection='Mercator', color_by_network=True, 
-                  legend_dict={}, **kwargs):
+                  edgecolor='k', projection='Mercator', resolution='110m',
+                  color_by_network=True, legend_dict={}, **kwargs)
         Maps the seismic receivers for which data are available        
         
         
@@ -318,7 +318,7 @@ def get_times_and_coords(self, files):
         return times, coords
 
             
-    def prepare_data(self):
+    def prepare_data(self, recompute=False):
         """ 
         Saves to disk the geographic coordinates and the starting and ending 
         time associated with each continuous seismogram. These are saved to
@@ -353,11 +353,24 @@ def prepare_data(self):
         considered in the calculation of phase velocities both the horizontal
         components (N and E, e.g. BHN and BHE) should be available. These will
         be rotated so as to analyse the transverse (T, for Love) and radial (R,
-        for radially-polarized Rayleigh waves) components.
+        for radially-polarized Rayleigh waves) components
+        
+        
+        Parameters
+        ----------
+        recompute : bool
+            If True, the station coordinates and times will be removed from
+            disk and recalculated. Otherwise (default), if they are present,
+            they will be loaded into memory, avoiding any computation. This
+            parameter should be set to True whenever one has added files to
+            the source directory
         """
         
         savecoords = os.path.join(self.savedir, 'stations.pickle')
         savetimes = os.path.join(self.savedir, 'timespans.pickle')
+        if recompute:
+            remove_file(savecoords)
+            remove_file(savetimes)
         if not os.path.exists(savecoords) or not os.path.exists(savetimes):
             times, coords = self.get_times_and_coords(self.files)
             save_pickle(savecoords, coords)
@@ -720,7 +733,8 @@ def display_progress(no_files, done, verbose=False):
     
     def plot_stations(self, ax=None, show=True, oceans_color='water', 
                       lands_color='land', edgecolor='k', projection='Mercator',
-                      color_by_network=True, legend_dict={}, **kwargs):
+                      resolution='110m', color_by_network=True, legend_dict={}, 
+                      **kwargs):
         """ Maps the seismic receivers for which data are available
         
         Parameters
@@ -749,6 +763,12 @@ def plot_stations(self, ax=None, show=True, oceans_color='water',
             (Visit the cartopy website for a list of valid projection names.)
             If ax is not None, `projection` is ignored. Default is 'Mercator'
             
+        resolution : str
+            Resolution of the Earth features displayed in the figure. Passed to
+            cartopy.feature.NaturalEarthFeature. Valid arguments are '110m',
+            '50m', '10m'. Default is '110m'
+
+            
         color_by_network : bool
             If True, each seismic network will have a different color in the
             resulting map, and a legend will be displayed. Otherwise, all
@@ -773,6 +793,7 @@ def plot_stations(self, ax=None, show=True, oceans_color='water',
                              lands_color=lands_color, 
                              edgecolor=edgecolor, 
                              projection=projection,
+                             resolution=resolution,
                              color_by_network=color_by_network, 
                              legend_dict=legend_dict,
                              **kwargs)      
diff --git a/seislib/clib/setup_all.py b/seislib/clib/setup_all.py
index ec6ff05..9ff2660 100644
--- a/seislib/clib/setup_all.py
+++ b/seislib/clib/setup_all.py
@@ -22,11 +22,27 @@
 from Cython.Build import cythonize
 import numpy
 import os
+import platform
+
 
 names = ['_math', '_spherical_geometry', '_tomography']
 sources = ['%s.pyx'%name for name in names]
 extra_compile_args = ["-O3", "-ffast-math", "-march=native", "-fopenmp" ]
-extra_link_args=['-fopenmp']
+
+platform_name = platform.system()
+
+
+if platform_name.lower() == 'darwin':
+    versions = os.listdir('/usr/local/Cellar/gcc/')
+    version = max(versions, key=lambda i: int(i.split('.')[0]))
+    version_int = version.split('.')[0]
+    path = '/usr/local/Cellar/gcc/%s/lib/gcc/%s'%(version, version_int)
+    os.environ['CC'] = 'gcc-%s'%version_int
+    os.environ['CXX'] = 'g++-%s'%version_int
+    extra_link_args=['-Wl,-rpath,%s'%path]
+
+else:
+    extra_link_args=['-fopenmp']
 
 
 for name, source in zip(names, sources):
diff --git a/seislib/eq/eq_velocity.py b/seislib/eq/eq_velocity.py
index 36c8d9e..dd308e5 100644
--- a/seislib/eq/eq_velocity.py
+++ b/seislib/eq/eq_velocity.py
@@ -16,13 +16,12 @@
 from scipy.interpolate import interp1d
 from scipy import signal
 import matplotlib.pyplot as plt
-import cartopy.crs as ccrs
 from obspy import read, read_events
 from obspy.geodetics.base import gps2dist_azimuth
 from seislib.utils import adapt_timespan, zeropad, bandpass_gaussian
 from seislib.utils import load_pickle, save_pickle, remove_file
 from seislib.utils import gaussian, skewed_normal, next_power_of_2
-from seislib.plotting import add_earth_features
+from seislib.plotting import plot_stations, plot_events
 from seislib.exceptions import DispersionCurveException
 
 
@@ -73,6 +72,10 @@ class EQVelocity:
         Writes to disk the geographic coordinates of the seismic receivers and 
         triplets of epicenters-receivers to be used for retrieving the
         dispersion curves
+       
+    get_events_used()
+        Retrieves the events id for which triplets of epicenter-receivers are
+        available to extract dispersion measurements
         
     delete_unused_files()    
         Deletes every file in the data directory which is not useful for 
@@ -95,12 +98,16 @@ class EQVelocity:
         at the specified period(s). (No extrapolation is made.)
         
     plot_stations(ax=None, show=True, oceans_color='water', lands_color='land', 
-                  edgecolor='k', projection='Mercator', color_by_network=True, 
-                  **kwargs):
+                  edgecolor='k', projection='Mercator', resolution='110m',
+                  color_by_network=True, legend_dict={}, **kwargs)
         Maps the seismic receivers for which data are available        
         
+    plot_events(ax=None, show=True, oceans_color='water', lands_color='land', 
+                edgecolor='k', projection='Mercator', resolution='110m', 
+                min_size=5, max_size=150, **kwargs):
+        Maps the seismic events used to obtain dispersion measurements
         
-    Class Methods
+        Class Methods
     -------------
     lie_on_same_gc(stla1, stlo1, stla2, stlo2, evla, evlo, azimuth_tolerance=5, 
                    distmin=None, distmax=None):
@@ -145,22 +152,29 @@ class EQVelocity:
                         distmax=2500, 
                         min_no_events=8)
     
-    which will extracts the geographic coordinates of each seismic receivers 
-    from the header of the sac files and prepare the triplets of epicenters-
-    receivers that will be used to obtain the dispersion curves.
+    which will extracts the geographic coordinates of each seismic receivers
+    from the header of the sac files, those of the seismic events, and prepare 
+    the triplets of epicenters-receivers that will be used to obtain the 
+    dispersion curves.
     
     ---------------------------------------------------------------------------        
-    NOTE: The coordinates are saved at /path/to/eq_velocity/Z/stations.pickle,
-    and are stored in a dictionary object where each key corresponds to a 
-    station code ($network_code.$station_code) and each value is a tuple 
-    containing latitude and longitude of the station. 
-    
-    For example:
+    NOTE: The geographic coordinates of the seismic receivers are saved at 
+    /$self.savedir/stations.pickle, and are stored in a dictionary object 
+    where each key corresponds to a station code ($network_code.$station_code) 
+    and each value is a tuple containing latitude and longitude of the station. 
         
+    For example:
+            
         { net1.sta1 : (lat1, lon1),
           net1.sta2 : (lat2, lon2),
           net2.sta3 : (lat3, lon3)
-          }
+            }
+            
+    The geographic coordinates of the events are saved at 
+    /$self.savedir/events.pickle, and have a similar structure to the above:
+    each key corresponds to an event origin time (in 
+    obspy.UTCDateTime.timestamp format), and each value is (lat, lon, mag) of 
+    the earthquake, where mag is magnitude.
         
     The triplets of epicenter-receivers are saved at 
     /path/to/eq_velocity/Z/triplets.pickle, and are stored in a dictionary object 
@@ -183,16 +197,20 @@ class EQVelocity:
                                     '1239825695.33']
           }
         
-    Note that each event in the list corresponds to a sub-directory of
-    the source data `src`. Also note that, as opposed to the above example, the 
+    Note that each event origin time corresponds to a sub-directory of
+    the source data (`src`). Also note that, as opposed to the above example, the 
     length of the lists of events will be at least 8, since `min_no_events` in 
     `prepare_data` was 8.
     ---------------------------------------------------------------------------
     
-    After the data have been "prepared", the receivers available can be 
-    displayed by typing
+    After the data have been "prepared", the receivers and events available can 
+    be displayed by typing
         
         eq.plot_stations()
+        
+    and
+    
+        eq.plot_events()
     
     Now we can calculate the dispersion curves automatically, provided that we
     pass a reference curve (i.e., ndarray of shape (n, 2), where the 1st column 
@@ -294,8 +312,8 @@ def __str__(self):
     def get_coords_and_triplets(self, events, azimuth_tolerance=5, distmin=None,
                                 distmax=None):
         """ 
-        Retrieves stations information and triplets of epicenter-receivers
-        to be used to calculate the phase velocities
+        Retrieves stations and events information, and the triplets of 
+        epicenter-receivers to be used to calculate the phase velocities
         
         Parameters
         ----------
@@ -321,7 +339,16 @@ def get_coords_and_triplets(self, events, azimuth_tolerance=5, distmin=None,
                   net1.sta2 : (lat2, lon2),
                   net2.sta3 : (lat3, lon3)
                   }
-                
+        
+        events_info : dict
+            each key corresponds to an event origin time and each value is a 
+            tuple containing latitude, longitude, and magnitude of the event.
+                    
+                    { '1222701570.84' : (lat1, lon1, mag1),
+                      '1237486660.74' : (lat2, lon2, mag2),
+                      '1239825695.33' : (lat3, lon3, mag3)
+                      }        
+                    
         triplets : dict
             each key is a tuple corresponding to a station pair and each value 
             is a list of events that are aligned on the same great circle path
@@ -348,13 +375,15 @@ def get_event_coords_from_xml(eventfile):
             event = read_events(eventfile)[0]
             evla = event.origins[0].latitude
             evlo = event.origins[0].longitude
-            return evla, evlo
+            mag = event.magnitudes[0].mag
+            return evla, evlo, mag
         
         def get_event_coords_from_sac(sacfile):
             tr = read(sacfile)[0]
             evla = tr.stats.sac.evla
             evlo = tr.stats.sac.evlo
-            return evla, evlo
+            mag = tr.stats.sac.mag
+            return evla, evlo, mag
         
         def get_station_coords(evdir, sacfile):
             nonlocal stations
@@ -368,6 +397,7 @@ def get_station_coords(evdir, sacfile):
         
         
         stations = {}
+        events_info = {}
         triplets = defaultdict(list)
         no_events = len(events)
         iprint = no_events//10 if no_events>10 else 1
@@ -380,11 +410,14 @@ def get_station_coords(evdir, sacfile):
             sacfiles = sorted([i for i in os.listdir(evdir) \
                                if i.endswith('sac') and i[-5]==self.component])
             try:
-                evla, evlo = get_event_coords_from_xml(os.path.join(evdir, 
-                                                                    '%s.xml'%event))
+                evla, evlo, mag = get_event_coords_from_xml(
+                        os.path.join(evdir, '%s.xml'%event)
+                        )
             except FileNotFoundError:
-                evla, evlo = get_event_coords_from_sac(os.path.join(evdir, 
-                                                                    sacfiles[0]))
+                evla, evlo, mag = get_event_coords_from_sac(
+                        os.path.join(evdir, sacfiles[0])
+                        )
+            events_info[event] = (evla, evlo, mag)
             for sac1, sac2 in it.combinations(sacfiles, 2):
                 sta1, (stla1, stlo1) = get_station_coords(evdir, sac1)
                 sta2, (stla2, stlo2) = get_station_coords(evdir, sac2)
@@ -398,7 +431,8 @@ def get_station_coords(evdir, sacfile):
                                        distmin=distmin,
                                        distmax=distmax):
                     triplets[(sta1, sta2)].append(event)
-        return stations, triplets
+                    
+        return stations, events_info, triplets
                 
                
     @classmethod
@@ -450,9 +484,9 @@ def lie_on_same_gc(cls, stla1, stlo1, stla2, stlo2, evla, evlo,
     def prepare_data(self, azimuth_tolerance=5, distmin=None, distmax=None, 
                      min_no_events=5, recompute=False, delete_unused_files=False):
         """
-        Saves to disk the geographic coordinates of the seismic receivers and 
-        triplets of epicenters-receivers to be used for retrieving the
-        dispersion curves.
+        Saves to disk the geographic coordinates of the seismic receivers and of
+        the seismic events, along with the triplets of epicenters-receivers to 
+        be used for retrieving the dispersion curves.
         
         Parameters
         ----------
@@ -490,10 +524,10 @@ def prepare_data(self, azimuth_tolerance=5, distmin=None, distmax=None,
             
         Note
         ----
-        The geographic coordinates are saved at /$self.savedir/stations.pickle,
-        and are stored in a dictionary object where each key corresponds to a 
-        station code ($network_code.$station_code) and each value is a tuple 
-        containing latitude and longitude of the station. 
+        The geographic coordinates of the seismic receivers are saved at 
+        /$self.savedir/stations.pickle, and are stored in a dictionary object 
+        where each key corresponds to a station code ($network_code.$station_code) 
+        and each value is a tuple containing latitude and longitude of the station. 
         
         For example:
             
@@ -502,6 +536,12 @@ def prepare_data(self, azimuth_tolerance=5, distmin=None, distmax=None,
               net2.sta3 : (lat3, lon3)
               }
             
+        The geographic coordinates of the events are saved at 
+        /$self.savedir/events.pickle, and have a similar structure to the above:
+        each key corresponds to an event origin time (in 
+        obspy.UTCDateTime.timestamp format), and each value is (lat, lon, mag) 
+        of the epicenter, where mag is the magnitude of the event.
+
         The triplets of epicenter-receivers are saved at 
         /$self.savedir/triplets.pickle, and are stored in a dictionary object 
         where each key is a tuple corresponding to a station pair and each value 
@@ -532,29 +572,59 @@ def prepare_data(self, azimuth_tolerance=5, distmin=None, distmax=None,
             of phase velocity, GJI
         """
                     
-        savecoords = os.path.join(self.savedir, 'stations.pickle')
-        savetriplets = os.path.join(self.savedir, 'triplets.pickle')
+        save_stations = os.path.join(self.savedir, 'stations.pickle')
+        save_events = os.path.join(self.savedir, 'events.pickle')
+        save_triplets = os.path.join(self.savedir, 'triplets.pickle')
         if recompute:
-            remove_file(savecoords)
-            remove_file(savetriplets)
-        if not os.path.exists(savecoords) or not os.path.exists(savetriplets):
-            coords, triplets = self.get_coords_and_triplets(
+            remove_file(save_stations)
+            remove_file(save_events)
+            remove_file(save_triplets)
+        exist_stations = os.path.exists(save_stations)
+        exist_events = os.path.exists(save_events)
+        exist_triplets = os.path.exists(save_triplets)
+        if not exist_stations or not exist_events or not exist_triplets:
+            stations, events_info, triplets = self.get_coords_and_triplets(
                     self.events, 
                     azimuth_tolerance=azimuth_tolerance,
                     distmin=distmin,
                     distmax=distmax
                     )
             triplets = {k: v for k, v in triplets.items() if len(v)>=min_no_events}
-            save_pickle(savecoords, coords)
-            save_pickle(savetriplets, triplets)
+            save_pickle(save_stations, stations)
+            save_pickle(save_events, events_info)
+            save_pickle(save_triplets, triplets)
         else:
-            coords = load_pickle(savecoords)
-            triplets = load_pickle(savetriplets)
+            stations = load_pickle(save_stations)
+            events_info = load_pickle(save_events)
+            triplets = load_pickle(save_triplets)
         self.min_no_events = min_no_events
-        self.stations = coords
+        self.stations = stations
+        self.events_info = events_info
         self.triplets = triplets
             
     
+    def get_events_used(self):
+        """ 
+        Retrieves the events id for which triplets of epicenter-receivers are
+        available to extract dispersion measurements
+        
+        Returns
+        -------
+        events_used : dict
+            Dictionary object where each key corresponds to an event (origin
+            time in obspy.UTCDateTime.timestamp format, i.e., the name of the
+            respective directory in self.src), and the associated values 
+            include all the station codes that exploit that event to extract
+            a dispersion measurement
+        
+        """
+        events_used = defaultdict(set)
+        for stapair, evlist in self.triplets.items():
+            for event in evlist:
+                events_used[event] = events_used[event].union(stapair)
+        return events_used
+    
+    
     def delete_unused_files(self):
         """ 
         Deletes every file in the data directory which is not useful for 
@@ -585,10 +655,7 @@ def remove_files(folder, files):
         
         deleted_files = 0
         deleted_folders = 0
-        events_used = defaultdict(set)
-        for stapair, evlist in self.triplets.items():
-            for event in evlist:
-                events_used[event] = events_used[event].union(stapair)
+        events_used = self.get_events_used()
         for event in sorted(events_used):
             folder = os.path.join(self.src, event)
             sacfiles = set([i for i in os.listdir(folder) if i.endswith('.sac')])
@@ -862,10 +929,9 @@ def display_progress(no_files, done, verbose=False):
     
     def plot_stations(self, ax=None, show=True, oceans_color='water', 
                       lands_color='land', edgecolor='k', projection='Mercator',
-                      color_by_network=True, **kwargs):
-        """ 
-        Maps the seismic receivers for which data triplets of epicenter-receivers
-        are available.
+                      resolution='110m', color_by_network=True, legend_dict={}, 
+                      **kwargs):
+        """ Maps the seismic receivers for which data are available
         
         Parameters
         ----------
@@ -893,61 +959,119 @@ def plot_stations(self, ax=None, show=True, oceans_color='water',
             (Visit the cartopy website for a list of valid projection names.)
             If ax is not None, `projection` is ignored. Default is 'Mercator'
             
+        resolution : str
+            Resolution of the Earth features displayed in the figure. Passed to
+            cartopy.feature.NaturalEarthFeature. Valid arguments are '110m',
+            '50m', '10m'. Default is '110m'
+            
         color_by_network : bool
             If True, each seismic network will have a different color in the
             resulting map, and a legend will be displayed. Otherwise, all
             stations will have the same color. Default is True
-            
-        kwargs : dict
-            Dictionary of additional key words that will be passed to the
-            matplotlib.pyplot.plot method.        
+        
+        legend_dict : dict, optional
+            Dictionary of keyword arguments passed to matplotlib.pyplot.legend
+        
+        kwargs : 
+            Additional keyword arguments passed to matplotlib.pyplot.plot 
             
             
         Returns
         -------
         If `show` is True, None, else `ax`, i.e. the GeoAxesSubplot
         """
-
-        def get_coords(stations, color_by_network):
-            codes, coords = zip(*[(k, v) for k, v in stations.items() \
-                                  if '_' not in k])
-            if not color_by_network:
-                yield np.array(coords), None
-            else:
-                networks = sorted(set(code.split('.')[0] for code in codes))
-                for net in networks:
-                    idx = [i for i, code in enumerate(codes) if code.startswith(net)]
-                    yield np.array(coords)[idx], net
-            
-        if ax is None:
-            projection = eval('ccrs.%s()'%projection)
-            fig = plt.figure()
-            ax = fig.add_subplot(1, 1, 1, projection=projection)
-            add_earth_features(ax, 
-                               oceans_color=oceans_color, 
-                               edgecolor=edgecolor,
-                               lands_color=lands_color)
-        transform = ccrs.PlateCarree()
-        marker = '^' if kwargs.get('marker', None) is None else kwargs['marker']
-        coords = np.array([i for i in self.stations.values()])
-        latmin, latmax = np.min(coords[:,0]), np.max(coords[:,0])
-        lonmin, lonmax = np.min(coords[:,1]), np.max(coords[:,1])
-        dlat = (latmax - latmin) * 0.1
-        dlon = (lonmax - lonmin) * 0.1
-        for coords, net in get_coords(self.stations, color_by_network):
-            ax.plot(*coords.T[::-1], marker=marker, lw=0, transform=transform,
-                    label=net, **kwargs) 
-        ax.set_extent((lonmin-dlon, lonmax+dlon, latmin-dlat, latmax+dlat), 
-                      transform)
-        if color_by_network:
-            ax.legend()
-        if show:
-            plt.show()
-        else:
-            return ax
+        
+        return plot_stations(stations=self.stations,
+                             ax=ax, 
+                             show=show, 
+                             oceans_color=oceans_color, 
+                             lands_color=lands_color, 
+                             edgecolor=edgecolor, 
+                             projection=projection,
+                             resolution=resolution,
+                             color_by_network=color_by_network, 
+                             legend_dict=legend_dict,
+                             **kwargs)      
                 
         
+    def plot_events(self, ax=None, show=True, oceans_color='water', lands_color='land', 
+                    edgecolor='k', projection='Mercator', resolution='110m', 
+                    min_size=5, max_size=150, legend_dict={}, **kwargs):
+        """ Creates a map of epicenters
+        
+        Parameters
+        ----------
+        lat, lon : ndarray of shape (n,)
+            Latitude and longitude of the epicenters
+            
+        mag : ndarray of shape(n,), optional
+            If not given, the size of each on the map will be constant
+        
+        ax : cartopy.mpl.geoaxes.GeoAxesSubplot
+            If not None, the receivers are plotted on the GeoAxesSubplot instance. 
+            Otherwise, a new figure and GeoAxesSubplot instance is created
+            
+        show : bool
+            If True, the plot is shown. Otherwise, a GeoAxesSubplot instance is
+            returned. Default is True
+            
+        oceans_color, lands_color : str
+            Color of oceans and lands. The arguments are ignored if ax is not
+            None. Otherwise, they are passed to cartopy.feature.GSHHSFeature 
+            (to the argument 'facecolor'). Defaults are 'water' and 'land'
+            
+        edgecolor : str
+            Color of the boundaries between, e.g., lakes and land. The argument 
+            is ignored if ax is not None. Otherwise, it is passed to 
+            cartopy.feature.GSHHSFeature (to the argument 'edgecolor'). Default
+            is 'k' (black)
+            
+        projection : str
+            Name of the geographic projection used to create the GeoAxesSubplot.
+            (Visit the cartopy website for a list of valid projection names.)
+            If ax is not None, `projection` is ignored. Default is 'Mercator'
+        
+        resolution : str
+            Resolution of the Earth features displayed in the figure. Passed to
+            cartopy.feature.NaturalEarthFeature. Valid arguments are '110m',
+            '50m', '10m'. Default is '110m'
+        
+        min_size, max_size : int, float
+            Minimum and maximum size of the epicenters on the map. These are used
+            to interpolate all magnitudes associated with each event, so as to
+            scale them appropriately on the map. (The final "sizes" are passed to
+            the argument `s` of matplotlib.pyplot.scatter)
+                
+        legend_dict : dict
+            Keyword arguments passed to matplotlib.pyplot.legend
+            
+        kwargs : 
+            Additional keyword arguments passed to matplotlib.pyplot.scatter
+            
+            
+        Returns
+        -------
+        If `show` is True, None, else `ax`, i.e. the GeoAxesSubplot
+        """
         
+        events_used = self.get_events_used().keys()
+        events_info = np.array([j for i, j in self.events_info.items() \
+                                if i in events_used])
+        lat, lon, mag = events_info.T
+        return plot_events(lat=lat,
+                           lon=lon,
+                           mag=mag,
+                           ax=ax, 
+                           show=show, 
+                           oceans_color=oceans_color, 
+                           lands_color=lands_color, 
+                           edgecolor=edgecolor, 
+                           projection=projection,
+                           resolution=resolution, 
+                           min_size=min_size,
+                           max_size=max_size,
+                           legend_dict=legend_dict, 
+                           **kwargs)
         
 ###############################################################################
         
diff --git a/seislib/plotting/plotting.py b/seislib/plotting/plotting.py
index 4256464..576ffc1 100644
--- a/seislib/plotting/plotting.py
+++ b/seislib/plotting/plotting.py
@@ -6,6 +6,7 @@
 @email2: fabrizio.magrini90@gmail.com
 """
 
+from collections.abc import Iterable
 import numpy as np
 from scipy.ndimage.filters import gaussian_filter
 import matplotlib.pyplot as plt
@@ -307,8 +308,9 @@ def contour(mesh, c, ax, smoothing=None, **kwargs):
 
 def plot_stations(stations, ax=None, show=True, oceans_color='water', 
                   lands_color='land', edgecolor='k', projection='Mercator',
-                  color_by_network=True, legend_dict={}, **kwargs):
-    """ Maps the seismic receivers for which data are available
+                  resolution='110m', color_by_network=True, legend_dict={}, 
+                  **kwargs):
+    """ Creates a maps of seismic receivers
     
     Parameters
     ----------
@@ -348,7 +350,12 @@ def plot_stations(stations, ax=None, show=True, oceans_color='water',
         Name of the geographic projection used to create the GeoAxesSubplot.
         (Visit the cartopy website for a list of valid projection names.)
         If ax is not None, `projection` is ignored. Default is 'Mercator'
-        
+    
+    resolution : str
+        Resolution of the Earth features displayed in the figure. Passed to
+        cartopy.feature.NaturalEarthFeature. Valid arguments are '110m',
+        '50m', '10m'. Default is '110m'
+    
     color_by_network : bool
         If True, each seismic network will have a different color in the
         resulting map, and a legend will be displayed. Otherwise, all
@@ -382,8 +389,8 @@ def get_map_boundaries(stations):
         coords = np.array([i for i in stations.values()])
         latmin, latmax = np.min(coords[:,0]), np.max(coords[:,0])
         lonmin, lonmax = np.min(coords[:,1]), np.max(coords[:,1])
-        dlat = (latmax - latmin) * 0.02
-        dlon = (lonmax - lonmin) * 0.02
+        dlat = (latmax - latmin) * 0.03
+        dlon = (lonmax - lonmin) * 0.03
         lonmin = lonmin-dlon if lonmin-dlon > -180 else lonmin
         lonmax = lonmax+dlon if lonmax+dlon < 180 else lonmax
         latmin = latmin-dlat if latmin-dlat > -90 else latmin
@@ -395,13 +402,15 @@ def get_map_boundaries(stations):
         projection = eval('ccrs.%s()'%projection)
         fig = plt.figure()
         ax = fig.add_subplot(1, 1, 1, projection=projection)
+        transform = ccrs.PlateCarree()
+        map_boundaries = get_map_boundaries(stations)
         add_earth_features(ax, 
+                           scale=resolution,
                            oceans_color=oceans_color, 
                            edgecolor=edgecolor,
                            lands_color=lands_color)
-    transform = ccrs.PlateCarree()
-    marker = '^' if kwargs.get('marker', None) is None else kwargs['marker']
-    map_boundaries = get_map_boundaries(stations)
+    
+    marker = kwargs.pop('marker', '^')
     for coords, net in get_coords(stations, color_by_network):
         ax.plot(*coords.T[::-1], marker=marker, lw=0, transform=transform,
                 label=net, **kwargs) 
@@ -414,11 +423,130 @@ def get_map_boundaries(stations):
         return ax
 
 
+def plot_events(lat, lon, mag=None, ax=None, show=True, oceans_color='water', 
+                lands_color='land', edgecolor='k', projection='Mercator',
+                resolution='110m', min_size=5, max_size=200, legend_dict={}, 
+                **kwargs):
+    """ Creates a map of epicenters
+    
+    Parameters
+    ----------
+    lat, lon : ndarray of shape (n,)
+        Latitude and longitude of the epicenters
+        
+    mag : ndarray of shape(n,), optional
+        If not given, the size of each on the map will be constant
+    
+    ax : cartopy.mpl.geoaxes.GeoAxesSubplot
+        If not None, the receivers are plotted on the GeoAxesSubplot instance. 
+        Otherwise, a new figure and GeoAxesSubplot instance is created
+        
+    show : bool
+        If True, the plot is shown. Otherwise, a GeoAxesSubplot instance is
+        returned. Default is True
+        
+    oceans_color, lands_color : str
+        Color of oceans and lands. The arguments are ignored if ax is not
+        None. Otherwise, they are passed to cartopy.feature.GSHHSFeature 
+        (to the argument 'facecolor'). Defaults are 'water' and 'land'
+        
+    edgecolor : str
+        Color of the boundaries between, e.g., lakes and land. The argument 
+        is ignored if ax is not None. Otherwise, it is passed to 
+        cartopy.feature.GSHHSFeature (to the argument 'edgecolor'). Default
+        is 'k' (black)
+        
+    projection : str
+        Name of the geographic projection used to create the GeoAxesSubplot.
+        (Visit the cartopy website for a list of valid projection names.)
+        If ax is not None, `projection` is ignored. Default is 'Mercator'
+    
+    resolution : str
+        Resolution of the Earth features displayed in the figure. Passed to
+        cartopy.feature.NaturalEarthFeature. Valid arguments are '110m',
+        '50m', '10m'. Default is '110m'
+    
+    min_size, max_size : int, float
+        Minimum and maximum size of the epicenters on the map. These are used
+        to interpolate all magnitudes associated with each event, so as to
+        scale them appropriately on the map. (The final "sizes" are passed to
+        the argument `s` of matplotlib.pyplot.scatter)
+            
+    legend_dict : dict
+        Keyword arguments passed to matplotlib.pyplot.legend
+        
+    kwargs : 
+        Additional keyword arguments passed to matplotlib.pyplot.scatter
+        
+        
+    Returns
+    -------
+    If `show` is True, None, else `ax`, i.e. the GeoAxesSubplot
+    """
+    
+    def get_map_boundaries(lat, lon):
+        latmin, latmax = np.min(lat), np.max(lat)
+        lonmin, lonmax = np.min(lon), np.max(lon)
+        dlat = (latmax - latmin) * 0.03
+        dlon = (lonmax - lonmin) * 0.03
+        lonmin = lonmin-dlon if lonmin-dlon > -180 else lonmin
+        lonmax = lonmax+dlon if lonmax+dlon < 180 else lonmax
+        latmin = latmin-dlat if latmin-dlat > -90 else latmin
+        latmax = latmax+dlat if latmax+dlat < 90 else latmax
+        return (lonmin, lonmax, latmin, latmax)
+    
+    def get_markers_size(mag, kwargs):
+        if mag is None:
+            return kwargs.pop('s', None)
+        x = np.linspace(min(mag), max(mag))
+        y = np.geomspace(min_size, max_size)
+        return np.interp(mag, x, y)
+    
+    def get_integer_magnitudes(mag):
+        magmin, magmax = min(mag), max(mag)
+        return np.round(np.geomspace(magmin, magmax, 4), 1)
+            
+    if ax is None:
+        projection = eval('ccrs.%s()'%projection)
+        fig = plt.figure()
+        ax = fig.add_subplot(1, 1, 1, projection=projection)
+        transform = ccrs.PlateCarree()
+        map_boundaries = get_map_boundaries(lat, lon)
+        ax.set_extent(map_boundaries, transform)
+        add_earth_features(ax, 
+                           scale=resolution,
+                           oceans_color=oceans_color, 
+                           edgecolor=edgecolor,
+                           lands_color=lands_color)
+        
+    marker = kwargs.pop('marker', '*')
+    s = kwargs.pop('s', None)
+    if s is None:
+        s = get_markers_size(mag, kwargs)
+    if isinstance(s, Iterable):
+        mags_legend = get_integer_magnitudes(mag)
+        idx = [np.argmin(np.abs(mag - mag_i)) for mag_i in mags_legend]
+        idx_all = np.setdiff1d(range(len(mag)), idx)
+        ax.scatter(lon[idx_all], lat[idx_all], marker=marker, transform=transform, 
+                   s=s[idx_all], zorder=100, **kwargs) 
+        for i, mag_legend in zip(idx, mags_legend):
+            ax.scatter(lon[i], lat[i], marker=marker, transform=transform, 
+                       s=s[i], label=mag_legend, zorder=100, **kwargs) 
+        ax.legend(**legend_dict)
+    else:
+        ax.scatter(lon, lat, marker=marker, transform=transform, s=s, zorder=100,
+                   **kwargs) 
+    if show:
+        plt.show()
+    else:
+        return ax
+    
+
 def plot_rays(data_coords, ax=None, show=True, stations_color='r', 
               paths_color='k', oceans_color='water', lands_color='land', 
               edgecolor='k', stations_alpha=None, paths_alpha=0.3, 
-              projection='Mercator', map_boundaries=None, bound_map=True, 
-              paths_width=0.2, **kwargs):
+              projection='Mercator', resolution='110m', map_boundaries=None, 
+              bound_map=True, paths_width=0.2, **kwargs):
     """ 
     Utility function to display the great-circle paths associated with pairs
     of data coordinates
@@ -463,6 +591,11 @@ def plot_rays(data_coords, ax=None, show=True, stations_color='r',
         Name of the geographic projection used to create the GeoAxesSubplot.
         (Visit the cartopy website for a list of valid projection names.)
         If ax is not None, `projection` is ignored. Default is 'Mercator'
+
+    resolution : str
+        Resolution of the Earth features displayed in the figure. Passed to
+        cartopy.feature.NaturalEarthFeature. Valid arguments are '110m',
+        '50m', '10m'. Default is '110m'
         
     map_boundaries : list or tuple of floats, shape (4,), optional
         Lonmin, lonmax, latmin, latmax (in degrees) defining the extent of
@@ -484,8 +617,8 @@ def get_map_boundaries(data):
         lons = np.concatenate((data[:,1], data[:, 3]))
         latmin, latmax = np.nanmin(lats), np.nanmax(lats)
         lonmin, lonmax = np.nanmin(lons), np.nanmax(lons)
-        dlon = (lonmax - lonmin) * 0.02
-        dlat = (latmax - latmin) * 0.02
+        dlon = (lonmax - lonmin) * 0.03
+        dlat = (latmax - latmin) * 0.03
         lonmin = lonmin-dlon if lonmin-dlon > -180 else lonmin
         lonmax = lonmax+dlon if lonmax+dlon < 180 else lonmax
         latmin = latmin-dlat if latmin-dlat > -90 else latmin
@@ -497,6 +630,7 @@ def get_map_boundaries(data):
         fig = plt.figure()
         ax = fig.add_subplot(1, 1, 1, projection=projection)
         add_earth_features(ax, 
+                           scale=resolution,
                            oceans_color=oceans_color, 
                            edgecolor=edgecolor,
                            lands_color=lands_color)
@@ -544,7 +678,7 @@ def get_map_boundaries(data):
 def plot_colored_rays(data_coords, c, ax=None, show=True, cmap=scm.roma,
                       vmin=None, vmax=None, stations_color='k', 
                       oceans_color='lightgrey', lands_color='w', edgecolor='k', 
-                      stations_alpha=None, paths_alpha=None, 
+                      stations_alpha=None, paths_alpha=None, resolution='110m',
                       projection='Mercator', map_boundaries=None, bound_map=True, 
                       paths_width=1, colorbar=True, cbar_dict={}, **kwargs):
     """ 
@@ -605,6 +739,11 @@ def plot_colored_rays(data_coords, c, ax=None, show=True, cmap=scm.roma,
         (Visit the cartopy website for a list of valid projection names.)
         If ax is not None, `projection` is ignored. Default is 'Mercator'
         
+    resolution : str
+        Resolution of the Earth features displayed in the figure. Passed to
+        cartopy.feature.NaturalEarthFeature. Valid arguments are '110m',
+        '50m', '10m'. Default is '110m'
+        
     map_boundaries : list or tuple of floats, shape (4,), optional
         Lonmin, lonmax, latmin, latmax (in degrees) defining the extent of
         the map
@@ -635,8 +774,8 @@ def get_map_boundaries(data):
         lons = np.concatenate((data[:,1], data[:, 3]))
         latmin, latmax = np.nanmin(lats), np.nanmax(lats)
         lonmin, lonmax = np.nanmin(lons), np.nanmax(lons)
-        dlon = (lonmax - lonmin) * 0.02
-        dlat = (latmax - latmin) * 0.02
+        dlon = (lonmax - lonmin) * 0.03
+        dlat = (latmax - latmin) * 0.03
         lonmin = lonmin-dlon if lonmin-dlon > -180 else lonmin
         lonmax = lonmax+dlon if lonmax+dlon < 180 else lonmax
         latmin = latmin-dlat if latmin-dlat > -90 else latmin
@@ -665,6 +804,7 @@ def create_colourbar(ax, cmap, norm, **kwargs):
         fig = plt.figure()
         ax = fig.add_subplot(1, 1, 1, projection=projection)
         add_earth_features(ax, 
+                           scale=resolution,
                            oceans_color=oceans_color, 
                            edgecolor=edgecolor,
                            lands_color=lands_color)
@@ -716,7 +856,7 @@ def create_colourbar(ax, cmap, norm, **kwargs):
 
 def plot_map(mesh, c, ax=None, projection='Mercator', map_boundaries=None, 
              bound_map=True, colorbar=True, show=True, style='colormesh', 
-             add_features=False, cbar_dict={}, **kwargs):
+             add_features=False, resolution='110m', cbar_dict={}, **kwargs):
     """ 
     Utility function to display the lateral variations of some quantity on 
     a equal-area grid
@@ -763,6 +903,11 @@ def plot_map(mesh, c, ax=None, projection='Mercator', map_boundaries=None,
         If True, natural Earth features will be added to the GeoAxesSubplot.
         Default is False. If `ax` is None, it is automatically set to True
     
+    resolution : str
+        Resolution of the Earth features displayed in the figure. Passed to
+        cartopy.feature.NaturalEarthFeature. Valid arguments are '110m',
+        '50m', '10m'. Default is '110m'
+    
     cbar_dict : dict
         Keyword arguments passed to matplotlib.colorbar.ColorbarBase
     
@@ -781,8 +926,8 @@ def plot_map(mesh, c, ax=None, projection='Mercator', map_boundaries=None,
     def get_map_boundaries(mesh):
         latmin, lonmin = np.min(mesh, axis=0)[::2]
         latmax, lonmax = np.max(mesh, axis=0)[1::2]
-        dlon = (lonmax - lonmin) * 0.02
-        dlat = (latmax - latmin) * 0.02
+        dlon = (lonmax - lonmin) * 0.03
+        dlat = (latmax - latmin) * 0.03
         
         lonmin = lonmin-dlon if lonmin-dlon > -180 else lonmin
         lonmax = lonmax+dlon if lonmax+dlon < 180 else lonmax
@@ -829,6 +974,7 @@ def set_colorbar_aspect(cb, kwargs):
         ax = fig.add_subplot(1, 1, 1, projection=projection)        
     if add_features:
         add_earth_features(ax, 
+                           scale=resolution,
                            oceans_color='none', 
                            edgecolor='k',
                            lands_color='none')
diff --git a/seislib/tomography/tomography.py b/seislib/tomography/tomography.py
index d552a8b..bb35dcf 100644
--- a/seislib/tomography/tomography.py
+++ b/seislib/tomography/tomography.py
@@ -107,14 +107,14 @@ class SeismicTomography:
     plot_rays(ax=None, show=True, stations_color='r', gc_paths_color='k',
               oceans_color='water', lands_color='land', edgecolor='k',
               alpha_stations=None, alpha_paths=0.3, projection='Mercator', 
-              map_boundaries=None, bound_map=True, **kwargs)
+              resolution='110m', map_boundaries=None, bound_map=True, **kwargs)
         Utility function to display the great circle paths associated with pairs
         of data coordinates
         
     plot_colored_rays(ax=None, show=True, cmap=scm.roma, vmin=None, vmax=None, 
                       stations_color='k', oceans_color='lightgrey', 
                       lands_color='w', edgecolor='k', stations_alpha=None, 
-                      paths_alpha=None, projection='Mercator', 
+                      paths_alpha=None, projection='Mercator', resolution='110m',
                       map_boundaries=None, bound_map=True, paths_width=1, 
                       colorbar=True, cbar_dict={}, **kwargs):
         Utility function to display the great-circle paths associated with pairs
@@ -143,7 +143,7 @@ class SeismicTomography:
 
     plot_map(mesh, c, ax=None, projection='Mercator', map_boundaries=None, 
              bound_map=True, colorbar=True, show=True, style='colormesh', 
-             add_features=False, cbar_dict={}, **kwargs)
+             add_features=False, resolution='110m', cbar_dict={}, **kwargs)
         Utility function to display the lateral variations of some quantity on 
         a equal-area grid
         
@@ -1341,7 +1341,7 @@ def spike_test(self, x0, y0, sigma_x, sigma_y, latmin=None, latmax=None,
     @classmethod
     def plot_map(cls, mesh, c, ax=None, projection='Mercator', map_boundaries=None, 
                  bound_map=True, colorbar=True, show=True, style='colormesh', 
-                 add_features=False, cbar_dict={}, **kwargs):
+                 add_features=False, resolution='110m', cbar_dict={}, **kwargs):
         """ 
         Utility function to display the lateral variations of some quantity on 
         a equal-area grid
@@ -1388,6 +1388,11 @@ def plot_map(cls, mesh, c, ax=None, projection='Mercator', map_boundaries=None,
             If True, natural Earth features will be added to the GeoAxesSubplot.
             Default is False. If `ax` is None, it is automatically set to True
             
+        resolution : str
+            Resolution of the Earth features displayed in the figure. Passed to
+            cartopy.feature.NaturalEarthFeature. Valid arguments are '110m',
+            '50m', '10m'. Default is '110m'
+            
         cbar_dict : dict
             Keyword arguments passed to matplotlib.colorbar.ColorbarBase
             
@@ -1413,6 +1418,7 @@ def plot_map(cls, mesh, c, ax=None, projection='Mercator', map_boundaries=None,
                                  show=show, 
                                  style=style, 
                                  add_features=add_features, 
+                                 resolution=resolution,
                                  cbar_dict=cbar_dict,
                                  **kwargs)
 
@@ -1517,8 +1523,8 @@ def contour(cls, mesh, c, ax, smoothing=None, **kwargs):
     def plot_rays(self, ax=None, show=True, stations_color='r', 
                   paths_color='k', oceans_color='water', lands_color='land', 
                   edgecolor='k', stations_alpha=None, paths_alpha=0.3, 
-                  projection='Mercator', map_boundaries=None, bound_map=True, 
-                  paths_width=0.2, **kwargs):
+                  projection='Mercator', resolution='110m', map_boundaries=None, 
+                  bound_map=True, paths_width=0.2, **kwargs):
         """ 
         Utility function to display the great-circle paths associated with pairs
         of data coordinates
@@ -1560,7 +1566,12 @@ def plot_rays(self, ax=None, show=True, stations_color='r',
             Name of the geographic projection used to create the GeoAxesSubplot.
             (Visit the cartopy website for a list of valid projection names.)
             If ax is not None, `projection` is ignored. Default is 'Mercator'
-            
+     
+        resolution : str
+            Resolution of the Earth features displayed in the figure. Passed to
+            cartopy.feature.NaturalEarthFeature. Valid arguments are '110m',
+            '50m', '10m'. Default is '110m'
+        
         map_boundaries : list or tuple of floats, shape (4,), optional
             Lonmin, lonmax, latmin, latmax (in degrees) defining the extent of
             the map
@@ -1586,6 +1597,7 @@ def plot_rays(self, ax=None, show=True, stations_color='r',
                                   stations_alpha=stations_alpha, 
                                   paths_alpha=paths_alpha, 
                                   projection=projection, 
+                                  resolution=resolution,
                                   map_boundaries=map_boundaries, 
                                   bound_map=bound_map, 
                                   paths_width=paths_width, 
@@ -1595,7 +1607,7 @@ def plot_rays(self, ax=None, show=True, stations_color='r',
     def plot_colored_rays(self, ax=None, show=True, cmap=scm.roma, vmin=None, 
                           vmax=None, stations_color='k', oceans_color='lightgrey', 
                           lands_color='w', edgecolor='k', stations_alpha=None, 
-                          paths_alpha=None, projection='Mercator', 
+                          paths_alpha=None, projection='Mercator', resolution='110m',
                           map_boundaries=None, bound_map=True, paths_width=1, 
                           colorbar=True, cbar_dict={}, **kwargs):
         """ 
@@ -1656,6 +1668,11 @@ def plot_colored_rays(self, ax=None, show=True, cmap=scm.roma, vmin=None,
             (Visit the cartopy website for a list of valid projection names.)
             If ax is not None, `projection` is ignored. Default is 'Mercator'
             
+        resolution : str
+            Resolution of the Earth features displayed in the figure. Passed to
+            cartopy.feature.NaturalEarthFeature. Valid arguments are '110m',
+            '50m', '10m'. Default is '110m'
+            
         map_boundaries : list or tuple of floats, shape (4,), optional
             Lonmin, lonmax, latmin, latmax (in degrees) defining the extent of
             the map
@@ -1695,6 +1712,7 @@ def plot_colored_rays(self, ax=None, show=True, cmap=scm.roma, vmin=None,
                                           stations_alpha=stations_alpha, 
                                           paths_alpha=paths_alpha, 
                                           projection=projection, 
+                                          resolution=resolution,
                                           map_boundaries=map_boundaries, 
                                           bound_map=bound_map, 
                                           paths_width=paths_width, 
diff --git a/setup.cfg b/setup.cfg
new file mode 100644
index 0000000..b88034e
--- /dev/null
+++ b/setup.cfg
@@ -0,0 +1,2 @@
+[metadata]
+description-file = README.md
diff --git a/setup.py b/setup.py
index d252709..26c2c76 100644
--- a/setup.py
+++ b/setup.py
@@ -53,6 +53,7 @@ def generate_cython():
         author_email="fabrizio.magrini90@gmail.com",
         license="MIT",
         packages=find_packages(),
+        python_requires=">=3.6",
         keywords="Seismic Imaging, Surface Waves, Seismic Ambient Noise, Earthquake Seismology, Tomographic Inversion",
         configuration=configuration,
         install_requires=['obspy>=1.1.0',
@@ -131,4 +132,4 @@ def parse_setuppy_commands():
     # Don't import numpy in other cases - non-build actions are required to succeed without NumPy,
     # for example when pip is used to install this when NumPy is not yet present in the system.
 
-    setup(**pkg_metadata)
\ No newline at end of file
+    setup(**pkg_metadata)