diff --git a/bin/geomag.py b/bin/geomag.py
index 67c34d75b4dc5cf1c77f95659b8784b481f3e4b1..46a1689d3015b01aa87c282049fc3d85f8cbaf47 100755
--- a/bin/geomag.py
+++ b/bin/geomag.py
@@ -17,6 +17,7 @@ import geomagio.pcdcp as pcdcp
from geomagio import Algorithm, \
Controller
from geomagio.XYZAlgorithm import XYZAlgorithm
+from geomagio.DeltaFAlgorithm import DeltaFAlgorithm
from obspy.core import UTCDateTime
@@ -142,6 +143,8 @@ def main():
if args.xyz is not None:
algorithm = XYZAlgorithm(informat=args.xyz[0],
outformat=args.xyz[1])
+ elif args.deltaf is not None:
+ algorithm = DeltaFAlgorithm(informat=args.deltaf)
else:
# TODO get smarter on inchannels/outchannels since input doesn't always
# need to use the --inchannels argument, but might (as in iaga2002),
@@ -277,6 +280,9 @@ def parse_args():
choices=['geo', 'mag', 'obs', 'obsd'],
help='Enter the geomagnetic orientation(s) you want to read from' +
' and to respectfully.')
+ algorithm_group.add_argument('--deltaf',
+ choices=['geo', 'obs', 'obsd'],
+ help='Enter the geomagnetic orientation you want to read from')
return parser.parse_args()
diff --git a/geomagio/ChannelConverter.py b/geomagio/ChannelConverter.py
index 8181e42cd2c96f78d8eda7e6703450f777ceb398..1058f6c9c059f3a66361277fa01faa68c23220ea 100644
--- a/geomagio/ChannelConverter.py
+++ b/geomagio/ChannelConverter.py
@@ -378,6 +378,44 @@ def get_obs_h_from_mag(h, d, d0=0):
return numpy.multiply(h, numpy.cos(obs_d))
+def get_deltaf(fv, fs):
+ """gets the deltaf value given the scalar F and computed F values.
+
+ Parameters
+ ----------
+ fv: array_like
+ the F vector computed from the observatory instruments
+ fs: array_like
+ the measured F value
+ """
+ return numpy.subtract(fv, fs)
+
+
+def get_computed_f_using_squares(x, y, z):
+ """gets the computed f value
+
+ Parameters
+ ----------
+ x: array_like
+ the x component from the observatory
+ y: array_like
+ the y component from the observatory
+ z: array_like
+ the z component from the observatory
+ Notes
+ -----
+ This works for geographic coordinates, or observatory coordinates.
+ ie x, y, z or h, e, z
+ We're using variables x,y,z to represent generic cartisian coordinates.
+ """
+ x2 = numpy.square(x)
+ y2 = numpy.square(y)
+ z2 = numpy.square(z)
+ fv = numpy.add(x2, y2)
+ fv = numpy.add(fv, z2)
+ return numpy.sqrt(fv)
+
+
def get_radians_from_minutes(m):
"""gets the radian value given the decimal value
Parameters
diff --git a/geomagio/ChannelConverter_test.py b/geomagio/ChannelConverter_test.py
index 00fe4bb6e5f98beae82402d369ff52b2f0a1277a..0c33b4b04b6f5cec4ad957a58eb36b412f204694 100644
--- a/geomagio/ChannelConverter_test.py
+++ b/geomagio/ChannelConverter_test.py
@@ -45,6 +45,21 @@ class ChannelConverterTest:
length to be cos(angle) and the opposite length to be sin(angle)
"""
+ def test_get_computed_f_using_sqaures(self):
+ """geomag.ChannelConverterTest.test_get_computed_f_using_sqaures
+
+ Computed f is the combination of the 3 vector components from
+ either the geographic coordinate system (X,Y,Z), or the observatory
+ coordinate system (h,e,z).
+ """
+ h = 2
+ e = 2
+ z = 2
+ fv = channel.get_computed_f_using_squares(h, e, z)
+ fs = math.sqrt(12)
+ assert_almost_equal(fv, fs, 8, 'Expect fv to almost equal sqrt(12)',
+ True)
+
def test_get_geo_from_obs(self):
"""geomag.ChannelConverterTest.test_get_geo_from_obs()
diff --git a/geomagio/DeltaFAlgorithm.py b/geomagio/DeltaFAlgorithm.py
new file mode 100644
index 0000000000000000000000000000000000000000..d0619dee63ae90534ec4034063c383219cc90947
--- /dev/null
+++ b/geomagio/DeltaFAlgorithm.py
@@ -0,0 +1,70 @@
+"""Algorithm that creates a deltaf
+
+"""
+
+from Algorithm import Algorithm
+from AlgorithmException import AlgorithmException
+import StreamConverter as StreamConverter
+
+# List of channels by geomagnetic observatory orientation.
+# geo represents a geographic north/south orientation
+# obs represents the sensor orientation aligned close to the mag orientation
+# obsd is the same as obs, but with D(declination) instead of E (e/w vector)
+CHANNELS = {
+ 'geo': ['X', 'Y', 'Z', 'F'],
+ 'obs': ['H', 'E', 'Z', 'F'],
+ 'obsd': ['H', 'D', 'Z', 'F']
+}
+
+
+class DeltaFAlgorithm(Algorithm):
+ """Algorithm for getting Delta F.
+
+ Parameters
+ ----------
+ informat: str
+ the code that represents the incoming data form that the Algorithm
+ will be converting from.
+ """
+
+ def __init__(self, informat=None):
+ Algorithm.__init__(self, inchannels=CHANNELS[informat],
+ outchannels=['G'])
+ self.informat = informat
+
+ def check_stream(self, timeseries):
+ """checks a stream to make certain all the required channels
+ exist.
+
+ Parameters
+ ----------
+ timeseries: obspy.core.Stream
+ stream to be checked.
+ """
+ for channel in self._inchannels:
+ if len(timeseries.select(channel=channel)) == 0:
+ raise AlgorithmException(
+ 'Channel %s not found in input' % channel)
+
+ def process(self, timeseries):
+ """converts a timeseries stream into a different coordinate system
+
+ Parameters
+ ----------
+ informat: string
+ indicates the input coordinate system.
+
+ Returns
+ -------
+ out_stream: obspy.core.Stream
+ new stream object containing the converted coordinates.
+ """
+ self.check_stream(timeseries)
+ out_stream = None
+
+ if self.informat == 'geo':
+ out_stream = StreamConverter.get_deltaf_from_geo(timeseries)
+ elif self.informat == 'obs' or self.informat == 'obsd':
+ out_stream = StreamConverter.get_deltaf_from_obs(timeseries)
+
+ return out_stream
diff --git a/geomagio/StreamConverter.py b/geomagio/StreamConverter.py
index 199ad55bbb7d096b7ec53477838fc581e6563ba5..e352ebb53501564078e693e16cfd1f4e8eb40b88 100644
--- a/geomagio/StreamConverter.py
+++ b/geomagio/StreamConverter.py
@@ -54,6 +54,52 @@ def get_geo_from_obs(obs):
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.
@@ -166,7 +212,7 @@ def get_obs_from_obs(obs, include_e=False, include_d=False):
Parameters
----------
stream: obspy.core.Stream
- stream containing the observatory components H, D of E, Z, and F.
+ stream containing the observatory components H, D or E, Z, and F.
include_e: boolean
whether to include the e component
include_d: boolean
diff --git a/geomagio/__init__.py b/geomagio/__init__.py
index 551b00c5272c3087a932b4a42e31a59530b26485..74c7d9ea771a19bcb02cd7068c234e885425d0f8 100644
--- a/geomagio/__init__.py
+++ b/geomagio/__init__.py
@@ -18,6 +18,7 @@ __all__ = [
'AlgorithmException',
'ChannelConverter',
'Controller',
+ 'DeltaFAlgorithm',
'StreamConverter',
'TimeseriesFactory',
'TimeseriesFactoryException',
diff --git a/geomagio/edge/EdgeFactory.py b/geomagio/edge/EdgeFactory.py
index 795eb2ebfda0ce4f5c73a059f05c760ef1c693b1..06305985d14d0a804695501bfcd56c252ce66a84 100644
--- a/geomagio/edge/EdgeFactory.py
+++ b/geomagio/edge/EdgeFactory.py
@@ -342,6 +342,8 @@ class EdgeFactory(TimeseriesFactory):
edge_channel = edge_interval_code + 'VH'
elif channel == 'Z':
edge_channel = edge_interval_code + 'VZ'
+ elif channel == 'G':
+ edge_channel = edge_interval_code + 'SG'
else:
raise TimeseriesFactoryException(
'Unexpected channel code "%s"' % channel)