Download data from the seismological resarch network Rhineland lower palatinate

seigerdown.py

#!/usr/bin/env pyrocko-python

import sys
import logging
import tempfile
import math
import os.path as op
import shutil
try:
    from urllib.error import HTTPError
except:
    from urllib2 import HTTPError
import glob
import pipes
from optparse import OptionParser
from collections import defaultdict

import numpy as num

from pyrocko import trace, util, io, cake, catalog, automap, pile, model
from pyrocko import orthodrome, weeding
from pyrocko.client import fdsn
from pyrocko.io import resp, enhanced_sacpz as epz, stationxml
import subprocess

km = 1000.

g_sites_available = sorted(fdsn.g_site_abbr.keys())

geofon = catalog.Geofon()
usgs = catalog.USGS(catalog=None)

tfade_factor = 1.0
ffade_factors = 0.5, 1.5

fdsn.g_timeout = 60.


class starfill(object):
    def __getitem__(self, k):
        return '*'


def nice_seconds_floor(s):
    nice = [1., 10., 60., 600., 3600., 3.*3600., 12*3600., 24*3600., 48*3600.]
    p = s
    for x in nice:
        if s < x:
            return p

        p = x

    return s


def get_events(time_range, region=None, catalog=geofon, **kwargs):
    if not region:
        return catalog.get_events(time_range, **kwargs)

    events = []
    for (west, east, south, north) in automap.split_region(region):
        events.extend(
            catalog.get_events(
                time_range=time_range,
                lonmin=west,
                lonmax=east,
                latmin=south,
                latmax=north, **kwargs))

    return events


def cut_n_dump(traces, win, out_path):
    otraces = []
    for tr in traces:
        try:
            otr = tr.chop(win[0], win[1], inplace=False)
            otraces.append(otr)
        except trace.NoData:
            pass

    return io.save(otraces, out_path)


aliases = {
    '2013_landau_1': '2013-04-07 07:59:00',
}


def get_events_by_name_or_date(event_names_or_dates, catalog=geofon):
    stimes = []
    for sev in event_names_or_dates:
        if sev in aliases:
            if isinstance(aliases[sev], str):
                stimes.append(aliases[sev])
            else:
                stimes.extend(aliases[sev])
        else:
            stimes.append(sev)

    events_out = []
    for stime in stimes:
        if op.isfile(stime):
            events_out.extend(model.Event.load_catalog(stime))
        elif stime.startswith('gfz'):
            event = geofon.get_event(stime)
            events_out.append(event)
        else:
            t = util.str_to_time(stime)
            events = get_events(time_range=(t-60., t+60.), catalog=catalog)
            events.sort(key=lambda ev: abs(ev.time - t))
            event = events[0]
            events_out.append(event)

    return events_out


class NoArrival(Exception):
    pass


class PhaseWindow(object):

    def __init__(self, model, phases, omin, omax):
        self.model = model
        self.phases = phases
        self.omin = omin
        self.omax = omax

    def __call__(self, time, distance, depth):
        for ray in self.model.arrivals(
                phases=self.phases,
                zstart=depth,
                distances=[distance*cake.m2d]):

            return time + ray.t + self.omin, time + ray.t + self.omax

        raise NoArrival


class VelocityWindow(object):
    def __init__(self, vmin, vmax=None, tpad=0.0):
        self.vmin = vmin
        self.vmax = vmax
        self.tpad = tpad

    def __call__(self, time, distance, depth):
        ttmax = (depth + distance) / self.vmin
        if self.vmax is not None:
            ttmin = (depth + distance) / self.vmax
        else:
            ttmin = 0.0

        return time + ttmin - self.tpad, time + ttmax + self.tpad


class FixedWindow(object):
    def __init__(self, tmin, tmax):
        self.tmin = tmin
        self.tmax = tmax

    def __call__(self, time, distance, depth):
        return self.tmin, self.tmax


def dump_commandline(argv, fn):
    s = ' '.join([pipes.quote(x) for x in argv])
    with open(fn, 'w') as f:
        f.write(s)
        f.write('\n')


g_user_credentials = {}
g_auth_tokens = {}


def get_user_credentials(site):
    user, passwd = g_user_credentials.get(site, (None, None))
    token = g_auth_tokens.get(site, None)
    return dict(user=user, passwd=passwd, token=token)


program_name = 'seigerdown'
description = '''
Download waveforms from FDSN web services and prepare for seiger
'''.strip()

