From b47ff27594e9cda5c14d29a5ef639fd96b233cc8 Mon Sep 17 00:00:00 2001
From: Jeremy Fee <jmfee@usgs.gov>
Date: Fri, 19 Dec 2014 15:34:45 -0700
Subject: [PATCH] Updated main to work with updated factory, added algorithm
interface, streamconverter
---
src/python/geomag/Algorithm.py | 20 ++++
src/python/geomag/StreamConverter.py | 162 +++++++++++++++++++++++++++
src/python/main.py | 7 +-
3 files changed, 186 insertions(+), 3 deletions(-)
create mode 100644 src/python/geomag/Algorithm.py
create mode 100644 src/python/geomag/StreamConverter.py
diff --git a/src/python/geomag/Algorithm.py b/src/python/geomag/Algorithm.py
new file mode 100644
index 00000000..1c05eaae
--- /dev/null
+++ b/src/python/geomag/Algorithm.py
@@ -0,0 +1,20 @@
+
+
+class Algorithm(object):
+ """An algorithm processes a stream of timeseries to produce new timeseries.
+ """
+
+ def process(self, stream):
+ """Process a chunk of data.
+
+ Parameters
+ ----------
+ stream : obspy.core.Stream
+ input data
+
+ Returns
+ -------
+ obspy.core.Stream
+ resulting data
+ """
+ return stream.copy()
diff --git a/src/python/geomag/StreamConverter.py b/src/python/geomag/StreamConverter.py
new file mode 100644
index 00000000..8b4f9cbe
--- /dev/null
+++ b/src/python/geomag/StreamConverter.py
@@ -0,0 +1,162 @@
+"""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 obspy.core
+import ChannelConverter
+
+
+def get_geo_from_mag(mag):
+ pass
+
+
+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_mag_from_geo(geo):
+ pass
+
+
+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]
+ d = __get_obs_d_from_obs(obs)
+ z = obs.select(channel='Z')[0]
+ f = obs.select(channel='F')[0]
+ obs_h = h.data
+ obs_d = d.data
+ d0 = ChannelConverter.get_radians_from_minutes(
+ numpy.float64(d.stats['DECBAS']) / 10)
+ (mag_h, mag_d) = ChannelConverter.get_mag_from_obs(obs_h, obs_d, d0)
+ return obspy.core.Stream((
+ __get_trace('H', h.stats, mag_h),
+ __get_trace('D', d.stats, mag_d),
+ z, f))
+
+
+def get_obs_from_geo(geo, include_e=False):
+ """Convert a stream to geographic coordinate system.
+
+ Parameters
+ ----------
+ stream : obspy.core.Stream
+ stream containing geographic components X, Y, Z, and F.
+ include_e : boolean
+ whether to also include the observatory E 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))
+
+
+def get_obs_from_mag(mag, include_e=False):
+ # TODO
+ obs = obspy.core.Stream()
+ if include_e:
+ obs += __get_obs_e_from_obs(obs)
+ return obs
+
+
+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
diff --git a/src/python/main.py b/src/python/main.py
index abc365af..daacb0cd 100644
--- a/src/python/main.py
+++ b/src/python/main.py
@@ -15,10 +15,11 @@ def main():
import pprint
iaga_dir = path.normpath(path.join(script_dir, '../../etc/iaga2002'))
factory = iaga2002.IAGA2002Factory('file://' + iaga_dir +
- '/%(OBS)s/%(interval)s%(type)s/%(obs)s%(ymd)s%(t)s%(i)s.%(i)s')
- timeseries = factory.get_timeseries('BOU',
- UTCDateTime('2014-11-01'), UTCDateTime('2014-11-02'),
+ '/%(OBS)s/%(interval)s%(type)s/%(obs)s%(ymd)s%(t)s%(i)s.%(i)s',
+ observatory='BOU', channels=('H', 'D', 'Z', 'F'),
interval='minute', type='variation')
+ timeseries = factory.get_timeseries(
+ UTCDateTime('2014-11-01'), UTCDateTime('2014-11-02'))
print timeseries
--
GitLab