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