Merge pull request #35 from hasimpson-usgs/deltaf

Deltaf

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()
 

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

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()
 

geomagio/DeltaFAlgorithm.py

+"""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

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

geomagio/__init__.py

@@ -18,6 +18,7 @@ __all__ = [
     'AlgorithmException',
     'ChannelConverter',
     'Controller',
+    'DeltaFAlgorithm',
     'StreamConverter',
     'TimeseriesFactory',
     'TimeseriesFactoryException',

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)