diff --git a/docs/algorithms/Average.md b/docs/algorithms/Average.md
new file mode 100644
index 0000000000000000000000000000000000000000..9a1b0fb098fcecb38b997c1f1d26bb4d2ff701ef
--- /dev/null
+++ b/docs/algorithms/Average.md
@@ -0,0 +1,38 @@
+Average Algorithm
+-----------------
+
+Algorithm Theoretical Basis for "Average Algorithm"
+
+## Summary
+Averages multiple observatory data for the same channel.
+Used mainly to average disturbance data in order to find
+the general disturbance in the magnetic field. The algorithm
+takes data from multiple observatories with the same channel
+and returns a stream of the averaged data.
+
+## Background and Motivation
+The averaging function is used to smooth out the plots and
+find a combined disturbance of the magnetic field that
+encompasses the Earth. This can be used to determine the
+overall effect of magnetic storms over a large area. The
+algorithm is also used to average the Solar Quiet response
+to find the daily change in the magnetic field. The algorithm
+can be used to find the average of any channel over multiple
+observation stations.
+
+## Math and Theory
+Multiple data streams are averaged using a numpy function
+(numpy.mean) that takes multiple ndarrays as an argument
+and averages them into one array. A latitude correction
+can be applied based on the different observation locations.
+The correction is really just a weighting value based on a
+0-1 scale in order to put more validity in some stations.
+
+## Practical Considerations
+The averaging function can be called from geomag.py and stating
+the '--algorithm average' option or by calling the average.py
+script which automatically chooses the averaging option to the
+algorithm. Only one channel at a time can be run through the
+algorithm. Any input that geomag.py can handle can be used
+such as the edge server or a file url input this is also true
+for any output.
diff --git a/geomagio/__init__.py b/geomagio/__init__.py
index e896583b7b50762cf9858d9032c35d3c050a25b4..d871a122884b9940dd200d5ac436c9b2414b5d50 100644
--- a/geomagio/__init__.py
+++ b/geomagio/__init__.py
@@ -25,6 +25,5 @@ __all__ = [
'TimeseriesFactoryException',
'TimeseriesUtility',
'Util',
- 'Url',
- 'XYZAlgorithm'
+ 'Url'
]
diff --git a/geomagio/algorithm/AverageAlgorithm.py b/geomagio/algorithm/AverageAlgorithm.py
new file mode 100644
index 0000000000000000000000000000000000000000..f030fd0e357008f1158bde2cb00fe792f132b058
--- /dev/null
+++ b/geomagio/algorithm/AverageAlgorithm.py
@@ -0,0 +1,210 @@
+"""Algorithm that creates an averaged Dst.
+
+"""
+from __future__ import absolute_import
+
+from .Algorithm import Algorithm
+from .AlgorithmException import AlgorithmException
+from ..ObservatoryMetadata import ObservatoryMetadata
+import numpy
+import obspy.core
+
+
+# Possible correction factors.
+# Defaults to 1.0 if station not found in list.
+CORR = {
+ 'HON': 1.0,
+ 'SJG': 1.0,
+ 'HER': 1.0,
+ 'KAK': 1.0,
+ 'GUA': 1.0
+}
+
+
+class AverageAlgorithm(Algorithm):
+ """Algorithm that creates an averaged Dst.
+
+ Parameters
+ ----------
+
+ """
+
+ def __init__(self):
+ Algorithm.__init__(self)
+ self._npts = -1
+ self._stt = -1
+ self._stats = None
+ self.observatoryMetadata = ObservatoryMetadata()
+
+ def check_stream(self, timeseries):
+ """checks a stream to make certain the required data
+ exists.
+
+ Parameters
+ ----------
+ timeseries: obspy.core.Stream
+ stream to be checked.
+ """
+
+ # A stream produced by EdgeFactory should always pass these checks.
+
+ # must have only one channel for each observatory
+ if len(timeseries) != len(self.observatories):
+ raise AlgorithmException(
+ 'Expected data for %d stations, received %d \n'
+ 'Only 1 channel may be averaged at one time'
+ % (len(self.observatories), len(timeseries)))
+
+ # timeseries starttime and number of samples must match
+ for ts in timeseries:
+ # grab 1st set of stats to use in output.
+ # Its values will be good if checks pass.
+ if self._stats is None:
+ self._stats = ts.stats
+
+ if self._npts == -1:
+ self._npts = ts.stats.npts
+ if ts.stats.npts != self._npts:
+ raise AlgorithmException(
+ 'Received timeseries have different lengths')
+
+ if self._stt == -1:
+ self._stt = ts.stats.starttime
+ if ts.stats.starttime != self._stt:
+ raise AlgorithmException(
+ 'Received timeseries have different starttimes')
+
+ def process(self, timeseries):
+ """averages a channel across multiple stations
+
+ Parameters
+ ----------
+
+ Returns
+ -------
+ out_stream:
+ new stream object containing the averaged values.
+ """
+
+ self.check_stream(timeseries)
+
+ # initialize array for data to be appended
+ combined = []
+ # loop over stations
+ for obsy in self.observatories:
+
+ # lookup latitude correction factor, default = 1.0
+ latcorr = 1.0
+ if obsy in CORR:
+ latcorr = CORR[obsy]
+
+ # create array of data for each station
+ # and take into account correction factor
+ ts = timeseries.select(station=obsy)[0]
+ combined.append(ts.data * latcorr)
+
+ # after looping over stations, compute average
+ dst_tot = numpy.mean(combined, axis=0)
+
+ # Create a stream from the trace function
+ stream = obspy.core.Stream((
+ get_trace(self.outchannel, self._stats, dst_tot), ))
+
+ # set the full metadata for the USGS station used for averaged
+ # data sets
+ self.set_metadata(
+ stream=stream,
+ observatory='USGS',
+ channel=self.outchannel,
+ type=stream[0].stats.data_type,
+ interval=timeseries[0].stats.data_interval)
+
+ # return averaged values as a stream
+ return stream
+
+ def set_metadata(self, stream, observatory, channel, type, interval):
+ """set metadata for a given stream/channel
+ Parameters
+ ----------
+ observatory : str
+ observatory code
+ channel : str
+ edge channel code {MVH, MVE, MVD, ...}
+ type : str
+ data type {definitive, quasi-definitive, variation}
+ interval : str
+ interval length {minute, second}
+ """
+
+ for trace in stream:
+ self.observatoryMetadata.set_metadata(trace.stats, observatory,
+ channel, type, interval)
+
+ @classmethod
+ def add_arguments(cls, parser):
+ """Add command line arguments to argparse parser.
+
+ Parameters
+ ----------
+ parser: ArgumentParser
+ command line argument parser
+ """
+ parser.add_argument('--average-observatory-scale',
+ default=(None,),
+ help='Scale factor for observatories specified with ' +
+ '--observatory argument',
+ nargs='*',
+ type=float)
+
+ def configure(self, arguments):
+ """Configure algorithm using comand line arguments.
+
+ Parameters
+ ----------
+ arguments: Namespace
+ parsed command line arguments
+ """
+
+ self.observatories = arguments.observatory
+ if len(arguments.outchannels) > 1:
+ raise AlgorithmException(
+ 'Only 1 channel can be specified')
+ if arguments.outchannels:
+ self.outchannel = arguments.outchannels[0]
+ else:
+ self.outchannel = 'MSD'
+ self.scales = arguments.average_observatory_scale
+ if self.scales[0] is not None:
+ if len(self.observatories) != len(self.scales):
+ raise AlgorithmException(
+ 'Mismatch between observatories and scale factors')
+ else:
+ for (i, obs) in enumerate(self.observatories):
+ CORR[obs] = self.scales[i]
+
+
+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
+ stats.station = 'USGS'
+ stats.network = 'NT'
+ stats.location = 'R0'
+
+ return obspy.core.Trace(data, stats)
diff --git a/geomagio/algorithm/__init__.py b/geomagio/algorithm/__init__.py
index 2a6ac517a2a7ac0758f017f848407342b5b5376f..9b41729b35b2f010e7079fee586717af3cc71c06 100644
--- a/geomagio/algorithm/__init__.py
+++ b/geomagio/algorithm/__init__.py
@@ -8,6 +8,7 @@ from .Algorithm import Algorithm
from .AlgorithmException import AlgorithmException
# algorithms
from .AdjustedAlgorithm import AdjustedAlgorithm
+from .AverageAlgorithm import AverageAlgorithm
from .DeltaFAlgorithm import DeltaFAlgorithm
from .SqDistAlgorithm import SqDistAlgorithm
from .XYZAlgorithm import XYZAlgorithm
@@ -17,6 +18,7 @@ from .XYZAlgorithm import XYZAlgorithm
algorithms = {
'identity': Algorithm,
'adjusted': AdjustedAlgorithm,
+ 'average': AverageAlgorithm,
'deltaf': DeltaFAlgorithm,
'sqdist': SqDistAlgorithm,
'xyz': XYZAlgorithm
@@ -29,6 +31,7 @@ __all__ = [
'AlgorithmException',
# algorithms
'AdjustedAlgorithm',
+ 'AverageAlgorithm',
'DeltaFAlgorithm',
'SqDistAlgorithm',
'XYZAlgorithm'
diff --git a/setup.py b/setup.py
index edb67644460431478a7ec44ed1be080979f79b99..d74bc0f5f276d4e3d0e4fc45049efde86ef8d630 100644
--- a/setup.py
+++ b/setup.py
@@ -26,6 +26,6 @@ setup(
],
scripts=[
'bin/geomag.py',
- 'bin/make_cal.py'
+ 'bin/make_cal.py',
]
)