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