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