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