StreamConverter.py

"""Convert streams from one coordinate system to another.

We use three coordinate systems.
Geo: Based on Geographic North.  X, Y, Z, F
     X is north, Y is east, Z is down
Obs: Based on the observatories orientaion. H, E, Z, F, d0
Mag: Based on Magnetic North. H, D, Z, F
"""

import numpy
import obspy.core
import ChannelConverter


def get_geo_from_mag(mag):
    """Convert a stream to geographic coordinate system.

    Parameters
    ----------
    stream : obspy.core.Stream
        stream containing observatory components H, D, Z, and F.

    Returns
    -------
    obspy.core.Stream
        new stream object containing geographic components X, Y, Z, and F.
    """
    h = mag.select(channel='H')[0]
    d = mag.select(channel='D')[0]
    z = mag.select(channel='Z')[0]
    f = mag.select(channel='F')[0]
    mag_h = h.data
    mag_d = d.data
    (geo_x, geo_y) = ChannelConverter.get_geo_from_mag(mag_h, mag_d)
    return obspy.core.Stream((
        __get_trace('X', h.stats, geo_x),
        __get_trace('Y', d.stats, geo_y),
        z, f))


def get_geo_from_obs(obs):
    """Convert a stream to geographic coordinate system.

    Parameters
    ----------
    stream : obspy.core.Stream
        stream containing observatory components H, D or E, Z, and F.

    Returns
    -------
    obspy.core.Stream
        new stream object containing geographic components X, Y, Z, and F.
    """
    return get_geo_from_mag(get_mag_from_obs(obs))


def get_deltaf_from_geo(geo):
    """Get deltaf given geographic coordinate values

    Parameters
    ----------
    obs: obspy.core.Stream
        stream containing the observatory components H, D or E, Z, and F.

    Returns
    -------
    obspy.core.Stream
        stream object containing delta f values
    """
    x = geo.select(channel='X')[0]
    y = geo.select(channel='Y')[0]
    z = geo.select(channel='Z')[0]
    fs = geo.select(channel='F')[0]
    fv = ChannelConverter.get_computed_f_using_squares(x, y, z)
    G = ChannelConverter.get_deltaf(fv, fs)
    return obspy.core.Stream((
            __get_trace('G', x.stats, G), ))


def get_deltaf_from_obs(obs):
    """Get deltaf given observatory coordinate values

    Parameters
    ----------
    obs: obspy.core.Stream
        stream containing the observatory components H, D or E, Z, and F.

    Returns
    -------
    obspy.core.Stream
        stream object containing delta f values
    """
    h = obs.select(channel='H')[0]
    z = obs.select(channel='Z')[0]
    fs = obs.select(channel='F')[0]
    e = __get_obs_e_from_obs(obs)
    fv = ChannelConverter.get_computed_f_using_squares(h, e, z)
    G = ChannelConverter.get_deltaf(fv, fs)
    return obspy.core.Stream((
            __get_trace('G', h.stats, G), ))


def get_mag_from_geo(geo):
    """Convert a stream to magnetic coordinate system.

    Parameters
    ----------
    geo: obspy.core.Stream
        stream containing observatory components X, Y, Z, and F.

    Returns
    -------
    obspy.core.Stream
        new stream object containing magnetic components H, D, Z, and F.
    """
    x = geo.select(channel='X')[0]
    y = geo.select(channel='Y')[0]
    z = geo.select(channel='Z')[0]
    f = geo.select(channel='F')[0]
    geo_x = x.data
    geo_y = y.data
    (mag_h, mag_d) = ChannelConverter.get_mag_from_geo(geo_x, geo_y)
    return obspy.core.Stream((
            __get_trace('H', x.stats, mag_h),
            __get_trace('D', y.stats, mag_d),
            z, f))


def get_mag_from_obs(obs):
    """Convert a stream to magnetic coordinate system.

    Parameters
    ----------
    obs : obspy.core.Stream
        stream containing observatory components H, D or E, Z, and F.

    Returns
    -------
    obspy.core.Stream
        new stream object containing magnetic components H, D, Z, and F.
    """
    h = obs.select(channel='H')[0]
    e = __get_obs_e_from_obs(obs)
    z = obs.select(channel='Z')[0]
    f = obs.select(channel='F')[0]
    obs_h = h.data
    obs_e = e.data
    d0 = ChannelConverter.get_radians_from_minutes(
            numpy.float64(e.stats.declination_base) / 10)
    (mag_h, mag_d) = ChannelConverter.get_mag_from_obs(obs_h, obs_e, d0)
    return obspy.core.Stream((
            __get_trace('H', h.stats, mag_h),
            __get_trace('D', e.stats, mag_d),
            z, f))


