schain Commit - r1006:8363978b5396 · Jicamarca Repository

22/08/2017

ebocanegra -

r1006:8363978b5396

parent child

Context file:

r1006:8363978b5396

Expand all files

.gitignore +3 -6

-            <<<<<<< HEAD
             # Byte-compiled / optimized / DLL files
             __pycache__/
             *.py[cod]
             *$py.class
             # C extensions
             *.so
             # Distribution / packaging
             .Python
             env/
             build/
             develop-eggs/
             dist/
             downloads/
             eggs/
             .eggs/
             lib/
             lib64/
             parts/
             sdist/
             var/
             wheels/
             *.egg-info/
             .installed.cfg
             *.egg
             # PyInstaller
             #  Usually these files are written by a python script from a template
             #  before PyInstaller builds the exe, so as to inject date/other infos into it.
             *.manifest
             *.spec
             # Installer logs
             pip-log.txt
             pip-delete-this-directory.txt
             # Unit test / coverage reports
             htmlcov/
             .tox/
             .coverage
             .coverage.*
             .cache
             nosetests.xml
             coverage.xml
             *,cover
             .hypothesis/
             # Translations
             *.mo
             *.pot
             # Django stuff:
             *.log
             local_settings.py
             # Flask stuff:
             instance/
             .webassets-cache
             # Scrapy stuff:
             .scrapy
             # Sphinx documentation
             docs/_build/
             # PyBuilder
             target/
             # Jupyter Notebook
             .ipynb_checkpoints
             # pyenv
             .python-version
             # celery beat schedule file
             celerybeat-schedule
             # SageMath parsed files
             *.sage.py
             # dotenv
             .env
             # virtualenv
             .venv
             venv/
             ENV/
             # Spyder project settings
             .spyderproject
             .spyproject
             # Rope project settings
             .ropeproject
             # mkdocs documentation
             /site
             # eclipse
             .project
             .pydevproject
-            =======
-            build/
-            dist/
-            schainpy.egg-info/
             .svn/
-            >>>>>>> 08c4507d6c3c48f6c52326d5dedfa1972fb26356
+            *.png
+            *.pyc
+            schainpy/scripts

schainpy/admin.pyc binary modified 0 0

NO CONTENT: modified file, binary diff hidden

schainpy/controller.pyc binary modified 0 0

NO CONTENT: modified file, binary diff hidden

schainpy/model/data/jrodata.pyc binary modified 0 0

	1	NO CONTENT: modified file, binary diff hidden			NO CONTENT: modified file, binary diff hidden

NO CONTENT: modified file, binary diff hidden

schainpy/model/data/jroheaderIO.pyc binary modified 0 0

	1	NO CONTENT: modified file, binary diff hidden			NO CONTENT: modified file, binary diff hidden

NO CONTENT: modified file, binary diff hidden

schainpy/model/graphics/jroplot_correlation.pyc binary modified 0 0

NO CONTENT: modified file, binary diff hidden

schainpy/model/graphics/jroplot_heispectra.pyc binary modified 0 0

NO CONTENT: modified file, binary diff hidden

schainpy/model/graphics/jroplot_parameters.pyc binary modified 0 0

NO CONTENT: modified file, binary diff hidden

schainpy/model/graphics/jroplot_spectra.pyc binary modified 0 0

NO CONTENT: modified file, binary diff hidden

schainpy/model/graphics/jroplot_voltage.pyc binary modified 0 0

NO CONTENT: modified file, binary diff hidden

schainpy/model/graphics/mpldriver.pyc binary modified 0 0

NO CONTENT: modified file, binary diff hidden

schainpy/model/io/__init__.py +1 -1

             '''
             $Author: murco $
             $Id: JRODataIO.py 169 2012-11-19 21:57:03Z murco $
             '''
             from jroIO_voltage import *
             from jroIO_spectra import *
             from jroIO_heispectra import *
             from jroIO_usrp import *
             from jroIO_kamisr import *
             from jroIO_param import *
             from jroIO_hf import *
             from jroIO_bltr import *
             from jroIO_mira35c import *
-            #from io_bltr_block import *
+            from io_bltr_block import *

schainpy/model/io/__init__.pyc binary modified 0 0

NO CONTENT: modified file, binary diff hidden

schainpy/model/io/io_bltr_block.py +8 -6

             '''
             Created on Nov 9, 2016
             @author: roj- LouVD
             '''
             import numpy
             import os.path
             import sys
             import time
             import datetime
             from sys import path
             from os.path import dirname
             from mimify import HeaderFile
             from numpy import size, asarray
             from datetime import datetime
             from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation
             from schainpy.model.data.jrodata import Parameters
             from schainpy.model.data.jroheaderIO import RadarControllerHeader, SystemHeader
             from schainpy.model.graphics.jroplot_parameters import WindProfilerPlot
             from schainpy.model.io.jroIO_base import *
             import schainpy
             #import madrigal
             #import madrigal.cedar
             #from madrigal.cedar import MadrigalCatalogRecord
             import warnings
             from time import gmtime
             from math import floor
             warnings.simplefilter("error")
             from numpy.lib.nanfunctions import nansum
             warnings.simplefilter('ignore', FutureWarning)
             class testBLTRReader(ProcessingUnit):
-                def __init__(self):
+                def __init__(self, **kwargs):
                     path = None
                     startDate = None
                     endDate = None
                     startTime = None
                     endTime = None
                     startTime = None
                     endTime = None
                     isConfig = False
                     dataOut = None
                     walk = None
                     ext = 'swwma'
                     fileList = []
                     fileIndex = -1
                     timezone = None
                     filename = None
                     timearray = None
                     height = None
                     snr_ref = None
                     zon_ref = None
                     ver_ref = None
                     mer_ref = None
                     nmodes = None
                     nchannels = None
                     nranges = None
                     year = None
                     month = None
                     day = None
                     lat = None
                     lon = None
                     siteFile = None
-                    ProcessingUnit.__init__(self)
+                    ProcessingUnit.__init__(self , **kwargs)
                     self.dataOut = self.createObjByDefault()
                     self.imode = 0
                     self.counter_records = 0
                     self.isConfig = False
                     self.flagNoMoreFiles = 0
                     self.buffer = None
                 def createObjByDefault(self):
                     dataObj = Parameters()
                     return dataObj
                 def info(self):
                     '''
                     Experience information
                     '''
                     self.hoy = datetime.datetime.now()
                     place = 'Jicamarca Radio Observatory'
                     signalchainweb='http://jro-dev.igp.gob.pe:3000/projects/signal-chain/wiki/Manual_de_Desarrollador'
                     print '{} at {}'.format(self.hoy,place)
                     print 'Boundary Layer and Tropospheric Radar (BLTR) script, Wind velocities and SNR from *.sswma files'
                     print '{}  \n'.format(signalchainweb)
-                def run(self, path, startDate, endDate, ext, startTime, endTime):
+                def run(self, path, startDate, endDate, ext, startTime, endTime, queue=None):
                     if not(self.isConfig):
                         self.setup(path, startDate, endDate, ext)
                         self.isConfig = True
                     self.getData()
                 def setup(self,
                           path=None,
                           startDate=None,
                           endDate=None,
                           ext=None,
                           startTime=datetime.time(0, 0, 0),
                           endTime=datetime.time(23, 59, 59),
                           timezone=0):
                     self.info()
                     self.path = path
                     if self.path == None:
                         raise ValueError, "The path is not valid"
                     if ext == None:
                         ext = self.ext
                     self.searchFiles(self.path, startDate, endDate, ext)
                     self.timezone = timezone
                     self.ext = ext
                     self.fileIndex = -1
                     if not(self.fileList):
                         raise  Warning, "There is no files matching these date in the folder: %s. \n Check 'startDate' and 'endDate' "%(path)
                     if not(self.setNextFile()):
                         print 'not next file'
                         if (startDate!=None) and (endDate!=None):
                             print "No files in range: %s - %s" %(datetime.datetime.combine(startDate,startTime).ctime(), datetime.datetime.combine(endDate,endTime).ctime())
                         elif startDate != None:
                             print "No files in range: %s" %(datetime.datetime.combine(startDate,startTime).ctime())
                         else:
                             print "No files"
                         sys.exit(-1)
                 def searchFiles(self, path, startDate, endDate, ext=None):
                     '''
                      Searching for BLTR rawdata file in path
                      Creating a list of file to proces included in [startDate,endDate]
                      Input:
                          path - Path to find BLTR rawdata files
                          startDate - Select file from this date
                          enDate - Select file until this date
                          ext - Extension of the file to read
                     '''
                     fullpath = path
                     foldercounter = 0
                     print 'Searching file in %s ' % (fullpath)
                     fileList0 = glob.glob1(fullpath, "*%s" % ext)
                     fileList0.sort()
                     self.fileList = []
                     self.dateFileList = []
                     for thisFile in fileList0:
                         year = thisFile[-14:-10]
                         if not isNumber(year):
                             continue
                         month = thisFile[-10:-8]
                         if not isNumber(month):
                             continue
                         day = thisFile[-8:-6]
                         if not isNumber(day):
                             continue
                         year, month, day = int(year), int(month), int(day)
                         dateFile = datetime.date(year, month, day)
                         if not ((startDate <= dateFile) and (endDate > dateFile)):
                             continue
                         self.fileList.append(thisFile)
                         self.dateFileList.append(dateFile)
                     return 1
                 def setNextFile(self):
                     idFile = self.fileIndex
                     while (True):
                         idFile += 1
                         if idFile >= len(self.fileList):
                             print '\nNo more files in the folder'
                             print 'Total number of file(s) read : {}'.format(self.fileIndex + 1)
                             print  'Time of processing : {}'.format(datetime.datetime.now()- self.hoy)
                             self.flagNoMoreFiles = 1
                             return 0
                         if self.isConfig: print '------------------------[Next File]---------------------------'
                         filename = os.path.join(self.path, self.fileList[idFile])
                         self.Open(filename)
                         print '\n[Setting file] (%s) ...' % self.fileList[idFile]
                         break
                     self.flagIsNewFile =0
                     self.fileIndex = idFile
                     self.filename = filename
                     print 'File:',self.filename
                     return 1
                 def readDataBlock(self):
                     self.readHeader()
                     self.dataRecords(0)
                     print '[New Record] record: {} /{} // file {}/{}'.format(self.counter_records,self.nrecords,self.fileIndex+1,len(self.fileList))
                     self.setDataBuffer()
                     self.flagIsNewBlock = 1
                     if self.counter_records > self.nrecords:
                         self.flagIsNewFile = 1
                         return 0
                     return 1
                 def setDataBuffer(self):
                     '''
                     Storing data from one block
                     '''
                     self.t = datetime.datetime(self.year, self.month, self.day)
                     self.doy = time.localtime(time.mktime(self.t.timetuple())).tm_yday
                     self.buffer = numpy.squeeze(numpy.array([[self.one_snr],[self.one_zonal],[self.one_vertical],[self.one_meridional],
                                                 [self.time],[self.height],[self.fileIndex],
                                                 [self.year],[self.month],[self.day],[self.t],[self.doy]]))
                     self.dataOut.time1 = self.time1
                 def Open(self, filename):
                     '''
                     Opening BLTR rawdata file defined by filename
                     Inputs:
                         filename    - Full path name of BLTR rawdata file
                     '''
                     [dir, name] = os.path.split(filename)
                     strFile = name.split('.')
                     self.siteFile = strFile[0]  # 'peru2' ---> Piura  -   'peru1' ---> Huancayo or Porcuya
                     self.filename = filename
                     if os.path.isfile(self.filename) == False:
                         print 'File do not exist. Check "filename"'
                         sys.exit(0)
                     self.h_file = numpy.dtype([
                                       ('FMN', '<u4'),
                                       ('nrec', '<u4'),
                                       ('fr_offset', '<u4'),
                                       ('id', '<u4'),
                                       ('site', 'u1', (32,))
                                       ])
                     self.pointer = open(self.filename, 'rb')  # rb : Read Binary
+                    print self.filename
                     self.header_file = numpy.fromfile(self.pointer, self.h_file, 1)
+                    print self.header_file
                     self.nrecords = self.header_file['nrec'][0]
                     self.sizeOfFile = os.path.getsize(self.filename)
                     self.time = numpy.zeros([2, self.nrecords], dtype='u4')
                     self.counter_records = 0
                     self.count = 0
                     self.flag_initialArray = False
                     self.year = 0
                     self.month = 0
                     self.day = 0
                 def hasNotDataInBuffer(self):
                     if self.buffer == None:
                         return 1
                     return 0
                 def getData(self):
                     '''
                     Storing data from databuffer to dataOut object
                     '''
                     if self.flagNoMoreFiles==1:
                         self.dataOut.flagNoData = True
                         print 'No file left to process'
                         return 0
                     self.flagIsNewBlock = 0
                     if self.hasNotDataInBuffer():
                         if self.flagIsNewFile==0:
                             self.readNextBlock()
                             '''RETURN A BLOCK OF DATA'''
                             if self.flagNoMoreFiles==0:
                                 self.dataOut.data_SNR = self.buffer[0]
                                 self.dataOut.time = self.buffer[4]
                                 self.dataOut.height = self.height
                                 self.dataOut.height= self.height
                                 self.dataOut.data_output = numpy.squeeze(numpy.array([[self.buffer[1]],
                                                                                     [self.buffer[3]],
                                                                                     [self.buffer[2]]]))
             #
                                 self.dataOut.day, self.dataOut.month, self.dataOut.year = self.buffer[9], self.buffer[8], self.buffer[7]
                                 self.dataOut.utctimeInit = self.time1
                                 self.dataOut.utctime = self.dataOut.utctimeInit
                                 self.dataOut.counter_records = self.counter_records
                                 self.dataOut.nrecords = self.nrecords
                                 self.setHeader()
                                 self.buffer = None
                                 self.dataOut.flagNoData = False
                 def readNextBlock(self):
                     if not(self.setNewBlock()):
                         return 0
                     if not(self.readDataBlock()):
                         return 0
                     if self.flagIsNewFile:
                         self.setNextFile()
                     return 1
                 def setNewBlock(self):
                     if self.pointer==None:
                         return 0
                     if self.flagIsNewFile:
                         return 1
                     if self.counter_records < self.nrecords:
                         return 1
                     if not(self.setNextFile()):
                         return 0
                     return 1
                 def readHeader(self):
                     '''
                     RecordHeader of BLTR rawdata file
                     '''
                     if self.pointer.tell() == self.sizeOfFile:
                         print 'End of File'
                         return
                     self.h_rec1 = numpy.dtype([
                                               ('rmn', '<u4'),
                                               ('rcounter', '<u4'),
                                               ('nr_offset', '<u4'),
                                               ('tr_offset', '<u4'),
                                               ('time', '<u4'),
                                               ('time_msec', '<u4'),
                                               ('tag', 'u1', (32,)),
                                               ('comments', 'u1', (32,)),
                                               ('lat', '<f4'),
                                               ('lon', '<f4'),
                                               ('gps_status', '<u4'),
                                               ('freq', '<u4'),
                                               ('freq0', '<u4'),
                                               ('nchan', '<u4'),
                                               ('delta_r', '<u4'),
                                               ('nranges', '<u4'),
                                               ('r0', '<u4'),
                                               ('prf', '<u4'),
                                               ('ncoh', '<u4'),
                                               ('npoints', '<u4'),
                                               ('polarization', '<i4'),
                                               ('rx_filter', '<u4'),
                                               ('nmodes', '<u4'),
                                               ('dmode_index', '<u4'),
                                               ('dmode_rngcorr', '<u4'),
                                               ('nrxs', '<u4'),
                                               ('acf_length', '<u4'),
                                               ('acf_lags', '<u4'),
                                               ('sea_to_atmos', '<f4'),
                                               ('sea_notch', '<u4'),
                                               ('lh_sea', '<u4'),
                                               ('hh_sea', '<u4'),
                                               ('nbins_sea', '<u4'),
                                               ('min_snr', '<f4'),
                                               ('min_cc', '<f4'),
                                               ('max_time_diff', '<f4')
                                               ])
                     self.header_rec1 = numpy.fromfile(self.pointer, self.h_rec1, 1)
                     self.lat = self.header_rec1['lat'][0]
                     self.lon = self.header_rec1['lon'][0]
                     self.nchannels = self.header_rec1['nchan'][0] / 2
                     self.kchan = self.header_rec1['nrxs'][0]
                     self.nranges = self.header_rec1['nranges'][0]
                     self.deltha = self.header_rec1['delta_r'][0]
                     self.correction = self.header_rec1['dmode_rngcorr'][0]
                     self.nmodes = self.header_rec1['nmodes'][0]
                     self.imode = self.header_rec1['dmode_index'][0]
                     self.h_rec2 = numpy.dtype([
                                       ('antenna_coord', 'f4', (2, self.nchannels)),
                                       ('rx_gains', 'u4', (self.nchannels,)),
                                       ('rx_analysis', 'u4', (self.nchannels,))
                                       ])
                     self.header_rec2 = numpy.fromfile(self.pointer, self.h_rec2, 1)  # header rec2
                     self.antenna = self.header_rec2['antenna_coord']
                     self.rx_gains = self.header_rec2['rx_gains']
                     self.d_rec = numpy.dtype ([
                                     ('range', '<u4'),
                                     ('status', '<u4'),
                                     ('zonal', '<f4'),
                                     ('meridional', '<f4'),
                                     ('vertical', '<f4'),
                                     ('zonal_a', '<f4'),
                                     ('meridional_a', '<f4'),
                                     ('corrected_fading', '<f4'),  # seconds
                                     ('uncorrected_fading', '<f4'),  # seconds
                                     ('time_diff', '<f4'),
                                     ('major_axis', '<f4'),
                                     ('axial_ratio', '<f4'),
                                     ('orientation', '<f4'),
                                     ('sea_power', '<u4'),
                                     ('sea_algorithm', '<u4'),
                                     ('rx_saturation', 'u4', (self.nchannels,)),
                                     ('chan_offset', 'u4', (2 * self.nchannels,)),
                                     ('rx_amp', 'u4', (self.nchannels,)),
                                     ('rx_snr', 'f4', (self.nchannels,)),
                                     ('cross_snr', 'f4', (self.kchan,)),
                                     ('sea_power_relative', 'f4', (self.kchan,))
                             ])
                     # Memory allocation
                     if not(self.flag_initialArray):
                         self.height = numpy.zeros([2, self.nranges], dtype='f4') + numpy.nan
                         self.p_zonal = numpy.zeros([self.nrecords, self.nranges, 2], dtype='f4') + numpy.nan
                         self.p_meridional = numpy.zeros([self.nrecords, self.nranges, 2], dtype='f4') + numpy.nan
                         self.p_vertical = numpy.zeros([self.nrecords, self.nranges, 2], dtype='f4') + numpy.nan
                         self.p_snr = numpy.zeros([self.nrecords, self.nranges, self.kchan, 2], dtype='f4') + numpy.nan
                         self.flag_initialArray = True
                     self.time[self.imode, self.count] = self.header_rec1['time'][0]
                     self.time1 = self.header_rec1['time'][0]
                     tseconds = self.header_rec1['time'][0]
                     local_t1 = time.localtime(tseconds)
                     self.year = local_t1.tm_year
                     self.month = local_t1.tm_mon
                     self.day = local_t1.tm_mday
                     self.t = datetime.datetime(self.year, self.month, self.day)
                 def setHeader(self):
                     '''
                     Saving metada on dataOut object
                     '''
                     self.dataOut.type = 'Parameters'
                     self.dataOut.useLocalTime = False
-                    self.dataOut.outputInterval = 157
+                    # self.dataOut.outputInterval = 157
+                    self.dataOut.paramInterval = 157
                     self.dataOut.timezone = self.timezone
                     self.dataOut.site = self.siteFile
                     self.dataOut.nrecords = self.nrecords
                     self.dataOut.sizeOfFile = self.sizeOfFile
                     self.dataOut.lat = self.lat
                     self.dataOut.lon = self.lon
                     self.dataOut.nchannels = self.nchannels
                     self.dataOut.kchan = self.kchan
                     self.dataOut.nranges = self.nranges
                     self.dataOut.deltha = self.deltha
                     self.dataOut.correction = self.correction
                     self.dataOut.nmodes = self.nmodes
                     self.dataOut.imode = self.imode
                     self.dataOut.antenna = self.antenna
                     self.dataOut.rx_gains = self.rx_gains
                 def dataRecords(self, status_value):
                     '''
                     Reading and filtering data block record of BLTR rawdata file, filtering is according to status_value.
                     Input:
                         status_value - Array data is set to NAN for values that are not equal to status_value
                     '''
                     data_rec = numpy.fromfile(self.pointer, self.d_rec, self.nranges)
                     status = []
                     zonal = []
                     meridional = []
                     vertical = []
                     rx_snr = []
                     index = 0
                     for rec in data_rec:
                         status.append(rec['status'])
                         zonal.append(rec['zonal'])
                         meridional.append(rec['meridional'])
                         vertical.append(rec['vertical'])
                         self.height[self.imode, index] = (rec['range'] - self.correction) / 1000.
                         numpy.seterr(all='ignore')
                         index = index + 1
                         rx_snr.append(rec['rx_snr'])
                     status = numpy.array(status, dtype='int')
                     zonal = numpy.array(zonal, dtype='float')
                     meridional = numpy.array(meridional, dtype='float')
                     vertical = numpy.array(vertical, dtype='float')
                     rx_snr = numpy.array(rx_snr, dtype='float')
                     rx_snr = rx_snr.reshape((self.nranges, self.nchannels))
                     # FILTERING DATA
                     stvalue = status_value
                     zonal[numpy.where(zonal == -9999.)] = numpy.nan
                     zonal[numpy.where(status != stvalue)] = numpy.nan
                     self.p_zonal[self.count, :, self.imode] = zonal
                     self.one_zonal= self.p_zonal[self.count, :, :]
                     meridional[numpy.where(meridional == -9999.)] = numpy.nan
                     meridional[numpy.where(status != stvalue)] = numpy.nan
                     self.p_meridional[self.count, :, self.imode] = meridional
                     self.one_meridional = self.p_meridional[self.count, :, :]
                     vertical[numpy.where(vertical == -9999.)] = numpy.nan
                     vertical[numpy.where(status != stvalue)] = numpy.nan
                     self.p_vertical[self.count, :, self.imode] = vertical
                     self.one_vertical = self.p_vertical[self.count, :, :]
                     rx_snr[numpy.where(rx_snr == -9999.)] = numpy.nan
                     rx_snr[numpy.where(status != stvalue), :] = numpy.nan
                     for k in range(self.kchan):
                         self.p_snr[self.count, :, k, self.imode] = numpy.power(10, rx_snr[:, k] / 10)
                     self.one_snr = self.p_snr[self.count, :, :, :]
                     if self.nmodes == 2:
                         self.count = self.count + self.imode
                     else:
                         self.count = self.count + 1
                     self.imode +=1
                     self.counter_records = self.counter_records + 1
                     self.zon_ref = self.p_zonal
                     self.ver_ref = self.p_vertical
                     self.mer_ref = self.p_meridional
                     self.snr_ref = self.p_snr
                 def Close (self):
                     '''
                     Closing BLTR rawdata file
                     '''
                     if self.pointer.tell() == self.sizeOfFile:
                         self.pointer.close()
                         return
             class testBLTRWriter(Operation):
                 def __init__(self):
                     Operation.__init__(self)
                     self.dataOut = Parameters()
                     self.path = None
                     self.dataOut = None
                     self.flagIsNewFile=1
                     self.ext = ".hdf5"
                     return
                 def run(self, dataOut, path , modetowrite,**kwargs):
                     if self.flagIsNewFile:
                         flagdata = self.setup(dataOut, path, modetowrite)
                     self.putData()
                     return
                 def setup(self, dataOut, path, modetowrite):
                     '''
                     Recovering data to write in new *.hdf5 file
                     Inputs:
                         modew -- mode to write (1 or 2)
                         path --  destination path
                     '''
                     self.im = modetowrite-1
                     if  self.im!=0 and self.im!=1:
                         raise ValueError, 'Check "modetowrite" value. Must be egual to 1 or 2, "{}" is not valid. '.format(modetowrite)
                     self.dataOut = dataOut
                     self.nmodes = self.dataOut.nmodes
                     self.nchannels = self.dataOut.nchannels
                     self.lat = self.dataOut.lat
                     self.lon = self.dataOut.lon
                     self.hcm = 3
                     self.thisDate = self.dataOut.utctimeInit
                     self.year = self.dataOut.year
                     self.month = self.dataOut.month
                     self.day = self.dataOut.day
                     self.path = path
                     self.flagIsNewFile = 0
                     return 1
                 def setFile(self):
                     '''
                     - Determining the file name for each mode of operation
                         kinst - Kind of Instrument (mnemotic)
                         kindat - Kind of Data (mnemotic)
                     - Creating a cedarObject
                     '''
                     lat_piura = -5.17
                     lat_huancayo = -12.04
                     lat_porcuya = -5.8
                     if '%2.2f' % self.lat == '%2.2f' % lat_piura:
                             self.instMnemonic = 'pbr'
                     elif '%2.2f' % self.lat == '%2.2f' % lat_huancayo:
                             self.instMnemonic = 'hbr'
                     elif '%2.2f' % self.lat == '%2.2f' % lat_porcuya:
                             self.instMnemonic = 'obr'
                     else: raise Warning, "The site of file read doesn't match any site known. Only file from Huancayo, Piura and Porcuya can be processed.\n Check the file "
                     mode = ['_mode1','_mode2']
                     self.hdf5filename = '%s%4.4d%2.2d%2.2d%s%s' % (self.instMnemonic,
                                                         self.year,
                                                         self.month,
                                                         self.day,
                                                         mode[self.im],
                                                         self.ext)
                     self.fullname=os.path.join(self.path,self.hdf5filename)
                     if os.path.isfile(self.fullname) :
                         print "Destination path '%s' already exists. Previous file deleted. " %self.fullname
                         os.remove(self.fullname)
                          # Identify kinst and kindat
                     InstName = self.hdf5filename[0:3]
                     KinstList = [1000, 1001, 1002]
                     KinstId = {'pbr':0, 'hbr':1, 'obr':2}  # pbr:piura, hbr:huancayo, obr:porcuya
                     KindatList = [1600, 1601]  # mode 1, mode 2
                     self.type = KinstId[InstName]
                     self.kinst = KinstList[self.type]
                     self.kindat = KindatList[self.im]
                     try:
                         self.cedarObj = madrigal.cedar.MadrigalCedarFile(self.fullname, True)
                     except ValueError, message:
                         print '[Error]: Impossible to create a cedar object with "madrigal.cedar.MadrigalCedarFile" '
                         return
                     return 1
                 def writeBlock(self):
                     '''
                     - Selecting mode of operation:
                         bltr high resolution mode 1 - Low Atmosphere (0 - 3km) //  bltr high resolution mode 2 - High Atmosphere (0 - 10km)
                         msnr - Average Signal Noise Ratio in dB
                         hcm - 3 km
                     - Filling the cedarObject by a block: each array data entry is assigned a code that defines the parameter to write to the file
                         GDLATR - Reference geod latitude (deg)
                         GDLONR - Reference geographic longitude (deg)
                         GDLAT2 - Geodetic latitude of second inst (deg)
                         GLON2 - Geographic longitude of second inst (deg)
                         GDALT - Geodetic altitude (height) (km)
                         SNL - Log10 (signal to noise ratio)
                         VN1P2 - Neutral wind in direction 1 (eastward) (m/s), ie zonal wind
                         VN2P2 - Neutral wind in direction 2 (northward) (m/s), ie meridional wind
                         EL2 - Ending elevation angle (deg), ie vertical wind
                         Other parameters: /madrigal3/metadata/parcodes.tab
                     '''
                     self.z_zon = self.dataOut.data_output[0,:,:]
                     self.z_mer =self.dataOut.data_output[1,:,:]
                     self.z_ver = self.dataOut.data_output[2,:,:]
                     if self.im == 0:
                         h_select = numpy.where(numpy.bitwise_and(self.dataOut.height[0, :] >= 0., self.dataOut.height[0, :] <= self.hcm, numpy.isfinite(self.dataOut.height[0, :])))
                     else:
                         h_select = numpy.where(numpy.bitwise_and(self.dataOut.height[0, :] >= 0., self.dataOut.height[0, :] < 20, numpy.isfinite(self.dataOut.height[0, :])))
                     ht = h_select[0]
                     self.o_height = self.dataOut.height[self.im, ht]
                     self.o_zon = self.z_zon[ht, self.im]
                     self.o_mer = self.z_mer[ht, self.im]
                     self.o_ver = self.z_ver[ht, self.im]
                     o_snr = self.dataOut.data_SNR[ :, :, self.im]
                     o_snr = o_snr[ht, :]
                     ndiv = numpy.nansum((numpy.isfinite(o_snr)), 1)
                     ndiv = ndiv.astype(float)
                     sel_div = numpy.where(ndiv == 0.)
                     ndiv[sel_div] = numpy.nan
                     if self.nchannels > 1:
                         msnr = numpy.nansum(o_snr, axis=1)
                     else:
                         msnr = o_snr
                     try:
                         self.msnr = 10 * numpy.log10(msnr / ndiv)
                     except ZeroDivisionError:
                         self.msnr = 10 * numpy.log10(msnr /1)
                         print 'Number of division (ndiv) egal to 1 by default. Check SNR'
                     time_t = time.gmtime(self.dataOut.time1)
                     year = time_t.tm_year
                     month = time_t.tm_mon
                     day = time_t.tm_mday
                     hour = time_t.tm_hour
                     minute = time_t.tm_min
                     second = time_t.tm_sec
                     timedate_0 = datetime.datetime(year, month, day, hour, minute, second)
                        # 1d parameters
                     GDLATR = self.lat
                     GDLONR = self.lon
                     GDLAT2 = self.lat
                     GLON2 = self.lon
                      # 2d parameters
                     GDALT = self.o_height
                     SNL = self.msnr
                     VN1P2 = self.o_zon
                     VN2P2 = self.o_mer
                     EL2 = self.o_ver
                     NROW = len(self.o_height)
                     startTime = timedate_0
                     endTime = startTime
                     self.dataRec = madrigal.cedar.MadrigalDataRecord(self.kinst,
                                                                  self.kindat,
                                                                  startTime.year,
                                                                  startTime.month,
                                                                  startTime.day,
                                                                  startTime.hour,
                                                                  startTime.minute,
                                                                  startTime.second,
 ,
                                                                  endTime.year,
                                                                  endTime.month,
                                                                  endTime.day,
                                                                  endTime.hour,
                                                                  endTime.minute,
                                                                  endTime.second,
 ,
                                                                  ('gdlatr', 'gdlonr', 'gdlat2', 'glon2'),
                                                                  ('gdalt', 'snl', 'vn1p2', 'vn2p2', 'el2'),
                                                                  NROW, ind2DList=['gdalt'])
                     # Setting 1d values
                     self.dataRec.set1D('gdlatr', GDLATR)
                     self.dataRec.set1D('gdlonr', GDLONR)
                     self.dataRec.set1D('gdlat2', GDLAT2)
                     self.dataRec.set1D('glon2', GLON2)
                     # Setting 2d values
                     for n in range(self.o_height.shape[0]):
                         self.dataRec.set2D('gdalt', n, GDALT[n])
                         self.dataRec.set2D('snl', n, SNL[n])
                         self.dataRec.set2D('vn1p2', n, VN1P2[n])
                         self.dataRec.set2D('vn2p2', n, VN2P2[n])
                         self.dataRec.set2D('el2', n, EL2[n])
                       # Appending new data record
                     '''
                     [MADRIGAL3]There are two ways to write to a MadrigalCedarFile.  Either this method (write) is called after all the
                     records have been appended to the MadrigalCedarFile, or dump is called after a certain number of records are appended,
                     and then at the end dump is called a final time if there were any records not yet dumped, followed by addArray.
                     '''
                     self.cedarObj.append(self.dataRec)
                     print '    [Writing] records {} (mode {}).'.format(self.dataOut.counter_records,self.im+1)
                     self.cedarObj.dump()
                 def setHeader(self):
                     '''
                     - Creating self.catHeadObj
                     - Adding information catalog
                     - Writing file header
                     '''
                     self.catHeadObj = madrigal.cedar.CatalogHeaderCreator(self.fullname)
                     kindatDesc, comments, analyst, history, principleInvestigator = self._info_BLTR()
                     self.catHeadObj.createCatalog(principleInvestigator="Jarjar",
                                               expPurpose='characterize the atmospheric dynamics in this region where frequently it happens the  El Nino',
                                               sciRemarks="http://madrigal3.haystack.mit.edu/static/CEDARMadrigalHdf5Format.pdf")
                     self.catHeadObj.createHeader(kindatDesc, analyst, comments, history)
                     self.catHeadObj.write()
                     print '[File created] path: %s' % (self.fullname)
                 def putData(self):
                     if self.dataOut.flagNoData:
                         return 0
                     if self.dataOut.counter_records == 1:
                         self.setFile()
                         print '[Writing] Setting new hdf5 file for the mode {}'.format(self.im+1)
                     if self.dataOut.counter_records <= self.dataOut.nrecords:
                         self.writeBlock()
                     if self.dataOut.counter_records == self.dataOut.nrecords:
                         self.cedarObj.addArray()
                         self.setHeader()
                         self.flagIsNewFile = 1
                 def _info_BLTR(self):
                     kindatDesc = '''--This header is for KINDAT = %d''' % self.kindat
                     history = None
                     analyst = '''Jarjar'''
                     principleInvestigator = '''
                         Jarjar
                         Radio Observatorio de Jicamarca
                         Instituto Geofisico del Peru
                         '''
                     if self.type == 1:
                         comments = '''
                         --These data are provided by two Boundary Layer and Tropospheric Radar (BLTR) deployed at two different locations at Peru(GMT-5), one of them at Piura(5.17 S, 80.64W) and another located at Huancayo (12.04 S, 75.32 W).
                         --The purpose of conducting these observations is to measure wind in the differents levels of height, this radar makes measurements the  Zonal(U), Meridional(V) and Vertical(W) wind velocities component in northcoast from Peru. And the main purpose of these mensurations is to characterize the atmospheric dynamics in this region where frequently it happens the  'El Nino Phenomenon'
                         --In Kindat = 1600, contains information of wind velocities component since 0 Km to 3 Km.
                         --In Kindat = 1601, contains information of wind velocities component since 0 Km to 10 Km.
                         --The Huancayo-BLTR is a VHF Profiler Radar System is a 3 channel coherent receiver pulsed radar utilising state-of-the-art software and computing techniques to acquire, decode, and translate signals obtained from partial reflection echoes in the troposphere, lower stratosphere and mesosphere.   It uses an array of three horizontal spaced and vertically directed receiving antennas. The data is recorded thirty seconds, averaged to one minute mean values of  Height, Zonal, Meridional and Vertical wind.
                         --The Huancayo-BLTR was installed in January 2010. This instrument was designed and constructed by Genesis Soft Pty. Ltd. Is constituted by three groups of spaced antennas (distributed) forming an isosceles triangle.
                         Station _______ Geographic Coord ______ Geomagnetic Coord
                         _______________ Latitude _ Longitude __ Latitude _ Longitude
                         Huancayo (HUA) __12.04 S ___ 75.32 W _____ -12.05 ____ 352.85
                         Piura (PIU) _____ 5.17 S ___ 80.64 W ______ 5.18 ____ 350.93
                         WIND OBSERVATIONS
                         --To obtain wind the BLTR uses Spaced Antenna technique (e.g., Briggs 1984). The scatter and reflection it still provided by variations in the refractive index as in the Doppler method(Gage and Basley,1978; Balsley and Gage 1982; Larsen and Rottger 1982), but instead of using the Doppler shift to derive the velocity components, the cross-correlation between signals in an array of three horizontally spaced and vertically directed receiving antennas is used.
                         ......................................................................
                         For more information, consult the following references:
                         - Balsley, B. B., and K. S. Gage., On the use of radars for operational wind profiling, Bull. Amer. Meteor.Soc.,63, 1009-1018, 1982.
                         - Briggs, B. H., The analysis of spaced sensor data by correations techniques, Handbook for MAP, Vol. 13, SCOTEP Secretariat, University of Illinois, Urbana, 166-186, 1984.
                         - Gage, K. S., and B.B. Balsley., Doppler radar probing of the clear atmosphere, Bull. Amer. Meteor.Soc., 59, 1074-1093, 1978.
                         - Larsen, M. F., The Spaced Antenna Technique for Radar Wind Profiling, Journal of Atm. and Ocean. Technology. , Vol.6, 920-937, 1989.
                         - Larsen, M. F., A method for single radar voracity measurements?, Handbook for MAP,SCOSTEP Secretariat, University of the Illinois, Urban, in press, 1989.
                         ......................................................................
                         ACKNOWLEDGEMENTS:
                         --The Piura and Huancayo BLTR are part of the network of instruments operated by the Jicamarca Radio Observatory.
                         --The Jicamarca Radio Observatory is a facility of the Instituto Geofisico del Peru operated with support from the NSF Cooperative Agreement ATM-0432565 through Cornell University
                         ......................................................................
                         Further questions and comments should be addressed to:
                             Radio Observatorio de Jicamarca
                             Instituto Geofisico del Peru
                             Lima, Peru
                             Web URL: http://jro.igp.gob.pe
                         ......................................................................
                 '''
                     return kindatDesc, comments, analyst, history, principleInvestigator

