This notebook demonstrates some DaViTPy tools for working with timeseries data. It includes filtering and visualization tools.
The DaViTPy signal is used by creating a VTSig object, which includes built-in methods for keeping track of all steps in a signal processing algorithm. As the signal is processed, all preceeding variants of the signal are saved within the object and a history is created. This makes it very easy to see what has been done with the data at any step along the way.
Written by N.A. Frissell, 9 July 2013.
In [1]:
#Import the modules we need.
import numpy as np
import scipy as sp
import matplotlib.pyplot as mp
import datetime
import copy
import pydarn
import utils
import pydarn.proc.signal as signal
In [2]:
#This cell just sets some preferences to make the plots look nice.
font = {'family' : 'sans-serif',
'weight' : 'bold',
'size' : 22}
matplotlib.rc('font', **font)
figure = {'figsize' : '15, 12',
'dpi' : 600}
matplotlib.rc('figure', **figure)
We will start with some GOES10 magnetometer data since it is regularly sampled and easy to load from the provided sample data file.
In [3]:
#Metadata is used to control aspects of signal processing and plotting throughout the routine, as well as a way
#of keeping track of different data properties and processing done to the data.
#Global metadata will control parameters pertaining to all signals being worked on.
#Here, we define the time period we want any processing we do to apply to.
signal.globalMetaData_add(validTimes = [datetime.datetime(2008,07,14,16,00), datetime.datetime(2008,07,15,00,00)])
In [4]:
#Read data from file and sort out variables.
satfile = 'signal-sample_data-goes10.txt'
dataIn = np.genfromtxt(satfile,comments='#',usecols=(0,1,2),dtype=[('date','S10'),('time','S5'),('btot','f8')])
timeStr = map(' '.join, zip(dataIn['date'],dataIn['time']))
datetimeObjList = [ datetime.datetime.strptime(x,'%Y-%m-%d %H:%M') for x in timeStr ]
Next we instantiate vt sig object. All you need is a list of Python datetime.datetime objects and the corresponding data.
In [5]:
#Instantiate vt sig object
goes10 = signal.sig(datetimeObjList,dataIn['btot'])
In [6]:
#Now lets look to see where the data is actually stored...
goes10.raw.dtv #Datetime vector (dtv)
Out[6]:
In [7]:
goes10.raw.data #Data vector
Out[7]:
Now we set the object metadata. This not only allows you to keep track of information about the date, but is also used by the plotting and signal processing routines.
goes10.metadata applies to all versions of the signal stored in the goes10 object. goes10.[dataSetName].metadata applies just to that data set. Anything set in goes10.[dataSetName].metadata will always override anything set in goes10.metadata. You can use the goes10.[dataSetName].getAllMetaData() method to see what metadata is actually being applied to a particular dataset.
In [8]:
goes10.metadata['title'] = 'GOES 10'
goes10.metadata['ylabel'] = 'B$_{tot}$ [nT]'
goes10.metadata['xlabel'] = 'Time [UT]\n14 July 2008'
#goes10.metadata['xmin'] = datetime.datetime(2008,07,14,18,00)
#goes10.metadata['xmax'] = datetime.datetime(2008,07,15,00,00)
#goes10.raw.metadata['ymin'] = 100.
#goes10.raw.metadata['ymax'] = 125.
In [9]:
#Now let's plot our raw data. This plot shows all of the GOES10 data which has been loaded into the object.
#The white region is the region defined by the validTimes in the global metadata.
goes10.raw.plot()
Let's run some processing now. Remember, each step does not discard the data from the previous steps.
In [10]:
#Apply a linear detrend through the data.
signal.detrend(goes10.raw)
In [11]:
#Here is the result of the detrending. Note that only the validTimes period is shown now. This is because
#the detrending was applied only to the validTimes period.
goes10.detrended.plot()
In [12]:
#Global metadata can also be used to define characteristics of a filter.
signal.globalMetaData_add(filter_numtaps = 101) #This is the length of the filter. A larger number is a better filter, but slower computationally.
signal.globalMetaData_add(filter_cutoff_high = 0.004) #Frequencies are in Hz.
signal.globalMetaData_add(filter_cutoff_low = 0.001)
In [13]:
#Now filter the data with filter we defined in the global metadata.
goes10Filt = signal.filter(goes10)
In [14]:
#Let's look at the current history of our signal object.
goes10.active.history
Out[14]:
In [15]:
#We can also see all of the metadata set for a particular dataObject.
goes10.active.getAllMetaData()
Out[15]:
In [16]:
#Let's see the transfer function of the filter we used.
goes10Filt.plotTransferFunction()
In [17]:
#And now the filtered result. Because the filter causes a delay in the signal and doesn't get to the end of
#the signal, part of the filtered signal should be discard. Appropriate new validTimes have been set, and the
#gray bars at the beginning and end mark the part of the signal to be discarded.
#Note that the validTimes here only apply to this filtered signal (goes10.filtered.metadata['validTimes']), and
#not the global metadata. When multiple validTimes are in place for a particular object, the most restrictive set of
#validTimes is always used.
goes10.filtered.plot()
In [18]:
#We can also look at the spectrum of any signal that is regularly sampled in time.
goes10.filtered.plotfft()
In [19]:
#Once you plot the FFT (or invoke the fft() method), you now have access to the spectrum
goes10.filtered.freqVec #Here are the frequency bins in Hz.
Out[19]:
In [20]:
goes10.filtered.spectrum #Here is the actual spectrum
Out[20]:
In [21]:
#We can also compare the spectra of two different signals.
a = signal.oplotfft([goes10.detrended,goes10.filtered])
Now lets try looking at some SuperDARN data from the same period. SuperDARN can be easily accessed through Virginia Tech's servers, as shown below. Working with this data is a little trickier because it is not necessarily regularly sampled.
In [22]:
#Set some parameters regarding the event of interest. Let's use a time series from beam 7, rangegate 0.
sTime = datetime.datetime(2008,7,14)
eTime = datetime.datetime(2008,7,14,23)
fileType = 'fitacf'
radar = 'rkn'
beam = 7
gate = 0
In [23]:
#Use this cell to load in data straight from DaViTPy. See the radDataRead notebook for
#additional inspiration on reading SuperDARN data.
myPtr = pydarn.sdio.radDataOpen(sTime,radar,eTime=eTime,bmnum=beam,fileType=fileType)
datetimeObjList = []
data = []
time_ii = copy.copy(sTime)
myBeam = 1
while time_ii < eTime:
myBeam = pydarn.sdio.radDataReadRec(myPtr)
if myBeam == None:
break
time_ii = myBeam.time
if not gate in myBeam.fit.slist:
continue
datetimeObjList.append(time_ii)
inx = int(np.where(np.array(myBeam.fit.slist) == gate)[0])
data.append(myBeam.fit.v[inx])
In [24]:
#Instantiate vt sig object.
rkn = signal.sig(datetimeObjList,data)
In [25]:
#Set Rankin Inlet Radar Metadata.
rkn.metadata['title'] = ' '.join([radar.upper(),'Beam:',str(beam),'Gate:',str(gate)])
rkn.metadata['ylabel'] = 'Velocity [m/s]'
rkn.metadata['xlabel'] = 'Time [UT]\n14 July 2008'
#rkn.metadata['lineStyle'] = '.-'
#rkn.active.metadata['xmin'] = datetime.datetime(2008,07,14,20,30)
#rkn.active.metadata['xmax'] = datetime.datetime(2008,07,14,23,30)
#rkn.metadata['validTimes'] = [datetime.datetime(2008,07,14,01,00), datetime.datetime(2008,07,15,01,00)]
In [26]:
#Here is a way to change the valid times for the raw dataset. Just uncomment the lines below and give it a try!
#valid = [rkn.raw.dtv[0], rkn.raw.dtv[-1]]
#rkn.raw.metadata['validTimes'] = valid
In [27]:
rkn.raw.plot()
In [28]:
#This method will tell us the sampling period of a dataset.
#If it is not regularly sampled, it will give us a warning.
rkn.raw.samplePeriod()
Out[28]:
In [29]:
#We can do a linear interpolation the signal to make it regularly sampled.
#This does need to be done with care... we will compare the raw and interpolated signal
#in the next cell.
signal.interpolate(rkn.raw,samplePeriod=3.)
In [30]:
#Now compare the two signals. You can see that the interpolated version misses parts of the raw signal.
signal.oplot([rkn.raw,rkn.active],xmin=datetime.datetime(2008,7,14,17))
Out[30]:
In [31]:
#Again, detrending and filtering. We are overriding the global metadata set before by using the numtaps keyword.
signal.detrend(rkn.active)
rknFilt = signal.filter(rkn.detrended,numtaps=1001)
In [32]:
#Print out some information about the filter being used.
rknFilt.comment
Out[32]:
In [33]:
#Plot the transfer function. The worN keyword computes the transfer function at a higher resolution than the
#default 512/(2*Pi). See plotTransferFunction docstring for details.
rknFilt.plotTransferFunction(xmax=0.009,worN=5000)
In [34]:
#Look at the filtered data.
rkn.filtered.plot(xmin=datetime.datetime(2008,7,14,17))
In [35]:
#Compare the filtered with the non-filetered spectrum.
a = signal.oplotfft([rkn.detrended,rkn.filtered],xmax=0.01)
We can now compare the Rankin Inlet data with the satellite data.
In [36]:
#The truncate method removes any part of the signal that is outside of the validTimes.
goes10.filtered.truncate()
rkn.filtered.truncate()
Out[36]:
In [37]:
#Now we overplot the time series. There is still some gray since the data sets do not overlap perfectly in time.
oplot = signal.oplot([goes10,rkn],normalize=True,legend_size=10)
In [38]:
#We can also compare the spectra of the two signals.
oplot = signal.oplotfft([goes10,rkn],normalize=True,legend_size=12,xmax=0.010)
In [39]:
#Finally, it is possible to place everything onto an identical time grid. The most restrictive range of validtimes and the highest time resolution is used.
signal.commonDtv([goes10.active,rkn.active])
In [40]:
#Overplot two signals. The amplitude of goes10 is small compare to rkn, so it is hard to really make anything out here.
#Note in the command below, no dataset is specified, so the active one is used by default.
signal.oplot([goes10,rkn])
Out[40]:
In [41]:
#Don't forget we can look at the history of an object!
goes10.active.history
Out[41]:
In [42]:
rkn.active.history
Out[42]:
In [43]:
#xcor = signal.xcor(rkn,goes10)
In [44]:
#xcor.xcor.plot()