##// END OF EJS Templates
schain
RSS
Changes in Correlations Processing unit
Julio Valdez -
r854:23b4ca740232
parent child
Show More
Show file before | Show file after | Hide comments
@@ -1,246 +1,178
1 import numpy
1 import numpy
2
2
3 from jroproc_base import ProcessingUnit, Operation
3 from jroproc_base import ProcessingUnit, Operation
4 from schainpy.model.data.jrodata import Correlation
4 from schainpy.model.data.jrodata import Correlation, hildebrand_sekhon
5
5
6 class CorrelationProc(ProcessingUnit):
6 class CorrelationProc(ProcessingUnit):
7
7
8 pairsList = None
9
10 data_cf = None
11
8 def __init__(self):
12 def __init__(self):
9
13
10 ProcessingUnit.__init__(self)
14 ProcessingUnit.__init__(self)
11
15
12 self.objectDict = {}
16 self.objectDict = {}
13 self.buffer = None
17 self.buffer = None
14 self.firstdatatime = None
18 self.firstdatatime = None
15 self.profIndex = 0
19 self.profIndex = 0
16 self.dataOut = Correlation()
20 self.dataOut = Correlation()
17
21
18 def __updateObjFromVoltage(self):
22 def __updateObjFromVoltage(self):
19
23
20 self.dataOut.timeZone = self.dataIn.timeZone
24 self.dataOut.timeZone = self.dataIn.timeZone
21 self.dataOut.dstFlag = self.dataIn.dstFlag
25 self.dataOut.dstFlag = self.dataIn.dstFlag
22 self.dataOut.errorCount = self.dataIn.errorCount
26 self.dataOut.errorCount = self.dataIn.errorCount
23 self.dataOut.useLocalTime = self.dataIn.useLocalTime
27 self.dataOut.useLocalTime = self.dataIn.useLocalTime
24
28
25 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
29 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
26 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
30 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
27 self.dataOut.channelList = self.dataIn.channelList
31 self.dataOut.channelList = self.dataIn.channelList
28 self.dataOut.heightList = self.dataIn.heightList
32 self.dataOut.heightList = self.dataIn.heightList
29 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
33 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
30 # self.dataOut.nHeights = self.dataIn.nHeights
34 # self.dataOut.nHeights = self.dataIn.nHeights
31 # self.dataOut.nChannels = self.dataIn.nChannels
35 # self.dataOut.nChannels = self.dataIn.nChannels
32 self.dataOut.nBaud = self.dataIn.nBaud
36 self.dataOut.nBaud = self.dataIn.nBaud
33 self.dataOut.nCode = self.dataIn.nCode
37 self.dataOut.nCode = self.dataIn.nCode
34 self.dataOut.code = self.dataIn.code
38 self.dataOut.code = self.dataIn.code
35 # self.dataOut.nProfiles = self.dataOut.nFFTPoints
39 # self.dataOut.nProfiles = self.dataOut.nFFTPoints
36 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
40 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
37 self.dataOut.utctime = self.firstdatatime
41 self.dataOut.utctime = self.firstdatatime
38 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
42 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
39 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
43 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
40 self.dataOut.nCohInt = self.dataIn.nCohInt
44 self.dataOut.nCohInt = self.dataIn.nCohInt
41 # self.dataOut.nIncohInt = 1
45 # self.dataOut.nIncohInt = 1
42 self.dataOut.ippSeconds = self.dataIn.ippSeconds
46 self.dataOut.ippSeconds = self.dataIn.ippSeconds
47 self.dataOut.nProfiles = self.dataIn.nProfiles
48 self.dataOut.utctime = self.dataIn.utctime
43 # self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
49 # self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
44
50
45 # self.dataOut.timeInterval = self.dataIn.timeInterval*self.dataOut.nPoints
51 # self.dataOut.timeInterval = self.dataIn.timeInterval*self.dataOut.nPoints
46
52
47
53
48 def removeDC(self, jspectra):
54 def removeDC(self, jspectra):
49
55
50 nChannel = jspectra.shape[0]
56 nChannel = jspectra.shape[0]
51
57
52 for i in range(nChannel):
58 for i in range(nChannel):
53 jspectra_tmp = jspectra[i,:,:]
59 jspectra_tmp = jspectra[i,:,:]
54 jspectra_DC = numpy.mean(jspectra_tmp,axis = 0)
60 jspectra_DC = numpy.mean(jspectra_tmp,axis = 0)
55
61
56 jspectra_tmp = jspectra_tmp - jspectra_DC
62 jspectra_tmp = jspectra_tmp - jspectra_DC
57 jspectra[i,:,:] = jspectra_tmp
63 jspectra[i,:,:] = jspectra_tmp
58
64
59 return jspectra
65 return jspectra
60
66
61
67
62 def removeNoise(self, mode = 2):
68 def removeNoise(self, mode = 2):
63 indR = numpy.where(self.dataOut.lagR == 0)[0][0]
69 indR = numpy.where(self.dataOut.lagR == 0)[0][0]
64 indT = numpy.where(self.dataOut.lagT == 0)[0][0]
70 indT = numpy.where(self.dataOut.lagT == 0)[0][0]
65
71
66 jspectra = self.dataOut.data_corr[:,:,indR,:]
72 jspectra = self.dataOut.data_corr[:,:,indR,:]
67
73
68 num_chan = jspectra.shape[0]
74 num_chan = jspectra.shape[0]
69 num_hei = jspectra.shape[2]
75 num_hei = jspectra.shape[2]
70
76
71 freq_dc = indT
77 freq_dc = indT
72 ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
78 ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
73
79
74 NPot = self.dataOut.getNoise(mode)
80 NPot = self.dataOut.getNoise(mode)
75 jspectra[:,freq_dc,:] = jspectra[:,freq_dc,:] - NPot
81 jspectra[:,freq_dc,:] = jspectra[:,freq_dc,:] - NPot
76 SPot = jspectra[:,freq_dc,:]
82 SPot = jspectra[:,freq_dc,:]
77 pairsAutoCorr = self.dataOut.getPairsAutoCorr()
83 pairsAutoCorr = self.dataOut.getPairsAutoCorr()
78 # self.dataOut.signalPotency = SPot
84 # self.dataOut.signalPotency = SPot
79 self.dataOut.noise = NPot
85 self.dataOut.noise = NPot
80 self.dataOut.SNR = (SPot/NPot)[pairsAutoCorr]
86 self.dataOut.SNR = (SPot/NPot)[pairsAutoCorr]
81 self.dataOut.data_corr[:,:,indR,:] = jspectra
87 self.dataOut.data_corr[:,:,indR,:] = jspectra
82
88
83 return 1
89 return 1
84
90
85
91 def run(self, lags=None, mode = 'time', pairsList=None, fullBuffer=False, nAvg = 1, removeDC = False, splitCF=False):
86 def calculateNormFactor(self):
87
88 pairsList = self.dataOut.pairsList
89 pairsAutoCorr = self.dataOut.pairsAutoCorr
90 nHeights = self.dataOut.nHeights
91 nPairs = len(pairsList)
92 normFactor = numpy.zeros((nPairs,nHeights))
93
94 indR = numpy.where(self.dataOut.lagR == 0)[0][0]
95 indT = numpy.where(self.dataOut.lagT == 0)[0][0]
96
97 for l in range(len(pairsList)):
98 firstChannel = pairsList[l][0]
99 secondChannel = pairsList[l][1]
100
101 AC1 = pairsAutoCorr[firstChannel]
102 AC2 = pairsAutoCorr[secondChannel]
103
104 if (AC1 >= 0 and AC2 >= 0):
105
106 data1 = numpy.abs(self.dataOut.data_corr[AC1,:,indR,:])
107 data2 = numpy.abs(self.dataOut.data_corr[AC2,:,indR,:])
108 maxim1 = data1.max(axis = 0)
109 maxim2 = data1.max(axis = 0)
110 maxim = numpy.sqrt(maxim1*maxim2)
111 else:
112 #In case there is no autocorrelation for the pair
113 data = numpy.abs(self.dataOut.data_corr[l,:,indR,:])
114 maxim = numpy.max(data, axis = 0)
115
116 normFactor[l,:] = maxim
117
118 self.dataOut.normFactor = normFactor
119
120 return 1
121
122 def run(self, lagT=None, lagR=None, pairsList=None,
123 nPoints=None, nAvg=None, bufferSize=None,
124 fullT = False, fullR = False, removeDC = False):
125
92
126 self.dataOut.flagNoData = True
93 self.dataOut.flagNoData = True
127
94
128 if self.dataIn.type == "Correlation":
95 if self.dataIn.type == "Correlation":
129
96
130 self.dataOut.copy(self.dataIn)
97 self.dataOut.copy(self.dataIn)
131
98
132 return
99 return
133
100
134 if self.dataIn.type == "Voltage":
101 if self.dataIn.type == "Voltage":
135
102
136 if pairsList == None:
103 nChannels = self.dataIn.nChannels
137 pairsList = [numpy.array([0,0])]
104 nProfiles = self.dataIn.nProfiles
138
105 nHeights = self.dataIn.nHeights
139 if nPoints == None:
106 data_pre = self.dataIn.data
140 nPoints = 128
141 #------------------------------------------------------------
142 #Condicionales para calcular Correlaciones en Tiempo y Rango
143 if fullT:
144 lagT = numpy.arange(nPoints*2 - 1) - nPoints + 1
145 elif lagT == None:
146 lagT = numpy.array([0])
147 else:
148 lagT = numpy.array(lagT)
149
150 if fullR:
151 lagR = numpy.arange(self.dataOut.nHeights)
152 elif lagR == None:
153 lagR = numpy.array([0])
154 #-------------------------------------------------------------
155
156 if nAvg == None:
157 nAvg = 1
158
159 if bufferSize == None:
160 bufferSize = 0
161
162 deltaH = self.dataIn.heightList[1] - self.dataIn.heightList[0]
163 self.dataOut.lagR = numpy.round(numpy.array(lagR)/deltaH)
164 self.dataOut.pairsList = pairsList
165 self.dataOut.nPoints = nPoints
166 # channels = numpy.sort(list(set(list(itertools.chain.from_iterable(pairsList)))))
167
107
168 if self.buffer is None:
108 #--------------- Remover DC ------------
109 if removeDC:
110 data_pre = self.removeDC(data_pre)
169
111
170 self.buffer = numpy.zeros((self.dataIn.nChannels,self.dataIn.nProfiles,self.dataIn.nHeights),dtype='complex')
112 #---------------------------------------------
171
113 # pairsList = list(ccfList)
172
114 # for i in acfList:
173 self.buffer[:,self.profIndex,:] = self.dataIn.data.copy()[:,:]
115 # pairsList.append((i,i))
116 #
117 # ccf_pairs = numpy.arange(len(ccfList))
118 # acf_pairs = numpy.arange(len(ccfList),len(pairsList))
119 self.__updateObjFromVoltage()
120 #----------------------------------------------------------------------
121 #Creating temporal buffers
122 if fullBuffer:
123 tmp = numpy.zeros((len(pairsList), len(lags), nProfiles, nHeights), dtype = 'complex')*numpy.nan
124 elif mode == 'time':
125 if lags == None:
126 lags = numpy.arange(-nProfiles+1, nProfiles)
127 tmp = numpy.zeros((len(pairsList), len(lags), nHeights),dtype='complex')
128 elif mode == 'height':
129 if lags == None:
130 lags = numpy.arange(-nHeights+1, nHeights)
131 tmp = numpy.zeros(len(pairsList), (len(lags), nProfiles),dtype='complex')
174
132
175 self.profIndex += 1
133 #For loop
176
134 for l in range(len(pairsList)):
177 if self.firstdatatime == None:
178
135
179 self.firstdatatime = self.dataIn.utctime
136 ch0 = pairsList[l][0]
180
137 ch1 = pairsList[l][1]
181 if self.profIndex == nPoints:
138
182
139 for i in range(len(lags)):
183 tmp = self.buffer[:,0:nPoints,:]
140 idx = lags[i]
184 self.buffer = None
185 self.buffer = tmp
186
187 #--------------- Remover DC ------------
188 if removeDC:
189 self.buffer = self.removeDC(self.buffer)
190 #---------------------------------------------
191 self.dataOut.data_volts = self.buffer
192 self.__updateObjFromVoltage()
193 self.dataOut.data_corr = numpy.zeros((len(pairsList),
194 len(lagT),len(lagR),
195 self.dataIn.nHeights),
196 dtype='complex')
197
198 for l in range(len(pairsList)):
199
200 firstChannel = pairsList[l][0]
201 secondChannel = pairsList[l][1]
202
203 tmp = None
204 tmp = numpy.zeros((len(lagT),len(lagR),self.dataIn.nHeights),dtype='complex')
205
141
206 for t in range(len(lagT)):
142 if idx >= 0:
143 if mode == 'time':
144 ccf0 = data_pre[ch0,:nProfiles-idx,:]*numpy.conj(data_pre[ch1,idx:,:]) #time
145 else:
146 ccf0 = data_pre[ch0,:,nHeights-idx]*numpy.conj(data_pre[ch1,:,idx:]) #heights
147 else:
148 if mode == 'time':
149 ccf0 = data_pre[ch0,-idx:,:]*numpy.conj(data_pre[ch1,:nProfiles+idx,:]) #time
150 else:
151 ccf0 = data_pre[ch0,:,-idx:]*numpy.conj(data_pre[ch1,:,:nHeights+idx]) #heights
207
152
208 for r in range(len(lagR)):
153 if fullBuffer:
209
154 tmp[l,i,:ccf0.shape[0],:] = ccf0
210 idxT = lagT[t]
155 else:
211 idxR = lagR[r]
156 tmp[l,i,:] = numpy.sum(ccf0, axis=0)
212
157
213 if idxT >= 0:
158 #-----------------------------------------------------------------
214 vStacked = numpy.vstack((self.buffer[secondChannel,idxT:,:],
159 if fullBuffer:
215 numpy.zeros((idxT,self.dataIn.nHeights),dtype='complex')))
160 tmp = numpy.sum(numpy.reshape(tmp,(tmp.shape[0],tmp.shape[1],tmp.shape[2]/nAvg,nAvg,tmp.shape[3])),axis=3)
216 else:
161 self.dataOut.nAvg = nAvg
217 vStacked = numpy.vstack((numpy.zeros((-idxT,self.dataIn.nHeights),dtype='complex'),
162
218 self.buffer[secondChannel,:(nPoints + idxT),:]))
163 self.dataOut.data_cf = tmp
219
164 self.dataOut.mode = mode
220 if idxR >= 0:
165 self.dataOut.nLags = len(lags)
221 hStacked = numpy.hstack((vStacked[:,idxR:],numpy.zeros((nPoints,idxR),dtype='complex')))
166 self.dataOut.pairsList = pairsList
222 else:
167 self.dataOut.nPairs = len(pairsList)
223 hStacked = numpy.hstack((numpy.zeros((nPoints,-idxR),dtype='complex'),vStacked[:,(self.dataOut.nHeights + idxR)]))
168
224
169 #Se Calcula los factores de Normalizacion
225
170 if mode == 'time':
226 tmp[t,r,:] = numpy.sum((numpy.conjugate(self.buffer[firstChannel,:,:])*hStacked),axis=0)
171 delta = self.dataIn.ippSeconds*self.dataIn.nCohInt
227
172 else:
228
173 delta = self.dataIn.heightList[1] - self.dataIn.heightList[0]
229 hStacked = None
174 self.dataOut.lagRange = numpy.array(lags)*delta
230 vStacked = None
175 # self.dataOut.nCohInt = self.dataIn.nCohInt*nAvg
231
176 self.dataOut.flagNoData = False
232 self.dataOut.data_corr[l,:,:,:] = tmp[:,:,:]
177 a = self.dataOut.normFactor
233
178 return No newline at end of file
234 #Se Calcula los factores de Normalizacion
235 self.dataOut.pairsAutoCorr = self.dataOut.getPairsAutoCorr()
236 self.dataOut.lagT = lagT*self.dataIn.ippSeconds*self.dataIn.nCohInt
237 self.dataOut.lagR = lagR
238
239 self.calculateNormFactor()
240
241 self.dataOut.flagNoData = False
242 self.buffer = None
243 self.firstdatatime = None
244 self.profIndex = 0
245
246 return No newline at end of file
Show file before | Show file after | Hide comments
@@ -1,2722 +1,2739
1 import numpy
1 import numpy
2 import math
2 import math
3 from scipy import optimize, interpolate, signal, stats, ndimage
3 from scipy import optimize, interpolate, signal, stats, ndimage
4 import re
4 import re
5 import datetime
5 import datetime
6 import copy
6 import copy
7 import sys
7 import sys
8 import importlib
8 import importlib
9 import itertools
9 import itertools
10
10
11 from jroproc_base import ProcessingUnit, Operation
11 from jroproc_base import ProcessingUnit, Operation
12 from schainpy.model.data.jrodata import Parameters, hildebrand_sekhon
12 from schainpy.model.data.jrodata import Parameters, hildebrand_sekhon
13
13
14
14
15 class ParametersProc(ProcessingUnit):
15 class ParametersProc(ProcessingUnit):
16
16
17 nSeconds = None
17 nSeconds = None
18
18
19 def __init__(self):
19 def __init__(self):
20 ProcessingUnit.__init__(self)
20 ProcessingUnit.__init__(self)
21
21
22 # self.objectDict = {}
22 # self.objectDict = {}
23 self.buffer = None
23 self.buffer = None
24 self.firstdatatime = None
24 self.firstdatatime = None
25 self.profIndex = 0
25 self.profIndex = 0
26 self.dataOut = Parameters()
26 self.dataOut = Parameters()
27
27
28 def __updateObjFromInput(self):
28 def __updateObjFromInput(self):
29
29
30 self.dataOut.inputUnit = self.dataIn.type
30 self.dataOut.inputUnit = self.dataIn.type
31
31
32 self.dataOut.timeZone = self.dataIn.timeZone
32 self.dataOut.timeZone = self.dataIn.timeZone
33 self.dataOut.dstFlag = self.dataIn.dstFlag
33 self.dataOut.dstFlag = self.dataIn.dstFlag
34 self.dataOut.errorCount = self.dataIn.errorCount
34 self.dataOut.errorCount = self.dataIn.errorCount
35 self.dataOut.useLocalTime = self.dataIn.useLocalTime
35 self.dataOut.useLocalTime = self.dataIn.useLocalTime
36
36
37 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
37 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
38 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
38 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
39 self.dataOut.channelList = self.dataIn.channelList
39 self.dataOut.channelList = self.dataIn.channelList
40 self.dataOut.heightList = self.dataIn.heightList
40 self.dataOut.heightList = self.dataIn.heightList
41 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
41 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
42 # self.dataOut.nHeights = self.dataIn.nHeights
42 # self.dataOut.nHeights = self.dataIn.nHeights
43 # self.dataOut.nChannels = self.dataIn.nChannels
43 # self.dataOut.nChannels = self.dataIn.nChannels
44 self.dataOut.nBaud = self.dataIn.nBaud
44 self.dataOut.nBaud = self.dataIn.nBaud
45 self.dataOut.nCode = self.dataIn.nCode
45 self.dataOut.nCode = self.dataIn.nCode
46 self.dataOut.code = self.dataIn.code
46 self.dataOut.code = self.dataIn.code
47 # self.dataOut.nProfiles = self.dataOut.nFFTPoints
47 # self.dataOut.nProfiles = self.dataOut.nFFTPoints
48 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
48 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
49 # self.dataOut.utctime = self.firstdatatime
49 # self.dataOut.utctime = self.firstdatatime
50 self.dataOut.utctime = self.dataIn.utctime
50 self.dataOut.utctime = self.dataIn.utctime
51 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
51 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
52 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
52 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
53 self.dataOut.nCohInt = self.dataIn.nCohInt
53 self.dataOut.nCohInt = self.dataIn.nCohInt
54 # self.dataOut.nIncohInt = 1
54 # self.dataOut.nIncohInt = 1
55 self.dataOut.ippSeconds = self.dataIn.ippSeconds
55 self.dataOut.ippSeconds = self.dataIn.ippSeconds
56 # self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
56 # self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
57 self.dataOut.timeInterval = self.dataIn.timeInterval
57 self.dataOut.timeInterval = self.dataIn.timeInterval
58 self.dataOut.heightList = self.dataIn.getHeiRange()
58 self.dataOut.heightList = self.dataIn.getHeiRange()
59 self.dataOut.frequency = self.dataIn.frequency
59 self.dataOut.frequency = self.dataIn.frequency
60 self.dataOut.noise = self.dataIn.noise
60 self.dataOut.noise = self.dataIn.noise
61
61
62 def run(self):
62 def run(self):
63
63
64 #---------------------- Voltage Data ---------------------------
64 #---------------------- Voltage Data ---------------------------
65
65
66 if self.dataIn.type == "Voltage":
66 if self.dataIn.type == "Voltage":
67
67
68 self.__updateObjFromInput()
68 self.__updateObjFromInput()
69 self.dataOut.data_pre = self.dataIn.data.copy()
69 self.dataOut.data_pre = self.dataIn.data.copy()
70 self.dataOut.flagNoData = False
70 self.dataOut.flagNoData = False
71 self.dataOut.utctimeInit = self.dataIn.utctime
71 self.dataOut.utctimeInit = self.dataIn.utctime
72 self.dataOut.paramInterval = self.dataIn.nProfiles*self.dataIn.nCohInt*self.dataIn.ippSeconds
72 self.dataOut.paramInterval = self.dataIn.nProfiles*self.dataIn.nCohInt*self.dataIn.ippSeconds
73 return
73 return
74
74
75 #---------------------- Spectra Data ---------------------------
75 #---------------------- Spectra Data ---------------------------
76
76
77 if self.dataIn.type == "Spectra":
77 if self.dataIn.type == "Spectra":
78
78
79 self.dataOut.data_pre = (self.dataIn.data_spc,self.dataIn.data_cspc)
79 self.dataOut.data_pre = (self.dataIn.data_spc,self.dataIn.data_cspc)
80 self.dataOut.abscissaList = self.dataIn.getVelRange(1)
80 self.dataOut.abscissaList = self.dataIn.getVelRange(1)
81 self.dataOut.noise = self.dataIn.getNoise()
81 self.dataOut.noise = self.dataIn.getNoise()
82 self.dataOut.normFactor = self.dataIn.normFactor
82 self.dataOut.normFactor = self.dataIn.normFactor
83 self.dataOut.groupList = self.dataIn.pairsList
83 self.dataOut.groupList = self.dataIn.pairsList
84 self.dataOut.flagNoData = False
84 self.dataOut.flagNoData = False
85
85
86 #---------------------- Correlation Data ---------------------------
86 #---------------------- Correlation Data ---------------------------
87
87
88 if self.dataIn.type == "Correlation":
88 if self.dataIn.type == "Correlation":
89 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.dataIn.splitFunctions()
89
90
90 if self.dataIn.data_ccf is not None:
91 self.dataOut.data_pre = (self.dataIn.data_cf[acf_ind,:], self.dataIn.data_cf[ccf_ind,:,:])
91 self.dataOut.data_pre = (self.dataIn.data_acf,self.dataIn.data_ccf)
92 self.dataOut.normFactor = (self.dataIn.normFactor[acf_ind,:], self.dataIn.normFactor[ccf_ind,:])
92 else:
93 self.dataOut.groupList = (acf_pairs, ccf_pairs)
93 self.dataOut.data_pre = self.dataIn.data_acf.copy()
94
94
95 self.dataOut.abscissaList = self.dataIn.lagRange
95 self.dataOut.abscissaList = self.dataIn.lagRange
96 self.dataOut.noise = self.dataIn.noise
96 self.dataOut.noise = self.dataIn.noise
97 self.dataOut.normFactor = self.dataIn.normFactor
98 self.dataOut.data_SNR = self.dataIn.SNR
97 self.dataOut.data_SNR = self.dataIn.SNR
99 self.dataOut.groupList = self.dataIn.pairsList
100 self.dataOut.flagNoData = False
98 self.dataOut.flagNoData = False
101 self.dataOut.nAvg = self.dataIn.nAvg
99 self.dataOut.nAvg = self.dataIn.nAvg
102
100
103 #---------------------- Parameters Data ---------------------------
101 #---------------------- Parameters Data ---------------------------
104
102
105 if self.dataIn.type == "Parameters":
103 if self.dataIn.type == "Parameters":
106 self.dataOut.copy(self.dataIn)
104 self.dataOut.copy(self.dataIn)
107 self.dataOut.flagNoData = False
105 self.dataOut.flagNoData = False
108
106
109 return True
107 return True
110
108
111 self.__updateObjFromInput()
109 self.__updateObjFromInput()
112 self.dataOut.utctimeInit = self.dataIn.utctime
110 self.dataOut.utctimeInit = self.dataIn.utctime
113 self.dataOut.paramInterval = self.dataIn.timeInterval
111 self.dataOut.paramInterval = self.dataIn.timeInterval
114
112
115 return
113 return
116
114
117 class SpectralMoments(Operation):
115 class SpectralMoments(Operation):
118
116
119 '''
117 '''
120 Function SpectralMoments()
118 Function SpectralMoments()
121
119
122 Calculates moments (power, mean, standard deviation) and SNR of the signal
120 Calculates moments (power, mean, standard deviation) and SNR of the signal
123
121
124 Type of dataIn: Spectra
122 Type of dataIn: Spectra
125
123
124 Configuration Parameters:
125
126 dirCosx : Cosine director in X axis
127 dirCosy : Cosine director in Y axis
128
129 elevation :
130 azimuth :
131
126 Input:
132 Input:
127 channelList : simple channel list to select e.g. [2,3,7]
133 channelList : simple channel list to select e.g. [2,3,7]
128 self.dataOut.data_pre : Spectral data
134 self.dataOut.data_pre : Spectral data
129 self.dataOut.abscissaList : List of frequencies
135 self.dataOut.abscissaList : List of frequencies
130 self.dataOut.noise : Noise level per channel
136 self.dataOut.noise : Noise level per channel
131
137
132 Affected:
138 Affected:
133 self.dataOut.data_param : Parameters per channel
139 self.dataOut.data_param : Parameters per channel
134 self.dataOut.data_SNR : SNR per channel
140 self.dataOut.data_SNR : SNR per channel
135
141
136 '''
142 '''
137
143
138 def run(self, dataOut, channelList=None):
144 def run(self, dataOut):
139
145
140 dataOut.data_pre = dataOut.data_pre[0]
146 dataOut.data_pre = dataOut.data_pre[0]
141 data = dataOut.data_pre
147 data = dataOut.data_pre
142 absc = dataOut.abscissaList[:-1]
148 absc = dataOut.abscissaList[:-1]
143 noise = dataOut.noise
149 noise = dataOut.noise
144
150 nChannel = data.shape[0]
145 data_param = numpy.zeros((data.shape[0], 4, data.shape[2]))
151 data_param = numpy.zeros((nChannel, 4, data.shape[2]))
146
152
147 if channelList== None:
153 for ind in range(nChannel):
148 channelList = dataOut.channelList
149 dataOut.channelList = channelList
150
151 for ind in channelList:
152 data_param[ind,:,:] = self.__calculateMoments(data[ind,:,:], absc, noise[ind])
154 data_param[ind,:,:] = self.__calculateMoments(data[ind,:,:], absc, noise[ind])
153
155
154 dataOut.data_param = data_param[:,1:,:]
156 dataOut.data_param = data_param[:,1:,:]
155 dataOut.data_SNR = data_param[:,0]
157 dataOut.data_SNR = data_param[:,0]
156 return
158 return
157
159
158 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):
160 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):
159
161
160 if (nicoh == None): nicoh = 1
162 if (nicoh == None): nicoh = 1
161 if (graph == None): graph = 0
163 if (graph == None): graph = 0
162 if (smooth == None): smooth = 0
164 if (smooth == None): smooth = 0
163 elif (self.smooth < 3): smooth = 0
165 elif (self.smooth < 3): smooth = 0
164
166
165 if (type1 == None): type1 = 0
167 if (type1 == None): type1 = 0
166 if (fwindow == None): fwindow = numpy.zeros(oldfreq.size) + 1
168 if (fwindow == None): fwindow = numpy.zeros(oldfreq.size) + 1
167 if (snrth == None): snrth = -3
169 if (snrth == None): snrth = -3
168 if (dc == None): dc = 0
170 if (dc == None): dc = 0
169 if (aliasing == None): aliasing = 0
171 if (aliasing == None): aliasing = 0
170 if (oldfd == None): oldfd = 0
172 if (oldfd == None): oldfd = 0
171 if (wwauto == None): wwauto = 0
173 if (wwauto == None): wwauto = 0
172
174
173 if (n0 < 1.e-20): n0 = 1.e-20
175 if (n0 < 1.e-20): n0 = 1.e-20
174
176
175 freq = oldfreq
177 freq = oldfreq
176 vec_power = numpy.zeros(oldspec.shape[1])
178 vec_power = numpy.zeros(oldspec.shape[1])
177 vec_fd = numpy.zeros(oldspec.shape[1])
179 vec_fd = numpy.zeros(oldspec.shape[1])
178 vec_w = numpy.zeros(oldspec.shape[1])
180 vec_w = numpy.zeros(oldspec.shape[1])
179 vec_snr = numpy.zeros(oldspec.shape[1])
181 vec_snr = numpy.zeros(oldspec.shape[1])
180
182
181 for ind in range(oldspec.shape[1]):
183 for ind in range(oldspec.shape[1]):
182
184
183 spec = oldspec[:,ind]
185 spec = oldspec[:,ind]
184 aux = spec*fwindow
186 aux = spec*fwindow
185 max_spec = aux.max()
187 max_spec = aux.max()
186 m = list(aux).index(max_spec)
188 m = list(aux).index(max_spec)
187
189
188 #Smooth
190 #Smooth
189 if (smooth == 0): spec2 = spec
191 if (smooth == 0): spec2 = spec
190 else: spec2 = scipy.ndimage.filters.uniform_filter1d(spec,size=smooth)
192 else: spec2 = scipy.ndimage.filters.uniform_filter1d(spec,size=smooth)
191
193
192 # Calculo de Momentos
194 # Calculo de Momentos
193 bb = spec2[range(m,spec2.size)]
195 bb = spec2[range(m,spec2.size)]
194 bb = (bb<n0).nonzero()
196 bb = (bb<n0).nonzero()
195 bb = bb[0]
197 bb = bb[0]
196
198
197 ss = spec2[range(0,m + 1)]
199 ss = spec2[range(0,m + 1)]
198 ss = (ss<n0).nonzero()
200 ss = (ss<n0).nonzero()
199 ss = ss[0]
201 ss = ss[0]
200
202
201 if (bb.size == 0):
203 if (bb.size == 0):
202 bb0 = spec.size - 1 - m
204 bb0 = spec.size - 1 - m
203 else:
205 else:
204 bb0 = bb[0] - 1
206 bb0 = bb[0] - 1
205 if (bb0 < 0):
207 if (bb0 < 0):
206 bb0 = 0
208 bb0 = 0
207
209
208 if (ss.size == 0): ss1 = 1
210 if (ss.size == 0): ss1 = 1
209 else: ss1 = max(ss) + 1
211 else: ss1 = max(ss) + 1
210
212
211 if (ss1 > m): ss1 = m
213 if (ss1 > m): ss1 = m
212
214
213 valid = numpy.asarray(range(int(m + bb0 - ss1 + 1))) + ss1
215 valid = numpy.asarray(range(int(m + bb0 - ss1 + 1))) + ss1
214 power = ((spec2[valid] - n0)*fwindow[valid]).sum()
216 power = ((spec2[valid] - n0)*fwindow[valid]).sum()
215 fd = ((spec2[valid]- n0)*freq[valid]*fwindow[valid]).sum()/power
217 fd = ((spec2[valid]- n0)*freq[valid]*fwindow[valid]).sum()/power
216 w = math.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum()/power)
218 w = math.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum()/power)
217 snr = (spec2.mean()-n0)/n0
219 snr = (spec2.mean()-n0)/n0
218
220
219 if (snr < 1.e-20) :
221 if (snr < 1.e-20) :
220 snr = 1.e-20
222 snr = 1.e-20
221
223
222 vec_power[ind] = power
224 vec_power[ind] = power
223 vec_fd[ind] = fd
225 vec_fd[ind] = fd
224 vec_w[ind] = w
226 vec_w[ind] = w
225 vec_snr[ind] = snr
227 vec_snr[ind] = snr
226
228
227 moments = numpy.vstack((vec_snr, vec_power, vec_fd, vec_w))
229 moments = numpy.vstack((vec_snr, vec_power, vec_fd, vec_w))
228 return moments
230 return moments
229
231
230 #------------------ Get SA Parameters --------------------------
232 #------------------ Get SA Parameters --------------------------
231
233
232 def GetSAParameters(self):
234 def GetSAParameters(self):
233 #SA en frecuencia
235 #SA en frecuencia
234 pairslist = self.dataOut.groupList
236 pairslist = self.dataOut.groupList
235 num_pairs = len(pairslist)
237 num_pairs = len(pairslist)
236
238
237 vel = self.dataOut.abscissaList
239 vel = self.dataOut.abscissaList
238 spectra = self.dataOut.data_pre
240 spectra = self.dataOut.data_pre
239 cspectra = self.dataIn.data_cspc
241 cspectra = self.dataIn.data_cspc
240 delta_v = vel[1] - vel[0]
242 delta_v = vel[1] - vel[0]
241
243
242 #Calculating the power spectrum
244 #Calculating the power spectrum
243 spc_pow = numpy.sum(spectra, 3)*delta_v
245 spc_pow = numpy.sum(spectra, 3)*delta_v
244 #Normalizing Spectra
246 #Normalizing Spectra
245 norm_spectra = spectra/spc_pow
247 norm_spectra = spectra/spc_pow
246 #Calculating the norm_spectra at peak
248 #Calculating the norm_spectra at peak
247 max_spectra = numpy.max(norm_spectra, 3)
249 max_spectra = numpy.max(norm_spectra, 3)
248
250
249 #Normalizing Cross Spectra
251 #Normalizing Cross Spectra
250 norm_cspectra = numpy.zeros(cspectra.shape)
252 norm_cspectra = numpy.zeros(cspectra.shape)
251
253
252 for i in range(num_chan):
254 for i in range(num_chan):
253 norm_cspectra[i,:,:] = cspectra[i,:,:]/numpy.sqrt(spc_pow[pairslist[i][0],:]*spc_pow[pairslist[i][1],:])
255 norm_cspectra[i,:,:] = cspectra[i,:,:]/numpy.sqrt(spc_pow[pairslist[i][0],:]*spc_pow[pairslist[i][1],:])
254
256
255 max_cspectra = numpy.max(norm_cspectra,2)
257 max_cspectra = numpy.max(norm_cspectra,2)
256 max_cspectra_index = numpy.argmax(norm_cspectra, 2)
258 max_cspectra_index = numpy.argmax(norm_cspectra, 2)
257
259
258 for i in range(num_pairs):
260 for i in range(num_pairs):
259 cspc_par[i,:,:] = __calculateMoments(norm_cspectra)
261 cspc_par[i,:,:] = __calculateMoments(norm_cspectra)
260 #------------------- Get Lags ----------------------------------
262 #------------------- Get Lags ----------------------------------
261
263
262 class SALags(Operation):
264 class SALags(Operation):
263 '''
265 '''
264 Function GetMoments()
266 Function GetMoments()
265
267
266 Input:
268 Input:
267 self.dataOut.data_pre
269 self.dataOut.data_pre
268 self.dataOut.abscissaList
270 self.dataOut.abscissaList
269 self.dataOut.noise
271 self.dataOut.noise
270 self.dataOut.normFactor
272 self.dataOut.normFactor
271 self.dataOut.data_SNR
273 self.dataOut.data_SNR
272 self.dataOut.groupList
274 self.dataOut.groupList
273 self.dataOut.nChannels
275 self.dataOut.nChannels
274
276
275 Affected:
277 Affected:
276 self.dataOut.data_param
278 self.dataOut.data_param
277
279
278 '''
280 '''
279 def run(self, dataOut):
281 def run(self, dataOut):
280 data = dataOut.data_pre
281 data_acf = dataOut.data_pre[0]
282 data_acf = dataOut.data_pre[0]
282 data_ccf = dataOut.data_pre[1]
283 data_ccf = dataOut.data_pre[1]
284 normFactor_acf = dataOut.normFactor[0]
285 normFactor_ccf = dataOut.normFactor[1]
286 pairs_acf = dataOut.groupList[0]
287 pairs_ccf = dataOut.groupList[1]
283
288
284 normFactor = dataOut.normFactor
285 nHeights = dataOut.nHeights
289 nHeights = dataOut.nHeights
286 absc = dataOut.abscissaList[:-1]
290 absc = dataOut.abscissaList
287 noise = dataOut.noise
291 noise = dataOut.noise
288 SNR = dataOut.data_SNR
292 SNR = dataOut.data_SNR
289 # pairsList = dataOut.groupList
290 nChannels = dataOut.nChannels
293 nChannels = dataOut.nChannels
294 # pairsList = dataOut.groupList
291 # pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)
295 # pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)
292 dataOut.data_param = numpy.zeros((len(pairsCrossCorr)*2 + 1, nHeights))
296
293
297 for l in range(len(pairs_acf)):
294 dataNorm = numpy.abs(data)
298 data_acf[l,:,:] = data_acf[l,:,:]/normFactor_acf[l,:]
295 for l in range(len(pairsList)):
299
296 dataNorm[l,:,:] = dataNorm[l,:,:]/normFactor[l,:]
300 for l in range(len(pairs_ccf)):
301 data_ccf[l,:,:] = data_ccf[l,:,:]/normFactor_ccf[l,:]
297
302
298 self.dataOut.data_param[:-1,:] = self.__calculateTaus(dataNorm, pairsCrossCorr, pairsAutoCorr, absc)
303 dataOut.data_param = numpy.zeros((len(pairs_ccf)*2 + 1, nHeights))
299 self.dataOut.data_param[-1,:] = self.__calculateLag1Phase(data_acf, absc)
304 dataOut.data_param[:-1,:] = self.__calculateTaus(data_acf, data_ccf, absc)
305 dataOut.data_param[-1,:] = self.__calculateLag1Phase(data_acf, absc)
300 return
306 return
301
307
302 # def __getPairsAutoCorr(self, pairsList, nChannels):
308 # def __getPairsAutoCorr(self, pairsList, nChannels):
303 #
309 #
304 # pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
310 # pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
305 #
311 #
306 # for l in range(len(pairsList)):
312 # for l in range(len(pairsList)):
307 # firstChannel = pairsList[l][0]
313 # firstChannel = pairsList[l][0]
308 # secondChannel = pairsList[l][1]
314 # secondChannel = pairsList[l][1]
309 #
315 #
310 # #Obteniendo pares de Autocorrelacion
316 # #Obteniendo pares de Autocorrelacion
311 # if firstChannel == secondChannel:
317 # if firstChannel == secondChannel:
312 # pairsAutoCorr[firstChannel] = int(l)
318 # pairsAutoCorr[firstChannel] = int(l)
313 #
319 #
314 # pairsAutoCorr = pairsAutoCorr.astype(int)
320 # pairsAutoCorr = pairsAutoCorr.astype(int)
315 #
321 #
316 # pairsCrossCorr = range(len(pairsList))
322 # pairsCrossCorr = range(len(pairsList))
317 # pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
323 # pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
318 #
324 #
319 # return pairsAutoCorr, pairsCrossCorr
325 # return pairsAutoCorr, pairsCrossCorr
320
326
321 def __calculateTaus(self, data, pairsCrossCorr, pairsAutoCorr, lagTRange):
327 def __calculateTaus(self, data_acf, data_ccf, lagRange):
322
328
323 Pt0 = data.shape[1]/2
329 lag0 = data_acf.shape[1]/2
324 #Funcion de Autocorrelacion
330 #Funcion de Autocorrelacion
325 dataAutoCorr = stats.nanmean(data[pairsAutoCorr,:,:], axis = 0)
331 mean_acf = stats.nanmean(data_acf, axis = 0)
326
332
327 #Obtencion Indice de TauCross
333 #Obtencion Indice de TauCross
328 indCross = data[pairsCrossCorr,:,:].argmax(axis = 1)
334 ind_ccf = data_ccf.argmax(axis = 1)
329 #Obtencion Indice de TauAuto
335 #Obtencion Indice de TauAuto
330 indAuto = numpy.zeros(indCross.shape,dtype = 'int')
336 ind_acf = numpy.zeros(ind_ccf.shape,dtype = 'int')
331 CCValue = data[pairsCrossCorr,Pt0,:]
337 ccf_lag0 = data_ccf[:,lag0,:]
332 for i in range(pairsCrossCorr.size):
338
333 indAuto[i,:] = numpy.abs(dataAutoCorr - CCValue[i,:]).argmin(axis = 0)
339 for i in range(ccf_lag0.shape[0]):
340 ind_acf[i,:] = numpy.abs(mean_acf - ccf_lag0[i,:]).argmin(axis = 0)
334
341
335 #Obtencion de TauCross y TauAuto
342 #Obtencion de TauCross y TauAuto
336 tauCross = lagTRange[indCross]
343 tau_ccf = lagRange[ind_ccf]
337 tauAuto = lagTRange[indAuto]
344 tau_acf = lagRange[ind_acf]
338
345
339 Nan1, Nan2 = numpy.where(tauCross == lagTRange[0])
346 Nan1, Nan2 = numpy.where(tau_ccf == lagRange[0])
340
347
341 tauCross[Nan1,Nan2] = numpy.nan
348 tau_ccf[Nan1,Nan2] = numpy.nan
342 tauAuto[Nan1,Nan2] = numpy.nan
349 tau_acf[Nan1,Nan2] = numpy.nan
343 tau = numpy.vstack((tauCross,tauAuto))
350 tau = numpy.vstack((tau_ccf,tau_acf))
344
351
345 return tau
352 return tau
346
353
347 def __calculateLag1Phase(self, data, lagTRange):
354 def __calculateLag1Phase(self, data, lagTRange):
348 data1 = stats.nanmean(data, axis = 0)
355 data1 = stats.nanmean(data, axis = 0)
349 lag1 = numpy.where(lagTRange == 0)[0][0] + 1
356 lag1 = numpy.where(lagTRange == 0)[0][0] + 1
350
357
351 phase = numpy.angle(data1[lag1,:])
358 phase = numpy.angle(data1[lag1,:])
352
359
353 return phase
360 return phase
354
361
355 class SpectralFitting(Operation):
362 class SpectralFitting(Operation):
356 '''
363 '''
357 Function GetMoments()
364 Function GetMoments()
358
365
359 Input:
366 Input:
360 Output:
367 Output:
361 Variables modified:
368 Variables modified:
362 '''
369 '''
363
370
364 def run(self, dataOut, getSNR = True, path=None, file=None, groupList=None):
371 def run(self, dataOut, getSNR = True, path=None, file=None, groupList=None):
365
372
366
373
367 if path != None:
374 if path != None:
368 sys.path.append(path)
375 sys.path.append(path)
369 self.dataOut.library = importlib.import_module(file)
376 self.dataOut.library = importlib.import_module(file)
370
377
371 #To be inserted as a parameter
378 #To be inserted as a parameter
372 groupArray = numpy.array(groupList)
379 groupArray = numpy.array(groupList)
373 # groupArray = numpy.array([[0,1],[2,3]])
380 # groupArray = numpy.array([[0,1],[2,3]])
374 self.dataOut.groupList = groupArray
381 self.dataOut.groupList = groupArray
375
382
376 nGroups = groupArray.shape[0]
383 nGroups = groupArray.shape[0]
377 nChannels = self.dataIn.nChannels
384 nChannels = self.dataIn.nChannels
378 nHeights=self.dataIn.heightList.size
385 nHeights=self.dataIn.heightList.size
379
386
380 #Parameters Array
387 #Parameters Array
381 self.dataOut.data_param = None
388 self.dataOut.data_param = None
382
389
383 #Set constants
390 #Set constants
384 constants = self.dataOut.library.setConstants(self.dataIn)
391 constants = self.dataOut.library.setConstants(self.dataIn)
385 self.dataOut.constants = constants
392 self.dataOut.constants = constants
386 M = self.dataIn.normFactor
393 M = self.dataIn.normFactor
387 N = self.dataIn.nFFTPoints
394 N = self.dataIn.nFFTPoints
388 ippSeconds = self.dataIn.ippSeconds
395 ippSeconds = self.dataIn.ippSeconds
389 K = self.dataIn.nIncohInt
396 K = self.dataIn.nIncohInt
390 pairsArray = numpy.array(self.dataIn.pairsList)
397 pairsArray = numpy.array(self.dataIn.pairsList)
391
398
392 #List of possible combinations
399 #List of possible combinations
393 listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)
400 listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)
394 indCross = numpy.zeros(len(list(listComb)), dtype = 'int')
401 indCross = numpy.zeros(len(list(listComb)), dtype = 'int')
395
402
396 if getSNR:
403 if getSNR:
397 listChannels = groupArray.reshape((groupArray.size))
404 listChannels = groupArray.reshape((groupArray.size))
398 listChannels.sort()
405 listChannels.sort()
399 noise = self.dataIn.getNoise()
406 noise = self.dataIn.getNoise()
400 self.dataOut.data_SNR = self.__getSNR(self.dataIn.data_spc[listChannels,:,:], noise[listChannels])
407 self.dataOut.data_SNR = self.__getSNR(self.dataIn.data_spc[listChannels,:,:], noise[listChannels])
401
408
402 for i in range(nGroups):
409 for i in range(nGroups):
403 coord = groupArray[i,:]
410 coord = groupArray[i,:]
404
411
405 #Input data array
412 #Input data array
406 data = self.dataIn.data_spc[coord,:,:]/(M*N)
413 data = self.dataIn.data_spc[coord,:,:]/(M*N)
407 data = data.reshape((data.shape[0]*data.shape[1],data.shape[2]))
414 data = data.reshape((data.shape[0]*data.shape[1],data.shape[2]))
408
415
409 #Cross Spectra data array for Covariance Matrixes
416 #Cross Spectra data array for Covariance Matrixes
410 ind = 0
417 ind = 0
411 for pairs in listComb:
418 for pairs in listComb:
412 pairsSel = numpy.array([coord[x],coord[y]])
419 pairsSel = numpy.array([coord[x],coord[y]])
413 indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0])
420 indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0])
414 ind += 1
421 ind += 1
415 dataCross = self.dataIn.data_cspc[indCross,:,:]/(M*N)
422 dataCross = self.dataIn.data_cspc[indCross,:,:]/(M*N)
416 dataCross = dataCross**2/K
423 dataCross = dataCross**2/K
417
424
418 for h in range(nHeights):
425 for h in range(nHeights):
419 # print self.dataOut.heightList[h]
426 # print self.dataOut.heightList[h]
420
427
421 #Input
428 #Input
422 d = data[:,h]
429 d = data[:,h]
423
430
424 #Covariance Matrix
431 #Covariance Matrix
425 D = numpy.diag(d**2/K)
432 D = numpy.diag(d**2/K)
426 ind = 0
433 ind = 0
427 for pairs in listComb:
434 for pairs in listComb:
428 #Coordinates in Covariance Matrix
435 #Coordinates in Covariance Matrix
429 x = pairs[0]
436 x = pairs[0]
430 y = pairs[1]
437 y = pairs[1]
431 #Channel Index
438 #Channel Index
432 S12 = dataCross[ind,:,h]
439 S12 = dataCross[ind,:,h]
433 D12 = numpy.diag(S12)
440 D12 = numpy.diag(S12)
434 #Completing Covariance Matrix with Cross Spectras
441 #Completing Covariance Matrix with Cross Spectras
435 D[x*N:(x+1)*N,y*N:(y+1)*N] = D12
442 D[x*N:(x+1)*N,y*N:(y+1)*N] = D12
436 D[y*N:(y+1)*N,x*N:(x+1)*N] = D12
443 D[y*N:(y+1)*N,x*N:(x+1)*N] = D12
437 ind += 1
444 ind += 1
438 Dinv=numpy.linalg.inv(D)
445 Dinv=numpy.linalg.inv(D)
439 L=numpy.linalg.cholesky(Dinv)
446 L=numpy.linalg.cholesky(Dinv)
440 LT=L.T
447 LT=L.T
441
448
442 dp = numpy.dot(LT,d)
449 dp = numpy.dot(LT,d)
443
450
444 #Initial values
451 #Initial values
445 data_spc = self.dataIn.data_spc[coord,:,h]
452 data_spc = self.dataIn.data_spc[coord,:,h]
446
453
447 if (h>0)and(error1[3]<5):
454 if (h>0)and(error1[3]<5):
448 p0 = self.dataOut.data_param[i,:,h-1]
455 p0 = self.dataOut.data_param[i,:,h-1]
449 else:
456 else:
450 p0 = numpy.array(self.dataOut.library.initialValuesFunction(data_spc, constants, i))
457 p0 = numpy.array(self.dataOut.library.initialValuesFunction(data_spc, constants, i))
451
458
452 try:
459 try:
453 #Least Squares
460 #Least Squares
454 minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True)
461 minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True)
455 # minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))
462 # minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))
456 #Chi square error
463 #Chi square error
457 error0 = numpy.sum(infodict['fvec']**2)/(2*N)
464 error0 = numpy.sum(infodict['fvec']**2)/(2*N)
458 #Error with Jacobian
465 #Error with Jacobian
459 error1 = self.dataOut.library.errorFunction(minp,constants,LT)
466 error1 = self.dataOut.library.errorFunction(minp,constants,LT)
460 except:
467 except:
461 minp = p0*numpy.nan
468 minp = p0*numpy.nan
462 error0 = numpy.nan
469 error0 = numpy.nan
463 error1 = p0*numpy.nan
470 error1 = p0*numpy.nan
464
471
465 #Save
472 #Save
466 if self.dataOut.data_param == None:
473 if self.dataOut.data_param == None:
467 self.dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan
474 self.dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan
468 self.dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan
475 self.dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan
469
476
470 self.dataOut.data_error[i,:,h] = numpy.hstack((error0,error1))
477 self.dataOut.data_error[i,:,h] = numpy.hstack((error0,error1))
471 self.dataOut.data_param[i,:,h] = minp
478 self.dataOut.data_param[i,:,h] = minp
472 return
479 return
473
480
474 def __residFunction(self, p, dp, LT, constants):
481 def __residFunction(self, p, dp, LT, constants):
475
482
476 fm = self.dataOut.library.modelFunction(p, constants)
483 fm = self.dataOut.library.modelFunction(p, constants)
477 fmp=numpy.dot(LT,fm)
484 fmp=numpy.dot(LT,fm)
478
485
479 return dp-fmp
486 return dp-fmp
480
487
481 def __getSNR(self, z, noise):
488 def __getSNR(self, z, noise):
482
489
483 avg = numpy.average(z, axis=1)
490 avg = numpy.average(z, axis=1)
484 SNR = (avg.T-noise)/noise
491 SNR = (avg.T-noise)/noise
485 SNR = SNR.T
492 SNR = SNR.T
486 return SNR
493 return SNR
487
494
488 def __chisq(p,chindex,hindex):
495 def __chisq(p,chindex,hindex):
489 #similar to Resid but calculates CHI**2
496 #similar to Resid but calculates CHI**2
490 [LT,d,fm]=setupLTdfm(p,chindex,hindex)
497 [LT,d,fm]=setupLTdfm(p,chindex,hindex)
491 dp=numpy.dot(LT,d)
498 dp=numpy.dot(LT,d)
492 fmp=numpy.dot(LT,fm)
499 fmp=numpy.dot(LT,fm)
493 chisq=numpy.dot((dp-fmp).T,(dp-fmp))
500 chisq=numpy.dot((dp-fmp).T,(dp-fmp))
494 return chisq
501 return chisq
495
502
496 class WindProfiler(Operation):
503 class WindProfiler(Operation):
497
504
498 __isConfig = False
505 __isConfig = False
499
506
500 __initime = None
507 __initime = None
501 __lastdatatime = None
508 __lastdatatime = None
502 __integrationtime = None
509 __integrationtime = None
503
510
504 __buffer = None
511 __buffer = None
505
512
506 __dataReady = False
513 __dataReady = False
507
514
508 __firstdata = None
515 __firstdata = None
509
516
510 n = None
517 n = None
511
518
512 def __init__(self):
519 def __init__(self):
513 Operation.__init__(self)
520 Operation.__init__(self)
514
521
515 def __calculateCosDir(self, elev, azim):
522 def __calculateCosDir(self, elev, azim):
516 zen = (90 - elev)*numpy.pi/180
523 zen = (90 - elev)*numpy.pi/180
517 azim = azim*numpy.pi/180
524 azim = azim*numpy.pi/180
518 cosDirX = numpy.sqrt((1-numpy.cos(zen)**2)/((1+numpy.tan(azim)**2)))
525 cosDirX = numpy.sqrt((1-numpy.cos(zen)**2)/((1+numpy.tan(azim)**2)))
519 cosDirY = numpy.sqrt(1-numpy.cos(zen)**2-cosDirX**2)
526 cosDirY = numpy.sqrt(1-numpy.cos(zen)**2-cosDirX**2)
520
527
521 signX = numpy.sign(numpy.cos(azim))
528 signX = numpy.sign(numpy.cos(azim))
522 signY = numpy.sign(numpy.sin(azim))
529 signY = numpy.sign(numpy.sin(azim))
523
530
524 cosDirX = numpy.copysign(cosDirX, signX)
531 cosDirX = numpy.copysign(cosDirX, signX)
525 cosDirY = numpy.copysign(cosDirY, signY)
532 cosDirY = numpy.copysign(cosDirY, signY)
526 return cosDirX, cosDirY
533 return cosDirX, cosDirY
527
534
528 def __calculateAngles(self, theta_x, theta_y, azimuth):
535 def __calculateAngles(self, theta_x, theta_y, azimuth):
529
536
530 dir_cosw = numpy.sqrt(1-theta_x**2-theta_y**2)
537 dir_cosw = numpy.sqrt(1-theta_x**2-theta_y**2)
531 zenith_arr = numpy.arccos(dir_cosw)
538 zenith_arr = numpy.arccos(dir_cosw)
532 azimuth_arr = numpy.arctan2(theta_x,theta_y) + azimuth*math.pi/180
539 azimuth_arr = numpy.arctan2(theta_x,theta_y) + azimuth*math.pi/180
533
540
534 dir_cosu = numpy.sin(azimuth_arr)*numpy.sin(zenith_arr)
541 dir_cosu = numpy.sin(azimuth_arr)*numpy.sin(zenith_arr)
535 dir_cosv = numpy.cos(azimuth_arr)*numpy.sin(zenith_arr)
542 dir_cosv = numpy.cos(azimuth_arr)*numpy.sin(zenith_arr)
536
543
537 return azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw
544 return azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw
538
545
539 def __calculateMatA(self, dir_cosu, dir_cosv, dir_cosw, horOnly):
546 def __calculateMatA(self, dir_cosu, dir_cosv, dir_cosw, horOnly):
540
547
541 #
548 #
542 if horOnly:
549 if horOnly:
543 A = numpy.c_[dir_cosu,dir_cosv]
550 A = numpy.c_[dir_cosu,dir_cosv]
544 else:
551 else:
545 A = numpy.c_[dir_cosu,dir_cosv,dir_cosw]
552 A = numpy.c_[dir_cosu,dir_cosv,dir_cosw]
546 A = numpy.asmatrix(A)
553 A = numpy.asmatrix(A)
547 A1 = numpy.linalg.inv(A.transpose()*A)*A.transpose()
554 A1 = numpy.linalg.inv(A.transpose()*A)*A.transpose()
548
555
549 return A1
556 return A1
550
557
551 def __correctValues(self, heiRang, phi, velRadial, SNR):
558 def __correctValues(self, heiRang, phi, velRadial, SNR):
552 listPhi = phi.tolist()
559 listPhi = phi.tolist()
553 maxid = listPhi.index(max(listPhi))
560 maxid = listPhi.index(max(listPhi))
554 minid = listPhi.index(min(listPhi))
561 minid = listPhi.index(min(listPhi))
555
562
556 rango = range(len(phi))
563 rango = range(len(phi))
557 # rango = numpy.delete(rango,maxid)
564 # rango = numpy.delete(rango,maxid)
558
565
559 heiRang1 = heiRang*math.cos(phi[maxid])
566 heiRang1 = heiRang*math.cos(phi[maxid])
560 heiRangAux = heiRang*math.cos(phi[minid])
567 heiRangAux = heiRang*math.cos(phi[minid])
561 indOut = (heiRang1 < heiRangAux[0]).nonzero()
568 indOut = (heiRang1 < heiRangAux[0]).nonzero()
562 heiRang1 = numpy.delete(heiRang1,indOut)
569 heiRang1 = numpy.delete(heiRang1,indOut)
563
570
564 velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
571 velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
565 SNR1 = numpy.zeros([len(phi),len(heiRang1)])
572 SNR1 = numpy.zeros([len(phi),len(heiRang1)])
566
573
567 for i in rango:
574 for i in rango:
568 x = heiRang*math.cos(phi[i])
575 x = heiRang*math.cos(phi[i])
569 y1 = velRadial[i,:]
576 y1 = velRadial[i,:]
570 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
577 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
571
578
572 x1 = heiRang1
579 x1 = heiRang1
573 y11 = f1(x1)
580 y11 = f1(x1)
574
581
575 y2 = SNR[i,:]
582 y2 = SNR[i,:]
576 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
583 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
577 y21 = f2(x1)
584 y21 = f2(x1)
578
585
579 velRadial1[i,:] = y11
586 velRadial1[i,:] = y11
580 SNR1[i,:] = y21
587 SNR1[i,:] = y21
581
588
582 return heiRang1, velRadial1, SNR1
589 return heiRang1, velRadial1, SNR1
583
590
584 def __calculateVelUVW(self, A, velRadial):
591 def __calculateVelUVW(self, A, velRadial):
585
592
586 #Operacion Matricial
593 #Operacion Matricial
587 # velUVW = numpy.zeros((velRadial.shape[1],3))
594 # velUVW = numpy.zeros((velRadial.shape[1],3))
588 # for ind in range(velRadial.shape[1]):
595 # for ind in range(velRadial.shape[1]):
589 # velUVW[ind,:] = numpy.dot(A,velRadial[:,ind])
596 # velUVW[ind,:] = numpy.dot(A,velRadial[:,ind])
590 # velUVW = velUVW.transpose()
597 # velUVW = velUVW.transpose()
591 velUVW = numpy.zeros((A.shape[0],velRadial.shape[1]))
598 velUVW = numpy.zeros((A.shape[0],velRadial.shape[1]))
592 velUVW[:,:] = numpy.dot(A,velRadial)
599 velUVW[:,:] = numpy.dot(A,velRadial)
593
600
594
601
595 return velUVW
602 return velUVW
596
603
597 # def techniqueDBS(self, velRadial0, dirCosx, disrCosy, azimuth, correct, horizontalOnly, heiRang, SNR0):
604 # def techniqueDBS(self, velRadial0, dirCosx, disrCosy, azimuth, correct, horizontalOnly, heiRang, SNR0):
598 def techniqueDBS(self, velRadial0, heiRang, SNR0, kwargs):
605
606 def techniqueDBS(self, kwargs):
599 """
607 """
600 Function that implements Doppler Beam Swinging (DBS) technique.
608 Function that implements Doppler Beam Swinging (DBS) technique.
601
609
602 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
610 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
603 Direction correction (if necessary), Ranges and SNR
611 Direction correction (if necessary), Ranges and SNR
604
612
605 Output: Winds estimation (Zonal, Meridional and Vertical)
613 Output: Winds estimation (Zonal, Meridional and Vertical)
606
614
607 Parameters affected: Winds, height range, SNR
615 Parameters affected: Winds, height range, SNR
608 """
616 """
617 velRadial0 = kwargs['velRadial']
618 heiRang = kwargs['heightList']
619 SNR0 = kwargs['SNR']
609
620
610 if kwargs.has_key('dirCosx') and kwargs.has_key('dirCosy'):
621 if kwargs.has_key('dirCosx') and kwargs.has_key('dirCosy'):
611 theta_x = numpy.array(kwargs['dirCosx'])
622 theta_x = numpy.array(kwargs['dirCosx'])
612 theta_y = numpy.array(kwargs['dirCosy'])
623 theta_y = numpy.array(kwargs['dirCosy'])
613 else:
624 else:
614 elev = numpy.array(kwargs['elevation'])
625 elev = numpy.array(kwargs['elevation'])
615 azim = numpy.array(kwargs['azimuth'])
626 azim = numpy.array(kwargs['azimuth'])
616 theta_x, theta_y = self.__calculateCosDir(elev, azim)
627 theta_x, theta_y = self.__calculateCosDir(elev, azim)
617 azimuth = kwargs['correctAzimuth']
628 azimuth = kwargs['correctAzimuth']
618 if kwargs.has_key('horizontalOnly'):
629 if kwargs.has_key('horizontalOnly'):
619 horizontalOnly = kwargs['horizontalOnly']
630 horizontalOnly = kwargs['horizontalOnly']
620 else: horizontalOnly = False
631 else: horizontalOnly = False
621 if kwargs.has_key('correctFactor'):
632 if kwargs.has_key('correctFactor'):
622 correctFactor = kwargs['correctFactor']
633 correctFactor = kwargs['correctFactor']
623 else: correctFactor = 1
634 else: correctFactor = 1
624 if kwargs.has_key('channelList'):
635 if kwargs.has_key('channelList'):
625 channelList = kwargs['channelList']
636 channelList = kwargs['channelList']
626 if len(channelList) == 2:
637 if len(channelList) == 2:
627 horizontalOnly = True
638 horizontalOnly = True
628 arrayChannel = numpy.array(channelList)
639 arrayChannel = numpy.array(channelList)
629 param = param[arrayChannel,:,:]
640 param = param[arrayChannel,:,:]
630 theta_x = theta_x[arrayChannel]
641 theta_x = theta_x[arrayChannel]
631 theta_y = theta_y[arrayChannel]
642 theta_y = theta_y[arrayChannel]
632
643
633 azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(theta_x, theta_y, azimuth)
644 azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(theta_x, theta_y, azimuth)
634 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correctFactor*velRadial0, SNR0)
645 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correctFactor*velRadial0, SNR0)
635 A = self.__calculateMatA(dir_cosu, dir_cosv, dir_cosw, horizontalOnly)
646 A = self.__calculateMatA(dir_cosu, dir_cosv, dir_cosw, horizontalOnly)
636
647
637 #Calculo de Componentes de la velocidad con DBS
648 #Calculo de Componentes de la velocidad con DBS
638 winds = self.__calculateVelUVW(A,velRadial1)
649 winds = self.__calculateVelUVW(A,velRadial1)
639
650
640 return winds, heiRang1, SNR1
651 return winds, heiRang1, SNR1
641
652
642 def __calculateDistance(self, posx, posy, pairsCrossCorr, pairsList, pairs, azimuth = None):
653 def __calculateDistance(self, posx, posy, pairs_ccf, azimuth = None):
643
654
655 nPairs = len(pairs_ccf)
644 posx = numpy.asarray(posx)
656 posx = numpy.asarray(posx)
645 posy = numpy.asarray(posy)
657 posy = numpy.asarray(posy)
646
658
647 #Rotacion Inversa para alinear con el azimuth
659 #Rotacion Inversa para alinear con el azimuth
648 if azimuth!= None:
660 if azimuth!= None:
649 azimuth = azimuth*math.pi/180
661 azimuth = azimuth*math.pi/180
650 posx1 = posx*math.cos(azimuth) + posy*math.sin(azimuth)
662 posx1 = posx*math.cos(azimuth) + posy*math.sin(azimuth)
651 posy1 = -posx*math.sin(azimuth) + posy*math.cos(azimuth)
663 posy1 = -posx*math.sin(azimuth) + posy*math.cos(azimuth)
652 else:
664 else:
653 posx1 = posx
665 posx1 = posx
654 posy1 = posy
666 posy1 = posy
655
667
656 #Calculo de Distancias
668 #Calculo de Distancias
657 distx = numpy.zeros(pairsCrossCorr.size)
669 distx = numpy.zeros(nPairs)
658 disty = numpy.zeros(pairsCrossCorr.size)
670 disty = numpy.zeros(nPairs)
659 dist = numpy.zeros(pairsCrossCorr.size)
671 dist = numpy.zeros(nPairs)
660 ang = numpy.zeros(pairsCrossCorr.size)
672 ang = numpy.zeros(nPairs)
661
673
662 for i in range(pairsCrossCorr.size):
674 for i in range(nPairs):
663 distx[i] = posx1[pairsList[pairsCrossCorr[i]][1]] - posx1[pairsList[pairsCrossCorr[i]][0]]
675 distx[i] = posx1[pairs_ccf[i][1]] - posx1[pairs_ccf[i][0]]
664 disty[i] = posy1[pairsList[pairsCrossCorr[i]][1]] - posy1[pairsList[pairsCrossCorr[i]][0]]
676 disty[i] = posy1[pairs_ccf[i][1]] - posy1[pairs_ccf[i][0]]
665 dist[i] = numpy.sqrt(distx[i]**2 + disty[i]**2)
677 dist[i] = numpy.sqrt(distx[i]**2 + disty[i]**2)
666 ang[i] = numpy.arctan2(disty[i],distx[i])
678 ang[i] = numpy.arctan2(disty[i],distx[i])
667 #Calculo de Matrices
668 nPairs = len(pairs)
669 ang1 = numpy.zeros((nPairs, 2, 1))
670 dist1 = numpy.zeros((nPairs, 2, 1))
671
679
672 for j in range(nPairs):
680 return distx, disty, dist, ang
673 dist1[j,0,0] = dist[pairs[j][0]]
681 #Calculo de Matrices
674 dist1[j,1,0] = dist[pairs[j][1]]
682 # nPairs = len(pairs)
675 ang1[j,0,0] = ang[pairs[j][0]]
683 # ang1 = numpy.zeros((nPairs, 2, 1))
676 ang1[j,1,0] = ang[pairs[j][1]]
684 # dist1 = numpy.zeros((nPairs, 2, 1))
677
685 #
678 return distx,disty, dist1,ang1
686 # for j in range(nPairs):
687 # dist1[j,0,0] = dist[pairs[j][0]]
688 # dist1[j,1,0] = dist[pairs[j][1]]
689 # ang1[j,0,0] = ang[pairs[j][0]]
690 # ang1[j,1,0] = ang[pairs[j][1]]
691 #
692 # return distx,disty, dist1,ang1
693
679
694
680 def __calculateVelVer(self, phase, lagTRange, _lambda):
695 def __calculateVelVer(self, phase, lagTRange, _lambda):
681
696
682 Ts = lagTRange[1] - lagTRange[0]
697 Ts = lagTRange[1] - lagTRange[0]
683 velW = -_lambda*phase/(4*math.pi*Ts)
698 velW = -_lambda*phase/(4*math.pi*Ts)
684
699
685 return velW
700 return velW
686
701
687 def __calculateVelHorDir(self, dist, tau1, tau2, ang):
702 def __calculateVelHorDir(self, dist, tau1, tau2, ang):
688 nPairs = tau1.shape[0]
703 nPairs = tau1.shape[0]
689 vel = numpy.zeros((nPairs,3,tau1.shape[2]))
704 nHeights = tau1.shape[1]
705 vel = numpy.zeros((nPairs,3,nHeights))
706 dist1 = numpy.reshape(dist, (dist.size,1))
690
707
691 angCos = numpy.cos(ang)
708 angCos = numpy.cos(ang)
692 angSin = numpy.sin(ang)
709 angSin = numpy.sin(ang)
693
710
694 vel0 = dist*tau1/(2*tau2**2)
711 vel0 = dist1*tau1/(2*tau2**2)
695 vel[:,0,:] = (vel0*angCos).sum(axis = 1)
712 vel[:,0,:] = (vel0*angCos).sum(axis = 1)
696 vel[:,1,:] = (vel0*angSin).sum(axis = 1)
713 vel[:,1,:] = (vel0*angSin).sum(axis = 1)
697
714
698 ind = numpy.where(numpy.isinf(vel))
715 ind = numpy.where(numpy.isinf(vel))
699 vel[ind] = numpy.nan
716 vel[ind] = numpy.nan
700
717
701 return vel
718 return vel
702
719
703 def __getPairsAutoCorr(self, pairsList, nChannels):
720 # def __getPairsAutoCorr(self, pairsList, nChannels):
704
721 #
705 pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
722 # pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
706
723 #
707 for l in range(len(pairsList)):
724 # for l in range(len(pairsList)):
708 firstChannel = pairsList[l][0]
725 # firstChannel = pairsList[l][0]
709 secondChannel = pairsList[l][1]
726 # secondChannel = pairsList[l][1]
710
727 #
711 #Obteniendo pares de Autocorrelacion
728 # #Obteniendo pares de Autocorrelacion
712 if firstChannel == secondChannel:
729 # if firstChannel == secondChannel:
713 pairsAutoCorr[firstChannel] = int(l)
730 # pairsAutoCorr[firstChannel] = int(l)
714
731 #
715 pairsAutoCorr = pairsAutoCorr.astype(int)
732 # pairsAutoCorr = pairsAutoCorr.astype(int)
716
733 #
717 pairsCrossCorr = range(len(pairsList))
734 # pairsCrossCorr = range(len(pairsList))
718 pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
735 # pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
719
736 #
720 return pairsAutoCorr, pairsCrossCorr
737 # return pairsAutoCorr, pairsCrossCorr
721
738
722 def techniqueSA(self, pairsSelected, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, lagTRange, correctFactor):
739 # def techniqueSA(self, pairsSelected, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, lagTRange, correctFactor):
740 def techniqueSA(self, kwargs):
741
723 """
742 """
724 Function that implements Spaced Antenna (SA) technique.
743 Function that implements Spaced Antenna (SA) technique.
725
744
726 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
745 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
727 Direction correction (if necessary), Ranges and SNR
746 Direction correction (if necessary), Ranges and SNR
728
747
729 Output: Winds estimation (Zonal, Meridional and Vertical)
748 Output: Winds estimation (Zonal, Meridional and Vertical)
730
749
731 Parameters affected: Winds
750 Parameters affected: Winds
732 """
751 """
733 #Cross Correlation pairs obtained
752 position_x = kwargs['positionX']
734 pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)
753 position_y = kwargs['positionY']
735 pairsArray = numpy.array(pairsList)[pairsCrossCorr]
754 azimuth = kwargs['azimuth']
736 pairsSelArray = numpy.array(pairsSelected)
755
737 pairs = []
756 if kwargs.has_key('correctFactor'):
757 correctFactor = kwargs['correctFactor']
758 else:
759 correctFactor = 1
760
761 groupList = kwargs['groupList']
762 pairs_ccf = groupList[1]
763 tau = kwargs['tau']
764 _lambda = kwargs['_lambda']
738
765
739 #Wind estimation pairs obtained
766 #Cross Correlation pairs obtained
740 for i in range(pairsSelArray.shape[0]/2):
767 # pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairssList, nChannels)
741 ind1 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i], axis = 1))[0][0]
768 # pairsArray = numpy.array(pairsList)[pairsCrossCorr]
742 ind2 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i + 1], axis = 1))[0][0]
769 # pairsSelArray = numpy.array(pairsSelected)
743 pairs.append((ind1,ind2))
770 # pairs = []
771 #
772 # #Wind estimation pairs obtained
773 # for i in range(pairsSelArray.shape[0]/2):
774 # ind1 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i], axis = 1))[0][0]
775 # ind2 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i + 1], axis = 1))[0][0]
776 # pairs.append((ind1,ind2))
744
777
745 indtau = tau.shape[0]/2
778 indtau = tau.shape[0]/2
746 tau1 = tau[:indtau,:]
779 tau1 = tau[:indtau,:]
747 tau2 = tau[indtau:-1,:]
780 tau2 = tau[indtau:-1,:]
748 tau1 = tau1[pairs,:]
781 # tau1 = tau1[pairs,:]
749 tau2 = tau2[pairs,:]
782 # tau2 = tau2[pairs,:]
750 phase1 = tau[-1,:]
783 phase1 = tau[-1,:]
751
784
752 #---------------------------------------------------------------------
785 #---------------------------------------------------------------------
753 #Metodo Directo
786 #Metodo Directo
754 distx, disty, dist, ang = self.__calculateDistance(position_x, position_y, pairsCrossCorr, pairsList, pairs,azimuth)
787 distx, disty, dist, ang = self.__calculateDistance(position_x, position_y, pairs_ccf,azimuth)
755 winds = self.__calculateVelHorDir(dist, tau1, tau2, ang)
788 winds = self.__calculateVelHorDir(dist, tau1, tau2, ang)
756 winds = stats.nanmean(winds, axis=0)
789 winds = stats.nanmean(winds, axis=0)
757 #---------------------------------------------------------------------
790 #---------------------------------------------------------------------
758 #Metodo General
791 #Metodo General
759 # distx, disty, dist = self.calculateDistance(position_x,position_y,pairsCrossCorr, pairsList, azimuth)
792 # distx, disty, dist = self.calculateDistance(position_x,position_y,pairsCrossCorr, pairsList, azimuth)
760 # #Calculo Coeficientes de Funcion de Correlacion
793 # #Calculo Coeficientes de Funcion de Correlacion
761 # F,G,A,B,H = self.calculateCoef(tau1,tau2,distx,disty,n)
794 # F,G,A,B,H = self.calculateCoef(tau1,tau2,distx,disty,n)
762 # #Calculo de Velocidades
795 # #Calculo de Velocidades
763 # winds = self.calculateVelUV(F,G,A,B,H)
796 # winds = self.calculateVelUV(F,G,A,B,H)
764
797
765 #---------------------------------------------------------------------
798 #---------------------------------------------------------------------
766 winds[2,:] = self.__calculateVelVer(phase1, lagTRange, _lambda)
799 winds[2,:] = self.__calculateVelVer(phase1, lagTRange, _lambda)
767 winds = correctFactor*winds
800 winds = correctFactor*winds
768 return winds
801 return winds
769
802
770 def __checkTime(self, currentTime, paramInterval, outputInterval):
803 def __checkTime(self, currentTime, paramInterval, outputInterval):
771
804
772 dataTime = currentTime + paramInterval
805 dataTime = currentTime + paramInterval
773 deltaTime = dataTime - self.__initime
806 deltaTime = dataTime - self.__initime
774
807
775 if deltaTime >= outputInterval or deltaTime < 0:
808 if deltaTime >= outputInterval or deltaTime < 0:
776 self.__dataReady = True
809 self.__dataReady = True
777 return
810 return
778
811
779 def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax):
812 def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax):
780 '''
813 '''
781 Function that implements winds estimation technique with detected meteors.
814 Function that implements winds estimation technique with detected meteors.
782
815
783 Input: Detected meteors, Minimum meteor quantity to wind estimation
816 Input: Detected meteors, Minimum meteor quantity to wind estimation
784
817
785 Output: Winds estimation (Zonal and Meridional)
818 Output: Winds estimation (Zonal and Meridional)
786
819
787 Parameters affected: Winds
820 Parameters affected: Winds
788 '''
821 '''
789 # print arrayMeteor.shape
822 # print arrayMeteor.shape
790 #Settings
823 #Settings
791 nInt = (heightMax - heightMin)/2
824 nInt = (heightMax - heightMin)/2
792 # print nInt
825 # print nInt
793 nInt = int(nInt)
826 nInt = int(nInt)
794 # print nInt
827 # print nInt
795 winds = numpy.zeros((2,nInt))*numpy.nan
828 winds = numpy.zeros((2,nInt))*numpy.nan
796
829
797 #Filter errors
830 #Filter errors
798 error = numpy.where(arrayMeteor[:,-1] == 0)[0]
831 error = numpy.where(arrayMeteor[:,-1] == 0)[0]
799 finalMeteor = arrayMeteor[error,:]
832 finalMeteor = arrayMeteor[error,:]
800
833
801 #Meteor Histogram
834 #Meteor Histogram
802 finalHeights = finalMeteor[:,2]
835 finalHeights = finalMeteor[:,2]
803 hist = numpy.histogram(finalHeights, bins = nInt, range = (heightMin,heightMax))
836 hist = numpy.histogram(finalHeights, bins = nInt, range = (heightMin,heightMax))
804 nMeteorsPerI = hist[0]
837 nMeteorsPerI = hist[0]
805 heightPerI = hist[1]
838 heightPerI = hist[1]
806
839
807 #Sort of meteors
840 #Sort of meteors
808 indSort = finalHeights.argsort()
841 indSort = finalHeights.argsort()
809 finalMeteor2 = finalMeteor[indSort,:]
842 finalMeteor2 = finalMeteor[indSort,:]
810
843
811 # Calculating winds
844 # Calculating winds
812 ind1 = 0
845 ind1 = 0
813 ind2 = 0
846 ind2 = 0
814
847
815 for i in range(nInt):
848 for i in range(nInt):
816 nMet = nMeteorsPerI[i]
849 nMet = nMeteorsPerI[i]
817 ind1 = ind2
850 ind1 = ind2
818 ind2 = ind1 + nMet
851 ind2 = ind1 + nMet
819
852
820 meteorAux = finalMeteor2[ind1:ind2,:]
853 meteorAux = finalMeteor2[ind1:ind2,:]
821
854
822 if meteorAux.shape[0] >= meteorThresh:
855 if meteorAux.shape[0] >= meteorThresh:
823 vel = meteorAux[:, 6]
856 vel = meteorAux[:, 6]
824 zen = meteorAux[:, 4]*numpy.pi/180
857 zen = meteorAux[:, 4]*numpy.pi/180
825 azim = meteorAux[:, 3]*numpy.pi/180
858 azim = meteorAux[:, 3]*numpy.pi/180
826
859
827 n = numpy.cos(zen)
860 n = numpy.cos(zen)
828 # m = (1 - n**2)/(1 - numpy.tan(azim)**2)
861 # m = (1 - n**2)/(1 - numpy.tan(azim)**2)
829 # l = m*numpy.tan(azim)
862 # l = m*numpy.tan(azim)
830 l = numpy.sin(zen)*numpy.sin(azim)
863 l = numpy.sin(zen)*numpy.sin(azim)
831 m = numpy.sin(zen)*numpy.cos(azim)
864 m = numpy.sin(zen)*numpy.cos(azim)
832
865
833 A = numpy.vstack((l, m)).transpose()
866 A = numpy.vstack((l, m)).transpose()
834 A1 = numpy.dot(numpy.linalg.inv( numpy.dot(A.transpose(),A) ),A.transpose())
867 A1 = numpy.dot(numpy.linalg.inv( numpy.dot(A.transpose(),A) ),A.transpose())
835 windsAux = numpy.dot(A1, vel)
868 windsAux = numpy.dot(A1, vel)
836
869
837 winds[0,i] = windsAux[0]
870 winds[0,i] = windsAux[0]
838 winds[1,i] = windsAux[1]
871 winds[1,i] = windsAux[1]
839
872
840 return winds, heightPerI[:-1]
873 return winds, heightPerI[:-1]
841
874
842 def techniqueNSM_SA(self, **kwargs):
875 def techniqueNSM_SA(self, **kwargs):
843 metArray = kwargs['metArray']
876 metArray = kwargs['metArray']
844 heightList = kwargs['heightList']
877 heightList = kwargs['heightList']
845 timeList = kwargs['timeList']
878 timeList = kwargs['timeList']
846
879
847 rx_location = kwargs['rx_location']
880 rx_location = kwargs['rx_location']
848 groupList = kwargs['groupList']
881 groupList = kwargs['groupList']
849 azimuth = kwargs['azimuth']
882 azimuth = kwargs['azimuth']
850 dfactor = kwargs['dfactor']
883 dfactor = kwargs['dfactor']
851 k = kwargs['k']
884 k = kwargs['k']
852
885
853 azimuth1, dist = self.__calculateAzimuth1(rx_location, groupList, azimuth)
886 azimuth1, dist = self.__calculateAzimuth1(rx_location, groupList, azimuth)
854 d = dist*dfactor
887 d = dist*dfactor
855 #Phase calculation
888 #Phase calculation
856 metArray1 = self.__getPhaseSlope(metArray, heightList, timeList)
889 metArray1 = self.__getPhaseSlope(metArray, heightList, timeList)
857
890
858 metArray1[:,-2] = metArray1[:,-2]*metArray1[:,2]*1000/(k*d[metArray1[:,1].astype(int)]) #angles into velocities
891 metArray1[:,-2] = metArray1[:,-2]*metArray1[:,2]*1000/(k*d[metArray1[:,1].astype(int)]) #angles into velocities
859
892
860 velEst = numpy.zeros((heightList.size,2))*numpy.nan
893 velEst = numpy.zeros((heightList.size,2))*numpy.nan
861 azimuth1 = azimuth1*numpy.pi/180
894 azimuth1 = azimuth1*numpy.pi/180
862
895
863 for i in range(heightList.size):
896 for i in range(heightList.size):
864 h = heightList[i]
897 h = heightList[i]
865 indH = numpy.where((metArray1[:,2] == h)&(numpy.abs(metArray1[:,-2]) < 100))[0]
898 indH = numpy.where((metArray1[:,2] == h)&(numpy.abs(metArray1[:,-2]) < 100))[0]
866 metHeight = metArray1[indH,:]
899 metHeight = metArray1[indH,:]
867 if metHeight.shape[0] >= 2:
900 if metHeight.shape[0] >= 2:
868 velAux = numpy.asmatrix(metHeight[:,-2]).T #Radial Velocities
901 velAux = numpy.asmatrix(metHeight[:,-2]).T #Radial Velocities
869 iazim = metHeight[:,1].astype(int)
902 iazim = metHeight[:,1].astype(int)
870 azimAux = numpy.asmatrix(azimuth1[iazim]).T #Azimuths
903 azimAux = numpy.asmatrix(azimuth1[iazim]).T #Azimuths
871 A = numpy.hstack((numpy.cos(azimAux),numpy.sin(azimAux)))
904 A = numpy.hstack((numpy.cos(azimAux),numpy.sin(azimAux)))
872 A = numpy.asmatrix(A)
905 A = numpy.asmatrix(A)
873 A1 = numpy.linalg.pinv(A.transpose()*A)*A.transpose()
906 A1 = numpy.linalg.pinv(A.transpose()*A)*A.transpose()
874 velHor = numpy.dot(A1,velAux)
907 velHor = numpy.dot(A1,velAux)
875
908
876 velEst[i,:] = numpy.squeeze(velHor)
909 velEst[i,:] = numpy.squeeze(velHor)
877 return velEst
910 return velEst
878
911
879 def __getPhaseSlope(self, metArray, heightList, timeList):
912 def __getPhaseSlope(self, metArray, heightList, timeList):
880 meteorList = []
913 meteorList = []
881 #utctime sec1 height SNR velRad ph0 ph1 ph2 coh0 coh1 coh2
914 #utctime sec1 height SNR velRad ph0 ph1 ph2 coh0 coh1 coh2
882 #Putting back together the meteor matrix
915 #Putting back together the meteor matrix
883 utctime = metArray[:,0]
916 utctime = metArray[:,0]
884 uniqueTime = numpy.unique(utctime)
917 uniqueTime = numpy.unique(utctime)
885
918
886 phaseDerThresh = 0.5
919 phaseDerThresh = 0.5
887 ippSeconds = timeList[1] - timeList[0]
920 ippSeconds = timeList[1] - timeList[0]
888 sec = numpy.where(timeList>1)[0][0]
921 sec = numpy.where(timeList>1)[0][0]
889 nPairs = metArray.shape[1] - 6
922 nPairs = metArray.shape[1] - 6
890 nHeights = len(heightList)
923 nHeights = len(heightList)
891
924
892 for t in uniqueTime:
925 for t in uniqueTime:
893 metArray1 = metArray[utctime==t,:]
926 metArray1 = metArray[utctime==t,:]
894 # phaseDerThresh = numpy.pi/4 #reducir Phase thresh
927 # phaseDerThresh = numpy.pi/4 #reducir Phase thresh
895 tmet = metArray1[:,1].astype(int)
928 tmet = metArray1[:,1].astype(int)
896 hmet = metArray1[:,2].astype(int)
929 hmet = metArray1[:,2].astype(int)
897
930
898 metPhase = numpy.zeros((nPairs, heightList.size, timeList.size - 1))
931 metPhase = numpy.zeros((nPairs, heightList.size, timeList.size - 1))
899 metPhase[:,:] = numpy.nan
932 metPhase[:,:] = numpy.nan
900 metPhase[:,hmet,tmet] = metArray1[:,6:].T
933 metPhase[:,hmet,tmet] = metArray1[:,6:].T
901
934
902 #Delete short trails
935 #Delete short trails
903 metBool = ~numpy.isnan(metPhase[0,:,:])
936 metBool = ~numpy.isnan(metPhase[0,:,:])
904 heightVect = numpy.sum(metBool, axis = 1)
937 heightVect = numpy.sum(metBool, axis = 1)
905 metBool[heightVect<sec,:] = False
938 metBool[heightVect<sec,:] = False
906 metPhase[:,heightVect<sec,:] = numpy.nan
939 metPhase[:,heightVect<sec,:] = numpy.nan
907
940
908 #Derivative
941 #Derivative
909 metDer = numpy.abs(metPhase[:,:,1:] - metPhase[:,:,:-1])
942 metDer = numpy.abs(metPhase[:,:,1:] - metPhase[:,:,:-1])
910 phDerAux = numpy.dstack((numpy.full((nPairs,nHeights,1), False, dtype=bool),metDer > phaseDerThresh))
943 phDerAux = numpy.dstack((numpy.full((nPairs,nHeights,1), False, dtype=bool),metDer > phaseDerThresh))
911 metPhase[phDerAux] = numpy.nan
944 metPhase[phDerAux] = numpy.nan
912
945
913 #--------------------------METEOR DETECTION -----------------------------------------
946 #--------------------------METEOR DETECTION -----------------------------------------
914 indMet = numpy.where(numpy.any(metBool,axis=1))[0]
947 indMet = numpy.where(numpy.any(metBool,axis=1))[0]
915
948
916 for p in numpy.arange(nPairs):
949 for p in numpy.arange(nPairs):
917 phase = metPhase[p,:,:]
950 phase = metPhase[p,:,:]
918 phDer = metDer[p,:,:]
951 phDer = metDer[p,:,:]
919
952
920 for h in indMet:
953 for h in indMet:
921 height = heightList[h]
954 height = heightList[h]
922 phase1 = phase[h,:] #82
955 phase1 = phase[h,:] #82
923 phDer1 = phDer[h,:]
956 phDer1 = phDer[h,:]
924
957
925 phase1[~numpy.isnan(phase1)] = numpy.unwrap(phase1[~numpy.isnan(phase1)]) #Unwrap
958 phase1[~numpy.isnan(phase1)] = numpy.unwrap(phase1[~numpy.isnan(phase1)]) #Unwrap
926
959
927 indValid = numpy.where(~numpy.isnan(phase1))[0]
960 indValid = numpy.where(~numpy.isnan(phase1))[0]
928 initMet = indValid[0]
961 initMet = indValid[0]
929 endMet = 0
962 endMet = 0
930
963
931 for i in range(len(indValid)-1):
964 for i in range(len(indValid)-1):
932
965
933 #Time difference
966 #Time difference
934 inow = indValid[i]
967 inow = indValid[i]
935 inext = indValid[i+1]
968 inext = indValid[i+1]
936 idiff = inext - inow
969 idiff = inext - inow
937 #Phase difference
970 #Phase difference
938 phDiff = numpy.abs(phase1[inext] - phase1[inow])
971 phDiff = numpy.abs(phase1[inext] - phase1[inow])
939
972
940 if idiff>sec or phDiff>numpy.pi/4 or inext==indValid[-1]: #End of Meteor
973 if idiff>sec or phDiff>numpy.pi/4 or inext==indValid[-1]: #End of Meteor
941 sizeTrail = inow - initMet + 1
974 sizeTrail = inow - initMet + 1
942 if sizeTrail>3*sec: #Too short meteors
975 if sizeTrail>3*sec: #Too short meteors
943 x = numpy.arange(initMet,inow+1)*ippSeconds
976 x = numpy.arange(initMet,inow+1)*ippSeconds
944 y = phase1[initMet:inow+1]
977 y = phase1[initMet:inow+1]
945 ynnan = ~numpy.isnan(y)
978 ynnan = ~numpy.isnan(y)
946 x = x[ynnan]
979 x = x[ynnan]
947 y = y[ynnan]
980 y = y[ynnan]
948 slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
981 slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
949 ylin = x*slope + intercept
982 ylin = x*slope + intercept
950 rsq = r_value**2
983 rsq = r_value**2
951 if rsq > 0.5:
984 if rsq > 0.5:
952 vel = slope#*height*1000/(k*d)
985 vel = slope#*height*1000/(k*d)
953 estAux = numpy.array([utctime,p,height, vel, rsq])
986 estAux = numpy.array([utctime,p,height, vel, rsq])
954 meteorList.append(estAux)
987 meteorList.append(estAux)
955 initMet = inext
988 initMet = inext
956 metArray2 = numpy.array(meteorList)
989 metArray2 = numpy.array(meteorList)
957
990
958 return metArray2
991 return metArray2
959
992
960 def __calculateAzimuth1(self, rx_location, pairslist, azimuth0):
993 def __calculateAzimuth1(self, rx_location, pairslist, azimuth0):
961
994
962 azimuth1 = numpy.zeros(len(pairslist))
995 azimuth1 = numpy.zeros(len(pairslist))
963 dist = numpy.zeros(len(pairslist))
996 dist = numpy.zeros(len(pairslist))
964
997
965 for i in range(len(rx_location)):
998 for i in range(len(rx_location)):
966 ch0 = pairslist[i][0]
999 ch0 = pairslist[i][0]
967 ch1 = pairslist[i][1]
1000 ch1 = pairslist[i][1]
968
1001
969 diffX = rx_location[ch0][0] - rx_location[ch1][0]
1002 diffX = rx_location[ch0][0] - rx_location[ch1][0]
970 diffY = rx_location[ch0][1] - rx_location[ch1][1]
1003 diffY = rx_location[ch0][1] - rx_location[ch1][1]
971 azimuth1[i] = numpy.arctan2(diffY,diffX)*180/numpy.pi
1004 azimuth1[i] = numpy.arctan2(diffY,diffX)*180/numpy.pi
972 dist[i] = numpy.sqrt(diffX**2 + diffY**2)
1005 dist[i] = numpy.sqrt(diffX**2 + diffY**2)
973
1006
974 azimuth1 -= azimuth0
1007 azimuth1 -= azimuth0
975 return azimuth1, dist
1008 return azimuth1, dist
976
1009
977 def techniqueNSM_DBS(self, **kwargs):
1010 def techniqueNSM_DBS(self, **kwargs):
978 metArray = kwargs['metArray']
1011 metArray = kwargs['metArray']
979 heightList = kwargs['heightList']
1012 heightList = kwargs['heightList']
980 timeList = kwargs['timeList']
1013 timeList = kwargs['timeList']
981 zenithList = kwargs['zenithList']
1014 zenithList = kwargs['zenithList']
982 nChan = numpy.max(cmet) + 1
1015 nChan = numpy.max(cmet) + 1
983 nHeights = len(heightList)
1016 nHeights = len(heightList)
984
1017
985 utctime = metArray[:,0]
1018 utctime = metArray[:,0]
986 cmet = metArray[:,1]
1019 cmet = metArray[:,1]
987 hmet = metArray1[:,3].astype(int)
1020 hmet = metArray1[:,3].astype(int)
988 h1met = heightList[hmet]*zenithList[cmet]
1021 h1met = heightList[hmet]*zenithList[cmet]
989 vmet = metArray1[:,5]
1022 vmet = metArray1[:,5]
990
1023
991 for i in range(nHeights - 1):
1024 for i in range(nHeights - 1):
992 hmin = heightList[i]
1025 hmin = heightList[i]
993 hmax = heightList[i + 1]
1026 hmax = heightList[i + 1]
994
1027
995 vthisH = vmet[(h1met>=hmin) & (h1met<hmax)]
1028 vthisH = vmet[(h1met>=hmin) & (h1met<hmax)]
996
1029
997
1030
998
1031
999 return data_output
1032 return data_output
1000
1033
1001 def run(self, dataOut, technique, **kwargs):
1034 def run(self, dataOut, technique, **kwargs):
1002
1035
1003 param = dataOut.data_param
1036 param = dataOut.data_param
1004 if dataOut.abscissaList != None:
1037 if dataOut.abscissaList != None:
1005 absc = dataOut.abscissaList[:-1]
1038 absc = dataOut.abscissaList[:-1]
1006 noise = dataOut.noise
1039 noise = dataOut.noise
1007 heightList = dataOut.heightList
1040 heightList = dataOut.heightList
1008 SNR = dataOut.data_SNR
1041 SNR = dataOut.data_SNR
1009
1042
1010 if technique == 'DBS':
1043 if technique == 'DBS':
1011 # if kwargs.has_key('dirCosx') and kwargs.has_key('dirCosy'):
1044
1012 # theta_x = numpy.array(kwargs['dirCosx'])
1045 kwargs['velRadial'] = param[:,1,:] #Radial velocity
1013 # theta_y = numpy.array(kwargs['dirCosy'])
1046 kwargs['heightList'] = heightList
1014 # else:
1047 kwargs['SNR'] = SNR
1015 # elev = numpy.array(kwargs['elevation'])
1048
1016 # azim = numpy.array(kwargs['azimuth'])
1049 dataOut.data_output, dataOut.heightList, dataOut.data_SNR = self.techniqueDBS(kwargs) #DBS Function
1017 # theta_x, theta_y = self.__calculateCosDir(elev, azim)
1018 # azimuth = kwargs['correctAzimuth']
1019 # if kwargs.has_key('horizontalOnly'):
1020 # horizontalOnly = kwargs['horizontalOnly']
1021 # else: horizontalOnly = False
1022 # if kwargs.has_key('correctFactor'):
1023 # correctFactor = kwargs['correctFactor']
1024 # else: correctFactor = 1
1025 # if kwargs.has_key('channelList'):
1026 # channelList = kwargs['channelList']
1027 # if len(channelList) == 2:
1028 # horizontalOnly = True
1029 # arrayChannel = numpy.array(channelList)
1030 # param = param[arrayChannel,:,:]
1031 # theta_x = theta_x[arrayChannel]
1032 # theta_y = theta_y[arrayChannel]
1033
1034 velRadial0 = param[:,1,:] #Radial velocity
1035 # dataOut.data_output, dataOut.heightList, dataOut.data_SNR = self.techniqueDBS(velRadial0, theta_x, theta_y, azimuth, correctFactor, horizontalOnly, heightList, SNR) #DBS Function
1036 dataOut.data_output, dataOut.heightList, dataOut.data_SNR = self.techniqueDBS(velRadial0, heightList, SNR, kwargs) #DBS Function
1037 dataOut.utctimeInit = dataOut.utctime
1050 dataOut.utctimeInit = dataOut.utctime
1038 dataOut.outputInterval = dataOut.paramInterval
1051 dataOut.outputInterval = dataOut.paramInterval
1039
1052
1040 elif technique == 'SA':
1053 elif technique == 'SA':
1041
1054
1042 #Parameters
1055 #Parameters
1043 position_x = kwargs['positionX']
1056 # position_x = kwargs['positionX']
1044 position_y = kwargs['positionY']
1057 # position_y = kwargs['positionY']
1045 azimuth = kwargs['azimuth']
1058 # azimuth = kwargs['azimuth']
1046
1059 #
1047 if kwargs.has_key('crosspairsList'):
1060 # if kwargs.has_key('crosspairsList'):
1048 pairs = kwargs['crosspairsList']
1061 # pairs = kwargs['crosspairsList']
1049 else:
1062 # else:
1050 pairs = None
1063 # pairs = None
1051
1064 #
1052 if kwargs.has_key('correctFactor'):
1065 # if kwargs.has_key('correctFactor'):
1053 correctFactor = kwargs['correctFactor']
1066 # correctFactor = kwargs['correctFactor']
1054 else:
1067 # else:
1055 correctFactor = 1
1068 # correctFactor = 1
1056
1069
1057 tau = dataOut.data_param
1070 # tau = dataOut.data_param
1058 _lambda = dataOut.C/dataOut.frequency
1071 # _lambda = dataOut.C/dataOut.frequency
1059 pairsList = dataOut.groupList
1072 # pairsList = dataOut.groupList
1060 nChannels = dataOut.nChannels
1073 # nChannels = dataOut.nChannels
1061
1074
1062 dataOut.data_output = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor)
1075 kwargs['groupList'] = dataOut.groupList
1076 kwargs['tau'] = dataOut.data_param
1077 kwargs['_lambda'] = dataOut.C/dataOut.frequency
1078 # dataOut.data_output = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor)
1079 dataOut.data_output = self.techniqueSA(kwargs)
1063 dataOut.utctimeInit = dataOut.utctime
1080 dataOut.utctimeInit = dataOut.utctime
1064 dataOut.outputInterval = dataOut.timeInterval
1081 dataOut.outputInterval = dataOut.timeInterval
1065
1082
1066 elif technique == 'Meteors':
1083 elif technique == 'Meteors':
1067 dataOut.flagNoData = True
1084 dataOut.flagNoData = True
1068 self.__dataReady = False
1085 self.__dataReady = False
1069
1086
1070 if kwargs.has_key('nHours'):
1087 if kwargs.has_key('nHours'):
1071 nHours = kwargs['nHours']
1088 nHours = kwargs['nHours']
1072 else:
1089 else:
1073 nHours = 1
1090 nHours = 1
1074
1091
1075 if kwargs.has_key('meteorsPerBin'):
1092 if kwargs.has_key('meteorsPerBin'):
1076 meteorThresh = kwargs['meteorsPerBin']
1093 meteorThresh = kwargs['meteorsPerBin']
1077 else:
1094 else:
1078 meteorThresh = 6
1095 meteorThresh = 6
1079
1096
1080 if kwargs.has_key('hmin'):
1097 if kwargs.has_key('hmin'):
1081 hmin = kwargs['hmin']
1098 hmin = kwargs['hmin']
1082 else: hmin = 70
1099 else: hmin = 70
1083 if kwargs.has_key('hmax'):
1100 if kwargs.has_key('hmax'):
1084 hmax = kwargs['hmax']
1101 hmax = kwargs['hmax']
1085 else: hmax = 110
1102 else: hmax = 110
1086
1103
1087 dataOut.outputInterval = nHours*3600
1104 dataOut.outputInterval = nHours*3600
1088
1105
1089 if self.__isConfig == False:
1106 if self.__isConfig == False:
1090 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
1107 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
1091 #Get Initial LTC time
1108 #Get Initial LTC time
1092 self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
1109 self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
1093 self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
1110 self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
1094
1111
1095 self.__isConfig = True
1112 self.__isConfig = True
1096
1113
1097 if self.__buffer == None:
1114 if self.__buffer == None:
1098 self.__buffer = dataOut.data_param
1115 self.__buffer = dataOut.data_param
1099 self.__firstdata = copy.copy(dataOut)
1116 self.__firstdata = copy.copy(dataOut)
1100
1117
1101 else:
1118 else:
1102 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
1119 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
1103
1120
1104 self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
1121 self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
1105
1122
1106 if self.__dataReady:
1123 if self.__dataReady:
1107 dataOut.utctimeInit = self.__initime
1124 dataOut.utctimeInit = self.__initime
1108
1125
1109 self.__initime += dataOut.outputInterval #to erase time offset
1126 self.__initime += dataOut.outputInterval #to erase time offset
1110
1127
1111 dataOut.data_output, dataOut.heightList = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax)
1128 dataOut.data_output, dataOut.heightList = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax)
1112 dataOut.flagNoData = False
1129 dataOut.flagNoData = False
1113 self.__buffer = None
1130 self.__buffer = None
1114
1131
1115 elif technique == 'Meteors1':
1132 elif technique == 'Meteors1':
1116 dataOut.flagNoData = True
1133 dataOut.flagNoData = True
1117 self.__dataReady = False
1134 self.__dataReady = False
1118
1135
1119 if kwargs.has_key('nMins'):
1136 if kwargs.has_key('nMins'):
1120 nMins = kwargs['nMins']
1137 nMins = kwargs['nMins']
1121 else: nMins = 20
1138 else: nMins = 20
1122 if kwargs.has_key('rx_location'):
1139 if kwargs.has_key('rx_location'):
1123 rx_location = kwargs['rx_location']
1140 rx_location = kwargs['rx_location']
1124 else: rx_location = [(0,1),(1,1),(1,0)]
1141 else: rx_location = [(0,1),(1,1),(1,0)]
1125 if kwargs.has_key('azimuth'):
1142 if kwargs.has_key('azimuth'):
1126 azimuth = kwargs['azimuth']
1143 azimuth = kwargs['azimuth']
1127 else: azimuth = 51
1144 else: azimuth = 51
1128 if kwargs.has_key('dfactor'):
1145 if kwargs.has_key('dfactor'):
1129 dfactor = kwargs['dfactor']
1146 dfactor = kwargs['dfactor']
1130 if kwargs.has_key('mode'):
1147 if kwargs.has_key('mode'):
1131 mode = kwargs['mode']
1148 mode = kwargs['mode']
1132 else: mode = 'SA'
1149 else: mode = 'SA'
1133
1150
1134 #Borrar luego esto
1151 #Borrar luego esto
1135 if dataOut.groupList == None:
1152 if dataOut.groupList == None:
1136 dataOut.groupList = [(0,1),(0,2),(1,2)]
1153 dataOut.groupList = [(0,1),(0,2),(1,2)]
1137 groupList = dataOut.groupList
1154 groupList = dataOut.groupList
1138 C = 3e8
1155 C = 3e8
1139 freq = 50e6
1156 freq = 50e6
1140 lamb = C/freq
1157 lamb = C/freq
1141 k = 2*numpy.pi/lamb
1158 k = 2*numpy.pi/lamb
1142
1159
1143 timeList = dataOut.abscissaList
1160 timeList = dataOut.abscissaList
1144 heightList = dataOut.heightList
1161 heightList = dataOut.heightList
1145
1162
1146 if self.__isConfig == False:
1163 if self.__isConfig == False:
1147 dataOut.outputInterval = nMins*60
1164 dataOut.outputInterval = nMins*60
1148 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
1165 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
1149 #Get Initial LTC time
1166 #Get Initial LTC time
1150 initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
1167 initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
1151 minuteAux = initime.minute
1168 minuteAux = initime.minute
1152 minuteNew = int(numpy.floor(minuteAux/nMins)*nMins)
1169 minuteNew = int(numpy.floor(minuteAux/nMins)*nMins)
1153 self.__initime = (initime.replace(minute = minuteNew, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
1170 self.__initime = (initime.replace(minute = minuteNew, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
1154
1171
1155 self.__isConfig = True
1172 self.__isConfig = True
1156
1173
1157 if self.__buffer == None:
1174 if self.__buffer == None:
1158 self.__buffer = dataOut.data_param
1175 self.__buffer = dataOut.data_param
1159 self.__firstdata = copy.copy(dataOut)
1176 self.__firstdata = copy.copy(dataOut)
1160
1177
1161 else:
1178 else:
1162 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
1179 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
1163
1180
1164 self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
1181 self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
1165
1182
1166 if self.__dataReady:
1183 if self.__dataReady:
1167 dataOut.utctimeInit = self.__initime
1184 dataOut.utctimeInit = self.__initime
1168 self.__initime += dataOut.outputInterval #to erase time offset
1185 self.__initime += dataOut.outputInterval #to erase time offset
1169
1186
1170 metArray = self.__buffer
1187 metArray = self.__buffer
1171 if mode == 'SA':
1188 if mode == 'SA':
1172 dataOut.data_output = self.techniqueNSM_SA(rx_location=rx_location, groupList=groupList, azimuth=azimuth, dfactor=dfactor, k=k,metArray=metArray, heightList=heightList,timeList=timeList)
1189 dataOut.data_output = self.techniqueNSM_SA(rx_location=rx_location, groupList=groupList, azimuth=azimuth, dfactor=dfactor, k=k,metArray=metArray, heightList=heightList,timeList=timeList)
1173 elif mode == 'DBS':
1190 elif mode == 'DBS':
1174 dataOut.data_output = self.techniqueNSM_DBS(metArray=metArray,heightList=heightList,timeList=timeList)
1191 dataOut.data_output = self.techniqueNSM_DBS(metArray=metArray,heightList=heightList,timeList=timeList)
1175 dataOut.data_output = dataOut.data_output.T
1192 dataOut.data_output = dataOut.data_output.T
1176 dataOut.flagNoData = False
1193 dataOut.flagNoData = False
1177 self.__buffer = None
1194 self.__buffer = None
1178
1195
1179 return
1196 return
1180
1197
1181 class EWDriftsEstimation(Operation):
1198 class EWDriftsEstimation(Operation):
1182
1199
1183 def __init__(self):
1200 def __init__(self):
1184 Operation.__init__(self)
1201 Operation.__init__(self)
1185
1202
1186 def __correctValues(self, heiRang, phi, velRadial, SNR):
1203 def __correctValues(self, heiRang, phi, velRadial, SNR):
1187 listPhi = phi.tolist()
1204 listPhi = phi.tolist()
1188 maxid = listPhi.index(max(listPhi))
1205 maxid = listPhi.index(max(listPhi))
1189 minid = listPhi.index(min(listPhi))
1206 minid = listPhi.index(min(listPhi))
1190
1207
1191 rango = range(len(phi))
1208 rango = range(len(phi))
1192 # rango = numpy.delete(rango,maxid)
1209 # rango = numpy.delete(rango,maxid)
1193
1210
1194 heiRang1 = heiRang*math.cos(phi[maxid])
1211 heiRang1 = heiRang*math.cos(phi[maxid])
1195 heiRangAux = heiRang*math.cos(phi[minid])
1212 heiRangAux = heiRang*math.cos(phi[minid])
1196 indOut = (heiRang1 < heiRangAux[0]).nonzero()
1213 indOut = (heiRang1 < heiRangAux[0]).nonzero()
1197 heiRang1 = numpy.delete(heiRang1,indOut)
1214 heiRang1 = numpy.delete(heiRang1,indOut)
1198
1215
1199 velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
1216 velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
1200 SNR1 = numpy.zeros([len(phi),len(heiRang1)])
1217 SNR1 = numpy.zeros([len(phi),len(heiRang1)])
1201
1218
1202 for i in rango:
1219 for i in rango:
1203 x = heiRang*math.cos(phi[i])
1220 x = heiRang*math.cos(phi[i])
1204 y1 = velRadial[i,:]
1221 y1 = velRadial[i,:]
1205 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
1222 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
1206
1223
1207 x1 = heiRang1
1224 x1 = heiRang1
1208 y11 = f1(x1)
1225 y11 = f1(x1)
1209
1226
1210 y2 = SNR[i,:]
1227 y2 = SNR[i,:]
1211 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
1228 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
1212 y21 = f2(x1)
1229 y21 = f2(x1)
1213
1230
1214 velRadial1[i,:] = y11
1231 velRadial1[i,:] = y11
1215 SNR1[i,:] = y21
1232 SNR1[i,:] = y21
1216
1233
1217 return heiRang1, velRadial1, SNR1
1234 return heiRang1, velRadial1, SNR1
1218
1235
1219 def run(self, dataOut, zenith, zenithCorrection):
1236 def run(self, dataOut, zenith, zenithCorrection):
1220 heiRang = dataOut.heightList
1237 heiRang = dataOut.heightList
1221 velRadial = dataOut.data_param[:,3,:]
1238 velRadial = dataOut.data_param[:,3,:]
1222 SNR = dataOut.data_SNR
1239 SNR = dataOut.data_SNR
1223
1240
1224 zenith = numpy.array(zenith)
1241 zenith = numpy.array(zenith)
1225 zenith -= zenithCorrection
1242 zenith -= zenithCorrection
1226 zenith *= numpy.pi/180
1243 zenith *= numpy.pi/180
1227
1244
1228 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR)
1245 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR)
1229
1246
1230 alp = zenith[0]
1247 alp = zenith[0]
1231 bet = zenith[1]
1248 bet = zenith[1]
1232
1249
1233 w_w = velRadial1[0,:]
1250 w_w = velRadial1[0,:]
1234 w_e = velRadial1[1,:]
1251 w_e = velRadial1[1,:]
1235
1252
1236 w = (w_w*numpy.sin(bet) - w_e*numpy.sin(alp))/(numpy.cos(alp)*numpy.sin(bet) - numpy.cos(bet)*numpy.sin(alp))
1253 w = (w_w*numpy.sin(bet) - w_e*numpy.sin(alp))/(numpy.cos(alp)*numpy.sin(bet) - numpy.cos(bet)*numpy.sin(alp))
1237 u = (w_w*numpy.cos(bet) - w_e*numpy.cos(alp))/(numpy.sin(alp)*numpy.cos(bet) - numpy.sin(bet)*numpy.cos(alp))
1254 u = (w_w*numpy.cos(bet) - w_e*numpy.cos(alp))/(numpy.sin(alp)*numpy.cos(bet) - numpy.sin(bet)*numpy.cos(alp))
1238
1255
1239 winds = numpy.vstack((u,w))
1256 winds = numpy.vstack((u,w))
1240
1257
1241 dataOut.heightList = heiRang1
1258 dataOut.heightList = heiRang1
1242 dataOut.data_output = winds
1259 dataOut.data_output = winds
1243 dataOut.data_SNR = SNR1
1260 dataOut.data_SNR = SNR1
1244
1261
1245 dataOut.utctimeInit = dataOut.utctime
1262 dataOut.utctimeInit = dataOut.utctime
1246 dataOut.outputInterval = dataOut.timeInterval
1263 dataOut.outputInterval = dataOut.timeInterval
1247 return
1264 return
1248
1265
1249 #--------------- Non Specular Meteor ----------------
1266 #--------------- Non Specular Meteor ----------------
1250
1267
1251 class NonSpecularMeteorDetection(Operation):
1268 class NonSpecularMeteorDetection(Operation):
1252
1269
1253 def run(self, mode, SNRthresh=8, phaseDerThresh=0.5, cohThresh=0.8, allData = False):
1270 def run(self, mode, SNRthresh=8, phaseDerThresh=0.5, cohThresh=0.8, allData = False):
1254 data_acf = self.dataOut.data_pre[0]
1271 data_acf = self.dataOut.data_pre[0]
1255 data_ccf = self.dataOut.data_pre[1]
1272 data_ccf = self.dataOut.data_pre[1]
1256
1273
1257 lamb = self.dataOut.C/self.dataOut.frequency
1274 lamb = self.dataOut.C/self.dataOut.frequency
1258 tSamp = self.dataOut.ippSeconds*self.dataOut.nCohInt
1275 tSamp = self.dataOut.ippSeconds*self.dataOut.nCohInt
1259 paramInterval = self.dataOut.paramInterval
1276 paramInterval = self.dataOut.paramInterval
1260
1277
1261 nChannels = data_acf.shape[0]
1278 nChannels = data_acf.shape[0]
1262 nLags = data_acf.shape[1]
1279 nLags = data_acf.shape[1]
1263 nProfiles = data_acf.shape[2]
1280 nProfiles = data_acf.shape[2]
1264 nHeights = self.dataOut.nHeights
1281 nHeights = self.dataOut.nHeights
1265 nCohInt = self.dataOut.nCohInt
1282 nCohInt = self.dataOut.nCohInt
1266 sec = numpy.round(nProfiles/self.dataOut.paramInterval)
1283 sec = numpy.round(nProfiles/self.dataOut.paramInterval)
1267 heightList = self.dataOut.heightList
1284 heightList = self.dataOut.heightList
1268 ippSeconds = self.dataOut.ippSeconds*self.dataOut.nCohInt*self.dataOut.nAvg
1285 ippSeconds = self.dataOut.ippSeconds*self.dataOut.nCohInt*self.dataOut.nAvg
1269 utctime = self.dataOut.utctime
1286 utctime = self.dataOut.utctime
1270
1287
1271 self.dataOut.abscissaList = numpy.arange(0,paramInterval+ippSeconds,ippSeconds)
1288 self.dataOut.abscissaList = numpy.arange(0,paramInterval+ippSeconds,ippSeconds)
1272
1289
1273 #------------------------ SNR --------------------------------------
1290 #------------------------ SNR --------------------------------------
1274 power = data_acf[:,0,:,:].real
1291 power = data_acf[:,0,:,:].real
1275 noise = numpy.zeros(nChannels)
1292 noise = numpy.zeros(nChannels)
1276 SNR = numpy.zeros(power.shape)
1293 SNR = numpy.zeros(power.shape)
1277 for i in range(nChannels):
1294 for i in range(nChannels):
1278 noise[i] = hildebrand_sekhon(power[i,:], nCohInt)
1295 noise[i] = hildebrand_sekhon(power[i,:], nCohInt)
1279 SNR[i] = (power[i]-noise[i])/noise[i]
1296 SNR[i] = (power[i]-noise[i])/noise[i]
1280 SNRm = numpy.nanmean(SNR, axis = 0)
1297 SNRm = numpy.nanmean(SNR, axis = 0)
1281 SNRdB = 10*numpy.log10(SNR)
1298 SNRdB = 10*numpy.log10(SNR)
1282
1299
1283 if mode == 'SA':
1300 if mode == 'SA':
1284 nPairs = data_ccf.shape[0]
1301 nPairs = data_ccf.shape[0]
1285 #---------------------- Coherence and Phase --------------------------
1302 #---------------------- Coherence and Phase --------------------------
1286 phase = numpy.zeros(data_ccf[:,0,:,:].shape)
1303 phase = numpy.zeros(data_ccf[:,0,:,:].shape)
1287 # phase1 = numpy.copy(phase)
1304 # phase1 = numpy.copy(phase)
1288 coh1 = numpy.zeros(data_ccf[:,0,:,:].shape)
1305 coh1 = numpy.zeros(data_ccf[:,0,:,:].shape)
1289
1306
1290 for p in range(nPairs):
1307 for p in range(nPairs):
1291 ch0 = self.dataOut.groupList[p][0]
1308 ch0 = self.dataOut.groupList[p][0]
1292 ch1 = self.dataOut.groupList[p][1]
1309 ch1 = self.dataOut.groupList[p][1]
1293 ccf = data_ccf[p,0,:,:]/numpy.sqrt(data_acf[ch0,0,:,:]*data_acf[ch1,0,:,:])
1310 ccf = data_ccf[p,0,:,:]/numpy.sqrt(data_acf[ch0,0,:,:]*data_acf[ch1,0,:,:])
1294 phase[p,:,:] = ndimage.median_filter(numpy.angle(ccf), size = (5,1)) #median filter
1311 phase[p,:,:] = ndimage.median_filter(numpy.angle(ccf), size = (5,1)) #median filter
1295 # phase1[p,:,:] = numpy.angle(ccf) #median filter
1312 # phase1[p,:,:] = numpy.angle(ccf) #median filter
1296 coh1[p,:,:] = ndimage.median_filter(numpy.abs(ccf), 5) #median filter
1313 coh1[p,:,:] = ndimage.median_filter(numpy.abs(ccf), 5) #median filter
1297 # coh1[p,:,:] = numpy.abs(ccf) #median filter
1314 # coh1[p,:,:] = numpy.abs(ccf) #median filter
1298 coh = numpy.nanmax(coh1, axis = 0)
1315 coh = numpy.nanmax(coh1, axis = 0)
1299 # struc = numpy.ones((5,1))
1316 # struc = numpy.ones((5,1))
1300 # coh = ndimage.morphology.grey_dilation(coh, size=(10,1))
1317 # coh = ndimage.morphology.grey_dilation(coh, size=(10,1))
1301 #---------------------- Radial Velocity ----------------------------
1318 #---------------------- Radial Velocity ----------------------------
1302 phaseAux = numpy.mean(numpy.angle(data_acf[:,1,:,:]), axis = 0)
1319 phaseAux = numpy.mean(numpy.angle(data_acf[:,1,:,:]), axis = 0)
1303 velRad = phaseAux*lamb/(4*numpy.pi*tSamp)
1320 velRad = phaseAux*lamb/(4*numpy.pi*tSamp)
1304
1321
1305 if allData:
1322 if allData:
1306 boolMetFin = ~numpy.isnan(SNRm)
1323 boolMetFin = ~numpy.isnan(SNRm)
1307 # coh[:-1,:] = numpy.nanmean(numpy.abs(phase[:,1:,:] - phase[:,:-1,:]),axis=0)
1324 # coh[:-1,:] = numpy.nanmean(numpy.abs(phase[:,1:,:] - phase[:,:-1,:]),axis=0)
1308 else:
1325 else:
1309 #------------------------ Meteor mask ---------------------------------
1326 #------------------------ Meteor mask ---------------------------------
1310 # #SNR mask
1327 # #SNR mask
1311 # boolMet = (SNRdB>SNRthresh)#|(~numpy.isnan(SNRdB))
1328 # boolMet = (SNRdB>SNRthresh)#|(~numpy.isnan(SNRdB))
1312 #
1329 #
1313 # #Erase small objects
1330 # #Erase small objects
1314 # boolMet1 = self.__erase_small(boolMet, 2*sec, 5)
1331 # boolMet1 = self.__erase_small(boolMet, 2*sec, 5)
1315 #
1332 #
1316 # auxEEJ = numpy.sum(boolMet1,axis=0)
1333 # auxEEJ = numpy.sum(boolMet1,axis=0)
1317 # indOver = auxEEJ>nProfiles*0.8 #Use this later
1334 # indOver = auxEEJ>nProfiles*0.8 #Use this later
1318 # indEEJ = numpy.where(indOver)[0]
1335 # indEEJ = numpy.where(indOver)[0]
1319 # indNEEJ = numpy.where(~indOver)[0]
1336 # indNEEJ = numpy.where(~indOver)[0]
1320 #
1337 #
1321 # boolMetFin = boolMet1
1338 # boolMetFin = boolMet1
1322 #
1339 #
1323 # if indEEJ.size > 0:
1340 # if indEEJ.size > 0:
1324 # boolMet1[:,indEEJ] = False #Erase heights with EEJ
1341 # boolMet1[:,indEEJ] = False #Erase heights with EEJ
1325 #
1342 #
1326 # boolMet2 = coh > cohThresh
1343 # boolMet2 = coh > cohThresh
1327 # boolMet2 = self.__erase_small(boolMet2, 2*sec,5)
1344 # boolMet2 = self.__erase_small(boolMet2, 2*sec,5)
1328 #
1345 #
1329 # #Final Meteor mask
1346 # #Final Meteor mask
1330 # boolMetFin = boolMet1|boolMet2
1347 # boolMetFin = boolMet1|boolMet2
1331
1348
1332 #Coherence mask
1349 #Coherence mask
1333 boolMet1 = coh > 0.75
1350 boolMet1 = coh > 0.75
1334 struc = numpy.ones((30,1))
1351 struc = numpy.ones((30,1))
1335 boolMet1 = ndimage.morphology.binary_dilation(boolMet1, structure=struc)
1352 boolMet1 = ndimage.morphology.binary_dilation(boolMet1, structure=struc)
1336
1353
1337 #Derivative mask
1354 #Derivative mask
1338 derPhase = numpy.nanmean(numpy.abs(phase[:,1:,:] - phase[:,:-1,:]),axis=0)
1355 derPhase = numpy.nanmean(numpy.abs(phase[:,1:,:] - phase[:,:-1,:]),axis=0)
1339 boolMet2 = derPhase < 0.2
1356 boolMet2 = derPhase < 0.2
1340 # boolMet2 = ndimage.morphology.binary_opening(boolMet2)
1357 # boolMet2 = ndimage.morphology.binary_opening(boolMet2)
1341 # boolMet2 = ndimage.morphology.binary_closing(boolMet2, structure = numpy.ones((10,1)))
1358 # boolMet2 = ndimage.morphology.binary_closing(boolMet2, structure = numpy.ones((10,1)))
1342 boolMet2 = ndimage.median_filter(boolMet2,size=5)
1359 boolMet2 = ndimage.median_filter(boolMet2,size=5)
1343 boolMet2 = numpy.vstack((boolMet2,numpy.full((1,nHeights), True, dtype=bool)))
1360 boolMet2 = numpy.vstack((boolMet2,numpy.full((1,nHeights), True, dtype=bool)))
1344 # #Final mask
1361 # #Final mask
1345 # boolMetFin = boolMet2
1362 # boolMetFin = boolMet2
1346 boolMetFin = boolMet1&boolMet2
1363 boolMetFin = boolMet1&boolMet2
1347 # boolMetFin = ndimage.morphology.binary_dilation(boolMetFin)
1364 # boolMetFin = ndimage.morphology.binary_dilation(boolMetFin)
1348 #Creating data_param
1365 #Creating data_param
1349 coordMet = numpy.where(boolMetFin)
1366 coordMet = numpy.where(boolMetFin)
1350
1367
1351 tmet = coordMet[0]
1368 tmet = coordMet[0]
1352 hmet = coordMet[1]
1369 hmet = coordMet[1]
1353
1370
1354 data_param = numpy.zeros((tmet.size, 6 + nPairs))
1371 data_param = numpy.zeros((tmet.size, 6 + nPairs))
1355 data_param[:,0] = utctime
1372 data_param[:,0] = utctime
1356 data_param[:,1] = tmet
1373 data_param[:,1] = tmet
1357 data_param[:,2] = hmet
1374 data_param[:,2] = hmet
1358 data_param[:,3] = SNRm[tmet,hmet]
1375 data_param[:,3] = SNRm[tmet,hmet]
1359 data_param[:,4] = velRad[tmet,hmet]
1376 data_param[:,4] = velRad[tmet,hmet]
1360 data_param[:,5] = coh[tmet,hmet]
1377 data_param[:,5] = coh[tmet,hmet]
1361 data_param[:,6:] = phase[:,tmet,hmet].T
1378 data_param[:,6:] = phase[:,tmet,hmet].T
1362
1379
1363 elif mode == 'DBS':
1380 elif mode == 'DBS':
1364 self.dataOut.groupList = numpy.arange(nChannels)
1381 self.dataOut.groupList = numpy.arange(nChannels)
1365
1382
1366 #Radial Velocities
1383 #Radial Velocities
1367 # phase = numpy.angle(data_acf[:,1,:,:])
1384 # phase = numpy.angle(data_acf[:,1,:,:])
1368 phase = ndimage.median_filter(numpy.angle(data_acf[:,1,:,:]), size = (1,5,1))
1385 phase = ndimage.median_filter(numpy.angle(data_acf[:,1,:,:]), size = (1,5,1))
1369 velRad = phase*lamb/(4*numpy.pi*tSamp)
1386 velRad = phase*lamb/(4*numpy.pi*tSamp)
1370
1387
1371 #Spectral width
1388 #Spectral width
1372 acf1 = ndimage.median_filter(numpy.abs(data_acf[:,1,:,:]), size = (1,5,1))
1389 acf1 = ndimage.median_filter(numpy.abs(data_acf[:,1,:,:]), size = (1,5,1))
1373 acf2 = ndimage.median_filter(numpy.abs(data_acf[:,2,:,:]), size = (1,5,1))
1390 acf2 = ndimage.median_filter(numpy.abs(data_acf[:,2,:,:]), size = (1,5,1))
1374
1391
1375 spcWidth = (lamb/(2*numpy.sqrt(6)*numpy.pi*tSamp))*numpy.sqrt(numpy.log(acf1/acf2))
1392 spcWidth = (lamb/(2*numpy.sqrt(6)*numpy.pi*tSamp))*numpy.sqrt(numpy.log(acf1/acf2))
1376 # velRad = ndimage.median_filter(velRad, size = (1,5,1))
1393 # velRad = ndimage.median_filter(velRad, size = (1,5,1))
1377 if allData:
1394 if allData:
1378 boolMetFin = ~numpy.isnan(SNRdB)
1395 boolMetFin = ~numpy.isnan(SNRdB)
1379 else:
1396 else:
1380 #SNR
1397 #SNR
1381 boolMet1 = (SNRdB>SNRthresh) #SNR mask
1398 boolMet1 = (SNRdB>SNRthresh) #SNR mask
1382 boolMet1 = ndimage.median_filter(boolMet1, size=(1,5,5))
1399 boolMet1 = ndimage.median_filter(boolMet1, size=(1,5,5))
1383
1400
1384 #Radial velocity
1401 #Radial velocity
1385 boolMet2 = numpy.abs(velRad) < 30
1402 boolMet2 = numpy.abs(velRad) < 30
1386 boolMet2 = ndimage.median_filter(boolMet2, (1,5,5))
1403 boolMet2 = ndimage.median_filter(boolMet2, (1,5,5))
1387
1404
1388 #Spectral Width
1405 #Spectral Width
1389 boolMet3 = spcWidth < 30
1406 boolMet3 = spcWidth < 30
1390 boolMet3 = ndimage.median_filter(boolMet3, (1,5,5))
1407 boolMet3 = ndimage.median_filter(boolMet3, (1,5,5))
1391 # boolMetFin = self.__erase_small(boolMet1, 10,5)
1408 # boolMetFin = self.__erase_small(boolMet1, 10,5)
1392 boolMetFin = boolMet1&boolMet2&boolMet3
1409 boolMetFin = boolMet1&boolMet2&boolMet3
1393
1410
1394 #Creating data_param
1411 #Creating data_param
1395 coordMet = numpy.where(boolMetFin)
1412 coordMet = numpy.where(boolMetFin)
1396
1413
1397 cmet = coordMet[0]
1414 cmet = coordMet[0]
1398 tmet = coordMet[1]
1415 tmet = coordMet[1]
1399 hmet = coordMet[2]
1416 hmet = coordMet[2]
1400
1417
1401 data_param = numpy.zeros((tmet.size, 7))
1418 data_param = numpy.zeros((tmet.size, 7))
1402 data_param[:,0] = utctime
1419 data_param[:,0] = utctime
1403 data_param[:,1] = cmet
1420 data_param[:,1] = cmet
1404 data_param[:,2] = tmet
1421 data_param[:,2] = tmet
1405 data_param[:,3] = hmet
1422 data_param[:,3] = hmet
1406 data_param[:,4] = SNR[cmet,tmet,hmet].T
1423 data_param[:,4] = SNR[cmet,tmet,hmet].T
1407 data_param[:,5] = velRad[cmet,tmet,hmet].T
1424 data_param[:,5] = velRad[cmet,tmet,hmet].T
1408 data_param[:,6] = spcWidth[cmet,tmet,hmet].T
1425 data_param[:,6] = spcWidth[cmet,tmet,hmet].T
1409
1426
1410 # self.dataOut.data_param = data_int
1427 # self.dataOut.data_param = data_int
1411 if len(data_param) == 0:
1428 if len(data_param) == 0:
1412 self.dataOut.flagNoData = True
1429 self.dataOut.flagNoData = True
1413 else:
1430 else:
1414 self.dataOut.data_param = data_param
1431 self.dataOut.data_param = data_param
1415
1432
1416 def __erase_small(self, binArray, threshX, threshY):
1433 def __erase_small(self, binArray, threshX, threshY):
1417 labarray, numfeat = ndimage.measurements.label(binArray)
1434 labarray, numfeat = ndimage.measurements.label(binArray)
1418 binArray1 = numpy.copy(binArray)
1435 binArray1 = numpy.copy(binArray)
1419
1436
1420 for i in range(1,numfeat + 1):
1437 for i in range(1,numfeat + 1):
1421 auxBin = (labarray==i)
1438 auxBin = (labarray==i)
1422 auxSize = auxBin.sum()
1439 auxSize = auxBin.sum()
1423
1440
1424 x,y = numpy.where(auxBin)
1441 x,y = numpy.where(auxBin)
1425 widthX = x.max() - x.min()
1442 widthX = x.max() - x.min()
1426 widthY = y.max() - y.min()
1443 widthY = y.max() - y.min()
1427
1444
1428 #width X: 3 seg -> 12.5*3
1445 #width X: 3 seg -> 12.5*3
1429 #width Y:
1446 #width Y:
1430
1447
1431 if (auxSize < 50) or (widthX < threshX) or (widthY < threshY):
1448 if (auxSize < 50) or (widthX < threshX) or (widthY < threshY):
1432 binArray1[auxBin] = False
1449 binArray1[auxBin] = False
1433
1450
1434 return binArray1
1451 return binArray1
1435
1452
1436 #--------------- Specular Meteor ----------------
1453 #--------------- Specular Meteor ----------------
1437
1454
1438 class MeteorDetection(Operation):
1455 class SMDetection(Operation):
1439 '''
1456 '''
1440 Function DetectMeteors()
1457 Function DetectMeteors()
1441 Project developed with paper:
1458 Project developed with paper:
1442 HOLDSWORTH ET AL. 2004
1459 HOLDSWORTH ET AL. 2004
1443
1460
1444 Input:
1461 Input:
1445 self.dataOut.data_pre
1462 self.dataOut.data_pre
1446
1463
1447 centerReceiverIndex: From the channels, which is the center receiver
1464 centerReceiverIndex: From the channels, which is the center receiver
1448
1465
1449 hei_ref: Height reference for the Beacon signal extraction
1466 hei_ref: Height reference for the Beacon signal extraction
1450 tauindex:
1467 tauindex:
1451 predefinedPhaseShifts: Predefined phase offset for the voltge signals
1468 predefinedPhaseShifts: Predefined phase offset for the voltge signals
1452
1469
1453 cohDetection: Whether to user Coherent detection or not
1470 cohDetection: Whether to user Coherent detection or not
1454 cohDet_timeStep: Coherent Detection calculation time step
1471 cohDet_timeStep: Coherent Detection calculation time step
1455 cohDet_thresh: Coherent Detection phase threshold to correct phases
1472 cohDet_thresh: Coherent Detection phase threshold to correct phases
1456
1473
1457 noise_timeStep: Noise calculation time step
1474 noise_timeStep: Noise calculation time step
1458 noise_multiple: Noise multiple to define signal threshold
1475 noise_multiple: Noise multiple to define signal threshold
1459
1476
1460 multDet_timeLimit: Multiple Detection Removal time limit in seconds
1477 multDet_timeLimit: Multiple Detection Removal time limit in seconds
1461 multDet_rangeLimit: Multiple Detection Removal range limit in km
1478 multDet_rangeLimit: Multiple Detection Removal range limit in km
1462
1479
1463 phaseThresh: Maximum phase difference between receiver to be consider a meteor
1480 phaseThresh: Maximum phase difference between receiver to be consider a meteor
1464 SNRThresh: Minimum SNR threshold of the meteor signal to be consider a meteor
1481 SNRThresh: Minimum SNR threshold of the meteor signal to be consider a meteor
1465
1482
1466 hmin: Minimum Height of the meteor to use it in the further wind estimations
1483 hmin: Minimum Height of the meteor to use it in the further wind estimations
1467 hmax: Maximum Height of the meteor to use it in the further wind estimations
1484 hmax: Maximum Height of the meteor to use it in the further wind estimations
1468 azimuth: Azimuth angle correction
1485 azimuth: Azimuth angle correction
1469
1486
1470 Affected:
1487 Affected:
1471 self.dataOut.data_param
1488 self.dataOut.data_param
1472
1489
1473 Rejection Criteria (Errors):
1490 Rejection Criteria (Errors):
1474 0: No error; analysis OK
1491 0: No error; analysis OK
1475 1: SNR < SNR threshold
1492 1: SNR < SNR threshold
1476 2: angle of arrival (AOA) ambiguously determined
1493 2: angle of arrival (AOA) ambiguously determined
1477 3: AOA estimate not feasible
1494 3: AOA estimate not feasible
1478 4: Large difference in AOAs obtained from different antenna baselines
1495 4: Large difference in AOAs obtained from different antenna baselines
1479 5: echo at start or end of time series
1496 5: echo at start or end of time series
1480 6: echo less than 5 examples long; too short for analysis
1497 6: echo less than 5 examples long; too short for analysis
1481 7: echo rise exceeds 0.3s
1498 7: echo rise exceeds 0.3s
1482 8: echo decay time less than twice rise time
1499 8: echo decay time less than twice rise time
1483 9: large power level before echo
1500 9: large power level before echo
1484 10: large power level after echo
1501 10: large power level after echo
1485 11: poor fit to amplitude for estimation of decay time
1502 11: poor fit to amplitude for estimation of decay time
1486 12: poor fit to CCF phase variation for estimation of radial drift velocity
1503 12: poor fit to CCF phase variation for estimation of radial drift velocity
1487 13: height unresolvable echo: not valid height within 70 to 110 km
1504 13: height unresolvable echo: not valid height within 70 to 110 km
1488 14: height ambiguous echo: more then one possible height within 70 to 110 km
1505 14: height ambiguous echo: more then one possible height within 70 to 110 km
1489 15: radial drift velocity or projected horizontal velocity exceeds 200 m/s
1506 15: radial drift velocity or projected horizontal velocity exceeds 200 m/s
1490 16: oscilatory echo, indicating event most likely not an underdense echo
1507 16: oscilatory echo, indicating event most likely not an underdense echo
1491
1508
1492 17: phase difference in meteor Reestimation
1509 17: phase difference in meteor Reestimation
1493
1510
1494 Data Storage:
1511 Data Storage:
1495 Meteors for Wind Estimation (8):
1512 Meteors for Wind Estimation (8):
1496 Utc Time | Range Height
1513 Utc Time | Range Height
1497 Azimuth Zenith errorCosDir
1514 Azimuth Zenith errorCosDir
1498 VelRad errorVelRad
1515 VelRad errorVelRad
1499 Phase0 Phase1 Phase2 Phase3
1516 Phase0 Phase1 Phase2 Phase3
1500 TypeError
1517 TypeError
1501
1518
1502 '''
1519 '''
1503
1520
1504 def run(self, dataOut, hei_ref = None, tauindex = 0,
1521 def run(self, dataOut, hei_ref = None, tauindex = 0,
1505 phaseOffsets = None,
1522 phaseOffsets = None,
1506 cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25,
1523 cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25,
1507 noise_timeStep = 4, noise_multiple = 4,
1524 noise_timeStep = 4, noise_multiple = 4,
1508 multDet_timeLimit = 1, multDet_rangeLimit = 3,
1525 multDet_timeLimit = 1, multDet_rangeLimit = 3,
1509 phaseThresh = 20, SNRThresh = 5,
1526 phaseThresh = 20, SNRThresh = 5,
1510 hmin = 50, hmax=150, azimuth = 0,
1527 hmin = 50, hmax=150, azimuth = 0,
1511 channelPositions = None) :
1528 channelPositions = None) :
1512
1529
1513
1530
1514 #Getting Pairslist
1531 #Getting Pairslist
1515 if channelPositions == None:
1532 if channelPositions == None:
1516 # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
1533 # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
1517 channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella
1534 channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella
1518 meteorOps = MeteorOperations()
1535 meteorOps = SMOperations()
1519 pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
1536 pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
1520
1521 #Get Beacon signal
1522 newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
1523
1524 if hei_ref != None:
1525 newheis = numpy.where(self.dataOut.heightList>hei_ref)
1526
1527 heiRang = dataOut.getHeiRange()
1537 heiRang = dataOut.getHeiRange()
1538 #Get Beacon signal - No Beacon signal anymore
1539 # newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
1540 #
1541 # if hei_ref != None:
1542 # newheis = numpy.where(self.dataOut.heightList>hei_ref)
1543 #
1544
1528
1545
1529 #****************REMOVING HARDWARE PHASE DIFFERENCES***************
1546 #****************REMOVING HARDWARE PHASE DIFFERENCES***************
1530 # see if the user put in pre defined phase shifts
1547 # see if the user put in pre defined phase shifts
1531 voltsPShift = self.dataOut.data_pre.copy()
1548 voltsPShift = dataOut.data_pre.copy()
1532
1549
1533 # if predefinedPhaseShifts != None:
1550 # if predefinedPhaseShifts != None:
1534 # hardwarePhaseShifts = numpy.array(predefinedPhaseShifts)*numpy.pi/180
1551 # hardwarePhaseShifts = numpy.array(predefinedPhaseShifts)*numpy.pi/180
1535 #
1552 #
1536 # # elif beaconPhaseShifts:
1553 # # elif beaconPhaseShifts:
1537 # # #get hardware phase shifts using beacon signal
1554 # # #get hardware phase shifts using beacon signal
1538 # # hardwarePhaseShifts = self.__getHardwarePhaseDiff(self.dataOut.data_pre, pairslist, newheis, 10)
1555 # # hardwarePhaseShifts = self.__getHardwarePhaseDiff(self.dataOut.data_pre, pairslist, newheis, 10)
1539 # # hardwarePhaseShifts = numpy.insert(hardwarePhaseShifts,centerReceiverIndex,0)
1556 # # hardwarePhaseShifts = numpy.insert(hardwarePhaseShifts,centerReceiverIndex,0)
1540 #
1557 #
1541 # else:
1558 # else:
1542 # hardwarePhaseShifts = numpy.zeros(5)
1559 # hardwarePhaseShifts = numpy.zeros(5)
1543 #
1560 #
1544 # voltsPShift = numpy.zeros((self.dataOut.data_pre.shape[0],self.dataOut.data_pre.shape[1],self.dataOut.data_pre.shape[2]), dtype = 'complex')
1561 # voltsPShift = numpy.zeros((self.dataOut.data_pre.shape[0],self.dataOut.data_pre.shape[1],self.dataOut.data_pre.shape[2]), dtype = 'complex')
1545 # for i in range(self.dataOut.data_pre.shape[0]):
1562 # for i in range(self.dataOut.data_pre.shape[0]):
1546 # voltsPShift[i,:,:] = self.__shiftPhase(self.dataOut.data_pre[i,:,:], hardwarePhaseShifts[i])
1563 # voltsPShift[i,:,:] = self.__shiftPhase(self.dataOut.data_pre[i,:,:], hardwarePhaseShifts[i])
1547
1564
1548 #******************END OF REMOVING HARDWARE PHASE DIFFERENCES*********
1565 #******************END OF REMOVING HARDWARE PHASE DIFFERENCES*********
1549
1566
1550 #Remove DC
1567 #Remove DC
1551 voltsDC = numpy.mean(voltsPShift,1)
1568 voltsDC = numpy.mean(voltsPShift,1)
1552 voltsDC = numpy.mean(voltsDC,1)
1569 voltsDC = numpy.mean(voltsDC,1)
1553 for i in range(voltsDC.shape[0]):
1570 for i in range(voltsDC.shape[0]):
1554 voltsPShift[i] = voltsPShift[i] - voltsDC[i]
1571 voltsPShift[i] = voltsPShift[i] - voltsDC[i]
1555
1572
1556 #Don't considerate last heights, theyre used to calculate Hardware Phase Shift
1573 #Don't considerate last heights, theyre used to calculate Hardware Phase Shift
1557 voltsPShift = voltsPShift[:,:,:newheis[0][0]]
1574 # voltsPShift = voltsPShift[:,:,:newheis[0][0]]
1558
1575
1559 #************ FIND POWER OF DATA W/COH OR NON COH DETECTION (3.4) **********
1576 #************ FIND POWER OF DATA W/COH OR NON COH DETECTION (3.4) **********
1560 #Coherent Detection
1577 #Coherent Detection
1561 if cohDetection:
1578 if cohDetection:
1562 #use coherent detection to get the net power
1579 #use coherent detection to get the net power
1563 cohDet_thresh = cohDet_thresh*numpy.pi/180
1580 cohDet_thresh = cohDet_thresh*numpy.pi/180
1564 voltsPShift = self.__coherentDetection(voltsPShift, cohDet_timeStep, self.dataOut.timeInterval, pairslist0, cohDet_thresh)
1581 voltsPShift = self.__coherentDetection(voltsPShift, cohDet_timeStep, dataOut.timeInterval, pairslist0, cohDet_thresh)
1565
1582
1566 #Non-coherent detection!
1583 #Non-coherent detection!
1567 powerNet = numpy.nansum(numpy.abs(voltsPShift[:,:,:])**2,0)
1584 powerNet = numpy.nansum(numpy.abs(voltsPShift[:,:,:])**2,0)
1568 #********** END OF COH/NON-COH POWER CALCULATION**********************
1585 #********** END OF COH/NON-COH POWER CALCULATION**********************
1569
1586
1570 #********** FIND THE NOISE LEVEL AND POSSIBLE METEORS ****************
1587 #********** FIND THE NOISE LEVEL AND POSSIBLE METEORS ****************
1571 #Get noise
1588 #Get noise
1572 noise, noise1 = self.__getNoise(powerNet, noise_timeStep, self.dataOut.timeInterval)
1589 noise, noise1 = self.__getNoise(powerNet, noise_timeStep, dataOut.timeInterval)
1573 # noise = self.getNoise1(powerNet, noise_timeStep, self.dataOut.timeInterval)
1590 # noise = self.getNoise1(powerNet, noise_timeStep, self.dataOut.timeInterval)
1574 #Get signal threshold
1591 #Get signal threshold
1575 signalThresh = noise_multiple*noise
1592 signalThresh = noise_multiple*noise
1576 #Meteor echoes detection
1593 #Meteor echoes detection
1577 listMeteors = self.__findMeteors(powerNet, signalThresh)
1594 listMeteors = self.__findMeteors(powerNet, signalThresh)
1578 #******* END OF NOISE LEVEL AND POSSIBLE METEORS CACULATION **********
1595 #******* END OF NOISE LEVEL AND POSSIBLE METEORS CACULATION **********
1579
1596
1580 #************** REMOVE MULTIPLE DETECTIONS (3.5) ***************************
1597 #************** REMOVE MULTIPLE DETECTIONS (3.5) ***************************
1581 #Parameters
1598 #Parameters
1582 heiRange = self.dataOut.getHeiRange()
1599 heiRange = dataOut.getHeiRange()
1583 rangeInterval = heiRange[1] - heiRange[0]
1600 rangeInterval = heiRange[1] - heiRange[0]
1584 rangeLimit = multDet_rangeLimit/rangeInterval
1601 rangeLimit = multDet_rangeLimit/rangeInterval
1585 timeLimit = multDet_timeLimit/self.dataOut.timeInterval
1602 timeLimit = multDet_timeLimit/dataOut.timeInterval
1586 #Multiple detection removals
1603 #Multiple detection removals
1587 listMeteors1 = self.__removeMultipleDetections(listMeteors, rangeLimit, timeLimit)
1604 listMeteors1 = self.__removeMultipleDetections(listMeteors, rangeLimit, timeLimit)
1588 #************ END OF REMOVE MULTIPLE DETECTIONS **********************
1605 #************ END OF REMOVE MULTIPLE DETECTIONS **********************
1589
1606
1590 #********************* METEOR REESTIMATION (3.7, 3.8, 3.9, 3.10) ********************
1607 #********************* METEOR REESTIMATION (3.7, 3.8, 3.9, 3.10) ********************
1591 #Parameters
1608 #Parameters
1592 phaseThresh = phaseThresh*numpy.pi/180
1609 phaseThresh = phaseThresh*numpy.pi/180
1593 thresh = [phaseThresh, noise_multiple, SNRThresh]
1610 thresh = [phaseThresh, noise_multiple, SNRThresh]
1594 #Meteor reestimation (Errors N 1, 6, 12, 17)
1611 #Meteor reestimation (Errors N 1, 6, 12, 17)
1595 listMeteors2, listMeteorsPower, listMeteorsVolts = self.__meteorReestimation(listMeteors1, voltsPShift, pairslist0, thresh, noise, self.dataOut.timeInterval, self.dataOut.frequency)
1612 listMeteors2, listMeteorsPower, listMeteorsVolts = self.__meteorReestimation(listMeteors1, voltsPShift, pairslist0, thresh, noise, dataOut.timeInterval, dataOut.frequency)
1596 # listMeteors2, listMeteorsPower, listMeteorsVolts = self.meteorReestimation3(listMeteors2, listMeteorsPower, listMeteorsVolts, voltsPShift, pairslist, thresh, noise)
1613 # listMeteors2, listMeteorsPower, listMeteorsVolts = self.meteorReestimation3(listMeteors2, listMeteorsPower, listMeteorsVolts, voltsPShift, pairslist, thresh, noise)
1597 #Estimation of decay times (Errors N 7, 8, 11)
1614 #Estimation of decay times (Errors N 7, 8, 11)
1598 listMeteors3 = self.__estimateDecayTime(listMeteors2, listMeteorsPower, self.dataOut.timeInterval, self.dataOut.frequency)
1615 listMeteors3 = self.__estimateDecayTime(listMeteors2, listMeteorsPower, dataOut.timeInterval, dataOut.frequency)
1599 #******************* END OF METEOR REESTIMATION *******************
1616 #******************* END OF METEOR REESTIMATION *******************
1600
1617
1601 #********************* METEOR PARAMETERS CALCULATION (3.11, 3.12, 3.13) **************************
1618 #********************* METEOR PARAMETERS CALCULATION (3.11, 3.12, 3.13) **************************
1602 #Calculating Radial Velocity (Error N 15)
1619 #Calculating Radial Velocity (Error N 15)
1603 radialStdThresh = 10
1620 radialStdThresh = 10
1604 listMeteors4 = self.__getRadialVelocity(listMeteors3, listMeteorsVolts, radialStdThresh, pairslist0, self.dataOut.timeInterval)
1621 listMeteors4 = self.__getRadialVelocity(listMeteors3, listMeteorsVolts, radialStdThresh, pairslist0, dataOut.timeInterval)
1605
1622
1606 if len(listMeteors4) > 0:
1623 if len(listMeteors4) > 0:
1607 #Setting New Array
1624 #Setting New Array
1608 date = self.dataOut.utctime
1625 date = dataOut.utctime
1609 arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang)
1626 arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang)
1610
1627
1611 #Correcting phase offset
1628 #Correcting phase offset
1612 if phaseOffsets != None:
1629 if phaseOffsets != None:
1613 phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180
1630 phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180
1614 arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets)
1631 arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets)
1615
1632
1616 #Second Pairslist
1633 #Second Pairslist
1617 pairsList = []
1634 pairsList = []
1618 pairx = (0,1)
1635 pairx = (0,1)
1619 pairy = (2,3)
1636 pairy = (2,3)
1620 pairsList.append(pairx)
1637 pairsList.append(pairx)
1621 pairsList.append(pairy)
1638 pairsList.append(pairy)
1622
1639
1623 jph = numpy.array([0,0,0,0])
1640 jph = numpy.array([0,0,0,0])
1624 h = (hmin,hmax)
1641 h = (hmin,hmax)
1625 arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, jph)
1642 arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, jph)
1626
1643
1627 # #Calculate AOA (Error N 3, 4)
1644 # #Calculate AOA (Error N 3, 4)
1628 # #JONES ET AL. 1998
1645 # #JONES ET AL. 1998
1629 # error = arrayParameters[:,-1]
1646 # error = arrayParameters[:,-1]
1630 # AOAthresh = numpy.pi/8
1647 # AOAthresh = numpy.pi/8
1631 # phases = -arrayParameters[:,9:13]
1648 # phases = -arrayParameters[:,9:13]
1632 # arrayParameters[:,4:7], arrayParameters[:,-1] = meteorOps.getAOA(phases, pairsList, error, AOAthresh, azimuth)
1649 # arrayParameters[:,4:7], arrayParameters[:,-1] = meteorOps.getAOA(phases, pairsList, error, AOAthresh, azimuth)
1633 #
1650 #
1634 # #Calculate Heights (Error N 13 and 14)
1651 # #Calculate Heights (Error N 13 and 14)
1635 # error = arrayParameters[:,-1]
1652 # error = arrayParameters[:,-1]
1636 # Ranges = arrayParameters[:,2]
1653 # Ranges = arrayParameters[:,2]
1637 # zenith = arrayParameters[:,5]
1654 # zenith = arrayParameters[:,5]
1638 # arrayParameters[:,3], arrayParameters[:,-1] = meteorOps.getHeights(Ranges, zenith, error, hmin, hmax)
1655 # arrayParameters[:,3], arrayParameters[:,-1] = meteorOps.getHeights(Ranges, zenith, error, hmin, hmax)
1639 # error = arrayParameters[:,-1]
1656 # error = arrayParameters[:,-1]
1640 #********************* END OF PARAMETERS CALCULATION **************************
1657 #********************* END OF PARAMETERS CALCULATION **************************
1641
1658
1642 #***************************+ PASS DATA TO NEXT STEP **********************
1659 #***************************+ PASS DATA TO NEXT STEP **********************
1643 # arrayFinal = arrayParameters.reshape((1,arrayParameters.shape[0],arrayParameters.shape[1]))
1660 # arrayFinal = arrayParameters.reshape((1,arrayParameters.shape[0],arrayParameters.shape[1]))
1644 self.dataOut.data_param = arrayParameters
1661 dataOut.data_param = arrayParameters
1645
1662
1646 if arrayParameters == None:
1663 if arrayParameters == None:
1647 self.dataOut.flagNoData = True
1664 dataOut.flagNoData = True
1648 else:
1665 else:
1649 self.dataOut.flagNoData = True
1666 dataOut.flagNoData = True
1650
1667
1651 return
1668 return
1652
1669
1653 def __getHardwarePhaseDiff(self, voltage0, pairslist, newheis, n):
1670 def __getHardwarePhaseDiff(self, voltage0, pairslist, newheis, n):
1654
1671
1655 minIndex = min(newheis[0])
1672 minIndex = min(newheis[0])
1656 maxIndex = max(newheis[0])
1673 maxIndex = max(newheis[0])
1657
1674
1658 voltage = voltage0[:,:,minIndex:maxIndex+1]
1675 voltage = voltage0[:,:,minIndex:maxIndex+1]
1659 nLength = voltage.shape[1]/n
1676 nLength = voltage.shape[1]/n
1660 nMin = 0
1677 nMin = 0
1661 nMax = 0
1678 nMax = 0
1662 phaseOffset = numpy.zeros((len(pairslist),n))
1679 phaseOffset = numpy.zeros((len(pairslist),n))
1663
1680
1664 for i in range(n):
1681 for i in range(n):
1665 nMax += nLength
1682 nMax += nLength
1666 phaseCCF = -numpy.angle(self.__calculateCCF(voltage[:,nMin:nMax,:], pairslist, [0]))
1683 phaseCCF = -numpy.angle(self.__calculateCCF(voltage[:,nMin:nMax,:], pairslist, [0]))
1667 phaseCCF = numpy.mean(phaseCCF, axis = 2)
1684 phaseCCF = numpy.mean(phaseCCF, axis = 2)
1668 phaseOffset[:,i] = phaseCCF.transpose()
1685 phaseOffset[:,i] = phaseCCF.transpose()
1669 nMin = nMax
1686 nMin = nMax
1670 # phaseDiff, phaseArrival = self.estimatePhaseDifference(voltage, pairslist)
1687 # phaseDiff, phaseArrival = self.estimatePhaseDifference(voltage, pairslist)
1671
1688
1672 #Remove Outliers
1689 #Remove Outliers
1673 factor = 2
1690 factor = 2
1674 wt = phaseOffset - signal.medfilt(phaseOffset,(1,5))
1691 wt = phaseOffset - signal.medfilt(phaseOffset,(1,5))
1675 dw = numpy.std(wt,axis = 1)
1692 dw = numpy.std(wt,axis = 1)
1676 dw = dw.reshape((dw.size,1))
1693 dw = dw.reshape((dw.size,1))
1677 ind = numpy.where(numpy.logical_or(wt>dw*factor,wt<-dw*factor))
1694 ind = numpy.where(numpy.logical_or(wt>dw*factor,wt<-dw*factor))
1678 phaseOffset[ind] = numpy.nan
1695 phaseOffset[ind] = numpy.nan
1679 phaseOffset = stats.nanmean(phaseOffset, axis=1)
1696 phaseOffset = stats.nanmean(phaseOffset, axis=1)
1680
1697
1681 return phaseOffset
1698 return phaseOffset
1682
1699
1683 def __shiftPhase(self, data, phaseShift):
1700 def __shiftPhase(self, data, phaseShift):
1684 #this will shift the phase of a complex number
1701 #this will shift the phase of a complex number
1685 dataShifted = numpy.abs(data) * numpy.exp((numpy.angle(data)+phaseShift)*1j)
1702 dataShifted = numpy.abs(data) * numpy.exp((numpy.angle(data)+phaseShift)*1j)
1686 return dataShifted
1703 return dataShifted
1687
1704
1688 def __estimatePhaseDifference(self, array, pairslist):
1705 def __estimatePhaseDifference(self, array, pairslist):
1689 nChannel = array.shape[0]
1706 nChannel = array.shape[0]
1690 nHeights = array.shape[2]
1707 nHeights = array.shape[2]
1691 numPairs = len(pairslist)
1708 numPairs = len(pairslist)
1692 # phaseCCF = numpy.zeros((nChannel, 5, nHeights))
1709 # phaseCCF = numpy.zeros((nChannel, 5, nHeights))
1693 phaseCCF = numpy.angle(self.__calculateCCF(array, pairslist, [-2,-1,0,1,2]))
1710 phaseCCF = numpy.angle(self.__calculateCCF(array, pairslist, [-2,-1,0,1,2]))
1694
1711
1695 #Correct phases
1712 #Correct phases
1696 derPhaseCCF = phaseCCF[:,1:,:] - phaseCCF[:,0:-1,:]
1713 derPhaseCCF = phaseCCF[:,1:,:] - phaseCCF[:,0:-1,:]
1697 indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
1714 indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
1698
1715
1699 if indDer[0].shape[0] > 0:
1716 if indDer[0].shape[0] > 0:
1700 for i in range(indDer[0].shape[0]):
1717 for i in range(indDer[0].shape[0]):
1701 signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i],indDer[2][i]])
1718 signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i],indDer[2][i]])
1702 phaseCCF[indDer[0][i],indDer[1][i]+1:,:] += signo*2*numpy.pi
1719 phaseCCF[indDer[0][i],indDer[1][i]+1:,:] += signo*2*numpy.pi
1703
1720
1704 # for j in range(numSides):
1721 # for j in range(numSides):
1705 # phaseCCFAux = self.calculateCCF(arrayCenter, arraySides[j,:,:], [-2,1,0,1,2])
1722 # phaseCCFAux = self.calculateCCF(arrayCenter, arraySides[j,:,:], [-2,1,0,1,2])
1706 # phaseCCF[j,:,:] = numpy.angle(phaseCCFAux)
1723 # phaseCCF[j,:,:] = numpy.angle(phaseCCFAux)
1707 #
1724 #
1708 #Linear
1725 #Linear
1709 phaseInt = numpy.zeros((numPairs,1))
1726 phaseInt = numpy.zeros((numPairs,1))
1710 angAllCCF = phaseCCF[:,[0,1,3,4],0]
1727 angAllCCF = phaseCCF[:,[0,1,3,4],0]
1711 for j in range(numPairs):
1728 for j in range(numPairs):
1712 fit = stats.linregress([-2,-1,1,2],angAllCCF[j,:])
1729 fit = stats.linregress([-2,-1,1,2],angAllCCF[j,:])
1713 phaseInt[j] = fit[1]
1730 phaseInt[j] = fit[1]
1714 #Phase Differences
1731 #Phase Differences
1715 phaseDiff = phaseInt - phaseCCF[:,2,:]
1732 phaseDiff = phaseInt - phaseCCF[:,2,:]
1716 phaseArrival = phaseInt.reshape(phaseInt.size)
1733 phaseArrival = phaseInt.reshape(phaseInt.size)
1717
1734
1718 #Dealias
1735 #Dealias
1719 phaseArrival = numpy.angle(numpy.exp(1j*phaseArrival))
1736 phaseArrival = numpy.angle(numpy.exp(1j*phaseArrival))
1720 # indAlias = numpy.where(phaseArrival > numpy.pi)
1737 # indAlias = numpy.where(phaseArrival > numpy.pi)
1721 # phaseArrival[indAlias] -= 2*numpy.pi
1738 # phaseArrival[indAlias] -= 2*numpy.pi
1722 # indAlias = numpy.where(phaseArrival < -numpy.pi)
1739 # indAlias = numpy.where(phaseArrival < -numpy.pi)
1723 # phaseArrival[indAlias] += 2*numpy.pi
1740 # phaseArrival[indAlias] += 2*numpy.pi
1724
1741
1725 return phaseDiff, phaseArrival
1742 return phaseDiff, phaseArrival
1726
1743
1727 def __coherentDetection(self, volts, timeSegment, timeInterval, pairslist, thresh):
1744 def __coherentDetection(self, volts, timeSegment, timeInterval, pairslist, thresh):
1728 #this function will run the coherent detection used in Holdworth et al. 2004 and return the net power
1745 #this function will run the coherent detection used in Holdworth et al. 2004 and return the net power
1729 #find the phase shifts of each channel over 1 second intervals
1746 #find the phase shifts of each channel over 1 second intervals
1730 #only look at ranges below the beacon signal
1747 #only look at ranges below the beacon signal
1731 numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
1748 numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
1732 numBlocks = int(volts.shape[1]/numProfPerBlock)
1749 numBlocks = int(volts.shape[1]/numProfPerBlock)
1733 numHeights = volts.shape[2]
1750 numHeights = volts.shape[2]
1734 nChannel = volts.shape[0]
1751 nChannel = volts.shape[0]
1735 voltsCohDet = volts.copy()
1752 voltsCohDet = volts.copy()
1736
1753
1737 pairsarray = numpy.array(pairslist)
1754 pairsarray = numpy.array(pairslist)
1738 indSides = pairsarray[:,1]
1755 indSides = pairsarray[:,1]
1739 # indSides = numpy.array(range(nChannel))
1756 # indSides = numpy.array(range(nChannel))
1740 # indSides = numpy.delete(indSides, indCenter)
1757 # indSides = numpy.delete(indSides, indCenter)
1741 #
1758 #
1742 # listCenter = numpy.array_split(volts[indCenter,:,:], numBlocks, 0)
1759 # listCenter = numpy.array_split(volts[indCenter,:,:], numBlocks, 0)
1743 listBlocks = numpy.array_split(volts, numBlocks, 1)
1760 listBlocks = numpy.array_split(volts, numBlocks, 1)
1744
1761
1745 startInd = 0
1762 startInd = 0
1746 endInd = 0
1763 endInd = 0
1747
1764
1748 for i in range(numBlocks):
1765 for i in range(numBlocks):
1749 startInd = endInd
1766 startInd = endInd
1750 endInd = endInd + listBlocks[i].shape[1]
1767 endInd = endInd + listBlocks[i].shape[1]
1751
1768
1752 arrayBlock = listBlocks[i]
1769 arrayBlock = listBlocks[i]
1753 # arrayBlockCenter = listCenter[i]
1770 # arrayBlockCenter = listCenter[i]
1754
1771
1755 #Estimate the Phase Difference
1772 #Estimate the Phase Difference
1756 phaseDiff, aux = self.__estimatePhaseDifference(arrayBlock, pairslist)
1773 phaseDiff, aux = self.__estimatePhaseDifference(arrayBlock, pairslist)
1757 #Phase Difference RMS
1774 #Phase Difference RMS
1758 arrayPhaseRMS = numpy.abs(phaseDiff)
1775 arrayPhaseRMS = numpy.abs(phaseDiff)
1759 phaseRMSaux = numpy.sum(arrayPhaseRMS < thresh,0)
1776 phaseRMSaux = numpy.sum(arrayPhaseRMS < thresh,0)
1760 indPhase = numpy.where(phaseRMSaux==4)
1777 indPhase = numpy.where(phaseRMSaux==4)
1761 #Shifting
1778 #Shifting
1762 if indPhase[0].shape[0] > 0:
1779 if indPhase[0].shape[0] > 0:
1763 for j in range(indSides.size):
1780 for j in range(indSides.size):
1764 arrayBlock[indSides[j],:,indPhase] = self.__shiftPhase(arrayBlock[indSides[j],:,indPhase], phaseDiff[j,indPhase].transpose())
1781 arrayBlock[indSides[j],:,indPhase] = self.__shiftPhase(arrayBlock[indSides[j],:,indPhase], phaseDiff[j,indPhase].transpose())
1765 voltsCohDet[:,startInd:endInd,:] = arrayBlock
1782 voltsCohDet[:,startInd:endInd,:] = arrayBlock
1766
1783
1767 return voltsCohDet
1784 return voltsCohDet
1768
1785
1769 def __calculateCCF(self, volts, pairslist ,laglist):
1786 def __calculateCCF(self, volts, pairslist ,laglist):
1770
1787
1771 nHeights = volts.shape[2]
1788 nHeights = volts.shape[2]
1772 nPoints = volts.shape[1]
1789 nPoints = volts.shape[1]
1773 voltsCCF = numpy.zeros((len(pairslist), len(laglist), nHeights),dtype = 'complex')
1790 voltsCCF = numpy.zeros((len(pairslist), len(laglist), nHeights),dtype = 'complex')
1774
1791
1775 for i in range(len(pairslist)):
1792 for i in range(len(pairslist)):
1776 volts1 = volts[pairslist[i][0]]
1793 volts1 = volts[pairslist[i][0]]
1777 volts2 = volts[pairslist[i][1]]
1794 volts2 = volts[pairslist[i][1]]
1778
1795
1779 for t in range(len(laglist)):
1796 for t in range(len(laglist)):
1780 idxT = laglist[t]
1797 idxT = laglist[t]
1781 if idxT >= 0:
1798 if idxT >= 0:
1782 vStacked = numpy.vstack((volts2[idxT:,:],
1799 vStacked = numpy.vstack((volts2[idxT:,:],
1783 numpy.zeros((idxT, nHeights),dtype='complex')))
1800 numpy.zeros((idxT, nHeights),dtype='complex')))
1784 else:
1801 else:
1785 vStacked = numpy.vstack((numpy.zeros((-idxT, nHeights),dtype='complex'),
1802 vStacked = numpy.vstack((numpy.zeros((-idxT, nHeights),dtype='complex'),
1786 volts2[:(nPoints + idxT),:]))
1803 volts2[:(nPoints + idxT),:]))
1787 voltsCCF[i,t,:] = numpy.sum((numpy.conjugate(volts1)*vStacked),axis=0)
1804 voltsCCF[i,t,:] = numpy.sum((numpy.conjugate(volts1)*vStacked),axis=0)
1788
1805
1789 vStacked = None
1806 vStacked = None
1790 return voltsCCF
1807 return voltsCCF
1791
1808
1792 def __getNoise(self, power, timeSegment, timeInterval):
1809 def __getNoise(self, power, timeSegment, timeInterval):
1793 numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
1810 numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
1794 numBlocks = int(power.shape[0]/numProfPerBlock)
1811 numBlocks = int(power.shape[0]/numProfPerBlock)
1795 numHeights = power.shape[1]
1812 numHeights = power.shape[1]
1796
1813
1797 listPower = numpy.array_split(power, numBlocks, 0)
1814 listPower = numpy.array_split(power, numBlocks, 0)
1798 noise = numpy.zeros((power.shape[0], power.shape[1]))
1815 noise = numpy.zeros((power.shape[0], power.shape[1]))
1799 noise1 = numpy.zeros((power.shape[0], power.shape[1]))
1816 noise1 = numpy.zeros((power.shape[0], power.shape[1]))
1800
1817
1801 startInd = 0
1818 startInd = 0
1802 endInd = 0
1819 endInd = 0
1803
1820
1804 for i in range(numBlocks): #split por canal
1821 for i in range(numBlocks): #split por canal
1805 startInd = endInd
1822 startInd = endInd
1806 endInd = endInd + listPower[i].shape[0]
1823 endInd = endInd + listPower[i].shape[0]
1807
1824
1808 arrayBlock = listPower[i]
1825 arrayBlock = listPower[i]
1809 noiseAux = numpy.mean(arrayBlock, 0)
1826 noiseAux = numpy.mean(arrayBlock, 0)
1810 # noiseAux = numpy.median(noiseAux)
1827 # noiseAux = numpy.median(noiseAux)
1811 # noiseAux = numpy.mean(arrayBlock)
1828 # noiseAux = numpy.mean(arrayBlock)
1812 noise[startInd:endInd,:] = noise[startInd:endInd,:] + noiseAux
1829 noise[startInd:endInd,:] = noise[startInd:endInd,:] + noiseAux
1813
1830
1814 noiseAux1 = numpy.mean(arrayBlock)
1831 noiseAux1 = numpy.mean(arrayBlock)
1815 noise1[startInd:endInd,:] = noise1[startInd:endInd,:] + noiseAux1
1832 noise1[startInd:endInd,:] = noise1[startInd:endInd,:] + noiseAux1
1816
1833
1817 return noise, noise1
1834 return noise, noise1
1818
1835
1819 def __findMeteors(self, power, thresh):
1836 def __findMeteors(self, power, thresh):
1820 nProf = power.shape[0]
1837 nProf = power.shape[0]
1821 nHeights = power.shape[1]
1838 nHeights = power.shape[1]
1822 listMeteors = []
1839 listMeteors = []
1823
1840
1824 for i in range(nHeights):
1841 for i in range(nHeights):
1825 powerAux = power[:,i]
1842 powerAux = power[:,i]
1826 threshAux = thresh[:,i]
1843 threshAux = thresh[:,i]
1827
1844
1828 indUPthresh = numpy.where(powerAux > threshAux)[0]
1845 indUPthresh = numpy.where(powerAux > threshAux)[0]
1829 indDNthresh = numpy.where(powerAux <= threshAux)[0]
1846 indDNthresh = numpy.where(powerAux <= threshAux)[0]
1830
1847
1831 j = 0
1848 j = 0
1832
1849
1833 while (j < indUPthresh.size - 2):
1850 while (j < indUPthresh.size - 2):
1834 if (indUPthresh[j + 2] == indUPthresh[j] + 2):
1851 if (indUPthresh[j + 2] == indUPthresh[j] + 2):
1835 indDNAux = numpy.where(indDNthresh > indUPthresh[j])
1852 indDNAux = numpy.where(indDNthresh > indUPthresh[j])
1836 indDNthresh = indDNthresh[indDNAux]
1853 indDNthresh = indDNthresh[indDNAux]
1837
1854
1838 if (indDNthresh.size > 0):
1855 if (indDNthresh.size > 0):
1839 indEnd = indDNthresh[0] - 1
1856 indEnd = indDNthresh[0] - 1
1840 indInit = indUPthresh[j]
1857 indInit = indUPthresh[j]
1841
1858
1842 meteor = powerAux[indInit:indEnd + 1]
1859 meteor = powerAux[indInit:indEnd + 1]
1843 indPeak = meteor.argmax() + indInit
1860 indPeak = meteor.argmax() + indInit
1844 FLA = sum(numpy.conj(meteor)*numpy.hstack((meteor[1:],0)))
1861 FLA = sum(numpy.conj(meteor)*numpy.hstack((meteor[1:],0)))
1845
1862
1846 listMeteors.append(numpy.array([i,indInit,indPeak,indEnd,FLA])) #CHEQUEAR!!!!!
1863 listMeteors.append(numpy.array([i,indInit,indPeak,indEnd,FLA])) #CHEQUEAR!!!!!
1847 j = numpy.where(indUPthresh == indEnd)[0] + 1
1864 j = numpy.where(indUPthresh == indEnd)[0] + 1
1848 else: j+=1
1865 else: j+=1
1849 else: j+=1
1866 else: j+=1
1850
1867
1851 return listMeteors
1868 return listMeteors
1852
1869
1853 def __removeMultipleDetections(self,listMeteors, rangeLimit, timeLimit):
1870 def __removeMultipleDetections(self,listMeteors, rangeLimit, timeLimit):
1854
1871
1855 arrayMeteors = numpy.asarray(listMeteors)
1872 arrayMeteors = numpy.asarray(listMeteors)
1856 listMeteors1 = []
1873 listMeteors1 = []
1857
1874
1858 while arrayMeteors.shape[0] > 0:
1875 while arrayMeteors.shape[0] > 0:
1859 FLAs = arrayMeteors[:,4]
1876 FLAs = arrayMeteors[:,4]
1860 maxFLA = FLAs.argmax()
1877 maxFLA = FLAs.argmax()
1861 listMeteors1.append(arrayMeteors[maxFLA,:])
1878 listMeteors1.append(arrayMeteors[maxFLA,:])
1862
1879
1863 MeteorInitTime = arrayMeteors[maxFLA,1]
1880 MeteorInitTime = arrayMeteors[maxFLA,1]
1864 MeteorEndTime = arrayMeteors[maxFLA,3]
1881 MeteorEndTime = arrayMeteors[maxFLA,3]
1865 MeteorHeight = arrayMeteors[maxFLA,0]
1882 MeteorHeight = arrayMeteors[maxFLA,0]
1866
1883
1867 #Check neighborhood
1884 #Check neighborhood
1868 maxHeightIndex = MeteorHeight + rangeLimit
1885 maxHeightIndex = MeteorHeight + rangeLimit
1869 minHeightIndex = MeteorHeight - rangeLimit
1886 minHeightIndex = MeteorHeight - rangeLimit
1870 minTimeIndex = MeteorInitTime - timeLimit
1887 minTimeIndex = MeteorInitTime - timeLimit
1871 maxTimeIndex = MeteorEndTime + timeLimit
1888 maxTimeIndex = MeteorEndTime + timeLimit
1872
1889
1873 #Check Heights
1890 #Check Heights
1874 indHeight = numpy.logical_and(arrayMeteors[:,0] >= minHeightIndex, arrayMeteors[:,0] <= maxHeightIndex)
1891 indHeight = numpy.logical_and(arrayMeteors[:,0] >= minHeightIndex, arrayMeteors[:,0] <= maxHeightIndex)
1875 indTime = numpy.logical_and(arrayMeteors[:,3] >= minTimeIndex, arrayMeteors[:,1] <= maxTimeIndex)
1892 indTime = numpy.logical_and(arrayMeteors[:,3] >= minTimeIndex, arrayMeteors[:,1] <= maxTimeIndex)
1876 indBoth = numpy.where(numpy.logical_and(indTime,indHeight))
1893 indBoth = numpy.where(numpy.logical_and(indTime,indHeight))
1877
1894
1878 arrayMeteors = numpy.delete(arrayMeteors, indBoth, axis = 0)
1895 arrayMeteors = numpy.delete(arrayMeteors, indBoth, axis = 0)
1879
1896
1880 return listMeteors1
1897 return listMeteors1
1881
1898
1882 def __meteorReestimation(self, listMeteors, volts, pairslist, thresh, noise, timeInterval,frequency):
1899 def __meteorReestimation(self, listMeteors, volts, pairslist, thresh, noise, timeInterval,frequency):
1883 numHeights = volts.shape[2]
1900 numHeights = volts.shape[2]
1884 nChannel = volts.shape[0]
1901 nChannel = volts.shape[0]
1885
1902
1886 thresholdPhase = thresh[0]
1903 thresholdPhase = thresh[0]
1887 thresholdNoise = thresh[1]
1904 thresholdNoise = thresh[1]
1888 thresholdDB = float(thresh[2])
1905 thresholdDB = float(thresh[2])
1889
1906
1890 thresholdDB1 = 10**(thresholdDB/10)
1907 thresholdDB1 = 10**(thresholdDB/10)
1891 pairsarray = numpy.array(pairslist)
1908 pairsarray = numpy.array(pairslist)
1892 indSides = pairsarray[:,1]
1909 indSides = pairsarray[:,1]
1893
1910
1894 pairslist1 = list(pairslist)
1911 pairslist1 = list(pairslist)
1895 pairslist1.append((0,1))
1912 pairslist1.append((0,1))
1896 pairslist1.append((3,4))
1913 pairslist1.append((3,4))
1897
1914
1898 listMeteors1 = []
1915 listMeteors1 = []
1899 listPowerSeries = []
1916 listPowerSeries = []
1900 listVoltageSeries = []
1917 listVoltageSeries = []
1901 #volts has the war data
1918 #volts has the war data
1902
1919
1903 if frequency == 30e6:
1920 if frequency == 30e6:
1904 timeLag = 45*10**-3
1921 timeLag = 45*10**-3
1905 else:
1922 else:
1906 timeLag = 15*10**-3
1923 timeLag = 15*10**-3
1907 lag = numpy.ceil(timeLag/timeInterval)
1924 lag = numpy.ceil(timeLag/timeInterval)
1908
1925
1909 for i in range(len(listMeteors)):
1926 for i in range(len(listMeteors)):
1910
1927
1911 ###################### 3.6 - 3.7 PARAMETERS REESTIMATION #########################
1928 ###################### 3.6 - 3.7 PARAMETERS REESTIMATION #########################
1912 meteorAux = numpy.zeros(16)
1929 meteorAux = numpy.zeros(16)
1913
1930
1914 #Loading meteor Data (mHeight, mStart, mPeak, mEnd)
1931 #Loading meteor Data (mHeight, mStart, mPeak, mEnd)
1915 mHeight = listMeteors[i][0]
1932 mHeight = listMeteors[i][0]
1916 mStart = listMeteors[i][1]
1933 mStart = listMeteors[i][1]
1917 mPeak = listMeteors[i][2]
1934 mPeak = listMeteors[i][2]
1918 mEnd = listMeteors[i][3]
1935 mEnd = listMeteors[i][3]
1919
1936
1920 #get the volt data between the start and end times of the meteor
1937 #get the volt data between the start and end times of the meteor
1921 meteorVolts = volts[:,mStart:mEnd+1,mHeight]
1938 meteorVolts = volts[:,mStart:mEnd+1,mHeight]
1922 meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
1939 meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
1923
1940
1924 #3.6. Phase Difference estimation
1941 #3.6. Phase Difference estimation
1925 phaseDiff, aux = self.__estimatePhaseDifference(meteorVolts, pairslist)
1942 phaseDiff, aux = self.__estimatePhaseDifference(meteorVolts, pairslist)
1926
1943
1927 #3.7. Phase difference removal & meteor start, peak and end times reestimated
1944 #3.7. Phase difference removal & meteor start, peak and end times reestimated
1928 #meteorVolts0.- all Channels, all Profiles
1945 #meteorVolts0.- all Channels, all Profiles
1929 meteorVolts0 = volts[:,:,mHeight]
1946 meteorVolts0 = volts[:,:,mHeight]
1930 meteorThresh = noise[:,mHeight]*thresholdNoise
1947 meteorThresh = noise[:,mHeight]*thresholdNoise
1931 meteorNoise = noise[:,mHeight]
1948 meteorNoise = noise[:,mHeight]
1932 meteorVolts0[indSides,:] = self.__shiftPhase(meteorVolts0[indSides,:], phaseDiff) #Phase Shifting
1949 meteorVolts0[indSides,:] = self.__shiftPhase(meteorVolts0[indSides,:], phaseDiff) #Phase Shifting
1933 powerNet0 = numpy.nansum(numpy.abs(meteorVolts0)**2, axis = 0) #Power
1950 powerNet0 = numpy.nansum(numpy.abs(meteorVolts0)**2, axis = 0) #Power
1934
1951
1935 #Times reestimation
1952 #Times reestimation
1936 mStart1 = numpy.where(powerNet0[:mPeak] < meteorThresh[:mPeak])[0]
1953 mStart1 = numpy.where(powerNet0[:mPeak] < meteorThresh[:mPeak])[0]
1937 if mStart1.size > 0:
1954 if mStart1.size > 0:
1938 mStart1 = mStart1[-1] + 1
1955 mStart1 = mStart1[-1] + 1
1939
1956
1940 else:
1957 else:
1941 mStart1 = mPeak
1958 mStart1 = mPeak
1942
1959
1943 mEnd1 = numpy.where(powerNet0[mPeak:] < meteorThresh[mPeak:])[0][0] + mPeak - 1
1960 mEnd1 = numpy.where(powerNet0[mPeak:] < meteorThresh[mPeak:])[0][0] + mPeak - 1
1944 mEndDecayTime1 = numpy.where(powerNet0[mPeak:] < meteorNoise[mPeak:])[0]
1961 mEndDecayTime1 = numpy.where(powerNet0[mPeak:] < meteorNoise[mPeak:])[0]
1945 if mEndDecayTime1.size == 0:
1962 if mEndDecayTime1.size == 0:
1946 mEndDecayTime1 = powerNet0.size
1963 mEndDecayTime1 = powerNet0.size
1947 else:
1964 else:
1948 mEndDecayTime1 = mEndDecayTime1[0] + mPeak - 1
1965 mEndDecayTime1 = mEndDecayTime1[0] + mPeak - 1
1949 # mPeak1 = meteorVolts0[mStart1:mEnd1 + 1].argmax()
1966 # mPeak1 = meteorVolts0[mStart1:mEnd1 + 1].argmax()
1950
1967
1951 #meteorVolts1.- all Channels, from start to end
1968 #meteorVolts1.- all Channels, from start to end
1952 meteorVolts1 = meteorVolts0[:,mStart1:mEnd1 + 1]
1969 meteorVolts1 = meteorVolts0[:,mStart1:mEnd1 + 1]
1953 meteorVolts2 = meteorVolts0[:,mPeak + lag:mEnd1 + 1]
1970 meteorVolts2 = meteorVolts0[:,mPeak + lag:mEnd1 + 1]
1954 if meteorVolts2.shape[1] == 0:
1971 if meteorVolts2.shape[1] == 0:
1955 meteorVolts2 = meteorVolts0[:,mPeak:mEnd1 + 1]
1972 meteorVolts2 = meteorVolts0[:,mPeak:mEnd1 + 1]
1956 meteorVolts1 = meteorVolts1.reshape(meteorVolts1.shape[0], meteorVolts1.shape[1], 1)
1973 meteorVolts1 = meteorVolts1.reshape(meteorVolts1.shape[0], meteorVolts1.shape[1], 1)
1957 meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1], 1)
1974 meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1], 1)
1958 ##################### END PARAMETERS REESTIMATION #########################
1975 ##################### END PARAMETERS REESTIMATION #########################
1959
1976
1960 ##################### 3.8 PHASE DIFFERENCE REESTIMATION ########################
1977 ##################### 3.8 PHASE DIFFERENCE REESTIMATION ########################
1961 # if mEnd1 - mStart1 > 4: #Error Number 6: echo less than 5 samples long; too short for analysis
1978 # if mEnd1 - mStart1 > 4: #Error Number 6: echo less than 5 samples long; too short for analysis
1962 if meteorVolts2.shape[1] > 0:
1979 if meteorVolts2.shape[1] > 0:
1963 #Phase Difference re-estimation
1980 #Phase Difference re-estimation
1964 phaseDiff1, phaseDiffint = self.__estimatePhaseDifference(meteorVolts2, pairslist1) #Phase Difference Estimation
1981 phaseDiff1, phaseDiffint = self.__estimatePhaseDifference(meteorVolts2, pairslist1) #Phase Difference Estimation
1965 # phaseDiff1, phaseDiffint = self.estimatePhaseDifference(meteorVolts2, pairslist)
1982 # phaseDiff1, phaseDiffint = self.estimatePhaseDifference(meteorVolts2, pairslist)
1966 meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1])
1983 meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1])
1967 phaseDiff11 = numpy.reshape(phaseDiff1, (phaseDiff1.shape[0],1))
1984 phaseDiff11 = numpy.reshape(phaseDiff1, (phaseDiff1.shape[0],1))
1968 meteorVolts2[indSides,:] = self.__shiftPhase(meteorVolts2[indSides,:], phaseDiff11[0:4]) #Phase Shifting
1985 meteorVolts2[indSides,:] = self.__shiftPhase(meteorVolts2[indSides,:], phaseDiff11[0:4]) #Phase Shifting
1969
1986
1970 #Phase Difference RMS
1987 #Phase Difference RMS
1971 phaseRMS1 = numpy.sqrt(numpy.mean(numpy.square(phaseDiff1)))
1988 phaseRMS1 = numpy.sqrt(numpy.mean(numpy.square(phaseDiff1)))
1972 powerNet1 = numpy.nansum(numpy.abs(meteorVolts1[:,:])**2,0)
1989 powerNet1 = numpy.nansum(numpy.abs(meteorVolts1[:,:])**2,0)
1973 #Data from Meteor
1990 #Data from Meteor
1974 mPeak1 = powerNet1.argmax() + mStart1
1991 mPeak1 = powerNet1.argmax() + mStart1
1975 mPeakPower1 = powerNet1.max()
1992 mPeakPower1 = powerNet1.max()
1976 noiseAux = sum(noise[mStart1:mEnd1 + 1,mHeight])
1993 noiseAux = sum(noise[mStart1:mEnd1 + 1,mHeight])
1977 mSNR1 = (sum(powerNet1)-noiseAux)/noiseAux
1994 mSNR1 = (sum(powerNet1)-noiseAux)/noiseAux
1978 Meteor1 = numpy.array([mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1])
1995 Meteor1 = numpy.array([mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1])
1979 Meteor1 = numpy.hstack((Meteor1,phaseDiffint))
1996 Meteor1 = numpy.hstack((Meteor1,phaseDiffint))
1980 PowerSeries = powerNet0[mStart1:mEndDecayTime1 + 1]
1997 PowerSeries = powerNet0[mStart1:mEndDecayTime1 + 1]
1981 #Vectorize
1998 #Vectorize
1982 meteorAux[0:7] = [mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1]
1999 meteorAux[0:7] = [mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1]
1983 meteorAux[7:11] = phaseDiffint[0:4]
2000 meteorAux[7:11] = phaseDiffint[0:4]
1984
2001
1985 #Rejection Criterions
2002 #Rejection Criterions
1986 if phaseRMS1 > thresholdPhase: #Error Number 17: Phase variation
2003 if phaseRMS1 > thresholdPhase: #Error Number 17: Phase variation
1987 meteorAux[-1] = 17
2004 meteorAux[-1] = 17
1988 elif mSNR1 < thresholdDB1: #Error Number 1: SNR < threshold dB
2005 elif mSNR1 < thresholdDB1: #Error Number 1: SNR < threshold dB
1989 meteorAux[-1] = 1
2006 meteorAux[-1] = 1
1990
2007
1991
2008
1992 else:
2009 else:
1993 meteorAux[0:4] = [mHeight, mStart, mPeak, mEnd]
2010 meteorAux[0:4] = [mHeight, mStart, mPeak, mEnd]
1994 meteorAux[-1] = 6 #Error Number 6: echo less than 5 samples long; too short for analysis
2011 meteorAux[-1] = 6 #Error Number 6: echo less than 5 samples long; too short for analysis
1995 PowerSeries = 0
2012 PowerSeries = 0
1996
2013
1997 listMeteors1.append(meteorAux)
2014 listMeteors1.append(meteorAux)
1998 listPowerSeries.append(PowerSeries)
2015 listPowerSeries.append(PowerSeries)
1999 listVoltageSeries.append(meteorVolts1)
2016 listVoltageSeries.append(meteorVolts1)
2000
2017
2001 return listMeteors1, listPowerSeries, listVoltageSeries
2018 return listMeteors1, listPowerSeries, listVoltageSeries
2002
2019
2003 def __estimateDecayTime(self, listMeteors, listPower, timeInterval, frequency):
2020 def __estimateDecayTime(self, listMeteors, listPower, timeInterval, frequency):
2004
2021
2005 threshError = 10
2022 threshError = 10
2006 #Depending if it is 30 or 50 MHz
2023 #Depending if it is 30 or 50 MHz
2007 if frequency == 30e6:
2024 if frequency == 30e6:
2008 timeLag = 45*10**-3
2025 timeLag = 45*10**-3
2009 else:
2026 else:
2010 timeLag = 15*10**-3
2027 timeLag = 15*10**-3
2011 lag = numpy.ceil(timeLag/timeInterval)
2028 lag = numpy.ceil(timeLag/timeInterval)
2012
2029
2013 listMeteors1 = []
2030 listMeteors1 = []
2014
2031
2015 for i in range(len(listMeteors)):
2032 for i in range(len(listMeteors)):
2016 meteorPower = listPower[i]
2033 meteorPower = listPower[i]
2017 meteorAux = listMeteors[i]
2034 meteorAux = listMeteors[i]
2018
2035
2019 if meteorAux[-1] == 0:
2036 if meteorAux[-1] == 0:
2020
2037
2021 try:
2038 try:
2022 indmax = meteorPower.argmax()
2039 indmax = meteorPower.argmax()
2023 indlag = indmax + lag
2040 indlag = indmax + lag
2024
2041
2025 y = meteorPower[indlag:]
2042 y = meteorPower[indlag:]
2026 x = numpy.arange(0, y.size)*timeLag
2043 x = numpy.arange(0, y.size)*timeLag
2027
2044
2028 #first guess
2045 #first guess
2029 a = y[0]
2046 a = y[0]
2030 tau = timeLag
2047 tau = timeLag
2031 #exponential fit
2048 #exponential fit
2032 popt, pcov = optimize.curve_fit(self.__exponential_function, x, y, p0 = [a, tau])
2049 popt, pcov = optimize.curve_fit(self.__exponential_function, x, y, p0 = [a, tau])
2033 y1 = self.__exponential_function(x, *popt)
2050 y1 = self.__exponential_function(x, *popt)
2034 #error estimation
2051 #error estimation
2035 error = sum((y - y1)**2)/(numpy.var(y)*(y.size - popt.size))
2052 error = sum((y - y1)**2)/(numpy.var(y)*(y.size - popt.size))
2036
2053
2037 decayTime = popt[1]
2054 decayTime = popt[1]
2038 riseTime = indmax*timeInterval
2055 riseTime = indmax*timeInterval
2039 meteorAux[11:13] = [decayTime, error]
2056 meteorAux[11:13] = [decayTime, error]
2040
2057
2041 #Table items 7, 8 and 11
2058 #Table items 7, 8 and 11
2042 if (riseTime > 0.3): #Number 7: Echo rise exceeds 0.3s
2059 if (riseTime > 0.3): #Number 7: Echo rise exceeds 0.3s
2043 meteorAux[-1] = 7
2060 meteorAux[-1] = 7
2044 elif (decayTime < 2*riseTime) : #Number 8: Echo decay time less than than twice rise time
2061 elif (decayTime < 2*riseTime) : #Number 8: Echo decay time less than than twice rise time
2045 meteorAux[-1] = 8
2062 meteorAux[-1] = 8
2046 if (error > threshError): #Number 11: Poor fit to amplitude for estimation of decay time
2063 if (error > threshError): #Number 11: Poor fit to amplitude for estimation of decay time
2047 meteorAux[-1] = 11
2064 meteorAux[-1] = 11
2048
2065
2049
2066
2050 except:
2067 except:
2051 meteorAux[-1] = 11
2068 meteorAux[-1] = 11
2052
2069
2053
2070
2054 listMeteors1.append(meteorAux)
2071 listMeteors1.append(meteorAux)
2055
2072
2056 return listMeteors1
2073 return listMeteors1
2057
2074
2058 #Exponential Function
2075 #Exponential Function
2059
2076
2060 def __exponential_function(self, x, a, tau):
2077 def __exponential_function(self, x, a, tau):
2061 y = a*numpy.exp(-x/tau)
2078 y = a*numpy.exp(-x/tau)
2062 return y
2079 return y
2063
2080
2064 def __getRadialVelocity(self, listMeteors, listVolts, radialStdThresh, pairslist, timeInterval):
2081 def __getRadialVelocity(self, listMeteors, listVolts, radialStdThresh, pairslist, timeInterval):
2065
2082
2066 pairslist1 = list(pairslist)
2083 pairslist1 = list(pairslist)
2067 pairslist1.append((0,1))
2084 pairslist1.append((0,1))
2068 pairslist1.append((3,4))
2085 pairslist1.append((3,4))
2069 numPairs = len(pairslist1)
2086 numPairs = len(pairslist1)
2070 #Time Lag
2087 #Time Lag
2071 timeLag = 45*10**-3
2088 timeLag = 45*10**-3
2072 c = 3e8
2089 c = 3e8
2073 lag = numpy.ceil(timeLag/timeInterval)
2090 lag = numpy.ceil(timeLag/timeInterval)
2074 freq = 30e6
2091 freq = 30e6
2075
2092
2076 listMeteors1 = []
2093 listMeteors1 = []
2077
2094
2078 for i in range(len(listMeteors)):
2095 for i in range(len(listMeteors)):
2079 meteorAux = listMeteors[i]
2096 meteorAux = listMeteors[i]
2080 if meteorAux[-1] == 0:
2097 if meteorAux[-1] == 0:
2081 mStart = listMeteors[i][1]
2098 mStart = listMeteors[i][1]
2082 mPeak = listMeteors[i][2]
2099 mPeak = listMeteors[i][2]
2083 mLag = mPeak - mStart + lag
2100 mLag = mPeak - mStart + lag
2084
2101
2085 #get the volt data between the start and end times of the meteor
2102 #get the volt data between the start and end times of the meteor
2086 meteorVolts = listVolts[i]
2103 meteorVolts = listVolts[i]
2087 meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
2104 meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
2088
2105
2089 #Get CCF
2106 #Get CCF
2090 allCCFs = self.__calculateCCF(meteorVolts, pairslist1, [-2,-1,0,1,2])
2107 allCCFs = self.__calculateCCF(meteorVolts, pairslist1, [-2,-1,0,1,2])
2091
2108
2092 #Method 2
2109 #Method 2
2093 slopes = numpy.zeros(numPairs)
2110 slopes = numpy.zeros(numPairs)
2094 time = numpy.array([-2,-1,1,2])*timeInterval
2111 time = numpy.array([-2,-1,1,2])*timeInterval
2095 angAllCCF = numpy.angle(allCCFs[:,[0,1,3,4],0])
2112 angAllCCF = numpy.angle(allCCFs[:,[0,1,3,4],0])
2096
2113
2097 #Correct phases
2114 #Correct phases
2098 derPhaseCCF = angAllCCF[:,1:] - angAllCCF[:,0:-1]
2115 derPhaseCCF = angAllCCF[:,1:] - angAllCCF[:,0:-1]
2099 indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
2116 indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
2100
2117
2101 if indDer[0].shape[0] > 0:
2118 if indDer[0].shape[0] > 0:
2102 for i in range(indDer[0].shape[0]):
2119 for i in range(indDer[0].shape[0]):
2103 signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i]])
2120 signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i]])
2104 angAllCCF[indDer[0][i],indDer[1][i]+1:] += signo*2*numpy.pi
2121 angAllCCF[indDer[0][i],indDer[1][i]+1:] += signo*2*numpy.pi
2105
2122
2106 # fit = scipy.stats.linregress(numpy.array([-2,-1,1,2])*timeInterval, numpy.array([phaseLagN2s[i],phaseLagN1s[i],phaseLag1s[i],phaseLag2s[i]]))
2123 # fit = scipy.stats.linregress(numpy.array([-2,-1,1,2])*timeInterval, numpy.array([phaseLagN2s[i],phaseLagN1s[i],phaseLag1s[i],phaseLag2s[i]]))
2107 for j in range(numPairs):
2124 for j in range(numPairs):
2108 fit = stats.linregress(time, angAllCCF[j,:])
2125 fit = stats.linregress(time, angAllCCF[j,:])
2109 slopes[j] = fit[0]
2126 slopes[j] = fit[0]
2110
2127
2111 #Remove Outlier
2128 #Remove Outlier
2112 # indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
2129 # indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
2113 # slopes = numpy.delete(slopes,indOut)
2130 # slopes = numpy.delete(slopes,indOut)
2114 # indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
2131 # indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
2115 # slopes = numpy.delete(slopes,indOut)
2132 # slopes = numpy.delete(slopes,indOut)
2116
2133
2117 radialVelocity = -numpy.mean(slopes)*(0.25/numpy.pi)*(c/freq)
2134 radialVelocity = -numpy.mean(slopes)*(0.25/numpy.pi)*(c/freq)
2118 radialError = numpy.std(slopes)*(0.25/numpy.pi)*(c/freq)
2135 radialError = numpy.std(slopes)*(0.25/numpy.pi)*(c/freq)
2119 meteorAux[-2] = radialError
2136 meteorAux[-2] = radialError
2120 meteorAux[-3] = radialVelocity
2137 meteorAux[-3] = radialVelocity
2121
2138
2122 #Setting Error
2139 #Setting Error
2123 #Number 15: Radial Drift velocity or projected horizontal velocity exceeds 200 m/s
2140 #Number 15: Radial Drift velocity or projected horizontal velocity exceeds 200 m/s
2124 if numpy.abs(radialVelocity) > 200:
2141 if numpy.abs(radialVelocity) > 200:
2125 meteorAux[-1] = 15
2142 meteorAux[-1] = 15
2126 #Number 12: Poor fit to CCF variation for estimation of radial drift velocity
2143 #Number 12: Poor fit to CCF variation for estimation of radial drift velocity
2127 elif radialError > radialStdThresh:
2144 elif radialError > radialStdThresh:
2128 meteorAux[-1] = 12
2145 meteorAux[-1] = 12
2129
2146
2130 listMeteors1.append(meteorAux)
2147 listMeteors1.append(meteorAux)
2131 return listMeteors1
2148 return listMeteors1
2132
2149
2133 def __setNewArrays(self, listMeteors, date, heiRang):
2150 def __setNewArrays(self, listMeteors, date, heiRang):
2134
2151
2135 #New arrays
2152 #New arrays
2136 arrayMeteors = numpy.array(listMeteors)
2153 arrayMeteors = numpy.array(listMeteors)
2137 arrayParameters = numpy.zeros((len(listMeteors), 13))
2154 arrayParameters = numpy.zeros((len(listMeteors), 13))
2138
2155
2139 #Date inclusion
2156 #Date inclusion
2140 # date = re.findall(r'\((.*?)\)', date)
2157 # date = re.findall(r'\((.*?)\)', date)
2141 # date = date[0].split(',')
2158 # date = date[0].split(',')
2142 # date = map(int, date)
2159 # date = map(int, date)
2143 #
2160 #
2144 # if len(date)<6:
2161 # if len(date)<6:
2145 # date.append(0)
2162 # date.append(0)
2146 #
2163 #
2147 # date = [date[0]*10000 + date[1]*100 + date[2], date[3]*10000 + date[4]*100 + date[5]]
2164 # date = [date[0]*10000 + date[1]*100 + date[2], date[3]*10000 + date[4]*100 + date[5]]
2148 # arrayDate = numpy.tile(date, (len(listMeteors), 1))
2165 # arrayDate = numpy.tile(date, (len(listMeteors), 1))
2149 arrayDate = numpy.tile(date, (len(listMeteors)))
2166 arrayDate = numpy.tile(date, (len(listMeteors)))
2150
2167
2151 #Meteor array
2168 #Meteor array
2152 # arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)]
2169 # arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)]
2153 # arrayMeteors = numpy.hstack((arrayDate, arrayMeteors))
2170 # arrayMeteors = numpy.hstack((arrayDate, arrayMeteors))
2154
2171
2155 #Parameters Array
2172 #Parameters Array
2156 arrayParameters[:,0] = arrayDate #Date
2173 arrayParameters[:,0] = arrayDate #Date
2157 arrayParameters[:,1] = heiRang[arrayMeteors[:,0].astype(int)] #Range
2174 arrayParameters[:,1] = heiRang[arrayMeteors[:,0].astype(int)] #Range
2158 arrayParameters[:,6:8] = arrayMeteors[:,-3:-1] #Radial velocity and its error
2175 arrayParameters[:,6:8] = arrayMeteors[:,-3:-1] #Radial velocity and its error
2159 arrayParameters[:,8:12] = arrayMeteors[:,7:11] #Phases
2176 arrayParameters[:,8:12] = arrayMeteors[:,7:11] #Phases
2160 arrayParameters[:,-1] = arrayMeteors[:,-1] #Error
2177 arrayParameters[:,-1] = arrayMeteors[:,-1] #Error
2161
2178
2162
2179
2163 return arrayParameters
2180 return arrayParameters
2164
2181
2165 class CorrectMeteorPhases(Operation):
2182 class CorrectSMPhases(Operation):
2166
2183
2167 def run(self, dataOut, phaseOffsets, hmin = 50, hmax = 150, azimuth = 45, channelPositions = None):
2184 def run(self, dataOut, phaseOffsets, hmin = 50, hmax = 150, azimuth = 45, channelPositions = None):
2168
2185
2169 arrayParameters = dataOut.data_param
2186 arrayParameters = dataOut.data_param
2170 pairsList = []
2187 pairsList = []
2171 pairx = (0,1)
2188 pairx = (0,1)
2172 pairy = (2,3)
2189 pairy = (2,3)
2173 pairsList.append(pairx)
2190 pairsList.append(pairx)
2174 pairsList.append(pairy)
2191 pairsList.append(pairy)
2175 jph = numpy.zeros(4)
2192 jph = numpy.zeros(4)
2176
2193
2177 phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180
2194 phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180
2178 # arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets)
2195 # arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets)
2179 arrayParameters[:,8:12] = numpy.angle(numpy.exp(1j*(arrayParameters[:,8:12] + phaseOffsets)))
2196 arrayParameters[:,8:12] = numpy.angle(numpy.exp(1j*(arrayParameters[:,8:12] + phaseOffsets)))
2180
2197
2181 meteorOps = MeteorOperations()
2198 meteorOps = SMOperations()
2182 if channelPositions == None:
2199 if channelPositions == None:
2183 # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
2200 # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
2184 channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella
2201 channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella
2185
2202
2186 pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
2203 pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
2187 h = (hmin,hmax)
2204 h = (hmin,hmax)
2188
2205
2189 arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, jph)
2206 arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, jph)
2190
2207
2191 dataOut.data_param = arrayParameters
2208 dataOut.data_param = arrayParameters
2192 return
2209 return
2193
2210
2194 class PhaseCalibration(Operation):
2211 class SMPhaseCalibration(Operation):
2195
2212
2196 __buffer = None
2213 __buffer = None
2197
2214
2198 __initime = None
2215 __initime = None
2199
2216
2200 __dataReady = False
2217 __dataReady = False
2201
2218
2202 __isConfig = False
2219 __isConfig = False
2203
2220
2204 def __checkTime(self, currentTime, initTime, paramInterval, outputInterval):
2221 def __checkTime(self, currentTime, initTime, paramInterval, outputInterval):
2205
2222
2206 dataTime = currentTime + paramInterval
2223 dataTime = currentTime + paramInterval
2207 deltaTime = dataTime - initTime
2224 deltaTime = dataTime - initTime
2208
2225
2209 if deltaTime >= outputInterval or deltaTime < 0:
2226 if deltaTime >= outputInterval or deltaTime < 0:
2210 return True
2227 return True
2211
2228
2212 return False
2229 return False
2213
2230
2214 def __getGammas(self, pairs, d, phases):
2231 def __getGammas(self, pairs, d, phases):
2215 gammas = numpy.zeros(2)
2232 gammas = numpy.zeros(2)
2216
2233
2217 for i in range(len(pairs)):
2234 for i in range(len(pairs)):
2218
2235
2219 pairi = pairs[i]
2236 pairi = pairs[i]
2220
2237
2221 phip3 = phases[:,pairi[1]]
2238 phip3 = phases[:,pairi[1]]
2222 d3 = d[pairi[1]]
2239 d3 = d[pairi[1]]
2223 phip2 = phases[:,pairi[0]]
2240 phip2 = phases[:,pairi[0]]
2224 d2 = d[pairi[0]]
2241 d2 = d[pairi[0]]
2225 #Calculating gamma
2242 #Calculating gamma
2226 # jdcos = alp1/(k*d1)
2243 # jdcos = alp1/(k*d1)
2227 # jgamma = numpy.angle(numpy.exp(1j*(d0*alp1/d1 - alp0)))
2244 # jgamma = numpy.angle(numpy.exp(1j*(d0*alp1/d1 - alp0)))
2228 jgamma = -phip2*d3/d2 - phip3
2245 jgamma = -phip2*d3/d2 - phip3
2229 jgamma = numpy.angle(numpy.exp(1j*jgamma))
2246 jgamma = numpy.angle(numpy.exp(1j*jgamma))
2230 # jgamma[jgamma>numpy.pi] -= 2*numpy.pi
2247 # jgamma[jgamma>numpy.pi] -= 2*numpy.pi
2231 # jgamma[jgamma<-numpy.pi] += 2*numpy.pi
2248 # jgamma[jgamma<-numpy.pi] += 2*numpy.pi
2232
2249
2233 #Revised distribution
2250 #Revised distribution
2234 jgammaArray = numpy.hstack((jgamma,jgamma+0.5*numpy.pi,jgamma-0.5*numpy.pi))
2251 jgammaArray = numpy.hstack((jgamma,jgamma+0.5*numpy.pi,jgamma-0.5*numpy.pi))
2235
2252
2236 #Histogram
2253 #Histogram
2237 nBins = 64.0
2254 nBins = 64.0
2238 rmin = -0.5*numpy.pi
2255 rmin = -0.5*numpy.pi
2239 rmax = 0.5*numpy.pi
2256 rmax = 0.5*numpy.pi
2240 phaseHisto = numpy.histogram(jgammaArray, bins=nBins, range=(rmin,rmax))
2257 phaseHisto = numpy.histogram(jgammaArray, bins=nBins, range=(rmin,rmax))
2241
2258
2242 meteorsY = phaseHisto[0]
2259 meteorsY = phaseHisto[0]
2243 phasesX = phaseHisto[1][:-1]
2260 phasesX = phaseHisto[1][:-1]
2244 width = phasesX[1] - phasesX[0]
2261 width = phasesX[1] - phasesX[0]
2245 phasesX += width/2
2262 phasesX += width/2
2246
2263
2247 #Gaussian aproximation
2264 #Gaussian aproximation
2248 bpeak = meteorsY.argmax()
2265 bpeak = meteorsY.argmax()
2249 peak = meteorsY.max()
2266 peak = meteorsY.max()
2250 jmin = bpeak - 5
2267 jmin = bpeak - 5
2251 jmax = bpeak + 5 + 1
2268 jmax = bpeak + 5 + 1
2252
2269
2253 if jmin<0:
2270 if jmin<0:
2254 jmin = 0
2271 jmin = 0
2255 jmax = 6
2272 jmax = 6
2256 elif jmax > meteorsY.size:
2273 elif jmax > meteorsY.size:
2257 jmin = meteorsY.size - 6
2274 jmin = meteorsY.size - 6
2258 jmax = meteorsY.size
2275 jmax = meteorsY.size
2259
2276
2260 x0 = numpy.array([peak,bpeak,50])
2277 x0 = numpy.array([peak,bpeak,50])
2261 coeff = optimize.leastsq(self.__residualFunction, x0, args=(meteorsY[jmin:jmax], phasesX[jmin:jmax]))
2278 coeff = optimize.leastsq(self.__residualFunction, x0, args=(meteorsY[jmin:jmax], phasesX[jmin:jmax]))
2262
2279
2263 #Gammas
2280 #Gammas
2264 gammas[i] = coeff[0][1]
2281 gammas[i] = coeff[0][1]
2265
2282
2266 return gammas
2283 return gammas
2267
2284
2268 def __residualFunction(self, coeffs, y, t):
2285 def __residualFunction(self, coeffs, y, t):
2269
2286
2270 return y - self.__gauss_function(t, coeffs)
2287 return y - self.__gauss_function(t, coeffs)
2271
2288
2272 def __gauss_function(self, t, coeffs):
2289 def __gauss_function(self, t, coeffs):
2273
2290
2274 return coeffs[0]*numpy.exp(-0.5*((t - coeffs[1]) / coeffs[2])**2)
2291 return coeffs[0]*numpy.exp(-0.5*((t - coeffs[1]) / coeffs[2])**2)
2275
2292
2276 def __getPhases(self, azimuth, h, pairsList, d, gammas, meteorsArray):
2293 def __getPhases(self, azimuth, h, pairsList, d, gammas, meteorsArray):
2277 meteorOps = MeteorOperations()
2294 meteorOps = SMOperations()
2278 nchan = 4
2295 nchan = 4
2279 pairx = pairsList[0]
2296 pairx = pairsList[0]
2280 pairy = pairsList[1]
2297 pairy = pairsList[1]
2281 center_xangle = 0
2298 center_xangle = 0
2282 center_yangle = 0
2299 center_yangle = 0
2283 range_angle = numpy.array([10*numpy.pi,numpy.pi,numpy.pi/2,numpy.pi/4])
2300 range_angle = numpy.array([10*numpy.pi,numpy.pi,numpy.pi/2,numpy.pi/4])
2284 ntimes = len(range_angle)
2301 ntimes = len(range_angle)
2285
2302
2286 nstepsx = 20.0
2303 nstepsx = 20.0
2287 nstepsy = 20.0
2304 nstepsy = 20.0
2288
2305
2289 for iz in range(ntimes):
2306 for iz in range(ntimes):
2290 min_xangle = -range_angle[iz]/2 + center_xangle
2307 min_xangle = -range_angle[iz]/2 + center_xangle
2291 max_xangle = range_angle[iz]/2 + center_xangle
2308 max_xangle = range_angle[iz]/2 + center_xangle
2292 min_yangle = -range_angle[iz]/2 + center_yangle
2309 min_yangle = -range_angle[iz]/2 + center_yangle
2293 max_yangle = range_angle[iz]/2 + center_yangle
2310 max_yangle = range_angle[iz]/2 + center_yangle
2294
2311
2295 inc_x = (max_xangle-min_xangle)/nstepsx
2312 inc_x = (max_xangle-min_xangle)/nstepsx
2296 inc_y = (max_yangle-min_yangle)/nstepsy
2313 inc_y = (max_yangle-min_yangle)/nstepsy
2297
2314
2298 alpha_y = numpy.arange(nstepsy)*inc_y + min_yangle
2315 alpha_y = numpy.arange(nstepsy)*inc_y + min_yangle
2299 alpha_x = numpy.arange(nstepsx)*inc_x + min_xangle
2316 alpha_x = numpy.arange(nstepsx)*inc_x + min_xangle
2300 penalty = numpy.zeros((nstepsx,nstepsy))
2317 penalty = numpy.zeros((nstepsx,nstepsy))
2301 jph_array = numpy.zeros((nchan,nstepsx,nstepsy))
2318 jph_array = numpy.zeros((nchan,nstepsx,nstepsy))
2302 jph = numpy.zeros(nchan)
2319 jph = numpy.zeros(nchan)
2303
2320
2304 # Iterations looking for the offset
2321 # Iterations looking for the offset
2305 for iy in range(int(nstepsy)):
2322 for iy in range(int(nstepsy)):
2306 for ix in range(int(nstepsx)):
2323 for ix in range(int(nstepsx)):
2307 jph[pairy[1]] = alpha_y[iy]
2324 jph[pairy[1]] = alpha_y[iy]
2308 jph[pairy[0]] = -gammas[1] - alpha_y[iy]*d[pairy[1]]/d[pairy[0]]
2325 jph[pairy[0]] = -gammas[1] - alpha_y[iy]*d[pairy[1]]/d[pairy[0]]
2309
2326
2310 jph[pairx[1]] = alpha_x[ix]
2327 jph[pairx[1]] = alpha_x[ix]
2311 jph[pairx[0]] = -gammas[0] - alpha_x[ix]*d[pairx[1]]/d[pairx[0]]
2328 jph[pairx[0]] = -gammas[0] - alpha_x[ix]*d[pairx[1]]/d[pairx[0]]
2312
2329
2313 jph_array[:,ix,iy] = jph
2330 jph_array[:,ix,iy] = jph
2314
2331
2315 meteorsArray1 = meteorOps.getMeteorParams(meteorsArray, azimuth, h, pairsList, d, jph)
2332 meteorsArray1 = meteorOps.getMeteorParams(meteorsArray, azimuth, h, pairsList, d, jph)
2316 error = meteorsArray1[:,-1]
2333 error = meteorsArray1[:,-1]
2317 ind1 = numpy.where(error==0)[0]
2334 ind1 = numpy.where(error==0)[0]
2318 penalty[ix,iy] = ind1.size
2335 penalty[ix,iy] = ind1.size
2319
2336
2320 i,j = numpy.unravel_index(penalty.argmax(), penalty.shape)
2337 i,j = numpy.unravel_index(penalty.argmax(), penalty.shape)
2321 phOffset = jph_array[:,i,j]
2338 phOffset = jph_array[:,i,j]
2322
2339
2323 center_xangle = phOffset[pairx[1]]
2340 center_xangle = phOffset[pairx[1]]
2324 center_yangle = phOffset[pairy[1]]
2341 center_yangle = phOffset[pairy[1]]
2325
2342
2326 phOffset = numpy.angle(numpy.exp(1j*jph_array[:,i,j]))
2343 phOffset = numpy.angle(numpy.exp(1j*jph_array[:,i,j]))
2327 phOffset = phOffset*180/numpy.pi
2344 phOffset = phOffset*180/numpy.pi
2328 return phOffset
2345 return phOffset
2329
2346
2330
2347
2331 def run(self, dataOut, hmin, hmax, channelPositions=None, nHours = 1):
2348 def run(self, dataOut, hmin, hmax, channelPositions=None, nHours = 1):
2332
2349
2333 dataOut.flagNoData = True
2350 dataOut.flagNoData = True
2334 self.__dataReady = False
2351 self.__dataReady = False
2335 dataOut.outputInterval = nHours*3600
2352 dataOut.outputInterval = nHours*3600
2336
2353
2337 if self.__isConfig == False:
2354 if self.__isConfig == False:
2338 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
2355 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
2339 #Get Initial LTC time
2356 #Get Initial LTC time
2340 self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
2357 self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
2341 self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
2358 self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
2342
2359
2343 self.__isConfig = True
2360 self.__isConfig = True
2344
2361
2345 if self.__buffer == None:
2362 if self.__buffer == None:
2346 self.__buffer = dataOut.data_param.copy()
2363 self.__buffer = dataOut.data_param.copy()
2347
2364
2348 else:
2365 else:
2349 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
2366 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
2350
2367
2351 self.__dataReady = self.__checkTime(dataOut.utctime, self.__initime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
2368 self.__dataReady = self.__checkTime(dataOut.utctime, self.__initime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
2352
2369
2353 if self.__dataReady:
2370 if self.__dataReady:
2354 dataOut.utctimeInit = self.__initime
2371 dataOut.utctimeInit = self.__initime
2355 self.__initime += dataOut.outputInterval #to erase time offset
2372 self.__initime += dataOut.outputInterval #to erase time offset
2356
2373
2357 freq = dataOut.frequency
2374 freq = dataOut.frequency
2358 c = dataOut.C #m/s
2375 c = dataOut.C #m/s
2359 lamb = c/freq
2376 lamb = c/freq
2360 k = 2*numpy.pi/lamb
2377 k = 2*numpy.pi/lamb
2361 azimuth = 0
2378 azimuth = 0
2362 h = (hmin, hmax)
2379 h = (hmin, hmax)
2363 pairs = ((0,1),(2,3))
2380 pairs = ((0,1),(2,3))
2364
2381
2365 if channelPositions == None:
2382 if channelPositions == None:
2366 # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
2383 # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
2367 channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella
2384 channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella
2368 meteorOps = MeteorOperations()
2385 meteorOps = SMOperations()
2369 pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
2386 pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
2370
2387
2371 # distances1 = [-distances[0]*lamb, distances[1]*lamb, -distances[2]*lamb, distances[3]*lamb]
2388 # distances1 = [-distances[0]*lamb, distances[1]*lamb, -distances[2]*lamb, distances[3]*lamb]
2372
2389
2373 meteorsArray = self.__buffer
2390 meteorsArray = self.__buffer
2374 error = meteorsArray[:,-1]
2391 error = meteorsArray[:,-1]
2375 boolError = (error==0)|(error==3)|(error==4)|(error==13)|(error==14)
2392 boolError = (error==0)|(error==3)|(error==4)|(error==13)|(error==14)
2376 ind1 = numpy.where(boolError)[0]
2393 ind1 = numpy.where(boolError)[0]
2377 meteorsArray = meteorsArray[ind1,:]
2394 meteorsArray = meteorsArray[ind1,:]
2378 meteorsArray[:,-1] = 0
2395 meteorsArray[:,-1] = 0
2379 phases = meteorsArray[:,8:12]
2396 phases = meteorsArray[:,8:12]
2380
2397
2381 #Calculate Gammas
2398 #Calculate Gammas
2382 gammas = self.__getGammas(pairs, distances, phases)
2399 gammas = self.__getGammas(pairs, distances, phases)
2383 # gammas = numpy.array([-21.70409463,45.76935864])*numpy.pi/180
2400 # gammas = numpy.array([-21.70409463,45.76935864])*numpy.pi/180
2384 #Calculate Phases
2401 #Calculate Phases
2385 phasesOff = self.__getPhases(azimuth, h, pairs, distances, gammas, meteorsArray)
2402 phasesOff = self.__getPhases(azimuth, h, pairs, distances, gammas, meteorsArray)
2386 phasesOff = phasesOff.reshape((1,phasesOff.size))
2403 phasesOff = phasesOff.reshape((1,phasesOff.size))
2387 dataOut.data_output = -phasesOff
2404 dataOut.data_output = -phasesOff
2388 dataOut.flagNoData = False
2405 dataOut.flagNoData = False
2389 self.__buffer = None
2406 self.__buffer = None
2390
2407
2391
2408
2392 return
2409 return
2393
2410
2394 class MeteorOperations():
2411 class SMOperations():
2395
2412
2396 def __init__(self):
2413 def __init__(self):
2397
2414
2398 return
2415 return
2399
2416
2400 def getMeteorParams(self, arrayParameters0, azimuth, h, pairsList, distances, jph):
2417 def getMeteorParams(self, arrayParameters0, azimuth, h, pairsList, distances, jph):
2401
2418
2402 arrayParameters = arrayParameters0.copy()
2419 arrayParameters = arrayParameters0.copy()
2403 hmin = h[0]
2420 hmin = h[0]
2404 hmax = h[1]
2421 hmax = h[1]
2405
2422
2406 #Calculate AOA (Error N 3, 4)
2423 #Calculate AOA (Error N 3, 4)
2407 #JONES ET AL. 1998
2424 #JONES ET AL. 1998
2408 AOAthresh = numpy.pi/8
2425 AOAthresh = numpy.pi/8
2409 error = arrayParameters[:,-1]
2426 error = arrayParameters[:,-1]
2410 phases = -arrayParameters[:,8:12] + jph
2427 phases = -arrayParameters[:,8:12] + jph
2411 # phases = numpy.unwrap(phases)
2428 # phases = numpy.unwrap(phases)
2412 arrayParameters[:,3:6], arrayParameters[:,-1] = self.__getAOA(phases, pairsList, distances, error, AOAthresh, azimuth)
2429 arrayParameters[:,3:6], arrayParameters[:,-1] = self.__getAOA(phases, pairsList, distances, error, AOAthresh, azimuth)
2413
2430
2414 #Calculate Heights (Error N 13 and 14)
2431 #Calculate Heights (Error N 13 and 14)
2415 error = arrayParameters[:,-1]
2432 error = arrayParameters[:,-1]
2416 Ranges = arrayParameters[:,1]
2433 Ranges = arrayParameters[:,1]
2417 zenith = arrayParameters[:,4]
2434 zenith = arrayParameters[:,4]
2418 arrayParameters[:,2], arrayParameters[:,-1] = self.__getHeights(Ranges, zenith, error, hmin, hmax)
2435 arrayParameters[:,2], arrayParameters[:,-1] = self.__getHeights(Ranges, zenith, error, hmin, hmax)
2419
2436
2420 #----------------------- Get Final data ------------------------------------
2437 #----------------------- Get Final data ------------------------------------
2421 # error = arrayParameters[:,-1]
2438 # error = arrayParameters[:,-1]
2422 # ind1 = numpy.where(error==0)[0]
2439 # ind1 = numpy.where(error==0)[0]
2423 # arrayParameters = arrayParameters[ind1,:]
2440 # arrayParameters = arrayParameters[ind1,:]
2424
2441
2425 return arrayParameters
2442 return arrayParameters
2426
2443
2427 def __getAOA(self, phases, pairsList, directions, error, AOAthresh, azimuth):
2444 def __getAOA(self, phases, pairsList, directions, error, AOAthresh, azimuth):
2428
2445
2429 arrayAOA = numpy.zeros((phases.shape[0],3))
2446 arrayAOA = numpy.zeros((phases.shape[0],3))
2430 cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList,directions)
2447 cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList,directions)
2431
2448
2432 arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
2449 arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
2433 cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
2450 cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
2434 arrayAOA[:,2] = cosDirError
2451 arrayAOA[:,2] = cosDirError
2435
2452
2436 azimuthAngle = arrayAOA[:,0]
2453 azimuthAngle = arrayAOA[:,0]
2437 zenithAngle = arrayAOA[:,1]
2454 zenithAngle = arrayAOA[:,1]
2438
2455
2439 #Setting Error
2456 #Setting Error
2440 indError = numpy.where(numpy.logical_or(error == 3, error == 4))[0]
2457 indError = numpy.where(numpy.logical_or(error == 3, error == 4))[0]
2441 error[indError] = 0
2458 error[indError] = 0
2442 #Number 3: AOA not fesible
2459 #Number 3: AOA not fesible
2443 indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
2460 indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
2444 error[indInvalid] = 3
2461 error[indInvalid] = 3
2445 #Number 4: Large difference in AOAs obtained from different antenna baselines
2462 #Number 4: Large difference in AOAs obtained from different antenna baselines
2446 indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
2463 indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
2447 error[indInvalid] = 4
2464 error[indInvalid] = 4
2448 return arrayAOA, error
2465 return arrayAOA, error
2449
2466
2450 def __getDirectionCosines(self, arrayPhase, pairsList, distances):
2467 def __getDirectionCosines(self, arrayPhase, pairsList, distances):
2451
2468
2452 #Initializing some variables
2469 #Initializing some variables
2453 ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
2470 ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
2454 ang_aux = ang_aux.reshape(1,ang_aux.size)
2471 ang_aux = ang_aux.reshape(1,ang_aux.size)
2455
2472
2456 cosdir = numpy.zeros((arrayPhase.shape[0],2))
2473 cosdir = numpy.zeros((arrayPhase.shape[0],2))
2457 cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
2474 cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
2458
2475
2459
2476
2460 for i in range(2):
2477 for i in range(2):
2461 ph0 = arrayPhase[:,pairsList[i][0]]
2478 ph0 = arrayPhase[:,pairsList[i][0]]
2462 ph1 = arrayPhase[:,pairsList[i][1]]
2479 ph1 = arrayPhase[:,pairsList[i][1]]
2463 d0 = distances[pairsList[i][0]]
2480 d0 = distances[pairsList[i][0]]
2464 d1 = distances[pairsList[i][1]]
2481 d1 = distances[pairsList[i][1]]
2465
2482
2466 ph0_aux = ph0 + ph1
2483 ph0_aux = ph0 + ph1
2467 ph0_aux = numpy.angle(numpy.exp(1j*ph0_aux))
2484 ph0_aux = numpy.angle(numpy.exp(1j*ph0_aux))
2468 # ph0_aux[ph0_aux > numpy.pi] -= 2*numpy.pi
2485 # ph0_aux[ph0_aux > numpy.pi] -= 2*numpy.pi
2469 # ph0_aux[ph0_aux < -numpy.pi] += 2*numpy.pi
2486 # ph0_aux[ph0_aux < -numpy.pi] += 2*numpy.pi
2470 #First Estimation
2487 #First Estimation
2471 cosdir0[:,i] = (ph0_aux)/(2*numpy.pi*(d0 - d1))
2488 cosdir0[:,i] = (ph0_aux)/(2*numpy.pi*(d0 - d1))
2472
2489
2473 #Most-Accurate Second Estimation
2490 #Most-Accurate Second Estimation
2474 phi1_aux = ph0 - ph1
2491 phi1_aux = ph0 - ph1
2475 phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
2492 phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
2476 #Direction Cosine 1
2493 #Direction Cosine 1
2477 cosdir1 = (phi1_aux + ang_aux)/(2*numpy.pi*(d0 + d1))
2494 cosdir1 = (phi1_aux + ang_aux)/(2*numpy.pi*(d0 + d1))
2478
2495
2479 #Searching the correct Direction Cosine
2496 #Searching the correct Direction Cosine
2480 cosdir0_aux = cosdir0[:,i]
2497 cosdir0_aux = cosdir0[:,i]
2481 cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
2498 cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
2482 #Minimum Distance
2499 #Minimum Distance
2483 cosDiff = (cosdir1 - cosdir0_aux)**2
2500 cosDiff = (cosdir1 - cosdir0_aux)**2
2484 indcos = cosDiff.argmin(axis = 1)
2501 indcos = cosDiff.argmin(axis = 1)
2485 #Saving Value obtained
2502 #Saving Value obtained
2486 cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
2503 cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
2487
2504
2488 return cosdir0, cosdir
2505 return cosdir0, cosdir
2489
2506
2490 def __calculateAOA(self, cosdir, azimuth):
2507 def __calculateAOA(self, cosdir, azimuth):
2491 cosdirX = cosdir[:,0]
2508 cosdirX = cosdir[:,0]
2492 cosdirY = cosdir[:,1]
2509 cosdirY = cosdir[:,1]
2493
2510
2494 zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
2511 zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
2495 azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth#0 deg north, 90 deg east
2512 azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth#0 deg north, 90 deg east
2496 angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
2513 angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
2497
2514
2498 return angles
2515 return angles
2499
2516
2500 def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
2517 def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
2501
2518
2502 Ramb = 375 #Ramb = c/(2*PRF)
2519 Ramb = 375 #Ramb = c/(2*PRF)
2503 Re = 6371 #Earth Radius
2520 Re = 6371 #Earth Radius
2504 heights = numpy.zeros(Ranges.shape)
2521 heights = numpy.zeros(Ranges.shape)
2505
2522
2506 R_aux = numpy.array([0,1,2])*Ramb
2523 R_aux = numpy.array([0,1,2])*Ramb
2507 R_aux = R_aux.reshape(1,R_aux.size)
2524 R_aux = R_aux.reshape(1,R_aux.size)
2508
2525
2509 Ranges = Ranges.reshape(Ranges.size,1)
2526 Ranges = Ranges.reshape(Ranges.size,1)
2510
2527
2511 Ri = Ranges + R_aux
2528 Ri = Ranges + R_aux
2512 hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re
2529 hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re
2513
2530
2514 #Check if there is a height between 70 and 110 km
2531 #Check if there is a height between 70 and 110 km
2515 h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
2532 h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
2516 ind_h = numpy.where(h_bool == 1)[0]
2533 ind_h = numpy.where(h_bool == 1)[0]
2517
2534
2518 hCorr = hi[ind_h, :]
2535 hCorr = hi[ind_h, :]
2519 ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
2536 ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
2520
2537
2521 hCorr = hi[ind_hCorr]
2538 hCorr = hi[ind_hCorr]
2522 heights[ind_h] = hCorr
2539 heights[ind_h] = hCorr
2523
2540
2524 #Setting Error
2541 #Setting Error
2525 #Number 13: Height unresolvable echo: not valid height within 70 to 110 km
2542 #Number 13: Height unresolvable echo: not valid height within 70 to 110 km
2526 #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
2543 #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
2527 indError = numpy.where(numpy.logical_or(error == 13, error == 14))[0]
2544 indError = numpy.where(numpy.logical_or(error == 13, error == 14))[0]
2528 error[indError] = 0
2545 error[indError] = 0
2529 indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
2546 indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
2530 error[indInvalid2] = 14
2547 error[indInvalid2] = 14
2531 indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
2548 indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
2532 error[indInvalid1] = 13
2549 error[indInvalid1] = 13
2533
2550
2534 return heights, error
2551 return heights, error
2535
2552
2536 def getPhasePairs(self, channelPositions):
2553 def getPhasePairs(self, channelPositions):
2537 chanPos = numpy.array(channelPositions)
2554 chanPos = numpy.array(channelPositions)
2538 listOper = list(itertools.combinations(range(5),2))
2555 listOper = list(itertools.combinations(range(5),2))
2539
2556
2540 distances = numpy.zeros(4)
2557 distances = numpy.zeros(4)
2541 axisX = []
2558 axisX = []
2542 axisY = []
2559 axisY = []
2543 distX = numpy.zeros(3)
2560 distX = numpy.zeros(3)
2544 distY = numpy.zeros(3)
2561 distY = numpy.zeros(3)
2545 ix = 0
2562 ix = 0
2546 iy = 0
2563 iy = 0
2547
2564
2548 pairX = numpy.zeros((2,2))
2565 pairX = numpy.zeros((2,2))
2549 pairY = numpy.zeros((2,2))
2566 pairY = numpy.zeros((2,2))
2550
2567
2551 for i in range(len(listOper)):
2568 for i in range(len(listOper)):
2552 pairi = listOper[i]
2569 pairi = listOper[i]
2553
2570
2554 posDif = numpy.abs(chanPos[pairi[0],:] - chanPos[pairi[1],:])
2571 posDif = numpy.abs(chanPos[pairi[0],:] - chanPos[pairi[1],:])
2555
2572
2556 if posDif[0] == 0:
2573 if posDif[0] == 0:
2557 axisY.append(pairi)
2574 axisY.append(pairi)
2558 distY[iy] = posDif[1]
2575 distY[iy] = posDif[1]
2559 iy += 1
2576 iy += 1
2560 elif posDif[1] == 0:
2577 elif posDif[1] == 0:
2561 axisX.append(pairi)
2578 axisX.append(pairi)
2562 distX[ix] = posDif[0]
2579 distX[ix] = posDif[0]
2563 ix += 1
2580 ix += 1
2564
2581
2565 for i in range(2):
2582 for i in range(2):
2566 if i==0:
2583 if i==0:
2567 dist0 = distX
2584 dist0 = distX
2568 axis0 = axisX
2585 axis0 = axisX
2569 else:
2586 else:
2570 dist0 = distY
2587 dist0 = distY
2571 axis0 = axisY
2588 axis0 = axisY
2572
2589
2573 side = numpy.argsort(dist0)[:-1]
2590 side = numpy.argsort(dist0)[:-1]
2574 axis0 = numpy.array(axis0)[side,:]
2591 axis0 = numpy.array(axis0)[side,:]
2575 chanC = int(numpy.intersect1d(axis0[0,:], axis0[1,:])[0])
2592 chanC = int(numpy.intersect1d(axis0[0,:], axis0[1,:])[0])
2576 axis1 = numpy.unique(numpy.reshape(axis0,4))
2593 axis1 = numpy.unique(numpy.reshape(axis0,4))
2577 side = axis1[axis1 != chanC]
2594 side = axis1[axis1 != chanC]
2578 diff1 = chanPos[chanC,i] - chanPos[side[0],i]
2595 diff1 = chanPos[chanC,i] - chanPos[side[0],i]
2579 diff2 = chanPos[chanC,i] - chanPos[side[1],i]
2596 diff2 = chanPos[chanC,i] - chanPos[side[1],i]
2580 if diff1<0:
2597 if diff1<0:
2581 chan2 = side[0]
2598 chan2 = side[0]
2582 d2 = numpy.abs(diff1)
2599 d2 = numpy.abs(diff1)
2583 chan1 = side[1]
2600 chan1 = side[1]
2584 d1 = numpy.abs(diff2)
2601 d1 = numpy.abs(diff2)
2585 else:
2602 else:
2586 chan2 = side[1]
2603 chan2 = side[1]
2587 d2 = numpy.abs(diff2)
2604 d2 = numpy.abs(diff2)
2588 chan1 = side[0]
2605 chan1 = side[0]
2589 d1 = numpy.abs(diff1)
2606 d1 = numpy.abs(diff1)
2590
2607
2591 if i==0:
2608 if i==0:
2592 chanCX = chanC
2609 chanCX = chanC
2593 chan1X = chan1
2610 chan1X = chan1
2594 chan2X = chan2
2611 chan2X = chan2
2595 distances[0:2] = numpy.array([d1,d2])
2612 distances[0:2] = numpy.array([d1,d2])
2596 else:
2613 else:
2597 chanCY = chanC
2614 chanCY = chanC
2598 chan1Y = chan1
2615 chan1Y = chan1
2599 chan2Y = chan2
2616 chan2Y = chan2
2600 distances[2:4] = numpy.array([d1,d2])
2617 distances[2:4] = numpy.array([d1,d2])
2601 # axisXsides = numpy.reshape(axisX[ix,:],4)
2618 # axisXsides = numpy.reshape(axisX[ix,:],4)
2602 #
2619 #
2603 # channelCentX = int(numpy.intersect1d(pairX[0,:], pairX[1,:])[0])
2620 # channelCentX = int(numpy.intersect1d(pairX[0,:], pairX[1,:])[0])
2604 # channelCentY = int(numpy.intersect1d(pairY[0,:], pairY[1,:])[0])
2621 # channelCentY = int(numpy.intersect1d(pairY[0,:], pairY[1,:])[0])
2605 #
2622 #
2606 # ind25X = numpy.where(pairX[0,:] != channelCentX)[0][0]
2623 # ind25X = numpy.where(pairX[0,:] != channelCentX)[0][0]
2607 # ind20X = numpy.where(pairX[1,:] != channelCentX)[0][0]
2624 # ind20X = numpy.where(pairX[1,:] != channelCentX)[0][0]
2608 # channel25X = int(pairX[0,ind25X])
2625 # channel25X = int(pairX[0,ind25X])
2609 # channel20X = int(pairX[1,ind20X])
2626 # channel20X = int(pairX[1,ind20X])
2610 # ind25Y = numpy.where(pairY[0,:] != channelCentY)[0][0]
2627 # ind25Y = numpy.where(pairY[0,:] != channelCentY)[0][0]
2611 # ind20Y = numpy.where(pairY[1,:] != channelCentY)[0][0]
2628 # ind20Y = numpy.where(pairY[1,:] != channelCentY)[0][0]
2612 # channel25Y = int(pairY[0,ind25Y])
2629 # channel25Y = int(pairY[0,ind25Y])
2613 # channel20Y = int(pairY[1,ind20Y])
2630 # channel20Y = int(pairY[1,ind20Y])
2614
2631
2615 # pairslist = [(channelCentX, channel25X),(channelCentX, channel20X),(channelCentY,channel25Y),(channelCentY, channel20Y)]
2632 # pairslist = [(channelCentX, channel25X),(channelCentX, channel20X),(channelCentY,channel25Y),(channelCentY, channel20Y)]
2616 pairslist = [(chanCX, chan1X),(chanCX, chan2X),(chanCY,chan1Y),(chanCY, chan2Y)]
2633 pairslist = [(chanCX, chan1X),(chanCX, chan2X),(chanCY,chan1Y),(chanCY, chan2Y)]
2617
2634
2618 return pairslist, distances
2635 return pairslist, distances
2619 # def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth):
2636 # def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth):
2620 #
2637 #
2621 # arrayAOA = numpy.zeros((phases.shape[0],3))
2638 # arrayAOA = numpy.zeros((phases.shape[0],3))
2622 # cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList)
2639 # cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList)
2623 #
2640 #
2624 # arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
2641 # arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
2625 # cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
2642 # cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
2626 # arrayAOA[:,2] = cosDirError
2643 # arrayAOA[:,2] = cosDirError
2627 #
2644 #
2628 # azimuthAngle = arrayAOA[:,0]
2645 # azimuthAngle = arrayAOA[:,0]
2629 # zenithAngle = arrayAOA[:,1]
2646 # zenithAngle = arrayAOA[:,1]
2630 #
2647 #
2631 # #Setting Error
2648 # #Setting Error
2632 # #Number 3: AOA not fesible
2649 # #Number 3: AOA not fesible
2633 # indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
2650 # indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
2634 # error[indInvalid] = 3
2651 # error[indInvalid] = 3
2635 # #Number 4: Large difference in AOAs obtained from different antenna baselines
2652 # #Number 4: Large difference in AOAs obtained from different antenna baselines
2636 # indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
2653 # indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
2637 # error[indInvalid] = 4
2654 # error[indInvalid] = 4
2638 # return arrayAOA, error
2655 # return arrayAOA, error
2639 #
2656 #
2640 # def __getDirectionCosines(self, arrayPhase, pairsList):
2657 # def __getDirectionCosines(self, arrayPhase, pairsList):
2641 #
2658 #
2642 # #Initializing some variables
2659 # #Initializing some variables
2643 # ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
2660 # ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
2644 # ang_aux = ang_aux.reshape(1,ang_aux.size)
2661 # ang_aux = ang_aux.reshape(1,ang_aux.size)
2645 #
2662 #
2646 # cosdir = numpy.zeros((arrayPhase.shape[0],2))
2663 # cosdir = numpy.zeros((arrayPhase.shape[0],2))
2647 # cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
2664 # cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
2648 #
2665 #
2649 #
2666 #
2650 # for i in range(2):
2667 # for i in range(2):
2651 # #First Estimation
2668 # #First Estimation
2652 # phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]]
2669 # phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]]
2653 # #Dealias
2670 # #Dealias
2654 # indcsi = numpy.where(phi0_aux > numpy.pi)
2671 # indcsi = numpy.where(phi0_aux > numpy.pi)
2655 # phi0_aux[indcsi] -= 2*numpy.pi
2672 # phi0_aux[indcsi] -= 2*numpy.pi
2656 # indcsi = numpy.where(phi0_aux < -numpy.pi)
2673 # indcsi = numpy.where(phi0_aux < -numpy.pi)
2657 # phi0_aux[indcsi] += 2*numpy.pi
2674 # phi0_aux[indcsi] += 2*numpy.pi
2658 # #Direction Cosine 0
2675 # #Direction Cosine 0
2659 # cosdir0[:,i] = -(phi0_aux)/(2*numpy.pi*0.5)
2676 # cosdir0[:,i] = -(phi0_aux)/(2*numpy.pi*0.5)
2660 #
2677 #
2661 # #Most-Accurate Second Estimation
2678 # #Most-Accurate Second Estimation
2662 # phi1_aux = arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]]
2679 # phi1_aux = arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]]
2663 # phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
2680 # phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
2664 # #Direction Cosine 1
2681 # #Direction Cosine 1
2665 # cosdir1 = -(phi1_aux + ang_aux)/(2*numpy.pi*4.5)
2682 # cosdir1 = -(phi1_aux + ang_aux)/(2*numpy.pi*4.5)
2666 #
2683 #
2667 # #Searching the correct Direction Cosine
2684 # #Searching the correct Direction Cosine
2668 # cosdir0_aux = cosdir0[:,i]
2685 # cosdir0_aux = cosdir0[:,i]
2669 # cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
2686 # cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
2670 # #Minimum Distance
2687 # #Minimum Distance
2671 # cosDiff = (cosdir1 - cosdir0_aux)**2
2688 # cosDiff = (cosdir1 - cosdir0_aux)**2
2672 # indcos = cosDiff.argmin(axis = 1)
2689 # indcos = cosDiff.argmin(axis = 1)
2673 # #Saving Value obtained
2690 # #Saving Value obtained
2674 # cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
2691 # cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
2675 #
2692 #
2676 # return cosdir0, cosdir
2693 # return cosdir0, cosdir
2677 #
2694 #
2678 # def __calculateAOA(self, cosdir, azimuth):
2695 # def __calculateAOA(self, cosdir, azimuth):
2679 # cosdirX = cosdir[:,0]
2696 # cosdirX = cosdir[:,0]
2680 # cosdirY = cosdir[:,1]
2697 # cosdirY = cosdir[:,1]
2681 #
2698 #
2682 # zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
2699 # zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
2683 # azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth #0 deg north, 90 deg east
2700 # azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth #0 deg north, 90 deg east
2684 # angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
2701 # angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
2685 #
2702 #
2686 # return angles
2703 # return angles
2687 #
2704 #
2688 # def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
2705 # def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
2689 #
2706 #
2690 # Ramb = 375 #Ramb = c/(2*PRF)
2707 # Ramb = 375 #Ramb = c/(2*PRF)
2691 # Re = 6371 #Earth Radius
2708 # Re = 6371 #Earth Radius
2692 # heights = numpy.zeros(Ranges.shape)
2709 # heights = numpy.zeros(Ranges.shape)
2693 #
2710 #
2694 # R_aux = numpy.array([0,1,2])*Ramb
2711 # R_aux = numpy.array([0,1,2])*Ramb
2695 # R_aux = R_aux.reshape(1,R_aux.size)
2712 # R_aux = R_aux.reshape(1,R_aux.size)
2696 #
2713 #
2697 # Ranges = Ranges.reshape(Ranges.size,1)
2714 # Ranges = Ranges.reshape(Ranges.size,1)
2698 #
2715 #
2699 # Ri = Ranges + R_aux
2716 # Ri = Ranges + R_aux
2700 # hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re
2717 # hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re
2701 #
2718 #
2702 # #Check if there is a height between 70 and 110 km
2719 # #Check if there is a height between 70 and 110 km
2703 # h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
2720 # h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
2704 # ind_h = numpy.where(h_bool == 1)[0]
2721 # ind_h = numpy.where(h_bool == 1)[0]
2705 #
2722 #
2706 # hCorr = hi[ind_h, :]
2723 # hCorr = hi[ind_h, :]
2707 # ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
2724 # ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
2708 #
2725 #
2709 # hCorr = hi[ind_hCorr]
2726 # hCorr = hi[ind_hCorr]
2710 # heights[ind_h] = hCorr
2727 # heights[ind_h] = hCorr
2711 #
2728 #
2712 # #Setting Error
2729 # #Setting Error
2713 # #Number 13: Height unresolvable echo: not valid height within 70 to 110 km
2730 # #Number 13: Height unresolvable echo: not valid height within 70 to 110 km
2714 # #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
2731 # #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
2715 #
2732 #
2716 # indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
2733 # indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
2717 # error[indInvalid2] = 14
2734 # error[indInvalid2] = 14
2718 # indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
2735 # indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
2719 # error[indInvalid1] = 13
2736 # error[indInvalid1] = 13
2720 #
2737 #
2721 # return heights, error
2738 # return heights, error
2722 No newline at end of file
2739
General Comments 0
You need to be logged in to leave comments. Login now