schainpy/model/io/io_bltr_block.pyc binary modified 0 0

NO CONTENT: modified file, binary diff hidden

schainpy/model/io/jroIO_base.pyc binary modified 0 0

NO CONTENT: modified file, binary diff hidden

schainpy/model/io/jroIO_bltr.pyc binary modified 0 0

NO CONTENT: modified file, binary diff hidden

schainpy/model/io/jroIO_mira35c.pyc binary modified 0 0

NO CONTENT: modified file, binary diff hidden

schainpy/model/io/jroIO_param.pyc binary modified 0 0

NO CONTENT: modified file, binary diff hidden

schainpy/model/io/jroIO_spectra.pyc binary modified 0 0

NO CONTENT: modified file, binary diff hidden

schainpy/model/io/jroIO_voltage.pyc binary modified 0 0

NO CONTENT: modified file, binary diff hidden

schainpy/model/proc/__init__.py +1 -1

             '''
             $Author: murco $
             $Id: Processor.py 1 2012-11-12 18:56:07Z murco $
             '''
             from jroproc_voltage import *
             from jroproc_spectra import *
             from jroproc_heispectra import *
             from jroproc_amisr import *
             from jroproc_correlation import *
             from jroproc_parameters import *
             from jroproc_spectra_lags import *
             from jroproc_spectra_acf import *
-            #from jroproc_bltr import *
+            from jroproc_bltr import *

schainpy/model/proc/__init__.pyc binary modified 0 0

NO CONTENT: modified file, binary diff hidden

schainpy/model/proc/jroproc_base.pyc binary modified 0 0

NO CONTENT: modified file, binary diff hidden

schainpy/model/proc/jroproc_bltr.py +2 -2

             '''
             Created on Oct 24, 2016
             @author: roj- LouVD
             '''
             import numpy
             import copy
             import datetime
             import time
             from time import gmtime
             from jroproc_base import ProcessingUnit
             from schainpy.model.data.jrodata import Parameters
             from numpy import transpose
             from matplotlib import cm
             import matplotlib.pyplot as plt
             from matplotlib.mlab import griddata
             class BLTRProcess(ProcessingUnit):
                 isConfig = False
                 '''
                 Processing unit for BLTR rawdata
                 Inputs:
                     self.dataOut.nmodes - Number of operation modes
                     self.dataOut.nchannels - Number of channels
                     self.dataOut.nranges - Number of ranges
                     self.dataOut.data_SNR - SNR array
                     self.dataOut.data_output - Zonal, Vertical and Meridional velocity array
                     self.dataOut.height - Height array (km)
                     self.dataOut.time - Time array (seconds)
                     self.dataOut.fileIndex -Index of the file currently read
                     self.dataOut.lat - Latitude coordinate of BLTR location
                     self.dataOut.doy - Experiment doy (number of the day in the current year)
                     self.dataOut.month - Experiment month
                     self.dataOut.day - Experiment day
                     self.dataOut.year - Experiment year
                 '''
-                def __init__(self):
+                def __init__(self, **kwargs):
                     '''
                     Inputs: None
                     '''
-                    ProcessingUnit.__init__(self)
+                    ProcessingUnit.__init__(self, **kwargs)
                     self.dataOut = Parameters()
                     # Filters
                     snr_val = None
                     value = None
                     svalue2 = None
                     method = None
                     factor = None
                     filter = None
                     npoints = None
                     status_value = None
                     width = None
                     self.flagfirstmode = 0
                 def run (self):
                     if self.dataIn.type == "Parameters":
                         self.dataOut.copy(self.dataIn)
                 def TimeSelect(self):
                     '''
                     Selecting the time array according to the day of the experiment with a duration of 24 hours
                     '''
                     k1 = datetime.datetime(self.dataOut.year, self.dataOut.month, self.dataOut.day) - datetime.timedelta(hours=5)
                     k2 = datetime.datetime(self.dataOut.year, self.dataOut.month, self.dataOut.day) + datetime.timedelta(hours=25) - datetime.timedelta(hours=5)
                     limit_sec1 = time.mktime(k1.timetuple())
                     limit_sec2 = time.mktime(k2.timetuple())
                     valid_data = 0
                     doy = self.dataOut.doy
                     t1 = numpy.where(self.dataOut.time[0, :] >= limit_sec1)
                     t2 = numpy.where(self.dataOut.time[0, :] < limit_sec2)
                     time_select = []
                     for val_sec in t1[0]:
                         if val_sec in t2[0]:
                             time_select.append(val_sec)
                     time_select = numpy.array(time_select, dtype='int')
                     valid_data = valid_data + len(time_select)
                     if len(time_select) > 0:
                         self.f_timesec = self.dataOut.time[:, time_select]
                         snr = self.dataOut.data_SNR[time_select, :, :, :]
                         zon = self.dataOut.data_output[0][time_select, :, :]
                         mer = self.dataOut.data_output[1][time_select, :, :]
                         ver = self.dataOut.data_output[2][time_select, :, :]
                     if valid_data > 0:
                         self.timesec1 = self.f_timesec[0, :]
                         self.f_height = self.dataOut.height
                         self.f_zon = zon
                         self.f_mer = mer
                         self.f_ver = ver
                         self.f_snr = snr
                         self.f_timedate = []
                         self.f_time = []
                         for valuet in self.timesec1:
                             time_t = time.gmtime(valuet)
                             year = time_t.tm_year
                             month = time_t.tm_mon
                             day = time_t.tm_mday
                             hour = time_t.tm_hour
                             minute = time_t.tm_min
                             second = time_t.tm_sec
                             f_timedate_0 = datetime.datetime(year, month, day, hour, minute, second)
                             self.f_timedate.append(f_timedate_0)
                         return self.f_timedate, self.f_timesec, self.f_height, self.f_zon, self.f_mer, self.f_ver, self.f_snr
                     else:
                         self.f_timesec = None
                         self.f_timedate = None
                         self.f_height = None
                         self.f_zon = None
                         self.f_mer = None
                         self.f_ver = None
                         self.f_snr = None
                         print 'Invalid time'
                         return self.f_timedate, self.f_height, self.f_zon, self.f_mer, self.f_ver, self.f_snr
                 def SnrFilter(self, snr_val,modetofilter):
                     '''
                     Inputs: snr_val - Threshold value
                     '''
                     if modetofilter!=2 and modetofilter!=1 :
                         raise ValueError,'Mode to filter should be "1" or "2". {} is not valid, check "Modetofilter" value.'.format(modetofilter)
                     m = modetofilter-1
                     print '    SNR filter [mode {}]: SNR <= {}: data_output = NA'.format(modetofilter,snr_val)
                     for k in range(self.dataOut.nchannels):
                         for r in range(self.dataOut.nranges):
                             if self.dataOut.data_SNR[r,k,m] <= snr_val:
                                 self.dataOut.data_output[2][r,m] = numpy.nan
                                 self.dataOut.data_output[1][r,m] = numpy.nan
                                 self.dataOut.data_output[0][r,m] = numpy.nan
                 def OutliersFilter(self,modetofilter,svalue,svalue2,method,factor,filter,npoints):
                     '''
                     Inputs:
                         svalue    -  string to select array velocity
                         svalue2    -  string to choose axis filtering
                         method    - 0 for SMOOTH or 1 for MEDIAN
                         factor    - number used to set threshold
                         filter    - 1 for data filtering using the standard deviation criteria else 0
                         npoints    - number of points for mask filter
                     '''
                     if modetofilter!=2 and modetofilter!=1 :
                         raise ValueError,'Mode to filter should be "1" or "2". {} is not valid, check "Modetofilter" value.'.format(modetofilter)
                     m = modetofilter-1
                     print '    Outliers Filter [mode {}]: {} {} / threshold = {}'.format(modetofilter,svalue,svalue,factor)
                     npoints = 9
                     novalid = 0.1
                     if svalue == 'zonal':
                         value = self.dataOut.data_output[0]
                     elif svalue == 'meridional':
                         value = self.dataOut.data_output[1]
                     elif svalue == 'vertical':
                         value = self.dataOut.data_output[2]
                     else:
                         print 'value is not defined'
                         return
                     if svalue2 == 'inTime':
                         yaxis = self.dataOut.height
                         xaxis = numpy.array([[self.dataOut.time1],[self.dataOut.time1]])
                     elif svalue2 == 'inHeight':
                         yaxis = numpy.array([[self.dataOut.time1],[self.dataOut.time1]])
                         xaxis = self.dataOut.height
                     else:
                         print 'svalue2 is required, either inHeight or inTime'
                         return
                     output_array = value
                     value_temp = value[:,m]
                     error = numpy.zeros(len(self.dataOut.time[m,:]))
                     if svalue2 == 'inHeight':
                         value_temp = numpy.transpose(value_temp)
                         error = numpy.zeros(len(self.dataOut.height))
                     htemp = yaxis[m,:]
                     std = value_temp
                     for h in range(len(htemp)):
                         if filter: #standard deviation filtering
                             std[h] = numpy.std(value_temp[h],ddof = npoints)
                             value_temp[numpy.where(std[h] > 5),h] = numpy.nan
                             error[numpy.where(std[h] > 5)] = error[numpy.where(std[h] > 5)] + 1
                         nvalues_valid = len(numpy.where(numpy.isfinite(value_temp[h]))[0])
                         minvalid = novalid*len(xaxis[m,:])
                         if minvalid <= npoints:
                             minvalid = npoints
                         #only if valid values greater than the minimum required (10%)
                         if nvalues_valid > minvalid:
                             if method == 0:
                                 #SMOOTH
                                 w = value_temp[h] - self.Smooth(input=value_temp[h], width=npoints, edge_truncate=1)
                             if method == 1:
                                 #MEDIAN
                                 w = value_temp[h] - self.Median(input=value_temp[h], width = npoints)
                             dw = numpy.std(w[numpy.where(numpy.isfinite(w))],ddof = 1)
                             threshold = dw*factor
                             value_temp[numpy.where(w > threshold),h] = numpy.nan
                             value_temp[numpy.where(w < -1*threshold),h] = numpy.nan
                     #At the end
                     if svalue2 == 'inHeight':
                         value_temp = numpy.transpose(value_temp)
                     output_array[:,m] = value_temp
                     if svalue == 'zonal':
                         self.dataOut.data_output[0] = output_array
                     elif svalue == 'meridional':
                         self.dataOut.data_output[1] = output_array
                     elif svalue == 'vertical':
                         self.dataOut.data_output[2] = output_array
                     return self.dataOut.data_output
                 def Median(self,input,width):
                     '''
                     Inputs:
                         input - Velocity array
                         width - Number of points for mask filter
                     '''
                     if numpy.mod(width,2) == 1:
                         pc = int((width - 1) / 2)
                     cont = 0
                     output = []
                     for i in range(len(input)):
                         if i >= pc and i < len(input) - pc:
                             new2 = input[i-pc:i+pc+1]
                             temp = numpy.where(numpy.isfinite(new2))
                             new = new2[temp]
                             value = numpy.median(new)
                             output.append(value)
                     output = numpy.array(output)
                     output = numpy.hstack((input[0:pc],output))
                     output = numpy.hstack((output,input[-pc:len(input)]))
                     return output
                 def Smooth(self,input,width,edge_truncate = None):
                     '''
                     Inputs:
                         input - Velocity array
                         width - Number of points for mask filter
                         edge_truncate - 1 for truncate the convolution product else
                     '''
                     if numpy.mod(width,2) == 0:
                         real_width = width + 1
                         nzeros = width / 2
                     else:
                         real_width = width
                         nzeros = (width - 1) / 2
                     half_width = int(real_width)/2
                     length = len(input)
                     gate = numpy.ones(real_width,dtype='float')
                     norm_of_gate = numpy.sum(gate)
                     nan_process = 0
                     nan_id = numpy.where(numpy.isnan(input))
                     if len(nan_id[0]) > 0:
                         nan_process = 1
                         pb = numpy.zeros(len(input))
                         pb[nan_id] = 1.
                         input[nan_id] = 0.
                     if edge_truncate == True:
                         output = numpy.convolve(input/norm_of_gate,gate,mode='same')
                     elif edge_truncate == False or edge_truncate == None:
                         output = numpy.convolve(input/norm_of_gate,gate,mode='valid')
                         output = numpy.hstack((input[0:half_width],output))
                         output = numpy.hstack((output,input[len(input)-half_width:len(input)]))
                     if nan_process:
                         pb = numpy.convolve(pb/norm_of_gate,gate,mode='valid')
                         pb = numpy.hstack((numpy.zeros(half_width),pb))
                         pb = numpy.hstack((pb,numpy.zeros(half_width)))
                         output[numpy.where(pb > 0.9999)] = numpy.nan
                         input[nan_id] = numpy.nan
                     return output
                 def Average(self,aver=0,nhaver=1):
                     '''
                     Inputs:
                         aver - Indicates the time period over which is averaged or consensus data
                         nhaver - Indicates the decimation factor in heights
                     '''
                     nhpoints = 48
                     lat_piura = -5.17
                     lat_huancayo = -12.04
                     lat_porcuya = -5.8
                     if '%2.2f'%self.dataOut.lat == '%2.2f'%lat_piura:
                         hcm = 3.
                         if self.dataOut.year == 2003 :
                             if self.dataOut.doy >= 25 and self.dataOut.doy < 64:
                                 nhpoints = 12
                     elif '%2.2f'%self.dataOut.lat == '%2.2f'%lat_huancayo:
                         hcm = 3.
                         if self.dataOut.year == 2003 :
                             if self.dataOut.doy >= 25 and self.dataOut.doy < 64:
                                 nhpoints = 12
                     elif '%2.2f'%self.dataOut.lat == '%2.2f'%lat_porcuya:
                         hcm = 5.#2
                     pdata = 0.2
                     taver = [1,2,3,4,6,8,12,24]
                     t0 = 0
                     tf = 24
                     ntime =(tf-t0)/taver[aver]
                     ti = numpy.arange(ntime)
                     tf = numpy.arange(ntime) + taver[aver]
                     old_height = self.dataOut.heightList
                     if nhaver > 1:
                         num_hei = len(self.dataOut.heightList)/nhaver/self.dataOut.nmodes
                         deltha = 0.05*nhaver
                         minhvalid = pdata*nhaver
                         for im in range(self.dataOut.nmodes):
                             new_height = numpy.arange(num_hei)*deltha + self.dataOut.height[im,0] + deltha/2.
                     data_fHeigths_List = []
                     data_fZonal_List = []
                     data_fMeridional_List = []
                     data_fVertical_List = []
                     startDTList = []
                     for i in range(ntime):
                         height = old_height
                         start = datetime.datetime(self.dataOut.year,self.dataOut.month,self.dataOut.day) + datetime.timedelta(hours = int(ti[i])) - datetime.timedelta(hours = 5)
                         stop = datetime.datetime(self.dataOut.year,self.dataOut.month,self.dataOut.day) + datetime.timedelta(hours = int(tf[i])) - datetime.timedelta(hours = 5)
                         limit_sec1 = time.mktime(start.timetuple())
                         limit_sec2 = time.mktime(stop.timetuple())
                         t1 = numpy.where(self.f_timesec >= limit_sec1)
                         t2 = numpy.where(self.f_timesec < limit_sec2)
                         time_select = []
                         for val_sec in t1[0]:
                             if val_sec in t2[0]:
                                 time_select.append(val_sec)
                         time_select = numpy.array(time_select,dtype = 'int')
                         minvalid = numpy.ceil(pdata*nhpoints)
                         zon_aver = numpy.zeros([self.dataOut.nranges,self.dataOut.nmodes],dtype='f4') + numpy.nan
                         mer_aver = numpy.zeros([self.dataOut.nranges,self.dataOut.nmodes],dtype='f4') + numpy.nan
                         ver_aver = numpy.zeros([self.dataOut.nranges,self.dataOut.nmodes],dtype='f4') + numpy.nan
                         if nhaver > 1:
                             new_zon_aver = numpy.zeros([num_hei,self.dataOut.nmodes],dtype='f4') + numpy.nan
                             new_mer_aver = numpy.zeros([num_hei,self.dataOut.nmodes],dtype='f4') + numpy.nan
                             new_ver_aver = numpy.zeros([num_hei,self.dataOut.nmodes],dtype='f4') + numpy.nan
                         if len(time_select) > minvalid:
                             time_average = self.f_timesec[time_select]
                             for im in range(self.dataOut.nmodes):
                                 for ih in range(self.dataOut.nranges):
                                     if numpy.sum(numpy.isfinite(self.f_zon[time_select,ih,im])) >= minvalid:
                                         zon_aver[ih,im] = numpy.nansum(self.f_zon[time_select,ih,im]) / numpy.sum(numpy.isfinite(self.f_zon[time_select,ih,im]))
                                     if numpy.sum(numpy.isfinite(self.f_mer[time_select,ih,im])) >= minvalid:
                                         mer_aver[ih,im] = numpy.nansum(self.f_mer[time_select,ih,im]) / numpy.sum(numpy.isfinite(self.f_mer[time_select,ih,im]))
                                     if numpy.sum(numpy.isfinite(self.f_ver[time_select,ih,im])) >= minvalid:
                                         ver_aver[ih,im] = numpy.nansum(self.f_ver[time_select,ih,im]) / numpy.sum(numpy.isfinite(self.f_ver[time_select,ih,im]))
                                 if nhaver > 1:
                                     for ih in range(num_hei):
                                         hvalid = numpy.arange(nhaver) + nhaver*ih
                                         if numpy.sum(numpy.isfinite(zon_aver[hvalid,im])) >= minvalid:
                                             new_zon_aver[ih,im] = numpy.nansum(zon_aver[hvalid,im]) /  numpy.sum(numpy.isfinite(zon_aver[hvalid,im]))
                                         if numpy.sum(numpy.isfinite(mer_aver[hvalid,im])) >= minvalid:
                                             new_mer_aver[ih,im] = numpy.nansum(mer_aver[hvalid,im]) /  numpy.sum(numpy.isfinite(mer_aver[hvalid,im]))
                                         if numpy.sum(numpy.isfinite(ver_aver[hvalid,im])) >= minvalid:
                                             new_ver_aver[ih,im] = numpy.nansum(ver_aver[hvalid,im]) /  numpy.sum(numpy.isfinite(ver_aver[hvalid,im]))
                             if nhaver > 1:
                                 zon_aver = new_zon_aver
                                 mer_aver = new_mer_aver
                                 ver_aver = new_ver_aver
                                 height = new_height
                             tstart = time_average[0]
                             tend = time_average[-1]
                             startTime = time.gmtime(tstart)
                             year = startTime.tm_year
                             month = startTime.tm_mon
                             day = startTime.tm_mday
                             hour = startTime.tm_hour
                             minute = startTime.tm_min
                             second = startTime.tm_sec
                             startDTList.append(datetime.datetime(year,month,day,hour,minute,second))
                             o_height = numpy.array([])
                             o_zon_aver = numpy.array([])
                             o_mer_aver = numpy.array([])
                             o_ver_aver = numpy.array([])
                             if self.dataOut.nmodes > 1:
                                 for im in range(self.dataOut.nmodes):
                                     if im == 0:
                                         h_select = numpy.where(numpy.bitwise_and(height[0,:] >=0,height[0,:] <= hcm,numpy.isfinite(height[0,:])))
                                     else:
                                         h_select = numpy.where(numpy.bitwise_and(height[1,:] > hcm,height[1,:] < 20,numpy.isfinite(height[1,:])))
                                     ht = h_select[0]
                                     o_height = numpy.hstack((o_height,height[im,ht]))
                                     o_zon_aver = numpy.hstack((o_zon_aver,zon_aver[ht,im]))
                                     o_mer_aver = numpy.hstack((o_mer_aver,mer_aver[ht,im]))
                                     o_ver_aver = numpy.hstack((o_ver_aver,ver_aver[ht,im]))
                                 data_fHeigths_List.append(o_height)
                                 data_fZonal_List.append(o_zon_aver)
                                 data_fMeridional_List.append(o_mer_aver)
                                 data_fVertical_List.append(o_ver_aver)
                             else:
                                 h_select = numpy.where(numpy.bitwise_and(height[0,:] <= hcm,numpy.isfinite(height[0,:])))
                                 ht = h_select[0]
                                 o_height = numpy.hstack((o_height,height[im,ht]))
                                 o_zon_aver = numpy.hstack((o_zon_aver,zon_aver[ht,im]))
                                 o_mer_aver = numpy.hstack((o_mer_aver,mer_aver[ht,im]))
                                 o_ver_aver = numpy.hstack((o_ver_aver,ver_aver[ht,im]))
                                 data_fHeigths_List.append(o_height)
                                 data_fZonal_List.append(o_zon_aver)
                                 data_fMeridional_List.append(o_mer_aver)
                                 data_fVertical_List.append(o_ver_aver)
                     return startDTList, data_fHeigths_List, data_fZonal_List, data_fMeridional_List, data_fVertical_List
                 def prePlot(self,modeselect=None):
                     '''
                   Inputs:
                       self.dataOut.data_output - Zonal, Meridional and Vertical velocity array
                       self.dataOut.height - height array
                       self.dataOut.time  - Time array (seconds)
                       self.dataOut.data_SNR - SNR array
                       '''
                     m = modeselect -1
                     print '    [Plotting mode {}]'.format(modeselect)
                     if not (m ==1 or m==0):
                         raise IndexError("'Mode' must be egual to : 1 or 2")
             #
                     if self.flagfirstmode==0:
                         #copy of the data
                         self.data_output_copy = self.dataOut.data_output.copy()
                         self.data_height_copy = self.dataOut.height.copy()
                         self.data_time_copy = self.dataOut.time.copy()
                         self.data_SNR_copy = self.dataOut.data_SNR.copy()
                         self.flagfirstmode = 1
                     else:
                         self.dataOut.data_output = self.data_output_copy
                         self.dataOut.height = self.data_height_copy
                         self.dataOut.time = self.data_time_copy
                         self.dataOut.data_SNR = self.data_SNR_copy
                         self.flagfirstmode = 0
                     #select data for mode m
                     #self.dataOut.data_output = self.dataOut.data_output[:,:,m]
                     self.dataOut.heightList = self.dataOut.height[0,:]
                     data_SNR = self.dataOut.data_SNR[:,:,m]
                     self.dataOut.data_SNR= transpose(data_SNR)
                     if m==1 and self.dataOut.counter_records%2==0:
                         print '*********'
                         print 'MODO 2'
                         #print 'Zonal', self.dataOut.data_output[0]
                         #print 'Meridional', self.dataOut.data_output[1]
                         #print 'Vertical', self.dataOut.data_output[2]
                         print '*********'
                         Vx=self.dataOut.data_output[0,:,m]
                         Vy=self.dataOut.data_output[1,:,m]
                         Vmag=numpy.sqrt(Vx**2+Vy**2)
                         Vang=numpy.arctan2(Vy,Vx)
                         #print 'Vmag', Vmag
                         #print 'Vang', Vang
                         self.dataOut.data_output[0,:,m]=Vmag
                         self.dataOut.data_output[1,:,m]=Vang
                         prin= self.dataOut.data_output[0,:,m][~numpy.isnan(self.dataOut.data_output[0,:,m])]
                         print ' '
                         print 'VmagAverage',numpy.mean(prin)
                         print ' '
                     self.dataOut.data_output = self.dataOut.data_output[:,:,m]
   No newline at end of file

schainpy/model/proc/jroproc_bltr.pyc binary modified 0 0

NO CONTENT: modified file, binary diff hidden

schainpy/model/proc/jroproc_correlation.pyc binary modified 0 0

NO CONTENT: modified file, binary diff hidden

schainpy/model/proc/jroproc_parameters.py +66 -35

             import numpy
             import math
             from scipy import optimize, interpolate, signal, stats, ndimage
             import scipy
             import re
             import datetime
             import copy
             import sys
             import importlib
             import itertools
             from multiprocessing import Pool, TimeoutError
             from multiprocessing.pool import ThreadPool
             import copy_reg
             import cPickle
             import types
             from functools import partial
             import time
             #from sklearn.cluster import KMeans
             import matplotlib.pyplot as plt
             from scipy.optimize import fmin_l_bfgs_b #optimize with bounds on state papameters
             from jroproc_base import ProcessingUnit, Operation
             from schainpy.model.data.jrodata import Parameters, hildebrand_sekhon
             from scipy import asarray as ar,exp
             from scipy.optimize import curve_fit
             import warnings
             from numpy import NaN
