small changes in the plots and avoiding to loop for filters (could fa…

…il if filter doesn't exist, e.g when moved filter 01 to 91)

msnoise/api.py

@@ -303,7 +303,7 @@ def get_params(session):
 # FILTERS PART
 
 
-def get_filters(session, all=False):
+def get_filters(session, all=False, ref=None):
     """Get Filters from the database.
 
     :type session: :class:`sqlalchemy.orm.session.Session`
@@ -317,7 +317,9 @@ def get_filters(session, all=False):
     :rtype: list of :class:`~msnoise.msnoise_table_def.declare_tables.Filter`
     :returns: a list of Filter
     """
-
+    if ref:
+        filter = session.query(Filter).filter(Filter.ref == ref).first()
+        return filter
     if all:
         filters = session.query(Filter).all()
     else:

msnoise/plots/ccftime.py

@@ -104,12 +104,10 @@ def main(sta1, sta2, filterid, components, mov_stack=1, ampli=5, seismic=False,
             y2 = line*ampli + i
             plt.fill_between(t, y1, y2, where=y2 >= y1, facecolor='k',
                              interpolate=True)
-    low = high = 0.0
-    for filterdb in get_filters(db, all=True):
-        if filterid == filterdb.ref:
-            low = float(filterdb.low)
-            high = float(filterdb.high)
-            break
+
+    filter = get_filters(db, ref=filterid)
+    low = float(filter.low)
+    high = float(filter.high)
 
     plt.xlabel("Lag Time (s)")
     plt.axhline(0, lw=0.5, c='k')

msnoise/plots/distance.py

@@ -76,12 +76,9 @@ def main(filterid, components, ampli=1, show=True, outfile=None,
 
     plt.ylabel("Interstation Distance in km")
     plt.xlabel("Lag Time")
-    low = high = 0.0
-    for filterdb in get_filters(db, all=True):
-        if filterid == filterdb.ref:
-            low = float(filterdb.low)
-            high = float(filterdb.high)
-            break
+    filter = get_filters(db, ref=filterid)
+    low = float(filter.low)
+    high = float(filter.high)
     title = '%s, Filter %d (%.2f - %.2f Hz)' % \
             (components, filterid, low, high,)
     if refilter:

msnoise/plots/dvv.py

@@ -13,39 +13,12 @@
 
 """
 
-import matplotlib.gridspec as gridspec
 import matplotlib.pyplot as plt
-
 from matplotlib.dates import DateFormatter
 
 from ..api import *
 
 
-def wavg(group, dttname, errname):
-    d = group[dttname]
-    group[errname][group[errname] == 0] = 1e-6
-    w = 1. / group[errname]
-    wavg = (d * w).sum() / w.sum()
-    return wavg
-
-
-def wstd(group, dttname, errname):
-    d = group[dttname]
-    group[errname][group[errname] == 0] = 1e-6
-    w = 1. / group[errname]
-    wavg = (d * w).sum() / w.sum()
-    N = len(np.nonzero(w)[0])
-    wstd = np.sqrt(np.sum(w * (d - wavg) ** 2) / ((N - 1) * np.sum(w) / N))
-    return wstd
-
-
-def get_wavgwstd(data, dttname, errname):
-    grouped = data.groupby(level=0)
-    g = grouped.apply(wavg, dttname=dttname, errname=errname)
-    h = grouped.apply(wstd, dttname=dttname, errname=errname)
-    return g, h
-
-
 def main(mov_stack=None, dttname="M", components='ZZ', filterid=1,
          pairs=[], showALL=False, show=False, outfile=None):
     db = connect()
@@ -63,11 +36,9 @@ def main(mov_stack=None, dttname="M", components='ZZ', filterid=1,
         components = [components, ]
 
     low = high = 0.0
-    for filterdb in get_filters(db, all=True):
-        if filterid == filterdb.ref:
-            low = float(filterdb.low)
-            high = float(filterdb.high)
-            break
+    filter = get_filters(db, ref=filterid)
+    low = float(filter.low)
+    high = float(filter.high)
 
     fig, axes = plt.subplots(len(mov_stacks), 1, sharex=True, figsize=(12, 9))
 

msnoise/plots/interferogram.py

@@ -75,11 +75,9 @@ def main(sta1, sta2, filterid, components, mov_stack=1, show=True,
     # ax.xaxis.set_minor_locator(DayLocator())
     # ax.xaxis.set_minor_formatter(DateFormatter('%Y-%m-%d %H:%M'))
 
-    for filterdb in get_filters(db, all=True):
-        if filterid == filterdb.ref:
-            low = float(filterdb.low)
-            high = float(filterdb.high)
-            break
+    filter = get_filters(db, ref=filterid)
+    low = float(filter.low)
+    high = float(filter.high)
 
     if "ylim" in kwargs:
         plt.ylim(kwargs["ylim"][0],kwargs["ylim"][1])

msnoise/plots/spectime.py

@@ -44,7 +44,8 @@
 from obspy.signal.filter import bandpass
 
 from msnoise.api import build_movstack_datelist, connect, get_config, \
-    get_filters, get_results, check_stations_uniqueness
+    get_filters, get_results, check_stations_uniqueness, xr_get_ccf,\
+    get_t_axis
 
 
 def main(sta1, sta2, filterid, components, mov_stack=1, ampli=5, show=False,
@@ -54,7 +55,7 @@ def main(sta1, sta2, filterid, components, mov_stack=1, ampli=5, show=False,
     cc_sampling_rate = float(get_config(db, 'cc_sampling_rate'))
     start, end, datelist = build_movstack_datelist(db)
     base = mdates.date2num(start)
-
+    taxis = get_t_axis(db)
     # TODO: Height adjustment of the plot for large number of stacks.
     # Preferably interactive
     fig = plt.figure(figsize=(12, 9))
@@ -75,10 +76,16 @@ def main(sta1, sta2, filterid, components, mov_stack=1, ampli=5, show=False,
 
     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):
+    stack_total = xr_get_ccf(sta1, sta2, components, filterid, mov_stack, taxis)
+
+    # convert index to mdates
+    stack_total.index = mdates.date2num(stack_total.index.to_pydatetime())
+
+    if len(stack_total) == 0:
+        print("No CCF found for this request")
+        return
+    ax = plt.subplot(111)
+    for i, line in stack_total.iterrows():
         if np.all(np.isnan(line)):
             continue
 
@@ -89,13 +96,12 @@ def main(sta1, sta2, filterid, components, mov_stack=1, ampli=5, show=False,
         freq, line = prepare_abs_postitive_fft(line, cc_sampling_rate)
         line /= line.max()
 
-        ax.plot(freq, line * ampli + i + base, c='k', lw=1)
+        ax.plot(freq, line * ampli + i, c='k', lw=1)
+
+    filter = get_filters(db, ref=filterid)
+    low = float(filter.low)
+    high = float(filter.high)
 
-    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))

msnoise/s04_stack2.py

@@ -173,7 +173,7 @@ def main(stype, interval=1.0, loglevel="INFO"):
                     if mov_stack > len(dr.times):
                         print("not enough data for mov_stack=%i" % mov_stack)
                         continue
-
+                    # TODO avoid the 1D resampler here
                     xx = dr.resample(times='1D').mean().rolling(
                         times=mov_stack, min_periods=1).mean().dropna("times", how="all")
                     xr_save_ccf(sta1, sta2, components, filterid, mov_stack, taxis, xx, overwrite=True)

msnoise/s05compute_mwcs2.py

@@ -129,17 +129,6 @@ def main(loglevel="INFO"):
     goal_sampling_rate = params.cc_sampling_rate
     maxlag = params.maxlag
 
-    # First we reset all DTT jobs to "T"odo if the REF is new for a given pair
-    # for station1, station2 in get_station_pairs(db, used=True):
-    #     sta1 = "%s.%s" % (station1.net, station1.sta)
-    #     sta2 = "%s.%s" % (station2.net, station2.sta)
-    #     pair = "%s:%s" % (sta1, sta2)
-    #     if is_dtt_next_job(db, jobtype='DTT', ref=pair):
-    #         logger.info(
-    #             "We will recompute all MWCS based on the new REF for %s" % pair)
-    #         reset_dtt_jobs(db, pair)
-    #         update_job(db, "REF", pair, jobtype='DTT', flag='D')
-    # 
     logger.debug('Ready to compute')
     # Then we compute the jobs
     outfolders = []