Update

.gitignore

@@ -3,6 +3,8 @@ __pycache__/
 *.pyc
 *.pyd
 *.so
+*.cpp
+*.c
 
 # Setuptools/distutils distribution folder.
 dist/

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
 
 

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,16 +96,16 @@ 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        
         
     
     Class Methods
     -------------
     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)

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)      

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):