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