def get_obs_from_geo(geo, include_d=False):
    """Convert a stream to observatory coordinate system.

    Parameters
    ----------
    stream : obspy.core.Stream
        stream containing geographic components X, Y, Z, and F.
    include_d : boolean
        whether to also include the observatory D component.

    Returns
    -------
    obspy.core.Stream
        new stream object containing observatory components H, D, E, Z, and F.
    """
    return get_obs_from_mag(get_mag_from_geo(geo), include_d)


def get_obs_from_mag(mag, include_d=False):
    """Convert a stream to magnetic observatory coordinate system.

    Parameters
    ----------
    stream: obspy.core.Stream
        stream containing magnetic components H, D, Z, and F.
    include_d: boolean
        whether to also include the observatory D component
    Returns
    -------
    obspy.core.Stream
        new stream object containing observatory components H, D, E, Z, and F
    """
    h = mag.select(channel='H')[0]
    d = mag.select(channel='D')[0]
    z = mag.select(channel='Z')[0]
    f = mag.select(channel='F')[0]
    mag_h = h.data
    mag_d = d.data
    d0 = ChannelConverter.get_radians_from_minutes(
        numpy.float64(d.stats.declination_base) / 10)
    (obs_h, obs_e) = ChannelConverter.get_obs_from_mag(mag_h, mag_d, d0)

    traces = (
        __get_trace('H', h.stats, obs_h),
        __get_trace('E', d.stats, obs_e),
        z, f)
    if include_d:
        obs_d = ChannelConverter.get_obs_d_from_obs(obs_h, obs_e)
        traces = traces + (__get_trace('D', d.stats, obs_d),)
    return obspy.core.Stream(traces)


def get_obs_from_obs(obs, include_e=False, include_d=False):
    """Fill in the observatory parameters as requested

    Parameters
    ----------
    stream: obspy.core.Stream
        stream containing the observatory components H, D or E, Z, and F.
    include_e: boolean
        whether to include the e component
    include_d: boolean
        whether to include the d component

    Returns
    -------
    obspy.core.Stream
        new stream object containing observatory components H, D, E, Z, and F
    """
    h = obs.select(channel='H')[0]
    z = obs.select(channel='Z')[0]
    f = obs.select(channel='F')[0]
    traces = (h, z, f)
    if include_d:
        d = __get_obs_d_from_obs(obs)
        traces = traces + (d, )
    if include_e:
        e = __get_obs_e_from_obs(obs)
        traces = traces + (e, )
    return obspy.core.Stream(traces)


def __get_trace(channel, stats, data):
    """Utility to create a new trace object.

    Parameters
    ----------
    channel : str
        channel name.
    stats : obspy.core.Stats
        channel metadata to clone.
    data : numpy.array
        channel data.

    Returns
    -------
    obspy.core.Trace
        trace containing data and metadata.
    """
    stats = obspy.core.Stats(stats)
    stats.channel = channel
    return obspy.core.Trace(data, stats)


def __get_obs_d_from_obs(obs):
    """Get trace containing observatory D component.

    Returns D if found, otherwise computes D from H, E, D0.

    Parameters
    ----------
    obs : obspy.core.Stream
        observatory components (D) or (H, E).

    Returns
    -------
    obspy.core.Trace
        observatory component D.
    """
    try:
        d = obs.select(channel='D')[0]
    except:
        h = obs.select(channel='H')[0]
        e = obs.select(channel='E')[0]
        d = __get_trace('D', e.stats,
                ChannelConverter.get_obs_d_from_obs(h.data, e.data))
    return d


def __get_obs_e_from_obs(obs):
    """Get trace containing observatory E component.

    Returns E if found, otherwise computes E from H,D.

    Parameters
    ----------
    obs : obspy.core.Stream
        observatory components (E) or (H, D).

    Returns
    -------
    obspy.core.Trace
        observatory component E.
    """
    try:
        e = obs.select(channel='E')[0]
    except:
        h = obs.select(channel='H')[0]
        d = obs.select(channel='D')[0]
        e = __get_trace('E', d.stats,
                ChannelConverter.get_obs_e_from_obs(h.data, d.data))
    return e