logger = logging.getLogger('')

usage = '''
usage: seigerdown [options] [--] <YYYY-MM-DD> <HH:MM:SS> <lat> <lon> \\
                               <depth_km> <radius_km> <fmin_hz> \\
                               <sampling_rate_hz> \\
                               <eventname>

       seigerdown [options] [--] <YYYY-MM-DD> <HH:MM:SS> <radius_km> <fmin_hz> \\
                               <sampling_rate_hz> <eventname>

       seigerdown [options] [--] <catalog-eventname> <radius_km> <fmin_hz> \\
                               <sampling_rate_hz> <eventname>

       seigerdown [options] --window="<YYYY-MM-DD HH:MM:SS, YYYY-MM-DD HH:MM:\
SS>" \\
                               [--] <lat> <lon> <radius_km> <fmin_hz> \\
                               <sampling_rate_hz> <eventname>
'''.strip()


def main():
    parser = OptionParser(
        usage=usage,
        description=description)

    parser.add_option(
        '--force',
        dest='force',
        action='store_true',
        default=False,
        help='allow recreation of output <directory>')

    parser.add_option(
        '--debug',
        dest='debug',
        action='store_true',
        default=False,
        help='print debugging information to stderr')

    parser.add_option(
        '--dry-run',
        dest='dry_run',
        action='store_true',
        default=False,
        help='show available stations/channels and exit '
             '(do not download waveforms)')

    parser.add_option(
        '--continue',
        dest='continue_',
        action='store_true',
        default=False,
        help='continue download after a accident')

    parser.add_option(
        '--local-data',
        dest='local_data',
        action='append',
        help='add file/directory with local data')

    parser.add_option(
        '--local-stations',
        dest='local_stations',
        action='append',
        help='add local stations file')

    parser.add_option(
        '--local-responses-resp',
        dest='local_responses_resp',
        action='append',
        help='add file/directory with local responses in RESP format')

    parser.add_option(
        '--local-responses-pz',
        dest='local_responses_pz',
        action='append',
        help='add file/directory with local pole-zero responses')

    parser.add_option(
        '--local-responses-stationxml',
        dest='local_responses_stationxml',
        help='add file with local response information in StationXML format')

    parser.add_option(
        '--window',
        dest='window',
        default='full',
        help='set time window to choose [full, p, "<time-start>,<time-end>"'
             '] (time format is YYYY-MM-DD HH:MM:SS)')

    parser.add_option(
        '--out-components',
        choices=['enu', 'rtu'],
        dest='out_components',
        default='rtu',
        help='set output component orientations to radial-transverse-up [rtu] '
             '(default) or east-north-up [enu]')

    parser.add_option(
        '--padding-factor',
        type=float,
        default=3.0,
        dest='padding_factor',
        help='extend time window on either side, in multiples of 1/<fmin_hz> '
             '(default: 5)')

    parser.add_option(
        '--credentials',
        dest='user_credentials',
        action='append',
        default=[],
        metavar='SITE,USER,PASSWD',
        help='user credentials for specific site to access restricted data '
             '(this option can be repeated)')

    parser.add_option(
        '--token',
        dest='auth_tokens',
        metavar='SITE,FILENAME',
        action='append',
        default=[],
        help='user authentication token for specific site to access '
             'restricted data (this option can be repeated)')

    parser.add_option(
        '--sites',
        dest='sites',
        metavar='SITE1,SITE2,...',
        default='http://ws.gpi.kit.edu,bgr',
        help='sites to query (available: %s, default: "%%default"'
        % ', '.join(g_sites_available))

    parser.add_option(
        '--band-codes',
        dest='priority_band_code',
        metavar='V,L,M,B,H,S,E,...',
        default='V,L,M,B,H,E',
        help='select and prioritize band codes (default: %default)')

    parser.add_option(
        '--instrument-codes',
        dest='priority_instrument_code',
        metavar='H,L,G,...',
        default='H,L,O,',
        help='select and prioritize instrument codes (default: %default)')

    parser.add_option(
        '--radius-min',
        dest='radius_min',
        metavar='VALUE',
        default=0.0,
        type=float,
        help='minimum radius [km]')

    parser.add_option(
        '--tinc',
        dest='tinc',
        metavar='VALUE',
        default=3600.*12.,
        type=float,
        help='length of seperate saved files in s')

    parser.add_option(
        '--nstations-wanted',
        dest='nstations_wanted',
        metavar='N',
        type=int,
        help='number of stations to select initially')


    (options, args) = parser.parse_args(sys.argv[1:])
    if len(args) not in (9, 6, 5):
        parser.print_help()
        sys.exit(1)

    if options.debug:
        util.setup_logging(program_name, 'debug')
    else:
        util.setup_logging(program_name, 'info')

    if options.local_responses_pz and options.local_responses_resp:
        logger.critical('cannot use local responses in PZ and RESP '
                        'format at the same time')
        sys.exit(1)

    n_resp_opt = 0
    for resp_opt in (
            options.local_responses_pz,
            options.local_responses_resp,
            options.local_responses_stationxml):

        if resp_opt:
            n_resp_opt += 1

    if n_resp_opt > 1:
        logger.critical('can only handle local responses from either PZ or '
                        'RESP or StationXML. Cannot yet merge different '
                        'response formats.')
        sys.exit(1)

    if options.local_responses_resp and not options.local_stations:
        logger.critical('--local-responses-resp can only be used '
                        'when --stations is also given.')
        sys.exit(1)

    try:
        ename = ''
        magnitude = None
        mt = None
        if len(args) == 9:
            time = util.str_to_time(args[0] + ' ' + args[1])
            lat = float(args[2])
            lon = float(args[3])
            depth = float(args[4])*km
            iarg = 5

        elif len(args) == 6:
            if args[1].find(':') == -1:
                sname_or_date = None
                lat = float(args[0])
                lon = float(args[1])
                event = None
                time = None
            else:
                sname_or_date = args[0] + ' ' + args[1]

            iarg = 2

        elif len(args) == 5:
            sname_or_date = args[0]
            iarg = 1

        if len(args) in (6, 5) and sname_or_date is not None:
            events = get_events_by_name_or_date([sname_or_date],
                                                catalog=geofon)
            if len(events) == 0:
                logger.critical('no event found')
                sys.exit(1)
            elif len(events) > 1:
                logger.critical('more than one event found')
                sys.exit(1)

            event = events[0]
            time = event.time
            lat = event.lat
            lon = event.lon
            depth = event.depth
            ename = event.name
            magnitude = event.magnitude
            mt = event.moment_tensor

        radius = float(args[iarg])*km
        fmin = float(args[iarg+1])
        sample_rate = float(args[iarg+2])

        eventname = args[iarg+3]
        event_dir = op.join('data', 'events', eventname)
        output_dir = op.join(event_dir, 'waveforms')
    except:
        raise
        parser.print_help()
        sys.exit(1)

    if options.force and op.isdir(event_dir):
        if not options.continue_:
            shutil.rmtree(event_dir)

    if op.exists(event_dir) and not options.continue_:
        logger.critical(
            'directory "%s" exists. Delete it first or use the --force option'
            % event_dir)
        sys.exit(1)

    util.ensuredir(output_dir)

    if time is not None:
        event = model.Event(
            time=time, lat=lat, lon=lon, depth=depth, name=ename,
            magnitude=magnitude, moment_tensor=mt)

    if options.window == 'full':
        if event is None:
            logger.critical('need event for --window=full')
            sys.exit(1)

        low_velocity = 1500.
        timewindow = VelocityWindow(
            low_velocity, tpad=options.padding_factor/fmin)

        tmin, tmax = timewindow(time, radius, depth)

    elif options.window == 'p':
        if event is None:
            logger.critical('need event for --window=p')
            sys.exit(1)

        phases = list(map(cake.PhaseDef, 'P p'.split()))
        emod = cake.load_model()

        tpad = options.padding_factor / fmin
        timewindow = PhaseWindow(emod, phases, -tpad, tpad)

        arrivaltimes = []
        for dist in num.linspace(0, radius, 20):
            try:
                arrivaltimes.extend(timewindow(time, dist, depth))
            except NoArrival:
                pass

        if not arrivaltimes:
            logger.error('required phase arrival not found')
            sys.exit(1)

        tmin = min(arrivaltimes)
        tmax = max(arrivaltimes)

    else:
        try:
            stmin, stmax = options.window.split(',')
            tmin = util.str_to_time(stmin.strip())
            tmax = util.str_to_time(stmax.strip())

            timewindow = FixedWindow(tmin, tmax)

        except ValueError:
            logger.critical('invalid argument to --window: "%s"'
                            % options.window)
            sys.exit(1)

    if event is not None:
        event.name = eventname

    tlen = tmax - tmin
    tfade = tfade_factor / fmin

    tpad = tfade

    tmin -= tpad
    tmax += tpad

    priority_band_code = options.priority_band_code.split(',')
    for s in priority_band_code:
        if len(s) != 1:
            logger.critical('invalid band code: %s' % s)

    priority_instrument_code = options.priority_instrument_code.split(',')
    for s in priority_instrument_code:
        if len(s) != 1:
            logger.critical('invalid instrument code: %s' % s)

    station_query_conf = dict(
        latitude=lat,
        longitude=lon,
        minradius=options.radius_min*km*cake.m2d,
        maxradius=radius*cake.m2d,
        channel=','.join('?%s?' % s for s in priority_band_code))

    target_sample_rate = sample_rate

    fmax = target_sample_rate

    # target_sample_rate = None
    # priority_band_code = ['H', 'B', 'M', 'L', 'V', 'E', 'S']

    priority_units = ['M/S', 'M', 'M/S**2']

    output_units = 'M'

    sites = [x.strip() for x in options.sites.split(',') if x.strip()]
    tinc = options.tinc
  #  for site in sites:
   #     if site not in g_sites_available:
    #        logger.critical('unknown FDSN site: %s' % site)
     #       sys.exit(1)

    for s in options.user_credentials:
        try:
            site, user, passwd = s.split(',')
            g_user_credentials[site] = user, passwd
        except ValueError:
            logger.critical('invalid format for user credentials: "%s"' % s)
            sys.exit(1)

    for s in options.auth_tokens:
        try:
            site, token_filename = s.split(',')
            with open(token_filename, 'r') as f:
                g_auth_tokens[site] = f.read()
        except (ValueError, OSError, IOError):
            logger.critical('cannot get token from file: %s' % token_filename)
            sys.exit(1)

    fn_template0 = \
