New plot type: ccf_freq

msnoise/api.py

@@ -193,7 +193,7 @@ def update_config(session, name, value):
     """
 
     config = session.query(Config).filter(Config.name == name).first()
-    if "NULL" in value: 
+    if "NULL" in value:
         config.value = None
     else:
         config.value = value
@@ -996,7 +996,7 @@ def export_mseed(db, filename, pair, components, filterid, corr, ncorr=0,
     except:
         pass
     filename += ".MSEED"
-    
+
     if maxlag is None:
         maxlag = float(get_config(db, "maxlag"))
     if cc_sampling_rate is None:
@@ -1053,7 +1053,15 @@ def get_results(session, station1, station2, filterid, components, dates,
     i = 0
     base = os.path.join("STACKS", "%02i" % filterid,
                         "%03i_DAYS" % mov_stack, components,
-                        "%s_%s" % (station1, station2), "%s")
+                        "%s_%s" % (station1, station2))
+
+    if not os.path.exists(base):
+        raise NotADirectoryError('Path ' + base + ' does not exist. Check if '
+                                 'provided station names are following '
+                                 'NETW.STATION pattern.')
+
+    base = os.path.join(base, "%s")
+
     if export_format == "BOTH":
         base += ".MSEED"
         export_format = "MSEED"
@@ -1063,7 +1071,7 @@ def get_results(session, station1, station2, filterid, components, dates,
         base += ".MSEED"
     logging.debug("Reading files...")
     for j, date in enumerate(dates):
-        daystack = base % str(date)
+        daystack = base  % str(date)
         try:
             stack_data[j, :] = read(daystack, format=export_format)[0].data[:]
             i += 1
@@ -1187,7 +1195,7 @@ def build_ref_datelist(session):
     return start, end, datelist.tolist()
 
 
-def build_movstack_datelist(session):
+def build_movstack_datelist(session, startdate=None, enddate=None):
     """
     Creates a date array for the analyse period.
     The returned tuple contains a start and an end date, and a list of
@@ -1196,12 +1204,19 @@ def build_movstack_datelist(session):
     :type session: :class:`sqlalchemy.orm.session.Session`
     :param session: A :class:`~sqlalchemy.orm.session.Session` object, as
         obtained by :func:`connect`
-
-    :rtype: tuple
-    :returns: (start, end, datelist)
-    """
-    begin = get_config(session, "startdate")
-    end = get_config(session, "enddate")
+    :type startdate: str
+    :param startdate: a startdate for which method should create a date array.
+        Defaults to None, in this case method gets the value from database.
+    :type enddate: str
+    :param enddate: a enddate for which method should create a date array.
+        Defaults to None, in this case method gets the value from database.
+
+    :rtype: (datetime.date, datetime.date, list of datetime.date)
+    :returns: Returns startdate, enddate and a list of all days between those.
+    """
+    begin = startdate if startdate is not None else get_config(session,
+                                                               "startdate")
+    end = enddate if enddate is not None else get_config(session, "enddate")
     if begin[0] == '-':
         start = datetime.date.today() + datetime.timedelta(days=int(begin))
         end = datetime.date.today() + datetime.timedelta(days=int(end))
@@ -1433,7 +1448,7 @@ def make_same_length(st):
     # apply the mask to all traces
     for tr in st:
         tr.data.mask = mask
-    
+
     # remove the masks from the stream 
     st = st.split()
     return st

msnoise/msnoise_admin.py

@@ -182,23 +182,23 @@ def mwcs_step(form, field):
         if field.data > form.data['mwcs_wlen']:
             raise ValidationError("'mwcs_step' should be smaller or equal to"
                                   " 'mwcs_wlen'")
-    
+
     form_args = dict(
         mwcs_low=dict(validators=[mwcs_low]),
         mwcs_high=dict(validators=[mwcs_high]),
         high=dict(validators=[high]),
         mwcs_step=dict(validators=[mwcs_step]),
     )
-    
+
     column_list = ('ref', 'low', 'mwcs_low', 'mwcs_high', 'high',
                    'rms_threshold', 'mwcs_wlen', 'mwcs_step', 'used')
     form_columns = ('low', 'mwcs_low', 'mwcs_high', 'high',
                     'rms_threshold', 'mwcs_wlen', 'mwcs_step', 'used')
-    
+
     def __init__(self, session, **kwargs):
         # You can pass name and other parameters if you want to
         super(FilterView, self).__init__(Filter, session, **kwargs)
-        
+
     @action('used',
             lazy_gettext('Toggle Used'),
             lazy_gettext('Are you sure you want to update selected models?'))
@@ -211,7 +211,7 @@ def used(self, ids):
             else:
                 s.used = True
         self.session.commit()
-        return  
+        return
 
 
 MY_DEFAULT_FORMATTERS = dict(typefmt.BASE_FORMATTERS)
