#!/usr/bin/env python
# coding: utf-8

# # Signal Processing Notebook

# 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.
get_ipython().run_line_magic('pylab', 'inline')
import numpy as np
import scipy as sp
import matplotlib.pyplot as mp
import datetime
import copy

from davitpy import pydarn
from davitpy import utils
import davitpy.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)


# # Signal Processing Example - GOES10 Data

# We will start with some GOES10 magnetometer data since it is regularly sampled and easy to load from the provided sample data file.

# ## Set Global Metadata

# 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)])


# ## Load Data

# 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.
# 
# * The goes10 object will have a number of subobjects and methods.
# * goes10.raw is subobject that contains the raw data of datetimeObjList,data.
# * goes10.active should always point to the most processed version of the data, unless you manually set it to point to something else.
# * When you apply a processing algorthim to the goes10 object, it will use the "active" dataset, unless you specify otherwise. 

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


# In[7]:


goes10.raw.data #Data vector


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


# ## Process Data

# 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


# In[15]:


#We can also see all of the metadata set for a particular dataObject.
goes10.active.getAllMetaData()


# 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.


# In[20]:


goes10.filtered.spectrum #Here is the actual spectrum


# In[21]:


#We can also compare the spectra of two different signals.
a = signal.oplotfft([goes10.detrended,goes10.filtered])


# # Signal Processing Example - SuperDARN Rankin Inlet Data

# 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.

# ## Load Data

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


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


# 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


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


# # Overplotting Data

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


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


# # Common Dtv Grid

# 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])


# In[41]:


#Don't forget we can look at the history of an object!
goes10.active.history


# In[42]:


rkn.active.history


# ## In the future, we can have cross correlation tools, too.  These commands are almost working, but not quite...

# In[ ]:


#xcor = signal.xcor(rkn,goes10)


# In[ ]:


#xcor.xcor.plot()