+            from scipy.optimize.optimize import OptimizeWarning
             warnings.filterwarnings('ignore')
             SPEED_OF_LIGHT = 299792458
             '''solving pickling issue'''
             def _pickle_method(method):
                 func_name = method.im_func.__name__
                 obj = method.im_self
                 cls = method.im_class
                 return _unpickle_method, (func_name, obj, cls)
             def _unpickle_method(func_name, obj, cls):
                 for cls in cls.mro():
                     try:
                         func = cls.__dict__[func_name]
                     except KeyError:
                         pass
                     else:
                         break
                 return func.__get__(obj, cls)
             class ParametersProc(ProcessingUnit):
                 nSeconds = None
                 def __init__(self):
                     ProcessingUnit.__init__(self)
             #         self.objectDict = {}
                     self.buffer = None
                     self.firstdatatime = None
                     self.profIndex = 0
                     self.dataOut = Parameters()
                 def __updateObjFromInput(self):
                     self.dataOut.inputUnit = self.dataIn.type
                     self.dataOut.timeZone = self.dataIn.timeZone
                     self.dataOut.dstFlag = self.dataIn.dstFlag
                     self.dataOut.errorCount = self.dataIn.errorCount
                     self.dataOut.useLocalTime = self.dataIn.useLocalTime
                     self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
                     self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
                     self.dataOut.channelList = self.dataIn.channelList
                     self.dataOut.heightList = self.dataIn.heightList
                     self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
             #         self.dataOut.nHeights = self.dataIn.nHeights
             #         self.dataOut.nChannels = self.dataIn.nChannels
                     self.dataOut.nBaud = self.dataIn.nBaud
                     self.dataOut.nCode = self.dataIn.nCode
                     self.dataOut.code = self.dataIn.code
             #        self.dataOut.nProfiles = self.dataOut.nFFTPoints
                     self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
             #         self.dataOut.utctime = self.firstdatatime
                     self.dataOut.utctime = self.dataIn.utctime
                     self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
                     self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
                     self.dataOut.nCohInt = self.dataIn.nCohInt
             #        self.dataOut.nIncohInt = 1
                     self.dataOut.ippSeconds = self.dataIn.ippSeconds
             #        self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
                     self.dataOut.timeInterval1 = self.dataIn.timeInterval
                     self.dataOut.heightList = self.dataIn.getHeiRange()
                     self.dataOut.frequency = self.dataIn.frequency
                     self.dataOut.noise = self.dataIn.noise
                 def run(self):
                     #----------------------    Voltage Data    ---------------------------
                     if self.dataIn.type == "Voltage":
                         self.__updateObjFromInput()
                         self.dataOut.data_pre = self.dataIn.data.copy()
                         self.dataOut.flagNoData = False
                         self.dataOut.utctimeInit = self.dataIn.utctime
                         self.dataOut.paramInterval = self.dataIn.nProfiles*self.dataIn.nCohInt*self.dataIn.ippSeconds
                         return
                     #----------------------    Spectra Data    ---------------------------
                     if self.dataIn.type == "Spectra":
                         self.dataOut.data_pre = (self.dataIn.data_spc,self.dataIn.data_cspc)
                         print 'self.dataIn.data_spc', self.dataIn.data_spc.shape
                         self.dataOut.abscissaList = self.dataIn.getVelRange(1)
                         self.dataOut.spc_noise = self.dataIn.getNoise()
                         self.dataOut.spc_range = (self.dataIn.getFreqRange(1)/1000. , self.dataIn.getAcfRange(1) , self.dataIn.getVelRange(1) )
                         self.dataOut.normFactor = self.dataIn.normFactor
                         #self.dataOut.outputInterval = self.dataIn.outputInterval
                         self.dataOut.groupList = self.dataIn.pairsList
                         self.dataOut.flagNoData = False
                         #print 'datain chandist ',self.dataIn.ChanDist
                         if hasattr(self.dataIn, 'ChanDist'): #Distances of receiver channels
                             self.dataOut.ChanDist = self.dataIn.ChanDist
                         else: self.dataOut.ChanDist = None
                         print 'datain chandist ',self.dataOut.ChanDist
                         if hasattr(self.dataIn, 'VelRange'): #Velocities range
                             self.dataOut.VelRange = self.dataIn.VelRange
                         else: self.dataOut.VelRange = None
                         if hasattr(self.dataIn, 'RadarConst'): #Radar Constant
                             self.dataOut.RadarConst = self.dataIn.RadarConst
                         if hasattr(self.dataIn, 'NPW'): #NPW
                             self.dataOut.NPW = self.dataIn.NPW
                         if hasattr(self.dataIn, 'COFA'): #COFA
                             self.dataOut.COFA = self.dataIn.COFA
                     #----------------------    Correlation Data    ---------------------------
                     if self.dataIn.type == "Correlation":
                         acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.dataIn.splitFunctions()
                         self.dataOut.data_pre = (self.dataIn.data_cf[acf_ind,:], self.dataIn.data_cf[ccf_ind,:,:])
                         self.dataOut.normFactor = (self.dataIn.normFactor[acf_ind,:], self.dataIn.normFactor[ccf_ind,:])
                         self.dataOut.groupList = (acf_pairs, ccf_pairs)
                         self.dataOut.abscissaList = self.dataIn.lagRange
                         self.dataOut.noise = self.dataIn.noise
                         self.dataOut.data_SNR = self.dataIn.SNR
                         self.dataOut.flagNoData = False
                         self.dataOut.nAvg = self.dataIn.nAvg
                     #----------------------    Parameters Data    ---------------------------
                     if self.dataIn.type == "Parameters":
                         self.dataOut.copy(self.dataIn)
                         self.dataOut.flagNoData = False
                         return True
                     self.__updateObjFromInput()
                     self.dataOut.utctimeInit = self.dataIn.utctime
                     self.dataOut.paramInterval = self.dataIn.timeInterval
                     return
             def target(tups):
                 obj, args = tups
                 #print 'TARGETTT', obj, args
                 return obj.FitGau(args)
             class GaussianFit(Operation):
                 '''
                     Function that fit of one and two generalized gaussians (gg) based
                     on the PSD shape across an "power band" identified from a cumsum of
                     the measured spectrum - noise.
                     Input:
                         self.dataOut.data_pre    :    SelfSpectra
                     Output:
                         self.dataOut.GauSPC :    SPC_ch1, SPC_ch2
                 '''
                 def __init__(self, **kwargs):
                     Operation.__init__(self, **kwargs)
                     self.i=0
                 def run(self, dataOut, num_intg=7, pnoise=1., vel_arr=None, SNRlimit=-9): #num_intg: Incoherent integrations, pnoise: Noise, vel_arr: range of velocities, similar to the ftt points
                     """This routine will find a couple of generalized Gaussians to a power spectrum
                     input: spc
                     output:
                         Amplitude0,shift0,width0,p0,Amplitude1,shift1,width1,p1,noise
                     """
                     self.spc = dataOut.data_pre[0].copy()
                     print 'SelfSpectra Shape', numpy.asarray(self.spc).shape
                     #plt.figure(50)
                     #plt.subplot(121)
                     #plt.plot(self.spc,'k',label='spc(66)')
                     #plt.plot(xFrec,ySamples[1],'g',label='Ch1')
                     #plt.plot(xFrec,ySamples[2],'r',label='Ch2')
                     #plt.plot(xFrec,FitGauss,'yo:',label='fit')
                     #plt.legend()
                     #plt.title('DATOS A ALTURA DE 7500 METROS')
                     #plt.show()
                     self.Num_Hei = self.spc.shape[2]
                     #self.Num_Bin = len(self.spc)
                     self.Num_Bin = self.spc.shape[1]
                     self.Num_Chn = self.spc.shape[0]
                     Vrange = dataOut.abscissaList
                     #print 'self.spc2', numpy.asarray(self.spc).shape
                     GauSPC = numpy.empty([2,self.Num_Bin,self.Num_Hei])
                     SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei])
                     SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei])
                     SPC_ch1[:] = numpy.NaN
                     SPC_ch2[:] = numpy.NaN
                     start_time = time.time()
                     noise_ = dataOut.spc_noise[0].copy()
                     pool = Pool(processes=self.Num_Chn)
                     args = [(Vrange, Ch, pnoise, noise_, num_intg, SNRlimit) for Ch in range(self.Num_Chn)]
                     objs = [self for __ in range(self.Num_Chn)]
                     attrs = zip(objs, args)
                     gauSPC = pool.map(target, attrs)
                     dataOut.GauSPC = numpy.asarray(gauSPC)
             #         ret = []
             #         for n in range(self.Num_Chn):
             #             self.FitGau(args[n])
             #         dataOut.GauSPC = ret
             #         for ch in range(self.Num_Chn):
             #
             #             for ht in range(self.Num_Hei):
             #                 #print (numpy.asarray(self.spc).shape)
             #                 spc =  numpy.asarray(self.spc)[ch,:,ht]
             #
             #                 #############################################
             #                 # normalizing spc and noise
             #                 # This part differs from gg1
             #                 spc_norm_max = max(spc)
             #                 spc = spc / spc_norm_max
             #                 pnoise = pnoise / spc_norm_max
             #                 #############################################
             #
             #                 if abs(vel_arr[0])<15.0: # this switch is for spectra collected with different length IPP's
             #                     fatspectra=1.0
             #                 else:
             #                     fatspectra=0.5
             #
             #                 wnoise =  noise_ / spc_norm_max
             #                 #print 'wnoise', noise_, dataOut.spc_noise[0], wnoise
             #                 #wnoise,stdv,i_max,index =enoise(spc,num_intg) #noise estimate using Hildebrand Sekhon, only wnoise is used
             #                 #if wnoise>1.1*pnoise: # to be tested later
             #                 #    wnoise=pnoise
             #                 noisebl=wnoise*0.9; noisebh=wnoise*1.1
             #                 spc=spc-wnoise
             #
             #                 minx=numpy.argmin(spc)
             #                 spcs=numpy.roll(spc,-minx)
             #                 cum=numpy.cumsum(spcs)
             #                 tot_noise=wnoise * self.Num_Bin  #64;
             #                 #tot_signal=sum(cum[-5:])/5.; ''' How does this line work? '''
             #                 #snr=tot_signal/tot_noise
             #                 #snr=cum[-1]/tot_noise
             #
             #                 #print 'spc' , spcs[5:8] , 'tot_noise', tot_noise
             #
             #                 snr = sum(spcs)/tot_noise
             #                 snrdB=10.*numpy.log10(snr)
             #
             #                 #if snrdB < -9 :
             #                 #    snrdB = numpy.NaN
             #                 #    continue
             #
             #                 #print 'snr',snrdB # , sum(spcs) , tot_noise
             #
             #
             #                 #if snrdB<-18 or numpy.isnan(snrdB) or num_intg<4:
             #                 #    return [None,]*4,[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None
             #
             #                 cummax=max(cum); epsi=0.08*fatspectra # cumsum to narrow down the energy region
             #                 cumlo=cummax*epsi;
             #                 cumhi=cummax*(1-epsi)
             #                 powerindex=numpy.array(numpy.where(numpy.logical_and(cum>cumlo, cum<cumhi))[0])
             #
             #                 #if len(powerindex)==1:
             #                     ##return [numpy.mod(powerindex[0]+minx,64),None,None,None,],[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None
             #                     #return [numpy.mod(powerindex[0]+minx, self.Num_Bin ),None,None,None,],[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None
             #                 #elif len(powerindex)<4*fatspectra:
             #                     #return [None,]*4,[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None
             #
             #                 if len(powerindex) < 1:# case for powerindex 0
             #                     continue
             #                 powerlo=powerindex[0]
             #                 powerhi=powerindex[-1]
             #                 powerwidth=powerhi-powerlo
             #
             #                 firstpeak=powerlo+powerwidth/10.# first gaussian energy location
             #                 secondpeak=powerhi-powerwidth/10.#second gaussian energy location
             #                 midpeak=(firstpeak+secondpeak)/2.
             #                 firstamp=spcs[int(firstpeak)]
             #                 secondamp=spcs[int(secondpeak)]
             #                 midamp=spcs[int(midpeak)]
             #                 #x=numpy.spc.shape[1]
             #
             #                 #x=numpy.arange(64)
             #                 x=numpy.arange( self.Num_Bin )
             #                 y_data=spc+wnoise
             #
             #                 # single gaussian
             #                 #shift0=numpy.mod(midpeak+minx,64)
             #                 shift0=numpy.mod(midpeak+minx, self.Num_Bin )
             #                 width0=powerwidth/4.#Initialization entire power of spectrum divided by 4
             #                 power0=2.
             #                 amplitude0=midamp
             #                 state0=[shift0,width0,amplitude0,power0,wnoise]
             #                 #bnds=((0,63),(1,powerwidth),(0,None),(0.5,3.),(noisebl,noisebh))
             #                 bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth),(0,None),(0.5,3.),(noisebl,noisebh))
             #                 #bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth),(0,None),(0.5,3.),(0.1,0.5))
             #                 # bnds = range of fft, power width, amplitude, power, noise
             #                 lsq1=fmin_l_bfgs_b(self.misfit1,state0,args=(y_data,x,num_intg),bounds=bnds,approx_grad=True)
             #
             #                 chiSq1=lsq1[1];
             #                 jack1= self.y_jacobian1(x,lsq1[0])
             #
             #
             #                 try:
             #                     sigmas1=numpy.sqrt(chiSq1*numpy.diag(numpy.linalg.inv(numpy.dot(jack1.T,jack1))))
             #                 except:
             #                     std1=32.; sigmas1=numpy.ones(5)
             #                 else:
             #                     std1=sigmas1[0]
             #
             #
             #                 if fatspectra<1.0 and powerwidth<4:
             #                         choice=0
             #                         Amplitude0=lsq1[0][2]
             #                         shift0=lsq1[0][0]
             #                         width0=lsq1[0][1]
             #                         p0=lsq1[0][3]
             #                         Amplitude1=0.
             #                         shift1=0.
             #                         width1=0.
             #                         p1=0.
             #                         noise=lsq1[0][4]
             #                         #return (numpy.array([shift0,width0,Amplitude0,p0]),
             #                         #        numpy.array([shift1,width1,Amplitude1,p1]),noise,snrdB,chiSq1,6.,sigmas1,[None,]*9,choice)
             #
             #                 # two gaussians
             #                 #shift0=numpy.mod(firstpeak+minx,64); shift1=numpy.mod(secondpeak+minx,64)
             #                 shift0=numpy.mod(firstpeak+minx, self.Num_Bin );
             #                 shift1=numpy.mod(secondpeak+minx, self.Num_Bin )
             #                 width0=powerwidth/6.;
             #                 width1=width0
             #                 power0=2.;
             #                 power1=power0
             #                 amplitude0=firstamp;
             #                 amplitude1=secondamp
             #                 state0=[shift0,width0,amplitude0,power0,shift1,width1,amplitude1,power1,wnoise]
             #                 #bnds=((0,63),(1,powerwidth/2.),(0,None),(0.5,3.),(0,63),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
             #                 bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth/2.),(0,None),(0.5,3.),( 0,(self.Num_Bin-1)),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
             #                 #bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth/2.),(0,None),(0.5,3.),( 0,(self.Num_Bin-1)),(1,powerwidth/2.),(0,None),(0.5,3.),(0.1,0.5))
             #
             #                 lsq2=fmin_l_bfgs_b(self.misfit2,state0,args=(y_data,x,num_intg),bounds=bnds,approx_grad=True)
             #
             #
             #                 chiSq2=lsq2[1]; jack2=self.y_jacobian2(x,lsq2[0])
             #
             #
             #                 try:
             #                     sigmas2=numpy.sqrt(chiSq2*numpy.diag(numpy.linalg.inv(numpy.dot(jack2.T,jack2))))
             #                 except:
             #                     std2a=32.; std2b=32.; sigmas2=numpy.ones(9)
             #                 else:
             #                     std2a=sigmas2[0]; std2b=sigmas2[4]
             #
             #
             #
             #                 oneG=(chiSq1<5 and chiSq1/chiSq2<2.0) and (abs(lsq2[0][0]-lsq2[0][4])<(lsq2[0][1]+lsq2[0][5])/3. or abs(lsq2[0][0]-lsq2[0][4])<10)
             #
             #                 if snrdB>-9: # when SNR is strong pick the peak with least shift (LOS velocity) error
             #                     if oneG:
             #                         choice=0
             #                     else:
             #                         w1=lsq2[0][1]; w2=lsq2[0][5]
             #                         a1=lsq2[0][2]; a2=lsq2[0][6]
             #                         p1=lsq2[0][3]; p2=lsq2[0][7]
             #                         s1=(2**(1+1./p1))*scipy.special.gamma(1./p1)/p1; s2=(2**(1+1./p2))*scipy.special.gamma(1./p2)/p2;
             #                         gp1=a1*w1*s1; gp2=a2*w2*s2 # power content of each ggaussian with proper p scaling
             #
             #                         if gp1>gp2:
             #                             if a1>0.7*a2:
             #                                 choice=1
             #                             else:
             #                                 choice=2
             #                         elif gp2>gp1:
             #                             if a2>0.7*a1:
             #                                 choice=2
             #                             else:
             #                                 choice=1
             #                         else:
             #                             choice=numpy.argmax([a1,a2])+1
             #                             #else:
             #                             #choice=argmin([std2a,std2b])+1
             #
             #                 else: # with low SNR go to the most energetic peak
             #                     choice=numpy.argmax([lsq1[0][2]*lsq1[0][1],lsq2[0][2]*lsq2[0][1],lsq2[0][6]*lsq2[0][5]])
             #
             #                 #print 'choice',choice
             #
             #                 if choice==0: # pick the single gaussian fit
             #                     Amplitude0=lsq1[0][2]
             #                     shift0=lsq1[0][0]
             #                     width0=lsq1[0][1]
             #                     p0=lsq1[0][3]
             #                     Amplitude1=0.
             #                     shift1=0.
             #                     width1=0.
             #                     p1=0.
             #                     noise=lsq1[0][4]
             #                 elif choice==1: # take the first one of the 2 gaussians fitted
             #                     Amplitude0 = lsq2[0][2]
             #                     shift0     = lsq2[0][0]
             #                     width0     = lsq2[0][1]
             #                     p0         = lsq2[0][3]
             #                     Amplitude1 = lsq2[0][6]     # This is 0 in gg1
             #                     shift1     = lsq2[0][4]     # This is 0 in gg1
             #                     width1     = lsq2[0][5]     # This is 0 in gg1
             #                     p1         = lsq2[0][7]     # This is 0 in gg1
             #                     noise      = lsq2[0][8]
             #                 else: # the second one
             #                     Amplitude0 = lsq2[0][6]
             #                     shift0     = lsq2[0][4]
             #                     width0     = lsq2[0][5]
             #                     p0         = lsq2[0][7]
             #                     Amplitude1 = lsq2[0][2]     # This is 0 in gg1
             #                     shift1     = lsq2[0][0]     # This is 0 in gg1
             #                     width1     = lsq2[0][1]     # This is 0 in gg1
             #                     p1         = lsq2[0][3]     # This is 0 in gg1
             #                     noise      = lsq2[0][8]
             #
             #                 #print len(noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0)
             #                 SPC_ch1[:,ht] = noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0
             #                 SPC_ch2[:,ht] = noise + Amplitude1*numpy.exp(-0.5*(abs(x-shift1))/width1)**p1
             #                 #print 'SPC_ch1.shape',SPC_ch1.shape
             #                 #print 'SPC_ch2.shape',SPC_ch2.shape
             #                 #dataOut.data_param = SPC_ch1
             #                 GauSPC[0] = SPC_ch1
             #                 GauSPC[1] = SPC_ch2
             #                 #plt.gcf().clear()
             #                plt.figure(50+self.i)
             #                self.i=self.i+1
             #                 #plt.subplot(121)
             #                plt.plot(self.spc,'k')#,label='spc(66)')
             #                plt.plot(SPC_ch1[ch,ht],'b')#,label='gg1')
             #                 #plt.plot(SPC_ch2,'r')#,label='gg2')
             #                 #plt.plot(xFrec,ySamples[1],'g',label='Ch1')
             #                 #plt.plot(xFrec,ySamples[2],'r',label='Ch2')
             #                 #plt.plot(xFrec,FitGauss,'yo:',label='fit')
             #                plt.legend()
             #                plt.title('DATOS A ALTURA DE 7500 METROS')
             #                plt.show()
             #             print 'shift0', shift0
             #             print 'Amplitude0', Amplitude0
             #             print 'width0', width0
             #             print 'p0', p0
             #             print '========================'
             #             print 'shift1', shift1
             #             print 'Amplitude1', Amplitude1
             #             print 'width1', width1
             #             print 'p1', p1
             #             print 'noise', noise
             #             print 's_noise', wnoise
                     print '========================================================'
                     print 'total_time: ', time.time()-start_time
                             # re-normalizing spc and noise
                             # This part differs from gg1
                     ''' Parameters:
 . Amplitude
 . Shift
 . Width
 . Power
                            '''
                     ###############################################################################
                 def FitGau(self, X):
                     Vrange, ch, pnoise, noise_, num_intg, SNRlimit = X
                     #print 'VARSSSS', ch, pnoise, noise, num_intg
                     #print 'HEIGHTS', self.Num_Hei
                     GauSPC = []
                     SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei])
                     SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei])
                     SPC_ch1[:] = 0#numpy.NaN
                     SPC_ch2[:] = 0#numpy.NaN
                     for ht in range(self.Num_Hei):
                         #print (numpy.asarray(self.spc).shape)
                         #print 'TTTTT', ch , ht
                         #print self.spc.shape
                         spc =  numpy.asarray(self.spc)[ch,:,ht]
                         #############################################
                         # normalizing spc and noise
                         # This part differs from gg1
                         spc_norm_max = max(spc)
                         spc = spc / spc_norm_max
                         pnoise = pnoise / spc_norm_max
                         #############################################
                         fatspectra=1.0
                         wnoise = noise_ / spc_norm_max
                             #wnoise,stdv,i_max,index =enoise(spc,num_intg) #noise estimate using Hildebrand Sekhon, only wnoise is used
                             #if wnoise>1.1*pnoise: # to be tested later
                             #    wnoise=pnoise
                         noisebl=wnoise*0.9; noisebh=wnoise*1.1
                         spc=spc-wnoise
                         # print 'wnoise', noise_[0], spc_norm_max, wnoise
                         minx=numpy.argmin(spc)
                         spcs=numpy.roll(spc,-minx)
                         cum=numpy.cumsum(spcs)
                         tot_noise=wnoise * self.Num_Bin  #64;
                         #print 'spc' , spcs[5:8] , 'tot_noise', tot_noise
                         #tot_signal=sum(cum[-5:])/5.; ''' How does this line work? '''
                         #snr=tot_signal/tot_noise
                         #snr=cum[-1]/tot_noise
                         snr = sum(spcs)/tot_noise
                         snrdB=10.*numpy.log10(snr)
                         if snrdB < SNRlimit :
                             snr = numpy.NaN
                             SPC_ch1[:,ht] = 0#numpy.NaN
                             SPC_ch1[:,ht] = 0#numpy.NaN
                             GauSPC = (SPC_ch1,SPC_ch2)
                             continue
                         #print 'snr',snrdB #, sum(spcs) , tot_noise
                         #if snrdB<-18 or numpy.isnan(snrdB) or num_intg<4:
                         #    return [None,]*4,[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None
                         cummax=max(cum); epsi=0.08*fatspectra # cumsum to narrow down the energy region
                         cumlo=cummax*epsi;
                         cumhi=cummax*(1-epsi)
                         powerindex=numpy.array(numpy.where(numpy.logical_and(cum>cumlo, cum<cumhi))[0])
                         if len(powerindex) < 1:# case for powerindex 0
                             continue
                         powerlo=powerindex[0]
                         powerhi=powerindex[-1]
                         powerwidth=powerhi-powerlo
                         firstpeak=powerlo+powerwidth/10.# first gaussian energy location
                         secondpeak=powerhi-powerwidth/10.#second gaussian energy location
                         midpeak=(firstpeak+secondpeak)/2.
                         firstamp=spcs[int(firstpeak)]
                         secondamp=spcs[int(secondpeak)]
                         midamp=spcs[int(midpeak)]
                         x=numpy.arange( self.Num_Bin )
                         y_data=spc+wnoise
                         # single gaussian
                         shift0=numpy.mod(midpeak+minx, self.Num_Bin )
                         width0=powerwidth/4.#Initialization entire power of spectrum divided by 4
                         power0=2.
                         amplitude0=midamp
                         state0=[shift0,width0,amplitude0,power0,wnoise]
                         bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth),(0,None),(0.5,3.),(noisebl,noisebh))
                         lsq1=fmin_l_bfgs_b(self.misfit1,state0,args=(y_data,x,num_intg),bounds=bnds,approx_grad=True)
                         chiSq1=lsq1[1];
                         jack1= self.y_jacobian1(x,lsq1[0])
                         try:
                             sigmas1=numpy.sqrt(chiSq1*numpy.diag(numpy.linalg.inv(numpy.dot(jack1.T,jack1))))
                         except:
                             std1=32.; sigmas1=numpy.ones(5)
                         else:
                             std1=sigmas1[0]
                         if fatspectra<1.0 and powerwidth<4:
                                 choice=0
                                 Amplitude0=lsq1[0][2]
                                 shift0=lsq1[0][0]
                                 width0=lsq1[0][1]
                                 p0=lsq1[0][3]
                                 Amplitude1=0.
                                 shift1=0.
                                 width1=0.
                                 p1=0.
                                 noise=lsq1[0][4]
                                 #return (numpy.array([shift0,width0,Amplitude0,p0]),
                                 #        numpy.array([shift1,width1,Amplitude1,p1]),noise,snrdB,chiSq1,6.,sigmas1,[None,]*9,choice)
                         # two gaussians
                         #shift0=numpy.mod(firstpeak+minx,64); shift1=numpy.mod(secondpeak+minx,64)
                         shift0=numpy.mod(firstpeak+minx, self.Num_Bin );
                         shift1=numpy.mod(secondpeak+minx, self.Num_Bin )
                         width0=powerwidth/6.;
                         width1=width0
                         power0=2.;
                         power1=power0
                         amplitude0=firstamp;
                         amplitude1=secondamp
                         state0=[shift0,width0,amplitude0,power0,shift1,width1,amplitude1,power1,wnoise]
                         #bnds=((0,63),(1,powerwidth/2.),(0,None),(0.5,3.),(0,63),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
                         bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth/2.),(0,None),(0.5,3.),( 0,(self.Num_Bin-1)),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
                         #bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth/2.),(0,None),(0.5,3.),( 0,(self.Num_Bin-1)),(1,powerwidth/2.),(0,None),(0.5,3.),(0.1,0.5))
                         lsq2=fmin_l_bfgs_b(self.misfit2,state0,args=(y_data,x,num_intg),bounds=bnds,approx_grad=True)
                         chiSq2=lsq2[1]; jack2=self.y_jacobian2(x,lsq2[0])
                         try:
                             sigmas2=numpy.sqrt(chiSq2*numpy.diag(numpy.linalg.inv(numpy.dot(jack2.T,jack2))))
                         except:
                             std2a=32.; std2b=32.; sigmas2=numpy.ones(9)
                         else:
                             std2a=sigmas2[0]; std2b=sigmas2[4]
                         oneG=(chiSq1<5 and chiSq1/chiSq2<2.0) and (abs(lsq2[0][0]-lsq2[0][4])<(lsq2[0][1]+lsq2[0][5])/3. or abs(lsq2[0][0]-lsq2[0][4])<10)
                         if snrdB>-6: # when SNR is strong pick the peak with least shift (LOS velocity) error
                             if oneG:
                                 choice=0
                             else:
                                 w1=lsq2[0][1]; w2=lsq2[0][5]
                                 a1=lsq2[0][2]; a2=lsq2[0][6]
                                 p1=lsq2[0][3]; p2=lsq2[0][7]
                                 s1=(2**(1+1./p1))*scipy.special.gamma(1./p1)/p1;
                                 s2=(2**(1+1./p2))*scipy.special.gamma(1./p2)/p2;
                                 gp1=a1*w1*s1; gp2=a2*w2*s2 # power content of each ggaussian with proper p scaling
                                 if gp1>gp2:
                                     if a1>0.7*a2:
                                         choice=1
                                     else:
                                         choice=2
                                 elif gp2>gp1:
                                     if a2>0.7*a1:
                                         choice=2
                                     else:
                                         choice=1
                                 else:
                                     choice=numpy.argmax([a1,a2])+1
                                     #else:
                                     #choice=argmin([std2a,std2b])+1
                         else: # with low SNR go to the most energetic peak
                             choice=numpy.argmax([lsq1[0][2]*lsq1[0][1],lsq2[0][2]*lsq2[0][1],lsq2[0][6]*lsq2[0][5]])
                         shift0=lsq2[0][0]; vel0=Vrange[0] + shift0*(Vrange[1]-Vrange[0])
                         shift1=lsq2[0][4]; vel1=Vrange[0] + shift1*(Vrange[1]-Vrange[0])
                         max_vel = 20
                         #first peak will be 0, second peak will be 1
                         if vel0 > 0 and vel0 < max_vel : #first peak is in the correct range
                             shift0=lsq2[0][0]
                             width0=lsq2[0][1]
                             Amplitude0=lsq2[0][2]
                             p0=lsq2[0][3]
                             shift1=lsq2[0][4]
                             width1=lsq2[0][5]
                             Amplitude1=lsq2[0][6]
                             p1=lsq2[0][7]
                             noise=lsq2[0][8]
                         else:
                             shift1=lsq2[0][0]
                             width1=lsq2[0][1]
                             Amplitude1=lsq2[0][2]
                             p1=lsq2[0][3]
                             shift0=lsq2[0][4]
                             width0=lsq2[0][5]
                             Amplitude0=lsq2[0][6]
                             p0=lsq2[0][7]
                             noise=lsq2[0][8]
                         if Amplitude0<0.1: # in case the peak is noise
                             shift0,width0,Amplitude0,p0 = 4*[numpy.NaN]
                         if Amplitude1<0.1:
                             shift1,width1,Amplitude1,p1 = 4*[numpy.NaN]
             #             if choice==0: # pick the single gaussian fit
             #                 Amplitude0=lsq1[0][2]
             #                 shift0=lsq1[0][0]
             #                 width0=lsq1[0][1]
             #                 p0=lsq1[0][3]
             #                 Amplitude1=0.
             #                 shift1=0.
             #                 width1=0.
             #                 p1=0.
             #                 noise=lsq1[0][4]
             #             elif choice==1: # take the first one of the 2 gaussians fitted
             #                 Amplitude0 = lsq2[0][2]
             #                 shift0     = lsq2[0][0]
             #                 width0     = lsq2[0][1]
             #                 p0         = lsq2[0][3]
             #                 Amplitude1 = lsq2[0][6]     # This is 0 in gg1
             #                 shift1     = lsq2[0][4]     # This is 0 in gg1
             #                 width1     = lsq2[0][5]     # This is 0 in gg1
             #                 p1         = lsq2[0][7]     # This is 0 in gg1
             #                 noise      = lsq2[0][8]
             #             else: # the second one
             #                 Amplitude0 = lsq2[0][6]
             #                 shift0     = lsq2[0][4]
             #                 width0     = lsq2[0][5]
             #                 p0         = lsq2[0][7]
             #                 Amplitude1 = lsq2[0][2]     # This is 0 in gg1
             #                 shift1     = lsq2[0][0]     # This is 0 in gg1
             #                 width1     = lsq2[0][1]     # This is 0 in gg1
             #                 p1         = lsq2[0][3]     # This is 0 in gg1
             #                 noise      = lsq2[0][8]
                         #print len(noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0)
                         SPC_ch1[:,ht] = noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0
                         SPC_ch2[:,ht] = noise + Amplitude1*numpy.exp(-0.5*(abs(x-shift1))/width1)**p1
                         #print 'SPC_ch1.shape',SPC_ch1.shape
                         #print 'SPC_ch2.shape',SPC_ch2.shape
                         #dataOut.data_param = SPC_ch1
                         GauSPC = (SPC_ch1,SPC_ch2)
                         #GauSPC[1] = SPC_ch2
             #         print 'shift0', shift0
             #         print 'Amplitude0', Amplitude0
             #         print 'width0', width0
             #         print 'p0', p0
             #         print '========================'
             #         print 'shift1', shift1
             #         print 'Amplitude1', Amplitude1
             #         print 'width1', width1
             #         print 'p1', p1
             #         print 'noise', noise
             #         print 's_noise', wnoise
                     return GauSPC
                 def y_jacobian1(self,x,state): # This function is for further analysis of generalized Gaussians, it is not too importan for the signal discrimination.
                     y_model=self.y_model1(x,state)
                     s0,w0,a0,p0,n=state
                     e0=((x-s0)/w0)**2;
                     e0u=((x-s0-self.Num_Bin)/w0)**2;
                     e0d=((x-s0+self.Num_Bin)/w0)**2
                     m0=numpy.exp(-0.5*e0**(p0/2.));
                     m0u=numpy.exp(-0.5*e0u**(p0/2.));
                     m0d=numpy.exp(-0.5*e0d**(p0/2.))
                     JA=m0+m0u+m0d
                     JP=(-1/4.)*a0*m0*e0**(p0/2.)*numpy.log(e0)+(-1/4.)*a0*m0u*e0u**(p0/2.)*numpy.log(e0u)+(-1/4.)*a0*m0d*e0d**(p0/2.)*numpy.log(e0d)
                     JS=(p0/w0/2.)*a0*m0*e0**(p0/2.-1)*((x-s0)/w0)+(p0/w0/2.)*a0*m0u*e0u**(p0/2.-1)*((x-s0- self.Num_Bin )/w0)+(p0/w0/2.)*a0*m0d*e0d**(p0/2.-1)*((x-s0+ self.Num_Bin )/w0)
                     JW=(p0/w0/2.)*a0*m0*e0**(p0/2.-1)*((x-s0)/w0)**2+(p0/w0/2.)*a0*m0u*e0u**(p0/2.-1)*((x-s0- self.Num_Bin )/w0)**2+(p0/w0/2.)*a0*m0d*e0d**(p0/2.-1)*((x-s0+ self.Num_Bin )/w0)**2
                     jack1=numpy.sqrt(7)*numpy.array([JS/y_model,JW/y_model,JA/y_model,JP/y_model,1./y_model])
                     return jack1.T
                 def y_jacobian2(self,x,state):
                     y_model=self.y_model2(x,state)
                     s0,w0,a0,p0,s1,w1,a1,p1,n=state
                     e0=((x-s0)/w0)**2;
                     e0u=((x-s0- self.Num_Bin )/w0)**2;
                     e0d=((x-s0+ self.Num_Bin )/w0)**2
                     e1=((x-s1)/w1)**2;
                     e1u=((x-s1- self.Num_Bin )/w1)**2;
                     e1d=((x-s1+ self.Num_Bin )/w1)**2
                     m0=numpy.exp(-0.5*e0**(p0/2.));
                     m0u=numpy.exp(-0.5*e0u**(p0/2.));
                     m0d=numpy.exp(-0.5*e0d**(p0/2.))
                     m1=numpy.exp(-0.5*e1**(p1/2.));
                     m1u=numpy.exp(-0.5*e1u**(p1/2.));
                     m1d=numpy.exp(-0.5*e1d**(p1/2.))
                     JA=m0+m0u+m0d
                     JA1=m1+m1u+m1d
                     JP=(-1/4.)*a0*m0*e0**(p0/2.)*numpy.log(e0)+(-1/4.)*a0*m0u*e0u**(p0/2.)*numpy.log(e0u)+(-1/4.)*a0*m0d*e0d**(p0/2.)*numpy.log(e0d)
                     JP1=(-1/4.)*a1*m1*e1**(p1/2.)*numpy.log(e1)+(-1/4.)*a1*m1u*e1u**(p1/2.)*numpy.log(e1u)+(-1/4.)*a1*m1d*e1d**(p1/2.)*numpy.log(e1d)
                     JS=(p0/w0/2.)*a0*m0*e0**(p0/2.-1)*((x-s0)/w0)+(p0/w0/2.)*a0*m0u*e0u**(p0/2.-1)*((x-s0- self.Num_Bin )/w0)+(p0/w0/2.)*a0*m0d*e0d**(p0/2.-1)*((x-s0+ self.Num_Bin )/w0)
                     JS1=(p1/w1/2.)*a1*m1*e1**(p1/2.-1)*((x-s1)/w1)+(p1/w1/2.)*a1*m1u*e1u**(p1/2.-1)*((x-s1- self.Num_Bin )/w1)+(p1/w1/2.)*a1*m1d*e1d**(p1/2.-1)*((x-s1+ self.Num_Bin )/w1)
                     JW=(p0/w0/2.)*a0*m0*e0**(p0/2.-1)*((x-s0)/w0)**2+(p0/w0/2.)*a0*m0u*e0u**(p0/2.-1)*((x-s0- self.Num_Bin )/w0)**2+(p0/w0/2.)*a0*m0d*e0d**(p0/2.-1)*((x-s0+ self.Num_Bin )/w0)**2
                     JW1=(p1/w1/2.)*a1*m1*e1**(p1/2.-1)*((x-s1)/w1)**2+(p1/w1/2.)*a1*m1u*e1u**(p1/2.-1)*((x-s1- self.Num_Bin )/w1)**2+(p1/w1/2.)*a1*m1d*e1d**(p1/2.-1)*((x-s1+ self.Num_Bin )/w1)**2
                     jack2=numpy.sqrt(7)*numpy.array([JS/y_model,JW/y_model,JA/y_model,JP/y_model,JS1/y_model,JW1/y_model,JA1/y_model,JP1/y_model,1./y_model])
                     return jack2.T
                 def y_model1(self,x,state):
                     shift0,width0,amplitude0,power0,noise=state
                     model0=amplitude0*numpy.exp(-0.5*abs((x-shift0)/width0)**power0)
                     model0u=amplitude0*numpy.exp(-0.5*abs((x-shift0- self.Num_Bin )/width0)**power0)
                     model0d=amplitude0*numpy.exp(-0.5*abs((x-shift0+ self.Num_Bin )/width0)**power0)
                     return model0+model0u+model0d+noise
                 def y_model2(self,x,state): #Equation for two generalized Gaussians with Nyquist
                     shift0,width0,amplitude0,power0,shift1,width1,amplitude1,power1,noise=state
                     model0=amplitude0*numpy.exp(-0.5*abs((x-shift0)/width0)**power0)
                     model0u=amplitude0*numpy.exp(-0.5*abs((x-shift0- self.Num_Bin )/width0)**power0)
                     model0d=amplitude0*numpy.exp(-0.5*abs((x-shift0+ self.Num_Bin )/width0)**power0)
                     model1=amplitude1*numpy.exp(-0.5*abs((x-shift1)/width1)**power1)
                     model1u=amplitude1*numpy.exp(-0.5*abs((x-shift1- self.Num_Bin )/width1)**power1)
                     model1d=amplitude1*numpy.exp(-0.5*abs((x-shift1+ self.Num_Bin )/width1)**power1)
                     return model0+model0u+model0d+model1+model1u+model1d+noise
                 def misfit1(self,state,y_data,x,num_intg): # This function compares how close real data is with the model data, the close it is, the better it is.
                     return num_intg*sum((numpy.log(y_data)-numpy.log(self.y_model1(x,state)))**2)#/(64-5.) # /(64-5.) can be commented
                 def misfit2(self,state,y_data,x,num_intg):
                     return num_intg*sum((numpy.log(y_data)-numpy.log(self.y_model2(x,state)))**2)#/(64-9.)
             class PrecipitationProc(Operation):
                 '''
                      Operator that estimates Reflectivity factor (Z), and estimates rainfall Rate (R)
                      Input:
                         self.dataOut.data_pre    :    SelfSpectra
                      Output:
                         self.dataOut.data_output :    Reflectivity factor, rainfall Rate
                      Parameters affected:
                 '''
                 def run(self, dataOut, radar=None, Pt=None, Gt=None, Gr=None, Lambda=None, aL=None,
                         tauW=None, ThetaT=None, ThetaR=None, Km = 0.93, Altitude=None):
                     self.spc = dataOut.data_pre[0].copy()
                     self.Num_Hei = self.spc.shape[2]
                     self.Num_Bin = self.spc.shape[1]
                     self.Num_Chn = self.spc.shape[0]
                     Velrange = dataOut.abscissaList
                     if radar == "MIRA35C" :
                         Ze = self.dBZeMODE2(dataOut)
                     else:
                         self.Pt = Pt
                         self.Gt = Gt
                         self.Gr = Gr
                         self.Lambda = Lambda
                         self.aL = aL
                         self.tauW = tauW
                         self.ThetaT = ThetaT
                         self.ThetaR = ThetaR
                         RadarConstant = GetRadarConstant()
                         SPCmean = numpy.mean(self.spc,0)
                         ETA = numpy.zeros(self.Num_Hei)
                         Pr = numpy.sum(SPCmean,0)
                         #for R in range(self.Num_Hei):
                         #    ETA[R] = RadarConstant * Pr[R] * R**2 #Reflectivity (ETA)
                         D_range = numpy.zeros(self.Num_Hei)
                         EqSec = numpy.zeros(self.Num_Hei)
                         del_V = numpy.zeros(self.Num_Hei)
                         for R in range(self.Num_Hei):
                             ETA[R] = RadarConstant * Pr[R] * R**2 #Reflectivity (ETA)
                             h = R + Altitude #Range from ground to radar pulse altitude
                             del_V[R] = 1 + 3.68 * 10**-5 * h + 1.71 * 10**-9 * h**2    #Density change correction for velocity
                             D_range[R] = numpy.log( (9.65 - (Velrange[R]/del_V[R])) / 10.3 ) / -0.6 #Range of Diameter of drops related to velocity
                             SIGMA[R] = numpy.pi**5 / Lambda**4 * Km * D_range[R]**6     #Equivalent Section of drops (sigma)
                             N_dist[R] = ETA[R] / SIGMA[R]
                         Ze = (ETA * Lambda**4) / (numpy.pi * Km)
                         Z = numpy.sum( N_dist * D_range**6 )
                         RR = 6*10**-4*numpy.pi * numpy.sum( D_range**3 * N_dist * Velrange ) #Rainfall rate
                     RR = (Ze/200)**(1/1.6)
                     dBRR = 10*numpy.log10(RR)
                     dBZe = 10*numpy.log10(Ze)
                     dataOut.data_output = Ze
                     dataOut.data_param = numpy.ones([2,self.Num_Hei])
                     dataOut.channelList = [0,1]
                     print 'channelList', dataOut.channelList
                     dataOut.data_param[0]=dBZe
                     dataOut.data_param[1]=dBRR
                     print 'RR SHAPE', dBRR.shape
                     print 'Ze SHAPE', dBZe.shape
                     print 'dataOut.data_param SHAPE', dataOut.data_param.shape
                 def dBZeMODE2(self, dataOut): #    Processing for MIRA35C
                     NPW = dataOut.NPW
                     COFA = dataOut.COFA
                     SNR = numpy.array([self.spc[0,:,:] / NPW[0]]) #, self.spc[1,:,:] / NPW[1]])
                     RadarConst = dataOut.RadarConst
                     #frequency = 34.85*10**9
                     ETA = numpy.zeros(([self.Num_Chn ,self.Num_Hei]))
                     data_output = numpy.ones([self.Num_Chn , self.Num_Hei])*numpy.NaN
                     ETA = numpy.sum(SNR,1)
                     print 'ETA' , ETA
                     ETA = numpy.where(ETA is not 0. , ETA, numpy.NaN)
                     Ze = numpy.ones([self.Num_Chn, self.Num_Hei] )
                     for r in range(self.Num_Hei):
                         Ze[0,r] =  ( ETA[0,r] ) * COFA[0,r][0] * RadarConst * ((r/5000.)**2)
                         #Ze[1,r] =  ( ETA[1,r] ) * COFA[1,r][0] * RadarConst * ((r/5000.)**2)
                     return Ze
                 def GetRadarConstant(self):
                     """
                     Constants:
                     Pt:     Transmission Power               dB
                     Gt:     Transmission Gain                dB
                     Gr:     Reception Gain                   dB
                     Lambda: Wavelenght                       m
                     aL:      Attenuation loses               dB
                     tauW:   Width of transmission pulse      s
                     ThetaT: Transmission antenna bean angle  rad
                     ThetaR: Reception antenna beam angle     rad
                     """
                     Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) )
                     Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * TauW * numpy.pi * ThetaT * TheraR)
                     RadarConstant =  Numerator / Denominator
                     return RadarConstant
             class FullSpectralAnalysis(Operation):
                 """
                     Function that implements Full Spectral Analisys technique.
                     Input:
                         self.dataOut.data_pre    :    SelfSpectra and CrossSPectra data
                         self.dataOut.groupList   :    Pairlist of channels
                         self.dataOut.ChanDist    :    Physical distance between receivers
                     Output:
                         self.dataOut.data_output :    Zonal wind, Meridional wind and Vertical wind
                     Parameters affected:    Winds, height range, SNR
                 """
