rough script to make two helicorder and spectrogram plots, one each unfiltered and filtered

obspy-day_spec_plots.py

from obspy.core.event.catalog import Catalog
from datetime import datetime, timedelta
import pytz
import os
import shutil
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import numpy as np
from obspy import read
from obspy.core.utcdatetime import UTCDateTime
from obspy import read_inventory, read_events
from obspy.clients.fdsn.client import Client
from obspy.clients.fdsn.header import FDSNException

"""
This script is designed to work with a directory structure like this with sorted miniSEED files (NET/STA/NET.STA.LOC.CHA.D.YYYY.JJJ).
All miniSEED files must be in the directory tree below the INPATH variable.
INPATH directory structure should be as follows:

INPATH/
├── AM # <- network folder
│   └── R35E7 # <- station folder
│       ├── AM.R35E7.00.SHZ.D.2018.174 # <- miniSEED files
│       ├── AM.R35E7.00.SHZ.D.2018.175
│       ├── AM.R35E7.00.SHZ.D.2018.176
│       ├── AM.R35E7.00.SHZ.D.2018.177
│       └── ... etc
├── events
│   └── EVENTFILE.xml # <- this XML event catalog file is generated by another script. if it's not there, this script
│                     # downloads some recent events from CLIENT both globally above M6.0 and locally near LAT and LON
├── output # <- location of PNG image outputs
│   ├── AM.R35E7.00.SHZ.D.2019.174-heli.png
│   ├── AM.R35E7.00.SHZ.D.2019.174-helibp.png
│   ├── AM.R35E7.00.SHZ.D.2019.174-spec.png
│   ├── AM.R35E7.00.SHZ.D.2019.174-specbp.png
│   └── ... etc
└── proc # <- processing directory for temporary storage of miniSEED files (prevents write collisions with Earthworm)
"""

# station name
STA = 'R0000'
CHN = 'EHZ'
# latitude and longitude
LAT = 45.0
LON = -69.0
# client for catalog download (for options see https://docs.obspy.org/packages/obspy.clients.fdsn.html#basic-fdsn-client-usage)
CLIENT = 'USGS'
# time zone database name (see https://en.wikipedia.org/wiki/List_of_tz_database_time_zones)
TZ = 'America/New_York'

INPATH = '/home/user/seismic/' # the root data directory
OUTPATH = os.path.join(INPATH, 'output') # output directory
EVENTPATH = os.path.join(INPATH, 'events') # directory of events catalog to read
EVENTFILE = 'recentevents.xml' # name of event catalog file as noted above
STA_DEF = 'AM.%s.00.%s.D' % (STA,CHN) # station name string

