Squirrel API usage

Variant 1: quick and dirty script

Get hourly RMS values for a seismic recording

import numpy as np
from pyrocko import progress
from pyrocko.util import setup_logging, time_to_str, str_to_time_fillup as stt
from pyrocko.squirrel import Squirrel


setup_logging('squirrel_rms1.py', 'info')


def rms(data):
    return np.sqrt(np.sum(data**2))


# Time span (Hunga Tonga explosion was on 2022-01-15):
tmin = stt('2022-01-14')
tmax = stt('2022-01-17')

# Filter range:
fmin = 0.01
fmax = 0.05

# Enable progress bars:
progress.set_default_viewer('log')

# All data access will happen through a single Squirrel instance:
sq = Squirrel()

# Add local data directories:
# sq.add('data/2022')

# Add online data source:
# Enable access to BGR's FDSN web service and restrict to GR network and
# LH channels only:
sq.add_fdsn('bgr', query_args=dict(network='GR', channel='LH?'))

# Ensure meta-data is up to date for the selected time span:
sq.update(tmin=tmin, tmax=tmax)

# Allow waveform download for station BFO. This does not download anything
# yet, it just enables downloads later, when needed. Omit `tmin`, `tmax`,
# and `codes` to enable download of all everything.
sq.update_waveform_promises(tmin=tmin, tmax=tmax, codes='*.BFO.*.*')

# Iterate window-wise, with some overlap over the data:
for batch in sq.chopper_waveforms(
        tmin=tmin,
        tmax=tmax,
        codes='*.*.*.LHZ',
        tinc=3600.,              # One hour time windows (payload).
        tpad=1.0/fmin,           # Add padding to absorb filter artifacts.
        want_incomplete=False,   # Skip incomplete windows.
        snap_window=True):       # Start all windows at full hours.

    for tr in batch.traces:

        # Filtering will introduce some artifacts in the padding interval.
        tr.highpass(4, fmin)
        tr.lowpass(4, fmax)

        # Cut off the contaminated padding. Trim to the payload interval.
        tr.chop(batch.tmin, batch.tmax)

        # Print channel codes, time-mark and RMS value of the hour.
        print(tr.str_codes, time_to_str(tr.tmin), rms(tr.ydata))

Variant 2: command line app (simple)

import numpy as np
from pyrocko.util import time_to_str as tts
from pyrocko import squirrel


def rms(data):
    return np.sqrt(np.sum(data**2))


sq = squirrel.from_command(
    description='Report hourly RMS values.')

fmin = 0.01
fmax = 0.05

for batch in sq.chopper_waveforms(
        tinc=3600.,
        tpad=1.0/fmin,
        want_incomplete=False,
        snap_window=True):

    for tr in batch.traces:
        tr.highpass(4, fmin)
        tr.lowpass(4, fmax)
        tr.chop(batch.tmin, batch.tmax)
        print(tr.str_codes, tts(tr.tmin), rms(tr.ydata))

Variant 3a: command line app (tight integration, single-command)

import numpy as np
from pyrocko.util import time_to_str as tts
from pyrocko import squirrel


def rms(data):
    return np.sqrt(np.sum(data**2))


class RMSTool(squirrel.SquirrelCommand):

    def setup_subcommand(self, subparsers):
        return self.add_parser(
            subparsers, 'rms', help='Report hourly RMS values.')

    def setup(self, parser):
        self.add_selection_arguments(parser)

        parser.add_argument(
            '--fmin',
            dest='fmin',
            metavar='FLOAT',
            type=float,
            help='Corner of highpass [Hz].')

        parser.add_argument(
            '--fmax',
            dest='fmax',
            metavar='FLOAT',
            type=float,
            help='Corner of lowpass [Hz].')

    def call(self, parser, args):
        sq = self.squirrel_from_selection_arguments(args)

        fmin = args.fmin
        fmax = args.fmax

        for batch in sq.chopper_waveforms(
                tinc=3600.,
                tpad=1.0/fmin if fmin is not None else 0.0,
                want_incomplete=False,
                snap_window=True):

            for tr in batch.traces:

                if fmin is not None:
                    tr.highpass(4, fmin)

                if fmax is not None:
                    tr.lowpass(4, fmax)

                tr.chop(batch.tmin, batch.tmax)
                print(tr.str_codes, tts(tr.tmin), rms(tr.ydata))


squirrel.from_command(
    RMSTool(),
    description='Report hourly RMS values.')

Variant 3b: command line app (tight intergation, multi-command)

import numpy as np
from pyrocko.util import time_to_str as tts
from pyrocko import squirrel


def rms(data):
    return np.sqrt(np.sum(data**2))


class RMSTool(squirrel.SquirrelCommand):

    def setup_subcommand(self, subparsers):
        return self.add_parser(
            subparsers, 'rms', help='Report hourly RMS values.')

    def setup(self, parser):
        self.add_selection_arguments(parser)

        parser.add_argument(
            '--fmin',
            dest='fmin',
            metavar='FLOAT',
            type=float,
            help='Corner of highpass [Hz].')

        parser.add_argument(
            '--fmax',
            dest='fmax',
            metavar='FLOAT',
            type=float,
            help='Corner of lowpass [Hz].')

    def call(self, parser, args):
        sq = self.squirrel_from_selection_arguments(args)

        fmin = args.fmin
        fmax = args.fmax

        for batch in sq.chopper_waveforms(
                tinc=3600.,
                tpad=1.0/fmin if fmin is not None else 0.0,
                want_incomplete=False,
                snap_window=True):

            for tr in batch.traces:

                if fmin is not None:
                    tr.highpass(4, fmin)

                if fmax is not None:
                    tr.lowpass(4, fmax)

                tr.chop(batch.tmin, batch.tmax)
                print(tr.str_codes, tts(tr.tmin), rms(tr.ydata))


class MinMaxTool(squirrel.SquirrelCommand):

    def setup_subcommand(self, subparsers):
        return self.add_parser(
            subparsers, 'minmax', help='Print hourly min/max values.')

    def call(self, parser, args):
        self.fail('Not implemented yet!')


squirrel.from_command(
    subcommands=[RMSTool(), MinMaxTool()],
    description='My favourite tools.')

Variant 4: command line app (loose integration)

import argparse
import numpy as np
from pyrocko.util import time_to_str as tts
from pyrocko.squirrel.tool import common as squirrel_tool_common


def rms(data):
    return np.sqrt(np.sum(data**2))


parser = argparse.ArgumentParser()
squirrel_tool_common.add_selection_arguments(parser)

parser.add_argument(
    '--fmin',
    dest='fmin',
    metavar='FLOAT',
    type=float,
    help='Corner of highpass [Hz].')

parser.add_argument(
    '--fmax',
    dest='fmax',
    metavar='FLOAT',
    type=float,
    help='Corner of lowpass [Hz].')

args = parser.parse_args()

sq = squirrel_tool_common.squirrel_from_selection_arguments(args)

fmin = args.fmin
fmax = args.fmax

for batch in sq.chopper_waveforms(
        tinc=3600.,
        tpad=1.0/fmin if fmin is not None else 0.0,
        want_incomplete=False,
        snap_window=True):

    for tr in batch.traces:

        if fmin is not None:
            tr.highpass(4, fmin)

        if fmax is not None:
            tr.lowpass(4, fmax)

        tr.chop(batch.tmin, batch.tmax)
        print(tr.str_codes, tts(tr.tmin), rms(tr.ydata))