-                def run(self, dataOut, E01=None, E02=None, E12=None, N01=None, N02=None, N12=None):
+                def run(self, dataOut, E01=None, E02=None, E12=None, N01=None, N02=None, N12=None, SNRlimit=7):
                     spc = dataOut.data_pre[0].copy()
                     cspc = dataOut.data_pre[1].copy()
                     nChannel = spc.shape[0]
                     nProfiles = spc.shape[1]
                     nHeights = spc.shape[2]
                     pairsList = dataOut.groupList
                     if dataOut.ChanDist is not None :
                         ChanDist = dataOut.ChanDist
                     else:
                         ChanDist = numpy.array([[E01, N01],[E02,N02],[E12,N12]])
                     #print 'ChanDist', ChanDist
                     if dataOut.VelRange is not None:
                         VelRange= dataOut.VelRange
                     else:
                         VelRange= dataOut.abscissaList
                     ySamples=numpy.ones([nChannel,nProfiles])
                     phase=numpy.ones([nChannel,nProfiles])
                     CSPCSamples=numpy.ones([nChannel,nProfiles],dtype=numpy.complex_)
                     coherence=numpy.ones([nChannel,nProfiles])
                     PhaseSlope=numpy.ones(nChannel)
                     PhaseInter=numpy.ones(nChannel)
+                    dataSNR = dataOut.data_SNR
                     data = dataOut.data_pre
                     noise = dataOut.noise
                     print 'noise',noise
-                    SNRdB = 10*numpy.log10(dataOut.data_SNR)
+                    #SNRdB = 10*numpy.log10(dataOut.data_SNR)
-                    FirstMoment = []
+                    FirstMoment = numpy.average(dataOut.data_param[:,1,:],0)
-                    SNRdBMean = []
+                    #SNRdBMean = []
-                    for j in range(nHeights):
+                    #for j in range(nHeights):
-                        FirstMoment = numpy.append(FirstMoment,numpy.mean([dataOut.data_param[0,1,j],dataOut.data_param[1,1,j],dataOut.data_param[2,1,j]]))
+                    #    FirstMoment = numpy.append(FirstMoment,numpy.mean([dataOut.data_param[0,1,j],dataOut.data_param[1,1,j],dataOut.data_param[2,1,j]]))
-                        SNRdBMean = numpy.append(SNRdBMean,numpy.mean([SNRdB[0,j],SNRdB[1,j],SNRdB[2,j]]))
+                    #    SNRdBMean = numpy.append(SNRdBMean,numpy.mean([SNRdB[0,j],SNRdB[1,j],SNRdB[2,j]]))
                     data_output=numpy.ones([3,spc.shape[2]])*numpy.NaN
                     velocityX=[]
                     velocityY=[]
                     velocityV=[]
+                    dbSNR = 10*numpy.log10(dataSNR)
+                    dbSNR = numpy.average(dbSNR,0)
                     for Height in range(nHeights):
-                        [Vzon,Vmer,Vver, GaussCenter]= self.WindEstimation(spc, cspc, pairsList, ChanDist, Height, noise, VelRange)
+                        [Vzon,Vmer,Vver, GaussCenter]= self.WindEstimation(spc, cspc, pairsList, ChanDist, Height, noise, VelRange, dbSNR[Height], SNRlimit)
                         if abs(Vzon)<100. and abs(Vzon)> 0.:
                             velocityX=numpy.append(velocityX, Vzon)#Vmag
                         else:
                             print 'Vzon',Vzon
                             velocityX=numpy.append(velocityX, numpy.NaN)
                         if abs(Vmer)<100. and abs(Vmer) > 0.:
                             velocityY=numpy.append(velocityY, Vmer)#Vang
                         else:
                             print 'Vmer',Vmer
                             velocityY=numpy.append(velocityY, numpy.NaN)
-                        if abs(GaussCenter)<10:
+                        if dbSNR[Height] > SNRlimit:
-                            velocityV=numpy.append(velocityV, Vver)
+                            velocityV=numpy.append(velocityV, FirstMoment[Height])
                         else:
                             velocityV=numpy.append(velocityV, numpy.NaN)
                             #FirstMoment[Height]= numpy.NaN
             #             if SNRdBMean[Height]  <12:
             #                 FirstMoment[Height] = numpy.NaN
             #                 velocityX[Height] = numpy.NaN
             #                 velocityY[Height] = numpy.NaN
                     data_output[0]=numpy.array(velocityX)
                     data_output[1]=numpy.array(velocityY)
-                    data_output[2]=-FirstMoment
+                    data_output[2]=-velocityV#FirstMoment
                     print ' '
                     #print 'FirstMoment'
                     #print FirstMoment
                     print 'velocityX',data_output[0]
                     print ' '
                     print 'velocityY',data_output[1]
                     #print numpy.array(velocityY)
                     print ' '
                     #print 'SNR'
                     #print 10*numpy.log10(dataOut.data_SNR)
                     #print numpy.shape(10*numpy.log10(dataOut.data_SNR))
                     print ' '
                     dataOut.data_output=data_output
                     return
                 def moving_average(self,x, N=2):
                     return numpy.convolve(x, numpy.ones((N,))/N)[(N-1):]
                 def gaus(self,xSamples,a,x0,sigma):
                     return a*numpy.exp(-(xSamples-x0)**2/(2*sigma**2))
                 def Find(self,x,value):
                     for index in range(len(x)):
                         if x[index]==value:
                             return index
-                def WindEstimation(self, spc, cspc, pairsList, ChanDist, Height, noise, VelRange):
+                def WindEstimation(self, spc, cspc, pairsList, ChanDist, Height, noise, VelRange, dbSNR, SNRlimit):
                     ySamples=numpy.ones([spc.shape[0],spc.shape[1]])
                     phase=numpy.ones([spc.shape[0],spc.shape[1]])
                     CSPCSamples=numpy.ones([spc.shape[0],spc.shape[1]],dtype=numpy.complex_)
                     coherence=numpy.ones([spc.shape[0],spc.shape[1]])
                     PhaseSlope=numpy.ones(spc.shape[0])
                     PhaseInter=numpy.ones(spc.shape[0])
                     xFrec=VelRange
                     '''Getting Eij and Nij'''
                     E01=ChanDist[0][0]
                     N01=ChanDist[0][1]
                     E02=ChanDist[1][0]
                     N02=ChanDist[1][1]
                     E12=ChanDist[2][0]
                     N12=ChanDist[2][1]
                     z = spc.copy()
                     z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
                     for i in range(spc.shape[0]):
                         '''****** Line of Data SPC ******'''
                         zline=z[i,:,Height]
                         '''****** SPC is normalized ******'''
-                        FactNorm= zline.copy() / numpy.sum(zline.copy())
+                        FactNorm= (zline.copy()-noise[i]) / numpy.sum(zline.copy())
                         FactNorm= FactNorm/numpy.sum(FactNorm)
                         SmoothSPC=self.moving_average(FactNorm,N=3)
                         xSamples = ar(range(len(SmoothSPC)))
-                        ySamples[i] = SmoothSPC-noise[i]
+                        ySamples[i] = SmoothSPC
+                    #dbSNR=10*numpy.log10(dataSNR)
                     print ' '
                     print ' '
                     print ' '
-                    print 'SmoothSPC',SmoothSPC
+                    #print 'dataSNR', dbSNR.shape, dbSNR[0,40:120]
+                    print 'SmoothSPC', SmoothSPC.shape, SmoothSPC[0:20]
                     print 'noise',noise
-                    print'zline',zline
+                    print 'zline',zline.shape, zline[0:20]
-                    print'FactNorm',FactNorm
+                    print 'FactNorm',FactNorm.shape, FactNorm[0:20]
+                    print 'FactNorm suma', numpy.sum(FactNorm)
                     for i in range(spc.shape[0]):
                         '''****** Line of Data CSPC ******'''
                         cspcLine=cspc[i,:,Height].copy()
                         '''****** CSPC is normalized ******'''
                         chan_index0 = pairsList[i][0]
                         chan_index1 = pairsList[i][1]
-                        CSPCFactor= numpy.sum(ySamples[chan_index0]) * numpy.sum(ySamples[chan_index1])
+                        CSPCFactor= abs(numpy.sum(ySamples[chan_index0]) * numpy.sum(ySamples[chan_index1])) #
-                        CSPCNorm= cspcLine.copy() / numpy.sqrt(CSPCFactor)
+                        CSPCNorm = (cspcLine.copy() -noise[i]) / numpy.sqrt(CSPCFactor)
-                        CSPCSamples[i] = CSPCNorm-noise[i]
+                        CSPCSamples[i] = CSPCNorm
                         coherence[i] = numpy.abs(CSPCSamples[i]) / numpy.sqrt(CSPCFactor)
                         coherence[i]= self.moving_average(coherence[i],N=2)
                         phase[i] = self.moving_average( numpy.arctan2(CSPCSamples[i].imag, CSPCSamples[i].real),N=1)#*180/numpy.pi
-                    print 'CSPCSamples', CSPCSamples
+                    print 'cspcLine', cspcLine.shape, cspcLine[0:20]
+                    print 'CSPCFactor', CSPCFactor#, CSPCFactor[0:20]
+                    print numpy.sum(ySamples[chan_index0]), numpy.sum(ySamples[chan_index1]), -noise[i]
+                    print 'CSPCNorm', CSPCNorm.shape, CSPCNorm[0:20]
+                    print 'CSPCNorm suma', numpy.sum(CSPCNorm)
+                    print 'CSPCSamples', CSPCSamples.shape, CSPCSamples[0,0:20]
                     '''****** Getting fij width ******'''
                     yMean=[]
                     yMean2=[]
                     for j in range(len(ySamples[1])):
                         yMean=numpy.append(yMean,numpy.mean([ySamples[0,j],ySamples[1,j],ySamples[2,j]]))
                     '''******* Getting fitting Gaussian ******'''
                     meanGauss=sum(xSamples*yMean) / len(xSamples)
                     sigma=sum(yMean*(xSamples-meanGauss)**2) / len(xSamples)
                     print '****************************'
                     print 'len(xSamples): ',len(xSamples)
-                    print 'yMean: ', yMean
+                    print 'yMean: ', yMean.shape, yMean[0:20]
-                    print 'ySamples', ySamples
+                    print 'ySamples', ySamples.shape, ySamples[0,0:20]
-                    print 'xSamples: ',xSamples
+                    print 'xSamples: ',xSamples.shape, xSamples[0:20]
                     print 'meanGauss',meanGauss
                     print 'sigma',sigma
-                    if (abs(meanGauss/sigma**2) > 0.000000001):#0.00001) :
+                    #if (abs(meanGauss/sigma**2) > 0.0001) : #0.000000001):
+                    if dbSNR > SNRlimit :
                         try:
                             popt,pcov = curve_fit(self.gaus,xSamples,yMean,p0=[1,meanGauss,sigma])
                             if numpy.amax(popt)>numpy.amax(yMean)*0.3:
                                 FitGauss=self.gaus(xSamples,*popt)
                             else:
                                 FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
                                 print 'Verificador:     Dentro', Height
-                        except RuntimeError:
+                        except :#RuntimeError:
                             FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
                     else:
                         FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
                     Maximun=numpy.amax(yMean)
-                    eMinus1=Maximun*numpy.exp(-1)*0.8
+                    eMinus1=Maximun*numpy.exp(-1)#*0.8
                     HWpos=self.Find(FitGauss,min(FitGauss, key=lambda value:abs(value-eMinus1)))
                     HalfWidth= xFrec[HWpos]
                     GCpos=self.Find(FitGauss, numpy.amax(FitGauss))
                     Vpos=self.Find(FactNorm, numpy.amax(FactNorm))
                     #Vpos=FirstMoment[]
                     '''****** Getting Fij ******'''
                     GaussCenter=xFrec[GCpos]
                     if (GaussCenter<0 and HalfWidth>0) or (GaussCenter>0 and HalfWidth<0):
                         Fij=abs(GaussCenter)+abs(HalfWidth)+0.0000001
                     else:
                         Fij=abs(GaussCenter-HalfWidth)+0.0000001
                     '''****** Getting Frecuency range of significant data ******'''
                     Rangpos=self.Find(FitGauss,min(FitGauss, key=lambda value:abs(value-Maximun*0.10)))
                     if Rangpos<GCpos:
                         Range=numpy.array([Rangpos,2*GCpos-Rangpos])
-                    else:
+                    elif Rangpos< ( len(xFrec)- len(xFrec)*0.1):
                         Range=numpy.array([2*GCpos-Rangpos,Rangpos])
+                    else:
+                        Range = numpy.array([0,0])
+                    print ' '
+                    print 'GCpos',GCpos, ( len(xFrec)- len(xFrec)*0.1)
+                    print 'Rangpos',Rangpos
+                    print 'RANGE: ', Range
                     FrecRange=xFrec[Range[0]:Range[1]]
                     '''****** Getting SCPC Slope ******'''
                     for i in range(spc.shape[0]):
                         if len(FrecRange)>5 and len(FrecRange)<spc.shape[1]*0.5:
                             PhaseRange=self.moving_average(phase[i,Range[0]:Range[1]],N=3)
