##// END OF EJS Templates
schain
RSS
Correction to Spaced Antenna libraries
Julio Valdez -
r762:7c90b920ac51
parent child
Show More
Show file before | Show file after | Hide comments
@@ -1,246 +1,246
1 1 import numpy
2 2
3 3 from jroproc_base import ProcessingUnit, Operation
4 4 from schainpy.model.data.jrodata import Correlation
5 5
6 6 class CorrelationProc(ProcessingUnit):
7 7
8 8 def __init__(self):
9 9
10 10 ProcessingUnit.__init__(self)
11 11
12 12 self.objectDict = {}
13 13 self.buffer = None
14 14 self.firstdatatime = None
15 15 self.profIndex = 0
16 16 self.dataOut = Correlation()
17 17
18 18 def __updateObjFromVoltage(self):
19 19
20 20 self.dataOut.timeZone = self.dataIn.timeZone
21 21 self.dataOut.dstFlag = self.dataIn.dstFlag
22 22 self.dataOut.errorCount = self.dataIn.errorCount
23 23 self.dataOut.useLocalTime = self.dataIn.useLocalTime
24 24
25 25 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
26 26 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
27 27 self.dataOut.channelList = self.dataIn.channelList
28 28 self.dataOut.heightList = self.dataIn.heightList
29 29 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
30 30 # self.dataOut.nHeights = self.dataIn.nHeights
31 31 # self.dataOut.nChannels = self.dataIn.nChannels
32 32 self.dataOut.nBaud = self.dataIn.nBaud
33 33 self.dataOut.nCode = self.dataIn.nCode
34 34 self.dataOut.code = self.dataIn.code
35 35 # self.dataOut.nProfiles = self.dataOut.nFFTPoints
36 36 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
37 37 self.dataOut.utctime = self.firstdatatime
38 38 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
39 39 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
40 # self.dataOut.nCohInt = self.dataIn.nCohInt
40 self.dataOut.nCohInt = self.dataIn.nCohInt
41 41 # self.dataOut.nIncohInt = 1
42 42 self.dataOut.ippSeconds = self.dataIn.ippSeconds
43 43 # self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
44 44
45 45 # self.dataOut.timeInterval = self.dataIn.timeInterval*self.dataOut.nPoints
46 46
47 47
48 48 def removeDC(self, jspectra):
49 49
50 50 nChannel = jspectra.shape[0]
51 51
52 52 for i in range(nChannel):
53 53 jspectra_tmp = jspectra[i,:,:]
54 54 jspectra_DC = numpy.mean(jspectra_tmp,axis = 0)
55 55
56 56 jspectra_tmp = jspectra_tmp - jspectra_DC
57 57 jspectra[i,:,:] = jspectra_tmp
58 58
59 59 return jspectra
60 60
61 61
62 62 def removeNoise(self, mode = 2):
63 63 indR = numpy.where(self.dataOut.lagR == 0)[0][0]
64 64 indT = numpy.where(self.dataOut.lagT == 0)[0][0]
65 65
66 66 jspectra = self.dataOut.data_corr[:,:,indR,:]
67 67
68 68 num_chan = jspectra.shape[0]
69 69 num_hei = jspectra.shape[2]
70 70
71 71 freq_dc = indT
72 72 ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
73 73
74 74 NPot = self.dataOut.getNoise(mode)
75 75 jspectra[:,freq_dc,:] = jspectra[:,freq_dc,:] - NPot
76 76 SPot = jspectra[:,freq_dc,:]
77 77 pairsAutoCorr = self.dataOut.getPairsAutoCorr()
78 78 # self.dataOut.signalPotency = SPot
79 79 self.dataOut.noise = NPot
80 80 self.dataOut.SNR = (SPot/NPot)[pairsAutoCorr]
81 81 self.dataOut.data_corr[:,:,indR,:] = jspectra
82 82
83 83 return 1
84 84
85 85
86 86 def calculateNormFactor(self):
87 87
88 88 pairsList = self.dataOut.pairsList
89 89 pairsAutoCorr = self.dataOut.pairsAutoCorr
90 90 nHeights = self.dataOut.nHeights
91 91 nPairs = len(pairsList)
92 92 normFactor = numpy.zeros((nPairs,nHeights))
93 93
94 94 indR = numpy.where(self.dataOut.lagR == 0)[0][0]
95 95 indT = numpy.where(self.dataOut.lagT == 0)[0][0]
96 96
97 97 for l in range(len(pairsList)):
98 98 firstChannel = pairsList[l][0]
99 99 secondChannel = pairsList[l][1]
100 100
101 101 AC1 = pairsAutoCorr[firstChannel]
102 102 AC2 = pairsAutoCorr[secondChannel]
103 103
104 104 if (AC1 >= 0 and AC2 >= 0):
105 105
106 106 data1 = numpy.abs(self.dataOut.data_corr[AC1,:,indR,:])
107 107 data2 = numpy.abs(self.dataOut.data_corr[AC2,:,indR,:])
108 108 maxim1 = data1.max(axis = 0)
109 109 maxim2 = data1.max(axis = 0)
110 110 maxim = numpy.sqrt(maxim1*maxim2)
111 111 else:
112 112 #In case there is no autocorrelation for the pair
113 113 data = numpy.abs(self.dataOut.data_corr[l,:,indR,:])
114 114 maxim = numpy.max(data, axis = 0)
115 115
116 116 normFactor[l,:] = maxim
117 117
118 118 self.dataOut.normFactor = normFactor
119 119
120 120 return 1
121 121
122 122 def run(self, lagT=None, lagR=None, pairsList=None,
123 123 nPoints=None, nAvg=None, bufferSize=None,
124 124 fullT = False, fullR = False, removeDC = False):
125 125
126 126 self.dataOut.flagNoData = True
127 127
128 128 if self.dataIn.type == "Correlation":
129 129
130 130 self.dataOut.copy(self.dataIn)
131 131
132 132 return
133 133
134 134 if self.dataIn.type == "Voltage":
135 135
136 136 if pairsList == None:
137 137 pairsList = [numpy.array([0,0])]
138 138
139 139 if nPoints == None:
140 140 nPoints = 128
141 141 #------------------------------------------------------------
142 142 #Condicionales para calcular Correlaciones en Tiempo y Rango
143 143 if fullT:
144 144 lagT = numpy.arange(nPoints*2 - 1) - nPoints + 1
145 145 elif lagT == None:
146 146 lagT = numpy.array([0])
147 147 else:
148 148 lagT = numpy.array(lagT)
149 149
150 150 if fullR:
151 151 lagR = numpy.arange(self.dataOut.nHeights)
152 152 elif lagR == None:
153 153 lagR = numpy.array([0])
154 154 #-------------------------------------------------------------
155 155
156 156 if nAvg == None:
157 157 nAvg = 1
158 158
159 159 if bufferSize == None:
160 160 bufferSize = 0
161 161
162 162 deltaH = self.dataIn.heightList[1] - self.dataIn.heightList[0]
163 163 self.dataOut.lagR = numpy.round(numpy.array(lagR)/deltaH)
164 164 self.dataOut.pairsList = pairsList
165 165 self.dataOut.nPoints = nPoints
166 166 # channels = numpy.sort(list(set(list(itertools.chain.from_iterable(pairsList)))))
167 167
168 168 if self.buffer is None:
169 169
170 170 self.buffer = numpy.zeros((self.dataIn.nChannels,self.dataIn.nProfiles,self.dataIn.nHeights),dtype='complex')
171 171
172 172
173 173 self.buffer[:,self.profIndex,:] = self.dataIn.data.copy()[:,:]
174 174
175 175 self.profIndex += 1
176 176
177 177 if self.firstdatatime == None:
178 178
179 179 self.firstdatatime = self.dataIn.utctime
180 180
181 181 if self.profIndex == nPoints:
182 182
183 183 tmp = self.buffer[:,0:nPoints,:]
184 184 self.buffer = None
185 185 self.buffer = tmp
186 186
187 187 #--------------- Remover DC ------------
188 188 if removeDC:
189 189 self.buffer = self.removeDC(self.buffer)
190 190 #---------------------------------------------
191 191 self.dataOut.data_volts = self.buffer
192 192 self.__updateObjFromVoltage()
193 193 self.dataOut.data_corr = numpy.zeros((len(pairsList),
194 194 len(lagT),len(lagR),
195 195 self.dataIn.nHeights),
196 196 dtype='complex')
197 197
198 198 for l in range(len(pairsList)):
199 199
200 200 firstChannel = pairsList[l][0]
201 201 secondChannel = pairsList[l][1]
202 202
203 203 tmp = None
204 204 tmp = numpy.zeros((len(lagT),len(lagR),self.dataIn.nHeights),dtype='complex')
205 205
206 206 for t in range(len(lagT)):
207 207
208 208 for r in range(len(lagR)):
209 209
210 210 idxT = lagT[t]
211 211 idxR = lagR[r]
212 212
213 213 if idxT >= 0:
214 214 vStacked = numpy.vstack((self.buffer[secondChannel,idxT:,:],
215 215 numpy.zeros((idxT,self.dataIn.nHeights),dtype='complex')))
216 216 else:
217 217 vStacked = numpy.vstack((numpy.zeros((-idxT,self.dataIn.nHeights),dtype='complex'),
218 218 self.buffer[secondChannel,:(nPoints + idxT),:]))
219 219
220 220 if idxR >= 0:
221 221 hStacked = numpy.hstack((vStacked[:,idxR:],numpy.zeros((nPoints,idxR),dtype='complex')))
222 222 else:
223 223 hStacked = numpy.hstack((numpy.zeros((nPoints,-idxR),dtype='complex'),vStacked[:,(self.dataOut.nHeights + idxR)]))
224 224
225 225
226 226 tmp[t,r,:] = numpy.sum((numpy.conjugate(self.buffer[firstChannel,:,:])*hStacked),axis=0)
227 227
228 228
229 229 hStacked = None
230 230 vStacked = None
231 231
232 232 self.dataOut.data_corr[l,:,:,:] = tmp[:,:,:]
233 233
234 234 #Se Calcula los factores de Normalizacion
235 235 self.dataOut.pairsAutoCorr = self.dataOut.getPairsAutoCorr()
236 236 self.dataOut.lagT = lagT*self.dataIn.ippSeconds*self.dataIn.nCohInt
237 237 self.dataOut.lagR = lagR
238 238
239 239 self.calculateNormFactor()
240 240
241 241 self.dataOut.flagNoData = False
242 242 self.buffer = None
243 243 self.firstdatatime = None
244 244 self.profIndex = 0
245 245
246 246 return No newline at end of file
Show file before | Show file after | Hide comments
@@ -1,2144 +1,2154
1 1 import numpy
2 2 import math
3 3 from scipy import optimize
4 4 from scipy import interpolate
5 5 from scipy import signal
6 6 from scipy import stats
7 7 import re
8 8 import datetime
9 9 import copy
10 10 import sys
11 11 import importlib
12 12 import itertools
13 13
14 14 from jroproc_base import ProcessingUnit, Operation
15 15 from schainpy.model.data.jrodata import Parameters
16 16
17 17
18 18 class ParametersProc(ProcessingUnit):
19 19
20 20 nSeconds = None
21 21
22 22 def __init__(self):
23 23 ProcessingUnit.__init__(self)
24 24
25 25 # self.objectDict = {}
26 26 self.buffer = None
27 27 self.firstdatatime = None
28 28 self.profIndex = 0
29 29 self.dataOut = Parameters()
30 30
31 31 def __updateObjFromInput(self):
32 32
33 33 self.dataOut.inputUnit = self.dataIn.type
34 34
35 35 self.dataOut.timeZone = self.dataIn.timeZone
36 36 self.dataOut.dstFlag = self.dataIn.dstFlag
37 37 self.dataOut.errorCount = self.dataIn.errorCount
38 38 self.dataOut.useLocalTime = self.dataIn.useLocalTime
39 39
40 40 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
41 41 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
42 42 self.dataOut.channelList = self.dataIn.channelList
43 43 self.dataOut.heightList = self.dataIn.heightList
44 44 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
45 45 # self.dataOut.nHeights = self.dataIn.nHeights
46 46 # self.dataOut.nChannels = self.dataIn.nChannels
47 47 self.dataOut.nBaud = self.dataIn.nBaud
48 48 self.dataOut.nCode = self.dataIn.nCode
49 49 self.dataOut.code = self.dataIn.code
50 50 # self.dataOut.nProfiles = self.dataOut.nFFTPoints
51 51 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
52 52 self.dataOut.utctime = self.firstdatatime
53 53 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
54 54 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
55 55 # self.dataOut.nCohInt = self.dataIn.nCohInt
56 56 # self.dataOut.nIncohInt = 1
57 57 self.dataOut.ippSeconds = self.dataIn.ippSeconds
58 58 # self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
59 59 self.dataOut.timeInterval = self.dataIn.timeInterval
60 60 self.dataOut.heightList = self.dataIn.getHeiRange()
61 61 self.dataOut.frequency = self.dataIn.frequency
62 62
63 def run(self, nSeconds = None, nProfiles = None):
63 def run(self, nSeconds = 100, nProfiles = None):
64 64
65 65
66 66
67 67 if self.firstdatatime == None:
68 68 self.firstdatatime = self.dataIn.utctime
69 69
70 70 #---------------------- Voltage Data ---------------------------
71 71
72 72 if self.dataIn.type == "Voltage":
73 73 self.dataOut.flagNoData = True
74 if nSeconds != None:
74
75
76 if self.buffer == None:
75 77 self.nSeconds = nSeconds
76 78 self.nProfiles= int(numpy.floor(nSeconds/(self.dataIn.ippSeconds*self.dataIn.nCohInt)))
79
80 self.buffer = numpy.zeros((self.dataIn.nChannels,
81 self.nProfiles,
82 self.dataIn.nHeights),
83 dtype='complex')
77 84
78 if self.buffer is None:
79 self.buffer = numpy.zeros((self.dataIn.nChannels,
80 self.nProfiles,
81 self.dataIn.nHeights),
82 dtype='complex')
83
85 if self.profIndex == 7990:
86 a = 1
87
84 88 self.buffer[:,self.profIndex,:] = self.dataIn.data.copy()
85 89 self.profIndex += 1
86 90
87 91 if self.profIndex == self.nProfiles:
88 92
89 93 self.__updateObjFromInput()
90 94 self.dataOut.data_pre = self.buffer.copy()
91 95 self.dataOut.paramInterval = nSeconds
92 96 self.dataOut.flagNoData = False
93 97
94 98 self.buffer = None
95 99 self.firstdatatime = None
96 100 self.profIndex = 0
97 101 return
98 102
99 103 #---------------------- Spectra Data ---------------------------
100 104
101 105 if self.dataIn.type == "Spectra":
102 106 self.dataOut.data_pre = self.dataIn.data_spc.copy()
103 107 self.dataOut.abscissaList = self.dataIn.getVelRange(1)
104 108 self.dataOut.noise = self.dataIn.getNoise()
105 109 self.dataOut.normFactor = self.dataIn.normFactor
106 110 self.dataOut.groupList = self.dataIn.pairsList
107 111 self.dataOut.flagNoData = False
108 112
109 113 #---------------------- Correlation Data ---------------------------
110 114
111 115 if self.dataIn.type == "Correlation":
112 116 lagRRange = self.dataIn.lagR
113 117 indR = numpy.where(lagRRange == 0)[0][0]
114 118
115 119 self.dataOut.data_pre = self.dataIn.data_corr.copy()[:,:,indR,:]
116 120 self.dataOut.abscissaList = self.dataIn.getLagTRange(1)
117 121 self.dataOut.noise = self.dataIn.noise
118 122 self.dataOut.normFactor = self.dataIn.normFactor
119 123 self.dataOut.data_SNR = self.dataIn.SNR
120 124 self.dataOut.groupList = self.dataIn.pairsList
121 125 self.dataOut.flagNoData = False
122 126
123 127 #---------------------- Correlation Data ---------------------------
124 128
125 129 if self.dataIn.type == "Parameters":
126 130 self.dataOut.copy(self.dataIn)
127 131 self.dataOut.flagNoData = False
128 132
129 133 return True
130 134
131 135 self.__updateObjFromInput()
132 136 self.firstdatatime = None
133 137 self.dataOut.utctimeInit = self.dataIn.utctime
134 self.dataOut.outputInterval = self.dataIn.timeInterval
138 self.dataOut.paramInterval = self.dataIn.timeInterval
135 139
136 140 #------------------- Get Moments ----------------------------------
137 141 def GetMoments(self, channelList = None):
138 142 '''
139 143 Function GetMoments()
140 144
141 145 Input:
142 146 channelList : simple channel list to select e.g. [2,3,7]
143 147 self.dataOut.data_pre
144 148 self.dataOut.abscissaList
145 149 self.dataOut.noise
146 150
147 151 Affected:
148 152 self.dataOut.data_param
149 153 self.dataOut.data_SNR
150 154
151 155 '''
152 156 data = self.dataOut.data_pre
153 157 absc = self.dataOut.abscissaList[:-1]
154 158 noise = self.dataOut.noise
155 159
156 160 data_param = numpy.zeros((data.shape[0], 4, data.shape[2]))
157 161
158 162 if channelList== None:
159 163 channelList = self.dataIn.channelList
160 164 self.dataOut.channelList = channelList
161 165
162 166 for ind in channelList:
163 167 data_param[ind,:,:] = self.__calculateMoments(data[ind,:,:], absc, noise[ind])
164 168
165 169 self.dataOut.data_param = data_param[:,1:,:]
166 170 self.dataOut.data_SNR = data_param[:,0]
167 171 return
168 172
169 173 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):
170 174
171 175 if (nicoh == None): nicoh = 1
172 176 if (graph == None): graph = 0
173 177 if (smooth == None): smooth = 0
174 178 elif (self.smooth < 3): smooth = 0
175 179
176 180 if (type1 == None): type1 = 0
177 if (fwindow is None): fwindow = numpy.zeros(oldfreq.size) + 1
181 if (fwindow == None): fwindow = numpy.zeros(oldfreq.size) + 1
178 182 if (snrth == None): snrth = -3
179 183 if (dc == None): dc = 0
180 184 if (aliasing == None): aliasing = 0
181 185 if (oldfd == None): oldfd = 0
182 186 if (wwauto == None): wwauto = 0
183 187
184 188 if (n0 < 1.e-20): n0 = 1.e-20
185 189
186 190 freq = oldfreq
187 191 vec_power = numpy.zeros(oldspec.shape[1])
188 192 vec_fd = numpy.zeros(oldspec.shape[1])
189 193 vec_w = numpy.zeros(oldspec.shape[1])
190 194 vec_snr = numpy.zeros(oldspec.shape[1])
191 195
192 196 for ind in range(oldspec.shape[1]):
193 197
194 198 spec = oldspec[:,ind]
195 199 aux = spec*fwindow
196 200 max_spec = aux.max()
197 201 m = list(aux).index(max_spec)
198 202
199 203 #Smooth
200 204 if (smooth == 0): spec2 = spec
201 205 else: spec2 = scipy.ndimage.filters.uniform_filter1d(spec,size=smooth)
202 206
203 207 # Calculo de Momentos
204 208 bb = spec2[range(m,spec2.size)]
205 209 bb = (bb<n0).nonzero()
206 210 bb = bb[0]
207 211
208 212 ss = spec2[range(0,m + 1)]
209 213 ss = (ss<n0).nonzero()
210 214 ss = ss[0]
211 215
212 216 if (bb.size == 0):
213 217 bb0 = spec.size - 1 - m
214 218 else:
215 219 bb0 = bb[0] - 1
216 220 if (bb0 < 0):
217 221 bb0 = 0
218 222
219 223 if (ss.size == 0): ss1 = 1
220 224 else: ss1 = max(ss) + 1
221 225
222 226 if (ss1 > m): ss1 = m
223 227
224 228 valid = numpy.asarray(range(int(m + bb0 - ss1 + 1))) + ss1
225 229 power = ((spec2[valid] - n0)*fwindow[valid]).sum()
226 230 fd = ((spec2[valid]- n0)*freq[valid]*fwindow[valid]).sum()/power
227 231 w = math.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum()/power)
228 232 snr = (spec2.mean()-n0)/n0
229 233
230 234 if (snr < 1.e-20) :
231 235 snr = 1.e-20
232 236
233 237 vec_power[ind] = power
234 238 vec_fd[ind] = fd
235 239 vec_w[ind] = w
236 240 vec_snr[ind] = snr
237 241
238 242 moments = numpy.vstack((vec_snr, vec_power, vec_fd, vec_w))
239 243 return moments
240 244
241 245 #------------------ Get SA Parameters --------------------------
242 246
243 247 def GetSAParameters(self):
244 248 pairslist = self.dataOut.groupList
245 249 num_pairs = len(pairslist)
246 250
247 251 vel = self.dataOut.abscissaList
248 252 spectra = self.dataOut.data_pre
249 253 cspectra = self.dataIn.data_cspc
250 254 delta_v = vel[1] - vel[0]
251 255
252 256 #Calculating the power spectrum
253 257 spc_pow = numpy.sum(spectra, 3)*delta_v
254 258 #Normalizing Spectra
255 259 norm_spectra = spectra/spc_pow
256 260 #Calculating the norm_spectra at peak
257 261 max_spectra = numpy.max(norm_spectra, 3)
258 262
259 263 #Normalizing Cross Spectra
260 264 norm_cspectra = numpy.zeros(cspectra.shape)
261 265
262 266 for i in range(num_chan):
263 267 norm_cspectra[i,:,:] = cspectra[i,:,:]/numpy.sqrt(spc_pow[pairslist[i][0],:]*spc_pow[pairslist[i][1],:])
264 268
265 269 max_cspectra = numpy.max(norm_cspectra,2)
266 270 max_cspectra_index = numpy.argmax(norm_cspectra, 2)
267 271
268 272 for i in range(num_pairs):
269 273 cspc_par[i,:,:] = __calculateMoments(norm_cspectra)
270 274 #------------------- Get Lags ----------------------------------
271 275
272 276 def GetLags(self):
273 277 '''
274 278 Function GetMoments()
275 279
276 280 Input:
277 281 self.dataOut.data_pre
278 282 self.dataOut.abscissaList
279 283 self.dataOut.noise
280 284 self.dataOut.normFactor
281 285 self.dataOut.data_SNR
282 286 self.dataOut.groupList
283 287 self.dataOut.nChannels
284 288
285 289 Affected:
286 290 self.dataOut.data_param
287 291
288 292 '''
289 293
290 294 data = self.dataOut.data_pre
291 295 normFactor = self.dataOut.normFactor
292 296 nHeights = self.dataOut.nHeights
293 297 absc = self.dataOut.abscissaList[:-1]
294 298 noise = self.dataOut.noise
295 299 SNR = self.dataOut.data_SNR
296 300 pairsList = self.dataOut.groupList
297 301 nChannels = self.dataOut.nChannels
298 302 pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)
299 303 self.dataOut.data_param = numpy.zeros((len(pairsCrossCorr)*2 + 1, nHeights))
300 304
301 305 dataNorm = numpy.abs(data)
302 306 for l in range(len(pairsList)):
303 307 dataNorm[l,:,:] = dataNorm[l,:,:]/normFactor[l,:]
304 308
305 309 self.dataOut.data_param[:-1,:] = self.__calculateTaus(dataNorm, pairsCrossCorr, pairsAutoCorr, absc)
306 310 self.dataOut.data_param[-1,:] = self.__calculateLag1Phase(data, pairsAutoCorr, absc)
307 311 return
308 312
309 313 def __getPairsAutoCorr(self, pairsList, nChannels):
310 314
311 315 pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
312 316
313 317 for l in range(len(pairsList)):
314 318 firstChannel = pairsList[l][0]
315 319 secondChannel = pairsList[l][1]
316 320
317 321 #Obteniendo pares de Autocorrelacion
318 322 if firstChannel == secondChannel:
319 323 pairsAutoCorr[firstChannel] = int(l)
320 324
321 325 pairsAutoCorr = pairsAutoCorr.astype(int)
322 326
323 327 pairsCrossCorr = range(len(pairsList))
324 328 pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
325 329
326 330 return pairsAutoCorr, pairsCrossCorr
327 331
328 332 def __calculateTaus(self, data, pairsCrossCorr, pairsAutoCorr, lagTRange):
329 333
330 334 Pt0 = data.shape[1]/2
331 335 #Funcion de Autocorrelacion
332 336 dataAutoCorr = stats.nanmean(data[pairsAutoCorr,:,:], axis = 0)
333 337
334 338 #Obtencion Indice de TauCross
335 339 indCross = data[pairsCrossCorr,:,:].argmax(axis = 1)
336 340 #Obtencion Indice de TauAuto
337 341 indAuto = numpy.zeros(indCross.shape,dtype = 'int')
338 342 CCValue = data[pairsCrossCorr,Pt0,:]
339 343 for i in range(pairsCrossCorr.size):
340 344 indAuto[i,:] = numpy.abs(dataAutoCorr - CCValue[i,:]).argmin(axis = 0)
341 345
342 346 #Obtencion de TauCross y TauAuto
343 347 tauCross = lagTRange[indCross]
344 348 tauAuto = lagTRange[indAuto]
345 349
346 350 Nan1, Nan2 = numpy.where(tauCross == lagTRange[0])
347 351
348 352 tauCross[Nan1,Nan2] = numpy.nan
349 353 tauAuto[Nan1,Nan2] = numpy.nan
350 354 tau = numpy.vstack((tauCross,tauAuto))
351 355
352 356 return tau
353 357
354 358 def __calculateLag1Phase(self, data, pairs, lagTRange):
355 359 data1 = stats.nanmean(data[pairs,:,:], axis = 0)
356 360 lag1 = numpy.where(lagTRange == 0)[0][0] + 1
357 361
358 362 phase = numpy.angle(data1[lag1,:])
359 363
360 364 return phase
361 365 #------------------- Detect Meteors ------------------------------
362 366
363 367 def MeteorDetection(self, hei_ref = None, tauindex = 0,
364 368 predefinedPhaseShifts = None, centerReceiverIndex = 2, saveAll = False,
365 369 cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25,
366 370 noise_timeStep = 4, noise_multiple = 4,
367 371 multDet_timeLimit = 1, multDet_rangeLimit = 3,
368 372 phaseThresh = 20, SNRThresh = 8,
369 373 hmin = 70, hmax=110, azimuth = 0) :
370 374
371 375 '''
372 376 Function DetectMeteors()
373 377 Project developed with paper:
374 378 HOLDSWORTH ET AL. 2004
375 379
376 380 Input:
377 381 self.dataOut.data_pre
378 382
379 383 centerReceiverIndex: From the channels, which is the center receiver
380 384
381 385 hei_ref: Height reference for the Beacon signal extraction
382 386 tauindex:
383 387 predefinedPhaseShifts: Predefined phase offset for the voltge signals
384 388
385 389 cohDetection: Whether to user Coherent detection or not
386 390 cohDet_timeStep: Coherent Detection calculation time step
387 391 cohDet_thresh: Coherent Detection phase threshold to correct phases
388 392
389 393 noise_timeStep: Noise calculation time step
390 394 noise_multiple: Noise multiple to define signal threshold
391 395
392 396 multDet_timeLimit: Multiple Detection Removal time limit in seconds
393 397 multDet_rangeLimit: Multiple Detection Removal range limit in km
394 398
395 399 phaseThresh: Maximum phase difference between receiver to be consider a meteor
396 400 SNRThresh: Minimum SNR threshold of the meteor signal to be consider a meteor
397 401
398 402 hmin: Minimum Height of the meteor to use it in the further wind estimations
399 403 hmax: Maximum Height of the meteor to use it in the further wind estimations
400 404 azimuth: Azimuth angle correction
401 405
402 406 Affected:
403 407 self.dataOut.data_param
404 408
405 409 Rejection Criteria (Errors):
406 410 0: No error; analysis OK
407 411 1: SNR < SNR threshold
408 412 2: angle of arrival (AOA) ambiguously determined
409 413 3: AOA estimate not feasible
410 414 4: Large difference in AOAs obtained from different antenna baselines
411 415 5: echo at start or end of time series
412 416 6: echo less than 5 examples long; too short for analysis
413 417 7: echo rise exceeds 0.3s
414 418 8: echo decay time less than twice rise time
415 419 9: large power level before echo
416 420 10: large power level after echo
417 421 11: poor fit to amplitude for estimation of decay time
418 422 12: poor fit to CCF phase variation for estimation of radial drift velocity
419 423 13: height unresolvable echo: not valid height within 70 to 110 km
420 424 14: height ambiguous echo: more then one possible height within 70 to 110 km
421 425 15: radial drift velocity or projected horizontal velocity exceeds 200 m/s
422 426 16: oscilatory echo, indicating event most likely not an underdense echo
423 427
424 428 17: phase difference in meteor Reestimation
425 429
426 430 Data Storage:
427 431 Meteors for Wind Estimation (8):
428 432 Day Hour | Range Height
429 433 Azimuth Zenith errorCosDir
430 434 VelRad errorVelRad
431 435 TypeError
432 436
433 437 '''
434 438 #Get Beacon signal
435 439 newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
436 440
437 441 if hei_ref != None:
438 442 newheis = numpy.where(self.dataOut.heightList>hei_ref)
439 443
440 444 heiRang = self.dataOut.getHeiRange()
441 445 #Pairs List
442 446 pairslist = []
443 447 nChannel = self.dataOut.nChannels
444 448 for i in range(nChannel):
445 449 if i != centerReceiverIndex:
446 450 pairslist.append((centerReceiverIndex,i))
447 451
448 452 #****************REMOVING HARDWARE PHASE DIFFERENCES***************
449 453 # see if the user put in pre defined phase shifts
450 454 voltsPShift = self.dataOut.data_pre.copy()
451 455
452 456 if predefinedPhaseShifts != None:
453 457 hardwarePhaseShifts = numpy.array(predefinedPhaseShifts)*numpy.pi/180
454 458
455 459 elif beaconPhaseShifts:
456 460 #get hardware phase shifts using beacon signal
457 461 hardwarePhaseShifts = self.__getHardwarePhaseDiff(self.dataOut.data_pre, pairslist, newheis, 10)
458 462 hardwarePhaseShifts = numpy.insert(hardwarePhaseShifts,centerReceiverIndex,0)
459 463
460 464 else:
461 465 hardwarePhaseShifts = numpy.zeros(5)
462 466
463 467
464 468 voltsPShift = numpy.zeros((self.dataOut.data_pre.shape[0],self.dataOut.data_pre.shape[1],self.dataOut.data_pre.shape[2]), dtype = 'complex')
465 469 for i in range(self.dataOut.data_pre.shape[0]):
466 470 voltsPShift[i,:,:] = self.__shiftPhase(self.dataOut.data_pre[i,:,:], hardwarePhaseShifts[i])
467 471
468 472
469 473 #******************END OF REMOVING HARDWARE PHASE DIFFERENCES*********
470 474
471 475 #Remove DC
472 476 voltsDC = numpy.mean(voltsPShift,1)
473 477 voltsDC = numpy.mean(voltsDC,1)
474 478 for i in range(voltsDC.shape[0]):
475 479 voltsPShift[i] = voltsPShift[i] - voltsDC[i]
476 480
477 481 #Don't considerate last heights, theyre used to calculate Hardware Phase Shift
478 482 voltsPShift = voltsPShift[:,:,:newheis[0][0]]
479 483
480 484 #************ FIND POWER OF DATA W/COH OR NON COH DETECTION (3.4) **********
481 485 #Coherent Detection
482 486 if cohDetection:
483 487 #use coherent detection to get the net power
484 488 cohDet_thresh = cohDet_thresh*numpy.pi/180
485 489 voltsPShift = self.__coherentDetection(voltsPShift, cohDet_timeStep, self.dataOut.timeInterval, pairslist, cohDet_thresh)
486 490
487 491 #Non-coherent detection!
488 492 powerNet = numpy.nansum(numpy.abs(voltsPShift[:,:,:])**2,0)
489 493 #********** END OF COH/NON-COH POWER CALCULATION**********************
490 494
491 495 #********** FIND THE NOISE LEVEL AND POSSIBLE METEORS ****************
492 496 #Get noise
493 497 noise, noise1 = self.__getNoise(powerNet, noise_timeStep, self.dataOut.timeInterval)
494 498 # noise = self.getNoise1(powerNet, noise_timeStep, self.dataOut.timeInterval)
495 499 #Get signal threshold
496 500 signalThresh = noise_multiple*noise
497 501 #Meteor echoes detection
498 502 listMeteors = self.__findMeteors(powerNet, signalThresh)
499 503 #******* END OF NOISE LEVEL AND POSSIBLE METEORS CACULATION **********
500 504
501 505 #************** REMOVE MULTIPLE DETECTIONS (3.5) ***************************
502 506 #Parameters
503 507 heiRange = self.dataOut.getHeiRange()
504 508 rangeInterval = heiRange[1] - heiRange[0]
505 509 rangeLimit = multDet_rangeLimit/rangeInterval
506 510 timeLimit = multDet_timeLimit/self.dataOut.timeInterval
507 511 #Multiple detection removals
508 512 listMeteors1 = self.__removeMultipleDetections(listMeteors, rangeLimit, timeLimit)
509 513 #************ END OF REMOVE MULTIPLE DETECTIONS **********************
510 514
511 515 #********************* METEOR REESTIMATION (3.7, 3.8, 3.9, 3.10) ********************
512 516 #Parameters
513 517 phaseThresh = phaseThresh*numpy.pi/180
514 518 thresh = [phaseThresh, noise_multiple, SNRThresh]
515 519 #Meteor reestimation (Errors N 1, 6, 12, 17)
516 520 listMeteors2, listMeteorsPower, listMeteorsVolts = self.__meteorReestimation(listMeteors1, voltsPShift, pairslist, thresh, noise, self.dataOut.timeInterval, self.dataOut.frequency)
517 521 # listMeteors2, listMeteorsPower, listMeteorsVolts = self.meteorReestimation3(listMeteors2, listMeteorsPower, listMeteorsVolts, voltsPShift, pairslist, thresh, noise)
518 522 #Estimation of decay times (Errors N 7, 8, 11)
519 523 listMeteors3 = self.__estimateDecayTime(listMeteors2, listMeteorsPower, self.dataOut.timeInterval, self.dataOut.frequency)
520 524 #******************* END OF METEOR REESTIMATION *******************
521 525
522 526 #********************* METEOR PARAMETERS CALCULATION (3.11, 3.12, 3.13) **************************
523 527 #Calculating Radial Velocity (Error N 15)
524 528 radialStdThresh = 10
525 529 listMeteors4 = self.__getRadialVelocity(listMeteors3, listMeteorsVolts, radialStdThresh, pairslist, self.dataOut.timeInterval)
526 530
527 531 if len(listMeteors4) > 0:
528 532
529 533 pairsList = []
530 534 pairx = (0,3)
531 535 pairy = (1,2)
532 536 pairsList.append(pairx)
533 537 pairsList.append(pairy)
534 538
535 539 #Setting New Array
536 540 date = repr(self.dataOut.datatime)
537 541 arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang)
538 542
539 543 meteorOps = MeteorOperations()
540 544 jph = numpy.array([0,0,0,0])
541 545 h = (hmin,hmax)
542 546 arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, jph)
543 547
544 548 # #Calculate AOA (Error N 3, 4)
545 549 # #JONES ET AL. 1998
546 550 # error = arrayParameters[:,-1]
547 551 # AOAthresh = numpy.pi/8
548 552 # phases = -arrayParameters[:,9:13]
549 553 # arrayParameters[:,4:7], arrayParameters[:,-1] = meteorOps.getAOA(phases, pairsList, error, AOAthresh, azimuth)
550 554 #
551 555 # #Calculate Heights (Error N 13 and 14)
552 556 # error = arrayParameters[:,-1]
553 557 # Ranges = arrayParameters[:,2]
554 558 # zenith = arrayParameters[:,5]
555 559 # arrayParameters[:,3], arrayParameters[:,-1] = meteorOps.getHeights(Ranges, zenith, error, hmin, hmax)
556 560 # error = arrayParameters[:,-1]
557 561 #********************* END OF PARAMETERS CALCULATION **************************
558 562
559 563 #***************************+ PASS DATA TO NEXT STEP **********************
560 564 arrayFinal = arrayParameters.reshape((1,arrayParameters.shape[0],arrayParameters.shape[1]))
561 565 self.dataOut.data_param = arrayFinal
562
566
567 if arrayFinal == None:
568 self.dataOut.flagNoData = True
569
563 570 return
564 571
565 572 def __getHardwarePhaseDiff(self, voltage0, pairslist, newheis, n):
566 573
567 574 minIndex = min(newheis[0])
568 575 maxIndex = max(newheis[0])
569 576
570 577 voltage = voltage0[:,:,minIndex:maxIndex+1]
571 578 nLength = voltage.shape[1]/n
572 579 nMin = 0
573 580 nMax = 0
574 581 phaseOffset = numpy.zeros((len(pairslist),n))
575 582
576 583 for i in range(n):
577 584 nMax += nLength
578 585 phaseCCF = -numpy.angle(self.__calculateCCF(voltage[:,nMin:nMax,:], pairslist, [0]))
579 586 phaseCCF = numpy.mean(phaseCCF, axis = 2)
580 587 phaseOffset[:,i] = phaseCCF.transpose()
581 588 nMin = nMax
582 589 # phaseDiff, phaseArrival = self.estimatePhaseDifference(voltage, pairslist)
583 590
584 591 #Remove Outliers
585 592 factor = 2
586 593 wt = phaseOffset - signal.medfilt(phaseOffset,(1,5))
587 594 dw = numpy.std(wt,axis = 1)
588 595 dw = dw.reshape((dw.size,1))
589 596 ind = numpy.where(numpy.logical_or(wt>dw*factor,wt<-dw*factor))
590 597 phaseOffset[ind] = numpy.nan
591 598 phaseOffset = stats.nanmean(phaseOffset, axis=1)
592 599
593 600 return phaseOffset
594 601
595 602 def __shiftPhase(self, data, phaseShift):
596 603 #this will shift the phase of a complex number
597 604 dataShifted = numpy.abs(data) * numpy.exp((numpy.angle(data)+phaseShift)*1j)
598 605 return dataShifted
599 606
600 607 def __estimatePhaseDifference(self, array, pairslist):
601 608 nChannel = array.shape[0]
602 609 nHeights = array.shape[2]
603 610 numPairs = len(pairslist)
604 611 # phaseCCF = numpy.zeros((nChannel, 5, nHeights))
605 612 phaseCCF = numpy.angle(self.__calculateCCF(array, pairslist, [-2,-1,0,1,2]))
606 613
607 614 #Correct phases
608 615 derPhaseCCF = phaseCCF[:,1:,:] - phaseCCF[:,0:-1,:]
609 616 indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
610 617
611 618 if indDer[0].shape[0] > 0:
612 619 for i in range(indDer[0].shape[0]):
613 620 signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i],indDer[2][i]])
614 621 phaseCCF[indDer[0][i],indDer[1][i]+1:,:] += signo*2*numpy.pi
615 622
616 623 # for j in range(numSides):
617 624 # phaseCCFAux = self.calculateCCF(arrayCenter, arraySides[j,:,:], [-2,1,0,1,2])
618 625 # phaseCCF[j,:,:] = numpy.angle(phaseCCFAux)
619 626 #
620 627 #Linear
621 628 phaseInt = numpy.zeros((numPairs,1))
622 629 angAllCCF = phaseCCF[:,[0,1,3,4],0]
623 630 for j in range(numPairs):
624 631 fit = stats.linregress([-2,-1,1,2],angAllCCF[j,:])
625 632 phaseInt[j] = fit[1]
626 633 #Phase Differences
627 634 phaseDiff = phaseInt - phaseCCF[:,2,:]
628 635 phaseArrival = phaseInt.reshape(phaseInt.size)
629 636
630 637 #Dealias
631 638 indAlias = numpy.where(phaseArrival > numpy.pi)
632 639 phaseArrival[indAlias] -= 2*numpy.pi
633 640 indAlias = numpy.where(phaseArrival < -numpy.pi)
634 641 phaseArrival[indAlias] += 2*numpy.pi
635 642
636 643 return phaseDiff, phaseArrival
637 644
638 645 def __coherentDetection(self, volts, timeSegment, timeInterval, pairslist, thresh):
639 646 #this function will run the coherent detection used in Holdworth et al. 2004 and return the net power
640 647 #find the phase shifts of each channel over 1 second intervals
641 648 #only look at ranges below the beacon signal
642 649 numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
643 650 numBlocks = int(volts.shape[1]/numProfPerBlock)
644 651 numHeights = volts.shape[2]
645 652 nChannel = volts.shape[0]
646 653 voltsCohDet = volts.copy()
647 654
648 655 pairsarray = numpy.array(pairslist)
649 656 indSides = pairsarray[:,1]
650 657 # indSides = numpy.array(range(nChannel))
651 658 # indSides = numpy.delete(indSides, indCenter)
652 659 #
653 660 # listCenter = numpy.array_split(volts[indCenter,:,:], numBlocks, 0)
654 661 listBlocks = numpy.array_split(volts, numBlocks, 1)
655 662
656 663 startInd = 0
657 664 endInd = 0
658 665
659 666 for i in range(numBlocks):
660 667 startInd = endInd
661 668 endInd = endInd + listBlocks[i].shape[1]
662 669
663 670 arrayBlock = listBlocks[i]
664 671 # arrayBlockCenter = listCenter[i]
665 672
666 673 #Estimate the Phase Difference
667 674 phaseDiff, aux = self.__estimatePhaseDifference(arrayBlock, pairslist)
668 675 #Phase Difference RMS
669 676 arrayPhaseRMS = numpy.abs(phaseDiff)
670 677 phaseRMSaux = numpy.sum(arrayPhaseRMS < thresh,0)
671 678 indPhase = numpy.where(phaseRMSaux==4)
672 679 #Shifting
673 680 if indPhase[0].shape[0] > 0:
674 681 for j in range(indSides.size):
675 682 arrayBlock[indSides[j],:,indPhase] = self.__shiftPhase(arrayBlock[indSides[j],:,indPhase], phaseDiff[j,indPhase].transpose())
676 683 voltsCohDet[:,startInd:endInd,:] = arrayBlock
677 684
678 685 return voltsCohDet
679 686
680 687 def __calculateCCF(self, volts, pairslist ,laglist):
681 688
682 689 nHeights = volts.shape[2]
683 690 nPoints = volts.shape[1]
684 691 voltsCCF = numpy.zeros((len(pairslist), len(laglist), nHeights),dtype = 'complex')
685 692
686 693 for i in range(len(pairslist)):
687 694 volts1 = volts[pairslist[i][0]]
688 695 volts2 = volts[pairslist[i][1]]
689 696
690 697 for t in range(len(laglist)):
691 698 idxT = laglist[t]
692 699 if idxT >= 0:
693 700 vStacked = numpy.vstack((volts2[idxT:,:],
694 701 numpy.zeros((idxT, nHeights),dtype='complex')))
695 702 else:
696 703 vStacked = numpy.vstack((numpy.zeros((-idxT, nHeights),dtype='complex'),
697 704 volts2[:(nPoints + idxT),:]))
698 705 voltsCCF[i,t,:] = numpy.sum((numpy.conjugate(volts1)*vStacked),axis=0)
699 706
700 707 vStacked = None
701 708 return voltsCCF
702 709
703 710 def __getNoise(self, power, timeSegment, timeInterval):
704 711 numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
705 712 numBlocks = int(power.shape[0]/numProfPerBlock)
706 713 numHeights = power.shape[1]
707 714
708 715 listPower = numpy.array_split(power, numBlocks, 0)
709 716 noise = numpy.zeros((power.shape[0], power.shape[1]))
710 717 noise1 = numpy.zeros((power.shape[0], power.shape[1]))
711 718
712 719 startInd = 0
713 720 endInd = 0
714 721
715 722 for i in range(numBlocks): #split por canal
716 723 startInd = endInd
717 724 endInd = endInd + listPower[i].shape[0]
718 725
719 726 arrayBlock = listPower[i]
720 727 noiseAux = numpy.mean(arrayBlock, 0)
721 728 # noiseAux = numpy.median(noiseAux)
722 729 # noiseAux = numpy.mean(arrayBlock)
723 730 noise[startInd:endInd,:] = noise[startInd:endInd,:] + noiseAux
724 731
725 732 noiseAux1 = numpy.mean(arrayBlock)
726 733 noise1[startInd:endInd,:] = noise1[startInd:endInd,:] + noiseAux1
727 734
728 735 return noise, noise1
729 736
730 737 def __findMeteors(self, power, thresh):
731 738 nProf = power.shape[0]
732 739 nHeights = power.shape[1]
733 740 listMeteors = []
734 741
735 742 for i in range(nHeights):
736 743 powerAux = power[:,i]
737 744 threshAux = thresh[:,i]
738 745
739 746 indUPthresh = numpy.where(powerAux > threshAux)[0]
740 747 indDNthresh = numpy.where(powerAux <= threshAux)[0]
741 748
742 749 j = 0
743 750
744 751 while (j < indUPthresh.size - 2):
745 752 if (indUPthresh[j + 2] == indUPthresh[j] + 2):
746 753 indDNAux = numpy.where(indDNthresh > indUPthresh[j])
747 754 indDNthresh = indDNthresh[indDNAux]
748 755
749 756 if (indDNthresh.size > 0):
750 757 indEnd = indDNthresh[0] - 1
751 758 indInit = indUPthresh[j]
752 759
753 760 meteor = powerAux[indInit:indEnd + 1]
754 761 indPeak = meteor.argmax() + indInit
755 762 FLA = sum(numpy.conj(meteor)*numpy.hstack((meteor[1:],0)))
756 763
757 764 listMeteors.append(numpy.array([i,indInit,indPeak,indEnd,FLA])) #CHEQUEAR!!!!!
758 765 j = numpy.where(indUPthresh == indEnd)[0] + 1
759 766 else: j+=1
760 767 else: j+=1
761 768
762 769 return listMeteors
763 770
764 771 def __removeMultipleDetections(self,listMeteors, rangeLimit, timeLimit):
765 772
766 773 arrayMeteors = numpy.asarray(listMeteors)
767 774 listMeteors1 = []
768 775
769 776 while arrayMeteors.shape[0] > 0:
770 777 FLAs = arrayMeteors[:,4]
771 778 maxFLA = FLAs.argmax()
772 779 listMeteors1.append(arrayMeteors[maxFLA,:])
773 780
774 781 MeteorInitTime = arrayMeteors[maxFLA,1]
775 782 MeteorEndTime = arrayMeteors[maxFLA,3]
776 783 MeteorHeight = arrayMeteors[maxFLA,0]
777 784
778 785 #Check neighborhood
779 786 maxHeightIndex = MeteorHeight + rangeLimit
780 787 minHeightIndex = MeteorHeight - rangeLimit
781 788 minTimeIndex = MeteorInitTime - timeLimit
782 789 maxTimeIndex = MeteorEndTime + timeLimit
783 790
784 791 #Check Heights
785 792 indHeight = numpy.logical_and(arrayMeteors[:,0] >= minHeightIndex, arrayMeteors[:,0] <= maxHeightIndex)
786 793 indTime = numpy.logical_and(arrayMeteors[:,3] >= minTimeIndex, arrayMeteors[:,1] <= maxTimeIndex)
787 794 indBoth = numpy.where(numpy.logical_and(indTime,indHeight))
788 795
789 796 arrayMeteors = numpy.delete(arrayMeteors, indBoth, axis = 0)
790 797
791 798 return listMeteors1
792 799
793 800 def __meteorReestimation(self, listMeteors, volts, pairslist, thresh, noise, timeInterval,frequency):
794 801 numHeights = volts.shape[2]
795 802 nChannel = volts.shape[0]
796 803
797 804 thresholdPhase = thresh[0]
798 805 thresholdNoise = thresh[1]
799 806 thresholdDB = float(thresh[2])
800 807
801 808 thresholdDB1 = 10**(thresholdDB/10)
802 809 pairsarray = numpy.array(pairslist)
803 810 indSides = pairsarray[:,1]
804 811
805 812 pairslist1 = list(pairslist)
806 813 pairslist1.append((0,1))
807 814 pairslist1.append((3,4))
808 815
809 816 listMeteors1 = []
810 817 listPowerSeries = []
811 818 listVoltageSeries = []
812 819 #volts has the war data
813 820
814 821 if frequency == 30e6:
815 822 timeLag = 45*10**-3
816 823 else:
817 824 timeLag = 15*10**-3
818 825 lag = numpy.ceil(timeLag/timeInterval)
819 826
820 827 for i in range(len(listMeteors)):
821 828
822 829 ###################### 3.6 - 3.7 PARAMETERS REESTIMATION #########################
823 830 meteorAux = numpy.zeros(16)
824 831
825 832 #Loading meteor Data (mHeight, mStart, mPeak, mEnd)
826 833 mHeight = listMeteors[i][0]
827 834 mStart = listMeteors[i][1]
828 835 mPeak = listMeteors[i][2]
829 836 mEnd = listMeteors[i][3]
830 837
831 838 #get the volt data between the start and end times of the meteor
832 839 meteorVolts = volts[:,mStart:mEnd+1,mHeight]
833 840 meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
834 841
835 842 #3.6. Phase Difference estimation
836 843 phaseDiff, aux = self.__estimatePhaseDifference(meteorVolts, pairslist)
837 844
838 845 #3.7. Phase difference removal & meteor start, peak and end times reestimated
839 846 #meteorVolts0.- all Channels, all Profiles
840 847 meteorVolts0 = volts[:,:,mHeight]
841 848 meteorThresh = noise[:,mHeight]*thresholdNoise
842 849 meteorNoise = noise[:,mHeight]
843 850 meteorVolts0[indSides,:] = self.__shiftPhase(meteorVolts0[indSides,:], phaseDiff) #Phase Shifting
844 851 powerNet0 = numpy.nansum(numpy.abs(meteorVolts0)**2, axis = 0) #Power
845 852
846 853 #Times reestimation
847 854 mStart1 = numpy.where(powerNet0[:mPeak] < meteorThresh[:mPeak])[0]
848 855 if mStart1.size > 0:
849 856 mStart1 = mStart1[-1] + 1
850 857
851 858 else:
852 859 mStart1 = mPeak
853 860
854 861 mEnd1 = numpy.where(powerNet0[mPeak:] < meteorThresh[mPeak:])[0][0] + mPeak - 1
855 862 mEndDecayTime1 = numpy.where(powerNet0[mPeak:] < meteorNoise[mPeak:])[0]
856 863 if mEndDecayTime1.size == 0:
857 864 mEndDecayTime1 = powerNet0.size
858 865 else:
859 866 mEndDecayTime1 = mEndDecayTime1[0] + mPeak - 1
860 867 # mPeak1 = meteorVolts0[mStart1:mEnd1 + 1].argmax()
861 868
862 869 #meteorVolts1.- all Channels, from start to end
863 870 meteorVolts1 = meteorVolts0[:,mStart1:mEnd1 + 1]
864 871 meteorVolts2 = meteorVolts0[:,mPeak + lag:mEnd1 + 1]
865 872 if meteorVolts2.shape[1] == 0:
866 873 meteorVolts2 = meteorVolts0[:,mPeak:mEnd1 + 1]
867 874 meteorVolts1 = meteorVolts1.reshape(meteorVolts1.shape[0], meteorVolts1.shape[1], 1)
868 875 meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1], 1)
869 876 ##################### END PARAMETERS REESTIMATION #########################
870 877
871 878 ##################### 3.8 PHASE DIFFERENCE REESTIMATION ########################
872 879 # if mEnd1 - mStart1 > 4: #Error Number 6: echo less than 5 samples long; too short for analysis
873 880 if meteorVolts2.shape[1] > 0:
874 881 #Phase Difference re-estimation
875 882 phaseDiff1, phaseDiffint = self.__estimatePhaseDifference(meteorVolts2, pairslist1) #Phase Difference Estimation
876 883 # phaseDiff1, phaseDiffint = self.estimatePhaseDifference(meteorVolts2, pairslist)
877 884 meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1])
878 885 phaseDiff11 = numpy.reshape(phaseDiff1, (phaseDiff1.shape[0],1))
879 886 meteorVolts2[indSides,:] = self.__shiftPhase(meteorVolts2[indSides,:], phaseDiff11[0:4]) #Phase Shifting
880 887
881 888 #Phase Difference RMS
882 889 phaseRMS1 = numpy.sqrt(numpy.mean(numpy.square(phaseDiff1)))
883 890 powerNet1 = numpy.nansum(numpy.abs(meteorVolts1[:,:])**2,0)
884 891 #Data from Meteor
885 892 mPeak1 = powerNet1.argmax() + mStart1
886 893 mPeakPower1 = powerNet1.max()
887 894 noiseAux = sum(noise[mStart1:mEnd1 + 1,mHeight])
888 895 mSNR1 = (sum(powerNet1)-noiseAux)/noiseAux
889 896 Meteor1 = numpy.array([mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1])
890 897 Meteor1 = numpy.hstack((Meteor1,phaseDiffint))
891 898 PowerSeries = powerNet0[mStart1:mEndDecayTime1 + 1]
892 899 #Vectorize
893 900 meteorAux[0:7] = [mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1]
894 901 meteorAux[7:11] = phaseDiffint[0:4]
895 902
896 903 #Rejection Criterions
897 904 if phaseRMS1 > thresholdPhase: #Error Number 17: Phase variation
898 905 meteorAux[-1] = 17
899 906 elif mSNR1 < thresholdDB1: #Error Number 1: SNR < threshold dB
900 907 meteorAux[-1] = 1
901 908
902 909
903 910 else:
904 911 meteorAux[0:4] = [mHeight, mStart, mPeak, mEnd]
905 912 meteorAux[-1] = 6 #Error Number 6: echo less than 5 samples long; too short for analysis
906 913 PowerSeries = 0
907 914
908 915 listMeteors1.append(meteorAux)
909 916 listPowerSeries.append(PowerSeries)
910 917 listVoltageSeries.append(meteorVolts1)
911 918
912 919 return listMeteors1, listPowerSeries, listVoltageSeries
913 920
914 921 def __estimateDecayTime(self, listMeteors, listPower, timeInterval, frequency):
915 922
916 923 threshError = 10
917 924 #Depending if it is 30 or 50 MHz
918 925 if frequency == 30e6:
919 926 timeLag = 45*10**-3
920 927 else:
921 928 timeLag = 15*10**-3
922 929 lag = numpy.ceil(timeLag/timeInterval)
923 930
924 931 listMeteors1 = []
925 932
926 933 for i in range(len(listMeteors)):
927 934 meteorPower = listPower[i]
928 935 meteorAux = listMeteors[i]
929 936
930 937 if meteorAux[-1] == 0:
931 938
932 939 try:
933 940 indmax = meteorPower.argmax()
934 941 indlag = indmax + lag
935 942
936 943 y = meteorPower[indlag:]
937 944 x = numpy.arange(0, y.size)*timeLag
938 945
939 946 #first guess
940 947 a = y[0]
941 948 tau = timeLag
942 949 #exponential fit
943 950 popt, pcov = optimize.curve_fit(self.__exponential_function, x, y, p0 = [a, tau])
944 951 y1 = self.__exponential_function(x, *popt)
945 952 #error estimation
946 953 error = sum((y - y1)**2)/(numpy.var(y)*(y.size - popt.size))
947 954
948 955 decayTime = popt[1]
949 956 riseTime = indmax*timeInterval
950 957 meteorAux[11:13] = [decayTime, error]
951 958
952 959 #Table items 7, 8 and 11
953 960 if (riseTime > 0.3): #Number 7: Echo rise exceeds 0.3s
954 961 meteorAux[-1] = 7
955 962 elif (decayTime < 2*riseTime) : #Number 8: Echo decay time less than than twice rise time
956 963 meteorAux[-1] = 8
957 964 if (error > threshError): #Number 11: Poor fit to amplitude for estimation of decay time
958 965 meteorAux[-1] = 11
959 966
960 967
961 968 except:
962 969 meteorAux[-1] = 11
963 970
964 971
965 972 listMeteors1.append(meteorAux)
966 973
967 974 return listMeteors1
968 975
969 976 #Exponential Function
970 977
971 978 def __exponential_function(self, x, a, tau):
972 979 y = a*numpy.exp(-x/tau)
973 980 return y
974 981
975 982 def __getRadialVelocity(self, listMeteors, listVolts, radialStdThresh, pairslist, timeInterval):
976 983
977 984 pairslist1 = list(pairslist)
978 985 pairslist1.append((0,1))
979 986 pairslist1.append((3,4))
980 987 numPairs = len(pairslist1)
981 988 #Time Lag
982 989 timeLag = 45*10**-3
983 990 c = 3e8
984 991 lag = numpy.ceil(timeLag/timeInterval)
985 992 freq = 30e6
986 993
987 994 listMeteors1 = []
988 995
989 996 for i in range(len(listMeteors)):
990 997 meteorAux = listMeteors[i]
991 998 if meteorAux[-1] == 0:
992 999 mStart = listMeteors[i][1]
993 1000 mPeak = listMeteors[i][2]
994 1001 mLag = mPeak - mStart + lag
995 1002
996 1003 #get the volt data between the start and end times of the meteor
997 1004 meteorVolts = listVolts[i]
998 1005 meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
999 1006
1000 1007 #Get CCF
1001 1008 allCCFs = self.__calculateCCF(meteorVolts, pairslist1, [-2,-1,0,1,2])
1002 1009
1003 1010 #Method 2
1004 1011 slopes = numpy.zeros(numPairs)
1005 1012 time = numpy.array([-2,-1,1,2])*timeInterval
1006 1013 angAllCCF = numpy.angle(allCCFs[:,[0,1,3,4],0])
1007 1014
1008 1015 #Correct phases
1009 1016 derPhaseCCF = angAllCCF[:,1:] - angAllCCF[:,0:-1]
1010 1017 indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
1011 1018
1012 1019 if indDer[0].shape[0] > 0:
1013 1020 for i in range(indDer[0].shape[0]):
1014 1021 signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i]])
1015 1022 angAllCCF[indDer[0][i],indDer[1][i]+1:] += signo*2*numpy.pi
1016 1023
1017 1024 # fit = scipy.stats.linregress(numpy.array([-2,-1,1,2])*timeInterval, numpy.array([phaseLagN2s[i],phaseLagN1s[i],phaseLag1s[i],phaseLag2s[i]]))
1018 1025 for j in range(numPairs):
1019 1026 fit = stats.linregress(time, angAllCCF[j,:])
1020 1027 slopes[j] = fit[0]
1021 1028
1022 1029 #Remove Outlier
1023 1030 # indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
1024 1031 # slopes = numpy.delete(slopes,indOut)
1025 1032 # indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
1026 1033 # slopes = numpy.delete(slopes,indOut)
1027 1034
1028 1035 radialVelocity = -numpy.mean(slopes)*(0.25/numpy.pi)*(c/freq)
1029 1036 radialError = numpy.std(slopes)*(0.25/numpy.pi)*(c/freq)
1030 1037 meteorAux[-2] = radialError
1031 1038 meteorAux[-3] = radialVelocity
1032 1039
1033 1040 #Setting Error
1034 1041 #Number 15: Radial Drift velocity or projected horizontal velocity exceeds 200 m/s
1035 1042 if numpy.abs(radialVelocity) > 200:
1036 1043 meteorAux[-1] = 15
1037 1044 #Number 12: Poor fit to CCF variation for estimation of radial drift velocity
1038 1045 elif radialError > radialStdThresh:
1039 1046 meteorAux[-1] = 12
1040 1047
1041 1048 listMeteors1.append(meteorAux)
1042 1049 return listMeteors1
1043 1050
1044 1051 def __setNewArrays(self, listMeteors, date, heiRang):
1045 1052
1046 1053 #New arrays
1047 1054 arrayMeteors = numpy.array(listMeteors)
1048 1055 arrayParameters = numpy.zeros((len(listMeteors), 14))
1049 1056
1050 1057 #Date inclusion
1051 1058 date = re.findall(r'\((.*?)\)', date)
1052 1059 date = date[0].split(',')
1053 1060 date = map(int, date)
1061
1062 if len(date)<6:
1063 date.append(0)
1064
1054 1065 date = [date[0]*10000 + date[1]*100 + date[2], date[3]*10000 + date[4]*100 + date[5]]
1055 1066 arrayDate = numpy.tile(date, (len(listMeteors), 1))
1056 1067
1057 1068 #Meteor array
1058 1069 # arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)]
1059 1070 # arrayMeteors = numpy.hstack((arrayDate, arrayMeteors))
1060 1071
1061 1072 #Parameters Array
1062 1073 arrayParameters[:,:2] = arrayDate #Date
1063 1074 arrayParameters[:,2] = heiRang[arrayMeteors[:,0].astype(int)] #Range
1064 1075 arrayParameters[:,7:9] = arrayMeteors[:,-3:-1] #Radial velocity and its error
1065 1076 arrayParameters[:,9:13] = arrayMeteors[:,7:11] #Phases
1066 1077 arrayParameters[:,-1] = arrayMeteors[:,-1] #Error
1067 1078
1068 1079
1069 1080 return arrayParameters
1070 1081
1071 1082 def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth):
1072 1083
1073 1084 arrayAOA = numpy.zeros((phases.shape[0],3))
1074 1085 cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList)
1075 1086
1076 1087 arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
1077 1088 cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
1078 1089 arrayAOA[:,2] = cosDirError
1079 1090
1080 1091 azimuthAngle = arrayAOA[:,0]
1081 1092 zenithAngle = arrayAOA[:,1]
1082 1093
1083 1094 #Setting Error
1084 1095 #Number 3: AOA not fesible
1085 1096 indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
1086 1097 error[indInvalid] = 3
1087 1098 #Number 4: Large difference in AOAs obtained from different antenna baselines
1088 1099 indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
1089 1100 error[indInvalid] = 4
1090 1101 return arrayAOA, error
1091 1102
1092 1103 def __getDirectionCosines(self, arrayPhase, pairsList):
1093 1104
1094 1105 #Initializing some variables
1095 1106 ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
1096 1107 ang_aux = ang_aux.reshape(1,ang_aux.size)
1097 1108
1098 1109 cosdir = numpy.zeros((arrayPhase.shape[0],2))
1099 1110 cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
1100 1111
1101 1112
1102 1113 for i in range(2):
1103 1114 #First Estimation
1104 1115 phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]]
1105 1116 #Dealias
1106 1117 indcsi = numpy.where(phi0_aux > numpy.pi)
1107 1118 phi0_aux[indcsi] -= 2*numpy.pi
1108 1119 indcsi = numpy.where(phi0_aux < -numpy.pi)
1109 1120 phi0_aux[indcsi] += 2*numpy.pi
1110 1121 #Direction Cosine 0
1111 1122 cosdir0[:,i] = -(phi0_aux)/(2*numpy.pi*0.5)
1112 1123
1113 1124 #Most-Accurate Second Estimation
1114 1125 phi1_aux = arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]]
1115 1126 phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
1116 1127 #Direction Cosine 1
1117 1128 cosdir1 = -(phi1_aux + ang_aux)/(2*numpy.pi*4.5)
1118 1129
1119 1130 #Searching the correct Direction Cosine
1120 1131 cosdir0_aux = cosdir0[:,i]
1121 1132 cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
1122 1133 #Minimum Distance
1123 1134 cosDiff = (cosdir1 - cosdir0_aux)**2
1124 1135 indcos = cosDiff.argmin(axis = 1)
1125 1136 #Saving Value obtained
1126 1137 cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
1127 1138
1128 1139 return cosdir0, cosdir
1129 1140
1130 1141 def __calculateAOA(self, cosdir, azimuth):
1131 1142 cosdirX = cosdir[:,0]
1132 1143 cosdirY = cosdir[:,1]
1133 1144
1134 1145 zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
1135 1146 azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth #0 deg north, 90 deg east
1136 1147 angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
1137 1148
1138 1149 return angles
1139 1150
1140 1151 def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
1141 1152
1142 1153 Ramb = 375 #Ramb = c/(2*PRF)
1143 1154 Re = 6371 #Earth Radius
1144 1155 heights = numpy.zeros(Ranges.shape)
1145 1156
1146 1157 R_aux = numpy.array([0,1,2])*Ramb
1147 1158 R_aux = R_aux.reshape(1,R_aux.size)
1148 1159
1149 1160 Ranges = Ranges.reshape(Ranges.size,1)
1150 1161
1151 1162 Ri = Ranges + R_aux
1152 1163 hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re
1153 1164
1154 1165 #Check if there is a height between 70 and 110 km
1155 1166 h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
1156 1167 ind_h = numpy.where(h_bool == 1)[0]
1157 1168
1158 1169 hCorr = hi[ind_h, :]
1159 1170 ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
1160 1171
1161 1172 hCorr = hi[ind_hCorr]
1162 1173 heights[ind_h] = hCorr
1163 1174
1164 1175 #Setting Error
1165 1176 #Number 13: Height unresolvable echo: not valid height within 70 to 110 km
1166 1177 #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
1167 1178
1168 1179 indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
1169 1180 error[indInvalid2] = 14
1170 1181 indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
1171 1182 error[indInvalid1] = 13
1172 1183
1173 1184 return heights, error
1174 1185
1175 1186 def SpectralFitting(self, getSNR = True, path=None, file=None, groupList=None):
1176 1187
1177 1188 '''
1178 1189 Function GetMoments()
1179 1190
1180 1191 Input:
1181 1192 Output:
1182 1193 Variables modified:
1183 1194 '''
1184 1195 if path != None:
1185 1196 sys.path.append(path)
1186 1197 self.dataOut.library = importlib.import_module(file)
1187 1198
1188 1199 #To be inserted as a parameter
1189 1200 groupArray = numpy.array(groupList)
1190 1201 # groupArray = numpy.array([[0,1],[2,3]])
1191 1202 self.dataOut.groupList = groupArray
1192 1203
1193 1204 nGroups = groupArray.shape[0]
1194 1205 nChannels = self.dataIn.nChannels
1195 1206 nHeights=self.dataIn.heightList.size
1196 1207
1197 1208 #Parameters Array
1198 1209 self.dataOut.data_param = None
1199 1210
1200 1211 #Set constants
1201 1212 constants = self.dataOut.library.setConstants(self.dataIn)
1202 1213 self.dataOut.constants = constants
1203 1214 M = self.dataIn.normFactor
1204 1215 N = self.dataIn.nFFTPoints
1205 1216 ippSeconds = self.dataIn.ippSeconds
1206 1217 K = self.dataIn.nIncohInt
1207 1218 pairsArray = numpy.array(self.dataIn.pairsList)
1208 1219
1209 1220 #List of possible combinations
1210 1221 listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)
1211 1222 indCross = numpy.zeros(len(list(listComb)), dtype = 'int')
1212 1223
1213 1224 if getSNR:
1214 1225 listChannels = groupArray.reshape((groupArray.size))
1215 1226 listChannels.sort()
1216 1227 noise = self.dataIn.getNoise()
1217 1228 self.dataOut.data_SNR = self.__getSNR(self.dataIn.data_spc[listChannels,:,:], noise[listChannels])
1218 1229
1219 1230 for i in range(nGroups):
1220 1231 coord = groupArray[i,:]
1221 1232
1222 1233 #Input data array
1223 1234 data = self.dataIn.data_spc[coord,:,:]/(M*N)
1224 1235 data = data.reshape((data.shape[0]*data.shape[1],data.shape[2]))
1225 1236
1226 1237 #Cross Spectra data array for Covariance Matrixes
1227 1238 ind = 0
1228 1239 for pairs in listComb:
1229 1240 pairsSel = numpy.array([coord[x],coord[y]])
1230 1241 indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0])
1231 1242 ind += 1
1232 1243 dataCross = self.dataIn.data_cspc[indCross,:,:]/(M*N)
1233 1244 dataCross = dataCross**2/K
1234 1245
1235 1246 for h in range(nHeights):
1236 1247 # print self.dataOut.heightList[h]
1237 1248
1238 1249 #Input
1239 1250 d = data[:,h]
1240 1251
1241 1252 #Covariance Matrix
1242 1253 D = numpy.diag(d**2/K)
1243 1254 ind = 0
1244 1255 for pairs in listComb:
1245 1256 #Coordinates in Covariance Matrix
1246 1257 x = pairs[0]
1247 1258 y = pairs[1]
1248 1259 #Channel Index
1249 1260 S12 = dataCross[ind,:,h]
1250 1261 D12 = numpy.diag(S12)
1251 1262 #Completing Covariance Matrix with Cross Spectras
1252 1263 D[x*N:(x+1)*N,y*N:(y+1)*N] = D12
1253 1264 D[y*N:(y+1)*N,x*N:(x+1)*N] = D12
1254 1265 ind += 1
1255 1266 Dinv=numpy.linalg.inv(D)
1256 1267 L=numpy.linalg.cholesky(Dinv)
1257 1268 LT=L.T
1258 1269
1259 1270 dp = numpy.dot(LT,d)
1260 1271
1261 1272 #Initial values
1262 1273 data_spc = self.dataIn.data_spc[coord,:,h]
1263 1274
1264 1275 if (h>0)and(error1[3]<5):
1265 1276 p0 = self.dataOut.data_param[i,:,h-1]
1266 1277 else:
1267 1278 p0 = numpy.array(self.dataOut.library.initialValuesFunction(data_spc, constants, i))
1268 1279
1269 1280 try:
1270 1281 #Least Squares
1271 1282 minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True)
1272 1283 # minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))
1273 1284 #Chi square error
1274 1285 error0 = numpy.sum(infodict['fvec']**2)/(2*N)
1275 1286 #Error with Jacobian
1276 1287 error1 = self.dataOut.library.errorFunction(minp,constants,LT)
1277 1288 except:
1278 1289 minp = p0*numpy.nan
1279 1290 error0 = numpy.nan
1280 1291 error1 = p0*numpy.nan
1281 1292
1282 1293 #Save
1283 if self.dataOut.data_param is None:
1294 if self.dataOut.data_param == None:
1284 1295 self.dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan
1285 1296 self.dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan
1286 1297
1287 1298 self.dataOut.data_error[i,:,h] = numpy.hstack((error0,error1))
1288 1299 self.dataOut.data_param[i,:,h] = minp
1289 1300 return
1290 1301
1291 1302 def __residFunction(self, p, dp, LT, constants):
1292 1303
1293 1304 fm = self.dataOut.library.modelFunction(p, constants)
1294 1305 fmp=numpy.dot(LT,fm)
1295 1306
1296 1307 return dp-fmp
1297 1308
1298 1309 def __getSNR(self, z, noise):
1299 1310
1300 1311 avg = numpy.average(z, axis=1)
1301 1312 SNR = (avg.T-noise)/noise
1302 1313 SNR = SNR.T
1303 1314 return SNR
1304 1315
1305 1316 def __chisq(p,chindex,hindex):
1306 1317 #similar to Resid but calculates CHI**2
1307 1318 [LT,d,fm]=setupLTdfm(p,chindex,hindex)
1308 1319 dp=numpy.dot(LT,d)
1309 1320 fmp=numpy.dot(LT,fm)
1310 1321 chisq=numpy.dot((dp-fmp).T,(dp-fmp))
1311 1322 return chisq
1312 1323
1313 1324
1314 1325
1315 1326 class WindProfiler(Operation):
1316 1327
1317 1328 __isConfig = False
1318 1329
1319 1330 __initime = None
1320 1331 __lastdatatime = None
1321 1332 __integrationtime = None
1322 1333
1323 1334 __buffer = None
1324 1335
1325 1336 __dataReady = False
1326 1337
1327 1338 __firstdata = None
1328 1339
1329 1340 n = None
1330 1341
1331 1342 def __init__(self):
1332 1343 Operation.__init__(self)
1333 1344
1334 1345 def __calculateCosDir(self, elev, azim):
1335 1346 zen = (90 - elev)*numpy.pi/180
1336 1347 azim = azim*numpy.pi/180
1337 1348 cosDirX = numpy.sqrt((1-numpy.cos(zen)**2)/((1+numpy.tan(azim)**2)))
1338 1349 cosDirY = numpy.sqrt(1-numpy.cos(zen)**2-cosDirX**2)
1339 1350
1340 1351 signX = numpy.sign(numpy.cos(azim))
1341 1352 signY = numpy.sign(numpy.sin(azim))
1342 1353
1343 1354 cosDirX = numpy.copysign(cosDirX, signX)
1344 1355 cosDirY = numpy.copysign(cosDirY, signY)
1345 1356 return cosDirX, cosDirY
1346 1357
1347 1358 def __calculateAngles(self, theta_x, theta_y, azimuth):
1348 1359
1349 1360 dir_cosw = numpy.sqrt(1-theta_x**2-theta_y**2)
1350 1361 zenith_arr = numpy.arccos(dir_cosw)
1351 1362 azimuth_arr = numpy.arctan2(theta_x,theta_y) + azimuth*math.pi/180
1352 1363
1353 1364 dir_cosu = numpy.sin(azimuth_arr)*numpy.sin(zenith_arr)
1354 1365 dir_cosv = numpy.cos(azimuth_arr)*numpy.sin(zenith_arr)
1355 1366
1356 1367 return azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw
1357 1368
1358 1369 def __calculateMatA(self, dir_cosu, dir_cosv, dir_cosw, horOnly):
1359 1370
1360 1371 #
1361 1372 if horOnly:
1362 1373 A = numpy.c_[dir_cosu,dir_cosv]
1363 1374 else:
1364 1375 A = numpy.c_[dir_cosu,dir_cosv,dir_cosw]
1365 1376 A = numpy.asmatrix(A)
1366 1377 A1 = numpy.linalg.inv(A.transpose()*A)*A.transpose()
1367 1378
1368 1379 return A1
1369 1380
1370 1381 def __correctValues(self, heiRang, phi, velRadial, SNR):
1371 1382 listPhi = phi.tolist()
1372 1383 maxid = listPhi.index(max(listPhi))
1373 1384 minid = listPhi.index(min(listPhi))
1374 1385
1375 1386 rango = range(len(phi))
1376 1387 # rango = numpy.delete(rango,maxid)
1377 1388
1378 1389 heiRang1 = heiRang*math.cos(phi[maxid])
1379 1390 heiRangAux = heiRang*math.cos(phi[minid])
1380 1391 indOut = (heiRang1 < heiRangAux[0]).nonzero()
1381 1392 heiRang1 = numpy.delete(heiRang1,indOut)
1382 1393
1383 1394 velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
1384 1395 SNR1 = numpy.zeros([len(phi),len(heiRang1)])
1385 1396
1386 1397 for i in rango:
1387 1398 x = heiRang*math.cos(phi[i])
1388 1399 y1 = velRadial[i,:]
1389 1400 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
1390 1401
1391 1402 x1 = heiRang1
1392 1403 y11 = f1(x1)
1393 1404
1394 1405 y2 = SNR[i,:]
1395 1406 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
1396 1407 y21 = f2(x1)
1397 1408
1398 1409 velRadial1[i,:] = y11
1399 1410 SNR1[i,:] = y21
1400 1411
1401 1412 return heiRang1, velRadial1, SNR1
1402 1413
1403 1414 def __calculateVelUVW(self, A, velRadial):
1404 1415
1405 1416 #Operacion Matricial
1406 1417 # velUVW = numpy.zeros((velRadial.shape[1],3))
1407 1418 # for ind in range(velRadial.shape[1]):
1408 1419 # velUVW[ind,:] = numpy.dot(A,velRadial[:,ind])
1409 1420 # velUVW = velUVW.transpose()
1410 1421 velUVW = numpy.zeros((A.shape[0],velRadial.shape[1]))
1411 1422 velUVW[:,:] = numpy.dot(A,velRadial)
1412 1423
1413 1424
1414 1425 return velUVW
1415 1426
1416 1427 def techniqueDBS(self, velRadial0, dirCosx, disrCosy, azimuth, correct, horizontalOnly, heiRang, SNR0):
1417 1428 """
1418 1429 Function that implements Doppler Beam Swinging (DBS) technique.
1419 1430
1420 1431 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
1421 1432 Direction correction (if necessary), Ranges and SNR
1422 1433
1423 1434 Output: Winds estimation (Zonal, Meridional and Vertical)
1424 1435
1425 1436 Parameters affected: Winds, height range, SNR
1426 1437 """
1427 1438 azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(dirCosx, disrCosy, azimuth)
1428 1439 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correct*velRadial0, SNR0)
1429 1440 A = self.__calculateMatA(dir_cosu, dir_cosv, dir_cosw, horizontalOnly)
1430 1441
1431 1442 #Calculo de Componentes de la velocidad con DBS
1432 1443 winds = self.__calculateVelUVW(A,velRadial1)
1433 1444
1434 1445 return winds, heiRang1, SNR1
1435 1446
1436 1447 def __calculateDistance(self, posx, posy, pairsCrossCorr, pairsList, pairs, azimuth = None):
1437 1448
1438 1449 posx = numpy.asarray(posx)
1439 1450 posy = numpy.asarray(posy)
1440 1451
1441 1452 #Rotacion Inversa para alinear con el azimuth
1442 1453 if azimuth!= None:
1443 1454 azimuth = azimuth*math.pi/180
1444 1455 posx1 = posx*math.cos(azimuth) + posy*math.sin(azimuth)
1445 1456 posy1 = -posx*math.sin(azimuth) + posy*math.cos(azimuth)
1446 1457 else:
1447 1458 posx1 = posx
1448 1459 posy1 = posy
1449 1460
1450 1461 #Calculo de Distancias
1451 1462 distx = numpy.zeros(pairsCrossCorr.size)
1452 1463 disty = numpy.zeros(pairsCrossCorr.size)
1453 1464 dist = numpy.zeros(pairsCrossCorr.size)
1454 1465 ang = numpy.zeros(pairsCrossCorr.size)
1455 1466
1456 1467 for i in range(pairsCrossCorr.size):
1457 1468 distx[i] = posx1[pairsList[pairsCrossCorr[i]][1]] - posx1[pairsList[pairsCrossCorr[i]][0]]
1458 1469 disty[i] = posy1[pairsList[pairsCrossCorr[i]][1]] - posy1[pairsList[pairsCrossCorr[i]][0]]
1459 1470 dist[i] = numpy.sqrt(distx[i]**2 + disty[i]**2)
1460 1471 ang[i] = numpy.arctan2(disty[i],distx[i])
1461 1472 #Calculo de Matrices
1462 1473 nPairs = len(pairs)
1463 1474 ang1 = numpy.zeros((nPairs, 2, 1))
1464 1475 dist1 = numpy.zeros((nPairs, 2, 1))
1465 1476
1466 1477 for j in range(nPairs):
1467 1478 dist1[j,0,0] = dist[pairs[j][0]]
1468 1479 dist1[j,1,0] = dist[pairs[j][1]]
1469 1480 ang1[j,0,0] = ang[pairs[j][0]]
1470 1481 ang1[j,1,0] = ang[pairs[j][1]]
1471 1482
1472 1483 return distx,disty, dist1,ang1
1473 1484
1474 1485 def __calculateVelVer(self, phase, lagTRange, _lambda):
1475 1486
1476 1487 Ts = lagTRange[1] - lagTRange[0]
1477 1488 velW = -_lambda*phase/(4*math.pi*Ts)
1478 1489
1479 1490 return velW
1480 1491
1481 1492 def __calculateVelHorDir(self, dist, tau1, tau2, ang):
1482 1493 nPairs = tau1.shape[0]
1483 1494 vel = numpy.zeros((nPairs,3,tau1.shape[2]))
1484 1495
1485 1496 angCos = numpy.cos(ang)
1486 1497 angSin = numpy.sin(ang)
1487 1498
1488 1499 vel0 = dist*tau1/(2*tau2**2)
1489 1500 vel[:,0,:] = (vel0*angCos).sum(axis = 1)
1490 1501 vel[:,1,:] = (vel0*angSin).sum(axis = 1)
1491 1502
1492 1503 ind = numpy.where(numpy.isinf(vel))
1493 1504 vel[ind] = numpy.nan
1494 1505
1495 1506 return vel
1496 1507
1497 1508 def __getPairsAutoCorr(self, pairsList, nChannels):
1498 1509
1499 1510 pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
1500 1511
1501 1512 for l in range(len(pairsList)):
1502 1513 firstChannel = pairsList[l][0]
1503 1514 secondChannel = pairsList[l][1]
1504 1515
1505 1516 #Obteniendo pares de Autocorrelacion
1506 1517 if firstChannel == secondChannel:
1507 1518 pairsAutoCorr[firstChannel] = int(l)
1508 1519
1509 1520 pairsAutoCorr = pairsAutoCorr.astype(int)
1510 1521
1511 1522 pairsCrossCorr = range(len(pairsList))
1512 1523 pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
1513 1524
1514 1525 return pairsAutoCorr, pairsCrossCorr
1515 1526
1516 1527 def techniqueSA(self, pairsSelected, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, lagTRange, correctFactor):
1517 1528 """
1518 1529 Function that implements Spaced Antenna (SA) technique.
1519 1530
1520 1531 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
1521 1532 Direction correction (if necessary), Ranges and SNR
1522 1533
1523 1534 Output: Winds estimation (Zonal, Meridional and Vertical)
1524 1535
1525 1536 Parameters affected: Winds
1526 1537 """
1527 1538 #Cross Correlation pairs obtained
1528 1539 pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)
1529 1540 pairsArray = numpy.array(pairsList)[pairsCrossCorr]
1530 1541 pairsSelArray = numpy.array(pairsSelected)
1531 1542 pairs = []
1532 1543
1533 1544 #Wind estimation pairs obtained
1534 1545 for i in range(pairsSelArray.shape[0]/2):
1535 1546 ind1 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i], axis = 1))[0][0]
1536 1547 ind2 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i + 1], axis = 1))[0][0]
1537 1548 pairs.append((ind1,ind2))
1538 1549
1539 1550 indtau = tau.shape[0]/2
1540 1551 tau1 = tau[:indtau,:]
1541 1552 tau2 = tau[indtau:-1,:]
1542 1553 tau1 = tau1[pairs,:]
1543 1554 tau2 = tau2[pairs,:]
1544 1555 phase1 = tau[-1,:]
1545 1556
1546 1557 #---------------------------------------------------------------------
1547 1558 #Metodo Directo
1548 1559 distx, disty, dist, ang = self.__calculateDistance(position_x, position_y, pairsCrossCorr, pairsList, pairs,azimuth)
1549 1560 winds = self.__calculateVelHorDir(dist, tau1, tau2, ang)
1550 1561 winds = stats.nanmean(winds, axis=0)
1551 1562 #---------------------------------------------------------------------
1552 1563 #Metodo General
1553 1564 # distx, disty, dist = self.calculateDistance(position_x,position_y,pairsCrossCorr, pairsList, azimuth)
1554 1565 # #Calculo Coeficientes de Funcion de Correlacion
1555 1566 # F,G,A,B,H = self.calculateCoef(tau1,tau2,distx,disty,n)
1556 1567 # #Calculo de Velocidades
1557 1568 # winds = self.calculateVelUV(F,G,A,B,H)
1558 1569
1559 1570 #---------------------------------------------------------------------
1560 1571 winds[2,:] = self.__calculateVelVer(phase1, lagTRange, _lambda)
1561 1572 winds = correctFactor*winds
1562 1573 return winds
1563 1574
1564 1575 def __checkTime(self, currentTime, paramInterval, outputInterval):
1565 1576
1566 1577 dataTime = currentTime + paramInterval
1567 1578 deltaTime = dataTime - self.__initime
1568 1579
1569 1580 if deltaTime >= outputInterval or deltaTime < 0:
1570 1581 self.__dataReady = True
1571 1582 return
1572 1583
1573 1584 def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax):
1574 1585 '''
1575 1586 Function that implements winds estimation technique with detected meteors.
1576 1587
1577 1588 Input: Detected meteors, Minimum meteor quantity to wind estimation
1578 1589
1579 1590 Output: Winds estimation (Zonal and Meridional)
1580 1591
1581 1592 Parameters affected: Winds
1582 1593 '''
1583 print arrayMeteor.shape
1594 # print arrayMeteor.shape
1584 1595 #Settings
1585 1596 nInt = (heightMax - heightMin)/2
1597 # print nInt
1598 nInt = int(nInt)
1599 # print nInt
1586 1600 winds = numpy.zeros((2,nInt))*numpy.nan
1587 1601
1588 1602 #Filter errors
1589 1603 error = numpy.where(arrayMeteor[0,:,-1] == 0)[0]
1590 1604 finalMeteor = arrayMeteor[0,error,:]
1591 1605
1592 1606 #Meteor Histogram
1593 1607 finalHeights = finalMeteor[:,3]
1594 1608 hist = numpy.histogram(finalHeights, bins = nInt, range = (heightMin,heightMax))
1595 1609 nMeteorsPerI = hist[0]
1596 1610 heightPerI = hist[1]
1597 1611
1598 1612 #Sort of meteors
1599 1613 indSort = finalHeights.argsort()
1600 1614 finalMeteor2 = finalMeteor[indSort,:]
1601 1615
1602 1616 # Calculating winds
1603 1617 ind1 = 0
1604 1618 ind2 = 0
1605 1619
1606 1620 for i in range(nInt):
1607 1621 nMet = nMeteorsPerI[i]
1608 1622 ind1 = ind2
1609 1623 ind2 = ind1 + nMet
1610 1624
1611 1625 meteorAux = finalMeteor2[ind1:ind2,:]
1612 1626
1613 1627 if meteorAux.shape[0] >= meteorThresh:
1614 1628 vel = meteorAux[:, 7]
1615 1629 zen = meteorAux[:, 5]*numpy.pi/180
1616 1630 azim = meteorAux[:, 4]*numpy.pi/180
1617 1631
1618 1632 n = numpy.cos(zen)
1619 1633 # m = (1 - n**2)/(1 - numpy.tan(azim)**2)
1620 1634 # l = m*numpy.tan(azim)
1621 1635 l = numpy.sin(zen)*numpy.sin(azim)
1622 1636 m = numpy.sin(zen)*numpy.cos(azim)
1623 1637
1624 1638 A = numpy.vstack((l, m)).transpose()
1625 1639 A1 = numpy.dot(numpy.linalg.inv( numpy.dot(A.transpose(),A) ),A.transpose())
1626 1640 windsAux = numpy.dot(A1, vel)
1627 1641
1628 1642 winds[0,i] = windsAux[0]
1629 1643 winds[1,i] = windsAux[1]
1630 1644
1631 1645 return winds, heightPerI[:-1]
1632 1646
1633 1647 def run(self, dataOut, technique, **kwargs):
1634 1648
1635 1649 param = dataOut.data_param
1636 # if dataOut.abscissaList != None:
1637 # absc = dataOut.abscissaList[:-1]
1650 if dataOut.abscissaList != None:
1651 absc = dataOut.abscissaList[:-1]
1638 1652 noise = dataOut.noise
1639 1653 heightList = dataOut.heightList
1640 1654 SNR = dataOut.data_SNR
1641 1655
1642 1656 if technique == 'DBS':
1643 1657
1644 1658 if kwargs.has_key('dirCosx') and kwargs.has_key('dirCosy'):
1645 1659 theta_x = numpy.array(kwargs['dirCosx'])
1646 1660 theta_y = numpy.array(kwargs['dirCosy'])
1647 1661 else:
1648 1662 elev = numpy.array(kwargs['elevation'])
1649 1663 azim = numpy.array(kwargs['azimuth'])
1650 1664 theta_x, theta_y = self.__calculateCosDir(elev, azim)
1651 1665 azimuth = kwargs['correctAzimuth']
1652 1666 if kwargs.has_key('horizontalOnly'):
1653 1667 horizontalOnly = kwargs['horizontalOnly']
1654 1668 else: horizontalOnly = False
1655 1669 if kwargs.has_key('correctFactor'):
1656 1670 correctFactor = kwargs['correctFactor']
1657 1671 else: correctFactor = 1
1658 1672 if kwargs.has_key('channelList'):
1659 1673 channelList = kwargs['channelList']
1660 1674 if len(channelList) == 2:
1661 1675 horizontalOnly = True
1662 1676 arrayChannel = numpy.array(channelList)
1663 1677 param = param[arrayChannel,:,:]
1664 1678 theta_x = theta_x[arrayChannel]
1665 1679 theta_y = theta_y[arrayChannel]
1666 1680
1667 1681 velRadial0 = param[:,1,:] #Radial velocity
1668 1682 dataOut.data_output, dataOut.heightList, dataOut.data_SNR = self.techniqueDBS(velRadial0, theta_x, theta_y, azimuth, correctFactor, horizontalOnly, heightList, SNR) #DBS Function
1669 1683 dataOut.utctimeInit = dataOut.utctime
1670 1684 dataOut.outputInterval = dataOut.timeInterval
1671 1685
1672 1686 elif technique == 'SA':
1673 1687
1674 1688 #Parameters
1675 1689 position_x = kwargs['positionX']
1676 1690 position_y = kwargs['positionY']
1677 1691 azimuth = kwargs['azimuth']
1678 1692
1679 1693 if kwargs.has_key('crosspairsList'):
1680 1694 pairs = kwargs['crosspairsList']
1681 1695 else:
1682 1696 pairs = None
1683 1697
1684 1698 if kwargs.has_key('correctFactor'):
1685 1699 correctFactor = kwargs['correctFactor']
1686 1700 else:
1687 1701 correctFactor = 1
1688 1702
1689 1703 tau = dataOut.data_param
1690 1704 _lambda = dataOut.C/dataOut.frequency
1691 1705 pairsList = dataOut.groupList
1692 1706 nChannels = dataOut.nChannels
1693 1707
1694 1708 dataOut.data_output = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor)
1695 1709 dataOut.utctimeInit = dataOut.utctime
1696 1710 dataOut.outputInterval = dataOut.timeInterval
1697 1711
1698 1712 elif technique == 'Meteors':
1699 1713 dataOut.flagNoData = True
1700 1714 self.__dataReady = False
1701 1715
1702 1716 if kwargs.has_key('nHours'):
1703 1717 nHours = kwargs['nHours']
1704 1718 else:
1705 1719 nHours = 1
1706 1720
1707 1721 if kwargs.has_key('meteorsPerBin'):
1708 1722 meteorThresh = kwargs['meteorsPerBin']
1709 1723 else:
1710 1724 meteorThresh = 6
1711 1725
1712 1726 if kwargs.has_key('hmin'):
1713 1727 hmin = kwargs['hmin']
1714 1728 else: hmin = 70
1715 1729 if kwargs.has_key('hmax'):
1716 1730 hmax = kwargs['hmax']
1717 1731 else: hmax = 110
1718 1732
1719 1733 dataOut.outputInterval = nHours*3600
1720 1734
1721 1735 if self.__isConfig == False:
1722 1736 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
1723 1737 #Get Initial LTC time
1724 1738 self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
1725 1739 self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
1726 1740
1727 1741 self.__isConfig = True
1728 1742
1729 if self.__buffer is None:
1743 if self.__buffer == None:
1730 1744 self.__buffer = dataOut.data_param
1731 1745 self.__firstdata = copy.copy(dataOut)
1732 1746
1733 1747 else:
1734 1748 self.__buffer = numpy.hstack((self.__buffer, dataOut.data_param))
1735 1749
1736 1750 self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
1737 1751
1738 1752 if self.__dataReady:
1739 1753 dataOut.utctimeInit = self.__initime
1740 1754
1741 1755 self.__initime += dataOut.outputInterval #to erase time offset
1742 1756
1743 1757 dataOut.data_output, dataOut.heightList = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax)
1744 1758 dataOut.flagNoData = False
1745 1759 self.__buffer = None
1746 1760
1747 1761 return
1748 1762
1749 1763 class EWDriftsEstimation(Operation):
1750 1764
1751 1765
1752 1766 def __init__(self):
1753 1767 Operation.__init__(self)
1754 1768
1755 1769 def __correctValues(self, heiRang, phi, velRadial, SNR):
1756 1770 listPhi = phi.tolist()
1757 1771 maxid = listPhi.index(max(listPhi))
1758 1772 minid = listPhi.index(min(listPhi))
1759 1773
1760 1774 rango = range(len(phi))
1761 1775 # rango = numpy.delete(rango,maxid)
1762 1776
1763 1777 heiRang1 = heiRang*math.cos(phi[maxid])
1764 1778 heiRangAux = heiRang*math.cos(phi[minid])
1765 1779 indOut = (heiRang1 < heiRangAux[0]).nonzero()
1766 1780 heiRang1 = numpy.delete(heiRang1,indOut)
1767 1781
1768 1782 velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
1769 1783 SNR1 = numpy.zeros([len(phi),len(heiRang1)])
1770 1784
1771 1785 for i in rango:
1772 1786 x = heiRang*math.cos(phi[i])
1773 1787 y1 = velRadial[i,:]
1774 1788 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
1775 1789
1776 1790 x1 = heiRang1
1777 1791 y11 = f1(x1)
1778 1792
1779 1793 y2 = SNR[i,:]
1780 1794 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
1781 1795 y21 = f2(x1)
1782 1796
1783 1797 velRadial1[i,:] = y11
1784 1798 SNR1[i,:] = y21
1785 1799
1786 1800 return heiRang1, velRadial1, SNR1
1787 1801
1788 1802 def run(self, dataOut, zenith, zenithCorrection):
1789 1803 heiRang = dataOut.heightList
1790 1804 velRadial = dataOut.data_param[:,3,:]
1791 1805 SNR = dataOut.data_SNR
1792 1806
1793 1807 zenith = numpy.array(zenith)
1794 1808 zenith -= zenithCorrection
1795 1809 zenith *= numpy.pi/180
1796 1810
1797 1811 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR)
1798 1812
1799 1813 alp = zenith[0]
1800 1814 bet = zenith[1]
1801 1815
1802 1816 w_w = velRadial1[0,:]
1803 1817 w_e = velRadial1[1,:]
1804 1818
1805 1819 w = (w_w*numpy.sin(bet) - w_e*numpy.sin(alp))/(numpy.cos(alp)*numpy.sin(bet) - numpy.cos(bet)*numpy.sin(alp))
1806 1820 u = (w_w*numpy.cos(bet) - w_e*numpy.cos(alp))/(numpy.sin(alp)*numpy.cos(bet) - numpy.sin(bet)*numpy.cos(alp))
1807 1821
1808 1822 winds = numpy.vstack((u,w))
1809 1823
1810 1824 dataOut.heightList = heiRang1
1811 1825 dataOut.data_output = winds
1812 1826 dataOut.data_SNR = SNR1
1813 1827
1814 1828 dataOut.utctimeInit = dataOut.utctime
1815 1829 dataOut.outputInterval = dataOut.timeInterval
1816 1830 return
1817 1831
1818 1832 class PhaseCalibration(Operation):
1819 1833
1820 1834 __buffer = None
1821 1835
1822 1836 __initime = None
1823 1837
1824 1838 __dataReady = False
1825 1839
1826 1840 __isConfig = False
1827 1841
1828 1842 def __checkTime(self, currentTime, initTime, paramInterval, outputInterval):
1829 1843
1830 1844 dataTime = currentTime + paramInterval
1831 1845 deltaTime = dataTime - initTime
1832 1846
1833 1847 if deltaTime >= outputInterval or deltaTime < 0:
1834 1848 return True
1835 1849
1836 1850 return False
1837 1851
1838 1852 def __getGammas(self, pairs, k, d, phases):
1839 1853 gammas = numpy.zeros(2)
1840
1854
1841 1855 for i in range(len(pairs)):
1842 1856
1843 1857 pairi = pairs[i]
1844 1858
1845 1859 #Calculating gamma
1846 1860 jdcos = phases[:,pairi[1]]/(k*d[pairi[1]])
1847 1861 jgamma = numpy.angle(numpy.exp(1j*(k*d[pairi[0]]*jdcos - phases[:,pairi[0]])))
1848 1862
1849 1863 #Revised distribution
1850 1864 jgammaArray = numpy.hstack((jgamma,jgamma+0.5*numpy.pi,jgamma-0.5*numpy.pi))
1851 1865
1852 1866 #Histogram
1853 1867 nBins = 64.0
1854 1868 rmin = -0.5*numpy.pi
1855 1869 rmax = 0.5*numpy.pi
1856 1870 phaseHisto = numpy.histogram(jgammaArray, bins=nBins, range=(rmin,rmax))
1857 1871
1858 1872 meteorsY = phaseHisto[0]
1859 1873 phasesX = phaseHisto[1][:-1]
1860 1874 width = phasesX[1] - phasesX[0]
1861 1875 phasesX += width/2
1862
1876
1863 1877 #Gaussian aproximation
1864 1878 bpeak = meteorsY.argmax()
1865 1879 peak = meteorsY.max()
1866 1880 jmin = bpeak - 5
1867 1881 jmax = bpeak + 5 + 1
1868 1882
1869 1883 if jmin<0:
1870 1884 jmin = 0
1871 1885 jmax = 6
1872 1886 elif jmax > meteorsY.size:
1873 1887 jmin = meteorsY.size - 6
1874 1888 jmax = meteorsY.size
1875 1889
1876 1890 x0 = numpy.array([peak,bpeak,50])
1877 1891 coeff = optimize.leastsq(self.__residualFunction, x0, args=(meteorsY[jmin:jmax], phasesX[jmin:jmax]))
1878 1892
1879 1893 #Gammas
1880 1894 gammas[i] = coeff[0][1]
1881 # gammas[i] = bpeak
1882
1895
1883 1896 return gammas
1884 1897
1885 1898 def __residualFunction(self, coeffs, y, t):
1886 1899
1887 1900 return y - self.__gauss_function(t, coeffs)
1888 1901
1889 1902 def __gauss_function(self, t, coeffs):
1890 1903
1891 1904 return coeffs[0]*numpy.exp(-0.5*((t - coeffs[1]) / coeffs[2])**2)
1892 1905
1893 1906 def __getPhases(self, azimuth, h, pairsList, d, gammas, meteorsArray):
1894 1907 meteorOps = MeteorOperations()
1895 1908 nchan = 4
1896 1909 pairx = pairsList[0]
1897 1910 pairy = pairsList[1]
1898 1911 center_xangle = 0
1899 1912 center_yangle = 0
1900 range_angle = numpy.array([8*numpy.pi,numpy.pi,numpy.pi/2,numpy.pi/4])
1913 range_angle = numpy.array([10*numpy.pi,numpy.pi,numpy.pi/2,numpy.pi/4])
1901 1914 ntimes = len(range_angle)
1902 1915
1903 1916 nstepsx = 20.0
1904 1917 nstepsy = 20.0
1905 1918
1906 1919 for iz in range(ntimes):
1907 1920 min_xangle = -range_angle[iz]/2 + center_xangle
1908 1921 max_xangle = range_angle[iz]/2 + center_xangle
1909 1922 min_yangle = -range_angle[iz]/2 + center_yangle
1910 1923 max_yangle = range_angle[iz]/2 + center_yangle
1911 1924
1912 1925 inc_x = (max_xangle-min_xangle)/nstepsx
1913 1926 inc_y = (max_yangle-min_yangle)/nstepsy
1914 1927
1915 1928 alpha_y = numpy.arange(nstepsy)*inc_y + min_yangle
1916 1929 alpha_x = numpy.arange(nstepsx)*inc_x + min_xangle
1917 1930 penalty = numpy.zeros((nstepsx,nstepsy))
1918 1931 jph_array = numpy.zeros((nchan,nstepsx,nstepsy))
1919 1932 jph = numpy.zeros(nchan)
1920 1933
1921 1934 # Iterations looking for the offset
1922 1935 for iy in range(int(nstepsy)):
1923 1936 for ix in range(int(nstepsx)):
1924 1937 jph[pairy[1]] = alpha_y[iy]
1925 1938 jph[pairy[0]] = -gammas[1] + alpha_y[iy]*d[pairy[0]]/d[pairy[1]]
1926 1939
1927 1940 jph[pairx[1]] = alpha_x[ix]
1928 1941 jph[pairx[0]] = -gammas[0] + alpha_x[ix]*d[pairx[0]]/d[pairx[1]]
1929 1942
1930 1943 jph_array[:,ix,iy] = jph
1931 1944
1932 1945 meteorsArray1 = meteorOps.getMeteorParams(meteorsArray, azimuth, h, pairsList, jph)
1933 1946 error = meteorsArray1[:,-1]
1934 1947 ind1 = numpy.where(error==0)[0]
1935 1948 penalty[ix,iy] = ind1.size
1936 1949
1937 1950 i,j = numpy.unravel_index(penalty.argmax(), penalty.shape)
1938 1951 phOffset = jph_array[:,i,j]
1939 1952
1940 1953 center_xangle = phOffset[pairx[1]]
1941 1954 center_yangle = phOffset[pairy[1]]
1942 1955
1943 1956 phOffset = numpy.angle(numpy.exp(1j*jph_array[:,i,j]))
1944 1957 phOffset = phOffset*180/numpy.pi
1945 1958 return phOffset
1946 1959
1947 1960
1948 1961 def run(self, dataOut, pairs, distances, hmin, hmax, nHours = None):
1949 1962
1950 1963 dataOut.flagNoData = True
1951 1964 self.__dataReady = False
1952 1965
1953 1966 if nHours == None:
1954 1967 nHours = 1
1955 1968
1956 1969 dataOut.outputInterval = nHours*3600
1957 1970
1958 1971 if self.__isConfig == False:
1959 1972 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
1960 1973 #Get Initial LTC time
1961 1974 self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
1962 1975 self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
1963 1976
1964 1977 self.__isConfig = True
1965 1978
1966 if self.__buffer is None:
1979 if self.__buffer == None:
1967 1980 self.__buffer = dataOut.data_param.copy()
1968 1981
1969 1982 else:
1970 1983 self.__buffer = numpy.hstack((self.__buffer, dataOut.data_param))
1971 1984
1972 1985 self.__dataReady = self.__checkTime(dataOut.utctime, self.__initime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
1973 1986
1974 1987 if self.__dataReady:
1975 1988 dataOut.utctimeInit = self.__initime
1976 1989 self.__initime += dataOut.outputInterval #to erase time offset
1977 1990
1978 1991 freq = dataOut.frequency
1979 1992 c = dataOut.C #m/s
1980 1993 lamb = c/freq
1981 1994 k = 2*numpy.pi/lamb
1982 1995 azimuth = 0
1983 1996 h = (hmin, hmax)
1984 1997 pairsList = ((0,3),(1,2))
1985 1998
1986 1999 meteorsArray = self.__buffer[0,:,:]
1987 2000 error = meteorsArray[:,-1]
1988 2001 boolError = (error==0)|(error==3)|(error==4)|(error==13)|(error==14)
1989 2002 ind1 = numpy.where(boolError)[0]
1990 2003 meteorsArray = meteorsArray[ind1,:]
1991 2004 meteorsArray[:,-1] = 0
1992 2005 phases = meteorsArray[:,9:13]
1993 2006
1994 2007 #Calculate Gammas
1995 2008 gammas = self.__getGammas(pairs, k, distances, phases)
1996 2009 # gammas = numpy.array([-21.70409463,45.76935864])*numpy.pi/180
1997 2010 #Calculate Phases
1998 2011 phasesOff = self.__getPhases(azimuth, h, pairsList, distances, gammas, meteorsArray)
1999 2012 phasesOff = phasesOff.reshape((1,phasesOff.size))
2000 2013 dataOut.data_output = -phasesOff
2001 2014 dataOut.flagNoData = False
2002 2015 self.__buffer = None
2003 2016
2004 2017
2005 2018 return
2006 2019
2007 2020 class MeteorOperations():
2008 2021
2009 2022 def __init__(self):
2010 2023
2011 2024 return
2012 2025
2013 2026 def getMeteorParams(self, arrayParameters0, azimuth, h, pairsList, jph):
2014 2027
2015 2028 arrayParameters = arrayParameters0.copy()
2016 2029 hmin = h[0]
2017 2030 hmax = h[1]
2018 2031
2019 2032 #Calculate AOA (Error N 3, 4)
2020 2033 #JONES ET AL. 1998
2021 2034 AOAthresh = numpy.pi/8
2022 2035 error = arrayParameters[:,-1]
2023 2036 phases = -arrayParameters[:,9:13] + jph
2024 2037 arrayParameters[:,4:7], arrayParameters[:,-1] = self.__getAOA(phases, pairsList, error, AOAthresh, azimuth)
2025 2038
2026 2039 #Calculate Heights (Error N 13 and 14)
2027 2040 error = arrayParameters[:,-1]
2028 2041 Ranges = arrayParameters[:,2]
2029 2042 zenith = arrayParameters[:,5]
2030 2043 arrayParameters[:,3], arrayParameters[:,-1] = self.__getHeights(Ranges, zenith, error, hmin, hmax)
2031 2044
2032 2045 #----------------------- Get Final data ------------------------------------
2033 2046 # error = arrayParameters[:,-1]
2034 2047 # ind1 = numpy.where(error==0)[0]
2035 2048 # arrayParameters = arrayParameters[ind1,:]
2036 2049
2037 2050 return arrayParameters
2038 2051
2039 2052 def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth):
2040 2053
2041 2054 arrayAOA = numpy.zeros((phases.shape[0],3))
2042 2055 cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList)
2043 2056
2044 2057 arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
2045 2058 cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
2046 2059 arrayAOA[:,2] = cosDirError
2047 2060
2048 2061 azimuthAngle = arrayAOA[:,0]
2049 2062 zenithAngle = arrayAOA[:,1]
2050 2063
2051 2064 #Setting Error
2052 2065 #Number 3: AOA not fesible
2053 2066 indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
2054 2067 error[indInvalid] = 3
2055 2068 #Number 4: Large difference in AOAs obtained from different antenna baselines
2056 2069 indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
2057 2070 error[indInvalid] = 4
2058 2071 return arrayAOA, error
2059 2072
2060 2073 def __getDirectionCosines(self, arrayPhase, pairsList):
2061 2074
2062 2075 #Initializing some variables
2063 2076 ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
2064 2077 ang_aux = ang_aux.reshape(1,ang_aux.size)
2065 2078
2066 2079 cosdir = numpy.zeros((arrayPhase.shape[0],2))
2067 2080 cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
2068 2081
2069 2082
2070 2083 for i in range(2):
2071 2084 #First Estimation
2072 2085 phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]]
2073 2086 #Dealias
2074 2087 indcsi = numpy.where(phi0_aux > numpy.pi)
2075 2088 phi0_aux[indcsi] -= 2*numpy.pi
2076 2089 indcsi = numpy.where(phi0_aux < -numpy.pi)
2077 2090 phi0_aux[indcsi] += 2*numpy.pi
2078 2091 #Direction Cosine 0
2079 2092 cosdir0[:,i] = -(phi0_aux)/(2*numpy.pi*0.5)
2080 2093
2081 2094 #Most-Accurate Second Estimation
2082 2095 phi1_aux = arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]]
2083 2096 phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
2084 2097 #Direction Cosine 1
2085 2098 cosdir1 = -(phi1_aux + ang_aux)/(2*numpy.pi*4.5)
2086 2099
2087 2100 #Searching the correct Direction Cosine
2088 2101 cosdir0_aux = cosdir0[:,i]
2089 2102 cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
2090 2103 #Minimum Distance
2091 2104 cosDiff = (cosdir1 - cosdir0_aux)**2
2092 2105 indcos = cosDiff.argmin(axis = 1)
2093 2106 #Saving Value obtained
2094 2107 cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
2095 2108
2096 2109 return cosdir0, cosdir
2097 2110
2098 2111 def __calculateAOA(self, cosdir, azimuth):
2099 2112 cosdirX = cosdir[:,0]
2100 2113 cosdirY = cosdir[:,1]
2101 2114
2102 2115 zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
2103 2116 azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth #0 deg north, 90 deg east
2104 2117 angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
2105 2118
2106 2119 return angles
2107 2120
2108 2121 def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
2109 2122
2110 2123 Ramb = 375 #Ramb = c/(2*PRF)
2111 2124 Re = 6371 #Earth Radius
2112 2125 heights = numpy.zeros(Ranges.shape)
2113 2126
2114 2127 R_aux = numpy.array([0,1,2])*Ramb
2115 2128 R_aux = R_aux.reshape(1,R_aux.size)
2116 2129
2117 2130 Ranges = Ranges.reshape(Ranges.size,1)
2118 2131
2119 2132 Ri = Ranges + R_aux
2120 2133 hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re
2121 2134
2122 2135 #Check if there is a height between 70 and 110 km
2123 2136 h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
2124 2137 ind_h = numpy.where(h_bool == 1)[0]
2125 2138
2126 2139 hCorr = hi[ind_h, :]
2127 2140 ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
2128 2141
2129 2142 hCorr = hi[ind_hCorr]
2130 2143 heights[ind_h] = hCorr
2131 2144
2132 2145 #Setting Error
2133 2146 #Number 13: Height unresolvable echo: not valid height within 70 to 110 km
2134 2147 #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
2135 2148
2136 2149 indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
2137 2150 error[indInvalid2] = 14
2138 2151 indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
2139 2152 error[indInvalid1] = 13
2140 2153
2141 return heights, error
2142
2143
2144 No newline at end of file
2154 return heights, error No newline at end of file
General Comments 0
You need to be logged in to leave comments. Login now