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