-                            slope, intercept, r_value, p_value, std_err = stats.linregress(FrecRange,PhaseRange)
+                            print 'FrecRange', len(FrecRange) , FrecRange
-                            PhaseSlope[i]=slope
+                            print 'PhaseRange', len(PhaseRange), PhaseRange
-                            PhaseInter[i]=intercept
+                            print ' '
+                            if len(FrecRange) == len(PhaseRange):
+                                slope, intercept, r_value, p_value, std_err = stats.linregress(FrecRange,PhaseRange)
+                                PhaseSlope[i]=slope
+                                PhaseInter[i]=intercept
+                            else:
+                                PhaseSlope[i]=0
+                                PhaseInter[i]=0
                         else:
                             PhaseSlope[i]=0
                             PhaseInter[i]=0
                         '''Getting constant C'''
                         cC=(Fij*numpy.pi)**2
                         '''****** Getting constants F and G ******'''
                         MijEijNij=numpy.array([[E02,N02], [E12,N12]])
                         MijResult0=(-PhaseSlope[1]*cC) / (2*numpy.pi)
                         MijResult1=(-PhaseSlope[2]*cC) / (2*numpy.pi)
                         MijResults=numpy.array([MijResult0,MijResult1])
                         (cF,cG) = numpy.linalg.solve(MijEijNij, MijResults)
                         '''****** Getting constants A, B and H ******'''
                         W01=numpy.amax(coherence[0])
                         W02=numpy.amax(coherence[1])
                         W12=numpy.amax(coherence[2])
                         WijResult0=((cF*E01+cG*N01)**2)/cC - numpy.log(W01 / numpy.sqrt(numpy.pi/cC))
                         WijResult1=((cF*E02+cG*N02)**2)/cC - numpy.log(W02 / numpy.sqrt(numpy.pi/cC))
                         WijResult2=((cF*E12+cG*N12)**2)/cC - numpy.log(W12 / numpy.sqrt(numpy.pi/cC))
                         WijResults=numpy.array([WijResult0, WijResult1, WijResult2])
                         WijEijNij=numpy.array([ [E01**2, N01**2, 2*E01*N01] , [E02**2, N02**2, 2*E02*N02] , [E12**2, N12**2, 2*E12*N12] ])
                         (cA,cB,cH) = numpy.linalg.solve(WijEijNij, WijResults)
                         VxVy=numpy.array([[cA,cH],[cH,cB]])
                         VxVyResults=numpy.array([-cF,-cG])
                         (Vx,Vy) = numpy.linalg.solve(VxVy, VxVyResults)
                         Vzon = Vy
                         Vmer = Vx
                         Vmag=numpy.sqrt(Vzon**2+Vmer**2)
                         Vang=numpy.arctan2(Vmer,Vzon)
                         Vver=xFrec[Vpos]
                     print 'vzon y vmer', Vzon, Vmer
                     return Vzon, Vmer, Vver, GaussCenter
             class SpectralMoments(Operation):
                 '''
                     Function SpectralMoments()
                     Calculates moments (power, mean, standard deviation) and SNR of the signal
                     Type of dataIn:    Spectra
                     Configuration Parameters:
                         dirCosx    :     Cosine director in X axis
                         dirCosy    :     Cosine director in Y axis
                         elevation  :
                         azimuth    :
                     Input:
                         channelList    :    simple channel list to select e.g. [2,3,7]
                         self.dataOut.data_pre        :    Spectral data
                         self.dataOut.abscissaList    :    List of frequencies
                         self.dataOut.noise           :    Noise level per channel
                     Affected:
                         self.dataOut.data_param      :    Parameters per channel
                         self.dataOut.data_SNR        :    SNR per channel
                 '''
                 def run(self, dataOut):
                     #dataOut.data_pre = dataOut.data_pre[0]
                     data = dataOut.data_pre[0]
                     absc = dataOut.abscissaList[:-1]
                     noise = dataOut.noise
                     nChannel = data.shape[0]
                     data_param = numpy.zeros((nChannel, 4, data.shape[2]))
                     for ind in range(nChannel):
                         data_param[ind,:,:] = self.__calculateMoments( data[ind,:,:] , absc , noise[ind] )
                     dataOut.data_param = data_param[:,1:,:]
                     dataOut.data_SNR = data_param[:,0]
                     return
                 def __calculateMoments(self, oldspec, oldfreq, n0,
                                        nicoh = None, graph = None, smooth = None, type1 = None, fwindow = None, snrth = None, dc = None, aliasing = None, oldfd = None, wwauto = None):
                     if (nicoh == None): nicoh = 1
                     if (graph == None): graph = 0
                     if (smooth == None): smooth = 0
                     elif (self.smooth < 3): smooth = 0
                     if (type1 == None): type1 = 0
                     if (fwindow == None): fwindow = numpy.zeros(oldfreq.size) + 1
                     if (snrth == None): snrth = -3
                     if (dc == None): dc = 0
                     if (aliasing == None): aliasing = 0
                     if (oldfd == None): oldfd = 0
                     if (wwauto == None): wwauto = 0
                     if (n0 < 1.e-20):   n0 = 1.e-20
                     freq = oldfreq
                     vec_power = numpy.zeros(oldspec.shape[1])
                     vec_fd = numpy.zeros(oldspec.shape[1])
                     vec_w = numpy.zeros(oldspec.shape[1])
                     vec_snr = numpy.zeros(oldspec.shape[1])
                     for ind in range(oldspec.shape[1]):
                         spec = oldspec[:,ind]
                         aux = spec*fwindow
                         max_spec = aux.max()
                         m = list(aux).index(max_spec)
                         #Smooth
                         if (smooth == 0):   spec2 = spec
                         else:   spec2 = scipy.ndimage.filters.uniform_filter1d(spec,size=smooth)
                         #    Calculo de Momentos
                         bb = spec2[range(m,spec2.size)]
                         bb = (bb<n0).nonzero()
                         bb = bb[0]
                         ss = spec2[range(0,m + 1)]
                         ss = (ss<n0).nonzero()
                         ss = ss[0]
                         if (bb.size == 0):
                             bb0 = spec.size - 1 - m
                         else:
                             bb0 = bb[0] - 1
                             if (bb0 < 0):
                                 bb0 = 0
                         if (ss.size == 0):   ss1 = 1
                         else: ss1 = max(ss) + 1
                         if (ss1 > m):   ss1 = m
                         valid = numpy.asarray(range(int(m + bb0 - ss1 + 1))) + ss1
                         power = ((spec2[valid] - n0)*fwindow[valid]).sum()
                         fd = ((spec2[valid]- n0)*freq[valid]*fwindow[valid]).sum()/power
                         w = math.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum()/power)
                         snr = (spec2.mean()-n0)/n0
                         if (snr < 1.e-20) :
                             snr = 1.e-20
                         vec_power[ind] = power
                         vec_fd[ind] = fd
                         vec_w[ind] = w
                         vec_snr[ind] = snr
                     moments = numpy.vstack((vec_snr, vec_power, vec_fd, vec_w))
                     return moments
                 #------------------    Get SA Parameters    --------------------------
                 def GetSAParameters(self):
                     #SA en frecuencia
                     pairslist = self.dataOut.groupList
                     num_pairs = len(pairslist)
                     vel = self.dataOut.abscissaList
                     spectra = self.dataOut.data_pre
                     cspectra = self.dataIn.data_cspc
                     delta_v = vel[1] - vel[0]
                     #Calculating the power spectrum
                     spc_pow = numpy.sum(spectra, 3)*delta_v
                     #Normalizing Spectra
                     norm_spectra = spectra/spc_pow
                     #Calculating the norm_spectra at peak
                     max_spectra = numpy.max(norm_spectra, 3)
                     #Normalizing Cross Spectra
                     norm_cspectra = numpy.zeros(cspectra.shape)
                     for i in range(num_chan):
                         norm_cspectra[i,:,:] = cspectra[i,:,:]/numpy.sqrt(spc_pow[pairslist[i][0],:]*spc_pow[pairslist[i][1],:])
                     max_cspectra = numpy.max(norm_cspectra,2)
                     max_cspectra_index = numpy.argmax(norm_cspectra, 2)
                     for i in range(num_pairs):
                         cspc_par[i,:,:] = __calculateMoments(norm_cspectra)
                 #-------------------    Get Lags    ----------------------------------
             class SALags(Operation):
                 '''
                 Function GetMoments()
                 Input:
                     self.dataOut.data_pre
                     self.dataOut.abscissaList
                     self.dataOut.noise
                     self.dataOut.normFactor
                     self.dataOut.data_SNR
                     self.dataOut.groupList
                     self.dataOut.nChannels
                 Affected:
                     self.dataOut.data_param
                 '''
                 def run(self, dataOut):
                     data_acf = dataOut.data_pre[0]
                     data_ccf = dataOut.data_pre[1]
                     normFactor_acf = dataOut.normFactor[0]
                     normFactor_ccf = dataOut.normFactor[1]
                     pairs_acf = dataOut.groupList[0]
                     pairs_ccf = dataOut.groupList[1]
                     nHeights = dataOut.nHeights
                     absc = dataOut.abscissaList
                     noise = dataOut.noise
                     SNR = dataOut.data_SNR
                     nChannels = dataOut.nChannels
             #         pairsList = dataOut.groupList
             #         pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)
                     for l in range(len(pairs_acf)):
                         data_acf[l,:,:] = data_acf[l,:,:]/normFactor_acf[l,:]
                     for l in range(len(pairs_ccf)):
                         data_ccf[l,:,:] = data_ccf[l,:,:]/normFactor_ccf[l,:]
                     dataOut.data_param = numpy.zeros((len(pairs_ccf)*2 + 1, nHeights))
                     dataOut.data_param[:-1,:] = self.__calculateTaus(data_acf, data_ccf, absc)
                     dataOut.data_param[-1,:] = self.__calculateLag1Phase(data_acf, absc)
                     return
             #     def __getPairsAutoCorr(self, pairsList, nChannels):
             #
             #         pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
             #
             #         for l in range(len(pairsList)):
             #             firstChannel = pairsList[l][0]
             #             secondChannel = pairsList[l][1]
             #
             #             #Obteniendo pares de Autocorrelacion
             #             if firstChannel == secondChannel:
             #                 pairsAutoCorr[firstChannel] = int(l)
             #
             #         pairsAutoCorr = pairsAutoCorr.astype(int)
             #
             #         pairsCrossCorr = range(len(pairsList))
             #         pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
             #
             #         return pairsAutoCorr, pairsCrossCorr
                 def __calculateTaus(self, data_acf, data_ccf, lagRange):
                     lag0 = data_acf.shape[1]/2
                     #Funcion de Autocorrelacion
                     mean_acf = stats.nanmean(data_acf, axis = 0)
                     #Obtencion Indice de TauCross
                     ind_ccf = data_ccf.argmax(axis = 1)
                     #Obtencion Indice de TauAuto
                     ind_acf = numpy.zeros(ind_ccf.shape,dtype = 'int')
                     ccf_lag0 = data_ccf[:,lag0,:]
                     for i in range(ccf_lag0.shape[0]):
                         ind_acf[i,:] = numpy.abs(mean_acf - ccf_lag0[i,:]).argmin(axis = 0)
                     #Obtencion de TauCross y TauAuto
                     tau_ccf = lagRange[ind_ccf]
                     tau_acf  = lagRange[ind_acf]
                     Nan1, Nan2 = numpy.where(tau_ccf == lagRange[0])
                     tau_ccf[Nan1,Nan2] = numpy.nan
                     tau_acf[Nan1,Nan2] = numpy.nan
                     tau = numpy.vstack((tau_ccf,tau_acf))
                     return tau
                 def __calculateLag1Phase(self, data, lagTRange):
                     data1 = stats.nanmean(data, axis = 0)
                     lag1 = numpy.where(lagTRange == 0)[0][0] + 1
                     phase = numpy.angle(data1[lag1,:])
                     return phase
             class SpectralFitting(Operation):
                 '''
                     Function GetMoments()
                     Input:
                     Output:
                     Variables modified:
                 '''
                 def run(self, dataOut, getSNR = True, path=None, file=None, groupList=None):
                     if path != None:
                         sys.path.append(path)
                     self.dataOut.library = importlib.import_module(file)
                     #To be inserted as a parameter
                     groupArray = numpy.array(groupList)
             #         groupArray = numpy.array([[0,1],[2,3]])
                     self.dataOut.groupList = groupArray
                     nGroups = groupArray.shape[0]
                     nChannels = self.dataIn.nChannels
                     nHeights=self.dataIn.heightList.size
                     #Parameters Array
                     self.dataOut.data_param = None
                     #Set constants
                     constants = self.dataOut.library.setConstants(self.dataIn)
                     self.dataOut.constants = constants
                     M = self.dataIn.normFactor
                     N = self.dataIn.nFFTPoints
                     ippSeconds = self.dataIn.ippSeconds
                     K = self.dataIn.nIncohInt
                     pairsArray = numpy.array(self.dataIn.pairsList)
                     #List of possible combinations
                     listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)
                     indCross = numpy.zeros(len(list(listComb)), dtype = 'int')
                     if getSNR:
                         listChannels = groupArray.reshape((groupArray.size))
                         listChannels.sort()
                         noise = self.dataIn.getNoise()
                         self.dataOut.data_SNR = self.__getSNR(self.dataIn.data_spc[listChannels,:,:], noise[listChannels])
                     for i in range(nGroups):
                         coord = groupArray[i,:]
                         #Input data array
                         data = self.dataIn.data_spc[coord,:,:]/(M*N)
                         data = data.reshape((data.shape[0]*data.shape[1],data.shape[2]))
                         #Cross Spectra data array for Covariance Matrixes
                         ind = 0
                         for pairs in listComb:
                             pairsSel = numpy.array([coord[x],coord[y]])
                             indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0])
                             ind += 1
                         dataCross = self.dataIn.data_cspc[indCross,:,:]/(M*N)
                         dataCross = dataCross**2/K
                         for h in range(nHeights):
             #                 print self.dataOut.heightList[h]
                             #Input
                             d = data[:,h]
                             #Covariance Matrix
                             D = numpy.diag(d**2/K)
                             ind = 0
                             for pairs in listComb:
                                 #Coordinates in Covariance Matrix
                                 x = pairs[0]
                                 y = pairs[1]
                                 #Channel Index
                                 S12 = dataCross[ind,:,h]
                                 D12 = numpy.diag(S12)
                                 #Completing Covariance Matrix with Cross Spectras
                                 D[x*N:(x+1)*N,y*N:(y+1)*N] = D12
                                 D[y*N:(y+1)*N,x*N:(x+1)*N] = D12
                                 ind += 1
                             Dinv=numpy.linalg.inv(D)
                             L=numpy.linalg.cholesky(Dinv)
                             LT=L.T
                             dp = numpy.dot(LT,d)
                             #Initial values
                             data_spc = self.dataIn.data_spc[coord,:,h]
                             if (h>0)and(error1[3]<5):
                                 p0 = self.dataOut.data_param[i,:,h-1]
                             else:
                                 p0 = numpy.array(self.dataOut.library.initialValuesFunction(data_spc, constants, i))
                             try:
                                 #Least Squares
                                 minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True)
             #                   minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))
                                 #Chi square error
                                 error0 = numpy.sum(infodict['fvec']**2)/(2*N)
                                 #Error with Jacobian
                                 error1 = self.dataOut.library.errorFunction(minp,constants,LT)
                             except:
                                 minp = p0*numpy.nan
                                 error0 = numpy.nan
                                 error1 = p0*numpy.nan
                             #Save
                             if self.dataOut.data_param == None:
                                 self.dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan
                                 self.dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan
                             self.dataOut.data_error[i,:,h] = numpy.hstack((error0,error1))
                             self.dataOut.data_param[i,:,h] = minp
                     return
                 def __residFunction(self, p, dp, LT, constants):
                     fm = self.dataOut.library.modelFunction(p, constants)
                     fmp=numpy.dot(LT,fm)
                     return  dp-fmp
                 def __getSNR(self, z, noise):
                     avg = numpy.average(z, axis=1)
                     SNR = (avg.T-noise)/noise
                     SNR = SNR.T
                     return SNR
                 def __chisq(p,chindex,hindex):
                     #similar to Resid but calculates CHI**2
                     [LT,d,fm]=setupLTdfm(p,chindex,hindex)
                     dp=numpy.dot(LT,d)
                     fmp=numpy.dot(LT,fm)
                     chisq=numpy.dot((dp-fmp).T,(dp-fmp))
                     return chisq
             class WindProfiler(Operation):
                 __isConfig = False
                 __initime = None
                 __lastdatatime = None
                 __integrationtime = None
                 __buffer = None
                 __dataReady = False
                 __firstdata = None
                 n = None
                 def __init__(self):
                     Operation.__init__(self)
                 def __calculateCosDir(self, elev, azim):
                     zen = (90 - elev)*numpy.pi/180
                     azim = azim*numpy.pi/180
                     cosDirX = numpy.sqrt((1-numpy.cos(zen)**2)/((1+numpy.tan(azim)**2)))
                     cosDirY = numpy.sqrt(1-numpy.cos(zen)**2-cosDirX**2)
                     signX = numpy.sign(numpy.cos(azim))
                     signY = numpy.sign(numpy.sin(azim))
                     cosDirX = numpy.copysign(cosDirX, signX)
                     cosDirY = numpy.copysign(cosDirY, signY)
                     return cosDirX, cosDirY
                 def __calculateAngles(self, theta_x, theta_y, azimuth):
                     dir_cosw = numpy.sqrt(1-theta_x**2-theta_y**2)
                     zenith_arr = numpy.arccos(dir_cosw)
                     azimuth_arr = numpy.arctan2(theta_x,theta_y) + azimuth*math.pi/180
                     dir_cosu = numpy.sin(azimuth_arr)*numpy.sin(zenith_arr)
                     dir_cosv = numpy.cos(azimuth_arr)*numpy.sin(zenith_arr)
                     return azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw
                 def __calculateMatA(self, dir_cosu, dir_cosv, dir_cosw, horOnly):
             #
                     if horOnly:
                         A = numpy.c_[dir_cosu,dir_cosv]
                     else:
                         A = numpy.c_[dir_cosu,dir_cosv,dir_cosw]
                     A = numpy.asmatrix(A)
                     A1 = numpy.linalg.inv(A.transpose()*A)*A.transpose()
                     return A1
                 def __correctValues(self, heiRang, phi, velRadial, SNR):
                     listPhi = phi.tolist()
                     maxid = listPhi.index(max(listPhi))
                     minid = listPhi.index(min(listPhi))
                     rango = range(len(phi))
                #     rango = numpy.delete(rango,maxid)
                     heiRang1 = heiRang*math.cos(phi[maxid])
                     heiRangAux = heiRang*math.cos(phi[minid])
                     indOut = (heiRang1 < heiRangAux[0]).nonzero()
                     heiRang1 = numpy.delete(heiRang1,indOut)
                     velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
                     SNR1 = numpy.zeros([len(phi),len(heiRang1)])
                     for i in rango:
                         x = heiRang*math.cos(phi[i])
                         y1 = velRadial[i,:]
                         f1 = interpolate.interp1d(x,y1,kind = 'cubic')
                         x1 = heiRang1
                         y11 = f1(x1)
                         y2 = SNR[i,:]
                         f2 = interpolate.interp1d(x,y2,kind = 'cubic')
                         y21 = f2(x1)
                         velRadial1[i,:] = y11
                         SNR1[i,:] = y21
                     return heiRang1, velRadial1, SNR1
                 def __calculateVelUVW(self, A, velRadial):
                     #Operacion Matricial
             #         velUVW = numpy.zeros((velRadial.shape[1],3))
             #         for ind in range(velRadial.shape[1]):
             #             velUVW[ind,:] = numpy.dot(A,velRadial[:,ind])
             #         velUVW = velUVW.transpose()
                     velUVW = numpy.zeros((A.shape[0],velRadial.shape[1]))
                     velUVW[:,:] = numpy.dot(A,velRadial)
                     return velUVW
             #     def techniqueDBS(self, velRadial0, dirCosx, disrCosy, azimuth, correct, horizontalOnly, heiRang, SNR0):
                 def techniqueDBS(self, kwargs):
                     """
                     Function that implements Doppler Beam Swinging (DBS) technique.
                     Input:    Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
                                 Direction correction (if necessary), Ranges and SNR
                     Output:    Winds estimation (Zonal, Meridional and Vertical)
                     Parameters affected:    Winds, height range, SNR
                     """
                     velRadial0 = kwargs['velRadial']
                     heiRang = kwargs['heightList']
                     SNR0 = kwargs['SNR']
                     if kwargs.has_key('dirCosx') and kwargs.has_key('dirCosy'):
                         theta_x = numpy.array(kwargs['dirCosx'])
                         theta_y = numpy.array(kwargs['dirCosy'])
                     else:
                         elev = numpy.array(kwargs['elevation'])
                         azim = numpy.array(kwargs['azimuth'])
                         theta_x, theta_y = self.__calculateCosDir(elev, azim)
                     azimuth = kwargs['correctAzimuth']
                     if kwargs.has_key('horizontalOnly'):
                         horizontalOnly = kwargs['horizontalOnly']
                     else:   horizontalOnly = False
                     if kwargs.has_key('correctFactor'):
                         correctFactor = kwargs['correctFactor']
                     else:   correctFactor = 1
                     if kwargs.has_key('channelList'):
                         channelList = kwargs['channelList']
                         if len(channelList) == 2:
                             horizontalOnly = True
                         arrayChannel = numpy.array(channelList)
                         param = param[arrayChannel,:,:]
                         theta_x = theta_x[arrayChannel]
                         theta_y = theta_y[arrayChannel]
                     azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(theta_x, theta_y, azimuth)
                     heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correctFactor*velRadial0, SNR0)
                     A = self.__calculateMatA(dir_cosu, dir_cosv, dir_cosw, horizontalOnly)
                     #Calculo de Componentes de la velocidad con DBS
                     winds = self.__calculateVelUVW(A,velRadial1)
                     return winds, heiRang1, SNR1
                 def __calculateDistance(self, posx, posy, pairs_ccf, azimuth = None):
                     nPairs = len(pairs_ccf)
                     posx = numpy.asarray(posx)
                     posy = numpy.asarray(posy)
                     #Rotacion Inversa para alinear con el azimuth
                     if azimuth!= None:
                         azimuth = azimuth*math.pi/180
                         posx1 = posx*math.cos(azimuth) + posy*math.sin(azimuth)
                         posy1 = -posx*math.sin(azimuth) + posy*math.cos(azimuth)
                     else:
                         posx1 = posx
                         posy1 = posy
                     #Calculo de Distancias
                     distx = numpy.zeros(nPairs)
                     disty = numpy.zeros(nPairs)
                     dist = numpy.zeros(nPairs)
                     ang = numpy.zeros(nPairs)
                     for i in range(nPairs):
                         distx[i] = posx1[pairs_ccf[i][1]] - posx1[pairs_ccf[i][0]]
                         disty[i] = posy1[pairs_ccf[i][1]] - posy1[pairs_ccf[i][0]]
                         dist[i] = numpy.sqrt(distx[i]**2 + disty[i]**2)
                         ang[i] = numpy.arctan2(disty[i],distx[i])
                     return distx, disty, dist, ang
                     #Calculo de Matrices
             #         nPairs = len(pairs)
             #         ang1 = numpy.zeros((nPairs, 2, 1))
             #         dist1 = numpy.zeros((nPairs, 2, 1))
             #
             #         for j in range(nPairs):
             #             dist1[j,0,0] = dist[pairs[j][0]]
             #             dist1[j,1,0] = dist[pairs[j][1]]
             #             ang1[j,0,0] = ang[pairs[j][0]]
             #             ang1[j,1,0] = ang[pairs[j][1]]
             #
             #         return distx,disty, dist1,ang1
                 def __calculateVelVer(self, phase, lagTRange, _lambda):
                     Ts = lagTRange[1] - lagTRange[0]
                     velW = -_lambda*phase/(4*math.pi*Ts)
                     return velW
                 def __calculateVelHorDir(self, dist, tau1, tau2, ang):
                     nPairs = tau1.shape[0]
                     nHeights = tau1.shape[1]
                     vel = numpy.zeros((nPairs,3,nHeights))
                     dist1 = numpy.reshape(dist, (dist.size,1))
                     angCos = numpy.cos(ang)
                     angSin = numpy.sin(ang)
                     vel0 = dist1*tau1/(2*tau2**2)
                     vel[:,0,:] = (vel0*angCos).sum(axis = 1)
                     vel[:,1,:] = (vel0*angSin).sum(axis = 1)
                     ind = numpy.where(numpy.isinf(vel))
                     vel[ind] = numpy.nan
                     return vel
             #     def __getPairsAutoCorr(self, pairsList, nChannels):
             #
             #         pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
             #
             #         for l in range(len(pairsList)):
             #             firstChannel = pairsList[l][0]
             #             secondChannel = pairsList[l][1]
             #
             #             #Obteniendo pares de Autocorrelacion
             #             if firstChannel == secondChannel:
             #                 pairsAutoCorr[firstChannel] = int(l)
             #
             #         pairsAutoCorr = pairsAutoCorr.astype(int)
             #
             #         pairsCrossCorr = range(len(pairsList))
             #         pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
             #
             #         return pairsAutoCorr, pairsCrossCorr
             #     def techniqueSA(self, pairsSelected, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, lagTRange, correctFactor):
                 def techniqueSA(self, kwargs):
                     """
                     Function that implements Spaced Antenna (SA) technique.
                     Input:    Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
                                 Direction correction (if necessary), Ranges and SNR
                     Output:    Winds estimation (Zonal, Meridional and Vertical)
                     Parameters affected:    Winds
                     """
                     position_x = kwargs['positionX']
                     position_y = kwargs['positionY']
                     azimuth = kwargs['azimuth']
                     if kwargs.has_key('correctFactor'):
                         correctFactor = kwargs['correctFactor']
                     else:
                         correctFactor = 1
                     groupList = kwargs['groupList']
                     pairs_ccf = groupList[1]
                     tau = kwargs['tau']
                     _lambda = kwargs['_lambda']
                     #Cross Correlation pairs obtained
             #         pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairssList, nChannels)
             #         pairsArray = numpy.array(pairsList)[pairsCrossCorr]
             #         pairsSelArray = numpy.array(pairsSelected)
             #         pairs = []
             #
             #         #Wind estimation pairs obtained
             #         for i in range(pairsSelArray.shape[0]/2):
             #             ind1 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i], axis = 1))[0][0]
             #             ind2 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i + 1], axis = 1))[0][0]
             #             pairs.append((ind1,ind2))
                     indtau = tau.shape[0]/2
                     tau1 = tau[:indtau,:]
                     tau2 = tau[indtau:-1,:]
             #         tau1 = tau1[pairs,:]
             #         tau2 = tau2[pairs,:]
                     phase1 = tau[-1,:]
                     #---------------------------------------------------------------------
                     #Metodo Directo
                     distx, disty, dist, ang = self.__calculateDistance(position_x, position_y, pairs_ccf,azimuth)
                     winds = self.__calculateVelHorDir(dist, tau1, tau2, ang)
                     winds = stats.nanmean(winds, axis=0)
                     #---------------------------------------------------------------------
                     #Metodo General
             #         distx, disty, dist = self.calculateDistance(position_x,position_y,pairsCrossCorr, pairsList, azimuth)
             #         #Calculo Coeficientes de Funcion de Correlacion
             #         F,G,A,B,H = self.calculateCoef(tau1,tau2,distx,disty,n)
             #         #Calculo de Velocidades
             #         winds = self.calculateVelUV(F,G,A,B,H)
                     #---------------------------------------------------------------------
                     winds[2,:] = self.__calculateVelVer(phase1, lagTRange, _lambda)
                     winds = correctFactor*winds
                     return winds
                 def __checkTime(self, currentTime, paramInterval, outputInterval):
                     dataTime = currentTime + paramInterval
                     deltaTime = dataTime - self.__initime
                     if deltaTime >= outputInterval or deltaTime < 0:
                         self.__dataReady = True
                     return
                 def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax):
                     '''
                     Function that implements winds estimation technique with detected meteors.
                     Input:    Detected meteors, Minimum meteor quantity to wind estimation
                     Output:    Winds estimation (Zonal and Meridional)
                     Parameters affected:    Winds
                     '''
             #         print arrayMeteor.shape
                     #Settings
                     nInt = (heightMax - heightMin)/2
             #         print nInt
                     nInt = int(nInt)
             #         print nInt
                     winds = numpy.zeros((2,nInt))*numpy.nan
                     #Filter errors
                     error = numpy.where(arrayMeteor[:,-1] == 0)[0]
                     finalMeteor = arrayMeteor[error,:]
                     #Meteor Histogram
                     finalHeights = finalMeteor[:,2]
                     hist = numpy.histogram(finalHeights, bins = nInt, range = (heightMin,heightMax))
                     nMeteorsPerI = hist[0]
                     heightPerI = hist[1]
                     #Sort of meteors
                     indSort = finalHeights.argsort()
                     finalMeteor2 = finalMeteor[indSort,:]
                     #    Calculating winds
                     ind1 = 0
                     ind2 = 0
                     for i in range(nInt):
                         nMet = nMeteorsPerI[i]
                         ind1 = ind2
                         ind2 = ind1 + nMet
                         meteorAux = finalMeteor2[ind1:ind2,:]
                         if meteorAux.shape[0] >= meteorThresh:
                             vel = meteorAux[:, 6]
                             zen = meteorAux[:, 4]*numpy.pi/180
                             azim = meteorAux[:, 3]*numpy.pi/180
                             n = numpy.cos(zen)
                     #         m = (1 - n**2)/(1 - numpy.tan(azim)**2)
                     #         l = m*numpy.tan(azim)
                             l = numpy.sin(zen)*numpy.sin(azim)
                             m = numpy.sin(zen)*numpy.cos(azim)
                             A = numpy.vstack((l, m)).transpose()
                             A1 = numpy.dot(numpy.linalg.inv( numpy.dot(A.transpose(),A) ),A.transpose())
                             windsAux = numpy.dot(A1, vel)
                             winds[0,i] = windsAux[0]
                             winds[1,i] = windsAux[1]
                     return winds, heightPerI[:-1]
                 def techniqueNSM_SA(self, **kwargs):
                     metArray = kwargs['metArray']
                     heightList = kwargs['heightList']
                     timeList = kwargs['timeList']
                     rx_location = kwargs['rx_location']
                     groupList = kwargs['groupList']
                     azimuth = kwargs['azimuth']
                     dfactor = kwargs['dfactor']
                     k = kwargs['k']
                     azimuth1, dist = self.__calculateAzimuth1(rx_location, groupList, azimuth)
                     d = dist*dfactor
                     #Phase calculation
                     metArray1 = self.__getPhaseSlope(metArray, heightList, timeList)
                     metArray1[:,-2] = metArray1[:,-2]*metArray1[:,2]*1000/(k*d[metArray1[:,1].astype(int)]) #angles into velocities
                     velEst = numpy.zeros((heightList.size,2))*numpy.nan
                     azimuth1 = azimuth1*numpy.pi/180
                     for i in range(heightList.size):
                         h = heightList[i]
                         indH = numpy.where((metArray1[:,2] == h)&(numpy.abs(metArray1[:,-2]) < 100))[0]
                         metHeight = metArray1[indH,:]
                         if metHeight.shape[0] >= 2:
                             velAux = numpy.asmatrix(metHeight[:,-2]).T    #Radial Velocities
                             iazim = metHeight[:,1].astype(int)
                             azimAux = numpy.asmatrix(azimuth1[iazim]).T    #Azimuths
                             A = numpy.hstack((numpy.cos(azimAux),numpy.sin(azimAux)))
                             A = numpy.asmatrix(A)
                             A1 = numpy.linalg.pinv(A.transpose()*A)*A.transpose()
                             velHor = numpy.dot(A1,velAux)
                             velEst[i,:] = numpy.squeeze(velHor)
                     return velEst
                 def __getPhaseSlope(self, metArray, heightList, timeList):
                     meteorList = []
                     #utctime sec1 height SNR velRad ph0 ph1 ph2 coh0 coh1 coh2
                     #Putting back together the meteor matrix
                     utctime = metArray[:,0]
                     uniqueTime = numpy.unique(utctime)
                     phaseDerThresh = 0.5
                     ippSeconds = timeList[1] - timeList[0]
                     sec = numpy.where(timeList>1)[0][0]
                     nPairs = metArray.shape[1] - 6
                     nHeights = len(heightList)
                     for t in uniqueTime:
                         metArray1 = metArray[utctime==t,:]
             #         phaseDerThresh = numpy.pi/4 #reducir Phase thresh
                         tmet = metArray1[:,1].astype(int)
                         hmet = metArray1[:,2].astype(int)
                         metPhase = numpy.zeros((nPairs, heightList.size, timeList.size - 1))
                         metPhase[:,:] = numpy.nan
                         metPhase[:,hmet,tmet] = metArray1[:,6:].T
                         #Delete short trails
                         metBool = ~numpy.isnan(metPhase[0,:,:])
                         heightVect = numpy.sum(metBool, axis = 1)
                         metBool[heightVect<sec,:] = False
                         metPhase[:,heightVect<sec,:] = numpy.nan
                         #Derivative
                         metDer = numpy.abs(metPhase[:,:,1:] - metPhase[:,:,:-1])
                         phDerAux = numpy.dstack((numpy.full((nPairs,nHeights,1), False, dtype=bool),metDer > phaseDerThresh))
                         metPhase[phDerAux] = numpy.nan
                         #--------------------------METEOR DETECTION    -----------------------------------------
                         indMet = numpy.where(numpy.any(metBool,axis=1))[0]
                         for p in numpy.arange(nPairs):
                             phase = metPhase[p,:,:]
                             phDer = metDer[p,:,:]
                             for h in indMet:
                                 height = heightList[h]
                                 phase1 = phase[h,:] #82
                                 phDer1 = phDer[h,:]
                                 phase1[~numpy.isnan(phase1)] = numpy.unwrap(phase1[~numpy.isnan(phase1)])   #Unwrap
                                 indValid = numpy.where(~numpy.isnan(phase1))[0]
                                 initMet = indValid[0]
                                 endMet = 0
                                 for i in range(len(indValid)-1):
                                     #Time difference
                                     inow = indValid[i]
                                     inext = indValid[i+1]
                                     idiff = inext - inow
                                     #Phase difference
                                     phDiff = numpy.abs(phase1[inext] - phase1[inow])
                                     if idiff>sec or phDiff>numpy.pi/4 or inext==indValid[-1]:   #End of Meteor
                                         sizeTrail = inow - initMet + 1
                                         if sizeTrail>3*sec:  #Too short meteors
                                             x = numpy.arange(initMet,inow+1)*ippSeconds
                                             y = phase1[initMet:inow+1]
                                             ynnan = ~numpy.isnan(y)
                                             x = x[ynnan]
                                             y = y[ynnan]
                                             slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
                                             ylin = x*slope + intercept
                                             rsq = r_value**2
                                             if rsq > 0.5:
                                                 vel = slope#*height*1000/(k*d)
                                                 estAux = numpy.array([utctime,p,height, vel, rsq])
                                                 meteorList.append(estAux)
                                         initMet = inext
                     metArray2 = numpy.array(meteorList)
                     return metArray2
                 def __calculateAzimuth1(self, rx_location, pairslist, azimuth0):
                     azimuth1 = numpy.zeros(len(pairslist))
                     dist = numpy.zeros(len(pairslist))
                     for i in range(len(rx_location)):
                         ch0 = pairslist[i][0]
                         ch1 = pairslist[i][1]
                         diffX = rx_location[ch0][0] - rx_location[ch1][0]
                         diffY = rx_location[ch0][1] - rx_location[ch1][1]
                         azimuth1[i] = numpy.arctan2(diffY,diffX)*180/numpy.pi
                         dist[i] = numpy.sqrt(diffX**2 + diffY**2)
                     azimuth1 -= azimuth0
                     return azimuth1, dist
                 def techniqueNSM_DBS(self, **kwargs):
                     metArray = kwargs['metArray']
                     heightList = kwargs['heightList']
                     timeList = kwargs['timeList']
                     zenithList = kwargs['zenithList']
                     nChan = numpy.max(cmet) + 1
                     nHeights = len(heightList)
                     utctime = metArray[:,0]
                     cmet = metArray[:,1]
                     hmet = metArray1[:,3].astype(int)
                     h1met = heightList[hmet]*zenithList[cmet]
                     vmet = metArray1[:,5]
                     for i in range(nHeights - 1):
                         hmin = heightList[i]
                         hmax = heightList[i + 1]
                         vthisH = vmet[(h1met>=hmin) & (h1met<hmax)]
                     return data_output
                 def run(self, dataOut, technique, **kwargs):
                     param = dataOut.data_param
                     if dataOut.abscissaList != None:
                         absc = dataOut.abscissaList[:-1]
                     noise = dataOut.noise
                     heightList = dataOut.heightList
                     SNR = dataOut.data_SNR
                     if technique == 'DBS':
                         kwargs['velRadial'] = param[:,1,:] #Radial velocity
                         kwargs['heightList'] = heightList
                         kwargs['SNR'] = SNR
                         dataOut.data_output, dataOut.heightList, dataOut.data_SNR = self.techniqueDBS(kwargs) #DBS Function
                         dataOut.utctimeInit = dataOut.utctime
                         dataOut.outputInterval = dataOut.paramInterval
                     elif technique == 'SA':
                         #Parameters
             #             position_x = kwargs['positionX']
             #             position_y = kwargs['positionY']
             #             azimuth = kwargs['azimuth']
             #
             #             if kwargs.has_key('crosspairsList'):
             #                 pairs = kwargs['crosspairsList']
             #             else:
             #                 pairs = None
             #
             #             if kwargs.has_key('correctFactor'):
             #                 correctFactor = kwargs['correctFactor']
             #             else:
             #                 correctFactor = 1
             #             tau = dataOut.data_param
             #             _lambda = dataOut.C/dataOut.frequency
             #             pairsList = dataOut.groupList
             #             nChannels = dataOut.nChannels
                         kwargs['groupList'] = dataOut.groupList
                         kwargs['tau'] = dataOut.data_param
                         kwargs['_lambda'] = dataOut.C/dataOut.frequency
             #             dataOut.data_output = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor)
                         dataOut.data_output = self.techniqueSA(kwargs)
                         dataOut.utctimeInit = dataOut.utctime
                         dataOut.outputInterval = dataOut.timeInterval
                     elif technique == 'Meteors':
                         dataOut.flagNoData = True
                         self.__dataReady = False
                         if kwargs.has_key('nHours'):
                             nHours = kwargs['nHours']
                         else:
                             nHours = 1
                         if kwargs.has_key('meteorsPerBin'):
                             meteorThresh = kwargs['meteorsPerBin']
                         else:
                             meteorThresh = 6
                         if kwargs.has_key('hmin'):
                             hmin = kwargs['hmin']
                         else:   hmin = 70
                         if kwargs.has_key('hmax'):
                             hmax = kwargs['hmax']
                         else:   hmax = 110
                         dataOut.outputInterval = nHours*3600
                         if self.__isConfig == False:
             #                 self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
                             #Get Initial LTC time
                             self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
                             self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
                             self.__isConfig = True
                         if self.__buffer == None:
                             self.__buffer = dataOut.data_param
                             self.__firstdata = copy.copy(dataOut)
                         else:
                             self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
                         self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
                         if self.__dataReady:
                             dataOut.utctimeInit = self.__initime
                             self.__initime += dataOut.outputInterval #to erase time offset
                             dataOut.data_output, dataOut.heightList = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax)
                             dataOut.flagNoData = False
                             self.__buffer = None
                     elif technique == 'Meteors1':
                         dataOut.flagNoData = True
                         self.__dataReady = False
                         if kwargs.has_key('nMins'):
                             nMins = kwargs['nMins']
                         else: nMins = 20
                         if kwargs.has_key('rx_location'):
                             rx_location = kwargs['rx_location']
                         else: rx_location = [(0,1),(1,1),(1,0)]
                         if kwargs.has_key('azimuth'):
                             azimuth = kwargs['azimuth']
                         else: azimuth = 51
                         if kwargs.has_key('dfactor'):
                             dfactor = kwargs['dfactor']
                         if kwargs.has_key('mode'):
                             mode = kwargs['mode']
                         else: mode = 'SA'
                         #Borrar luego esto
                         if dataOut.groupList == None:
                             dataOut.groupList = [(0,1),(0,2),(1,2)]
                         groupList = dataOut.groupList
                         C = 3e8
                         freq = 50e6
                         lamb = C/freq
                         k = 2*numpy.pi/lamb
                         timeList = dataOut.abscissaList
                         heightList = dataOut.heightList
                         if self.__isConfig == False:
                             dataOut.outputInterval = nMins*60
             #                 self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
                             #Get Initial LTC time
                             initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
                             minuteAux = initime.minute
                             minuteNew = int(numpy.floor(minuteAux/nMins)*nMins)
                             self.__initime = (initime.replace(minute = minuteNew, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
                             self.__isConfig = True
                         if self.__buffer == None:
                             self.__buffer = dataOut.data_param
                             self.__firstdata = copy.copy(dataOut)
                         else:
                             self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
                         self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
                         if self.__dataReady:
                             dataOut.utctimeInit = self.__initime
                             self.__initime += dataOut.outputInterval #to erase time offset
                             metArray = self.__buffer
                             if mode == 'SA':
                                 dataOut.data_output = self.techniqueNSM_SA(rx_location=rx_location, groupList=groupList, azimuth=azimuth, dfactor=dfactor, k=k,metArray=metArray, heightList=heightList,timeList=timeList)
                             elif mode == 'DBS':
                                 dataOut.data_output = self.techniqueNSM_DBS(metArray=metArray,heightList=heightList,timeList=timeList)
                             dataOut.data_output = dataOut.data_output.T
                             dataOut.flagNoData = False
                             self.__buffer = None
                     return
             class EWDriftsEstimation(Operation):
                 def __init__(self):
                     Operation.__init__(self)
                 def __correctValues(self, heiRang, phi, velRadial, SNR):
                     listPhi = phi.tolist()
                     maxid = listPhi.index(max(listPhi))
                     minid = listPhi.index(min(listPhi))
                     rango = range(len(phi))
                #     rango = numpy.delete(rango,maxid)
                     heiRang1 = heiRang*math.cos(phi[maxid])
                     heiRangAux = heiRang*math.cos(phi[minid])
                     indOut = (heiRang1 < heiRangAux[0]).nonzero()
                     heiRang1 = numpy.delete(heiRang1,indOut)
                     velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
                     SNR1 = numpy.zeros([len(phi),len(heiRang1)])
                     for i in rango:
                         x = heiRang*math.cos(phi[i])
                         y1 = velRadial[i,:]
                         f1 = interpolate.interp1d(x,y1,kind = 'cubic')
                         x1 = heiRang1
                         y11 = f1(x1)
                         y2 = SNR[i,:]
                         f2 = interpolate.interp1d(x,y2,kind = 'cubic')
                         y21 = f2(x1)
                         velRadial1[i,:] = y11
                         SNR1[i,:] = y21
                     return heiRang1, velRadial1, SNR1
                 def run(self, dataOut, zenith, zenithCorrection):
                     heiRang = dataOut.heightList
                     velRadial = dataOut.data_param[:,3,:]
                     SNR = dataOut.data_SNR
                     zenith = numpy.array(zenith)
                     zenith -= zenithCorrection
                     zenith *= numpy.pi/180
                     heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR)
                     alp = zenith[0]
                     bet = zenith[1]
                     w_w = velRadial1[0,:]
                     w_e = velRadial1[1,:]
                     w = (w_w*numpy.sin(bet) - w_e*numpy.sin(alp))/(numpy.cos(alp)*numpy.sin(bet) - numpy.cos(bet)*numpy.sin(alp))
                     u = (w_w*numpy.cos(bet) - w_e*numpy.cos(alp))/(numpy.sin(alp)*numpy.cos(bet) - numpy.sin(bet)*numpy.cos(alp))
                     winds = numpy.vstack((u,w))
                     dataOut.heightList = heiRang1
                     dataOut.data_output = winds
                     dataOut.data_SNR = SNR1
                     dataOut.utctimeInit = dataOut.utctime
                     dataOut.outputInterval = dataOut.timeInterval
                     return
             #---------------    Non Specular Meteor    ----------------
             class NonSpecularMeteorDetection(Operation):
                 def run(self, mode, SNRthresh=8, phaseDerThresh=0.5, cohThresh=0.8, allData = False):
                     data_acf = self.dataOut.data_pre[0]
                     data_ccf = self.dataOut.data_pre[1]
                     lamb = self.dataOut.C/self.dataOut.frequency
                     tSamp = self.dataOut.ippSeconds*self.dataOut.nCohInt
                     paramInterval = self.dataOut.paramInterval
                     nChannels = data_acf.shape[0]
                     nLags = data_acf.shape[1]
                     nProfiles = data_acf.shape[2]
                     nHeights = self.dataOut.nHeights
                     nCohInt = self.dataOut.nCohInt
                     sec = numpy.round(nProfiles/self.dataOut.paramInterval)
                     heightList = self.dataOut.heightList
                     ippSeconds = self.dataOut.ippSeconds*self.dataOut.nCohInt*self.dataOut.nAvg
                     utctime = self.dataOut.utctime
                     self.dataOut.abscissaList = numpy.arange(0,paramInterval+ippSeconds,ippSeconds)
                     #------------------------    SNR    --------------------------------------
                     power = data_acf[:,0,:,:].real
                     noise = numpy.zeros(nChannels)
                     SNR = numpy.zeros(power.shape)
                     for i in range(nChannels):
                         noise[i] = hildebrand_sekhon(power[i,:], nCohInt)
                         SNR[i] = (power[i]-noise[i])/noise[i]
                     SNRm = numpy.nanmean(SNR, axis = 0)
                     SNRdB = 10*numpy.log10(SNR)
                     if mode == 'SA':
                         nPairs = data_ccf.shape[0]
                         #----------------------    Coherence and Phase   --------------------------
                         phase = numpy.zeros(data_ccf[:,0,:,:].shape)
             #             phase1 = numpy.copy(phase)
                         coh1 = numpy.zeros(data_ccf[:,0,:,:].shape)
                         for p in range(nPairs):
                             ch0 = self.dataOut.groupList[p][0]
                             ch1 = self.dataOut.groupList[p][1]
                             ccf = data_ccf[p,0,:,:]/numpy.sqrt(data_acf[ch0,0,:,:]*data_acf[ch1,0,:,:])
                             phase[p,:,:] = ndimage.median_filter(numpy.angle(ccf), size = (5,1)) #median filter
             #                 phase1[p,:,:] = numpy.angle(ccf) #median filter
                             coh1[p,:,:] = ndimage.median_filter(numpy.abs(ccf), 5) #median filter
             #                 coh1[p,:,:] = numpy.abs(ccf) #median filter
                         coh = numpy.nanmax(coh1, axis = 0)
             #             struc = numpy.ones((5,1))
             #             coh = ndimage.morphology.grey_dilation(coh, size=(10,1))
                         #----------------------    Radial Velocity    ----------------------------
                         phaseAux = numpy.mean(numpy.angle(data_acf[:,1,:,:]), axis = 0)
                         velRad = phaseAux*lamb/(4*numpy.pi*tSamp)
                         if allData:
                             boolMetFin = ~numpy.isnan(SNRm)
             #                 coh[:-1,:] = numpy.nanmean(numpy.abs(phase[:,1:,:] - phase[:,:-1,:]),axis=0)
                         else:
                             #------------------------    Meteor mask    ---------------------------------
             #                 #SNR mask
             #                 boolMet = (SNRdB>SNRthresh)#|(~numpy.isnan(SNRdB))
             #
             #                 #Erase small objects
             #                 boolMet1 = self.__erase_small(boolMet, 2*sec, 5)
             #
             #                 auxEEJ = numpy.sum(boolMet1,axis=0)
             #                 indOver = auxEEJ>nProfiles*0.8  #Use this later
             #                 indEEJ = numpy.where(indOver)[0]
             #                 indNEEJ = numpy.where(~indOver)[0]
             #
             #                 boolMetFin = boolMet1
             #
             #                 if indEEJ.size > 0:
             #                     boolMet1[:,indEEJ] = False  #Erase heights with EEJ
             #
             #                     boolMet2 = coh > cohThresh
             #                     boolMet2 = self.__erase_small(boolMet2, 2*sec,5)
             #
             #                     #Final Meteor mask
             #                     boolMetFin = boolMet1|boolMet2
                             #Coherence mask
                             boolMet1 = coh > 0.75
                             struc = numpy.ones((30,1))
                             boolMet1 = ndimage.morphology.binary_dilation(boolMet1, structure=struc)
                             #Derivative mask
                             derPhase = numpy.nanmean(numpy.abs(phase[:,1:,:] - phase[:,:-1,:]),axis=0)
                             boolMet2 = derPhase < 0.2
             #                 boolMet2 = ndimage.morphology.binary_opening(boolMet2)
             #                 boolMet2 = ndimage.morphology.binary_closing(boolMet2, structure = numpy.ones((10,1)))
                             boolMet2 = ndimage.median_filter(boolMet2,size=5)
                             boolMet2 = numpy.vstack((boolMet2,numpy.full((1,nHeights), True, dtype=bool)))
             #                 #Final mask
             #                 boolMetFin = boolMet2
                             boolMetFin = boolMet1&boolMet2
             #                 boolMetFin = ndimage.morphology.binary_dilation(boolMetFin)
                         #Creating data_param
                         coordMet = numpy.where(boolMetFin)
                         tmet = coordMet[0]
                         hmet = coordMet[1]
                         data_param = numpy.zeros((tmet.size, 6 + nPairs))
                         data_param[:,0] = utctime
                         data_param[:,1] = tmet
                         data_param[:,2] = hmet
                         data_param[:,3] = SNRm[tmet,hmet]
                         data_param[:,4] = velRad[tmet,hmet]
                         data_param[:,5] = coh[tmet,hmet]
                         data_param[:,6:] = phase[:,tmet,hmet].T
                     elif mode == 'DBS':
                         self.dataOut.groupList = numpy.arange(nChannels)
                         #Radial Velocities
             #             phase = numpy.angle(data_acf[:,1,:,:])
                         phase = ndimage.median_filter(numpy.angle(data_acf[:,1,:,:]), size = (1,5,1))
                         velRad = phase*lamb/(4*numpy.pi*tSamp)
                         #Spectral width
                         acf1 = ndimage.median_filter(numpy.abs(data_acf[:,1,:,:]), size = (1,5,1))
                         acf2 = ndimage.median_filter(numpy.abs(data_acf[:,2,:,:]), size = (1,5,1))
                         spcWidth = (lamb/(2*numpy.sqrt(6)*numpy.pi*tSamp))*numpy.sqrt(numpy.log(acf1/acf2))
             #             velRad = ndimage.median_filter(velRad, size = (1,5,1))
                         if allData:
                             boolMetFin = ~numpy.isnan(SNRdB)
                         else:
                             #SNR
                             boolMet1 = (SNRdB>SNRthresh) #SNR mask
                             boolMet1 = ndimage.median_filter(boolMet1, size=(1,5,5))
                             #Radial velocity
                             boolMet2 = numpy.abs(velRad) < 30
                             boolMet2 = ndimage.median_filter(boolMet2, (1,5,5))
                             #Spectral Width
                             boolMet3 = spcWidth < 30
                             boolMet3 = ndimage.median_filter(boolMet3, (1,5,5))
             #                 boolMetFin = self.__erase_small(boolMet1, 10,5)
                             boolMetFin = boolMet1&boolMet2&boolMet3
                         #Creating data_param
                         coordMet = numpy.where(boolMetFin)
                         cmet = coordMet[0]
                         tmet = coordMet[1]
                         hmet = coordMet[2]
                         data_param = numpy.zeros((tmet.size, 7))
                         data_param[:,0] = utctime
                         data_param[:,1] = cmet
                         data_param[:,2] = tmet
                         data_param[:,3] = hmet
                         data_param[:,4] = SNR[cmet,tmet,hmet].T
                         data_param[:,5] = velRad[cmet,tmet,hmet].T
                         data_param[:,6] = spcWidth[cmet,tmet,hmet].T
             #         self.dataOut.data_param = data_int
                     if len(data_param) == 0:
                         self.dataOut.flagNoData = True
                     else:
                         self.dataOut.data_param = data_param
                 def __erase_small(self, binArray, threshX, threshY):
                     labarray, numfeat = ndimage.measurements.label(binArray)
                     binArray1 = numpy.copy(binArray)
                     for i in range(1,numfeat + 1):
                         auxBin = (labarray==i)
                         auxSize = auxBin.sum()
                         x,y = numpy.where(auxBin)
                         widthX = x.max() - x.min()
                         widthY = y.max() - y.min()
                         #width X: 3 seg -> 12.5*3
                         #width Y:
                         if (auxSize < 50) or (widthX < threshX) or (widthY < threshY):
                             binArray1[auxBin] = False
                     return binArray1
             #---------------    Specular Meteor    ----------------
             class SMDetection(Operation):
                 '''
                     Function DetectMeteors()
                         Project developed with paper:
                         HOLDSWORTH ET AL. 2004
                     Input:
                         self.dataOut.data_pre
                         centerReceiverIndex:      From the channels, which is the center receiver
                         hei_ref:                  Height reference for the Beacon signal extraction
                         tauindex:
                         predefinedPhaseShifts:    Predefined phase offset for the voltge signals
                         cohDetection:             Whether to user Coherent detection or not
                         cohDet_timeStep:          Coherent Detection calculation time step
                         cohDet_thresh:            Coherent Detection phase threshold to correct phases
                         noise_timeStep:           Noise calculation time step
                         noise_multiple:           Noise multiple to define signal threshold
                         multDet_timeLimit:        Multiple Detection Removal time limit in seconds
                         multDet_rangeLimit:       Multiple Detection Removal range limit in km
                         phaseThresh:              Maximum phase difference between receiver to be consider a meteor
                         SNRThresh:                Minimum SNR threshold of the meteor signal to be consider a meteor
                         hmin:                     Minimum Height of the meteor to use it in the further wind estimations
                         hmax:                     Maximum Height of the meteor to use it in the further wind estimations
                         azimuth:                  Azimuth angle correction
                     Affected:
                         self.dataOut.data_param
                     Rejection Criteria (Errors):
 : No error; analysis OK
 : SNR < SNR threshold
 : angle of arrival (AOA) ambiguously determined
 : AOA estimate not feasible
 : Large difference in AOAs obtained from different antenna baselines
 : echo at start or end of time series
 : echo less than 5 examples long; too short for analysis
 : echo rise exceeds 0.3s
 : echo decay time less than twice rise time
 : large power level before echo
 : large power level after echo
 : poor fit to amplitude for estimation of decay time
 : poor fit to CCF phase variation for estimation of radial drift velocity
 : height unresolvable echo: not valid height within 70 to 110 km
 : height ambiguous echo: more then one possible height within 70 to 110 km
 : radial drift velocity or projected horizontal velocity exceeds 200 m/s
 : oscilatory echo, indicating event most likely not an underdense echo
 : phase difference in meteor Reestimation
                     Data Storage:
                         Meteors for Wind Estimation   (8):
                         Utc Time   |    Range    Height
                         Azimuth    Zenith    errorCosDir
                         VelRad    errorVelRad
                         Phase0 Phase1 Phase2 Phase3
                         TypeError
                      '''
                 def run(self, dataOut, hei_ref = None, tauindex = 0,
                                   phaseOffsets = None,
                                   cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25,
                                   noise_timeStep = 4, noise_multiple = 4,
                                   multDet_timeLimit = 1, multDet_rangeLimit = 3,
                                   phaseThresh = 20, SNRThresh = 5,
                                   hmin = 50, hmax=150, azimuth = 0,
                                   channelPositions = None) :
                     #Getting Pairslist
                     if channelPositions == None:
             #             channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)]   #T
                         channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)]   #Estrella
                     meteorOps = SMOperations()
                     pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
                     heiRang = dataOut.getHeiRange()
                     #Get Beacon signal - No Beacon signal anymore
             #         newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
             #
             #         if hei_ref != None:
             #             newheis = numpy.where(self.dataOut.heightList>hei_ref)
             #
                     #****************REMOVING HARDWARE PHASE DIFFERENCES***************
                     # see if the user put in pre defined phase shifts
                     voltsPShift = dataOut.data_pre.copy()
             #         if predefinedPhaseShifts != None:
             #             hardwarePhaseShifts = numpy.array(predefinedPhaseShifts)*numpy.pi/180
             #
             # #         elif beaconPhaseShifts:
             # #             #get hardware phase shifts using beacon signal
             # #             hardwarePhaseShifts = self.__getHardwarePhaseDiff(self.dataOut.data_pre, pairslist, newheis, 10)
             # #             hardwarePhaseShifts = numpy.insert(hardwarePhaseShifts,centerReceiverIndex,0)
             #
             #         else:
             #             hardwarePhaseShifts = numpy.zeros(5)
             #
             #         voltsPShift = numpy.zeros((self.dataOut.data_pre.shape[0],self.dataOut.data_pre.shape[1],self.dataOut.data_pre.shape[2]), dtype = 'complex')
             #         for i in range(self.dataOut.data_pre.shape[0]):
             #             voltsPShift[i,:,:] = self.__shiftPhase(self.dataOut.data_pre[i,:,:], hardwarePhaseShifts[i])
                     #******************END OF REMOVING HARDWARE PHASE DIFFERENCES*********
                     #Remove DC
                     voltsDC = numpy.mean(voltsPShift,1)
                     voltsDC = numpy.mean(voltsDC,1)
                     for i in range(voltsDC.shape[0]):
                         voltsPShift[i] = voltsPShift[i] - voltsDC[i]
                     #Don't considerate last heights, theyre used to calculate Hardware Phase Shift
             #         voltsPShift = voltsPShift[:,:,:newheis[0][0]]
                     #************ FIND POWER OF DATA W/COH OR NON COH DETECTION (3.4) **********
                     #Coherent Detection
                     if cohDetection:
                         #use coherent detection to get the net power
                         cohDet_thresh = cohDet_thresh*numpy.pi/180
                         voltsPShift = self.__coherentDetection(voltsPShift, cohDet_timeStep, dataOut.timeInterval, pairslist0, cohDet_thresh)
                     #Non-coherent detection!
                     powerNet = numpy.nansum(numpy.abs(voltsPShift[:,:,:])**2,0)
                     #********** END OF COH/NON-COH POWER CALCULATION**********************
                     #********** FIND THE NOISE LEVEL AND POSSIBLE METEORS ****************
                     #Get noise
                     noise, noise1 = self.__getNoise(powerNet, noise_timeStep, dataOut.timeInterval)
             #         noise = self.getNoise1(powerNet, noise_timeStep, self.dataOut.timeInterval)
                     #Get signal threshold
                     signalThresh = noise_multiple*noise
                     #Meteor echoes detection
                     listMeteors = self.__findMeteors(powerNet, signalThresh)
                     #******* END OF NOISE LEVEL AND POSSIBLE METEORS CACULATION **********
                     #************** REMOVE MULTIPLE DETECTIONS (3.5) ***************************
                     #Parameters
                     heiRange = dataOut.getHeiRange()
                     rangeInterval = heiRange[1] - heiRange[0]
                     rangeLimit = multDet_rangeLimit/rangeInterval
                     timeLimit = multDet_timeLimit/dataOut.timeInterval
                     #Multiple detection removals
                     listMeteors1 = self.__removeMultipleDetections(listMeteors, rangeLimit, timeLimit)
                     #************ END OF REMOVE MULTIPLE DETECTIONS **********************
                     #*********************     METEOR REESTIMATION  (3.7, 3.8, 3.9, 3.10)   ********************
                     #Parameters
                     phaseThresh = phaseThresh*numpy.pi/180
                     thresh = [phaseThresh, noise_multiple, SNRThresh]
                     #Meteor reestimation  (Errors N 1, 6, 12, 17)
                     listMeteors2, listMeteorsPower, listMeteorsVolts = self.__meteorReestimation(listMeteors1, voltsPShift, pairslist0, thresh, noise, dataOut.timeInterval, dataOut.frequency)
             #         listMeteors2, listMeteorsPower, listMeteorsVolts = self.meteorReestimation3(listMeteors2, listMeteorsPower, listMeteorsVolts, voltsPShift, pairslist, thresh, noise)
                     #Estimation of decay times (Errors N 7, 8, 11)
                     listMeteors3 = self.__estimateDecayTime(listMeteors2, listMeteorsPower, dataOut.timeInterval, dataOut.frequency)
                     #*******************     END OF METEOR REESTIMATION    *******************
                     #*********************    METEOR PARAMETERS CALCULATION (3.11, 3.12, 3.13)    **************************
                     #Calculating Radial Velocity (Error N 15)
                     radialStdThresh = 10
                     listMeteors4 = self.__getRadialVelocity(listMeteors3, listMeteorsVolts, radialStdThresh, pairslist0, dataOut.timeInterval)
                     if len(listMeteors4) > 0:
                         #Setting New Array
                         date = dataOut.utctime
                         arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang)
                         #Correcting phase offset
                         if phaseOffsets != None:
                             phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180
                             arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets)
                         #Second Pairslist
                         pairsList = []
                         pairx = (0,1)
                         pairy = (2,3)
                         pairsList.append(pairx)
                         pairsList.append(pairy)
                         jph = numpy.array([0,0,0,0])
                         h = (hmin,hmax)
                         arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, jph)
             #             #Calculate AOA (Error N 3, 4)
             #             #JONES ET AL. 1998
             #             error = arrayParameters[:,-1]
             #             AOAthresh = numpy.pi/8
             #             phases = -arrayParameters[:,9:13]
             #             arrayParameters[:,4:7], arrayParameters[:,-1] = meteorOps.getAOA(phases, pairsList, error, AOAthresh, azimuth)
             #
             #             #Calculate Heights (Error N 13 and 14)
             #             error = arrayParameters[:,-1]
             #             Ranges = arrayParameters[:,2]
             #             zenith = arrayParameters[:,5]
             #             arrayParameters[:,3], arrayParameters[:,-1] = meteorOps.getHeights(Ranges, zenith, error, hmin, hmax)
             #             error = arrayParameters[:,-1]
                     #*********************    END OF PARAMETERS CALCULATION    **************************
                     #***************************+     PASS DATA TO NEXT STEP    **********************
             #             arrayFinal = arrayParameters.reshape((1,arrayParameters.shape[0],arrayParameters.shape[1]))
                         dataOut.data_param = arrayParameters
                         if arrayParameters == None:
                             dataOut.flagNoData = True
                     else:
                         dataOut.flagNoData = True
                     return
                 def __getHardwarePhaseDiff(self, voltage0, pairslist, newheis, n):
                     minIndex = min(newheis[0])
                     maxIndex = max(newheis[0])
                     voltage = voltage0[:,:,minIndex:maxIndex+1]
                     nLength = voltage.shape[1]/n
                     nMin = 0
                     nMax = 0
                     phaseOffset = numpy.zeros((len(pairslist),n))
                     for i in range(n):
                         nMax += nLength
                         phaseCCF = -numpy.angle(self.__calculateCCF(voltage[:,nMin:nMax,:], pairslist, [0]))
                         phaseCCF = numpy.mean(phaseCCF, axis = 2)
                         phaseOffset[:,i] = phaseCCF.transpose()
                         nMin = nMax
             #         phaseDiff, phaseArrival = self.estimatePhaseDifference(voltage, pairslist)
                     #Remove Outliers
                     factor = 2
                     wt = phaseOffset - signal.medfilt(phaseOffset,(1,5))
                     dw = numpy.std(wt,axis = 1)
                     dw = dw.reshape((dw.size,1))
                     ind = numpy.where(numpy.logical_or(wt>dw*factor,wt<-dw*factor))
                     phaseOffset[ind] = numpy.nan
                     phaseOffset = stats.nanmean(phaseOffset, axis=1)
                     return phaseOffset
                 def __shiftPhase(self, data, phaseShift):
                     #this will shift the phase of a complex number
                     dataShifted = numpy.abs(data) * numpy.exp((numpy.angle(data)+phaseShift)*1j)
                     return dataShifted
                 def __estimatePhaseDifference(self, array, pairslist):
                     nChannel = array.shape[0]
                     nHeights = array.shape[2]
                     numPairs = len(pairslist)
             #         phaseCCF = numpy.zeros((nChannel, 5, nHeights))
                     phaseCCF = numpy.angle(self.__calculateCCF(array, pairslist, [-2,-1,0,1,2]))
                     #Correct phases
                     derPhaseCCF = phaseCCF[:,1:,:] - phaseCCF[:,0:-1,:]
                     indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
                     if indDer[0].shape[0] > 0:
                         for i in range(indDer[0].shape[0]):
                             signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i],indDer[2][i]])
                             phaseCCF[indDer[0][i],indDer[1][i]+1:,:] += signo*2*numpy.pi
             #         for j in range(numSides):
             #             phaseCCFAux = self.calculateCCF(arrayCenter, arraySides[j,:,:], [-2,1,0,1,2])
             #             phaseCCF[j,:,:] = numpy.angle(phaseCCFAux)
             #
                     #Linear
                     phaseInt = numpy.zeros((numPairs,1))
                     angAllCCF = phaseCCF[:,[0,1,3,4],0]
                     for j in range(numPairs):
                         fit = stats.linregress([-2,-1,1,2],angAllCCF[j,:])
                         phaseInt[j] = fit[1]
                     #Phase Differences
                     phaseDiff = phaseInt - phaseCCF[:,2,:]
                     phaseArrival = phaseInt.reshape(phaseInt.size)
                     #Dealias
                     phaseArrival = numpy.angle(numpy.exp(1j*phaseArrival))
             #         indAlias = numpy.where(phaseArrival > numpy.pi)
             #         phaseArrival[indAlias] -= 2*numpy.pi
             #         indAlias = numpy.where(phaseArrival < -numpy.pi)
             #         phaseArrival[indAlias] += 2*numpy.pi
                     return phaseDiff, phaseArrival
                 def __coherentDetection(self, volts, timeSegment, timeInterval, pairslist, thresh):
                     #this function will run the coherent detection used in Holdworth et al. 2004 and return the net power
                     #find the phase shifts of each channel over 1 second intervals
                     #only look at ranges below the beacon signal
                     numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
                     numBlocks = int(volts.shape[1]/numProfPerBlock)
                     numHeights = volts.shape[2]
                     nChannel = volts.shape[0]
                     voltsCohDet = volts.copy()
                     pairsarray = numpy.array(pairslist)
                     indSides = pairsarray[:,1]
             #         indSides = numpy.array(range(nChannel))
             #         indSides = numpy.delete(indSides, indCenter)
             #
             #         listCenter = numpy.array_split(volts[indCenter,:,:], numBlocks, 0)
                     listBlocks = numpy.array_split(volts, numBlocks, 1)
                     startInd = 0
                     endInd = 0
                     for i in range(numBlocks):
                         startInd = endInd
                         endInd = endInd + listBlocks[i].shape[1]
                         arrayBlock = listBlocks[i]
             #             arrayBlockCenter = listCenter[i]
                         #Estimate the Phase Difference
                         phaseDiff, aux = self.__estimatePhaseDifference(arrayBlock, pairslist)
                         #Phase Difference RMS
                         arrayPhaseRMS = numpy.abs(phaseDiff)
                         phaseRMSaux = numpy.sum(arrayPhaseRMS < thresh,0)
                         indPhase = numpy.where(phaseRMSaux==4)
                         #Shifting
                         if indPhase[0].shape[0] > 0:
                             for j in range(indSides.size):
                                 arrayBlock[indSides[j],:,indPhase] = self.__shiftPhase(arrayBlock[indSides[j],:,indPhase], phaseDiff[j,indPhase].transpose())
                             voltsCohDet[:,startInd:endInd,:] = arrayBlock
                     return voltsCohDet
                 def __calculateCCF(self, volts, pairslist ,laglist):
                     nHeights = volts.shape[2]
                     nPoints = volts.shape[1]
                     voltsCCF = numpy.zeros((len(pairslist), len(laglist), nHeights),dtype = 'complex')
                     for i in range(len(pairslist)):
                         volts1 = volts[pairslist[i][0]]
                         volts2 = volts[pairslist[i][1]]
                         for t in range(len(laglist)):
                             idxT = laglist[t]
                             if idxT >= 0:
                                 vStacked = numpy.vstack((volts2[idxT:,:],
                                                        numpy.zeros((idxT, nHeights),dtype='complex')))
                             else:
                                 vStacked = numpy.vstack((numpy.zeros((-idxT, nHeights),dtype='complex'),
                                                        volts2[:(nPoints + idxT),:]))
                             voltsCCF[i,t,:] = numpy.sum((numpy.conjugate(volts1)*vStacked),axis=0)
                             vStacked = None
                     return voltsCCF
                 def __getNoise(self, power, timeSegment, timeInterval):
                     numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
                     numBlocks = int(power.shape[0]/numProfPerBlock)
                     numHeights = power.shape[1]
                     listPower = numpy.array_split(power, numBlocks, 0)
                     noise = numpy.zeros((power.shape[0], power.shape[1]))
                     noise1 = numpy.zeros((power.shape[0], power.shape[1]))
                     startInd = 0
                     endInd = 0
                     for i in range(numBlocks):             #split por canal
                         startInd = endInd
                         endInd = endInd + listPower[i].shape[0]
                         arrayBlock = listPower[i]
                         noiseAux = numpy.mean(arrayBlock, 0)
             #             noiseAux = numpy.median(noiseAux)
             #             noiseAux = numpy.mean(arrayBlock)
                         noise[startInd:endInd,:] = noise[startInd:endInd,:] + noiseAux
                         noiseAux1 = numpy.mean(arrayBlock)
                         noise1[startInd:endInd,:] = noise1[startInd:endInd,:] + noiseAux1
                     return noise, noise1
                 def __findMeteors(self, power, thresh):
                     nProf = power.shape[0]
                     nHeights = power.shape[1]
                     listMeteors = []
                     for i in range(nHeights):
                         powerAux = power[:,i]
                         threshAux = thresh[:,i]
                         indUPthresh = numpy.where(powerAux > threshAux)[0]
                         indDNthresh = numpy.where(powerAux <= threshAux)[0]
                         j = 0
                         while (j < indUPthresh.size - 2):
                             if (indUPthresh[j + 2] == indUPthresh[j] + 2):
                                 indDNAux = numpy.where(indDNthresh > indUPthresh[j])
                                 indDNthresh = indDNthresh[indDNAux]
                                 if (indDNthresh.size > 0):
                                     indEnd = indDNthresh[0] - 1
                                     indInit = indUPthresh[j]
                                     meteor = powerAux[indInit:indEnd + 1]
                                     indPeak = meteor.argmax() + indInit
                                     FLA = sum(numpy.conj(meteor)*numpy.hstack((meteor[1:],0)))
                                     listMeteors.append(numpy.array([i,indInit,indPeak,indEnd,FLA])) #CHEQUEAR!!!!!
                                     j = numpy.where(indUPthresh == indEnd)[0] + 1
                                 else: j+=1
                             else: j+=1
                     return listMeteors
                 def __removeMultipleDetections(self,listMeteors, rangeLimit, timeLimit):
                     arrayMeteors = numpy.asarray(listMeteors)
                     listMeteors1 = []
                     while arrayMeteors.shape[0] > 0:
                         FLAs = arrayMeteors[:,4]
                         maxFLA = FLAs.argmax()
                         listMeteors1.append(arrayMeteors[maxFLA,:])
                         MeteorInitTime = arrayMeteors[maxFLA,1]
                         MeteorEndTime = arrayMeteors[maxFLA,3]
                         MeteorHeight = arrayMeteors[maxFLA,0]
                         #Check neighborhood
                         maxHeightIndex = MeteorHeight + rangeLimit
                         minHeightIndex = MeteorHeight - rangeLimit
                         minTimeIndex = MeteorInitTime - timeLimit
                         maxTimeIndex = MeteorEndTime + timeLimit
                         #Check Heights
                         indHeight = numpy.logical_and(arrayMeteors[:,0] >= minHeightIndex, arrayMeteors[:,0] <= maxHeightIndex)
                         indTime = numpy.logical_and(arrayMeteors[:,3] >= minTimeIndex, arrayMeteors[:,1] <= maxTimeIndex)
                         indBoth = numpy.where(numpy.logical_and(indTime,indHeight))
                         arrayMeteors = numpy.delete(arrayMeteors, indBoth, axis = 0)
                     return listMeteors1
                 def __meteorReestimation(self, listMeteors, volts, pairslist, thresh, noise, timeInterval,frequency):
                     numHeights = volts.shape[2]
                     nChannel = volts.shape[0]
                     thresholdPhase = thresh[0]
                     thresholdNoise = thresh[1]
                     thresholdDB = float(thresh[2])
                     thresholdDB1 = 10**(thresholdDB/10)
                     pairsarray = numpy.array(pairslist)
                     indSides = pairsarray[:,1]
                     pairslist1 = list(pairslist)
                     pairslist1.append((0,1))
                     pairslist1.append((3,4))
                     listMeteors1 = []
                     listPowerSeries = []
                     listVoltageSeries = []
                     #volts has the war data
                     if frequency == 30e6:
                         timeLag = 45*10**-3
                     else:
                         timeLag = 15*10**-3
                     lag = numpy.ceil(timeLag/timeInterval)
                     for i in range(len(listMeteors)):
                         ######################   3.6 - 3.7 PARAMETERS REESTIMATION    #########################
                         meteorAux = numpy.zeros(16)
                         #Loading meteor Data (mHeight, mStart, mPeak, mEnd)
                         mHeight = listMeteors[i][0]
                         mStart = listMeteors[i][1]
                         mPeak = listMeteors[i][2]
                         mEnd = listMeteors[i][3]
                         #get the volt data between the start and end times of the meteor
                         meteorVolts = volts[:,mStart:mEnd+1,mHeight]
                         meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
                         #3.6. Phase Difference estimation
                         phaseDiff, aux = self.__estimatePhaseDifference(meteorVolts, pairslist)
                         #3.7. Phase difference removal & meteor start, peak and end times reestimated
                         #meteorVolts0.- all Channels, all Profiles
                         meteorVolts0 = volts[:,:,mHeight]
                         meteorThresh = noise[:,mHeight]*thresholdNoise
                         meteorNoise = noise[:,mHeight]
                         meteorVolts0[indSides,:] = self.__shiftPhase(meteorVolts0[indSides,:], phaseDiff) #Phase Shifting
                         powerNet0 = numpy.nansum(numpy.abs(meteorVolts0)**2, axis = 0)  #Power
                         #Times reestimation
                         mStart1 = numpy.where(powerNet0[:mPeak] < meteorThresh[:mPeak])[0]
                         if mStart1.size > 0:
                             mStart1 = mStart1[-1] + 1
                         else:
                             mStart1 = mPeak
                         mEnd1 = numpy.where(powerNet0[mPeak:] < meteorThresh[mPeak:])[0][0] + mPeak - 1
                         mEndDecayTime1 = numpy.where(powerNet0[mPeak:] < meteorNoise[mPeak:])[0]
                         if mEndDecayTime1.size == 0:
                             mEndDecayTime1 = powerNet0.size
                         else:
                             mEndDecayTime1 = mEndDecayTime1[0] + mPeak - 1
             #             mPeak1 = meteorVolts0[mStart1:mEnd1 + 1].argmax()
                         #meteorVolts1.- all Channels, from start to end
                         meteorVolts1 = meteorVolts0[:,mStart1:mEnd1 + 1]
                         meteorVolts2 = meteorVolts0[:,mPeak + lag:mEnd1 + 1]
                         if meteorVolts2.shape[1] == 0:
                             meteorVolts2 = meteorVolts0[:,mPeak:mEnd1 + 1]
                         meteorVolts1 = meteorVolts1.reshape(meteorVolts1.shape[0], meteorVolts1.shape[1], 1)
                         meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1], 1)
                         #####################    END PARAMETERS REESTIMATION    #########################
                         #####################   3.8 PHASE DIFFERENCE REESTIMATION  ########################
             #             if mEnd1 - mStart1 > 4:       #Error Number 6: echo less than 5 samples long; too short for analysis
                         if meteorVolts2.shape[1] > 0:
                             #Phase Difference re-estimation
                             phaseDiff1, phaseDiffint = self.__estimatePhaseDifference(meteorVolts2, pairslist1)   #Phase Difference Estimation
             #                 phaseDiff1, phaseDiffint = self.estimatePhaseDifference(meteorVolts2, pairslist)
                             meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1])
                             phaseDiff11 = numpy.reshape(phaseDiff1, (phaseDiff1.shape[0],1))
                             meteorVolts2[indSides,:] = self.__shiftPhase(meteorVolts2[indSides,:], phaseDiff11[0:4])     #Phase Shifting
                             #Phase Difference RMS
                             phaseRMS1 = numpy.sqrt(numpy.mean(numpy.square(phaseDiff1)))
                             powerNet1 = numpy.nansum(numpy.abs(meteorVolts1[:,:])**2,0)
                             #Data from Meteor
                             mPeak1 = powerNet1.argmax() + mStart1
                             mPeakPower1 = powerNet1.max()
                             noiseAux = sum(noise[mStart1:mEnd1 + 1,mHeight])
                             mSNR1 = (sum(powerNet1)-noiseAux)/noiseAux
                             Meteor1 = numpy.array([mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1])
                             Meteor1 = numpy.hstack((Meteor1,phaseDiffint))
                             PowerSeries  = powerNet0[mStart1:mEndDecayTime1 + 1]
                             #Vectorize
                             meteorAux[0:7] = [mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1]
                             meteorAux[7:11] = phaseDiffint[0:4]
                             #Rejection Criterions
                             if phaseRMS1 > thresholdPhase:  #Error Number 17: Phase variation
                                 meteorAux[-1] = 17
                             elif mSNR1 < thresholdDB1:      #Error Number 1: SNR < threshold dB
                                 meteorAux[-1] = 1
                         else:
                             meteorAux[0:4] = [mHeight, mStart, mPeak, mEnd]
                             meteorAux[-1] = 6 #Error Number 6: echo less than 5 samples long; too short for analysis
                             PowerSeries = 0
                         listMeteors1.append(meteorAux)
                         listPowerSeries.append(PowerSeries)
                         listVoltageSeries.append(meteorVolts1)
                     return listMeteors1, listPowerSeries, listVoltageSeries
                 def __estimateDecayTime(self, listMeteors, listPower, timeInterval, frequency):
                     threshError = 10
                     #Depending if it is 30 or 50 MHz
                     if frequency == 30e6:
                         timeLag = 45*10**-3
                     else:
                         timeLag = 15*10**-3
                     lag = numpy.ceil(timeLag/timeInterval)
                     listMeteors1 = []
                     for i in range(len(listMeteors)):
                         meteorPower = listPower[i]
                         meteorAux = listMeteors[i]
                         if meteorAux[-1] == 0:
                             try:
                                 indmax = meteorPower.argmax()
                                 indlag = indmax + lag
                                 y = meteorPower[indlag:]
                                 x = numpy.arange(0, y.size)*timeLag
                                 #first guess
                                 a = y[0]
                                 tau = timeLag
                                 #exponential fit
                                 popt, pcov = optimize.curve_fit(self.__exponential_function, x, y, p0 = [a, tau])
                                 y1 = self.__exponential_function(x, *popt)
                                 #error estimation
                                 error = sum((y - y1)**2)/(numpy.var(y)*(y.size - popt.size))
                                 decayTime = popt[1]
                                 riseTime = indmax*timeInterval
                                 meteorAux[11:13] = [decayTime, error]
                                 #Table items 7, 8 and 11
                                 if (riseTime > 0.3):            #Number 7: Echo rise exceeds 0.3s
                                     meteorAux[-1] = 7
                                 elif (decayTime < 2*riseTime) : #Number 8: Echo decay time less than than twice rise time
                                     meteorAux[-1] = 8
                                 if (error > threshError):       #Number 11: Poor fit to amplitude for estimation of decay time
                                     meteorAux[-1] = 11
                             except:
                                 meteorAux[-1] = 11
                         listMeteors1.append(meteorAux)
                     return listMeteors1
                 #Exponential Function
                 def __exponential_function(self, x, a, tau):
                     y = a*numpy.exp(-x/tau)
                     return y
                 def __getRadialVelocity(self, listMeteors, listVolts, radialStdThresh, pairslist,  timeInterval):
                     pairslist1 = list(pairslist)
                     pairslist1.append((0,1))
                     pairslist1.append((3,4))
                     numPairs = len(pairslist1)
                     #Time Lag
                     timeLag = 45*10**-3
                     c = 3e8
                     lag = numpy.ceil(timeLag/timeInterval)
                     freq = 30e6
                     listMeteors1 = []
                     for i in range(len(listMeteors)):
                         meteorAux = listMeteors[i]
                         if meteorAux[-1] == 0:
                             mStart = listMeteors[i][1]
                             mPeak = listMeteors[i][2]
                             mLag = mPeak - mStart + lag
                             #get the volt data between the start and end times of the meteor
                             meteorVolts = listVolts[i]
                             meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
                             #Get CCF
                             allCCFs = self.__calculateCCF(meteorVolts, pairslist1, [-2,-1,0,1,2])
                             #Method 2
                             slopes = numpy.zeros(numPairs)
                             time = numpy.array([-2,-1,1,2])*timeInterval
                             angAllCCF = numpy.angle(allCCFs[:,[0,1,3,4],0])
                             #Correct phases
                             derPhaseCCF = angAllCCF[:,1:] - angAllCCF[:,0:-1]
                             indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
                             if indDer[0].shape[0] > 0:
                                 for i in range(indDer[0].shape[0]):
                                     signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i]])
                                     angAllCCF[indDer[0][i],indDer[1][i]+1:] += signo*2*numpy.pi
             #                     fit = scipy.stats.linregress(numpy.array([-2,-1,1,2])*timeInterval, numpy.array([phaseLagN2s[i],phaseLagN1s[i],phaseLag1s[i],phaseLag2s[i]]))
                             for j in range(numPairs):
                                 fit = stats.linregress(time, angAllCCF[j,:])
                                 slopes[j] = fit[0]
                             #Remove Outlier
             #                 indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
             #                 slopes = numpy.delete(slopes,indOut)
             #                 indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
             #                 slopes = numpy.delete(slopes,indOut)
                             radialVelocity = -numpy.mean(slopes)*(0.25/numpy.pi)*(c/freq)
                             radialError = numpy.std(slopes)*(0.25/numpy.pi)*(c/freq)
                             meteorAux[-2] = radialError
                             meteorAux[-3] = radialVelocity
                             #Setting Error
                             #Number 15: Radial Drift velocity or projected horizontal velocity exceeds 200 m/s
                             if numpy.abs(radialVelocity) > 200:
                                 meteorAux[-1] = 15
                             #Number 12: Poor fit to CCF variation for estimation of radial drift velocity
                             elif radialError > radialStdThresh:
                                 meteorAux[-1] = 12
                         listMeteors1.append(meteorAux)
                     return listMeteors1
                 def __setNewArrays(self, listMeteors, date, heiRang):
                     #New arrays
                     arrayMeteors = numpy.array(listMeteors)
                     arrayParameters = numpy.zeros((len(listMeteors), 13))
                     #Date inclusion
             #         date = re.findall(r'\((.*?)\)', date)
             #         date = date[0].split(',')
             #         date = map(int, date)
             #
             #         if len(date)<6:
             #             date.append(0)
             #
             #         date = [date[0]*10000 + date[1]*100 + date[2], date[3]*10000 + date[4]*100 + date[5]]
             #         arrayDate = numpy.tile(date, (len(listMeteors), 1))
                     arrayDate = numpy.tile(date, (len(listMeteors)))
                     #Meteor array
             #         arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)]
             #         arrayMeteors = numpy.hstack((arrayDate, arrayMeteors))
                     #Parameters Array
                     arrayParameters[:,0] = arrayDate #Date
                     arrayParameters[:,1] = heiRang[arrayMeteors[:,0].astype(int)] #Range
                     arrayParameters[:,6:8] = arrayMeteors[:,-3:-1] #Radial velocity and its error
                     arrayParameters[:,8:12] = arrayMeteors[:,7:11]  #Phases
                     arrayParameters[:,-1] = arrayMeteors[:,-1]  #Error
                     return arrayParameters
             class CorrectSMPhases(Operation):
                 def run(self, dataOut, phaseOffsets, hmin = 50, hmax = 150, azimuth = 45, channelPositions = None):
                     arrayParameters = dataOut.data_param
                     pairsList = []
                     pairx = (0,1)
                     pairy = (2,3)
                     pairsList.append(pairx)
                     pairsList.append(pairy)
                     jph = numpy.zeros(4)
                     phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180
                 #         arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets)
                     arrayParameters[:,8:12] = numpy.angle(numpy.exp(1j*(arrayParameters[:,8:12] + phaseOffsets)))
                     meteorOps = SMOperations()
                     if channelPositions == None:
                 #             channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)]   #T
                         channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)]   #Estrella
                     pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
                     h = (hmin,hmax)
                     arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, jph)
                     dataOut.data_param = arrayParameters
                     return
             class SMPhaseCalibration(Operation):
                 __buffer = None
                 __initime = None
                 __dataReady = False
                 __isConfig = False
                 def __checkTime(self, currentTime, initTime, paramInterval, outputInterval):
                     dataTime = currentTime + paramInterval
                     deltaTime = dataTime - initTime
                     if deltaTime >= outputInterval or deltaTime < 0:
                         return True
                     return False
                 def __getGammas(self, pairs, d, phases):
                     gammas = numpy.zeros(2)
                     for i in range(len(pairs)):
                         pairi = pairs[i]
                         phip3 = phases[:,pairi[1]]
                         d3 = d[pairi[1]]
                         phip2 = phases[:,pairi[0]]
                         d2 = d[pairi[0]]
                         #Calculating gamma
             #             jdcos = alp1/(k*d1)
             #             jgamma = numpy.angle(numpy.exp(1j*(d0*alp1/d1 - alp0)))
                         jgamma = -phip2*d3/d2 - phip3
                         jgamma = numpy.angle(numpy.exp(1j*jgamma))
             #             jgamma[jgamma>numpy.pi] -= 2*numpy.pi
             #             jgamma[jgamma<-numpy.pi] += 2*numpy.pi
                         #Revised distribution
                         jgammaArray = numpy.hstack((jgamma,jgamma+0.5*numpy.pi,jgamma-0.5*numpy.pi))
                         #Histogram
                         nBins = 64.0
                         rmin = -0.5*numpy.pi
                         rmax = 0.5*numpy.pi
                         phaseHisto = numpy.histogram(jgammaArray, bins=nBins, range=(rmin,rmax))
                         meteorsY = phaseHisto[0]
                         phasesX = phaseHisto[1][:-1]
                         width = phasesX[1] - phasesX[0]
                         phasesX += width/2
                         #Gaussian aproximation
                         bpeak = meteorsY.argmax()
                         peak = meteorsY.max()
                         jmin = bpeak - 5
                         jmax = bpeak + 5 + 1
                         if jmin<0:
                             jmin = 0
                             jmax = 6
                         elif jmax > meteorsY.size:
                             jmin = meteorsY.size - 6
                             jmax = meteorsY.size
                         x0 = numpy.array([peak,bpeak,50])
                         coeff = optimize.leastsq(self.__residualFunction, x0, args=(meteorsY[jmin:jmax], phasesX[jmin:jmax]))
                         #Gammas
                         gammas[i] = coeff[0][1]
                     return gammas
                 def __residualFunction(self, coeffs, y, t):
                     return y - self.__gauss_function(t, coeffs)
                 def __gauss_function(self, t, coeffs):
                     return coeffs[0]*numpy.exp(-0.5*((t - coeffs[1]) / coeffs[2])**2)
                 def __getPhases(self, azimuth, h, pairsList, d, gammas, meteorsArray):
                     meteorOps = SMOperations()
                     nchan = 4
                     pairx = pairsList[0]
                     pairy = pairsList[1]
                     center_xangle = 0
                     center_yangle = 0
                     range_angle = numpy.array([10*numpy.pi,numpy.pi,numpy.pi/2,numpy.pi/4])
                     ntimes = len(range_angle)
                     nstepsx = 20.0
                     nstepsy = 20.0
                     for iz in range(ntimes):
                         min_xangle = -range_angle[iz]/2 + center_xangle
                         max_xangle = range_angle[iz]/2 + center_xangle
                         min_yangle = -range_angle[iz]/2 + center_yangle
                         max_yangle = range_angle[iz]/2 + center_yangle
                         inc_x = (max_xangle-min_xangle)/nstepsx
                         inc_y = (max_yangle-min_yangle)/nstepsy
                         alpha_y = numpy.arange(nstepsy)*inc_y + min_yangle
                         alpha_x = numpy.arange(nstepsx)*inc_x + min_xangle
                         penalty = numpy.zeros((nstepsx,nstepsy))
                         jph_array = numpy.zeros((nchan,nstepsx,nstepsy))
                         jph = numpy.zeros(nchan)
                         # Iterations looking for the offset
                         for iy in range(int(nstepsy)):
                             for ix in range(int(nstepsx)):
                                 jph[pairy[1]] = alpha_y[iy]
                                 jph[pairy[0]] = -gammas[1] - alpha_y[iy]*d[pairy[1]]/d[pairy[0]]
                                 jph[pairx[1]] = alpha_x[ix]
                                 jph[pairx[0]] = -gammas[0] - alpha_x[ix]*d[pairx[1]]/d[pairx[0]]
                                 jph_array[:,ix,iy] = jph
                                 meteorsArray1 = meteorOps.getMeteorParams(meteorsArray, azimuth, h, pairsList, d, jph)
                                 error = meteorsArray1[:,-1]
                                 ind1 = numpy.where(error==0)[0]
                                 penalty[ix,iy] = ind1.size
                         i,j = numpy.unravel_index(penalty.argmax(), penalty.shape)
                         phOffset = jph_array[:,i,j]
                         center_xangle = phOffset[pairx[1]]
                         center_yangle = phOffset[pairy[1]]
                     phOffset = numpy.angle(numpy.exp(1j*jph_array[:,i,j]))
                     phOffset = phOffset*180/numpy.pi
                     return phOffset
                 def run(self, dataOut, hmin, hmax, channelPositions=None, nHours = 1):
                     dataOut.flagNoData = True
                     self.__dataReady = False
                     dataOut.outputInterval = nHours*3600
                     if self.__isConfig == False:
             #                 self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
                         #Get Initial LTC time
                         self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
                         self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
                         self.__isConfig = True
                     if self.__buffer == None:
                         self.__buffer = dataOut.data_param.copy()
                     else:
                         self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
                     self.__dataReady = self.__checkTime(dataOut.utctime, self.__initime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
                     if self.__dataReady:
                         dataOut.utctimeInit = self.__initime
                         self.__initime += dataOut.outputInterval #to erase time offset
                         freq = dataOut.frequency
                         c = dataOut.C #m/s
                         lamb = c/freq
                         k = 2*numpy.pi/lamb
                         azimuth = 0
                         h = (hmin, hmax)
                         pairs = ((0,1),(2,3))
                         if channelPositions == None:
             #             channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)]   #T
                             channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)]   #Estrella
                         meteorOps = SMOperations()
                         pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
             #             distances1 = [-distances[0]*lamb, distances[1]*lamb, -distances[2]*lamb, distances[3]*lamb]
                         meteorsArray = self.__buffer
                         error = meteorsArray[:,-1]
                         boolError = (error==0)|(error==3)|(error==4)|(error==13)|(error==14)
                         ind1 = numpy.where(boolError)[0]
                         meteorsArray = meteorsArray[ind1,:]
                         meteorsArray[:,-1] = 0
                         phases = meteorsArray[:,8:12]
                         #Calculate Gammas
                         gammas = self.__getGammas(pairs, distances, phases)
             #             gammas = numpy.array([-21.70409463,45.76935864])*numpy.pi/180
                         #Calculate Phases
                         phasesOff = self.__getPhases(azimuth, h, pairs, distances, gammas, meteorsArray)
                         phasesOff = phasesOff.reshape((1,phasesOff.size))
                         dataOut.data_output = -phasesOff
                         dataOut.flagNoData = False
                         self.__buffer = None
                     return
             class SMOperations():
                 def __init__(self):
                     return
                 def getMeteorParams(self, arrayParameters0, azimuth, h, pairsList, distances, jph):
                     arrayParameters = arrayParameters0.copy()
                     hmin = h[0]
                     hmax = h[1]
                     #Calculate AOA (Error N 3, 4)
                     #JONES ET AL. 1998
                     AOAthresh = numpy.pi/8
                     error = arrayParameters[:,-1]
                     phases = -arrayParameters[:,8:12] + jph
             #         phases = numpy.unwrap(phases)
                     arrayParameters[:,3:6], arrayParameters[:,-1] = self.__getAOA(phases, pairsList, distances, error, AOAthresh, azimuth)
                     #Calculate Heights (Error N 13 and 14)
                     error = arrayParameters[:,-1]
                     Ranges = arrayParameters[:,1]
                     zenith = arrayParameters[:,4]
                     arrayParameters[:,2], arrayParameters[:,-1] = self.__getHeights(Ranges, zenith, error, hmin, hmax)
                     #----------------------- Get Final data    ------------------------------------
             #         error = arrayParameters[:,-1]
             #         ind1 = numpy.where(error==0)[0]
             #         arrayParameters = arrayParameters[ind1,:]
                     return arrayParameters
                 def __getAOA(self, phases, pairsList, directions, error, AOAthresh, azimuth):
                     arrayAOA = numpy.zeros((phases.shape[0],3))
                     cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList,directions)
                     arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
                     cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
                     arrayAOA[:,2] = cosDirError
                     azimuthAngle = arrayAOA[:,0]
                     zenithAngle = arrayAOA[:,1]
                     #Setting Error
                     indError = numpy.where(numpy.logical_or(error == 3, error == 4))[0]
                     error[indError] = 0
                     #Number 3: AOA not fesible
                     indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
                     error[indInvalid] = 3
                     #Number 4: Large difference in AOAs obtained from different antenna baselines
                     indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
                     error[indInvalid] = 4
                     return arrayAOA, error
                 def __getDirectionCosines(self, arrayPhase, pairsList, distances):
                     #Initializing some variables
                     ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
                     ang_aux = ang_aux.reshape(1,ang_aux.size)
                     cosdir = numpy.zeros((arrayPhase.shape[0],2))
                     cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
                     for i in range(2):
                         ph0 = arrayPhase[:,pairsList[i][0]]
                         ph1 = arrayPhase[:,pairsList[i][1]]
                         d0 = distances[pairsList[i][0]]
                         d1 = distances[pairsList[i][1]]
                         ph0_aux = ph0 + ph1
                         ph0_aux = numpy.angle(numpy.exp(1j*ph0_aux))
             #             ph0_aux[ph0_aux > numpy.pi] -= 2*numpy.pi
             #             ph0_aux[ph0_aux < -numpy.pi] += 2*numpy.pi
                         #First Estimation
                         cosdir0[:,i] = (ph0_aux)/(2*numpy.pi*(d0 - d1))
                         #Most-Accurate Second Estimation
                         phi1_aux =  ph0 - ph1
                         phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
                         #Direction Cosine 1
                         cosdir1 = (phi1_aux + ang_aux)/(2*numpy.pi*(d0 + d1))
                         #Searching the correct Direction Cosine
                         cosdir0_aux = cosdir0[:,i]
                         cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
                         #Minimum Distance
                         cosDiff = (cosdir1 - cosdir0_aux)**2
                         indcos = cosDiff.argmin(axis = 1)
                         #Saving Value obtained
                         cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
                     return cosdir0, cosdir
                 def __calculateAOA(self, cosdir, azimuth):
                     cosdirX = cosdir[:,0]
                     cosdirY = cosdir[:,1]
                     zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
                     azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth#0 deg north, 90 deg east
                     angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
                     return angles
                 def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
                     Ramb = 375  #Ramb = c/(2*PRF)
                     Re = 6371   #Earth Radius
                     heights = numpy.zeros(Ranges.shape)
                     R_aux = numpy.array([0,1,2])*Ramb
                     R_aux = R_aux.reshape(1,R_aux.size)
                     Ranges = Ranges.reshape(Ranges.size,1)
                     Ri = Ranges + R_aux
                     hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re
                     #Check if there is a height between 70 and 110 km
                     h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
                     ind_h = numpy.where(h_bool == 1)[0]
                     hCorr = hi[ind_h, :]
                     ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
                     hCorr = hi[ind_hCorr]
                     heights[ind_h] = hCorr
                     #Setting Error
                     #Number 13: Height unresolvable echo: not valid height within 70 to 110 km
                     #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
                     indError = numpy.where(numpy.logical_or(error == 13, error == 14))[0]
                     error[indError] = 0
                     indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
                     error[indInvalid2] = 14
                     indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
                     error[indInvalid1] = 13
                     return heights, error
                 def getPhasePairs(self, channelPositions):
                     chanPos = numpy.array(channelPositions)
                     listOper = list(itertools.combinations(range(5),2))
                     distances = numpy.zeros(4)
                     axisX = []
                     axisY = []
                     distX = numpy.zeros(3)
                     distY = numpy.zeros(3)
                     ix = 0
                     iy = 0
                     pairX = numpy.zeros((2,2))
                     pairY = numpy.zeros((2,2))
                     for i in range(len(listOper)):
                         pairi = listOper[i]
                         posDif = numpy.abs(chanPos[pairi[0],:] - chanPos[pairi[1],:])
                         if posDif[0] == 0:
                             axisY.append(pairi)
                             distY[iy] = posDif[1]
                             iy += 1
                         elif posDif[1] == 0:
                             axisX.append(pairi)
                             distX[ix] = posDif[0]
                             ix += 1
                     for i in range(2):
                         if i==0:
                             dist0 = distX
                             axis0 = axisX
                         else:
                             dist0 = distY
                             axis0 = axisY
                         side = numpy.argsort(dist0)[:-1]
                         axis0 = numpy.array(axis0)[side,:]
                         chanC = int(numpy.intersect1d(axis0[0,:], axis0[1,:])[0])
                         axis1 = numpy.unique(numpy.reshape(axis0,4))
                         side = axis1[axis1 != chanC]
                         diff1 = chanPos[chanC,i] - chanPos[side[0],i]
                         diff2 = chanPos[chanC,i] - chanPos[side[1],i]
                         if diff1<0:
                             chan2 = side[0]
                             d2 = numpy.abs(diff1)
                             chan1 = side[1]
                             d1 = numpy.abs(diff2)
                         else:
                             chan2 = side[1]
                             d2 = numpy.abs(diff2)
                             chan1 = side[0]
                             d1 = numpy.abs(diff1)
                         if i==0:
                             chanCX = chanC
                             chan1X = chan1
                             chan2X = chan2
                             distances[0:2] = numpy.array([d1,d2])
                         else:
                             chanCY = chanC
                             chan1Y = chan1
                             chan2Y = chan2
                             distances[2:4] = numpy.array([d1,d2])
             #         axisXsides = numpy.reshape(axisX[ix,:],4)
             #
             #         channelCentX = int(numpy.intersect1d(pairX[0,:], pairX[1,:])[0])
             #         channelCentY = int(numpy.intersect1d(pairY[0,:], pairY[1,:])[0])
             #
             #         ind25X = numpy.where(pairX[0,:] != channelCentX)[0][0]
             #         ind20X = numpy.where(pairX[1,:] != channelCentX)[0][0]
             #         channel25X = int(pairX[0,ind25X])
             #         channel20X = int(pairX[1,ind20X])
             #         ind25Y = numpy.where(pairY[0,:] != channelCentY)[0][0]
             #         ind20Y = numpy.where(pairY[1,:] != channelCentY)[0][0]
             #         channel25Y = int(pairY[0,ind25Y])
             #         channel20Y = int(pairY[1,ind20Y])
             #         pairslist = [(channelCentX, channel25X),(channelCentX, channel20X),(channelCentY,channel25Y),(channelCentY, channel20Y)]
                     pairslist = [(chanCX, chan1X),(chanCX, chan2X),(chanCY,chan1Y),(chanCY, chan2Y)]
                     return pairslist, distances
             #     def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth):
             #
             #         arrayAOA = numpy.zeros((phases.shape[0],3))
             #         cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList)
             #
             #         arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
             #         cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
             #         arrayAOA[:,2] = cosDirError
             #
             #         azimuthAngle = arrayAOA[:,0]
             #         zenithAngle = arrayAOA[:,1]
             #
             #         #Setting Error
             #         #Number 3: AOA not fesible
             #         indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
             #         error[indInvalid] = 3
             #         #Number 4: Large difference in AOAs obtained from different antenna baselines
             #         indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
             #         error[indInvalid] = 4
             #         return arrayAOA, error
             #
             #     def __getDirectionCosines(self, arrayPhase, pairsList):
             #
             #         #Initializing some variables
             #         ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
             #         ang_aux = ang_aux.reshape(1,ang_aux.size)
             #
             #         cosdir = numpy.zeros((arrayPhase.shape[0],2))
             #         cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
             #
             #
             #         for i in range(2):
             #             #First Estimation
             #             phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]]
             #             #Dealias
             #             indcsi = numpy.where(phi0_aux > numpy.pi)
             #             phi0_aux[indcsi] -= 2*numpy.pi
             #             indcsi = numpy.where(phi0_aux < -numpy.pi)
             #             phi0_aux[indcsi] += 2*numpy.pi
             #             #Direction Cosine 0
             #             cosdir0[:,i] = -(phi0_aux)/(2*numpy.pi*0.5)
             #
             #             #Most-Accurate Second Estimation
             #             phi1_aux =  arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]]
             #             phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
             #             #Direction Cosine 1
             #             cosdir1 = -(phi1_aux + ang_aux)/(2*numpy.pi*4.5)
             #
             #             #Searching the correct Direction Cosine
             #             cosdir0_aux = cosdir0[:,i]
             #             cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
             #             #Minimum Distance
             #             cosDiff = (cosdir1 - cosdir0_aux)**2
             #             indcos = cosDiff.argmin(axis = 1)
             #             #Saving Value obtained
             #             cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
             #
             #         return cosdir0, cosdir
             #
             #     def __calculateAOA(self, cosdir, azimuth):
             #         cosdirX = cosdir[:,0]
             #         cosdirY = cosdir[:,1]
             #
             #         zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
             #         azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth #0 deg north, 90 deg east
             #         angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
             #
             #         return angles
             #
             #     def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
             #
             #         Ramb = 375  #Ramb = c/(2*PRF)
             #         Re = 6371   #Earth Radius
             #         heights = numpy.zeros(Ranges.shape)
             #
             #         R_aux = numpy.array([0,1,2])*Ramb
             #         R_aux = R_aux.reshape(1,R_aux.size)
             #
             #         Ranges = Ranges.reshape(Ranges.size,1)
             #
             #         Ri = Ranges + R_aux
             #         hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re
             #
             #         #Check if there is a height between 70 and 110 km
             #         h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
             #         ind_h = numpy.where(h_bool == 1)[0]
             #
             #         hCorr = hi[ind_h, :]
             #         ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
             #
             #         hCorr = hi[ind_hCorr]
             #         heights[ind_h] = hCorr
             #
             #         #Setting Error
             #         #Number 13: Height unresolvable echo: not valid height within 70 to 110 km
             #         #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
             #
             #         indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
             #         error[indInvalid2] = 14
             #         indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
             #         error[indInvalid1] = 13
             #
             #         return heights, error
   No newline at end of file