'data_%(network)s.%(station)s.%(location)s.%(channel)s_%(tmin)s.mseed'

    fn_template_raw = op.join(output_dir, 'raw',  fn_template0)
    fn_template_raw_folder = op.join(output_dir, 'raw/' , 'traces.mseed')
    fn_stations_raw = op.join(output_dir, 'stations.raw.txt')
    fn_template_rest = op.join(output_dir, 'rest',  fn_template0)
    fn_commandline = op.join(output_dir, 'seigerdown.command')

    ftap = (ffade_factors[0]*fmin, fmin, fmax, ffade_factors[1]*fmax)

    # chapter 1: download

    sxs = []
    for site in sites:
        try:
            extra_args = {
                'iris': dict(matchtimeseries=True),
            }.get(site, {})

            extra_args.update(station_query_conf)

            if site == 'geonet':
                extra_args.update(
                    starttime=tmin,
                    endtime=tmax)
            else:
                extra_args.update(
                    startbefore=tmax,
                    endafter=tmin,
                    includerestricted=(
                        site in g_user_credentials or site in g_auth_tokens))

            logger.info('downloading channel information (%s)' % site)
            sx = fdsn.station(
                site=site,
                format='text',
                level='channel',
                **extra_args)

        except fdsn.EmptyResult:
            logger.error('No stations matching given criteria. (%s)' % site)
            sx = None

        if sx is not None:
            sxs.append(sx)

    if all(sx is None for sx in sxs) and not options.local_data:
        sys.exit(1)

    nsl_to_sites = defaultdict(list)
    nsl_to_station = {}
    for sx, site in zip(sxs, sites):
        site_stations = sx.get_pyrocko_stations()
        for s in site_stations:
            nsl = s.nsl()
            nsl_to_sites[nsl].append(site)
            if nsl not in nsl_to_station:
                nsl_to_station[nsl] = s  # using first site with this station
    logger.info('number of stations found: %i' % len(nsl_to_station))

    # station weeding

    nsls_selected = None
    if options.nstations_wanted:
        stations_all = [
            nsl_to_station[nsl_] for nsl_ in sorted(nsl_to_station.keys())]

        for s in stations_all:
            s.set_event_relative_data(event)

        stations_selected = weeding.weed_stations(
            stations_all, options.nstations_wanted)[0]

        nsls_selected = set(s.nsl() for s in stations_selected)
        logger.info('number of stations selected: %i' % len(nsls_selected))

    have_data = set()

    if options.continue_:
        fns = glob.glob(fn_template_raw % starfill())
        p = pile.make_pile(fns)
    else:
        fns = []

    have_data_site = {}
    could_have_data_site = {}
    for site in sites:
        have_data_site[site] = set()
        could_have_data_site[site] = set()

    available_through = defaultdict(set)
    it = 0
    nt = int(math.ceil((tmax - tmin) / tinc))
    for it in range(nt):
        tmin_win = tmin + it * tinc
        tmax_win = min(tmin + (it + 1) * tinc, tmax)
        logger.info('time window %i/%i (%s - %s)' % (it+1, nt,
util.tts(tmin_win),
util.tts(tmax_win)))

        have_data_this_window = set()
        if options.continue_:
            trs_avail = p.all(tmin=tmin_win, tmax=tmax_win, load_data=False)
            for tr in trs_avail:
                have_data_this_window.add(tr.nslc_id)
        for site, sx in zip(sites, sxs):
            if sx is None:
                continue

            selection = []
            channels = sx.choose_channels(
                target_sample_rate=target_sample_rate,
                priority_band_code=priority_band_code,
                priority_units=priority_units,
                priority_instrument_code=priority_instrument_code,
                timespan=(tmin_win, tmax_win))

            for nslc in sorted(channels.keys()):
                if nsls_selected is not None and nslc[:3] not in nsls_selected:
                    continue

                could_have_data_site[site].add(nslc)

                if nslc not in have_data_this_window:
                    channel = channels[nslc]
                    if event:
                        lat_, lon_ = event.lat, event.lon
                    else:
                        lat_, lon_ = lat, lon

                    dist = orthodrome.distance_accurate50m_numpy(
                        lat_, lon_,
                        channel.latitude.value, channel.longitude.value)

                    if event:
                        depth_ = event.depth
                        time_ = event.time
                    else:
                        depth_ = None
                        time_ = None

                    tmin_, tmax_ = timewindow(time_, dist, depth_)

                    tmin_this = tmin_ - tpad
                    tmax_this = tmax_ + tpad

                    tmin_req = max(tmin_win, tmin_this)
                    tmax_req = min(tmax_win, tmax_this)

                    if channel.sample_rate:
                        deltat = 1.0 / channel.sample_rate.value
                    else:
                        deltat = 1.0

                    if tmin_req < tmax_req:
                        # extend time window by some samples because otherwise
                        # sometimes gaps are produced
                        selection.append(
                            nslc + (
                                tmin_req-deltat*10.0,
                                tmax_req+deltat*10.0))

            if options.dry_run:
                for (net, sta, loc, cha, tmin, tmax) in selection:
                    available_through[net, sta, loc, cha].add(site)

            else:
                neach = 100
                i = 0
                nbatches = ((len(selection)-1) // neach) + 1
                while i < len(selection):
                    selection_now = selection[i:i+neach]

                    f = tempfile.NamedTemporaryFile()
                    try:
                        sbatch = ''
                        if nbatches > 1:
                            sbatch = ' (batch %i/%i)' % (
                                (i//neach) + 1, nbatches)

                        logger.info('downloading data (%s)%s' % (site, sbatch))
                        data = fdsn.dataselect(
                            site=site, selection=selection_now,
                            **get_user_credentials(site))

                        while True:
                            buf = data.read(1024)
                            if not buf:
                                break
                            f.write(buf)

                        f.flush()

                        trs = io.load(f.name)
                        for tr in trs:
                            try:
                                tr.chop(tmin_win, tmax_win)
                                have_data.add(tr.nslc_id)
                                have_data_site[site].add(tr.nslc_id)
                            except trace.NoData:
                                pass

                        fns2 = io.save(trs, fn_template_raw)
                        io.save(trs, fn_template_raw_folder)
                        for fn in fns2:
                            if fn in fns:
                                logger.warn('overwriting file %s', fn)
                        fns.extend(fns2)

                    except fdsn.EmptyResult:
                        pass

                    except HTTPError:
                        logger.warn(
                            'an error occurred while downloading data '
                            'for channels \n  %s' % '\n  '.join(
                                '.'.join(x[:4]) for x in selection_now))

                    f.close()
                    i += neach

    if options.dry_run:
        nslcs = sorted(available_through.keys())

        all_channels = defaultdict(set)
        all_stations = defaultdict(set)

        def plural_s(x):
            return '' if x == 1 else 's'

        for nslc in nslcs:
            sites = tuple(sorted(available_through[nslc]))
            logger.info('selected: %s.%s.%s.%s from site%s %s' % (
                nslc + (plural_s(len(sites)), '+'.join(sites))))

            all_channels[sites].add(nslc)
            all_stations[sites].add(nslc[:3])

        nchannels_all = 0
        nstations_all = 0
        for sites in sorted(
                all_channels.keys(),
                key=lambda sites: (-len(sites), sites)):

            nchannels = len(all_channels[sites])
            nstations = len(all_stations[sites])
            nchannels_all += nchannels
            nstations_all += nstations
            logger.info(
                'selected (%s): %i channel%s (%i station%s)' % (
                    '+'.join(sites),
                    nchannels,
                    plural_s(nchannels),
                    nstations,
                    plural_s(nstations)))

        logger.info(
            'selected total: %i channel%s (%i station%s)' % (
                nchannels_all,
                plural_s(nchannels_all),
                nstations_all,
                plural_s(nstations_all)))

        logger.info('dry run done.')
        sys.exit(0)

    for nslc in have_data:
        # if we are in continue mode, we have to guess where the data came from
        if not any(nslc in have_data_site[site] for site in sites):
            for site in sites:
                if nslc in could_have_data_site[site]:
                    have_data_site[site].add(nslc)

    sxs = {}
    for site in sites:
        selection = []
        for nslc in sorted(have_data_site[site]):
            selection.append(nslc + (tmin-tpad, tmax+tpad))

        if selection:
            logger.info('downloading response information (%s)' % site)
            sxs[site] = fdsn.station(
                site=site, level='response', selection=selection)
            sited = site

            if site == "http://192.168.11.220:8080":
               sited= "bgr_internal"
            elif site == "http://ws.gpi.kit.edu":
               sited= "kit"
            sxs[site].dump_xml(
                filename=op.join(output_dir, 'stations.%s.xml' % sited))

    # chapter 1.5: inject local data

    if options.local_data:
        have_data_site['local'] = set()
        plocal = pile.make_pile(options.local_data, fileformat='detect')
        for traces in plocal.chopper_grouped(
                gather=lambda tr: tr.nslc_id,
                tmin=tmin,
                tmax=tmax,
                tinc=tinc):

            for tr in traces:
                if tr.nslc_id not in have_data:
                    fns.extend(io.save(traces, fn_template_raw))
                    have_data_site['local'].add(tr.nslc_id)
                    have_data.add(tr.nslc_id)

        sites.append('local')

    if options.local_responses_pz:
        sxs['local'] = epz.make_stationxml(
            epz.iload(options.local_responses_pz))

    if options.local_responses_resp:
        local_stations = []
        for fn in options.local_stations:
            local_stations.extend(
                model.load_stations(fn))

        sxs['local'] = resp.make_stationxml(
            local_stations, resp.iload(options.local_responses_resp))

    if options.local_responses_stationxml:
        sxs['local'] = stationxml.load_xml(
            filename=options.local_responses_stationxml)

    # chapter 1.6: dump raw data stations file

    nsl_to_station = {}
    for site in sites:
        if site in sxs:
            stations = sxs[site].get_pyrocko_stations(timespan=(tmin, tmax))
            for s in stations:
                nsl = s.nsl()
                if nsl not in nsl_to_station:
                    nsl_to_station[nsl] = s

    stations = [
        nsl_to_station[nsl_] for nsl_ in sorted(nsl_to_station.keys())]

    util.ensuredirs(fn_stations_raw)
    model.dump_stations(stations, fn_stations_raw)

    dump_commandline(sys.argv, fn_commandline)

    # chapter 2: restitution

    if not fns:
        logger.error('no data available')
        sys.exit(1)

    p = pile.make_pile(fns, show_progress=False)
    p.get_deltatmin()
    otinc = None
    if otinc is None:
        otinc = nice_seconds_floor(p.get_deltatmin() * 500000.)
    otinc = 3600.
    otmin = math.floor(p.tmin / otinc) * otinc
    otmax = math.ceil(p.tmax / otinc) * otinc
    otpad = tpad*2

    fns = []
    rest_traces_b = []
    win_b = None
    for traces_a in p.chopper_grouped(
            gather=lambda tr: tr.nslc_id,
            tmin=otmin,
            tmax=otmax,
            tinc=otinc,
            tpad=otpad):

        rest_traces_a = []
        win_a = None
        for tr in traces_a:
            win_a = tr.wmin, tr.wmax

            if win_b and win_b[0] >= win_a[0]:
                fns.extend(cut_n_dump(rest_traces_b, win_b, fn_template_rest))
                rest_traces_b = []
                win_b = None

            response = None
            failure = []
            for site in sites:
                try:
                    if site not in sxs:
                        continue
                    response = sxs[site].get_pyrocko_response(
                        tr.nslc_id,
                        timespan=(tr.tmin, tr.tmax),
                        fake_input_units=output_units)

                    break

                except stationxml.NoResponseInformation:
                    failure.append('%s: no response information' % site)

                except stationxml.MultipleResponseInformation:
                    failure.append('%s: multiple response information' % site)

            if response is None:
                failure = ', '.join(failure)

            else:
                failure = ''
                try:
                    rest_tr = tr.transfer(tfade, ftap, response, invert=True)
                    rest_traces_a.append(rest_tr)

                except (trace.TraceTooShort, trace.NoData):
                    failure = 'trace too short'

            if failure:
                logger.warn('failed to restitute trace %s.%s.%s.%s (%s)' %
                            (tr.nslc_id + (failure,)))

        if rest_traces_b:
            rest_traces = trace.degapper(rest_traces_b + rest_traces_a,
deoverlap='crossfade_cos')

            fns.extend(cut_n_dump(rest_traces, win_b, fn_template_rest))
            rest_traces_a = []
            if win_a:
                for tr in rest_traces:
                    try:
                        rest_traces_a.append(
                            tr.chop(win_a[0], win_a[1]+otpad,
                                    inplace=False))
                    except trace.NoData:
                        pass

        rest_traces_b = rest_traces_a
        win_b = win_a

    fns.extend(cut_n_dump(rest_traces_b, win_b, fn_template_rest))

    # chapter 3: rotated restituted traces for inspection

    if not event:
        sys.exit(0)

    fn_template1 = \
        'DISPL.%(network)s.%(station)s.%(location)s.%(channel)s'

    fn_waveforms = op.join(output_dir, 'prepared',  fn_template1)
    fn_stations = op.join(output_dir, 'stations.prepared.txt')
    fn_event = op.join(event_dir, 'event.txt')

    nsl_to_station = {}
    for site in sites:
        if site in sxs:
            stations = sxs[site].get_pyrocko_stations(timespan=(tmin, tmax))
            for s in stations:
                nsl = s.nsl()
                if nsl not in nsl_to_station:
                    nsl_to_station[nsl] = s

    p = pile.make_pile(fns, show_progress=False)

    deltat = None
    if sample_rate is not None:
        deltat = 1.0 / sample_rate

    used_stations = []
    for nsl, s in nsl_to_station.items():
        s.set_event_relative_data(event)
        traces = p.all(trace_selector=lambda tr: tr.nslc_id[:3] == nsl)

        keep = []
        for tr in traces:
            if deltat is not None:
                try:
                    tr.downsample_to(deltat, snap=True, allow_upsample_max=5)
                    keep.append(tr)
                except util.UnavailableDecimation as e:
                    logger.warn('Cannot downsample %s.%s.%s.%s: %s'
                                % (tr.nslc_id + (e,)))
                    continue

        if options.out_components == 'rtu':
            pios = s.guess_projections_to_rtu(out_channels=('R', 'T', 'Z'))
        elif options.out_components == 'enu':
            pios = s.guess_projections_to_enu(out_channels=('E', 'N', 'Z'))
        else:
            assert False

        for (proj, in_channels, out_channels) in pios:

            proc = trace.project(traces, proj, in_channels, out_channels)
            for tr in proc:
                for ch in out_channels:
                    if ch.name == tr.channel:
                        s.add_channel(ch)

            if proc:
                io.save(proc, fn_waveforms)
                used_stations.append(s)

    stations = list(used_stations)
    util.ensuredirs(fn_stations)
    model.dump_stations(stations, fn_stations)
    model.dump_events([event], fn_event)

    logger.info('prepared waveforms from %i stations' % len(stations))
    # for seiscomp format:
    # subprocess.run(['jackseis', '%s/traces_disp_%s_%s.mseed' % (folder, tmin, tmax), '--output=%s' % folder + '/%(wmin_year)s/%(network)s/%(station)s/%(channel)s.D/%(network)s.%(station)s..%(channel)s.D.%(wmin)s' , '--force'])

if __name__ == '__main__':
    main()