@@ -224,10 +224,10 @@ class StationView(ModelView):
     view_title = "Station Configuration"
     column_filters = ('net', 'used')
     column_type_formatters = MY_DEFAULT_FORMATTERS
-    
+
     def __init__(self, session, **kwargs):
         super(StationView, self).__init__(Station, session, **kwargs)
-    
+
     @action('used',
             lazy_gettext('Toggle Used'),
             lazy_gettext('Are you sure you want to update selected models?'))
@@ -240,7 +240,7 @@ def used(self, ids):
             else:
                 s.used = True
         self.session.commit()
-        return  
+        return
 
 
 class DataAvailabilityView(ModelView):
@@ -255,7 +255,7 @@ class DataAvailabilityView(ModelView):
     def __init__(self, session, **kwargs):
         super(DataAvailabilityView, self).__init__(DataAvailability, session,
                                                    **kwargs)
-    
+
     @action('modified',
             lazy_gettext('Mark as (M)odified'),
             lazy_gettext('Are you sure you want to update selected models?'))
@@ -271,7 +271,7 @@ def modified(self, ids):
                        '%(count)s models were successfully flagged (M)odified.',
               count, count=count))
         return
-    
+
 
 class JobView(ModelView):
     view_title = "Jobs"
@@ -284,7 +284,7 @@ class JobView(ModelView):
 
     def __init__(self, session, **kwargs):
         super(JobView, self).__init__(Job, session, **kwargs)
-    
+
     @action('todo',
             lazy_gettext('Mark as (T)odo'),
             lazy_gettext('Are you sure you want to update selected models?'))
@@ -295,7 +295,7 @@ def todo(self, ids):
             s.flag = 'T'
         self.session.commit()
         return
-    
+
     @action('done',
             lazy_gettext('Mark as (D)one'),
             lazy_gettext('Are you sure you want to update selected models?'))
@@ -306,7 +306,7 @@ def done(self, ids):
             s.flag = 'D'
         self.session.commit()
         return
-    
+
     @action('deletetype',
             lazy_gettext('Delete all Jobs of the same "Type"'),
             lazy_gettext('Are you sure you want to delete all those models?'))
@@ -318,7 +318,7 @@ def deletetype(self, ids):
         self.get_query().filter(Job.jobtype == type_to_delete).delete()
         self.session.commit()
         return
-    
+
     @action('massTodo',
             lazy_gettext('Mark all Jobs of the same Type as (T)odo'),
             lazy_gettext('Are you sure you want to update all those models?'))
@@ -327,7 +327,7 @@ def massTodo(self, ids):
         query = self.get_query().filter(model_pk.in_(ids))
         for s in query.all():
             type_to_delete = s.jobtype
-        
+
         for s in self.get_query().filter(Job.jobtype == type_to_delete).all():
             s.flag = 'T'
         self.session.commit()
@@ -564,7 +564,7 @@ def allresults():
     start, end, dates = build_ref_datelist(db)
     i, result = get_results(db,station1, station2, filterid, components, dates,
                             format=format)
-    
+
     data = {}
     if format == 'stack':
         if i != 0:
@@ -603,7 +603,7 @@ def new_jobsTRIG():
     o['count'] = count
     o = json.dumps(o)
     return flask.Response(o, mimetype='application/json')
-    
+
 
 @app.route('/admin/jobs_list.json')
 def joblists():
@@ -666,7 +666,7 @@ def main(port=5000):
     db.close()
 
     admin = Admin(app, template_mode='bootstrap2')
-    
+
     if "msnoise_brand" in os.environ:
         tmp = eval(os.environ["msnoise_brand"])
         name, logo = tmp.split("|")

msnoise/plots/ccffreq.py

