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