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