@@ -0,0 +1,133 @@
+"""
+This plot shows the cross-correlation functions (CCF) vs frequency. The
+parameters allow to plot the daily or the mov-stacked CCF. Filters and
+components are selectable too. The ``--ampli`` argument allows to increase the
+vertical scale of the CCFs. Passing ``--refilter`` allows to bandpass filter
+CCFs before plotting. Passing ``--startdate`` and ``--enddate`` parameters
+allows to specify which period of data should be plotted. By default the plot
+uses dates determined in database.
+
+.. include:: clickhelp/msnoise-plot-ccffreq.rst
+
+
+Example:
+
+``msnoise plot ccftime ID.KWUI ID.POSI`` will plot all defaults:
+
+.. image:: .static/ccffreq.png
+"""
+# plot ccffreq
+
+import datetime
+
+import matplotlib.dates as mdates
+import matplotlib.pyplot as plt
+import numpy as np
+from matplotlib.widgets import Cursor
+from obspy.signal.filter import bandpass
+
+from msnoise.api import build_movstack_datelist, connect, get_config, \
+    get_filters, get_results
+
+
+def main(sta1, sta2, filterid, components, mov_stack=1, ampli=5, show=False,
+         outfile=None, refilter=None, startdate=None, enddate=None):
+
+    db = connect()
+    cc_sampling_rate = float(get_config(db, 'cc_sampling_rate'))
+    start, end, datelist = build_movstack_datelist(db, startdate, enddate)
+    base = mdates.date2num(start)
+    sta1 = sta1.replace('.', '_')
+    sta2 = sta2.replace('.', '_')
+
+    # TODO: Height adjustment of the plot for large number of stacks.
+    # Preferably interactive
+    fig = plt.figure(figsize=(12, 9))
+
+    if refilter:
+        freqmin, freqmax = refilter.split(':')
+        freqmin = float(freqmin)
+        freqmax = float(freqmax)
+
+    if sta2 >= sta1:
+        pair = "%s:%s" % (sta1, sta2)
+
+        print("New Data for %s-%s-%i-%i" % (pair, components, filterid,
+                                            mov_stack))
+        nstack, stack_total = get_results(db, sta1, sta2, filterid, components,
+                                          datelist, mov_stack, format="matrix")
+        ax = fig.add_subplot(111)
+        for i, line in enumerate(stack_total):
+            if np.all(np.isnan(line)):
+                continue
+
+            freq, line = prepare_abs_postitive_fft(line, cc_sampling_rate)
+
+            if refilter:
+                line = bandpass(line, freqmin, freqmax, cc_sampling_rate,
+                                zerophase=True)
+
+            ax.plot(freq, line * ampli + i + base, c='k')
+
+        for filterdb in get_filters(db, all=True):
+            if filterid == filterdb.ref:
+                low = float(filterdb.low)
+                high = float(filterdb.high)
+                break
+
+        ax.set_ylim(start-datetime.timedelta(days=ampli),
+                    end+datetime.timedelta(days=ampli))
+        ax.yaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d'))
+
+        ax.set_xlabel("Frequency [Hz]")
+        ax.set_xscale('log')
+        ax.grid()
+
+        title = '%s : %s, %s, Filter %d (%.2f - %.2f Hz), Stack %d' %\
+                (sta1.replace('_', '.'), sta2.replace('_', '.'), components,
+                 filterid, low, high, mov_stack)
+        if refilter:
+            title += ", Re-filtered (%.2f - %.2f Hz)" % (freqmin, freqmax)
+        ax.set_title(title)
+
+        cursor = Cursor(ax, useblit=True, color='red', linewidth=1.2)
+
+        if outfile:
+            if outfile.startswith("?"):
+                pair = pair.replace(':', '-')
+                outfile = outfile.replace('?', '%s-%s-f%i-m%i' % (pair,
+                                                                  components,
+                                                                  filterid,
+                                                                  mov_stack))
+            outfile = "ccffreq_" + outfile
+            print("output to:", outfile)
+            fig.savefig(outfile)
+        if show:
+            fig.show()
+        else:
+            plt.close(fig)
+
+
+def prepare_abs_postitive_fft(line, sampling_rate):
+    """
+    Method that returns a positive part of FFT of provided signal along with
+    a corresponding frequency vector.
+
+    :type line: numpy.ndarray
+    :param line: Signal to calculate fft.
+    :type sampling_rate: float
+    :param sampling_rate: Sampling rate of provided signal
+
+    :rtype: tuple(numpy.ndarray, numpy.ndarray)
+    :return: Tuple of two arrays. One contains frequency vector for positive
+    part of FFT, second contains positive and absolute FFT of input array.
+    """
+    val = np.fft.fft(line)
+    val = np.abs(val)
+
+    freq = np.fft.fftfreq(len(line), (1/sampling_rate))
+    freq = [x for x in freq if x >= 0]
+
+    val = val[:len(freq)]
+
+    return freq, val

msnoise/plots/ccftime.py

@@ -27,12 +27,12 @@
 
 
 def main(sta1, sta2, filterid, components, mov_stack=1, ampli=5, seismic=False,
-         show=False, outfile=None, envelope=False, refilter=None):
+         show=False, outfile=None, envelope=False, refilter=None,  startdate=None, enddate=None):
     db = connect()
     maxlag = float(get_config(db, 'maxlag'))
     samples = get_maxlag_samples(db)
     cc_sampling_rate = float(get_config(db, 'cc_sampling_rate'))
-    start, end, datelist = build_movstack_datelist(db)
+    start, end, datelist = build_movstack_datelist(db, startdate, enddate)
     base = mdates.date2num(start) 
     plt.figure(figsize=(12, 9))
     sta1 = sta1.replace('.', '_')

msnoise/s000installer.py

@@ -30,7 +30,6 @@
     directory. It's not very safe, but until now we haven't thought of another
     solution.
 """
-import argparse
 import sys
 from getpass import getpass