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