schainpy/model/proc/jroproc_parameters.pyc binary modified 0 0

NO CONTENT: modified file, binary diff hidden

schainpy/model/proc/jroproc_spectra.pyc binary modified 0 0

NO CONTENT: modified file, binary diff hidden

schainpy/model/proc/jroproc_voltage.pyc binary modified 0 0

NO CONTENT: modified file, binary diff hidden

schainpy/scripts/schain.xml +1 -1

	@@ -1,1 +1,1
	1	<Project description="Segundo Test" id="191" name="test01"><ReadUnit datatype="VoltageReader" id="1911" inputId="0" name="VoltageReader"><Operation id="19111" name="run" priority="1" type="self"><Parameter format="str" id="191111" name="datatype" value="VoltageReader" /><Parameter format="str" id="191112" name="path" value="/home/erick/Documents/Data/Claire_Data/raw" /><Parameter format="date" id="191113" name="startDate" value="2017/07/26" /><Parameter format="date" id="191114" name="endDate" value="2017/10/28" /><Parameter format="time" id="191115" name="startTime" value="09:40:00" /><Parameter format="time" id="191116" name="endTime" value="23:59:00" /><Parameter format="int" id="191118" name="online" value="0" /><Parameter format="int" id="191119" name="walk" value="1" /></Operation><Operation id="19112" name="printNumberOfBlock" priority="2" type="self" /></ReadUnit><ProcUnit datatype="SpectraProc" id="1913" inputId="1912" name="SpectraProc"><Operation id="19131" name="run" priority="1" type="self"><Parameter format="int" id="191311" name="nFFTPoints" value="256" /><Parameter format="pairslist" id="191312" name="pairsList" value="(0,1),(0,2),(1,2)" /></Operation><Operation id="19132" name="removeDC" priority="2" type="self" /><Operation id="19133" name="IncohInt" priority="3" type="external"><Parameter format="float" id="191331" name="n" value="30" /></Operation><Operation id="19134" name="SpectraPlot" priority="4" type="external"><Parameter format="int" id="191341" name="id" value="11" /><Parameter format="str" id="191342" name="wintitle" value="SpectraPlot" /><Parameter format="str" id="191343" name="xaxis" value="velocity" /><Parameter format="int" id="191344" name="zmin" value="0" /><Parameter format="float" id="191345" name="ymin" value="1.5" /><Parameter format="int" id="191346" name="ymax" value="7" /><Parameter format="int" id="191347" name="zmax" value="40" /><Parameter format="int" id="191348" name="save" value="1" /></Operation><Operation id="19135" name="RTIPlot" priority="5" type="other"><Parameter format="int" id="191351" name="id" value="30" /><Parameter format="str" id="191352" name="wintitle" value="RTI" /><Parameter format="int" id="191353" name="zmin" value="0" /><Parameter format="int" id="191354" name="zmax" value="40" /><Parameter format="float" id="191355" name="ymin" value="1.5" /><Parameter format="int" id="191356" name="ymax" value="8" /><Parameter format="int" id="191357" name="showprofile" value="1" /><Parameter format="float" id="191358" name="xmin" value="0" /><Parameter format="float" id="191359" name="xmax" value="23.9" /><Parameter format="int" id="191360" name="save" value="1" /></Operation><Operation id="19136" name="CrossSpectraPlot" priority="6" type="other"><Parameter format="str" id="191361" name="phase_cmap" value="bwr" /><Parameter format="int" id="191362" name="id" value="2005" /><Parameter format="str" id="191363" name="wintitle" value="CrossSpectraPlot_ShortPulse" /><Parameter format="str" id="191364" name="xaxis" value="Velocity" /><Parameter format="float" id="191365" name="ymin" value="1.5" /><Parameter format="int" id="191366" name="ymax" value="7" /></Operation></ProcUnit><ProcUnit datatype="VoltageProc" id="1912" inputId="1911" name="VoltageProc"><Operation id="19121" name="run" priority="1" type="self" /><Operation id="19122" name="setRadarFrequency" priority="2" type="self"><Parameter format="float" id="191221" name="frequency" value="445.09e6" /></Operation><Operation id="19123" name="selectHeights" priority="3" type="self"><Parameter format="float" id="191231" name="minHei" value="0" /><Parameter format="float" id="191232" name="maxHei" value="64" /></Operation></ProcUnit><ProcUnit datatype="Parameters" id="1914" inputId="1913" name="ParametersProc"><Operation id="19141" name="run" priority="1" type="self" /><Operation id="19142" name="SpectralMoments" priority="2" type="other" /><Operation id="19143" name="FullSpectralAnalysis" priority="3" type="other"><Parameter format="float" id="191431" name="E01" value="1.5" /><Parameter format="float" id="191432" name="E02" value="1.5" /><Parameter format="float" id="191433" name="E12" value="0" /><Parameter format="float" id="191434" name="N01" value="0.875" /><Parameter format="float" id="191435" name="N02" value="-0.875" /><Parameter format="float" id="191436" name="N12" value="-1.75" /></Operation><Operation id="19144" name="WindProfilerPlot" priority="4" type="other"><Parameter format="int" id="191441" name="id" value="4" /><Parameter format="str" id="191442" name="wintitle" value="Wind Profiler" /><Parameter format="bool" id="191443" name="save" value="1" /><Parameter format="float" id="191444" name="ymin" value="1.5" /><Parameter format="int" id="191445" name="ymax" value="7" /><Parameter format="float" id="191446" name="zmin" value="-0.00" /><Parameter format="float" id="191447" name="zmax" value="1.1" /></Operation></ProcUnit></Project> No newline at end of file		1	<Project description="read bltr data sswma file" id="191" name="test1"><ReadUnit datatype="testBLTRReader" id="1911" inputId="0" name="testBLTRReader"><Operation id="19111" name="run" priority="1" type="self"><Parameter format="str" id="191111" name="datatype" value="testBLTRReader" /><Parameter format="str" id="191112" name="path" value="/media/erick/6F60F7113095A154/BLTR" /><Parameter format="date" id="191113" name="startDate" value="2017/01/17" /><Parameter format="date" id="191114" name="endDate" value="2018/01/01" /><Parameter format="time" id="191115" name="startTime" value="00:00:00" /><Parameter format="time" id="191116" name="endTime" value="23:59:59" /><Parameter format="str" id="191118" name="ext" value="sswma" /></Operation></ReadUnit><ProcUnit datatype="BLTRProcess" id="1912" inputId="1911" name="BLTRProcess"><Operation id="19121" name="run" priority="1" type="self" /><Operation id="19122" name="SnrFilter" priority="2" type="self"><Parameter format="float" id="191221" name="snr_val" value="-20" /><Parameter format="int" id="191222" name="modetofilter" value="2" /></Operation><Operation id="19123" name="OutliersFilter" priority="3" type="self"><Parameter format="str" id="191231" name="svalue" value="meridional" /><Parameter format="str" id="191232" name="svalue2" value="inTime" /><Parameter format="float" id="191233" name="method" value="0" /><Parameter format="float" id="191234" name="factor" value="2" /><Parameter format="float" id="191235" name="filter" value="0" /><Parameter format="float" id="191236" name="npoints" value="9" /><Parameter format="int" id="191237" name="modetofilter" value="2" /></Operation><Operation id="19124" name="OutliersFilter" priority="4" type="self"><Parameter format="str" id="191241" name="svalue" value="zonal" /><Parameter format="str" id="191242" name="svalue2" value="inTime" /><Parameter format="float" id="191243" name="method" value="0" /><Parameter format="float" id="191244" name="factor" value="2" /><Parameter format="float" id="191245" name="filter" value="0" /><Parameter format="float" id="191246" name="npoints" value="9" /><Parameter format="int" id="191247" name="modetofilter" value="2" /></Operation><Operation id="19125" name="OutliersFilter" priority="5" type="self"><Parameter format="str" id="191251" name="svalue" value="vertical" /><Parameter format="str" id="191252" name="svalue2" value="inHeight" /><Parameter format="float" id="191253" name="method" value="0" /><Parameter format="float" id="191254" name="factor" value="2" /><Parameter format="float" id="191255" name="filter" value="0" /><Parameter format="float" id="191256" name="npoints" value="9" /><Parameter format="int" id="191257" name="modetofilter" value="2" /></Operation><Operation id="19126" name="prePlot" priority="6" type="self"><Parameter format="int" id="191261" name="modeselect" value="1" /></Operation><Operation id="19127" name="WindProfilerPlot" priority="7" type="other"><Parameter format="int" id="191271" name="id" value="1" /><Parameter format="str" id="191272" name="wintitle" value="" /><Parameter format="intlist" id="191273" name="channelList" value="0" /><Parameter format="int" id="191274" name="SNRmin" value="-10" /><Parameter format="int" id="191275" name="SNRmax" value="50" /><Parameter format="float" id="191276" name="SNRthresh" value="0" /><Parameter format="float" id="191277" name="xmin" value="0" /><Parameter format="float" id="191278" name="xmax" value="24" /><Parameter format="float" id="191279" name="ymax" value="3" /><Parameter format="float" id="191280" name="zmin" value="-20" /><Parameter format="float" id="191281" name="zmax" value="20" /><Parameter format="float" id="191282" name="zmin_ver" value="-200" /><Parameter format="float" id="191283" name="zmax_ver" value="200" /></Operation><Operation id="19128" name="prePlot" priority="8" type="self"><Parameter format="int" id="191281" name="modeselect" value="2" /></Operation><Operation id="19129" name="WindProfilerPlot" priority="9" type="other"><Parameter format="int" id="191291" name="id" value="2" /><Parameter format="str" id="191292" name="wintitle" value="" /><Parameter format="bool" id="191293" name="save" value="1" /><Parameter format="str" id="191294" name="figpath" value="/media/erick/6F60F7113095A154/BLTR/" /><Parameter format="int" id="191295" name="SNRmin" value="-20" /><Parameter format="int" id="191296" name="SNRmax" value="40" /><Parameter format="float" id="191297" name="SNRthresh" value="0" /><Parameter format="float" id="191298" name="xmin" value="0" /><Parameter format="float" id="191299" name="xmax" value="24" /><Parameter format="float" id="191300" name="ymin" value="0" /><Parameter format="float" id="191301" name="ymax" value="10" /><Parameter format="float" id="191302" name="zmin" value="-4" /><Parameter format="float" id="191303" name="zmax" value="4" /><Parameter format="float" id="191304" name="zmin_ver" value="-200" /><Parameter format="float" id="191305" name="zmax_ver" value="200" /></Operation></ProcUnit></Project> No newline at end of file

