def amp_pick_event(event, st, respdir, chans=['Z'], var_wintype=True,
                   winlen=0.9, pre_pick=0.2, pre_filt=True, lowcut=1.0,
                   highcut=20.0, corners=4, min_snr=1.0, plot=False,
                   remove_old=False, ps_multiplier=0.34, velocity=False):
    """
    Pick amplitudes for local magnitude for a single event.

    Looks for maximum peak-to-trough amplitude for a channel in a stream, and
    picks this amplitude and period.  There are a few things it does
    internally to stabilise the result:

        1. Applies a given filter to the data - very necessary for small
        magnitude earthquakes;

        2. Keeps track of the poles and zeros of this filter and removes them
        from the picked amplitude;

        3. Picks the peak-to-trough amplitude, but records half of this: the
        specification for the local magnitude is to use a peak amplitude on
        a horizontal, however, with modern digital seismometers, the peak
        amplitude often has an additional, DC-shift applied to it, to
        stabilise this, and to remove possible issues with de-meaning data
        recorded during the wave-train of an event (e.g. the mean may not be
        the same as it would be for longer durations), we use half the
        peak-to-trough amplitude;

        4. Despite the original definition of local magnitude requiring the
        use of a horizontal channel, more recent work has shown that the
        vertical channels give more consistent magnitude estimations between
        stations, due to a reduction in site-amplification effects, we
        therefore use the vertical channels by default, but allow the user
        to chose which channels they deem appropriate;

        5. We do not specify that the maximum amplitude should be the
        S-phase: The original definition holds that the maximum body-wave
        amplitude should be used - while this is often the S-phase, we do not
        discriminate against the P-phase.  We do note that, unless the user
        takes care when assigning winlen and filters, they may end up with
        amplitude picks for surface waves;

        6. We use a variable window-length by default that takes into account
        P-S times if available, this is in an effort to include only the
        body waves.  When P-S times are not available we us the ps_multiplier
        variable, which defaults to 0.34 x hypocentral distance.

    :type event: obspy.core.event.event.Event
    :param event: Event to pick
    :type st: obspy.core.stream.Stream
    :param st: Stream associated with event
    :type respdir: str
    :param respdir: Path to the response information directory
    :type chans: list
    :param chans:
        List of the channels to pick on, defaults to ['Z'] - should just be
        the orientations, e.g. Z, 1, 2, N, E
    :type var_wintype: bool
    :param var_wintype:
        If True, the winlen will be multiplied by the P-S time if both P and
        S picks are available, otherwise it will be multiplied by the
        hypocentral distance*ps_multiplier, defaults to True
    :type winlen: float
    :param winlen:
        Length of window, see above parameter, if var_wintype is False then
        this will be in seconds, otherwise it is the multiplier to the
        p-s time, defaults to 0.9.
    :type pre_pick: float
    :param pre_pick:
        Time before the s-pick to start the cut window, defaults to 0.2.
    :type pre_filt: bool
    :param pre_filt: To apply a pre-filter or not, defaults to True
    :type lowcut: float
    :param lowcut: Lowcut in Hz for the pre-filter, defaults to 1.0
    :type highcut: float
    :param highcut: Highcut in Hz for the pre-filter, defaults to 20.0
    :type corners: int
    :param corners: Number of corners to use in the pre-filter
    :type min_snr: float
    :param min_snr:
        Minimum signal-to-noise ratio to allow a pick - see note below on
        signal-to-noise ratio calculation.
    :type plot: bool
    :param plot: Turn plotting on or off.
    :type remove_old: bool
    :param remove_old:
        If True, will remove old amplitude picks from event and overwrite
        with new picks. Defaults to False.
    :type ps_multiplier: float
    :param ps_multiplier:
        A p-s time multiplier of hypocentral distance - defaults to 0.34,
        based on p-s ratio of 1.68 and an S-velocity 0f 1.5km/s, deliberately
        chosen to be quite slow.
    :type velocity: bool
    :param velocity:
        Whether to make the pick in velocity space or not. Original definition
        of local magnitude used displacement of Wood-Anderson, MLv in seiscomp
        and Antelope uses a velocity measurement.

    :returns: Picked event
    :rtype: :class:`obspy.core.event.Event`

    .. Note::
        Signal-to-noise ratio is calculated using the filtered data by
        dividing the maximum amplitude in the signal window (pick window)
        by the normalized noise amplitude (taken from the whole window
        supplied).

    .. Warning::
        Works in place on data - will filter and remove response from data,
        you are recommended to give this function a copy of the data if you
        are using it in a loop.
    """
    # Convert these picks into a lists
    stations = []  # List of stations
    channels = []  # List of channels
    picktimes = []  # List of pick times
    picktypes = []  # List of pick types
    picks_out = []
    try:
        depth = _get_origin(event).depth
    except MatchFilterError:
        depth = 0
    if remove_old and event.amplitudes:
        for amp in event.amplitudes:
            # Find the pick and remove it too
            pick = [p for p in event.picks if p.resource_id == amp.pick_id][0]
            event.picks.remove(pick)
        event.amplitudes = []
    for pick in event.picks:
        if pick.phase_hint in ['P', 'S']:
            picks_out.append(pick)  # Need to be able to remove this if there
            # isn't data for a station!
            stations.append(pick.waveform_id.station_code)
            channels.append(pick.waveform_id.channel_code)
            picktimes.append(pick.time)
            picktypes.append(pick.phase_hint)
    if len(picktypes) == 0:
        warnings.warn('No P or S picks found')
    st.merge()  # merge the data, just in case!
    # For each station cut the window
    uniq_stas = list(set(stations))
    for sta in uniq_stas:
        for chan in chans:
            print('Working on ' + sta + ' ' + chan)
            tr = st.select(station=sta, channel='*' + chan)
            if not tr:
                warnings.warn(
                    'There is no station and channel match in the wavefile!')
                continue
            else:
                tr = tr[0]
            # Apply the pre-filter
            if pre_filt:
                try:
                    tr.split().detrend('simple').merge(fill_value=0)
                except:
                    print('Some issue splitting this one')
                    dummy = tr.split()
                    dummy.detrend('simple')
                    tr = dummy.merge(fill_value=0)
                try:
                    tr.filter('bandpass', freqmin=lowcut, freqmax=highcut,
                              corners=corners)
                except NotImplementedError:
                    print('For some reason trace is not continuous:')
                    print(tr)
                    continue
            # Find the response information
            resp_info = _find_resp(
                tr.stats.station, tr.stats.channel, tr.stats.network,
                tr.stats.starttime, tr.stats.delta, respdir)
            PAZ = []
            seedresp = []
            if resp_info and 'gain' in resp_info:
                PAZ = resp_info
            elif resp_info:
                seedresp = resp_info
            # Simulate a Wood Anderson Seismograph
            if PAZ and len(tr.data) > 10:
                # Set ten data points to be the minimum to pass
                tr = _sim_WA(tr, PAZ, None, 10, velocity=velocity)
            elif seedresp and len(tr.data) > 10:
                tr = _sim_WA(tr, None, seedresp, 10, velocity=velocity)
            elif len(tr.data) > 10:
                warnings.warn('No PAZ for ' + tr.stats.station + ' ' +
                              tr.stats.channel + ' at time: ' +
                              str(tr.stats.starttime))
                continue
            sta_picks = [i for i in range(len(stations))
                         if stations[i] == sta]
            pick_id = event.picks[sta_picks[0]].resource_id
            arrival = [arrival for arrival in event.origins[0].arrivals
                       if arrival.pick_id == pick_id][0]
            hypo_dist = np.sqrt(
                np.square(degrees2kilometers(arrival.distance)) +
                np.square(depth / 1000))
            if var_wintype and hypo_dist:
                if 'S' in [picktypes[i] for i in sta_picks] and\
                   'P' in [picktypes[i] for i in sta_picks]:
                    # If there is an S-pick we can use this :D
                    s_pick = [picktimes[i] for i in sta_picks
                              if picktypes[i] == 'S']
                    s_pick = min(s_pick)
                    p_pick = [picktimes[i] for i in sta_picks
                              if picktypes[i] == 'P']
                    p_pick = min(p_pick)
                    try:
                        tr.trim(starttime=s_pick - pre_pick,
                                endtime=s_pick + (s_pick - p_pick) * winlen)
                    except ValueError:
                        continue
                elif 'S' in [picktypes[i] for i in sta_picks]:
                    s_pick = [picktimes[i] for i in sta_picks
                              if picktypes[i] == 'S']
                    s_pick = min(s_pick)
                    p_modelled = s_pick - (hypo_dist * ps_multiplier)
                    try:
                        tr.trim(starttime=s_pick - pre_pick,
                                endtime=s_pick + (s_pick - p_modelled) *
                                winlen)
                    except ValueError:
                        continue
                else:
                    # In this case we only have a P pick
                    p_pick = [picktimes[i] for i in sta_picks
                              if picktypes[i] == 'P']
                    p_pick = min(p_pick)
                    s_modelled = p_pick + (hypo_dist * ps_multiplier)
                    print('P_pick=%s' % str(p_pick))
                    print('hypo_dist: %s' % str(hypo_dist))
                    print('S modelled=%s' % str(s_modelled))
                    try:
                        tr.trim(starttime=s_modelled - pre_pick,
                                endtime=s_modelled + (s_modelled - p_pick) *
                                winlen)
                        print(tr)
                    except ValueError:
                        continue
                # Work out the window length based on p-s time or distance
            elif 'S' in [picktypes[i] for i in sta_picks]:
                # If the window is fixed we still need to find the start time,
                # which can be based either on the S-pick (this elif), or
                # on the hypocentral distance and the P-pick

                # Take the minimum S-pick time if more than one S-pick is
                # available
                s_pick = [picktimes[i] for i in sta_picks
                          if picktypes[i] == 'S']
                s_pick = min(s_pick)
                try:
                    tr.trim(starttime=s_pick - pre_pick,
                            endtime=s_pick + winlen)
                except ValueError:
                    continue
            else:
                # In this case, there is no S-pick and the window length is
                # fixed we need to calculate an expected S_pick based on the
                # hypocentral distance, this will be quite hand-wavey as we
                # are not using any kind of velocity model.
                p_pick = [picktimes[i] for i in sta_picks
                          if picktypes[i] == 'P']
                print(picktimes)
                p_pick = min(p_pick)
                s_modelled = p_pick + hypo_dist * ps_multiplier
                try:
                    tr.trim(starttime=s_modelled - pre_pick,
                            endtime=s_modelled + winlen)
                except ValueError:
                    continue
            if len(tr.data) <= 10:
                warnings.warn('No data found for: ' + tr.stats.station)
                continue
            # Get the amplitude
            try:
                amplitude, period, delay = _max_p2t(tr.data, tr.stats.delta)
            except ValueError:
                print('No amplitude picked for tr %s' % str(tr))
                continue
            # Calculate the normalized noise amplitude
            noise_amplitude = np.sqrt(np.mean(np.square(tr.data)))
            if amplitude == 0.0:
                continue
            if amplitude / noise_amplitude < min_snr:
                print('Signal to noise ratio of %s is below threshold.' %
                      (amplitude / noise_amplitude))
                continue
            if plot:
                plt.plot(np.arange(len(tr.data)), tr.data, 'k')
                plt.scatter(tr.stats.sampling_rate * delay, amplitude / 2)
                plt.scatter(tr.stats.sampling_rate * (delay + period),
                            -amplitude / 2)
                plt.show()
            print('Amplitude picked: ' + str(amplitude))
            print('Signal-to-noise ratio is: %s' %
                  (amplitude / noise_amplitude))
            # Note, amplitude should be in meters at the moment!
            # Remove the pre-filter response
            if pre_filt:
                # Generate poles and zeros for the filter we used earlier: this
                # is how the filter is designed in the convenience methods of
                # filtering in obspy.
                z, p, k = iirfilter(
                    corners, [lowcut / (0.5 * tr.stats.sampling_rate),
                              highcut / (0.5 * tr.stats.sampling_rate)],
                    btype='band', ftype='butter', output='zpk')
                filt_paz = {'poles': list(p), 'zeros': list(z), 'gain': k,
                            'sensitivity': 1.0}
                amplitude /= (paz_2_amplitude_value_of_freq_resp(
                    filt_paz, 1 / period) * filt_paz['sensitivity'])
            if PAZ:
                amplitude /= 1000
            if seedresp:  # Seedresp method returns mm
                amplitude *= 1000000
            # Write out the half amplitude, approximately the peak amplitude as
            # used directly in magnitude calculations
            amplitude *= 0.5
            # Append an amplitude reading to the event
            _waveform_id = WaveformStreamID(
                station_code=tr.stats.station, channel_code=tr.stats.channel,
                network_code=tr.stats.network)
            pick_ind = len(event.picks)
            event.picks.append(Pick(
                waveform_id=_waveform_id, phase_hint='IAML',
                polarity='undecidable', time=tr.stats.starttime + delay,
                evaluation_mode='automatic'))
            if not velocity:
                event.amplitudes.append(Amplitude(
                    generic_amplitude=amplitude / 1e9, period=period,
                    pick_id=event.picks[pick_ind].resource_id,
                    waveform_id=event.picks[pick_ind].waveform_id, unit='m',
                    magnitude_hint='ML', type='AML', category='point'))
            else:
                event.amplitudes.append(Amplitude(
                    generic_amplitude=amplitude / 1e9, period=period,
                    pick_id=event.picks[pick_ind].resource_id,
                    waveform_id=event.picks[pick_ind].waveform_id, unit='m/s',
                    magnitude_hint='ML', type='AML', category='point'))
    return event