def processor(start, end, sta=STA_DEF, dbscale=False, filtmin=0.7, filtmax=2.0, inpath=INPATH, outpath=OUTPATH):
    """
    Inputs:
    start: obspy.core.utcdatetime.UTCDateTime (start of spectrogram)
    end: obspy.core.utcdatetime.UTCDateTime (end of spectrogram)
    sta: string ('AM.Rxxxx.00.xxx.D.')
    dbscale: boolean (False: normal spectrogram colormap; True: log colormap. Default: False)
    filtmin: float or int (high-pass filter value, Hz)
    filtmax: float or int (low-pass filter value, Hz)
    inpath: os.path or string (the location where miniSEED files are stored)
    outpath: os.path or string (the location where plots are generated)

    Returns:
    dict{
        'filtmax': float (low-pass filter value, Hz),
        'filtmin': float (high-pass filter value, Hz),
        'heli': str (location of helicorder plot on disk),
        'helibp': str (location of bandpassed helicorder plot on disk),
        'spec': str (location of spectrogram plot on disk),
        'specbp': str (location of bandpassed spectrogram plot on disk)
    }
    """
    day = start.strftime('%Y.%j')
    yday = (start - timedelta(days=1)).strftime('%Y.%j')
    daystart = UTCDateTime(start.year, start.month, start.day)
    dayend = daystart + timedelta(days=1)
    if dayend > datetime.now():
        now = UTCDateTime.now()
        mins = 0
        hourdelta = timedelta(hours=0)
        if 14 >= now.minute >= 0:
            mins = 15
        elif 29 >= now.minute >= 15:
            mins = 30
        elif 44 >= now.minute >= 30:
            mins = 45
        else:
            mins = 0
            hourdelta = timedelta(hours=1)
        now += hourdelta
        dayend = UTCDateTime(now.year, now.month, now.day, now.hour, mins)
        daystart = dayend - timedelta(days=1)

    if sta:
        stn = sta
        #sta = sta + '.D.'
    else:
        stn = str(STA_DEF[0:-2])
        #sta = STA_DEF
        sta = stn

    stc = stn.split('.')
    net = stc[0]
    sta = stc[1]
    loc = stc[2]
    ch = stc[3]

    fn = '%s.%s.%s.%s.%s' % (sta, net, loc, ch, day)
    yfn = '%s.%s.%s.%s.%s' % (sta, net, loc, ch, yday)

    inpath = os.path.join(inpath, stc[0], stc[1])
    if os.path.isdir(os.path.join(inpath, 'proc')):
        pass
    else:
        os.mkdir(os.path.join(inpath, 'proc'))

    shutil.copy2(os.path.join(inpath, fn), os.path.join(inpath, 'proc'))
    shutil.copy2(os.path.join(inpath, yfn), os.path.join(inpath, 'proc'))
    ypath = inpath
    inpath = os.path.join(inpath, 'proc')


    tz = int(datetime.now(pytz.timezone(TZ)).strftime('%z'))/100
    fmin = 0.1
    fmax = 25
    fminbp = filtmin
    fmaxbp = filtmax

    heli = os.path.join(outpath, stn + '.' + day + '-heli.png')
    helibp = os.path.join(outpath, stn + '.' + day + '-heli-band.png')
    dur, spec = '', ''

    st = read().clear()
    yst = st.copy()
    try:
        yst = read(os.path.join(ypath, yfn))
    except:
        print("error reading yesterday's miniSEED file. may be further errors...")

    try:
        st = read(os.path.join(inpath, fn))
        os.remove(os.path.join(inpath, fn)) # remove the temporary file
    except:
        print("error reading today's miniSEED file. may be further errors...")

    net = str(st[0].stats.network)
    sta = str(st[0].stats.station)
    loc = str(st[0].stats.location)
    ch = str(st[0].stats.channel)
    startt = str(st[0].stats.starttime)
    sr = str(st[0].stats.sampling_rate)

    try:
        st = yst + st
    except:
        pass
    #st.merge()
    try:
        st = st.slice(starttime=daystart, endtime=dayend)
    except:
        st = yst
        st = st.slice(starttime=daystart, endtime=dayend)

    sbp = st.copy()
    sbp = sbp.filter('bandpass', freqmin=fminbp, freqmax=fmaxbp, zerophase=True)
    spu = st.slice(starttime=start, endtime=end)
    sps = sbp.slice(starttime=start, endtime=end) # slice for bandpass spectrogram

    cat = Catalog()
    try:
        cat.extend(read_events(pathname_or_url=os.path.join(EVENTPATH, EVENTFILE), format='QUAKEML'))
    except:
        print('no events XML found. trying to download...')
    # get events
    client = Client(CLIENT)
    cat = Catalog()
    try:
        cat += client.get_events(starttime=daystart, endtime=dayend, latitude=LAT, longitude=LON, maxradius=10)
    except FDSNException:
        print('no local events in %s catalog' % CLIENT)
    try:
        cat += client.get_events(starttime=daystart, endtime=dayend, latitude=LAT, longitude=LON,
                    minradius=10, maxradius=15, minmagnitude=2.5)
    except FDSNException:
        print('no local events in %s catalog' % CLIENT)
    try:
        cat += client.get_events(starttime=daystart, endtime=dayend, minradius=10, minmagnitude=6.0)
    except FDSNException:
        print('no local events in %s catalog' % CLIENT)

    title = net + '.' + sta + '.' + loc + '.' + ch + ' - ' + startt + ' - rate: ' + sr

    st.plot(type="dayplot", size=(1600, 1200), title=title + 'Hz - band: 0-25Hz', vertical_scaling_range=2000,
        tick_format='%H:%M', outfile=heli, color=['k', 'r', 'b', 'g'], linewidth=0.3, time_offset=tz, events=cat)
    sbp.plot(type="dayplot", size=(1600, 1200), title=title + 'Hz - band: '+ str(fminbp) + '-' + str(fmaxbp) + 'Hz', vertical_scaling_range=200,
        tick_format='%H:%M', outfile=helibp, color=['k', 'r', 'b', 'g'], linewidth=0.3, time_offset=tz, events=cat)

    #st.plot(type="dayplot", title=net + '.' + sta + '.' + loc + '.' + ch + ' - ' + startt + ' - rate: ' + sr + 'Hz - band: 0-25Hz', vertical_scaling_range=8e3, outfile=heli, color=['k', 'r', 'b', 'g'], time_offset=tz, events={'min_magnitude': 6.5})
    #sbp.plot(type="dayplot", title=net + '.' + sta + '.' + loc + '.' + ch + ' - ' + startt + ' - rate: ' + sr + 'Hz - band: '+ str(fminbp) + '-' + str(fmaxbp) + 'Hz', vertical_scaling_range=7e2, outfile=helibp, color=['k', 'r', 'b', 'g'], time_offset=tz, events={'min_magnitude': 6.5})

    heli = outpath + os.path.split(heli)[1]
    helibp = outpath + os.path.split(helibp)[1]

    if end:
        dur = end - start


    sp = spu.detrend(type='constant')
    ss = sps.detrend(type='constant')

    startt = str(sp[0].stats.starttime)


    ## ------------------------- ##
    # make spectrogram figure 1
    fig = plt.figure(figsize=(16,6), dpi=100)
    ax1 = fig.add_axes([0.068, 0.75, 0.85, 0.2]) #[left bottom width height]
    ax2 = fig.add_axes([0.068, 0.1, 0.85, 0.6], sharex=ax1)
    ax3 = fig.add_axes([0.931, 0.1, 0.03, 0.6])

    # labels
    fig.suptitle(net + '.' + sta + '.' + loc + '.' + ch + ' - ' + startt + ' - samplerate: ' + sr + 'Hz - frequency band: 0-25 Hz')
    ax1.set_ylabel('Traces')
    ax2.set_xlabel('Time [s]')
    ax2.set_ylabel('Frequency [Hz]')
    ax3.set_ylabel('Energy density [dimensionless]') # doesn't work

    # make time vector
    t = np.arange(sp[0].stats.npts) / sp[0].stats.sampling_rate

    # plot waveform (top subfigure)
    ax1.plot(t, sp[0].data, 'k', linewidth=0.5)

    # plot spectrogram (bottom subfigure)
    fig = sp[0].spectrogram(show=False, axes=ax2, log=False, dbscale=dbscale, cmap='viridis')
    mappable = ax2.images[0]
    plt.colorbar(mappable=mappable, cax=ax3)

    ax2.set_ylim(fmin, fmax)

    spec = os.path.join(outpath, stn + '.' + start.strftime('%Y.%j') + "-spec.png")
    plt.savefig(spec)
    spec = outpath + os.path.split(spec)[1]


    ## ------------------------- ##
    # make spectrogram figure 2
    sfig2 = plt.figure(figsize=(16,4), dpi=100)
    ax1 = sfig2.add_axes([0.068, 0.600, 0.85, 0.3]) #[left bottom width height]
    ax2 = sfig2.add_axes([0.068, 0.115, 0.85, 0.4], sharex=ax1)
    ax3 = sfig2.add_axes([0.932, 0.115, 0.03, 0.4])

    # labels
    sfig2.suptitle(net + '.' + sta + '.' + loc + '.' + ch + ' - ' + startt + ' - samplerate: ' + sr + 'Hz - bandpass: ' + str(fminbp) + '-' + str(fmaxbp) + ' Hz')
    ax1.set_ylabel('Counts')
    ax2.set_xlabel('Time [s]')
    ax2.set_ylabel('Frequency [Hz]')
    ax3.set_ylabel('Energy density [dimensionless]') # doesn't work

    # make time vector
    t = np.arange(ss[0].stats.npts) / ss[0].stats.sampling_rate

    # plot waveform (top subfigure)
    ax1.plot(t, ss[0].data, 'k', linewidth=0.5)

    # plot spectrogram (bottom subfigure)
    sfig2 = ss[0].spectrogram(show=False, axes=ax2, log=False, dbscale=dbscale, cmap='viridis')
    mappable = ax2.images[0]
    plt.colorbar(mappable=mappable, cax=ax3)

    ax2.set_ylim(fminbp, fmaxbp)


    specbp = os.path.join(outpath, stn + '.' + start.strftime('%Y.%j') + "-spec-band.png")
    plt.savefig(specbp)
    specbp = outpath + os.path.split(specbp)[1]



    returns = {
        'filtmin': fminbp,
        'filtmax': fmaxbp,
        'heli': heli,
        'helibp': helibp,
        'spec': spec,
        'specbp': specbp,
    }
    return returns