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