schainpy/scripts/testBLTR_block.py +8 -7

             '''
             Created on Nov 09, 2016
             @author: roj- LouVD
             '''
             import os, sys
             path = os.path.split(os.getcwd())[0]
             path = os.path.split(path)[0]
             sys.path.insert(0, path)
             from schainpy.controller import Project
             filename = 'test1.xml'
             # path = '/home/jespinoza/workspace/data/bltr/'
-            path = '/home/erick/Documents/Data/BLTR_Data/sswma/'
+            path = '/media/erick/6F60F7113095A154/BLTR/'
             desc = "read bltr data sswma file"
-            figpath = '/home/erick/workspace'
+            figpath = '/media/erick/6F60F7113095A154/BLTR/'
             pathhdf5 = '/tmp/'
             controllerObj = Project()
             controllerObj.setup(id = '191', name='test1', description=desc)
             readUnitConfObj = controllerObj.addReadUnit(datatype='testBLTRReader',
                                                         path=path,
-                                                        startDate='2015/01/17',
+                                                        startDate='2017/01/17',
-                                                        endDate='2017/01/01',
+                                                        endDate='2018/01/01',
                                                         startTime='00:00:00',
                                                         endTime='23:59:59',
                                                         ext='sswma')
             procUnitConfObj1 = controllerObj.addProcUnit(datatype='BLTRProcess',
                                                          inputId=readUnitConfObj.getId())
             '''-------------------------------------------Processing--------------------------------------------'''
             '''MODE 1: LOW ATMOSPHERE: 0- 3 km'''
             # opObj10 = procUnitConfObj1.addOperation(name='SnrFilter')
             # opObj10.addParameter(name='snr_val', value='-10', format='float')
             # opObj10.addParameter(name='modetofilter', value='1', format='int')
             #
             # opObj10 = procUnitConfObj1.addOperation(name='OutliersFilter')
             # opObj10.addParameter(name='svalue', value='meridional', format='str')
             # opObj10.addParameter(name='svalue2', value='inTime', format='str')
             # opObj10.addParameter(name='method', value='0', format='float')
             # opObj10.addParameter(name='factor', value='1', format='float')
             # opObj10.addParameter(name='filter', value='0', format='float')
             # opObj10.addParameter(name='npoints', value='5', format='float')
             # opObj10.addParameter(name='modetofilter', value='1', format='int')
             # #
             # opObj10 = procUnitConfObj1.addOperation(name='OutliersFilter')
             # opObj10.addParameter(name='svalue', value='zonal', format='str')
             # opObj10.addParameter(name='svalue2', value='inTime', format='str')
             # opObj10.addParameter(name='method', value='0', format='float')
             # opObj10.addParameter(name='factor', value='1', format='float')
             # opObj10.addParameter(name='filter', value='0', format='float')
             # opObj10.addParameter(name='npoints', value='5', format='float')
             # opObj10.addParameter(name='modetofilter', value='1', format='int')
             # #
             # opObj10 = procUnitConfObj1.addOperation(name='OutliersFilter')
             # opObj10.addParameter(name='svalue', value='vertical', format='str')
             # opObj10.addParameter(name='svalue2', value='inHeight', format='str')
             # opObj10.addParameter(name='method', value='0', format='float')
             # opObj10.addParameter(name='factor', value='2', format='float')
             # opObj10.addParameter(name='filter', value='0', format='float')
             # opObj10.addParameter(name='npoints', value='9', format='float')
             # opObj10.addParameter(name='modetofilter', value='1', format='int')
             #
             ''' MODE 2: 0 - 10 km '''
             opObj10 = procUnitConfObj1.addOperation(name='SnrFilter')
             opObj10.addParameter(name='snr_val', value='-20', format='float')
             opObj10.addParameter(name='modetofilter', value='2', format='int')
             opObj10 = procUnitConfObj1.addOperation(name='OutliersFilter')
             opObj10.addParameter(name='svalue', value='meridional', format='str')
             opObj10.addParameter(name='svalue2', value='inTime', format='str')
             opObj10.addParameter(name='method', value='0', format='float')
             opObj10.addParameter(name='factor', value='2', format='float')
             opObj10.addParameter(name='filter', value='0', format='float')
             opObj10.addParameter(name='npoints', value='9', format='float')
             opObj10.addParameter(name='modetofilter', value='2', format='int')
             # #
             opObj10 = procUnitConfObj1.addOperation(name='OutliersFilter')
             opObj10.addParameter(name='svalue', value='zonal', format='str')
             opObj10.addParameter(name='svalue2', value='inTime', format='str')
             opObj10.addParameter(name='method', value='0', format='float')
             opObj10.addParameter(name='factor', value='2', format='float')
             opObj10.addParameter(name='filter', value='0', format='float')
             opObj10.addParameter(name='npoints', value='9', format='float')
             opObj10.addParameter(name='modetofilter', value='2', format='int')
             # #
             opObj10 = procUnitConfObj1.addOperation(name='OutliersFilter')
             opObj10.addParameter(name='svalue', value='vertical', format='str')
             opObj10.addParameter(name='svalue2', value='inHeight', format='str')
             opObj10.addParameter(name='method', value='0', format='float')
             opObj10.addParameter(name='factor', value='2', format='float')
             opObj10.addParameter(name='filter', value='0', format='float')
             opObj10.addParameter(name='npoints', value='9', format='float')
             opObj10.addParameter(name='modetofilter', value='2', format='int')
             # '''-----------------------------------------Writing-------------------------------------------'''
             #
             # # opObj10 =  procUnitConfObj1.addOperation(name='testBLTRWriter',optype='other')
             # # opObj10.addParameter(name='path', value = pathhdf5)
             # # opObj10.addParameter(name='modetowrite', value = '2',format='int')
             # #
             # # opObj10 =  procUnitConfObj1.addOperation(name='testBLTRWriter',optype='other')
             # # opObj10.addParameter(name='path', value = pathhdf5)
             # # opObj10.addParameter(name='modetowrite', value = '1',format='int')
             #
             # '''----------------------------------------Plotting--------------------------------------------'''
             #
             opObj10 = procUnitConfObj1.addOperation(name='prePlot')
             opObj10.addParameter(name='modeselect',value='1',format='int')
             # #
             opObj10 = procUnitConfObj1.addOperation(name='WindProfilerPlot', optype='other')
             opObj10.addParameter(name='id', value='1', format='int')
             opObj10.addParameter(name='wintitle', value='', format='str')
             opObj10.addParameter(name='channelList', value='0', format='intlist')
             #opObj10.addParameter(name='save', value='1', format='bool')
             #opObj10.addParameter(name='figpath', value=figpath, format='str')
             opObj10.addParameter(name='SNRmin', value='-10', format='int')
             opObj10.addParameter(name='SNRmax', value='50', format='int')
             opObj10.addParameter(name='SNRthresh', value='0', format='float')
             opObj10.addParameter(name='xmin', value='0', format='float')
             opObj10.addParameter(name='xmax', value='24', format='float')
             opObj10.addParameter(name='ymax', value='3', format='float')
             opObj10.addParameter(name='zmin', value='-20', format='float')
             opObj10.addParameter(name='zmax', value='20', format='float')
             opObj10.addParameter(name='zmin_ver', value='-200', format='float')
             opObj10.addParameter(name='zmax_ver', value='200', format='float')
             #opObj10.addParameter(name='showprofile', value='1', format='bool')
             #opObj10.addParameter(name='show', value='1', format='bool')
             opObj10 = procUnitConfObj1.addOperation(name='prePlot')
             opObj10.addParameter(name='modeselect',value='2',format='int')
             #
             opObj10 = procUnitConfObj1.addOperation(name='WindProfilerPlot', optype='other')
             opObj10.addParameter(name='id', value='2', format='int')
             opObj10.addParameter(name='wintitle', value='', format='str')
             #opObj10.addParameter(name='channelList', value='0', format='intlist')