Pick amplitudes for local magnitude for a single event.

    Looks for maximum peak-to-trough amplitude for a channel in a stream, and
    picks this amplitude and period.  There are a few things it does
    internally to stabilise the result:

        1. Applies a given filter to the data - very necessary for small
        magnitude earthquakes;

        2. Keeps track of the poles and zeros of this filter and removes them
        from the picked amplitude;

        3. Picks the peak-to-trough amplitude, but records half of this: the
        specification for the local magnitude is to use a peak amplitude on
        a horizontal, however, with modern digital seismometers, the peak
        amplitude often has an additional, DC-shift applied to it, to
        stabilise this, and to remove possible issues with de-meaning data
        recorded during the wave-train of an event (e.g. the mean may not be
        the same as it would be for longer durations), we use half the
        peak-to-trough amplitude;

        4. Despite the original definition of local magnitude requiring the
        use of a horizontal channel, more recent work has shown that the
        vertical channels give more consistent magnitude estimations between
        stations, due to a reduction in site-amplification effects, we
        therefore use the vertical channels by default, but allow the user
        to chose which channels they deem appropriate;

        5. We do not specify that the maximum amplitude should be the
        S-phase: The original definition holds that the maximum body-wave
        amplitude should be used - while this is often the S-phase, we do not
        discriminate against the P-phase.  We do note that, unless the user
        takes care when assigning winlen and filters, they may end up with
        amplitude picks for surface waves;

        6. We use a variable window-length by default that takes into account
        P-S times if available, this is in an effort to include only the
        body waves.  When P-S times are not available we us the ps_multiplier
        variable, which defaults to 0.34 x hypocentral distance.

    :type event: obspy.core.event.event.Event
    :param event: Event to pick
    :type st: obspy.core.stream.Stream
    :param st: Stream associated with event
    :type respdir: str
    :param respdir: Path to the response information directory
    :type chans: list
    :param chans:
        List of the channels to pick on, defaults to ['Z'] - should just be
        the orientations, e.g. Z, 1, 2, N, E
    :type var_wintype: bool
    :param var_wintype:
        If True, the winlen will be multiplied by the P-S time if both P and
        S picks are available, otherwise it will be multiplied by the
        hypocentral distance*ps_multiplier, defaults to True
    :type winlen: float
    :param winlen:
        Length of window, see above parameter, if var_wintype is False then
        this will be in seconds, otherwise it is the multiplier to the
        p-s time, defaults to 0.9.
    :type pre_pick: float
    :param pre_pick:
        Time before the s-pick to start the cut window, defaults to 0.2.
    :type pre_filt: bool
    :param pre_filt: To apply a pre-filter or not, defaults to True
    :type lowcut: float
    :param lowcut: Lowcut in Hz for the pre-filter, defaults to 1.0
    :type highcut: float
    :param highcut: Highcut in Hz for the pre-filter, defaults to 20.0
    :type corners: int
    :param corners: Number of corners to use in the pre-filter
    :type min_snr: float
    :param min_snr:
        Minimum signal-to-noise ratio to allow a pick - see note below on
        signal-to-noise ratio calculation.
    :type plot: bool
    :param plot: Turn plotting on or off.
    :type remove_old: bool
    :param remove_old:
        If True, will remove old amplitude picks from event and overwrite
        with new picks. Defaults to False.
    :type ps_multiplier: float
    :param ps_multiplier:
        A p-s time multiplier of hypocentral distance - defaults to 0.34,
        based on p-s ratio of 1.68 and an S-velocity 0f 1.5km/s, deliberately
        chosen to be quite slow.
    :type velocity: bool
    :param velocity:
        Whether to make the pick in velocity space or not. Original definition
        of local magnitude used displacement of Wood-Anderson, MLv in seiscomp
        and Antelope uses a velocity measurement.

    :returns: Picked event
    :rtype: :class:`obspy.core.event.Event`

    .. Note::
        Signal-to-noise ratio is calculated using the filtered data by
        dividing the maximum amplitude in the signal window (pick window)
        by the normalized noise amplitude (taken from the whole window
        supplied).

    .. Warning::
        Works in place on data - will filter and remove response from data,
        you are recommended to give this function a copy of the data if you
        are using it in a loop.