diff --git a/src/python/geomag/io/Timeseries.py b/src/python/geomag/io/Timeseries.py new file mode 100644 index 0000000000000000000000000000000000000000..3b3efd55e3eaa4b8d9e6d85ad12689530535d7fb --- /dev/null +++ b/src/python/geomag/io/Timeseries.py @@ -0,0 +1,139 @@ +import numpy +from obspy.core.utcdatetime import UTCDateTime +from TimeseriesException import TimeseriesException + + +class Timeseries(object): + """A sequence of samples at regular intervals. + + Attributes + ---------- + observatory : string + the observatory code + channels : dict + keys are channel names. + values are numpy arrays of timeseries values. + numpy.nan is used in arrays when values are missing. + starttime : UTCDateTime + time of the first sample. + endtime : UTCDateTime + time of the last sample. + metadata : array_like + array of metadata dict objects describing the timeseries. + metadata keys and values vary by factory. + when multiple metadata exist, starttime and endtime keys + describe metadata extents. + length : int + number of samples in each channel. + rate : float + sample rate in hertz. + """ + def __init__(self, observatory, channels, starttime, endtime, + metadata=None): + self.observatory = observatory + self.channels = channels + self.starttime = starttime + self.endtime = endtime + self.metadata = metadata or [] + self.length = len(channels[channels.keys()[0]]) + self.rate = (self.length - 1) / (endtime - starttime) + + def get_time(self, i): + """Get the time of a sample for the given index. + + Parameters + ---------- + i : int + index of the sample + + Returns + ------- + UTCDateTime + time of the i-th sample + """ + return UTCDateTime(self.starttime.timestamp + (i / self.rate)) + + def get_index(self, time, exact=True): + """Get the index of the sample for a given time. + + Parameters + ---------- + time : UTCDateTime + the given time. + exact : Boolean, optional + when True, return None if no sample exists at exactly ``time``. + when False, return time of nearest sample. + + Returns + ------- + int + index of sample for given time, or None. + """ + index = round((time - self.starttime) * self.rate) + if exact: + indexTime = self.get_time(index) + if indexTime != time: + return None + return index + + def __add__(self, other): + """Combine two timeseries objects. + + Parameters + ---------- + other : Timeseries + timeseries object to append to current object. + other must have the same observatory, same sampling rate, + same channels, and result in a continuous timeseries + (i.e. ``self.endtime + 1/rate == other.starttime``). + + Returns + ------- + Timeseries + new object with all samples of self and other. + + Raises + ------ + TimeseriesException + if the preconditions for ``other`` are not met. + """ + if self.observatory != other.observatory: + raise TimeseriesException('cannot add, different observatory') + if self.rate != other.rate: + raise TimeseriesException('cannot add, different rates') + if sorted(self.channels.keys()) != sorted(other.channels.keys()): + raise TimeseriesException('cannot add, different channels') + expected_start = self.get_time(self.length) + if expected_start != other.starttime: + print 'expected=', expected_start, 'start=', other.starttime + raise TimeseriesException('cannot add, non continuous timeseries') + merged_channels = {} + for key in self.channels.keys(): + merged_channels[key] = numpy.concatenate( + (self.channels[key], other.channels[key])) + return Timeseries(self.observatory, merged_channels, + self.starttime, other.endtime, self.metadata + other.metadata) + + def __len__(self): + """Get the number of samples in each channel. + + Returns + ------- + int + the number of samples in each channel. + """ + return self.length + + def __str__(self): + """Get a summary of this timeseries. + + Returns + ------- + str + summary of this timeseries + """ + return 'Timeseries ' + str(self.observatory) + \ + ' ' + str(self.channels.keys()) + \ + ' ' + str(self.length) + ' samples' + \ + ' at ' + str(self.rate) + ' hertz' + \ + ' (' +str(self.starttime) + ' - ' + str(self.endtime) + ')' diff --git a/src/python/geomag/io/TimeseriesException.py b/src/python/geomag/io/TimeseriesException.py new file mode 100644 index 0000000000000000000000000000000000000000..9fa2e718f71f0dcfccbddd8c882aa97a55c7ef9f --- /dev/null +++ b/src/python/geomag/io/TimeseriesException.py @@ -0,0 +1,7 @@ +""" +Exception thrown by various Timeseries methods. +""" + +class TimeseriesException(Exception): + """Exception thrown by various Timeseries methods.""" + pass diff --git a/src/python/geomag/io/TimeseriesFactory.py b/src/python/geomag/io/TimeseriesFactory.py new file mode 100644 index 0000000000000000000000000000000000000000..9e05d4f71cfd078ed39bc83ef3961cd12a5655f8 --- /dev/null +++ b/src/python/geomag/io/TimeseriesFactory.py @@ -0,0 +1,36 @@ +"""Abstract Timeseries Factory Interface.""" + + +class TimeseriesFactory(object): + """Base class for timeseries factories.""" + + def get_timeseries(self, observatory, starttime, endtime, + channels=('H', 'D', 'Z', 'F'), + type='variation', interval='minute'): + """Get timeseries data. + + Support for specific channels, types, and intervals varies + between factory and observatory. Subclasses should raise + TimeseriesFactoryException if the data is not available, or + if an error occurs accessing data. + + Parameters + ---------- + observatory : str + observatory code, usually 3 characters. + starttime : UTCDateTime + time of first sample in timeseries. + endtime : UTCDateTime + time of last sample in timeseries. + channels : array_like + list of channels to load. + type : {'definitive', 'provisional', 'quasi-definitive', 'variation'} + data type. + interval : {'daily', 'hourly', 'minute', 'monthly', 'second'} + data interval. + """ + raise NotImplementedError('"get_timeseries" not implemented') + + def put_timeseries(self, timeseries): + """Store timeseries data.""" + raise NotImplementedError('"put_timeseries" not implemented') diff --git a/src/python/geomag/io/TimeseriesFactoryException.py b/src/python/geomag/io/TimeseriesFactoryException.py new file mode 100644 index 0000000000000000000000000000000000000000..06acf4473255a4c8e478ca5e218b5b2959fc4596 --- /dev/null +++ b/src/python/geomag/io/TimeseriesFactoryException.py @@ -0,0 +1,8 @@ +""" +Base class for exceptions thrown by factories. +""" + + +class TimeseriesFactoryException(Exception): + """Base class for exceptions thrown by factories.""" + pass diff --git a/src/python/geomag/io/Timeseries_test.py b/src/python/geomag/io/Timeseries_test.py new file mode 100644 index 0000000000000000000000000000000000000000..efb19d994389b3c548878e82dd524372af6c36f8 --- /dev/null +++ b/src/python/geomag/io/Timeseries_test.py @@ -0,0 +1,141 @@ +"""Test the Timeseries class.""" + +from obspy.core.utcdatetime import UTCDateTime +from nose.tools import assert_equals +from Timeseries import Timeseries + + +def test_get_time(): + """geomag.io.Timeseries_test.test_get_time() + + Construct a timeseries object with 4 samples per channel, + and known starttime and endtime. + Verify get_time handles + - negative indices + - zero (starttime) + - length (endtime) + - beyond length + """ + timeseries = Timeseries('OBS', + { + 'H': [1, 2, 3, 4], + 'E': [2, 3, 4, 5], + 'Z': [3, 4, 5, 6], + 'F': [4, 5, 6, 7] + }, + UTCDateTime('2014-01-01T01:01:01'), + UTCDateTime('2014-01-01T01:01:04')) + assert_equals(UTCDateTime('2014-01-01T01:00:59'), + timeseries.get_time(-2)) + assert_equals(UTCDateTime('2014-01-01T01:01:01'), + timeseries.get_time(0)) + assert_equals(UTCDateTime('2014-01-01T01:01:05'), + timeseries.get_time(len(timeseries))) + assert_equals(UTCDateTime('2014-01-01T01:01:07'), + timeseries.get_time(6)) + +def test_get_index(): + """geomag.io.Timeseries_test.test_get_index() + + Construct a timeseries object with 4 samples per channel, + and known starttime and endtime. + Verify get_index handles + - times before starttime + - starttime + - endtime + - times after endtime + """ + timeseries = Timeseries('OBS', + { + 'H': [1, 2, 3, 4], + 'E': [2, 3, 4, 5], + 'Z': [3, 4, 5, 6], + 'F': [4, 5, 6, 7] + }, + UTCDateTime('2014-01-01T01:01:01'), + UTCDateTime('2014-01-01T01:01:04')) + assert_equals(-2, + timeseries.get_index(UTCDateTime('2014-01-01T01:00:59'), False)) + assert_equals(0, + timeseries.get_index(UTCDateTime('2014-01-01T01:01:01'), False)) + assert_equals(len(timeseries), + timeseries.get_index(UTCDateTime('2014-01-01T01:01:05'), False)) + assert_equals(6, + timeseries.get_index(UTCDateTime('2014-01-01T01:01:07'), False)) + +def test_get_index_exact(): + """geomag.io.Timeseries_test.test_get_time() + + Construct a timeseries object with 4 samples per channel, + and known starttime and endtime. + Verify get_time returns + - None when time is between samples + - index when time is at sample + """ + timeseries = Timeseries('OBS', + { + 'H': [1, 2, 3, 4], + 'E': [2, 3, 4, 5], + 'Z': [3, 4, 5, 6], + 'F': [4, 5, 6, 7] + }, + UTCDateTime('2014-01-01T01:01:01'), + UTCDateTime('2014-01-01T01:01:04')) + assert_equals(None, + timeseries.get_index(UTCDateTime('2014-01-01T01:00:59.5'), True)) + assert_equals(1, + timeseries.get_index(UTCDateTime('2014-01-01T01:01:02'), True)) + +def test_length(): + """geomag.io.Timeseries_test.test_length() + + Construct a timeseries object with 4 samples per channel. + Verify ``len(timeseries)`` returns 4 + """ + timeseries = Timeseries('OBS', + { + 'H': [1, 2, 3, 4], + 'E': [2, 3, 4, 5], + 'Z': [3, 4, 5, 6], + 'F': [4, 5, 6, 7] + }, + UTCDateTime('2014-01-01T01:01:01'), + UTCDateTime('2014-01-01T01:01:04')) + assert_equals(4, len(timeseries)) + + +def test_rate_seconds(): + """geomag.io.Timeseries_test.test_rate_seconds() + + Construct a timeseries with 4 samples, and starttime and endtime + at 01 and 04 seconds within the same minute. + Verify ``timeseries.rate`` is 1 (hertz). + """ + timeseries = Timeseries('OBS', + { + 'H': [1, 2, 3, 4], + 'E': [2, 3, 4, 5], + 'Z': [3, 4, 5, 6], + 'F': [4, 5, 6, 7] + }, + UTCDateTime('2014-01-01T01:01:01'), + UTCDateTime('2014-01-01T01:01:04')) + assert_equals(1, timeseries.rate) + +def test_rate_minutes(): + """geomag.io.Timeseries_test.test_rate_minutes() + + Construct a timeseries with 4 samples, and starttime and endtime + at 01 and 04 minutes within the same same. + Verify ``timeseries.rate`` is 1/60 (hertz). + """ + timeseries = Timeseries('OBS', + { + 'H': [1, 2, 3, 4], + 'E': [2, 3, 4, 5], + 'Z': [3, 4, 5, 6], + 'F': [4, 5, 6, 7] + }, + UTCDateTime('2014-01-01T01:01:00'), + UTCDateTime('2014-01-01T01:04:00')) + assert_equals(1.0 / 60.0, timeseries.rate)