-            #opObj10.addParameter(name='save', value='1', format='bool')
+            opObj10.addParameter(name='save', value='1', format='bool')
-            #opObj10.addParameter(name='figpath', value=figpath, format='str')
+            opObj10.addParameter(name='figpath', value=figpath, format='str')
             opObj10.addParameter(name='SNRmin', value='-20', format='int')
             opObj10.addParameter(name='SNRmax', value='40', format='int')
             opObj10.addParameter(name='SNRthresh', value='0', format='float')
             opObj10.addParameter(name='xmin', value='0', format='float')
             opObj10.addParameter(name='xmax', value='24', format='float')
-            #opObj10.addParameter(name='ymax', value='8', format='float')
+            opObj10.addParameter(name='ymin', value='0', format='float')
+            opObj10.addParameter(name='ymax', value='10', format='float')
             opObj10.addParameter(name='zmin', value='-4', format='float')
             opObj10.addParameter(name='zmax', value='4', format='float')
             opObj10.addParameter(name='zmin_ver', value='-200', format='float')
             opObj10.addParameter(name='zmax_ver', value='200', format='float')
             #opObj10.addParameter(name='showprofile', value='1', format='bool')
             #opObj10.addParameter(name='show', value='1', format='bool')
             # # print "Escribiendo el archivo XML"
             # controllerObj.writeXml(filename)
             # # print "Leyendo el archivo XML"
             # controllerObj.readXml(filename)
             # controllerObj.createObjects()
             # controllerObj.connectObjects()
             # controllerObj.run()
             controllerObj.start()

schainpy/scripts/testProcBLTR.py +11 -7

             #!/usr/bin/env python
             import os, sys
             # path = os.path.dirname(os.getcwd())
             # path = os.path.join(path, 'source')
             # sys.path.insert(0, '../')
             from schainpy.controller import Project
             xmin = '15.5'
             xmax = '24'
             desc = "ProcBLTR Test"
             filename = "ProcBLTR.xml"
+            figpath = '/media/erick/6F60F7113095A154/BLTR'
             controllerObj = Project()
             controllerObj.setup(id='191', name='test01', description=desc)
             readUnitConfObj = controllerObj.addReadUnit(datatype='BLTRReader',
-                                                            path='/home/erick/Documents/Data/BLTR_Data/fdt/',
+                                                            path='/media/erick/6F60F7113095A154/BLTR/',
                                                             endDate='2017/10/19',
                                                             startTime='13:00:00',
                                                             startDate='2016/11/8',
                                                             endTime='23:59:59',
                                                             online=0,
                                                             walk=0,
                                                             ReadMode='1')
                                                             # expLabel='')
             # opObj11 = readUnitConfObj.addOperation(name='printNumberOfBlock')
             procUnitConfObj1 = controllerObj.addProcUnit(datatype='Spectra', inputId=readUnitConfObj.getId())
             opObj11 = procUnitConfObj1.addOperation(name='IncohInt', optype='other')
-            opObj11.addParameter(name='n', value='2', format='float')
+            opObj11.addParameter(name='n', value='3', format='float')
             opObj10 = procUnitConfObj1.addOperation(name='removeDC')
             # opObj10 = procUnitConfObj1.addOperation(name='calcMag')
             # opObj11 = procUnitConfObj1.addOperation(name='SpectraPlot', optype='other')
             # opObj11.addParameter(name='id', value='21', format='int')
             # opObj11.addParameter(name='wintitle', value='SpectraCutPlot', format='str')
             # opObj11.addParameter(name='xaxis', value='frequency', format='str')
             # opObj11.addParameter(name='colormap', value='winter', format='str')
             # opObj11.addParameter(name='xmin', value='-0.005', format='float')
             # opObj11.addParameter(name='xmax', value='0.005', format='float')
             # #opObj10 = procUnitConfObj1.addOperation(name='selectChannels')
             # #opObj10.addParameter(name='channelList', value='0,1', format='intlist')
             # opObj11 = procUnitConfObj1.addOperation(name='SpectraPlot', optype='other')
             # opObj11.addParameter(name='id', value='21', format='int')
             # opObj11.addParameter(name='wintitle', value='SpectraPlot', format='str')
             # #opObj11.addParameter(name='xaxis', value='Velocity', format='str')
             # opObj11.addParameter(name='xaxis', value='velocity', format='str')
             # opObj11.addParameter(name='xmin', value='-0.005', format='float')
             # opObj11.addParameter(name='xmax', value='0.005', format='float')
             # opObj11.addParameter(name='ymin', value='225', format='float')
             # opObj11.addParameter(name='ymax', value='3000', format='float')
             # opObj11.addParameter(name='zmin', value='-100', format='int')
             # opObj11.addParameter(name='zmax', value='-65', format='int')
             # opObj11 = procUnitConfObj1.addOperation(name='RTIPlot', optype='other')
             # opObj11.addParameter(name='id', value='10', format='int')
             # opObj11.addParameter(name='wintitle', value='RTI', format='str')
              # opObj11.addParameter(name='ymin', value='0', format='float')
              # opObj11.addParameter(name='ymax', value='4000', format='float')
             # #opObj11.addParameter(name='zmin', value='-100', format='int')
             # #opObj11.addParameter(name='zmax', value='-70', format='int')
              # opObj11.addParameter(name='zmin', value='-90', format='int')
              # opObj11.addParameter(name='zmax', value='-40', format='int')
             # opObj11.addParameter(name='showprofile', value='1', format='int')
             # opObj11.addParameter(name='timerange', value=str(2*60*60), format='int')
             opObj11 = procUnitConfObj1.addOperation(name='CrossSpectraPlot', optype='other')
             procUnitConfObj1.addParameter(name='pairsList', value='(0,1),(0,2),(1,2)', format='pairsList')
             opObj11.addParameter(name='id', value='2005', format='int')
             opObj11.addParameter(name='wintitle', value='CrossSpectraPlot_ShortPulse', format='str')
             # opObj11.addParameter(name='exp_code', value='13', format='int')
             opObj11.addParameter(name='xaxis', value='Velocity', format='str')
             #opObj11.addParameter(name='xmin', value='-10', format='float')
             #opObj11.addParameter(name='xmax', value='10', format='float')
             #opObj11.addParameter(name='ymin', value='225', format='float')
             #opObj11.addParameter(name='ymax', value='3000', format='float')
             #opObj11.addParameter(name='phase_min', value='-4', format='int')
             #opObj11.addParameter(name='phase_max', value='4', format='int')
             # procUnitConfObj2 = controllerObj.addProcUnit(datatype='CorrelationProc', inputId=procUnitConfObj1.getId())
             # procUnitConfObj2.addParameter(name='pairsList', value='(0,1),(0,2),(1,2)', format='pairsList')
             procUnitConfObj2 = controllerObj.addProcUnit(datatype='Parameters', inputId=procUnitConfObj1.getId())
             opObj11 = procUnitConfObj2.addOperation(name='SpectralMoments', optype='other')
             opObj22 = procUnitConfObj2.addOperation(name='FullSpectralAnalysis', optype='other')
             #
             opObj22 = procUnitConfObj2.addOperation(name='WindProfilerPlot', optype='other')
             opObj22.addParameter(name='id', value='4', format='int')
             opObj22.addParameter(name='wintitle', value='Wind Profiler', format='str')
             opObj22.addParameter(name='save', value='1', format='bool')
             # opObj22.addParameter(name='figpath', value = '/home/erick/Pictures', format='str')
             opObj22.addParameter(name='zmin', value='-20', format='int')
             opObj22.addParameter(name='zmax', value='20', format='int')
             opObj22.addParameter(name='zmin_ver', value='-250', format='float')
             opObj22.addParameter(name='zmax_ver', value='250', format='float')
             opObj22.addParameter(name='SNRmin', value='-5', format='int')
             opObj22.addParameter(name='SNRmax', value='30', format='int')
             # opObj22.addParameter(name='SNRthresh', value='-3.5', format='float')
             opObj22.addParameter(name='xmin', value=0, format='float')
             opObj22.addParameter(name='xmax', value=24, format='float')
             opObj22.addParameter(name='ymin', value='225', format='float')
             #opObj22.addParameter(name='ymax', value='2000', format='float')
+            opObj22.addParameter(name='save', value='1', format='int')
+            opObj22.addParameter(name='figpath', value=figpath, format='str')
             # opObj11.addParameter(name='pairlist', value='(1,0),(0,2),(1,2)', format='pairsList')
             #opObj10 = procUnitConfObj1.addOperation(name='selectHeights')
             #opObj10.addParameter(name='minHei', value='225', format='float')
             #opObj10.addParameter(name='maxHei', value='1000', format='float')
             # opObj11 = procUnitConfObj1.addOperation(name='CoherenceMap', optype='other')
             # opObj11.addParameter(name='id', value='102', format='int')
             # opObj11.addParameter(name='wintitle', value='Coherence', format='str')
             # opObj11.addParameter(name='ymin', value='225', format='float')
             # opObj11.addParameter(name='ymax', value='4000', format='float')
             # opObj11.addParameter(name='phase_cmap', value='jet', format='str')
             # opObj11.addParameter(name='xmin', value='8.5', format='float')
             # opObj11.addParameter(name='xmax', value='9.5', format='float')
             # opObj11.addParameter(name='figpath', value=figpath, format='str')
             # opObj11.addParameter(name='save', value=1, format='bool')
             # opObj11.addParameter(name='pairsList', value='(1,0),(3,2)', format='pairsList')
             # opObj12 = procUnitConfObj1.addOperation(name='PublishData', optype='other')
             # opObj12.addParameter(name='zeromq', value=1, format='int')
             # opObj12.addParameter(name='verbose', value=0, format='bool')
             # opObj12.addParameter(name='server', value='erick2', format='str')
             controllerObj.start()

schainpy/scripts/wind004/wind_20161108_000201.png binary modified 0 0

NO CONTENT: modified file, binary diff hidden

General Comments 0

Write
Preview

You need to be logged in to leave comments. Login now

No TODOs yet

	Site-wide shortcuts
/	Use quick search box
g h	Goto home page
g g	Goto my private gists page
g G	Goto my public gists page
g 0-9	Goto bookmarked items from 0-9
n r	New repository page
n g	New gist page

	Repositories
g s	Goto summary page
g c	Goto changelog page
g f	Goto files page
g F	Goto files page with file search activated
g p	Goto pull requests page
g o	Goto repository settings
g O	Goto repository access permissions settings
t s	Toggle sidebar on some pages