##// END OF EJS Templates
schain
RSS
First Spectral Fitting and EW Drifts operative module inside Signal Chain TRUNK
Julio Valdez -
r513:19f921674eb5
parent child
Show More
Show file before | Show file after | Hide comments
@@ -0,0 +1,135
1 # DIAS 19 Y 20 FEB 2014
2 # Comprobacion de Resultados DBS con SA
3
4 import os, sys
5
6 path = os.path.split(os.getcwd())[0]
7 sys.path.append(path)
8
9 from controller import *
10
11 desc = "DBS Experiment Test"
12 filename = "DBStest.xml"
13
14 controllerObj = Project()
15
16 controllerObj.setup(id = '191', name='test01', description=desc)
17
18 #Experimentos
19
20 path = '/host/Jicamarca/EW_Drifts/d2012248'
21 pathFigure = '/home/propietario/workspace/Graficos/drifts'
22
23
24 path = "/home/soporte/Data/drifts"
25 pathFigure = '/home/soporte/workspace/Graficos/drifts/prueba'
26
27 xmin = 11.75
28 xmax = 14.75
29 #------------------------------------------------------------------------------------------------
30 readUnitConfObj = controllerObj.addReadUnit(datatype='VoltageReader',
31 path=path,
32 startDate='2012/01/01',
33 endDate='2012/12/31',
34 startTime='00:00:00',
35 endTime='23:59:59',
36 online=0,
37 walk=1)
38
39 opObj11 = readUnitConfObj.addOperation(name='printNumberOfBlock')
40
41 #--------------------------------------------------------------------------------------------------
42
43 procUnitConfObj0 = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())
44
45 opObj11 = procUnitConfObj0.addOperation(name='ProfileSelector', optype='other')
46 opObj11.addParameter(name='profileRangeList', value='0,127', format='intlist')
47
48 opObj11 = procUnitConfObj0.addOperation(name='filterByHeights')
49 opObj11.addParameter(name='window', value='3', format='int')
50
51 opObj11 = procUnitConfObj0.addOperation(name='Decoder', optype='other')
52 # opObj11.addParameter(name='code', value='1,-1', format='floatlist')
53 # opObj11.addParameter(name='nCode', value='2', format='int')
54 # opObj11.addParameter(name='nBaud', value='1', format='int')
55
56 procUnitConfObj1 = controllerObj.addProcUnit(datatype='SpectraProc', inputId=procUnitConfObj0.getId())
57 procUnitConfObj1.addParameter(name='nFFTPoints', value='128', format='int')
58 procUnitConfObj1.addParameter(name='nProfiles', value='128', format='int')
59 procUnitConfObj1.addParameter(name='pairsList', value='(0,1),(2,3)', format='pairsList')#,(2,3)
60
61 opObj11 = procUnitConfObj1.addOperation(name='selectHeights')
62 # # opObj11.addParameter(name='minHei', value='320.0', format='float')
63 # # opObj11.addParameter(name='maxHei', value='350.0', format='float')
64 opObj11.addParameter(name='minHei', value='200.0', format='float')
65 opObj11.addParameter(name='maxHei', value='600.0', format='float')
66
67 opObj11 = procUnitConfObj1.addOperation(name='selectChannels')
68 opObj11.addParameter(name='channelList', value='0,1,2,3', format='intlist')
69
70 opObj11 = procUnitConfObj1.addOperation(name='IncohInt', optype='other')
71 opObj11.addParameter(name='timeInterval', value='300.0', format='float')
72
73 opObj13 = procUnitConfObj1.addOperation(name='removeDC')
74
75 # opObj14 = procUnitConfObj1.addOperation(name='SpectraPlot', optype='other')
76 # opObj14.addParameter(name='id', value='1', format='int')
77 # # opObj14.addParameter(name='wintitle', value='Con interf', format='str')
78 # opObj14.addParameter(name='save', value='1', format='bool')
79 # opObj14.addParameter(name='figpath', value=pathFigure, format='str')
80 # # opObj14.addParameter(name='zmin', value='5', format='int')
81 # opObj14.addParameter(name='zmax', value='30', format='int')
82 #
83 # opObj12 = procUnitConfObj1.addOperation(name='RTIPlot', optype='other')
84 # opObj12.addParameter(name='id', value='2', format='int')
85 # opObj12.addParameter(name='wintitle', value='RTI Plot', format='str')
86 # opObj12.addParameter(name='save', value='1', format='bool')
87 # opObj12.addParameter(name='figpath', value = pathFigure, format='str')
88 # opObj12.addParameter(name='xmin', value=xmin, format='float')
89 # opObj12.addParameter(name='xmax', value=xmax, format='float')
90 # # opObj12.addParameter(name='zmin', value='5', format='int')
91 # opObj12.addParameter(name='zmax', value='30', format='int')
92
93 #--------------------------------------------------------------------------------------------------
94
95 procUnitConfObj2 = controllerObj.addProcUnit(datatype='ParametersProc', inputId=procUnitConfObj1.getId())
96 opObj20 = procUnitConfObj2.addOperation(name='SpectralFitting')
97 opObj20.addParameter(name='path', value='/home/soporte/workspace/RemoteSystemsTempFiles', format='str')
98 opObj20.addParameter(name='file', value='modelSpectralFitting', format='str')
99 opObj20.addParameter(name='groupList', value='(0,1),(2,3)',format='multiList')
100
101 opObj11 = procUnitConfObj2.addOperation(name='SpectralFittingPlot', optype='other')
102 opObj11.addParameter(name='id', value='3', format='int')
103 opObj11.addParameter(name='wintitle', value='DopplerPlot', format='str')
104 opObj11.addParameter(name='cutHeight', value='350', format='int')
105 opObj11.addParameter(name='fit', value='1', format='int')#1--True/include fit
106 opObj11.addParameter(name='save', value='1', format='bool')
107 opObj11.addParameter(name='figpath', value = pathFigure, format='str')
108
109 opObj11 = procUnitConfObj2.addOperation(name='EWDriftsEstimation', optype='other')
110 opObj11.addParameter(name='zenith', value='-3.80208,3.10658', format='floatlist')
111 opObj11.addParameter(name='zenithCorrection', value='0.183201', format='float')
112
113 opObj23 = procUnitConfObj2.addOperation(name='EWDriftsPlot', optype='other')
114 opObj23.addParameter(name='id', value='4', format='int')
115 opObj23.addParameter(name='wintitle', value='EW Drifts', format='str')
116 opObj23.addParameter(name='save', value='1', format='bool')
117 opObj23.addParameter(name='figpath', value = pathFigure, format='str')
118 opObj23.addParameter(name='zminZonal', value='-150', format='int')
119 opObj23.addParameter(name='zmaxZonal', value='150', format='int')
120 opObj23.addParameter(name='zminVertical', value='-30', format='float')
121 opObj23.addParameter(name='zmaxVertical', value='30', format='float')
122 opObj23.addParameter(name='SNR_1', value='1', format='bool')
123 opObj23.addParameter(name='SNRmax', value='5', format='int')
124 # opObj23.addParameter(name='SNRthresh', value='-50', format='float')
125 opObj23.addParameter(name='xmin', value=xmin, format='float')
126 opObj23.addParameter(name='xmax', value=xmax, format='float')
127 #--------------------------------------------------------------------------------------------------
128 print "Escribiendo el archivo XML"
129 controllerObj.writeXml(filename)
130 print "Leyendo el archivo XML"
131 controllerObj.readXml(filename)
132
133 controllerObj.createObjects()
134 controllerObj.connectObjects()
135 controllerObj.run() No newline at end of file
Show file before | Show file after | Hide comments
@@ -1,967 +1,983
1 '''
1 '''
2
2
3 $Author: murco $
3 $Author: murco $
4 $Id: JROData.py 173 2012-11-20 15:06:21Z murco $
4 $Id: JROData.py 173 2012-11-20 15:06:21Z murco $
5 '''
5 '''
6
6
7 import copy
7 import copy
8 import numpy
8 import numpy
9 import datetime
9 import datetime
10
10
11 from jroheaderIO import SystemHeader, RadarControllerHeader
11 from jroheaderIO import SystemHeader, RadarControllerHeader
12
12
13 def getNumpyDtype(dataTypeCode):
13 def getNumpyDtype(dataTypeCode):
14
14
15 if dataTypeCode == 0:
15 if dataTypeCode == 0:
16 numpyDtype = numpy.dtype([('real','<i1'),('imag','<i1')])
16 numpyDtype = numpy.dtype([('real','<i1'),('imag','<i1')])
17 elif dataTypeCode == 1:
17 elif dataTypeCode == 1:
18 numpyDtype = numpy.dtype([('real','<i2'),('imag','<i2')])
18 numpyDtype = numpy.dtype([('real','<i2'),('imag','<i2')])
19 elif dataTypeCode == 2:
19 elif dataTypeCode == 2:
20 numpyDtype = numpy.dtype([('real','<i4'),('imag','<i4')])
20 numpyDtype = numpy.dtype([('real','<i4'),('imag','<i4')])
21 elif dataTypeCode == 3:
21 elif dataTypeCode == 3:
22 numpyDtype = numpy.dtype([('real','<i8'),('imag','<i8')])
22 numpyDtype = numpy.dtype([('real','<i8'),('imag','<i8')])
23 elif dataTypeCode == 4:
23 elif dataTypeCode == 4:
24 numpyDtype = numpy.dtype([('real','<f4'),('imag','<f4')])
24 numpyDtype = numpy.dtype([('real','<f4'),('imag','<f4')])
25 elif dataTypeCode == 5:
25 elif dataTypeCode == 5:
26 numpyDtype = numpy.dtype([('real','<f8'),('imag','<f8')])
26 numpyDtype = numpy.dtype([('real','<f8'),('imag','<f8')])
27 else:
27 else:
28 raise ValueError, 'dataTypeCode was not defined'
28 raise ValueError, 'dataTypeCode was not defined'
29
29
30 return numpyDtype
30 return numpyDtype
31
31
32 def getDataTypeCode(numpyDtype):
32 def getDataTypeCode(numpyDtype):
33
33
34 if numpyDtype == numpy.dtype([('real','<i1'),('imag','<i1')]):
34 if numpyDtype == numpy.dtype([('real','<i1'),('imag','<i1')]):
35 datatype = 0
35 datatype = 0
36 elif numpyDtype == numpy.dtype([('real','<i2'),('imag','<i2')]):
36 elif numpyDtype == numpy.dtype([('real','<i2'),('imag','<i2')]):
37 datatype = 1
37 datatype = 1
38 elif numpyDtype == numpy.dtype([('real','<i4'),('imag','<i4')]):
38 elif numpyDtype == numpy.dtype([('real','<i4'),('imag','<i4')]):
39 datatype = 2
39 datatype = 2
40 elif numpyDtype == numpy.dtype([('real','<i8'),('imag','<i8')]):
40 elif numpyDtype == numpy.dtype([('real','<i8'),('imag','<i8')]):
41 datatype = 3
41 datatype = 3
42 elif numpyDtype == numpy.dtype([('real','<f4'),('imag','<f4')]):
42 elif numpyDtype == numpy.dtype([('real','<f4'),('imag','<f4')]):
43 datatype = 4
43 datatype = 4
44 elif numpyDtype == numpy.dtype([('real','<f8'),('imag','<f8')]):
44 elif numpyDtype == numpy.dtype([('real','<f8'),('imag','<f8')]):
45 datatype = 5
45 datatype = 5
46 else:
46 else:
47 datatype = None
47 datatype = None
48
48
49 return datatype
49 return datatype
50
50
51 def hildebrand_sekhon(data, navg):
51 def hildebrand_sekhon(data, navg):
52
52
53 data = data.copy()
53 data = data.copy()
54
54
55 sortdata = numpy.sort(data,axis=None)
55 sortdata = numpy.sort(data,axis=None)
56 lenOfData = len(sortdata)
56 lenOfData = len(sortdata)
57 nums_min = lenOfData/10
57 nums_min = lenOfData/10
58
58
59 if (lenOfData/10) > 2:
59 if (lenOfData/10) > 2:
60 nums_min = lenOfData/10
60 nums_min = lenOfData/10
61 else:
61 else:
62 nums_min = 2
62 nums_min = 2
63
63
64 sump = 0.
64 sump = 0.
65
65
66 sumq = 0.
66 sumq = 0.
67
67
68 j = 0
68 j = 0
69
69
70 cont = 1
70 cont = 1
71
71
72 while((cont==1)and(j<lenOfData)):
72 while((cont==1)and(j<lenOfData)):
73
73
74 sump += sortdata[j]
74 sump += sortdata[j]
75
75
76 sumq += sortdata[j]**2
76 sumq += sortdata[j]**2
77
77
78 j += 1
78 j += 1
79
79
80 if j > nums_min:
80 if j > nums_min:
81 rtest = float(j)/(j-1) + 1.0/navg
81 rtest = float(j)/(j-1) + 1.0/navg
82 if ((sumq*j) > (rtest*sump**2)):
82 if ((sumq*j) > (rtest*sump**2)):
83 j = j - 1
83 j = j - 1
84 sump = sump - sortdata[j]
84 sump = sump - sortdata[j]
85 sumq = sumq - sortdata[j]**2
85 sumq = sumq - sortdata[j]**2
86 cont = 0
86 cont = 0
87
87
88 lnoise = sump /j
88 lnoise = sump /j
89 stdv = numpy.sqrt((sumq - lnoise**2)/(j - 1))
89 stdv = numpy.sqrt((sumq - lnoise**2)/(j - 1))
90 return lnoise
90 return lnoise
91
91
92 class Beam:
92 class Beam:
93 def __init__(self):
93 def __init__(self):
94 self.codeList = []
94 self.codeList = []
95 self.azimuthList = []
95 self.azimuthList = []
96 self.zenithList = []
96 self.zenithList = []
97
97
98 class GenericData(object):
98 class GenericData(object):
99
99
100 flagNoData = True
100 flagNoData = True
101
101
102 def __init__(self):
102 def __init__(self):
103
103
104 raise ValueError, "This class has not been implemented"
104 raise ValueError, "This class has not been implemented"
105
105
106 def copy(self, inputObj=None):
106 def copy(self, inputObj=None):
107
107
108 if inputObj == None:
108 if inputObj == None:
109 return copy.deepcopy(self)
109 return copy.deepcopy(self)
110
110
111 for key in inputObj.__dict__.keys():
111 for key in inputObj.__dict__.keys():
112 self.__dict__[key] = inputObj.__dict__[key]
112 self.__dict__[key] = inputObj.__dict__[key]
113
113
114 def deepcopy(self):
114 def deepcopy(self):
115
115
116 return copy.deepcopy(self)
116 return copy.deepcopy(self)
117
117
118 def isEmpty(self):
118 def isEmpty(self):
119
119
120 return self.flagNoData
120 return self.flagNoData
121
121
122 class JROData(GenericData):
122 class JROData(GenericData):
123
123
124 # m_BasicHeader = BasicHeader()
124 # m_BasicHeader = BasicHeader()
125 # m_ProcessingHeader = ProcessingHeader()
125 # m_ProcessingHeader = ProcessingHeader()
126
126
127 systemHeaderObj = SystemHeader()
127 systemHeaderObj = SystemHeader()
128
128
129 radarControllerHeaderObj = RadarControllerHeader()
129 radarControllerHeaderObj = RadarControllerHeader()
130
130
131 # data = None
131 # data = None
132
132
133 type = None
133 type = None
134
134
135 datatype = None #dtype but in string
135 datatype = None #dtype but in string
136
136
137 # dtype = None
137 # dtype = None
138
138
139 # nChannels = None
139 # nChannels = None
140
140
141 # nHeights = None
141 # nHeights = None
142
142
143 nProfiles = None
143 nProfiles = None
144
144
145 heightList = None
145 heightList = None
146
146
147 channelList = None
147 channelList = None
148
148
149 flagTimeBlock = False
149 flagTimeBlock = False
150
150
151 useLocalTime = False
151 useLocalTime = False
152
152
153 utctime = None
153 utctime = None
154
154
155 timeZone = None
155 timeZone = None
156
156
157 dstFlag = None
157 dstFlag = None
158
158
159 errorCount = None
159 errorCount = None
160
160
161 blocksize = None
161 blocksize = None
162
162
163 nCode = None
163 nCode = None
164
164
165 nBaud = None
165 nBaud = None
166
166
167 code = None
167 code = None
168
168
169 flagDecodeData = False #asumo q la data no esta decodificada
169 flagDecodeData = False #asumo q la data no esta decodificada
170
170
171 flagDeflipData = False #asumo q la data no esta sin flip
171 flagDeflipData = False #asumo q la data no esta sin flip
172
172
173 flagShiftFFT = False
173 flagShiftFFT = False
174
174
175 # ippSeconds = None
175 # ippSeconds = None
176
176
177 timeInterval = None
177 timeInterval = None
178
178
179 nCohInt = None
179 nCohInt = None
180
180
181 noise = None
181 noise = None
182
182
183 windowOfFilter = 1
183 windowOfFilter = 1
184
184
185 #Speed of ligth
185 #Speed of ligth
186 C = 3e8
186 C = 3e8
187
187
188 frequency = 49.92e6
188 frequency = 49.92e6
189
189
190 realtime = False
190 realtime = False
191
191
192 beacon_heiIndexList = None
192 beacon_heiIndexList = None
193
193
194 last_block = None
194 last_block = None
195
195
196 blocknow = None
196 blocknow = None
197
197
198 azimuth = None
198 azimuth = None
199
199
200 zenith = None
200 zenith = None
201
201
202 beam = Beam()
202 beam = Beam()
203
203
204 def __init__(self):
204 def __init__(self):
205
205
206 raise ValueError, "This class has not been implemented"
206 raise ValueError, "This class has not been implemented"
207
207
208 def getNoise(self):
208 def getNoise(self):
209
209
210 raise ValueError, "Not implemented"
210 raise ValueError, "Not implemented"
211
211
212 def getNChannels(self):
212 def getNChannels(self):
213
213
214 return len(self.channelList)
214 return len(self.channelList)
215
215
216 def getChannelIndexList(self):
216 def getChannelIndexList(self):
217
217
218 return range(self.nChannels)
218 return range(self.nChannels)
219
219
220 def getNHeights(self):
220 def getNHeights(self):
221
221
222 return len(self.heightList)
222 return len(self.heightList)
223
223
224 def getHeiRange(self, extrapoints=0):
224 def getHeiRange(self, extrapoints=0):
225
225
226 heis = self.heightList
226 heis = self.heightList
227 # deltah = self.heightList[1] - self.heightList[0]
227 # deltah = self.heightList[1] - self.heightList[0]
228 #
228 #
229 # heis.append(self.heightList[-1])
229 # heis.append(self.heightList[-1])
230
230
231 return heis
231 return heis
232
232
233 def getltctime(self):
233 def getltctime(self):
234
234
235 if self.useLocalTime:
235 if self.useLocalTime:
236 return self.utctime - self.timeZone*60
236 return self.utctime - self.timeZone*60
237
237
238 return self.utctime
238 return self.utctime
239
239
240 def getDatatime(self):
240 def getDatatime(self):
241
241
242 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
242 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
243 return datatimeValue
243 return datatimeValue
244
244
245 def getTimeRange(self):
245 def getTimeRange(self):
246
246
247 datatime = []
247 datatime = []
248
248
249 datatime.append(self.ltctime)
249 datatime.append(self.ltctime)
250 datatime.append(self.ltctime + self.timeInterval)
250 datatime.append(self.ltctime + self.timeInterval)
251
251
252 datatime = numpy.array(datatime)
252 datatime = numpy.array(datatime)
253
253
254 return datatime
254 return datatime
255
255
256 def getFmax(self):
256 def getFmax(self):
257
257
258 PRF = 1./(self.ippSeconds * self.nCohInt)
258 PRF = 1./(self.ippSeconds * self.nCohInt)
259
259
260 fmax = PRF/2.
260 fmax = PRF/2.
261
261
262 return fmax
262 return fmax
263
263
264 def getVmax(self):
264 def getVmax(self):
265
265
266 _lambda = self.C/self.frequency
266 _lambda = self.C/self.frequency
267
267
268 vmax = self.getFmax() * _lambda
268 vmax = self.getFmax() * _lambda
269
269
270 return vmax
270 return vmax
271
271
272 def get_ippSeconds(self):
272 def get_ippSeconds(self):
273 '''
273 '''
274 '''
274 '''
275 return self.radarControllerHeaderObj.ippSeconds
275 return self.radarControllerHeaderObj.ippSeconds
276
276
277 def set_ippSeconds(self, ippSeconds):
277 def set_ippSeconds(self, ippSeconds):
278 '''
278 '''
279 '''
279 '''
280
280
281 self.radarControllerHeaderObj.ippSeconds = ippSeconds
281 self.radarControllerHeaderObj.ippSeconds = ippSeconds
282
282
283 return
283 return
284
284
285 def get_dtype(self):
285 def get_dtype(self):
286 '''
286 '''
287 '''
287 '''
288 return getNumpyDtype(self.datatype)
288 return getNumpyDtype(self.datatype)
289
289
290 def set_dtype(self, numpyDtype):
290 def set_dtype(self, numpyDtype):
291 '''
291 '''
292 '''
292 '''
293
293
294 self.datatype = getDataTypeCode(numpyDtype)
294 self.datatype = getDataTypeCode(numpyDtype)
295
295
296 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
296 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
297 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
297 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
298 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
298 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
299 #noise = property(getNoise, "I'm the 'nHeights' property.")
299 #noise = property(getNoise, "I'm the 'nHeights' property.")
300 datatime = property(getDatatime, "I'm the 'datatime' property")
300 datatime = property(getDatatime, "I'm the 'datatime' property")
301 ltctime = property(getltctime, "I'm the 'ltctime' property")
301 ltctime = property(getltctime, "I'm the 'ltctime' property")
302 ippSeconds = property(get_ippSeconds, set_ippSeconds)
302 ippSeconds = property(get_ippSeconds, set_ippSeconds)
303 dtype = property(get_dtype, set_dtype)
303 dtype = property(get_dtype, set_dtype)
304
304
305 class Voltage(JROData):
305 class Voltage(JROData):
306
306
307 #data es un numpy array de 2 dmensiones (canales, alturas)
307 #data es un numpy array de 2 dmensiones (canales, alturas)
308 data = None
308 data = None
309
309
310 def __init__(self):
310 def __init__(self):
311 '''
311 '''
312 Constructor
312 Constructor
313 '''
313 '''
314
314
315 self.radarControllerHeaderObj = RadarControllerHeader()
315 self.radarControllerHeaderObj = RadarControllerHeader()
316
316
317 self.systemHeaderObj = SystemHeader()
317 self.systemHeaderObj = SystemHeader()
318
318
319 self.type = "Voltage"
319 self.type = "Voltage"
320
320
321 self.data = None
321 self.data = None
322
322
323 # self.dtype = None
323 # self.dtype = None
324
324
325 # self.nChannels = 0
325 # self.nChannels = 0
326
326
327 # self.nHeights = 0
327 # self.nHeights = 0
328
328
329 self.nProfiles = None
329 self.nProfiles = None
330
330
331 self.heightList = None
331 self.heightList = None
332
332
333 self.channelList = None
333 self.channelList = None
334
334
335 # self.channelIndexList = None
335 # self.channelIndexList = None
336
336
337 self.flagNoData = True
337 self.flagNoData = True
338
338
339 self.flagTimeBlock = False
339 self.flagTimeBlock = False
340
340
341 self.utctime = None
341 self.utctime = None
342
342
343 self.timeZone = None
343 self.timeZone = None
344
344
345 self.dstFlag = None
345 self.dstFlag = None
346
346
347 self.errorCount = None
347 self.errorCount = None
348
348
349 self.nCohInt = None
349 self.nCohInt = None
350
350
351 self.blocksize = None
351 self.blocksize = None
352
352
353 self.flagDecodeData = False #asumo q la data no esta decodificada
353 self.flagDecodeData = False #asumo q la data no esta decodificada
354
354
355 self.flagDeflipData = False #asumo q la data no esta sin flip
355 self.flagDeflipData = False #asumo q la data no esta sin flip
356
356
357 self.flagShiftFFT = False
357 self.flagShiftFFT = False
358
358
359
359
360 def getNoisebyHildebrand(self):
360 def getNoisebyHildebrand(self):
361 """
361 """
362 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
362 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
363
363
364 Return:
364 Return:
365 noiselevel
365 noiselevel
366 """
366 """
367
367
368 for channel in range(self.nChannels):
368 for channel in range(self.nChannels):
369 daux = self.data_spc[channel,:,:]
369 daux = self.data_spc[channel,:,:]
370 self.noise[channel] = hildebrand_sekhon(daux, self.nCohInt)
370 self.noise[channel] = hildebrand_sekhon(daux, self.nCohInt)
371
371
372 return self.noise
372 return self.noise
373
373
374 def getNoise(self, type = 1):
374 def getNoise(self, type = 1):
375
375
376 self.noise = numpy.zeros(self.nChannels)
376 self.noise = numpy.zeros(self.nChannels)
377
377
378 if type == 1:
378 if type == 1:
379 noise = self.getNoisebyHildebrand()
379 noise = self.getNoisebyHildebrand()
380
380
381 return 10*numpy.log10(noise)
381 return 10*numpy.log10(noise)
382
382
383 noise = property(getNoise, "I'm the 'nHeights' property.")
383 noise = property(getNoise, "I'm the 'nHeights' property.")
384
384
385 class Spectra(JROData):
385 class Spectra(JROData):
386
386
387 #data es un numpy array de 2 dmensiones (canales, perfiles, alturas)
387 #data es un numpy array de 2 dmensiones (canales, perfiles, alturas)
388 data_spc = None
388 data_spc = None
389
389
390 #data es un numpy array de 2 dmensiones (canales, pares, alturas)
390 #data es un numpy array de 2 dmensiones (canales, pares, alturas)
391 data_cspc = None
391 data_cspc = None
392
392
393 #data es un numpy array de 2 dmensiones (canales, alturas)
393 #data es un numpy array de 2 dmensiones (canales, alturas)
394 data_dc = None
394 data_dc = None
395
395
396 nFFTPoints = None
396 nFFTPoints = None
397
397
398 # nPairs = None
398 # nPairs = None
399
399
400 pairsList = None
400 pairsList = None
401
401
402 nIncohInt = None
402 nIncohInt = None
403
403
404 wavelength = None #Necesario para cacular el rango de velocidad desde la frecuencia
404 wavelength = None #Necesario para cacular el rango de velocidad desde la frecuencia
405
405
406 nCohInt = None #se requiere para determinar el valor de timeInterval
406 nCohInt = None #se requiere para determinar el valor de timeInterval
407
407
408 ippFactor = None
408 ippFactor = None
409
409
410 def __init__(self):
410 def __init__(self):
411 '''
411 '''
412 Constructor
412 Constructor
413 '''
413 '''
414
414
415 self.radarControllerHeaderObj = RadarControllerHeader()
415 self.radarControllerHeaderObj = RadarControllerHeader()
416
416
417 self.systemHeaderObj = SystemHeader()
417 self.systemHeaderObj = SystemHeader()
418
418
419 self.type = "Spectra"
419 self.type = "Spectra"
420
420
421 # self.data = None
421 # self.data = None
422
422
423 # self.dtype = None
423 # self.dtype = None
424
424
425 # self.nChannels = 0
425 # self.nChannels = 0
426
426
427 # self.nHeights = 0
427 # self.nHeights = 0
428
428
429 self.nProfiles = None
429 self.nProfiles = None
430
430
431 self.heightList = None
431 self.heightList = None
432
432
433 self.channelList = None
433 self.channelList = None
434
434
435 # self.channelIndexList = None
435 # self.channelIndexList = None
436
436
437 self.pairsList = None
437 self.pairsList = None
438
438
439 self.flagNoData = True
439 self.flagNoData = True
440
440
441 self.flagTimeBlock = False
441 self.flagTimeBlock = False
442
442
443 self.utctime = None
443 self.utctime = None
444
444
445 self.nCohInt = None
445 self.nCohInt = None
446
446
447 self.nIncohInt = None
447 self.nIncohInt = None
448
448
449 self.blocksize = None
449 self.blocksize = None
450
450
451 self.nFFTPoints = None
451 self.nFFTPoints = None
452
452
453 self.wavelength = None
453 self.wavelength = None
454
454
455 self.flagDecodeData = False #asumo q la data no esta decodificada
455 self.flagDecodeData = False #asumo q la data no esta decodificada
456
456
457 self.flagDeflipData = False #asumo q la data no esta sin flip
457 self.flagDeflipData = False #asumo q la data no esta sin flip
458
458
459 self.flagShiftFFT = False
459 self.flagShiftFFT = False
460
460
461 self.ippFactor = 1
461 self.ippFactor = 1
462
462
463 #self.noise = None
463 #self.noise = None
464
464
465 self.beacon_heiIndexList = []
465 self.beacon_heiIndexList = []
466
466
467 self.noise_estimation = None
467 self.noise_estimation = None
468
468
469
469
470 def getNoisebyHildebrand(self):
470 def getNoisebyHildebrand(self):
471 """
471 """
472 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
472 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
473
473
474 Return:
474 Return:
475 noiselevel
475 noiselevel
476 """
476 """
477
477
478 noise = numpy.zeros(self.nChannels)
478 noise = numpy.zeros(self.nChannels)
479 for channel in range(self.nChannels):
479 for channel in range(self.nChannels):
480 daux = self.data_spc[channel,:,:]
480 daux = self.data_spc[channel,:,:]
481 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
481 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
482
482
483 return noise
483 return noise
484
484
485 def getNoise(self):
485 def getNoise(self):
486 if self.noise_estimation != None:
486 if self.noise_estimation != None:
487 return self.noise_estimation #this was estimated by getNoise Operation defined in jroproc_spectra.py
487 return self.noise_estimation #this was estimated by getNoise Operation defined in jroproc_spectra.py
488 else:
488 else:
489 noise = self.getNoisebyHildebrand()
489 noise = self.getNoisebyHildebrand()
490 return noise
490 return noise
491
491
492
492
493 def getFreqRange(self, extrapoints=0):
493 def getFreqRange(self, extrapoints=0):
494
494
495 deltafreq = self.getFmax() / (self.nFFTPoints*self.ippFactor)
495 deltafreq = self.getFmax() / (self.nFFTPoints*self.ippFactor)
496 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
496 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
497
497
498 return freqrange
498 return freqrange
499
499
500 def getVelRange(self, extrapoints=0):
500 def getVelRange(self, extrapoints=0):
501
501
502 deltav = self.getVmax() / (self.nFFTPoints*self.ippFactor)
502 deltav = self.getVmax() / (self.nFFTPoints*self.ippFactor)
503 velrange = deltav*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltav/2
503 velrange = deltav*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltav/2
504
504
505 return velrange
505 return velrange
506
506
507 def getNPairs(self):
507 def getNPairs(self):
508
508
509 return len(self.pairsList)
509 return len(self.pairsList)
510
510
511 def getPairsIndexList(self):
511 def getPairsIndexList(self):
512
512
513 return range(self.nPairs)
513 return range(self.nPairs)
514
514
515 def getNormFactor(self):
515 def getNormFactor(self):
516 pwcode = 1
516 pwcode = 1
517 if self.flagDecodeData:
517 if self.flagDecodeData:
518 pwcode = numpy.sum(self.code[0]**2)
518 pwcode = numpy.sum(self.code[0]**2)
519 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
519 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
520 normFactor = self.nProfiles*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
520 normFactor = self.nProfiles*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
521
521
522 return normFactor
522 return normFactor
523
523
524 def getFlagCspc(self):
524 def getFlagCspc(self):
525
525
526 if self.data_cspc == None:
526 if self.data_cspc == None:
527 return True
527 return True
528
528
529 return False
529 return False
530
530
531 def getFlagDc(self):
531 def getFlagDc(self):
532
532
533 if self.data_dc == None:
533 if self.data_dc == None:
534 return True
534 return True
535
535
536 return False
536 return False
537
537
538 nPairs = property(getNPairs, "I'm the 'nPairs' property.")
538 nPairs = property(getNPairs, "I'm the 'nPairs' property.")
539 pairsIndexList = property(getPairsIndexList, "I'm the 'pairsIndexList' property.")
539 pairsIndexList = property(getPairsIndexList, "I'm the 'pairsIndexList' property.")
540 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
540 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
541 flag_cspc = property(getFlagCspc)
541 flag_cspc = property(getFlagCspc)
542 flag_dc = property(getFlagDc)
542 flag_dc = property(getFlagDc)
543 noise = property(getNoise, "I'm the 'nHeights' property.")
543 noise = property(getNoise, "I'm the 'nHeights' property.")
544
544
545 class SpectraHeis(Spectra):
545 class SpectraHeis(Spectra):
546
546
547 data_spc = None
547 data_spc = None
548
548
549 data_cspc = None
549 data_cspc = None
550
550
551 data_dc = None
551 data_dc = None
552
552
553 nFFTPoints = None
553 nFFTPoints = None
554
554
555 # nPairs = None
555 # nPairs = None
556
556
557 pairsList = None
557 pairsList = None
558
558
559 nIncohInt = None
559 nIncohInt = None
560
560
561 def __init__(self):
561 def __init__(self):
562
562
563 self.radarControllerHeaderObj = RadarControllerHeader()
563 self.radarControllerHeaderObj = RadarControllerHeader()
564
564
565 self.systemHeaderObj = SystemHeader()
565 self.systemHeaderObj = SystemHeader()
566
566
567 self.type = "SpectraHeis"
567 self.type = "SpectraHeis"
568
568
569 # self.dtype = None
569 # self.dtype = None
570
570
571 # self.nChannels = 0
571 # self.nChannels = 0
572
572
573 # self.nHeights = 0
573 # self.nHeights = 0
574
574
575 self.nProfiles = None
575 self.nProfiles = None
576
576
577 self.heightList = None
577 self.heightList = None
578
578
579 self.channelList = None
579 self.channelList = None
580
580
581 # self.channelIndexList = None
581 # self.channelIndexList = None
582
582
583 self.flagNoData = True
583 self.flagNoData = True
584
584
585 self.flagTimeBlock = False
585 self.flagTimeBlock = False
586
586
587 # self.nPairs = 0
587 # self.nPairs = 0
588
588
589 self.utctime = None
589 self.utctime = None
590
590
591 self.blocksize = None
591 self.blocksize = None
592
592
593 def getNormFactor(self):
593 def getNormFactor(self):
594 pwcode = 1
594 pwcode = 1
595 if self.flagDecodeData:
595 if self.flagDecodeData:
596 pwcode = numpy.sum(self.code[0]**2)
596 pwcode = numpy.sum(self.code[0]**2)
597
597
598 normFactor = self.nIncohInt*self.nCohInt*pwcode
598 normFactor = self.nIncohInt*self.nCohInt*pwcode
599
599
600 return normFactor
600 return normFactor
601
601
602 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
602 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
603
603
604 class Fits:
604 class Fits:
605
605
606 heightList = None
606 heightList = None
607
607
608 channelList = None
608 channelList = None
609
609
610 flagNoData = True
610 flagNoData = True
611
611
612 flagTimeBlock = False
612 flagTimeBlock = False
613
613
614 useLocalTime = False
614 useLocalTime = False
615
615
616 utctime = None
616 utctime = None
617
617
618 timeZone = None
618 timeZone = None
619
619
620 # ippSeconds = None
620 # ippSeconds = None
621
621
622 timeInterval = None
622 timeInterval = None
623
623
624 nCohInt = None
624 nCohInt = None
625
625
626 nIncohInt = None
626 nIncohInt = None
627
627
628 noise = None
628 noise = None
629
629
630 windowOfFilter = 1
630 windowOfFilter = 1
631
631
632 #Speed of ligth
632 #Speed of ligth
633 C = 3e8
633 C = 3e8
634
634
635 frequency = 49.92e6
635 frequency = 49.92e6
636
636
637 realtime = False
637 realtime = False
638
638
639
639
640 def __init__(self):
640 def __init__(self):
641
641
642 self.type = "Fits"
642 self.type = "Fits"
643
643
644 self.nProfiles = None
644 self.nProfiles = None
645
645
646 self.heightList = None
646 self.heightList = None
647
647
648 self.channelList = None
648 self.channelList = None
649
649
650 # self.channelIndexList = None
650 # self.channelIndexList = None
651
651
652 self.flagNoData = True
652 self.flagNoData = True
653
653
654 self.utctime = None
654 self.utctime = None
655
655
656 self.nCohInt = None
656 self.nCohInt = None
657
657
658 self.nIncohInt = None
658 self.nIncohInt = None
659
659
660 self.useLocalTime = True
660 self.useLocalTime = True
661
661
662 # self.utctime = None
662 # self.utctime = None
663 # self.timeZone = None
663 # self.timeZone = None
664 # self.ltctime = None
664 # self.ltctime = None
665 # self.timeInterval = None
665 # self.timeInterval = None
666 # self.header = None
666 # self.header = None
667 # self.data_header = None
667 # self.data_header = None
668 # self.data = None
668 # self.data = None
669 # self.datatime = None
669 # self.datatime = None
670 # self.flagNoData = False
670 # self.flagNoData = False
671 # self.expName = ''
671 # self.expName = ''
672 # self.nChannels = None
672 # self.nChannels = None
673 # self.nSamples = None
673 # self.nSamples = None
674 # self.dataBlocksPerFile = None
674 # self.dataBlocksPerFile = None
675 # self.comments = ''
675 # self.comments = ''
676 #
676 #
677
677
678
678
679 def getltctime(self):
679 def getltctime(self):
680
680
681 if self.useLocalTime:
681 if self.useLocalTime:
682 return self.utctime - self.timeZone*60
682 return self.utctime - self.timeZone*60
683
683
684 return self.utctime
684 return self.utctime
685
685
686 def getDatatime(self):
686 def getDatatime(self):
687
687
688 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
688 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
689 return datatime
689 return datatime
690
690
691 def getTimeRange(self):
691 def getTimeRange(self):
692
692
693 datatime = []
693 datatime = []
694
694
695 datatime.append(self.ltctime)
695 datatime.append(self.ltctime)
696 datatime.append(self.ltctime + self.timeInterval)
696 datatime.append(self.ltctime + self.timeInterval)
697
697
698 datatime = numpy.array(datatime)
698 datatime = numpy.array(datatime)
699
699
700 return datatime
700 return datatime
701
701
702 def getHeiRange(self):
702 def getHeiRange(self):
703
703
704 heis = self.heightList
704 heis = self.heightList
705
705
706 return heis
706 return heis
707
707
708 def isEmpty(self):
708 def isEmpty(self):
709
709
710 return self.flagNoData
710 return self.flagNoData
711
711
712 def getNHeights(self):
712 def getNHeights(self):
713
713
714 return len(self.heightList)
714 return len(self.heightList)
715
715
716 def getNChannels(self):
716 def getNChannels(self):
717
717
718 return len(self.channelList)
718 return len(self.channelList)
719
719
720 def getChannelIndexList(self):
720 def getChannelIndexList(self):
721
721
722 return range(self.nChannels)
722 return range(self.nChannels)
723
723
724 def getNoise(self, type = 1):
724 def getNoise(self, type = 1):
725
725
726 self.noise = numpy.zeros(self.nChannels)
726 self.noise = numpy.zeros(self.nChannels)
727
727
728 if type == 1:
728 if type == 1:
729 noise = self.getNoisebyHildebrand()
729 noise = self.getNoisebyHildebrand()
730
730
731 if type == 2:
731 if type == 2:
732 noise = self.getNoisebySort()
732 noise = self.getNoisebySort()
733
733
734 if type == 3:
734 if type == 3:
735 noise = self.getNoisebyWindow()
735 noise = self.getNoisebyWindow()
736
736
737 return noise
737 return noise
738
738
739 datatime = property(getDatatime, "I'm the 'datatime' property")
739 datatime = property(getDatatime, "I'm the 'datatime' property")
740 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
740 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
741 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
741 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
742 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
742 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
743 noise = property(getNoise, "I'm the 'nHeights' property.")
743 noise = property(getNoise, "I'm the 'nHeights' property.")
744 datatime = property(getDatatime, "I'm the 'datatime' property")
744 datatime = property(getDatatime, "I'm the 'datatime' property")
745 ltctime = property(getltctime, "I'm the 'ltctime' property")
745 ltctime = property(getltctime, "I'm the 'ltctime' property")
746
746
747 class Correlation(JROData):
747 class Correlation(JROData):
748
748
749 noise = None
749 noise = None
750
750
751 SNR = None
751 SNR = None
752
752
753 pairsAutoCorr = None #Pairs of Autocorrelation
753 pairsAutoCorr = None #Pairs of Autocorrelation
754
754
755 #--------------------------------------------------
755 #--------------------------------------------------
756
756
757 data_corr = None
757 data_corr = None
758
758
759 data_volt = None
759 data_volt = None
760
760
761 lagT = None # each element value is a profileIndex
761 lagT = None # each element value is a profileIndex
762
762
763 lagR = None # each element value is in km
763 lagR = None # each element value is in km
764
764
765 pairsList = None
765 pairsList = None
766
766
767 calculateVelocity = None
767 calculateVelocity = None
768
768
769 nPoints = None
769 nPoints = None
770
770
771 nAvg = None
771 nAvg = None
772
772
773 bufferSize = None
773 bufferSize = None
774
774
775 def __init__(self):
775 def __init__(self):
776 '''
776 '''
777 Constructor
777 Constructor
778 '''
778 '''
779 self.radarControllerHeaderObj = RadarControllerHeader()
779 self.radarControllerHeaderObj = RadarControllerHeader()
780
780
781 self.systemHeaderObj = SystemHeader()
781 self.systemHeaderObj = SystemHeader()
782
782
783 self.type = "Correlation"
783 self.type = "Correlation"
784
784
785 self.data = None
785 self.data = None
786
786
787 self.dtype = None
787 self.dtype = None
788
788
789 self.nProfiles = None
789 self.nProfiles = None
790
790
791 self.heightList = None
791 self.heightList = None
792
792
793 self.channelList = None
793 self.channelList = None
794
794
795 self.flagNoData = True
795 self.flagNoData = True
796
796
797 self.flagTimeBlock = False
797 self.flagTimeBlock = False
798
798
799 self.utctime = None
799 self.utctime = None
800
800
801 self.timeZone = None
801 self.timeZone = None
802
802
803 self.dstFlag = None
803 self.dstFlag = None
804
804
805 self.errorCount = None
805 self.errorCount = None
806
806
807 self.blocksize = None
807 self.blocksize = None
808
808
809 self.flagDecodeData = False #asumo q la data no esta decodificada
809 self.flagDecodeData = False #asumo q la data no esta decodificada
810
810
811 self.flagDeflipData = False #asumo q la data no esta sin flip
811 self.flagDeflipData = False #asumo q la data no esta sin flip
812
812
813 self.pairsList = None
813 self.pairsList = None
814
814
815 self.nPoints = None
815 self.nPoints = None
816
816
817 def getLagTRange(self, extrapoints=0):
817 def getLagTRange(self, extrapoints=0):
818
818
819 lagTRange = self.lagT
819 lagTRange = self.lagT
820 diff = lagTRange[1] - lagTRange[0]
820 diff = lagTRange[1] - lagTRange[0]
821 extra = numpy.arange(1,extrapoints + 1)*diff + lagTRange[-1]
821 extra = numpy.arange(1,extrapoints + 1)*diff + lagTRange[-1]
822 lagTRange = numpy.hstack((lagTRange, extra))
822 lagTRange = numpy.hstack((lagTRange, extra))
823
823
824 return lagTRange
824 return lagTRange
825
825
826 def getLagRRange(self, extrapoints=0):
826 def getLagRRange(self, extrapoints=0):
827
827
828 return self.lagR
828 return self.lagR
829
829
830 def getPairsList(self):
830 def getPairsList(self):
831
831
832 return self.pairsList
832 return self.pairsList
833
833
834 def getCalculateVelocity(self):
834 def getCalculateVelocity(self):
835
835
836 return self.calculateVelocity
836 return self.calculateVelocity
837
837
838 def getNPoints(self):
838 def getNPoints(self):
839
839
840 return self.nPoints
840 return self.nPoints
841
841
842 def getNAvg(self):
842 def getNAvg(self):
843
843
844 return self.nAvg
844 return self.nAvg
845
845
846 def getBufferSize(self):
846 def getBufferSize(self):
847
847
848 return self.bufferSize
848 return self.bufferSize
849
849
850 def getPairsAutoCorr(self):
850 def getPairsAutoCorr(self):
851 pairsList = self.pairsList
851 pairsList = self.pairsList
852 pairsAutoCorr = numpy.zeros(self.nChannels, dtype = 'int')*numpy.nan
852 pairsAutoCorr = numpy.zeros(self.nChannels, dtype = 'int')*numpy.nan
853
853
854 for l in range(len(pairsList)):
854 for l in range(len(pairsList)):
855 firstChannel = pairsList[l][0]
855 firstChannel = pairsList[l][0]
856 secondChannel = pairsList[l][1]
856 secondChannel = pairsList[l][1]
857
857
858 #Obteniendo pares de Autocorrelacion
858 #Obteniendo pares de Autocorrelacion
859 if firstChannel == secondChannel:
859 if firstChannel == secondChannel:
860 pairsAutoCorr[firstChannel] = int(l)
860 pairsAutoCorr[firstChannel] = int(l)
861
861
862 pairsAutoCorr = pairsAutoCorr.astype(int)
862 pairsAutoCorr = pairsAutoCorr.astype(int)
863
863
864 return pairsAutoCorr
864 return pairsAutoCorr
865
865
866 def getNoise(self, mode = 2):
866 def getNoise(self, mode = 2):
867
867
868 indR = numpy.where(self.lagR == 0)[0][0]
868 indR = numpy.where(self.lagR == 0)[0][0]
869 indT = numpy.where(self.lagT == 0)[0][0]
869 indT = numpy.where(self.lagT == 0)[0][0]
870
870
871 jspectra0 = self.data_corr[:,:,indR,:]
871 jspectra0 = self.data_corr[:,:,indR,:]
872 jspectra = copy.copy(jspectra0)
872 jspectra = copy.copy(jspectra0)
873
873
874 num_chan = jspectra.shape[0]
874 num_chan = jspectra.shape[0]
875 num_hei = jspectra.shape[2]
875 num_hei = jspectra.shape[2]
876
876
877 freq_dc = jspectra.shape[1]/2
877 freq_dc = jspectra.shape[1]/2
878 ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
878 ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
879
879
880 if ind_vel[0]<0:
880 if ind_vel[0]<0:
881 ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
881 ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
882
882
883 if mode == 1:
883 if mode == 1:
884 jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
884 jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
885
885
886 if mode == 2:
886 if mode == 2:
887
887
888 vel = numpy.array([-2,-1,1,2])
888 vel = numpy.array([-2,-1,1,2])
889 xx = numpy.zeros([4,4])
889 xx = numpy.zeros([4,4])
890
890
891 for fil in range(4):
891 for fil in range(4):
892 xx[fil,:] = vel[fil]**numpy.asarray(range(4))
892 xx[fil,:] = vel[fil]**numpy.asarray(range(4))
893
893
894 xx_inv = numpy.linalg.inv(xx)
894 xx_inv = numpy.linalg.inv(xx)
895 xx_aux = xx_inv[0,:]
895 xx_aux = xx_inv[0,:]
896
896
897 for ich in range(num_chan):
897 for ich in range(num_chan):
898 yy = jspectra[ich,ind_vel,:]
898 yy = jspectra[ich,ind_vel,:]
899 jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy)
899 jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy)
900
900
901 junkid = jspectra[ich,freq_dc,:]<=0
901 junkid = jspectra[ich,freq_dc,:]<=0
902 cjunkid = sum(junkid)
902 cjunkid = sum(junkid)
903
903
904 if cjunkid.any():
904 if cjunkid.any():
905 jspectra[ich,freq_dc,junkid.nonzero()] = (jspectra[ich,ind_vel[1],junkid] + jspectra[ich,ind_vel[2],junkid])/2
905 jspectra[ich,freq_dc,junkid.nonzero()] = (jspectra[ich,ind_vel[1],junkid] + jspectra[ich,ind_vel[2],junkid])/2
906
906
907 noise = jspectra0[:,freq_dc,:] - jspectra[:,freq_dc,:]
907 noise = jspectra0[:,freq_dc,:] - jspectra[:,freq_dc,:]
908
908
909 return noise
909 return noise
910
910
911 # pairsList = property(getPairsList, "I'm the 'pairsList' property.")
911 # pairsList = property(getPairsList, "I'm the 'pairsList' property.")
912 # nPoints = property(getNPoints, "I'm the 'nPoints' property.")
912 # nPoints = property(getNPoints, "I'm the 'nPoints' property.")
913 calculateVelocity = property(getCalculateVelocity, "I'm the 'calculateVelocity' property.")
913 calculateVelocity = property(getCalculateVelocity, "I'm the 'calculateVelocity' property.")
914 nAvg = property(getNAvg, "I'm the 'nAvg' property.")
914 nAvg = property(getNAvg, "I'm the 'nAvg' property.")
915 bufferSize = property(getBufferSize, "I'm the 'bufferSize' property.")
915 bufferSize = property(getBufferSize, "I'm the 'bufferSize' property.")
916
916
917
917
918 class Parameters(JROData):
918 class Parameters(JROData):
919
919
920 #Information from previous data
921
920 inputUnit = None #Type of data to be processed
922 inputUnit = None #Type of data to be processed
921
923
922 operation = None #Type of operation to parametrize
924 operation = None #Type of operation to parametrize
923
925
926 normFactor = None #Normalization Factor
927
928 groupList = None #List of Pairs, Groups, etc
929
930 #Parameters
931
924 data_param = None #Parameters obtained
932 data_param = None #Parameters obtained
925
933
926 data_pre = None #Data Pre Parametrization
934 data_pre = None #Data Pre Parametrization
927
935
928 heightRange = None #Heights
936 heightRange = None #Heights
929
937
930 abscissaRange = None #Abscissa, can be velocities, lags or time
938 abscissaRange = None #Abscissa, can be velocities, lags or time
931
939
932 noise = None #Noise Potency
940 noise = None #Noise Potency
933
941
934 SNR = None #Signal to Noise Ratio
942 SNR = None #Signal to Noise Ratio
935
943
936 pairsList = None #List of Pairs for Cross correlations or Cross spectrum
937
938 initUtcTime = None #Initial UTC time
944 initUtcTime = None #Initial UTC time
939
945
940 paramInterval = None #Time interval to calculate Parameters in seconds
946 paramInterval = None #Time interval to calculate Parameters in seconds
941
947
942 windsInterval = None #Time interval to calculate Winds in seconds
948 #Fitting
949
950 constants = None
951
952 error = None
953
954 library = None
955
956 #Output signal
957
958 outputInterval = None #Time interval to calculate output signal in seconds
959
960 data_output = None #Out signal
943
961
944 normFactor = None #Normalization Factor
945
962
946 winds = None #Wind estimations
947
963
948 def __init__(self):
964 def __init__(self):
949 '''
965 '''
950 Constructor
966 Constructor
951 '''
967 '''
952 self.radarControllerHeaderObj = RadarControllerHeader()
968 self.radarControllerHeaderObj = RadarControllerHeader()
953
969
954 self.systemHeaderObj = SystemHeader()
970 self.systemHeaderObj = SystemHeader()
955
971
956 self.type = "Parameters"
972 self.type = "Parameters"
957
973
958 def getTimeRange1(self):
974 def getTimeRange1(self):
959
975
960 datatime = []
976 datatime = []
961
977
962 datatime.append(self.initUtcTime)
978 datatime.append(self.initUtcTime)
963 datatime.append(self.initUtcTime + self.windsInterval - 1)
979 datatime.append(self.initUtcTime + self.outputInterval - 1)
964
980
965 datatime = numpy.array(datatime)
981 datatime = numpy.array(datatime)
966
982
967 return datatime
983 return datatime
Show file before | Show file after | Hide comments
@@ -1,774 +1,1178
1 import os
1 import os
2 import datetime
2 import datetime
3 import numpy
3 import numpy
4
4
5 from figure import Figure, isRealtime
5 from figure import Figure, isRealtime
6
6
7 class MomentsPlot(Figure):
7 class MomentsPlot(Figure):
8
8
9 isConfig = None
9 isConfig = None
10 __nsubplots = None
10 __nsubplots = None
11
11
12 WIDTHPROF = None
12 WIDTHPROF = None
13 HEIGHTPROF = None
13 HEIGHTPROF = None
14 PREFIX = 'prm'
14 PREFIX = 'prm'
15
15
16 def __init__(self):
16 def __init__(self):
17
17
18 self.isConfig = False
18 self.isConfig = False
19 self.__nsubplots = 1
19 self.__nsubplots = 1
20
20
21 self.WIDTH = 280
21 self.WIDTH = 280
22 self.HEIGHT = 250
22 self.HEIGHT = 250
23 self.WIDTHPROF = 120
23 self.WIDTHPROF = 120
24 self.HEIGHTPROF = 0
24 self.HEIGHTPROF = 0
25 self.counter_imagwr = 0
25 self.counter_imagwr = 0
26
26
27 self.PLOT_CODE = 1
27 self.PLOT_CODE = 1
28 self.FTP_WEI = None
28 self.FTP_WEI = None
29 self.EXP_CODE = None
29 self.EXP_CODE = None
30 self.SUB_EXP_CODE = None
30 self.SUB_EXP_CODE = None
31 self.PLOT_POS = None
31 self.PLOT_POS = None
32
32
33 def getSubplots(self):
33 def getSubplots(self):
34
34
35 ncol = int(numpy.sqrt(self.nplots)+0.9)
35 ncol = int(numpy.sqrt(self.nplots)+0.9)
36 nrow = int(self.nplots*1./ncol + 0.9)
36 nrow = int(self.nplots*1./ncol + 0.9)
37
37
38 return nrow, ncol
38 return nrow, ncol
39
39
40 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
40 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
41
41
42 self.__showprofile = showprofile
42 self.__showprofile = showprofile
43 self.nplots = nplots
43 self.nplots = nplots
44
44
45 ncolspan = 1
45 ncolspan = 1
46 colspan = 1
46 colspan = 1
47 if showprofile:
47 if showprofile:
48 ncolspan = 3
48 ncolspan = 3
49 colspan = 2
49 colspan = 2
50 self.__nsubplots = 2
50 self.__nsubplots = 2
51
51
52 self.createFigure(id = id,
52 self.createFigure(id = id,
53 wintitle = wintitle,
53 wintitle = wintitle,
54 widthplot = self.WIDTH + self.WIDTHPROF,
54 widthplot = self.WIDTH + self.WIDTHPROF,
55 heightplot = self.HEIGHT + self.HEIGHTPROF,
55 heightplot = self.HEIGHT + self.HEIGHTPROF,
56 show=show)
56 show=show)
57
57
58 nrow, ncol = self.getSubplots()
58 nrow, ncol = self.getSubplots()
59
59
60 counter = 0
60 counter = 0
61 for y in range(nrow):
61 for y in range(nrow):
62 for x in range(ncol):
62 for x in range(ncol):
63
63
64 if counter >= self.nplots:
64 if counter >= self.nplots:
65 break
65 break
66
66
67 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
67 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
68
68
69 if showprofile:
69 if showprofile:
70 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
70 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
71
71
72 counter += 1
72 counter += 1
73
73
74 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=True,
74 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=True,
75 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
75 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
76 save=False, figpath='', figfile=None, show=True, ftp=False, wr_period=1,
76 save=False, figpath='', figfile=None, show=True, ftp=False, wr_period=1,
77 server=None, folder=None, username=None, password=None,
77 server=None, folder=None, username=None, password=None,
78 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False):
78 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False):
79
79
80 """
80 """
81
81
82 Input:
82 Input:
83 dataOut :
83 dataOut :
84 id :
84 id :
85 wintitle :
85 wintitle :
86 channelList :
86 channelList :
87 showProfile :
87 showProfile :
88 xmin : None,
88 xmin : None,
89 xmax : None,
89 xmax : None,
90 ymin : None,
90 ymin : None,
91 ymax : None,
91 ymax : None,
92 zmin : None,
92 zmin : None,
93 zmax : None
93 zmax : None
94 """
94 """
95
95
96 if dataOut.flagNoData:
96 if dataOut.flagNoData:
97 return None
97 return None
98
98
99 if realtime:
99 if realtime:
100 if not(isRealtime(utcdatatime = dataOut.utctime)):
100 if not(isRealtime(utcdatatime = dataOut.utctime)):
101 print 'Skipping this plot function'
101 print 'Skipping this plot function'
102 return
102 return
103
103
104 if channelList == None:
104 if channelList == None:
105 channelIndexList = dataOut.channelIndexList
105 channelIndexList = dataOut.channelIndexList
106 else:
106 else:
107 channelIndexList = []
107 channelIndexList = []
108 for channel in channelList:
108 for channel in channelList:
109 if channel not in dataOut.channelList:
109 if channel not in dataOut.channelList:
110 raise ValueError, "Channel %d is not in dataOut.channelList"
110 raise ValueError, "Channel %d is not in dataOut.channelList"
111 channelIndexList.append(dataOut.channelList.index(channel))
111 channelIndexList.append(dataOut.channelList.index(channel))
112
112
113 factor = dataOut.normFactor
113 factor = dataOut.normFactor
114 x = dataOut.abscissaRange
114 x = dataOut.abscissaRange
115 y = dataOut.heightRange
115 y = dataOut.heightRange
116
116
117 z = dataOut.data_pre[channelIndexList,:,:]/factor
117 z = dataOut.data_pre[channelIndexList,:,:]/factor
118 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
118 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
119 avg = numpy.average(z, axis=1)
119 avg = numpy.average(z, axis=1)
120 noise = dataOut.noise/factor
120 noise = dataOut.noise/factor
121
121
122 zdB = 10*numpy.log10(z)
122 zdB = 10*numpy.log10(z)
123 avgdB = 10*numpy.log10(avg)
123 avgdB = 10*numpy.log10(avg)
124 noisedB = 10*numpy.log10(noise)
124 noisedB = 10*numpy.log10(noise)
125
125
126 #thisDatetime = dataOut.datatime
126 #thisDatetime = dataOut.datatime
127 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
127 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
128 title = wintitle + " Parameters"
128 title = wintitle + " Parameters"
129 xlabel = "Velocity (m/s)"
129 xlabel = "Velocity (m/s)"
130 ylabel = "Range (Km)"
130 ylabel = "Range (Km)"
131
131
132 if not self.isConfig:
132 if not self.isConfig:
133
133
134 nplots = len(channelIndexList)
134 nplots = len(channelIndexList)
135
135
136 self.setup(id=id,
136 self.setup(id=id,
137 nplots=nplots,
137 nplots=nplots,
138 wintitle=wintitle,
138 wintitle=wintitle,
139 showprofile=showprofile,
139 showprofile=showprofile,
140 show=show)
140 show=show)
141
141
142 if xmin == None: xmin = numpy.nanmin(x)
142 if xmin == None: xmin = numpy.nanmin(x)
143 if xmax == None: xmax = numpy.nanmax(x)
143 if xmax == None: xmax = numpy.nanmax(x)
144 if ymin == None: ymin = numpy.nanmin(y)
144 if ymin == None: ymin = numpy.nanmin(y)
145 if ymax == None: ymax = numpy.nanmax(y)
145 if ymax == None: ymax = numpy.nanmax(y)
146 if zmin == None: zmin = numpy.nanmin(avgdB)*0.9
146 if zmin == None: zmin = numpy.nanmin(avgdB)*0.9
147 if zmax == None: zmax = numpy.nanmax(avgdB)*0.9
147 if zmax == None: zmax = numpy.nanmax(avgdB)*0.9
148
148
149 self.FTP_WEI = ftp_wei
149 self.FTP_WEI = ftp_wei
150 self.EXP_CODE = exp_code
150 self.EXP_CODE = exp_code
151 self.SUB_EXP_CODE = sub_exp_code
151 self.SUB_EXP_CODE = sub_exp_code
152 self.PLOT_POS = plot_pos
152 self.PLOT_POS = plot_pos
153
153
154 self.isConfig = True
154 self.isConfig = True
155
155
156 self.setWinTitle(title)
156 self.setWinTitle(title)
157
157
158 for i in range(self.nplots):
158 for i in range(self.nplots):
159 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
159 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
160 title = "Channel %d: %4.2fdB: %s" %(dataOut.channelList[i]+1, noisedB[i], str_datetime)
160 title = "Channel %d: %4.2fdB: %s" %(dataOut.channelList[i]+1, noisedB[i], str_datetime)
161 axes = self.axesList[i*self.__nsubplots]
161 axes = self.axesList[i*self.__nsubplots]
162 axes.pcolor(x, y, zdB[i,:,:],
162 axes.pcolor(x, y, zdB[i,:,:],
163 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
163 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
164 xlabel=xlabel, ylabel=ylabel, title=title,
164 xlabel=xlabel, ylabel=ylabel, title=title,
165 ticksize=9, cblabel='')
165 ticksize=9, cblabel='')
166 #Mean Line
166 #Mean Line
167 mean = dataOut.data_param[i, 1, :]
167 mean = dataOut.data_param[i, 1, :]
168 axes.addpline(mean, y, idline=0, color="black", linestyle="solid", lw=1)
168 axes.addpline(mean, y, idline=0, color="black", linestyle="solid", lw=1)
169
169
170 if self.__showprofile:
170 if self.__showprofile:
171 axes = self.axesList[i*self.__nsubplots +1]
171 axes = self.axesList[i*self.__nsubplots +1]
172 axes.pline(avgdB[i], y,
172 axes.pline(avgdB[i], y,
173 xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
173 xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
174 xlabel='dB', ylabel='', title='',
174 xlabel='dB', ylabel='', title='',
175 ytick_visible=False,
175 ytick_visible=False,
176 grid='x')
176 grid='x')
177
177
178 noiseline = numpy.repeat(noisedB[i], len(y))
178 noiseline = numpy.repeat(noisedB[i], len(y))
179 axes.addpline(noiseline, y, idline=1, color="black", linestyle="dashed", lw=2)
179 axes.addpline(noiseline, y, idline=1, color="black", linestyle="dashed", lw=2)
180
180
181 self.draw()
181 self.draw()
182
182
183 if figfile == None:
183 if figfile == None:
184 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
184 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
185 figfile = self.getFilename(name = str_datetime)
185 figfile = self.getFilename(name = str_datetime)
186
186
187 if figpath != '':
187 if figpath != '':
188 self.counter_imagwr += 1
188 self.counter_imagwr += 1
189 if (self.counter_imagwr>=wr_period):
189 if (self.counter_imagwr>=wr_period):
190 # store png plot to local folder
190 # store png plot to local folder
191 self.saveFigure(figpath, figfile)
191 self.saveFigure(figpath, figfile)
192 # store png plot to FTP server according to RT-Web format
192 # store png plot to FTP server according to RT-Web format
193 name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
193 name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
194 ftp_filename = os.path.join(figpath, name)
194 ftp_filename = os.path.join(figpath, name)
195 self.saveFigure(figpath, ftp_filename)
195 self.saveFigure(figpath, ftp_filename)
196 self.counter_imagwr = 0
196 self.counter_imagwr = 0
197
197
198 class SkyMapPlot(Figure):
198 class SkyMapPlot(Figure):
199
199
200 __isConfig = None
200 __isConfig = None
201 __nsubplots = None
201 __nsubplots = None
202
202
203 WIDTHPROF = None
203 WIDTHPROF = None
204 HEIGHTPROF = None
204 HEIGHTPROF = None
205 PREFIX = 'prm'
205 PREFIX = 'prm'
206
206
207 def __init__(self):
207 def __init__(self):
208
208
209 self.__isConfig = False
209 self.__isConfig = False
210 self.__nsubplots = 1
210 self.__nsubplots = 1
211
211
212 # self.WIDTH = 280
212 # self.WIDTH = 280
213 # self.HEIGHT = 250
213 # self.HEIGHT = 250
214 self.WIDTH = 600
214 self.WIDTH = 600
215 self.HEIGHT = 600
215 self.HEIGHT = 600
216 self.WIDTHPROF = 120
216 self.WIDTHPROF = 120
217 self.HEIGHTPROF = 0
217 self.HEIGHTPROF = 0
218 self.counter_imagwr = 0
218 self.counter_imagwr = 0
219
219
220 self.PLOT_CODE = 1
220 self.PLOT_CODE = 1
221 self.FTP_WEI = None
221 self.FTP_WEI = None
222 self.EXP_CODE = None
222 self.EXP_CODE = None
223 self.SUB_EXP_CODE = None
223 self.SUB_EXP_CODE = None
224 self.PLOT_POS = None
224 self.PLOT_POS = None
225
225
226 def getSubplots(self):
226 def getSubplots(self):
227
227
228 ncol = int(numpy.sqrt(self.nplots)+0.9)
228 ncol = int(numpy.sqrt(self.nplots)+0.9)
229 nrow = int(self.nplots*1./ncol + 0.9)
229 nrow = int(self.nplots*1./ncol + 0.9)
230
230
231 return nrow, ncol
231 return nrow, ncol
232
232
233 def setup(self, id, nplots, wintitle, showprofile=False, show=True):
233 def setup(self, id, nplots, wintitle, showprofile=False, show=True):
234
234
235 self.__showprofile = showprofile
235 self.__showprofile = showprofile
236 self.nplots = nplots
236 self.nplots = nplots
237
237
238 ncolspan = 1
238 ncolspan = 1
239 colspan = 1
239 colspan = 1
240
240
241 self.createFigure(id = id,
241 self.createFigure(id = id,
242 wintitle = wintitle,
242 wintitle = wintitle,
243 widthplot = self.WIDTH, #+ self.WIDTHPROF,
243 widthplot = self.WIDTH, #+ self.WIDTHPROF,
244 heightplot = self.HEIGHT,# + self.HEIGHTPROF,
244 heightplot = self.HEIGHT,# + self.HEIGHTPROF,
245 show=show)
245 show=show)
246
246
247 nrow, ncol = 1,1
247 nrow, ncol = 1,1
248 counter = 0
248 counter = 0
249 x = 0
249 x = 0
250 y = 0
250 y = 0
251 self.addAxes(1, 1, 0, 0, 1, 1, True)
251 self.addAxes(1, 1, 0, 0, 1, 1, True)
252
252
253 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=False,
253 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=False,
254 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
254 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
255 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
255 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
256 server=None, folder=None, username=None, password=None,
256 server=None, folder=None, username=None, password=None,
257 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False):
257 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False):
258
258
259 """
259 """
260
260
261 Input:
261 Input:
262 dataOut :
262 dataOut :
263 id :
263 id :
264 wintitle :
264 wintitle :
265 channelList :
265 channelList :
266 showProfile :
266 showProfile :
267 xmin : None,
267 xmin : None,
268 xmax : None,
268 xmax : None,
269 ymin : None,
269 ymin : None,
270 ymax : None,
270 ymax : None,
271 zmin : None,
271 zmin : None,
272 zmax : None
272 zmax : None
273 """
273 """
274
274
275 arrayParameters = dataOut.data_param
275 arrayParameters = dataOut.data_param
276 error = arrayParameters[:,-1]
276 error = arrayParameters[:,-1]
277 indValid = numpy.where(error == 0)[0]
277 indValid = numpy.where(error == 0)[0]
278 finalMeteor = arrayParameters[indValid,:]
278 finalMeteor = arrayParameters[indValid,:]
279 finalAzimuth = finalMeteor[:,4]
279 finalAzimuth = finalMeteor[:,4]
280 finalZenith = finalMeteor[:,5]
280 finalZenith = finalMeteor[:,5]
281
281
282 x = finalAzimuth*numpy.pi/180
282 x = finalAzimuth*numpy.pi/180
283 y = finalZenith
283 y = finalZenith
284
284
285
285
286 #thisDatetime = dataOut.datatime
286 #thisDatetime = dataOut.datatime
287 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
287 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
288 title = wintitle + " Parameters"
288 title = wintitle + " Parameters"
289 xlabel = "Zonal Zenith Angle (deg) "
289 xlabel = "Zonal Zenith Angle (deg) "
290 ylabel = "Meridional Zenith Angle (deg)"
290 ylabel = "Meridional Zenith Angle (deg)"
291
291
292 if not self.__isConfig:
292 if not self.__isConfig:
293
293
294 nplots = 1
294 nplots = 1
295
295
296 self.setup(id=id,
296 self.setup(id=id,
297 nplots=nplots,
297 nplots=nplots,
298 wintitle=wintitle,
298 wintitle=wintitle,
299 showprofile=showprofile,
299 showprofile=showprofile,
300 show=show)
300 show=show)
301
301
302 self.FTP_WEI = ftp_wei
302 self.FTP_WEI = ftp_wei
303 self.EXP_CODE = exp_code
303 self.EXP_CODE = exp_code
304 self.SUB_EXP_CODE = sub_exp_code
304 self.SUB_EXP_CODE = sub_exp_code
305 self.PLOT_POS = plot_pos
305 self.PLOT_POS = plot_pos
306 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
306 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
307 self.firstdate = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
307 self.firstdate = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
308 self.__isConfig = True
308 self.__isConfig = True
309
309
310 self.setWinTitle(title)
310 self.setWinTitle(title)
311
311
312 i = 0
312 i = 0
313 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
313 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
314
314
315 axes = self.axesList[i*self.__nsubplots]
315 axes = self.axesList[i*self.__nsubplots]
316 nevents = axes.x_buffer.shape[0] + x.shape[0]
316 nevents = axes.x_buffer.shape[0] + x.shape[0]
317 title = "Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n" %(self.firstdate,str_datetime,nevents)
317 title = "Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n" %(self.firstdate,str_datetime,nevents)
318 axes.polar(x, y,
318 axes.polar(x, y,
319 title=title, xlabel=xlabel, ylabel=ylabel,
319 title=title, xlabel=xlabel, ylabel=ylabel,
320 ticksize=9, cblabel='')
320 ticksize=9, cblabel='')
321
321
322 self.draw()
322 self.draw()
323
323
324 if save:
324 if save:
325
325
326 self.counter_imagwr += 1
326 self.counter_imagwr += 1
327 if (self.counter_imagwr==wr_period):
327 if (self.counter_imagwr==wr_period):
328
328
329 if figfile == None:
329 if figfile == None:
330 figfile = self.getFilename(name = self.name)
330 figfile = self.getFilename(name = self.name)
331 self.saveFigure(figpath, figfile)
331 self.saveFigure(figpath, figfile)
332
332
333 if ftp:
333 if ftp:
334 #provisionalmente envia archivos en el formato de la web en tiempo real
334 #provisionalmente envia archivos en el formato de la web en tiempo real
335 name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
335 name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
336 path = '%s%03d' %(self.PREFIX, self.id)
336 path = '%s%03d' %(self.PREFIX, self.id)
337 ftp_file = os.path.join(path,'ftp','%s.png'%name)
337 ftp_file = os.path.join(path,'ftp','%s.png'%name)
338 self.saveFigure(figpath, ftp_file)
338 self.saveFigure(figpath, ftp_file)
339 ftp_filename = os.path.join(figpath,ftp_file)
339 ftp_filename = os.path.join(figpath,ftp_file)
340
340
341
341
342 try:
342 try:
343 self.sendByFTP(ftp_filename, server, folder, username, password)
343 self.sendByFTP(ftp_filename, server, folder, username, password)
344 except:
344 except:
345 self.counter_imagwr = 0
345 self.counter_imagwr = 0
346 raise ValueError, 'Error FTP'
346 raise ValueError, 'Error FTP'
347
347
348 self.counter_imagwr = 0
348 self.counter_imagwr = 0
349
349
350
350
351 class WindProfilerPlot(Figure):
351 class WindProfilerPlot(Figure):
352
352
353 __isConfig = None
353 __isConfig = None
354 __nsubplots = None
354 __nsubplots = None
355
355
356 WIDTHPROF = None
356 WIDTHPROF = None
357 HEIGHTPROF = None
357 HEIGHTPROF = None
358 PREFIX = 'wind'
358 PREFIX = 'wind'
359
359
360 def __init__(self):
360 def __init__(self):
361
361
362 self.timerange = 2*60*60
362 self.timerange = 2*60*60
363 self.__isConfig = False
363 self.__isConfig = False
364 self.__nsubplots = 1
364 self.__nsubplots = 1
365
365
366 self.WIDTH = 800
366 self.WIDTH = 800
367 self.HEIGHT = 150
367 self.HEIGHT = 150
368 self.WIDTHPROF = 120
368 self.WIDTHPROF = 120
369 self.HEIGHTPROF = 0
369 self.HEIGHTPROF = 0
370 self.counter_imagwr = 0
370 self.counter_imagwr = 0
371
371
372 self.PLOT_CODE = 0
372 self.PLOT_CODE = 0
373 self.FTP_WEI = None
373 self.FTP_WEI = None
374 self.EXP_CODE = None
374 self.EXP_CODE = None
375 self.SUB_EXP_CODE = None
375 self.SUB_EXP_CODE = None
376 self.PLOT_POS = None
376 self.PLOT_POS = None
377 self.tmin = None
377 self.tmin = None
378 self.tmax = None
378 self.tmax = None
379
379
380 self.xmin = None
380 self.xmin = None
381 self.xmax = None
381 self.xmax = None
382
382
383 self.figfile = None
383 self.figfile = None
384
384
385 def getSubplots(self):
385 def getSubplots(self):
386
386
387 ncol = 1
387 ncol = 1
388 nrow = self.nplots
388 nrow = self.nplots
389
389
390 return nrow, ncol
390 return nrow, ncol
391
391
392 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
392 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
393
393
394 self.__showprofile = showprofile
394 self.__showprofile = showprofile
395 self.nplots = nplots
395 self.nplots = nplots
396
396
397 ncolspan = 1
397 ncolspan = 1
398 colspan = 1
398 colspan = 1
399
399
400 self.createFigure(id = id,
400 self.createFigure(id = id,
401 wintitle = wintitle,
401 wintitle = wintitle,
402 widthplot = self.WIDTH + self.WIDTHPROF,
402 widthplot = self.WIDTH + self.WIDTHPROF,
403 heightplot = self.HEIGHT + self.HEIGHTPROF,
403 heightplot = self.HEIGHT + self.HEIGHTPROF,
404 show=show)
404 show=show)
405
405
406 nrow, ncol = self.getSubplots()
406 nrow, ncol = self.getSubplots()
407
407
408 counter = 0
408 counter = 0
409 for y in range(nrow):
409 for y in range(nrow):
410 if counter >= self.nplots:
410 if counter >= self.nplots:
411 break
411 break
412
412
413 self.addAxes(nrow, ncol*ncolspan, y, 0, colspan, 1)
413 self.addAxes(nrow, ncol*ncolspan, y, 0, colspan, 1)
414 counter += 1
414 counter += 1
415
415
416 def run(self, dataOut, id, wintitle="", channelList=None,
416 def run(self, dataOut, id, wintitle="", channelList=None,
417 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
417 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
418 zmax_ver = None, zmin_ver = None, SNRmin = None, SNRmax = None,
418 zmax_ver = None, zmin_ver = None, SNRmin = None, SNRmax = None,
419 timerange=None, SNRthresh = None,
419 timerange=None, SNRthresh = None,
420 save=False, figpath='', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
420 save=False, figpath='', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
421 server=None, folder=None, username=None, password=None,
421 server=None, folder=None, username=None, password=None,
422 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
422 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
423 """
423 """
424
424
425 Input:
425 Input:
426 dataOut :
426 dataOut :
427 id :
427 id :
428 wintitle :
428 wintitle :
429 channelList :
429 channelList :
430 showProfile :
430 showProfile :
431 xmin : None,
431 xmin : None,
432 xmax : None,
432 xmax : None,
433 ymin : None,
433 ymin : None,
434 ymax : None,
434 ymax : None,
435 zmin : None,
435 zmin : None,
436 zmax : None
436 zmax : None
437 """
437 """
438
438
439 if channelList == None:
439 if channelList == None:
440 channelIndexList = dataOut.channelIndexList
440 channelIndexList = dataOut.channelIndexList
441 else:
441 else:
442 channelIndexList = []
442 channelIndexList = []
443 for channel in channelList:
443 for channel in channelList:
444 if channel not in dataOut.channelList:
444 if channel not in dataOut.channelList:
445 raise ValueError, "Channel %d is not in dataOut.channelList"
445 raise ValueError, "Channel %d is not in dataOut.channelList"
446 channelIndexList.append(dataOut.channelList.index(channel))
446 channelIndexList.append(dataOut.channelList.index(channel))
447
447
448 if timerange != None:
448 if timerange != None:
449 self.timerange = timerange
449 self.timerange = timerange
450
450
451 tmin = None
451 tmin = None
452 tmax = None
452 tmax = None
453
453
454 x = dataOut.getTimeRange1()
454 x = dataOut.getTimeRange1()
455 # y = dataOut.heightRange
455 # y = dataOut.heightRange
456 y = dataOut.heightRange
456 y = dataOut.heightRange
457
457
458 z = dataOut.winds.copy()
458 z = dataOut.data_output.copy()
459 nplots = z.shape[0] #Number of wind dimensions estimated
459 nplots = z.shape[0] #Number of wind dimensions estimated
460 nplotsw = nplots
460 nplotsw = nplots
461
461
462 #If there is a SNR function defined
462 #If there is a SNR function defined
463 if dataOut.SNR != None:
463 if dataOut.SNR != None:
464 nplots += 1
464 nplots += 1
465 SNR = dataOut.SNR
465 SNR = dataOut.SNR
466 SNRavg = numpy.average(SNR, axis=0)
466 SNRavg = numpy.average(SNR, axis=0)
467
467
468 SNRdB = 10*numpy.log10(SNR)
468 SNRdB = 10*numpy.log10(SNR)
469 SNRavgdB = 10*numpy.log10(SNRavg)
469 SNRavgdB = 10*numpy.log10(SNRavg)
470
470
471 if SNRthresh == None: SNRthresh = -5.0
471 if SNRthresh == None: SNRthresh = -5.0
472 ind = numpy.where(SNRavg < 10**(SNRthresh/10))[0]
472 ind = numpy.where(SNRavg < 10**(SNRthresh/10))[0]
473
473
474 for i in range(nplotsw):
474 for i in range(nplotsw):
475 z[i,ind] = numpy.nan
475 z[i,ind] = numpy.nan
476
476
477
477
478 showprofile = False
478 showprofile = False
479 # thisDatetime = dataOut.datatime
479 # thisDatetime = dataOut.datatime
480 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
480 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
481 title = wintitle + "Wind"
481 title = wintitle + "Wind"
482 xlabel = ""
482 xlabel = ""
483 ylabel = "Range (Km)"
483 ylabel = "Range (Km)"
484
484
485 if not self.__isConfig:
485 if not self.__isConfig:
486
486
487 self.setup(id=id,
487 self.setup(id=id,
488 nplots=nplots,
488 nplots=nplots,
489 wintitle=wintitle,
489 wintitle=wintitle,
490 showprofile=showprofile,
490 showprofile=showprofile,
491 show=show)
491 show=show)
492
492
493 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
493 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
494
494
495 if ymin == None: ymin = numpy.nanmin(y)
495 if ymin == None: ymin = numpy.nanmin(y)
496 if ymax == None: ymax = numpy.nanmax(y)
496 if ymax == None: ymax = numpy.nanmax(y)
497
497
498 if zmax == None: zmax = numpy.nanmax(abs(z[range(2),:]))
498 if zmax == None: zmax = numpy.nanmax(abs(z[range(2),:]))
499 #if numpy.isnan(zmax): zmax = 50
499 #if numpy.isnan(zmax): zmax = 50
500 if zmin == None: zmin = -zmax
500 if zmin == None: zmin = -zmax
501
501
502 if nplotsw == 3:
502 if nplotsw == 3:
503 if zmax_ver == None: zmax_ver = numpy.nanmax(abs(z[2,:]))
503 if zmax_ver == None: zmax_ver = numpy.nanmax(abs(z[2,:]))
504 if zmin_ver == None: zmin_ver = -zmax_ver
504 if zmin_ver == None: zmin_ver = -zmax_ver
505
505
506 if dataOut.SNR != None:
506 if dataOut.SNR != None:
507 if SNRmin == None: SNRmin = numpy.nanmin(SNRavgdB)
507 if SNRmin == None: SNRmin = numpy.nanmin(SNRavgdB)
508 if SNRmax == None: SNRmax = numpy.nanmax(SNRavgdB)
508 if SNRmax == None: SNRmax = numpy.nanmax(SNRavgdB)
509
509
510 self.FTP_WEI = ftp_wei
510 self.FTP_WEI = ftp_wei
511 self.EXP_CODE = exp_code
511 self.EXP_CODE = exp_code
512 self.SUB_EXP_CODE = sub_exp_code
512 self.SUB_EXP_CODE = sub_exp_code
513 self.PLOT_POS = plot_pos
513 self.PLOT_POS = plot_pos
514
514
515 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
515 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
516 self.__isConfig = True
516 self.__isConfig = True
517
517
518
518
519 self.setWinTitle(title)
519 self.setWinTitle(title)
520
520
521 if ((self.xmax - x[1]) < (x[1]-x[0])):
521 if ((self.xmax - x[1]) < (x[1]-x[0])):
522 x[1] = self.xmax
522 x[1] = self.xmax
523
523
524 strWind = ['Zonal', 'Meridional', 'Vertical']
524 strWind = ['Zonal', 'Meridional', 'Vertical']
525 strCb = ['Velocity (m/s)','Velocity (m/s)','Velocity (cm/s)']
525 strCb = ['Velocity (m/s)','Velocity (m/s)','Velocity (cm/s)']
526 zmaxVector = [zmax, zmax, zmax_ver]
526 zmaxVector = [zmax, zmax, zmax_ver]
527 zminVector = [zmin, zmin, zmin_ver]
527 zminVector = [zmin, zmin, zmin_ver]
528 windFactor = [1,1,100]
528 windFactor = [1,1,100]
529
529
530 for i in range(nplotsw):
530 for i in range(nplotsw):
531
531
532 title = "%s Wind: %s" %(strWind[i], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
532 title = "%s Wind: %s" %(strWind[i], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
533 axes = self.axesList[i*self.__nsubplots]
533 axes = self.axesList[i*self.__nsubplots]
534
534
535 z1 = z[i,:].reshape((1,-1))*windFactor[i]
535 z1 = z[i,:].reshape((1,-1))*windFactor[i]
536
536
537 axes.pcolorbuffer(x, y, z1,
537 axes.pcolorbuffer(x, y, z1,
538 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zminVector[i], zmax=zmaxVector[i],
538 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zminVector[i], zmax=zmaxVector[i],
539 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
539 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
540 ticksize=9, cblabel=strCb[i], cbsize="1%", colormap="RdBu_r" )
540 ticksize=9, cblabel=strCb[i], cbsize="1%", colormap="RdBu_r" )
541
541
542 if dataOut.SNR != None:
542 if dataOut.SNR != None:
543 i += 1
543 i += 1
544 title = "Signal Noise Ratio (SNR): %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
544 title = "Signal Noise Ratio (SNR): %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
545 axes = self.axesList[i*self.__nsubplots]
545 axes = self.axesList[i*self.__nsubplots]
546
546
547 SNRavgdB = SNRavgdB.reshape((1,-1))
547 SNRavgdB = SNRavgdB.reshape((1,-1))
548
548
549 axes.pcolorbuffer(x, y, SNRavgdB,
549 axes.pcolorbuffer(x, y, SNRavgdB,
550 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
550 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
551 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
551 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
552 ticksize=9, cblabel='', cbsize="1%", colormap="jet")
552 ticksize=9, cblabel='', cbsize="1%", colormap="jet")
553
553
554 self.draw()
554 self.draw()
555
555
556 if self.figfile == None:
556 if self.figfile == None:
557 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
557 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
558 self.figfile = self.getFilename(name = str_datetime)
558 self.figfile = self.getFilename(name = str_datetime)
559
559
560 if figpath != '':
560 if figpath != '':
561
561
562 self.counter_imagwr += 1
562 self.counter_imagwr += 1
563 if (self.counter_imagwr>=wr_period):
563 if (self.counter_imagwr>=wr_period):
564 # store png plot to local folder
564 # store png plot to local folder
565 self.saveFigure(figpath, self.figfile)
565 self.saveFigure(figpath, self.figfile)
566 # store png plot to FTP server according to RT-Web format
566 # store png plot to FTP server according to RT-Web format
567 name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
567 name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
568 ftp_filename = os.path.join(figpath, name)
568 ftp_filename = os.path.join(figpath, name)
569 self.saveFigure(figpath, ftp_filename)
569 self.saveFigure(figpath, ftp_filename)
570
570
571 self.counter_imagwr = 0
571 self.counter_imagwr = 0
572
572
573 if x[1] >= self.axesList[0].xmax:
573 if x[1] >= self.axesList[0].xmax:
574 self.counter_imagwr = wr_period
574 self.counter_imagwr = wr_period
575 self.__isConfig = False
575 self.__isConfig = False
576 self.figfile = None
576 self.figfile = None
577
577
578
578
579 class ParametersPlot(Figure):
579 class ParametersPlot(Figure):
580
580
581 __isConfig = None
581 __isConfig = None
582 __nsubplots = None
582 __nsubplots = None
583
583
584 WIDTHPROF = None
584 WIDTHPROF = None
585 HEIGHTPROF = None
585 HEIGHTPROF = None
586 PREFIX = 'prm'
586 PREFIX = 'prm'
587
587
588 def __init__(self):
588 def __init__(self):
589
589
590 self.timerange = 2*60*60
590 self.timerange = 2*60*60
591 self.__isConfig = False
591 self.__isConfig = False
592 self.__nsubplots = 1
592 self.__nsubplots = 1
593
593
594 self.WIDTH = 800
594 self.WIDTH = 800
595 self.HEIGHT = 150
595 self.HEIGHT = 150
596 self.WIDTHPROF = 120
596 self.WIDTHPROF = 120
597 self.HEIGHTPROF = 0
597 self.HEIGHTPROF = 0
598 self.counter_imagwr = 0
598 self.counter_imagwr = 0
599
599
600 self.PLOT_CODE = 0
600 self.PLOT_CODE = 0
601 self.FTP_WEI = None
601 self.FTP_WEI = None
602 self.EXP_CODE = None
602 self.EXP_CODE = None
603 self.SUB_EXP_CODE = None
603 self.SUB_EXP_CODE = None
604 self.PLOT_POS = None
604 self.PLOT_POS = None
605 self.tmin = None
605 self.tmin = None
606 self.tmax = None
606 self.tmax = None
607
607
608 self.xmin = None
608 self.xmin = None
609 self.xmax = None
609 self.xmax = None
610
610
611 self.figfile = None
611 self.figfile = None
612
612
613 def getSubplots(self):
613 def getSubplots(self):
614
614
615 ncol = 1
615 ncol = 1
616 nrow = self.nplots
616 nrow = self.nplots
617
617
618 return nrow, ncol
618 return nrow, ncol
619
619
620 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
620 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
621
621
622 self.__showprofile = showprofile
622 self.__showprofile = showprofile
623 self.nplots = nplots
623 self.nplots = nplots
624
624
625 ncolspan = 1
625 ncolspan = 1
626 colspan = 1
626 colspan = 1
627
627
628 self.createFigure(id = id,
628 self.createFigure(id = id,
629 wintitle = wintitle,
629 wintitle = wintitle,
630 widthplot = self.WIDTH + self.WIDTHPROF,
630 widthplot = self.WIDTH + self.WIDTHPROF,
631 heightplot = self.HEIGHT + self.HEIGHTPROF,
631 heightplot = self.HEIGHT + self.HEIGHTPROF,
632 show=show)
632 show=show)
633
633
634 nrow, ncol = self.getSubplots()
634 nrow, ncol = self.getSubplots()
635
635
636 counter = 0
636 counter = 0
637 for y in range(nrow):
637 for y in range(nrow):
638 for x in range(ncol):
638 for x in range(ncol):
639
639
640 if counter >= self.nplots:
640 if counter >= self.nplots:
641 break
641 break
642
642
643 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
643 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
644
644
645 if showprofile:
645 if showprofile:
646 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
646 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
647
647
648 counter += 1
648 counter += 1
649
649
650 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=False,
650 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=False,
651 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,timerange=None,
651 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,timerange=None,
652 SNRmin = None, SNRmax = None, SNRthresh = None, paramIndex = None, onlyPositive = False,
652 SNRmin = None, SNRmax = None, SNRthresh = None, paramIndex = None, onlyPositive = False,
653 zlabel = "", parameterName = "",
653 zlabel = "", parameterName = "",
654 save=False, figpath='', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
654 save=False, figpath='', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
655 server=None, folder=None, username=None, password=None,
655 server=None, folder=None, username=None, password=None,
656 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
656 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
657
657
658 """
658 """
659
659
660 Input:
660 Input:
661 dataOut :
661 dataOut :
662 id :
662 id :
663 wintitle :
663 wintitle :
664 channelList :
664 channelList :
665 showProfile :
665 showProfile :
666 xmin : None,
666 xmin : None,
667 xmax : None,
667 xmax : None,
668 ymin : None,
668 ymin : None,
669 ymax : None,
669 ymax : None,
670 zmin : None,
670 zmin : None,
671 zmax : None
671 zmax : None
672 """
672 """
673
673
674 if channelList == None:
674 if channelList == None:
675 channelIndexList = dataOut.channelIndexList
675 channelIndexList = dataOut.channelIndexList
676 else:
676 else:
677 channelIndexList = []
677 channelIndexList = []
678 for channel in channelList:
678 for channel in channelList:
679 if channel not in dataOut.channelList:
679 if channel not in dataOut.channelList:
680 raise ValueError, "Channel %d is not in dataOut.channelList"
680 raise ValueError, "Channel %d is not in dataOut.channelList"
681 channelIndexList.append(dataOut.channelList.index(channel))
681 channelIndexList.append(dataOut.channelList.index(channel))
682
682
683 if timerange != None:
683 if timerange != None:
684 self.timerange = timerange
684 self.timerange = timerange
685
685
686 #tmin = None
686 #tmin = None
687 #tmax = None
687 #tmax = None
688 if paramIndex == None:
688 if paramIndex == None:
689 paramIndex = 1
689 paramIndex = 1
690 x = dataOut.getTimeRange1()
690 x = dataOut.getTimeRange1()
691 y = dataOut.heightRange
691 y = dataOut.heightRange
692 z = dataOut.data_param[channelIndexList,paramIndex,:].copy()
692 z = dataOut.data_param[channelIndexList,paramIndex,:].copy()
693
693
694 zRange = dataOut.abscissaRange
694 zRange = dataOut.abscissaRange
695 nplots = z.shape[0] #Number of wind dimensions estimated
695 nplots = z.shape[0] #Number of wind dimensions estimated
696 # thisDatetime = dataOut.datatime
696 # thisDatetime = dataOut.datatime
697 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
697 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
698 title = wintitle + " Parameters Plot" #: %s" %(thisDatetime.strftime("%d-%b-%Y"))
698 title = wintitle + " Parameters Plot" #: %s" %(thisDatetime.strftime("%d-%b-%Y"))
699 xlabel = ""
699 xlabel = ""
700 ylabel = "Range (Km)"
700 ylabel = "Range (Km)"
701
701
702 if onlyPositive:
702 if onlyPositive:
703 colormap = "jet"
703 colormap = "jet"
704 zmin = 0
704 zmin = 0
705 else: colormap = "RdBu_r"
705 else: colormap = "RdBu_r"
706
706
707 if not self.__isConfig:
707 if not self.__isConfig:
708
708
709 self.setup(id=id,
709 self.setup(id=id,
710 nplots=nplots,
710 nplots=nplots,
711 wintitle=wintitle,
711 wintitle=wintitle,
712 showprofile=showprofile,
712 showprofile=showprofile,
713 show=show)
713 show=show)
714
714
715 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
715 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
716
716
717 if ymin == None: ymin = numpy.nanmin(y)
717 if ymin == None: ymin = numpy.nanmin(y)
718 if ymax == None: ymax = numpy.nanmax(y)
718 if ymax == None: ymax = numpy.nanmax(y)
719 if zmin == None: zmin = numpy.nanmin(zRange)
719 if zmin == None: zmin = numpy.nanmin(zRange)
720 if zmax == None: zmax = numpy.nanmax(zRange)
720 if zmax == None: zmax = numpy.nanmax(zRange)
721
721
722 if dataOut.SNR != None:
722 if dataOut.SNR != None:
723 if SNRmin == None: SNRmin = numpy.nanmin(SNRavgdB)
723 if SNRmin == None: SNRmin = numpy.nanmin(SNRavgdB)
724 if SNRmax == None: SNRmax = numpy.nanmax(SNRavgdB)
724 if SNRmax == None: SNRmax = numpy.nanmax(SNRavgdB)
725
725
726 self.FTP_WEI = ftp_wei
726 self.FTP_WEI = ftp_wei
727 self.EXP_CODE = exp_code
727 self.EXP_CODE = exp_code
728 self.SUB_EXP_CODE = sub_exp_code
728 self.SUB_EXP_CODE = sub_exp_code
729 self.PLOT_POS = plot_pos
729 self.PLOT_POS = plot_pos
730
730
731 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
731 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
732 self.__isConfig = True
732 self.__isConfig = True
733 self.figfile = figfile
733 self.figfile = figfile
734
734
735 self.setWinTitle(title)
735 self.setWinTitle(title)
736
736
737 if ((self.xmax - x[1]) < (x[1]-x[0])):
737 if ((self.xmax - x[1]) < (x[1]-x[0])):
738 x[1] = self.xmax
738 x[1] = self.xmax
739
739
740 for i in range(nplots):
740 for i in range(nplots):
741 title = "%s Channel %d: %s" %(parameterName, dataOut.channelList[i]+1, thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
741 title = "%s Channel %d: %s" %(parameterName, dataOut.channelList[i]+1, thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
742
742
743 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
743 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
744 title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith)
744 title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith)
745 axes = self.axesList[i*self.__nsubplots]
745 axes = self.axesList[i*self.__nsubplots]
746 z1 = z[i,:].reshape((1,-1))
746 z1 = z[i,:].reshape((1,-1))
747 axes.pcolorbuffer(x, y, z1,
747 axes.pcolorbuffer(x, y, z1,
748 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
748 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
749 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,colormap=colormap,
749 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,colormap=colormap,
750 ticksize=9, cblabel=zlabel, cbsize="1%")
750 ticksize=9, cblabel=zlabel, cbsize="1%")
751
751
752 self.draw()
752 self.draw()
753
753
754 if self.figfile == None:
754 if self.figfile == None:
755 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
755 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
756 self.figfile = self.getFilename(name = str_datetime)
756 self.figfile = self.getFilename(name = str_datetime)
757
757
758 if figpath != '':
758 if figpath != '':
759
759
760 self.counter_imagwr += 1
760 self.counter_imagwr += 1
761 if (self.counter_imagwr>=wr_period):
761 if (self.counter_imagwr>=wr_period):
762 # store png plot to local folder
762 # store png plot to local folder
763 self.saveFigure(figpath, self.figfile)
763 self.saveFigure(figpath, self.figfile)
764 # store png plot to FTP server according to RT-Web format
764 # store png plot to FTP server according to RT-Web format
765 name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
765 name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
766 ftp_filename = os.path.join(figpath, name)
766 ftp_filename = os.path.join(figpath, name)
767 self.saveFigure(figpath, ftp_filename)
767 self.saveFigure(figpath, ftp_filename)
768
768
769 self.counter_imagwr = 0
769 self.counter_imagwr = 0
770
770
771 if x[1] >= self.axesList[0].xmax:
771 if x[1] >= self.axesList[0].xmax:
772 self.counter_imagwr = wr_period
772 self.counter_imagwr = wr_period
773 self.__isConfig = False
773 self.__isConfig = False
774 self.figfile = None
775
776
777 class SpectralFittingPlot(Figure):
778
779 __isConfig = None
780 __nsubplots = None
781
782 WIDTHPROF = None
783 HEIGHTPROF = None
784 PREFIX = 'prm'
785
786
787 N = None
788 ippSeconds = None
789
790 def __init__(self):
791 self.__isConfig = False
792 self.__nsubplots = 1
793
794 self.WIDTH = 450
795 self.HEIGHT = 250
796 self.WIDTHPROF = 0
797 self.HEIGHTPROF = 0
798
799 def getSubplots(self):
800
801 ncol = int(numpy.sqrt(self.nplots)+0.9)
802 nrow = int(self.nplots*1./ncol + 0.9)
803
804 return nrow, ncol
805
806 def setup(self, id, nplots, wintitle, showprofile=False, show=True):
807
808 showprofile = False
809 self.__showprofile = showprofile
810 self.nplots = nplots
811
812 ncolspan = 5
813 colspan = 4
814 if showprofile:
815 ncolspan = 5
816 colspan = 4
817 self.__nsubplots = 2
818
819 self.createFigure(id = id,
820 wintitle = wintitle,
821 widthplot = self.WIDTH + self.WIDTHPROF,
822 heightplot = self.HEIGHT + self.HEIGHTPROF,
823 show=show)
824
825 nrow, ncol = self.getSubplots()
826
827 counter = 0
828 for y in range(nrow):
829 for x in range(ncol):
830
831 if counter >= self.nplots:
832 break
833
834 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
835
836 if showprofile:
837 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
838
839 counter += 1
840
841 def run(self, dataOut, id, cutHeight=None, fit=False, wintitle="", channelList=None, showprofile=True,
842 xmin=None, xmax=None, ymin=None, ymax=None,
843 save=False, figpath='./', figfile=None, show=True):
844
845 """
846
847 Input:
848 dataOut :
849 id :
850 wintitle :
851 channelList :
852 showProfile :
853 xmin : None,
854 xmax : None,
855 zmin : None,
856 zmax : None
857 """
858
859 if cutHeight==None:
860 h=270
861 heightindex = numpy.abs(cutHeight - dataOut.heightList).argmin()
862 cutHeight = dataOut.heightList[heightindex]
863
864 factor = dataOut.normFactor
865 x = dataOut.abscissaRange[:-1]
866 #y = dataOut.getHeiRange()
867
868 z = dataOut.data_pre[:,:,heightindex]/factor
869 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
870 avg = numpy.average(z, axis=1)
871 listChannels = z.shape[0]
872
873 #Reconstruct Function
874 if fit==True:
875 groupArray = dataOut.groupList
876 listChannels = groupArray.reshape((groupArray.size))
877 listChannels.sort()
878 spcFitLine = numpy.zeros(z.shape)
879 constants = dataOut.constants
880
881 nGroups = groupArray.shape[0]
882 nChannels = groupArray.shape[1]
883 nProfiles = z.shape[1]
884
885 for f in range(nGroups):
886 groupChann = groupArray[f,:]
887 p = dataOut.data_param[f,:,heightindex]
888 # p = numpy.array([ 89.343967,0.14036615,0.17086219,18.89835291,1.58388365,1.55099167])
889 fitLineAux = dataOut.library.modelFunction(p, constants)*nProfiles
890 fitLineAux = fitLineAux.reshape((nChannels,nProfiles))
891 spcFitLine[groupChann,:] = fitLineAux
892 # spcFitLine = spcFitLine/factor
893
894 z = z[listChannels,:]
895 spcFitLine = spcFitLine[listChannels,:]
896 spcFitLinedB = 10*numpy.log10(spcFitLine)
897
898 zdB = 10*numpy.log10(z)
899 #thisDatetime = dataOut.datatime
900 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
901 title = wintitle + " Doppler Spectra: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
902 xlabel = "Velocity (m/s)"
903 ylabel = "Spectrum"
904
905 if not self.__isConfig:
906
907 nplots = listChannels.size
908
909 self.setup(id=id,
910 nplots=nplots,
911 wintitle=wintitle,
912 showprofile=showprofile,
913 show=show)
914
915 if xmin == None: xmin = numpy.nanmin(x)
916 if xmax == None: xmax = numpy.nanmax(x)
917 if ymin == None: ymin = numpy.nanmin(zdB)
918 if ymax == None: ymax = numpy.nanmax(zdB)+2
919
920 self.__isConfig = True
921
922 self.setWinTitle(title)
923 for i in range(self.nplots):
924 # title = "Channel %d: %4.2fdB" %(dataOut.channelList[i]+1, noisedB[i])
925 title = "Height %4.1f km\nChannel %d:" %(cutHeight, listChannels[i]+1)
926 axes = self.axesList[i*self.__nsubplots]
927 if fit == False:
928 axes.pline(x, zdB[i,:],
929 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
930 xlabel=xlabel, ylabel=ylabel, title=title
931 )
932 if fit == True:
933 fitline=spcFitLinedB[i,:]
934 y=numpy.vstack([zdB[i,:],fitline] )
935 legendlabels=['Data','Fitting']
936 axes.pmultilineyaxis(x, y,
937 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
938 xlabel=xlabel, ylabel=ylabel, title=title,
939 legendlabels=legendlabels, marker=None,
940 linestyle='solid', grid='both')
941
942 self.draw()
943
944 if save:
945 date = thisDatetime.strftime("%Y%m%d_%H%M%S")
946 if figfile == None:
947 figfile = self.getFilename(name = date)
948
949 self.saveFigure(figpath, figfile)
950
951
952 class EWDriftsPlot(Figure):
953
954 __isConfig = None
955 __nsubplots = None
956
957 WIDTHPROF = None
958 HEIGHTPROF = None
959 PREFIX = 'drift'
960
961 def __init__(self):
962
963 self.timerange = 2*60*60
964 self.isConfig = False
965 self.__nsubplots = 1
966
967 self.WIDTH = 800
968 self.HEIGHT = 150
969 self.WIDTHPROF = 120
970 self.HEIGHTPROF = 0
971 self.counter_imagwr = 0
972
973 self.PLOT_CODE = 0
974 self.FTP_WEI = None
975 self.EXP_CODE = None
976 self.SUB_EXP_CODE = None
977 self.PLOT_POS = None
978 self.tmin = None
979 self.tmax = None
980
981 self.xmin = None
982 self.xmax = None
983
984 self.figfile = None
985
986 def getSubplots(self):
987
988 ncol = 1
989 nrow = self.nplots
990
991 return nrow, ncol
992
993 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
994
995 self.__showprofile = showprofile
996 self.nplots = nplots
997
998 ncolspan = 1
999 colspan = 1
1000
1001 self.createFigure(id = id,
1002 wintitle = wintitle,
1003 widthplot = self.WIDTH + self.WIDTHPROF,
1004 heightplot = self.HEIGHT + self.HEIGHTPROF,
1005 show=show)
1006
1007 nrow, ncol = self.getSubplots()
1008
1009 counter = 0
1010 for y in range(nrow):
1011 if counter >= self.nplots:
1012 break
1013
1014 self.addAxes(nrow, ncol*ncolspan, y, 0, colspan, 1)
1015 counter += 1
1016
1017 def run(self, dataOut, id, wintitle="", channelList=None,
1018 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
1019 zmaxVertical = None, zminVertical = None, zmaxZonal = None, zminZonal = None,
1020 timerange=None, SNRthresh = -numpy.inf, SNRmin = None, SNRmax = None, SNR_1 = False,
1021 save=False, figpath='', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
1022 server=None, folder=None, username=None, password=None,
1023 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1024 """
1025
1026 Input:
1027 dataOut :
1028 id :
1029 wintitle :
1030 channelList :
1031 showProfile :
1032 xmin : None,
1033 xmax : None,
1034 ymin : None,
1035 ymax : None,
1036 zmin : None,
1037 zmax : None
1038 """
1039
1040 if channelList == None:
1041 channelIndexList = dataOut.channelIndexList
1042 else:
1043 channelIndexList = []
1044 for channel in channelList:
1045 if channel not in dataOut.channelList:
1046 raise ValueError, "Channel %d is not in dataOut.channelList"
1047 channelIndexList.append(dataOut.channelList.index(channel))
1048
1049 if timerange != None:
1050 self.timerange = timerange
1051
1052 tmin = None
1053 tmax = None
1054
1055 x = dataOut.getTimeRange1()
1056 # y = dataOut.heightRange
1057 y = dataOut.heightList
1058
1059 z = dataOut.data_output
1060 nplots = z.shape[0] #Number of wind dimensions estimated
1061 nplotsw = nplots
1062
1063 #If there is a SNR function defined
1064 if dataOut.SNR != None:
1065 nplots += 1
1066 SNR = dataOut.SNR
1067
1068 if SNR_1:
1069 SNR += 1
1070
1071 SNRavg = numpy.average(SNR, axis=0)
1072
1073 SNRdB = 10*numpy.log10(SNR)
1074 SNRavgdB = 10*numpy.log10(SNRavg)
1075
1076 ind = numpy.where(SNRavg < 10**(SNRthresh/10))[0]
1077
1078 for i in range(nplotsw):
1079 z[i,ind] = numpy.nan
1080
1081
1082 showprofile = False
1083 # thisDatetime = dataOut.datatime
1084 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
1085 title = wintitle + " EW Drifts"
1086 xlabel = ""
1087 ylabel = "Height (Km)"
1088
1089 if not self.__isConfig:
1090
1091 self.setup(id=id,
1092 nplots=nplots,
1093 wintitle=wintitle,
1094 showprofile=showprofile,
1095 show=show)
1096
1097 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
1098
1099 if ymin == None: ymin = numpy.nanmin(y)
1100 if ymax == None: ymax = numpy.nanmax(y)
1101
1102 if zmaxZonal == None: zmaxZonal = numpy.nanmax(abs(z[0,:]))
1103 if zminZonal == None: zminZonal = -zmaxZonal
1104 if zmaxVertical == None: zmaxVertical = numpy.nanmax(abs(z[1,:]))
1105 if zminVertical == None: zminVertical = -zmaxVertical
1106
1107 if dataOut.SNR != None:
1108 if SNRmin == None: SNRmin = numpy.nanmin(SNRavgdB)
1109 if SNRmax == None: SNRmax = numpy.nanmax(SNRavgdB)
1110
1111 self.FTP_WEI = ftp_wei
1112 self.EXP_CODE = exp_code
1113 self.SUB_EXP_CODE = sub_exp_code
1114 self.PLOT_POS = plot_pos
1115
1116 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1117 self.__isConfig = True
1118
1119
1120 self.setWinTitle(title)
1121
1122 if ((self.xmax - x[1]) < (x[1]-x[0])):
1123 x[1] = self.xmax
1124
1125 strWind = ['Zonal','Vertical']
1126 strCb = 'Velocity (m/s)'
1127 zmaxVector = [zmaxZonal, zmaxVertical]
1128 zminVector = [zminZonal, zminVertical]
1129
1130 for i in range(nplotsw):
1131
1132 title = "%s Drifts: %s" %(strWind[i], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1133 axes = self.axesList[i*self.__nsubplots]
1134
1135 z1 = z[i,:].reshape((1,-1))
1136
1137 axes.pcolorbuffer(x, y, z1,
1138 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zminVector[i], zmax=zmaxVector[i],
1139 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1140 ticksize=9, cblabel=strCb, cbsize="1%", colormap="RdBu_r")
1141
1142 if dataOut.SNR != None:
1143 i += 1
1144 if SNR_1:
1145 title = "Signal Noise Ratio + 1 (SNR+1): %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1146 else:
1147 title = "Signal Noise Ratio (SNR): %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1148 axes = self.axesList[i*self.__nsubplots]
1149 SNRavgdB = SNRavgdB.reshape((1,-1))
1150
1151 axes.pcolorbuffer(x, y, SNRavgdB,
1152 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
1153 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1154 ticksize=9, cblabel='', cbsize="1%", colormap="jet")
1155
1156 self.draw()
1157
1158 if self.figfile == None:
1159 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
1160 self.figfile = self.getFilename(name = str_datetime)
1161
1162 if figpath != '':
1163
1164 self.counter_imagwr += 1
1165 if (self.counter_imagwr>=wr_period):
1166 # store png plot to local folder
1167 self.saveFigure(figpath, self.figfile)
1168 # store png plot to FTP server according to RT-Web format
1169 name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
1170 ftp_filename = os.path.join(figpath, name)
1171 self.saveFigure(figpath, ftp_filename)
1172
1173 self.counter_imagwr = 0
1174
1175 if x[1] >= self.axesList[0].xmax:
1176 self.counter_imagwr = wr_period
1177 self.__isConfig = False
774 self.figfile = None No newline at end of file
1178 self.figfile = None
Show file before | Show file after | Hide comments
@@ -1,427 +1,427
1 import numpy
1 import numpy
2 import datetime
2 import datetime
3 import sys
3 import sys
4 import matplotlib
4 import matplotlib
5
5
6 if 'linux' in sys.platform:
6 if 'linux' in sys.platform:
7 matplotlib.use("TKAgg")
7 matplotlib.use("TKAgg")
8
8
9 if 'darwin' in sys.platform:
9 if 'darwin' in sys.platform:
10 matplotlib.use("TKAgg")
10 matplotlib.use("TKAgg")
11
11
12 import matplotlib.pyplot
12 import matplotlib.pyplot
13
13
14 from mpl_toolkits.axes_grid1 import make_axes_locatable
14 from mpl_toolkits.axes_grid1 import make_axes_locatable
15 from matplotlib.ticker import *
15 from matplotlib.ticker import *
16
16
17 ###########################################
17 ###########################################
18 #Actualizacion de las funciones del driver
18 #Actualizacion de las funciones del driver
19 ###########################################
19 ###########################################
20
20
21 def createFigure(id, wintitle, width, height, facecolor="w", show=True):
21 def createFigure(id, wintitle, width, height, facecolor="w", show=True):
22
22
23 matplotlib.pyplot.ioff()
23 matplotlib.pyplot.ioff()
24 fig = matplotlib.pyplot.figure(num=id, facecolor=facecolor)
24 fig = matplotlib.pyplot.figure(num=id, facecolor=facecolor)
25 fig.canvas.manager.set_window_title(wintitle)
25 fig.canvas.manager.set_window_title(wintitle)
26 fig.canvas.manager.resize(width, height)
26 fig.canvas.manager.resize(width, height)
27 matplotlib.pyplot.ion()
27 matplotlib.pyplot.ion()
28 if show:
28 if show:
29 matplotlib.pyplot.show()
29 matplotlib.pyplot.show()
30
30
31 return fig
31 return fig
32
32
33 def closeFigure(show=True):
33 def closeFigure(show=True):
34
34
35 matplotlib.pyplot.ioff()
35 matplotlib.pyplot.ioff()
36 if show:
36 if show:
37 matplotlib.pyplot.show()
37 matplotlib.pyplot.show()
38
38
39 return
39 return
40
40
41 def saveFigure(fig, filename):
41 def saveFigure(fig, filename):
42
42
43 matplotlib.pyplot.ioff()
43 matplotlib.pyplot.ioff()
44 fig.savefig(filename)
44 fig.savefig(filename)
45 matplotlib.pyplot.ion()
45 matplotlib.pyplot.ion()
46
46
47 def setWinTitle(fig, title):
47 def setWinTitle(fig, title):
48
48
49 fig.canvas.manager.set_window_title(title)
49 fig.canvas.manager.set_window_title(title)
50
50
51 def setTitle(fig, title):
51 def setTitle(fig, title):
52
52
53 fig.suptitle(title)
53 fig.suptitle(title)
54
54
55 def createAxes(fig, nrow, ncol, xpos, ypos, colspan, rowspan, polar=False):
55 def createAxes(fig, nrow, ncol, xpos, ypos, colspan, rowspan, polar=False):
56
56
57 matplotlib.pyplot.ioff()
57 matplotlib.pyplot.ioff()
58 matplotlib.pyplot.figure(fig.number)
58 matplotlib.pyplot.figure(fig.number)
59 axes = matplotlib.pyplot.subplot2grid((nrow, ncol),
59 axes = matplotlib.pyplot.subplot2grid((nrow, ncol),
60 (xpos, ypos),
60 (xpos, ypos),
61 colspan=colspan,
61 colspan=colspan,
62 rowspan=rowspan,
62 rowspan=rowspan,
63 polar=polar)
63 polar=polar)
64
64
65 matplotlib.pyplot.ion()
65 matplotlib.pyplot.ion()
66 return axes
66 return axes
67
67
68 def setAxesText(ax, text):
68 def setAxesText(ax, text):
69
69
70 ax.annotate(text,
70 ax.annotate(text,
71 xy = (.1, .99),
71 xy = (.1, .99),
72 xycoords = 'figure fraction',
72 xycoords = 'figure fraction',
73 horizontalalignment = 'left',
73 horizontalalignment = 'left',
74 verticalalignment = 'top',
74 verticalalignment = 'top',
75 fontsize = 10)
75 fontsize = 10)
76
76
77 def printLabels(ax, xlabel, ylabel, title):
77 def printLabels(ax, xlabel, ylabel, title):
78
78
79 ax.set_xlabel(xlabel, size=11)
79 ax.set_xlabel(xlabel, size=11)
80 ax.set_ylabel(ylabel, size=11)
80 ax.set_ylabel(ylabel, size=11)
81 ax.set_title(title, size=8)
81 ax.set_title(title, size=8)
82
82
83 def createPline(ax, x, y, xmin, xmax, ymin, ymax, xlabel='', ylabel='', title='',
83 def createPline(ax, x, y, xmin, xmax, ymin, ymax, xlabel='', ylabel='', title='',
84 ticksize=9, xtick_visible=True, ytick_visible=True,
84 ticksize=9, xtick_visible=True, ytick_visible=True,
85 nxticks=4, nyticks=10,
85 nxticks=4, nyticks=10,
86 grid=None,color='blue'):
86 grid=None,color='blue'):
87
87
88 """
88 """
89
89
90 Input:
90 Input:
91 grid : None, 'both', 'x', 'y'
91 grid : None, 'both', 'x', 'y'
92 """
92 """
93
93
94 matplotlib.pyplot.ioff()
94 matplotlib.pyplot.ioff()
95
95
96 ax.set_xlim([xmin,xmax])
96 ax.set_xlim([xmin,xmax])
97 ax.set_ylim([ymin,ymax])
97 ax.set_ylim([ymin,ymax])
98
98
99 printLabels(ax, xlabel, ylabel, title)
99 printLabels(ax, xlabel, ylabel, title)
100
100
101 ######################################################
101 ######################################################
102 if (xmax-xmin)<=1:
102 if (xmax-xmin)<=1:
103 xtickspos = numpy.linspace(xmin,xmax,nxticks)
103 xtickspos = numpy.linspace(xmin,xmax,nxticks)
104 xtickspos = numpy.array([float("%.1f"%i) for i in xtickspos])
104 xtickspos = numpy.array([float("%.1f"%i) for i in xtickspos])
105 ax.set_xticks(xtickspos)
105 ax.set_xticks(xtickspos)
106 else:
106 else:
107 xtickspos = numpy.arange(nxticks)*int((xmax-xmin)/(nxticks)) + int(xmin)
107 xtickspos = numpy.arange(nxticks)*int((xmax-xmin)/(nxticks)) + int(xmin)
108 # xtickspos = numpy.arange(nxticks)*float(xmax-xmin)/float(nxticks) + int(xmin)
108 # xtickspos = numpy.arange(nxticks)*float(xmax-xmin)/float(nxticks) + int(xmin)
109 ax.set_xticks(xtickspos)
109 ax.set_xticks(xtickspos)
110
110
111 for tick in ax.get_xticklabels():
111 for tick in ax.get_xticklabels():
112 tick.set_visible(xtick_visible)
112 tick.set_visible(xtick_visible)
113
113
114 for tick in ax.xaxis.get_major_ticks():
114 for tick in ax.xaxis.get_major_ticks():
115 tick.label.set_fontsize(ticksize)
115 tick.label.set_fontsize(ticksize)
116
116
117 ######################################################
117 ######################################################
118 for tick in ax.get_yticklabels():
118 for tick in ax.get_yticklabels():
119 tick.set_visible(ytick_visible)
119 tick.set_visible(ytick_visible)
120
120
121 for tick in ax.yaxis.get_major_ticks():
121 for tick in ax.yaxis.get_major_ticks():
122 tick.label.set_fontsize(ticksize)
122 tick.label.set_fontsize(ticksize)
123
123
124 ax.plot(x, y, color=color)
124 ax.plot(x, y, color=color)
125 iplot = ax.lines[-1]
125 iplot = ax.lines[-1]
126
126
127 ######################################################
127 ######################################################
128 if '0.' in matplotlib.__version__[0:2]:
128 if '0.' in matplotlib.__version__[0:2]:
129 print "The matplotlib version has to be updated to 1.1 or newer"
129 print "The matplotlib version has to be updated to 1.1 or newer"
130 return iplot
130 return iplot
131
131
132 if '1.0.' in matplotlib.__version__[0:4]:
132 if '1.0.' in matplotlib.__version__[0:4]:
133 print "The matplotlib version has to be updated to 1.1 or newer"
133 print "The matplotlib version has to be updated to 1.1 or newer"
134 return iplot
134 return iplot
135
135
136 if grid != None:
136 if grid != None:
137 ax.grid(b=True, which='major', axis=grid)
137 ax.grid(b=True, which='major', axis=grid)
138
138
139 matplotlib.pyplot.tight_layout()
139 matplotlib.pyplot.tight_layout()
140
140
141 matplotlib.pyplot.ion()
141 matplotlib.pyplot.ion()
142
142
143 return iplot
143 return iplot
144
144
145 def set_linedata(ax, x, y, idline):
145 def set_linedata(ax, x, y, idline):
146
146
147 ax.lines[idline].set_data(x,y)
147 ax.lines[idline].set_data(x,y)
148
148
149 def pline(iplot, x, y, xlabel='', ylabel='', title=''):
149 def pline(iplot, x, y, xlabel='', ylabel='', title=''):
150
150
151 ax = iplot.get_axes()
151 ax = iplot.get_axes()
152
152
153 printLabels(ax, xlabel, ylabel, title)
153 printLabels(ax, xlabel, ylabel, title)
154
154
155 set_linedata(ax, x, y, idline=0)
155 set_linedata(ax, x, y, idline=0)
156
156
157 def addpline(ax, x, y, color, linestyle, lw):
157 def addpline(ax, x, y, color, linestyle, lw):
158
158
159 ax.plot(x,y,color=color,linestyle=linestyle,lw=lw)
159 ax.plot(x,y,color=color,linestyle=linestyle,lw=lw)
160
160
161
161
162 def createPcolor(ax, x, y, z, xmin, xmax, ymin, ymax, zmin, zmax,
162 def createPcolor(ax, x, y, z, xmin, xmax, ymin, ymax, zmin, zmax,
163 xlabel='', ylabel='', title='', ticksize = 9,
163 xlabel='', ylabel='', title='', ticksize = 9,
164 colormap='jet',cblabel='', cbsize="5%",
164 colormap='jet',cblabel='', cbsize="5%",
165 XAxisAsTime=False):
165 XAxisAsTime=False):
166
166
167 matplotlib.pyplot.ioff()
167 matplotlib.pyplot.ioff()
168
168
169 divider = make_axes_locatable(ax)
169 divider = make_axes_locatable(ax)
170 ax_cb = divider.new_horizontal(size=cbsize, pad=0.05)
170 ax_cb = divider.new_horizontal(size=cbsize, pad=0.05)
171 fig = ax.get_figure()
171 fig = ax.get_figure()
172 fig.add_axes(ax_cb)
172 fig.add_axes(ax_cb)
173
173
174 ax.set_xlim([xmin,xmax])
174 ax.set_xlim([xmin,xmax])
175 ax.set_ylim([ymin,ymax])
175 ax.set_ylim([ymin,ymax])
176
176
177 printLabels(ax, xlabel, ylabel, title)
177 printLabels(ax, xlabel, ylabel, title)
178
178
179 imesh = ax.pcolormesh(x,y,z.T, vmin=zmin, vmax=zmax, cmap=matplotlib.pyplot.get_cmap(colormap))
179 imesh = ax.pcolormesh(x,y,z.T, vmin=zmin, vmax=zmax, cmap=matplotlib.pyplot.get_cmap(colormap))
180 cb = matplotlib.pyplot.colorbar(imesh, cax=ax_cb)
180 cb = matplotlib.pyplot.colorbar(imesh, cax=ax_cb)
181 cb.set_label(cblabel)
181 cb.set_label(cblabel)
182
182
183 # for tl in ax_cb.get_yticklabels():
183 # for tl in ax_cb.get_yticklabels():
184 # tl.set_visible(True)
184 # tl.set_visible(True)
185
185
186 for tick in ax.yaxis.get_major_ticks():
186 for tick in ax.yaxis.get_major_ticks():
187 tick.label.set_fontsize(ticksize)
187 tick.label.set_fontsize(ticksize)
188
188
189 for tick in ax.xaxis.get_major_ticks():
189 for tick in ax.xaxis.get_major_ticks():
190 tick.label.set_fontsize(ticksize)
190 tick.label.set_fontsize(ticksize)
191
191
192 for tick in cb.ax.get_yticklabels():
192 for tick in cb.ax.get_yticklabels():
193 tick.set_fontsize(ticksize)
193 tick.set_fontsize(ticksize)
194
194
195 ax_cb.yaxis.tick_right()
195 ax_cb.yaxis.tick_right()
196
196
197 if '0.' in matplotlib.__version__[0:2]:
197 if '0.' in matplotlib.__version__[0:2]:
198 print "The matplotlib version has to be updated to 1.1 or newer"
198 print "The matplotlib version has to be updated to 1.1 or newer"
199 return imesh
199 return imesh
200
200
201 if '1.0.' in matplotlib.__version__[0:4]:
201 if '1.0.' in matplotlib.__version__[0:4]:
202 print "The matplotlib version has to be updated to 1.1 or newer"
202 print "The matplotlib version has to be updated to 1.1 or newer"
203 return imesh
203 return imesh
204
204
205 matplotlib.pyplot.tight_layout()
205 matplotlib.pyplot.tight_layout()
206
206
207 if XAxisAsTime:
207 if XAxisAsTime:
208
208
209 func = lambda x, pos: ('%s') %(datetime.datetime.utcfromtimestamp(x).strftime("%H:%M:%S"))
209 func = lambda x, pos: ('%s') %(datetime.datetime.utcfromtimestamp(x).strftime("%H:%M:%S"))
210 ax.xaxis.set_major_formatter(FuncFormatter(func))
210 ax.xaxis.set_major_formatter(FuncFormatter(func))
211 ax.xaxis.set_major_locator(LinearLocator(7))
211 ax.xaxis.set_major_locator(LinearLocator(7))
212
212
213 matplotlib.pyplot.ion()
213 matplotlib.pyplot.ion()
214 return imesh
214 return imesh
215
215
216 def pcolor(imesh, z, xlabel='', ylabel='', title=''):
216 def pcolor(imesh, z, xlabel='', ylabel='', title=''):
217
217
218 z = z.T
218 z = z.T
219
219
220 ax = imesh.get_axes()
220 ax = imesh.get_axes()
221
221
222 printLabels(ax, xlabel, ylabel, title)
222 printLabels(ax, xlabel, ylabel, title)
223
223
224 imesh.set_array(z.ravel())
224 imesh.set_array(z.ravel())
225
225
226 def addpcolor(ax, x, y, z, zmin, zmax, xlabel='', ylabel='', title='', colormap='jet'):
226 def addpcolor(ax, x, y, z, zmin, zmax, xlabel='', ylabel='', title='', colormap='jet'):
227
227
228 printLabels(ax, xlabel, ylabel, title)
228 printLabels(ax, xlabel, ylabel, title)
229
229
230 ax.pcolormesh(x,y,z.T,vmin=zmin,vmax=zmax, cmap=matplotlib.pyplot.get_cmap(colormap))
230 ax.pcolormesh(x,y,z.T,vmin=zmin,vmax=zmax, cmap=matplotlib.pyplot.get_cmap(colormap))
231
231
232 def addpcolorbuffer(ax, x, y, z, zmin, zmax, xlabel='', ylabel='', title='', colormap='jet'):
232 def addpcolorbuffer(ax, x, y, z, zmin, zmax, xlabel='', ylabel='', title='', colormap='jet'):
233
233
234 printLabels(ax, xlabel, ylabel, title)
234 printLabels(ax, xlabel, ylabel, title)
235
235
236 ax.collections.remove(ax.collections[0])
236 ax.collections.remove(ax.collections[0])
237
237
238 ax.pcolormesh(x,y,z.T,vmin=zmin,vmax=zmax, cmap=matplotlib.pyplot.get_cmap(colormap))
238 ax.pcolormesh(x,y,z.T,vmin=zmin,vmax=zmax, cmap=matplotlib.pyplot.get_cmap(colormap))
239
239
240 def createPmultiline(ax, x, y, xmin, xmax, ymin, ymax, xlabel='', ylabel='', title='', legendlabels=None,
240 def createPmultiline(ax, x, y, xmin, xmax, ymin, ymax, xlabel='', ylabel='', title='', legendlabels=None,
241 ticksize=9, xtick_visible=True, ytick_visible=True,
241 ticksize=9, xtick_visible=True, ytick_visible=True,
242 nxticks=4, nyticks=10,
242 nxticks=4, nyticks=10,
243 grid=None):
243 grid=None):
244
244
245 """
245 """
246
246
247 Input:
247 Input:
248 grid : None, 'both', 'x', 'y'
248 grid : None, 'both', 'x', 'y'
249 """
249 """
250
250
251 matplotlib.pyplot.ioff()
251 matplotlib.pyplot.ioff()
252
252
253 lines = ax.plot(x.T, y)
253 lines = ax.plot(x.T, y)
254 leg = ax.legend(lines, legendlabels, loc='upper right')
254 leg = ax.legend(lines, legendlabels, loc='upper right')
255 leg.get_frame().set_alpha(0.5)
255 leg.get_frame().set_alpha(0.5)
256 ax.set_xlim([xmin,xmax])
256 ax.set_xlim([xmin,xmax])
257 ax.set_ylim([ymin,ymax])
257 ax.set_ylim([ymin,ymax])
258 printLabels(ax, xlabel, ylabel, title)
258 printLabels(ax, xlabel, ylabel, title)
259
259
260 xtickspos = numpy.arange(nxticks)*int((xmax-xmin)/(nxticks)) + int(xmin)
260 xtickspos = numpy.arange(nxticks)*int((xmax-xmin)/(nxticks)) + int(xmin)
261 ax.set_xticks(xtickspos)
261 ax.set_xticks(xtickspos)
262
262
263 for tick in ax.get_xticklabels():
263 for tick in ax.get_xticklabels():
264 tick.set_visible(xtick_visible)
264 tick.set_visible(xtick_visible)
265
265
266 for tick in ax.xaxis.get_major_ticks():
266 for tick in ax.xaxis.get_major_ticks():
267 tick.label.set_fontsize(ticksize)
267 tick.label.set_fontsize(ticksize)
268
268
269 for tick in ax.get_yticklabels():
269 for tick in ax.get_yticklabels():
270 tick.set_visible(ytick_visible)
270 tick.set_visible(ytick_visible)
271
271
272 for tick in ax.yaxis.get_major_ticks():
272 for tick in ax.yaxis.get_major_ticks():
273 tick.label.set_fontsize(ticksize)
273 tick.label.set_fontsize(ticksize)
274
274
275 iplot = ax.lines[-1]
275 iplot = ax.lines[-1]
276
276
277 if '0.' in matplotlib.__version__[0:2]:
277 if '0.' in matplotlib.__version__[0:2]:
278 print "The matplotlib version has to be updated to 1.1 or newer"
278 print "The matplotlib version has to be updated to 1.1 or newer"
279 return iplot
279 return iplot
280
280
281 if '1.0.' in matplotlib.__version__[0:4]:
281 if '1.0.' in matplotlib.__version__[0:4]:
282 print "The matplotlib version has to be updated to 1.1 or newer"
282 print "The matplotlib version has to be updated to 1.1 or newer"
283 return iplot
283 return iplot
284
284
285 if grid != None:
285 if grid != None:
286 ax.grid(b=True, which='major', axis=grid)
286 ax.grid(b=True, which='major', axis=grid)
287
287
288 matplotlib.pyplot.tight_layout()
288 matplotlib.pyplot.tight_layout()
289
289
290 matplotlib.pyplot.ion()
290 matplotlib.pyplot.ion()
291
291
292 return iplot
292 return iplot
293
293
294
294
295 def pmultiline(iplot, x, y, xlabel='', ylabel='', title=''):
295 def pmultiline(iplot, x, y, xlabel='', ylabel='', title=''):
296
296
297 ax = iplot.get_axes()
297 ax = iplot.get_axes()
298
298
299 printLabels(ax, xlabel, ylabel, title)
299 printLabels(ax, xlabel, ylabel, title)
300
300
301 for i in range(len(ax.lines)):
301 for i in range(len(ax.lines)):
302 line = ax.lines[i]
302 line = ax.lines[i]
303 line.set_data(x[i,:],y)
303 line.set_data(x[i,:],y)
304
304
305 def createPmultilineYAxis(ax, x, y, xmin, xmax, ymin, ymax, xlabel='', ylabel='', title='', legendlabels=None,
305 def createPmultilineYAxis(ax, x, y, xmin, xmax, ymin, ymax, xlabel='', ylabel='', title='', legendlabels=None,
306 ticksize=9, xtick_visible=True, ytick_visible=True,
306 ticksize=9, xtick_visible=True, ytick_visible=True,
307 nxticks=4, nyticks=10, marker='.', markersize=10, linestyle="None",
307 nxticks=4, nyticks=10, marker='.', markersize=10, linestyle="None",
308 grid=None, XAxisAsTime=False):
308 grid=None, XAxisAsTime=False):
309
309
310 """
310 """
311
311
312 Input:
312 Input:
313 grid : None, 'both', 'x', 'y'
313 grid : None, 'both', 'x', 'y'
314 """
314 """
315
315
316 matplotlib.pyplot.ioff()
316 matplotlib.pyplot.ioff()
317
317
318 # lines = ax.plot(x, y.T, marker=marker,markersize=markersize,linestyle=linestyle)
318 # lines = ax.plot(x, y.T, marker=marker,markersize=markersize,linestyle=linestyle)
319 lines = ax.plot(x, y.T, linestyle='None', marker='.', markersize=markersize)
319 lines = ax.plot(x, y.T, linestyle=linestyle, marker=marker, markersize=markersize)
320 leg = ax.legend(lines, legendlabels, loc='upper left', bbox_to_anchor=(1.01, 1.00), numpoints=1, handlelength=1.5, \
320 leg = ax.legend(lines, legendlabels, loc='upper left', bbox_to_anchor=(1.01, 1.00), numpoints=1, handlelength=1.5, \
321 handletextpad=0.5, borderpad=0.5, labelspacing=0.5, borderaxespad=0.)
321 handletextpad=0.5, borderpad=0.5, labelspacing=0.5, borderaxespad=0.)
322
322
323 for label in leg.get_texts(): label.set_fontsize(9)
323 for label in leg.get_texts(): label.set_fontsize(9)
324
324
325 ax.set_xlim([xmin,xmax])
325 ax.set_xlim([xmin,xmax])
326 ax.set_ylim([ymin,ymax])
326 ax.set_ylim([ymin,ymax])
327 printLabels(ax, xlabel, ylabel, title)
327 printLabels(ax, xlabel, ylabel, title)
328
328
329 # xtickspos = numpy.arange(nxticks)*int((xmax-xmin)/(nxticks)) + int(xmin)
329 # xtickspos = numpy.arange(nxticks)*int((xmax-xmin)/(nxticks)) + int(xmin)
330 # ax.set_xticks(xtickspos)
330 # ax.set_xticks(xtickspos)
331
331
332 for tick in ax.get_xticklabels():
332 for tick in ax.get_xticklabels():
333 tick.set_visible(xtick_visible)
333 tick.set_visible(xtick_visible)
334
334
335 for tick in ax.xaxis.get_major_ticks():
335 for tick in ax.xaxis.get_major_ticks():
336 tick.label.set_fontsize(ticksize)
336 tick.label.set_fontsize(ticksize)
337
337
338 for tick in ax.get_yticklabels():
338 for tick in ax.get_yticklabels():
339 tick.set_visible(ytick_visible)
339 tick.set_visible(ytick_visible)
340
340
341 for tick in ax.yaxis.get_major_ticks():
341 for tick in ax.yaxis.get_major_ticks():
342 tick.label.set_fontsize(ticksize)
342 tick.label.set_fontsize(ticksize)
343
343
344 iplot = ax.lines[-1]
344 iplot = ax.lines[-1]
345
345
346 if '0.' in matplotlib.__version__[0:2]:
346 if '0.' in matplotlib.__version__[0:2]:
347 print "The matplotlib version has to be updated to 1.1 or newer"
347 print "The matplotlib version has to be updated to 1.1 or newer"
348 return iplot
348 return iplot
349
349
350 if '1.0.' in matplotlib.__version__[0:4]:
350 if '1.0.' in matplotlib.__version__[0:4]:
351 print "The matplotlib version has to be updated to 1.1 or newer"
351 print "The matplotlib version has to be updated to 1.1 or newer"
352 return iplot
352 return iplot
353
353
354 if grid != None:
354 if grid != None:
355 ax.grid(b=True, which='major', axis=grid)
355 ax.grid(b=True, which='major', axis=grid)
356
356
357 matplotlib.pyplot.tight_layout()
357 matplotlib.pyplot.tight_layout()
358
358
359 if XAxisAsTime:
359 if XAxisAsTime:
360
360
361 func = lambda x, pos: ('%s') %(datetime.datetime.utcfromtimestamp(x).strftime("%H:%M:%S"))
361 func = lambda x, pos: ('%s') %(datetime.datetime.utcfromtimestamp(x).strftime("%H:%M:%S"))
362 ax.xaxis.set_major_formatter(FuncFormatter(func))
362 ax.xaxis.set_major_formatter(FuncFormatter(func))
363 ax.xaxis.set_major_locator(LinearLocator(7))
363 ax.xaxis.set_major_locator(LinearLocator(7))
364
364
365 matplotlib.pyplot.ion()
365 matplotlib.pyplot.ion()
366
366
367 return iplot
367 return iplot
368
368
369 def pmultilineyaxis(iplot, x, y, xlabel='', ylabel='', title=''):
369 def pmultilineyaxis(iplot, x, y, xlabel='', ylabel='', title=''):
370
370
371 ax = iplot.get_axes()
371 ax = iplot.get_axes()
372
372
373 printLabels(ax, xlabel, ylabel, title)
373 printLabels(ax, xlabel, ylabel, title)
374
374
375 for i in range(len(ax.lines)):
375 for i in range(len(ax.lines)):
376 line = ax.lines[i]
376 line = ax.lines[i]
377 line.set_data(x,y[i,:])
377 line.set_data(x,y[i,:])
378
378
379 def createPolar(ax, x, y,
379 def createPolar(ax, x, y,
380 xlabel='', ylabel='', title='', ticksize = 9,
380 xlabel='', ylabel='', title='', ticksize = 9,
381 colormap='jet',cblabel='', cbsize="5%",
381 colormap='jet',cblabel='', cbsize="5%",
382 XAxisAsTime=False):
382 XAxisAsTime=False):
383
383
384 matplotlib.pyplot.ioff()
384 matplotlib.pyplot.ioff()
385
385
386 ax.plot(x,y,'bo', markersize=5)
386 ax.plot(x,y,'bo', markersize=5)
387 # ax.set_rmax(90)
387 # ax.set_rmax(90)
388 ax.set_ylim(0,90)
388 ax.set_ylim(0,90)
389 ax.set_yticks(numpy.arange(0,90,20))
389 ax.set_yticks(numpy.arange(0,90,20))
390 ax.text(0, -110, ylabel, rotation='vertical', va ='center', ha = 'center' ,size='11')
390 ax.text(0, -110, ylabel, rotation='vertical', va ='center', ha = 'center' ,size='11')
391 # ax.text(100, 100, 'example', ha='left', va='center', rotation='vertical')
391 # ax.text(100, 100, 'example', ha='left', va='center', rotation='vertical')
392 printLabels(ax, xlabel, '', title)
392 printLabels(ax, xlabel, '', title)
393 iplot = ax.lines[-1]
393 iplot = ax.lines[-1]
394
394
395 if '0.' in matplotlib.__version__[0:2]:
395 if '0.' in matplotlib.__version__[0:2]:
396 print "The matplotlib version has to be updated to 1.1 or newer"
396 print "The matplotlib version has to be updated to 1.1 or newer"
397 return iplot
397 return iplot
398
398
399 if '1.0.' in matplotlib.__version__[0:4]:
399 if '1.0.' in matplotlib.__version__[0:4]:
400 print "The matplotlib version has to be updated to 1.1 or newer"
400 print "The matplotlib version has to be updated to 1.1 or newer"
401 return iplot
401 return iplot
402
402
403 # if grid != None:
403 # if grid != None:
404 # ax.grid(b=True, which='major', axis=grid)
404 # ax.grid(b=True, which='major', axis=grid)
405
405
406 matplotlib.pyplot.tight_layout()
406 matplotlib.pyplot.tight_layout()
407
407
408 matplotlib.pyplot.ion()
408 matplotlib.pyplot.ion()
409
409
410
410
411 return iplot
411 return iplot
412
412
413 def polar(iplot, x, y, xlabel='', ylabel='', title=''):
413 def polar(iplot, x, y, xlabel='', ylabel='', title=''):
414
414
415 ax = iplot.get_axes()
415 ax = iplot.get_axes()
416
416
417 # ax.text(0, -110, ylabel, rotation='vertical', va ='center', ha = 'center',size='11')
417 # ax.text(0, -110, ylabel, rotation='vertical', va ='center', ha = 'center',size='11')
418 printLabels(ax, xlabel, '', title)
418 printLabels(ax, xlabel, '', title)
419
419
420 set_linedata(ax, x, y, idline=0)
420 set_linedata(ax, x, y, idline=0)
421
421
422 def draw(fig):
422 def draw(fig):
423
423
424 if type(fig) == 'int':
424 if type(fig) == 'int':
425 raise ValueError, "This parameter should be of tpye matplotlib figure"
425 raise ValueError, "This parameter should be of tpye matplotlib figure"
426
426
427 fig.canvas.draw()
427 fig.canvas.draw()
Show file before | Show file after | Hide comments
@@ -1,1539 +1,1749
1 import numpy
1 import numpy
2 import math
2 import math
3 from scipy import optimize
3 from scipy import optimize
4 from scipy import interpolate
4 from scipy import interpolate
5 from scipy import signal
5 from scipy import signal
6 from scipy import stats
6 from scipy import stats
7 import re
7 import re
8 import datetime
8 import datetime
9 import copy
9 import copy
10
10 import sys
11 import importlib
12 import itertools
11
13
12 from jroproc_base import ProcessingUnit, Operation
14 from jroproc_base import ProcessingUnit, Operation
13 from model.data.jrodata import Parameters
15 from model.data.jrodata import Parameters
14
16
15
17
16 class ParametersProc(ProcessingUnit):
18 class ParametersProc(ProcessingUnit):
17
19
18 nSeconds = None
20 nSeconds = None
19
21
20 def __init__(self):
22 def __init__(self):
21 ProcessingUnit.__init__(self)
23 ProcessingUnit.__init__(self)
22
24
23 self.objectDict = {}
25 self.objectDict = {}
24 self.buffer = None
26 self.buffer = None
25 self.firstdatatime = None
27 self.firstdatatime = None
26 self.profIndex = 0
28 self.profIndex = 0
27 self.dataOut = Parameters()
29 self.dataOut = Parameters()
28
30
29 def __updateObjFromInput(self):
31 def __updateObjFromInput(self):
30
32
31 self.dataOut.inputUnit = self.dataIn.type
33 self.dataOut.inputUnit = self.dataIn.type
32
34
33 self.dataOut.timeZone = self.dataIn.timeZone
35 self.dataOut.timeZone = self.dataIn.timeZone
34 self.dataOut.dstFlag = self.dataIn.dstFlag
36 self.dataOut.dstFlag = self.dataIn.dstFlag
35 self.dataOut.errorCount = self.dataIn.errorCount
37 self.dataOut.errorCount = self.dataIn.errorCount
36 self.dataOut.useLocalTime = self.dataIn.useLocalTime
38 self.dataOut.useLocalTime = self.dataIn.useLocalTime
37
39
38 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
40 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
39 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
41 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
40 self.dataOut.channelList = self.dataIn.channelList
42 self.dataOut.channelList = self.dataIn.channelList
41 self.dataOut.heightList = self.dataIn.heightList
43 self.dataOut.heightList = self.dataIn.heightList
42 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
44 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
43 # self.dataOut.nHeights = self.dataIn.nHeights
45 # self.dataOut.nHeights = self.dataIn.nHeights
44 # self.dataOut.nChannels = self.dataIn.nChannels
46 # self.dataOut.nChannels = self.dataIn.nChannels
45 self.dataOut.nBaud = self.dataIn.nBaud
47 self.dataOut.nBaud = self.dataIn.nBaud
46 self.dataOut.nCode = self.dataIn.nCode
48 self.dataOut.nCode = self.dataIn.nCode
47 self.dataOut.code = self.dataIn.code
49 self.dataOut.code = self.dataIn.code
48 # self.dataOut.nProfiles = self.dataOut.nFFTPoints
50 # self.dataOut.nProfiles = self.dataOut.nFFTPoints
49 self.dataOut.flagTimeBlock = self.dataIn.flagTimeBlock
51 self.dataOut.flagTimeBlock = self.dataIn.flagTimeBlock
50 self.dataOut.utctime = self.firstdatatime
52 self.dataOut.utctime = self.firstdatatime
51 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
53 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
52 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
54 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
53 # self.dataOut.nCohInt = self.dataIn.nCohInt
55 # self.dataOut.nCohInt = self.dataIn.nCohInt
54 # self.dataOut.nIncohInt = 1
56 # self.dataOut.nIncohInt = 1
55 self.dataOut.ippSeconds = self.dataIn.ippSeconds
57 self.dataOut.ippSeconds = self.dataIn.ippSeconds
56 # self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
58 # self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
57 self.dataOut.timeInterval = self.dataIn.timeInterval
59 self.dataOut.timeInterval = self.dataIn.timeInterval
58 self.dataOut.heightRange = self.dataIn.getHeiRange()
60 self.dataOut.heightRange = self.dataIn.getHeiRange()
59 self.dataOut.frequency = self.dataIn.frequency
61 self.dataOut.frequency = self.dataIn.frequency
60
62
61 def run(self, nSeconds = None, nProfiles = None):
63 def run(self, nSeconds = None, nProfiles = None):
62
64
63 self.dataOut.flagNoData = True
65
64
66
65 if self.firstdatatime == None:
67 if self.firstdatatime == None:
66 self.firstdatatime = self.dataIn.utctime
68 self.firstdatatime = self.dataIn.utctime
67
69
68 #---------------------- Voltage Data ---------------------------
70 #---------------------- Voltage Data ---------------------------
69
71
70 if self.dataIn.type == "Voltage":
72 if self.dataIn.type == "Voltage":
73 self.dataOut.flagNoData = True
71 if nSeconds != None:
74 if nSeconds != None:
72 self.nSeconds = nSeconds
75 self.nSeconds = nSeconds
73 self.nProfiles= int(numpy.floor(nSeconds/(self.dataIn.ippSeconds*self.dataIn.nCohInt)))
76 self.nProfiles= int(numpy.floor(nSeconds/(self.dataIn.ippSeconds*self.dataIn.nCohInt)))
74
77
75 if self.buffer == None:
78 if self.buffer == None:
76 self.buffer = numpy.zeros((self.dataIn.nChannels,
79 self.buffer = numpy.zeros((self.dataIn.nChannels,
77 self.nProfiles,
80 self.nProfiles,
78 self.dataIn.nHeights),
81 self.dataIn.nHeights),
79 dtype='complex')
82 dtype='complex')
80
83
81 self.buffer[:,self.profIndex,:] = self.dataIn.data.copy()
84 self.buffer[:,self.profIndex,:] = self.dataIn.data.copy()
82 self.profIndex += 1
85 self.profIndex += 1
83
86
84 if self.profIndex == self.nProfiles:
87 if self.profIndex == self.nProfiles:
85
88
86 self.__updateObjFromInput()
89 self.__updateObjFromInput()
87 self.dataOut.data_pre = self.buffer.copy()
90 self.dataOut.data_pre = self.buffer.copy()
88 self.dataOut.paramInterval = nSeconds
91 self.dataOut.paramInterval = nSeconds
89 self.dataOut.flagNoData = False
92 self.dataOut.flagNoData = False
90
93
91 self.buffer = None
94 self.buffer = None
92 self.firstdatatime = None
95 self.firstdatatime = None
93 self.profIndex = 0
96 self.profIndex = 0
94 return
97 return
95
98
96 #---------------------- Spectra Data ---------------------------
99 #---------------------- Spectra Data ---------------------------
97
100
98 if self.dataIn.type == "Spectra":
101 if self.dataIn.type == "Spectra":
99 self.dataOut.data_pre = self.dataIn.data_spc.copy()
102 self.dataOut.data_pre = self.dataIn.data_spc.copy()
100 self.dataOut.abscissaRange = self.dataIn.getVelRange(1)
103 self.dataOut.abscissaRange = self.dataIn.getVelRange(1)
101 self.dataOut.noise = self.dataIn.getNoise()
104 self.dataOut.noise = self.dataIn.getNoise()
102 self.dataOut.normFactor = self.dataIn.normFactor
105 self.dataOut.normFactor = self.dataIn.normFactor
106 self.dataOut.flagNoData = False
103
107
104 #---------------------- Correlation Data ---------------------------
108 #---------------------- Correlation Data ---------------------------
105
109
106 if self.dataIn.type == "Correlation":
110 if self.dataIn.type == "Correlation":
107 lagRRange = self.dataIn.lagR
111 lagRRange = self.dataIn.lagR
108 indR = numpy.where(lagRRange == 0)[0][0]
112 indR = numpy.where(lagRRange == 0)[0][0]
109
113
110 self.dataOut.data_pre = self.dataIn.data_corr.copy()[:,:,indR,:]
114 self.dataOut.data_pre = self.dataIn.data_corr.copy()[:,:,indR,:]
111 self.dataOut.abscissaRange = self.dataIn.getLagTRange(1)
115 self.dataOut.abscissaRange = self.dataIn.getLagTRange(1)
112 self.dataOut.noise = self.dataIn.noise
116 self.dataOut.noise = self.dataIn.noise
113 self.dataOut.normFactor = self.dataIn.normFactor
117 self.dataOut.normFactor = self.dataIn.normFactor
114 self.dataOut.SNR = self.dataIn.SNR
118 self.dataOut.SNR = self.dataIn.SNR
115 self.dataOut.pairsList = self.dataIn.pairsList
119 self.dataOut.groupList = self.dataIn.pairsList
120 self.dataOut.flagNoData = False
116
121
117
122
118 self.__updateObjFromInput()
123 self.__updateObjFromInput()
119 self.dataOut.flagNoData = False
120 self.firstdatatime = None
124 self.firstdatatime = None
121 self.dataOut.initUtcTime = self.dataIn.ltctime
125 self.dataOut.initUtcTime = self.dataIn.ltctime
122 self.dataOut.windsInterval = self.dataIn.timeInterval
126 self.dataOut.outputInterval = self.dataIn.timeInterval
123
127
124 #------------------- Get Moments ----------------------------------
128 #------------------- Get Moments ----------------------------------
125 def GetMoments(self, channelList = None):
129 def GetMoments(self, channelList = None):
126 '''
130 '''
127 Function GetMoments()
131 Function GetMoments()
128
132
129 Input:
133 Input:
130 channelList : simple channel list to select e.g. [2,3,7]
134 channelList : simple channel list to select e.g. [2,3,7]
131 self.dataOut.data_pre
135 self.dataOut.data_pre
132 self.dataOut.abscissaRange
136 self.dataOut.abscissaRange
133 self.dataOut.noise
137 self.dataOut.noise
134
138
135 Affected:
139 Affected:
136 self.dataOut.data_param
140 self.dataOut.data_param
137 self.dataOut.SNR
141 self.dataOut.SNR
138
142
139 '''
143 '''
140 data = self.dataOut.data_pre
144 data = self.dataOut.data_pre
141 absc = self.dataOut.abscissaRange[:-1]
145 absc = self.dataOut.abscissaRange[:-1]
142 noise = self.dataOut.noise
146 noise = self.dataOut.noise
143
147
144 data_param = numpy.zeros((data.shape[0], 4, data.shape[2]))
148 data_param = numpy.zeros((data.shape[0], 4, data.shape[2]))
145
149
146 if channelList== None:
150 if channelList== None:
147 channelList = self.dataIn.channelList
151 channelList = self.dataIn.channelList
148 self.dataOut.channelList = channelList
152 self.dataOut.channelList = channelList
149
153
150 for ind in channelList:
154 for ind in channelList:
151 data_param[ind,:,:] = self.__calculateMoments(data[ind,:,:], absc, noise[ind])
155 data_param[ind,:,:] = self.__calculateMoments(data[ind,:,:], absc, noise[ind])
152
156
153 self.dataOut.data_param = data_param[:,1:,:]
157 self.dataOut.data_param = data_param[:,1:,:]
154 self.dataOut.SNR = data_param[:,0]
158 self.dataOut.SNR = data_param[:,0]
155 return
159 return
156
160
157 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):
161 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):
158
162
159 if (nicoh == None): nicoh = 1
163 if (nicoh == None): nicoh = 1
160 if (graph == None): graph = 0
164 if (graph == None): graph = 0
161 if (smooth == None): smooth = 0
165 if (smooth == None): smooth = 0
162 elif (self.smooth < 3): smooth = 0
166 elif (self.smooth < 3): smooth = 0
163
167
164 if (type1 == None): type1 = 0
168 if (type1 == None): type1 = 0
165 if (fwindow == None): fwindow = numpy.zeros(oldfreq.size) + 1
169 if (fwindow == None): fwindow = numpy.zeros(oldfreq.size) + 1
166 if (snrth == None): snrth = -3
170 if (snrth == None): snrth = -3
167 if (dc == None): dc = 0
171 if (dc == None): dc = 0
168 if (aliasing == None): aliasing = 0
172 if (aliasing == None): aliasing = 0
169 if (oldfd == None): oldfd = 0
173 if (oldfd == None): oldfd = 0
170 if (wwauto == None): wwauto = 0
174 if (wwauto == None): wwauto = 0
171
175
172 if (n0 < 1.e-20): n0 = 1.e-20
176 if (n0 < 1.e-20): n0 = 1.e-20
173
177
174 freq = oldfreq
178 freq = oldfreq
175 vec_power = numpy.zeros(oldspec.shape[1])
179 vec_power = numpy.zeros(oldspec.shape[1])
176 vec_fd = numpy.zeros(oldspec.shape[1])
180 vec_fd = numpy.zeros(oldspec.shape[1])
177 vec_w = numpy.zeros(oldspec.shape[1])
181 vec_w = numpy.zeros(oldspec.shape[1])
178 vec_snr = numpy.zeros(oldspec.shape[1])
182 vec_snr = numpy.zeros(oldspec.shape[1])
179
183
180 for ind in range(oldspec.shape[1]):
184 for ind in range(oldspec.shape[1]):
181
185
182 spec = oldspec[:,ind]
186 spec = oldspec[:,ind]
183 aux = spec*fwindow
187 aux = spec*fwindow
184 max_spec = aux.max()
188 max_spec = aux.max()
185 m = list(aux).index(max_spec)
189 m = list(aux).index(max_spec)
186
190
187 #Smooth
191 #Smooth
188 if (smooth == 0): spec2 = spec
192 if (smooth == 0): spec2 = spec
189 else: spec2 = scipy.ndimage.filters.uniform_filter1d(spec,size=smooth)
193 else: spec2 = scipy.ndimage.filters.uniform_filter1d(spec,size=smooth)
190
194
191 # Calculo de Momentos
195 # Calculo de Momentos
192 bb = spec2[range(m,spec2.size)]
196 bb = spec2[range(m,spec2.size)]
193 bb = (bb<n0).nonzero()
197 bb = (bb<n0).nonzero()
194 bb = bb[0]
198 bb = bb[0]
195
199
196 ss = spec2[range(0,m + 1)]
200 ss = spec2[range(0,m + 1)]
197 ss = (ss<n0).nonzero()
201 ss = (ss<n0).nonzero()
198 ss = ss[0]
202 ss = ss[0]
199
203
200 if (bb.size == 0):
204 if (bb.size == 0):
201 bb0 = spec.size - 1 - m
205 bb0 = spec.size - 1 - m
202 else:
206 else:
203 bb0 = bb[0] - 1
207 bb0 = bb[0] - 1
204 if (bb0 < 0):
208 if (bb0 < 0):
205 bb0 = 0
209 bb0 = 0
206
210
207 if (ss.size == 0): ss1 = 1
211 if (ss.size == 0): ss1 = 1
208 else: ss1 = max(ss) + 1
212 else: ss1 = max(ss) + 1
209
213
210 if (ss1 > m): ss1 = m
214 if (ss1 > m): ss1 = m
211
215
212 valid = numpy.asarray(range(int(m + bb0 - ss1 + 1))) + ss1
216 valid = numpy.asarray(range(int(m + bb0 - ss1 + 1))) + ss1
213 power = ((spec2[valid] - n0)*fwindow[valid]).sum()
217 power = ((spec2[valid] - n0)*fwindow[valid]).sum()
214 fd = ((spec2[valid]- n0)*freq[valid]*fwindow[valid]).sum()/power
218 fd = ((spec2[valid]- n0)*freq[valid]*fwindow[valid]).sum()/power
215 w = math.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum()/power)
219 w = math.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum()/power)
216 snr = (spec2.mean()-n0)/n0
220 snr = (spec2.mean()-n0)/n0
217
221
218 if (snr < 1.e-20) :
222 if (snr < 1.e-20) :
219 snr = 1.e-20
223 snr = 1.e-20
220
224
221 vec_power[ind] = power
225 vec_power[ind] = power
222 vec_fd[ind] = fd
226 vec_fd[ind] = fd
223 vec_w[ind] = w
227 vec_w[ind] = w
224 vec_snr[ind] = snr
228 vec_snr[ind] = snr
225
229
226 moments = numpy.vstack((vec_snr, vec_power, vec_fd, vec_w))
230 moments = numpy.vstack((vec_snr, vec_power, vec_fd, vec_w))
227 return moments
231 return moments
228
232
229 #------------------- Get Lags ----------------------------------
233 #------------------- Get Lags ----------------------------------
230
234
231 def GetLags(self):
235 def GetLags(self):
232 '''
236 '''
233 Function GetMoments()
237 Function GetMoments()
234
238
235 Input:
239 Input:
236 self.dataOut.data_pre
240 self.dataOut.data_pre
237 self.dataOut.abscissaRange
241 self.dataOut.abscissaRange
238 self.dataOut.noise
242 self.dataOut.noise
239 self.dataOut.normFactor
243 self.dataOut.normFactor
240 self.dataOut.SNR
244 self.dataOut.SNR
241 self.dataOut.pairsList
245 self.dataOut.groupList
242 self.dataOut.nChannels
246 self.dataOut.nChannels
243
247
244 Affected:
248 Affected:
245 self.dataOut.data_param
249 self.dataOut.data_param
246
250
247 '''
251 '''
248 data = self.dataOut.data_pre
252 data = self.dataOut.data_pre
249 normFactor = self.dataOut.normFactor
253 normFactor = self.dataOut.normFactor
250 nHeights = self.dataOut.nHeights
254 nHeights = self.dataOut.nHeights
251 absc = self.dataOut.abscissaRange[:-1]
255 absc = self.dataOut.abscissaRange[:-1]
252 noise = self.dataOut.noise
256 noise = self.dataOut.noise
253 SNR = self.dataOut.SNR
257 SNR = self.dataOut.SNR
254 pairsList = self.dataOut.pairsList
258 pairsList = self.dataOut.groupList
255 nChannels = self.dataOut.nChannels
259 nChannels = self.dataOut.nChannels
256 pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)
260 pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)
257 self.dataOut.data_param = numpy.zeros((len(pairsCrossCorr)*2 + 1, nHeights))
261 self.dataOut.data_param = numpy.zeros((len(pairsCrossCorr)*2 + 1, nHeights))
258
262
259 dataNorm = numpy.abs(data)
263 dataNorm = numpy.abs(data)
260 for l in range(len(pairsList)):
264 for l in range(len(pairsList)):
261 dataNorm[l,:,:] = dataNorm[l,:,:]/normFactor[l,:]
265 dataNorm[l,:,:] = dataNorm[l,:,:]/normFactor[l,:]
262
266
263 self.dataOut.data_param[:-1,:] = self.__calculateTaus(dataNorm, pairsCrossCorr, pairsAutoCorr, absc)
267 self.dataOut.data_param[:-1,:] = self.__calculateTaus(dataNorm, pairsCrossCorr, pairsAutoCorr, absc)
264 self.dataOut.data_param[-1,:] = self.__calculateLag1Phase(data, pairsAutoCorr, absc)
268 self.dataOut.data_param[-1,:] = self.__calculateLag1Phase(data, pairsAutoCorr, absc)
265 return
269 return
266
270
267 def __getPairsAutoCorr(self, pairsList, nChannels):
271 def __getPairsAutoCorr(self, pairsList, nChannels):
268
272
269 pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
273 pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
270
274
271 for l in range(len(pairsList)):
275 for l in range(len(pairsList)):
272 firstChannel = pairsList[l][0]
276 firstChannel = pairsList[l][0]
273 secondChannel = pairsList[l][1]
277 secondChannel = pairsList[l][1]
274
278
275 #Obteniendo pares de Autocorrelacion
279 #Obteniendo pares de Autocorrelacion
276 if firstChannel == secondChannel:
280 if firstChannel == secondChannel:
277 pairsAutoCorr[firstChannel] = int(l)
281 pairsAutoCorr[firstChannel] = int(l)
278
282
279 pairsAutoCorr = pairsAutoCorr.astype(int)
283 pairsAutoCorr = pairsAutoCorr.astype(int)
280
284
281 pairsCrossCorr = range(len(pairsList))
285 pairsCrossCorr = range(len(pairsList))
282 pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
286 pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
283
287
284 return pairsAutoCorr, pairsCrossCorr
288 return pairsAutoCorr, pairsCrossCorr
285
289
286 def __calculateTaus(self, data, pairsCrossCorr, pairsAutoCorr, lagTRange):
290 def __calculateTaus(self, data, pairsCrossCorr, pairsAutoCorr, lagTRange):
287
291
288 Pt0 = data.shape[1]/2
292 Pt0 = data.shape[1]/2
289 #Funcion de Autocorrelacion
293 #Funcion de Autocorrelacion
290 dataAutoCorr = stats.nanmean(data[pairsAutoCorr,:,:], axis = 0)
294 dataAutoCorr = stats.nanmean(data[pairsAutoCorr,:,:], axis = 0)
291
295
292 #Obtencion Indice de TauCross
296 #Obtencion Indice de TauCross
293 indCross = data[pairsCrossCorr,:,:].argmax(axis = 1)
297 indCross = data[pairsCrossCorr,:,:].argmax(axis = 1)
294 #Obtencion Indice de TauAuto
298 #Obtencion Indice de TauAuto
295 indAuto = numpy.zeros(indCross.shape,dtype = 'int')
299 indAuto = numpy.zeros(indCross.shape,dtype = 'int')
296 CCValue = data[pairsCrossCorr,Pt0,:]
300 CCValue = data[pairsCrossCorr,Pt0,:]
297 for i in range(pairsCrossCorr.size):
301 for i in range(pairsCrossCorr.size):
298 indAuto[i,:] = numpy.abs(dataAutoCorr - CCValue[i,:]).argmin(axis = 0)
302 indAuto[i,:] = numpy.abs(dataAutoCorr - CCValue[i,:]).argmin(axis = 0)
299
303
300 #Obtencion de TauCross y TauAuto
304 #Obtencion de TauCross y TauAuto
301 tauCross = lagTRange[indCross]
305 tauCross = lagTRange[indCross]
302 tauAuto = lagTRange[indAuto]
306 tauAuto = lagTRange[indAuto]
303
307
304 Nan1, Nan2 = numpy.where(tauCross == lagTRange[0])
308 Nan1, Nan2 = numpy.where(tauCross == lagTRange[0])
305
309
306 tauCross[Nan1,Nan2] = numpy.nan
310 tauCross[Nan1,Nan2] = numpy.nan
307 tauAuto[Nan1,Nan2] = numpy.nan
311 tauAuto[Nan1,Nan2] = numpy.nan
308 tau = numpy.vstack((tauCross,tauAuto))
312 tau = numpy.vstack((tauCross,tauAuto))
309
313
310 return tau
314 return tau
311
315
312 def __calculateLag1Phase(self, data, pairs, lagTRange):
316 def __calculateLag1Phase(self, data, pairs, lagTRange):
313 data1 = stats.nanmean(data[pairs,:,:], axis = 0)
317 data1 = stats.nanmean(data[pairs,:,:], axis = 0)
314 lag1 = numpy.where(lagTRange == 0)[0][0] + 1
318 lag1 = numpy.where(lagTRange == 0)[0][0] + 1
315
319
316 phase = numpy.angle(data1[lag1,:])
320 phase = numpy.angle(data1[lag1,:])
317
321
318 return phase
322 return phase
319 #------------------- Detect Meteors ------------------------------
323 #------------------- Detect Meteors ------------------------------
320
324
321 def DetectMeteors(self, hei_ref = None, tauindex = 0,
325 def DetectMeteors(self, hei_ref = None, tauindex = 0,
322 predefinedPhaseShifts = None, centerReceiverIndex = 2,
326 predefinedPhaseShifts = None, centerReceiverIndex = 2,
323 cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25,
327 cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25,
324 noise_timeStep = 4, noise_multiple = 4,
328 noise_timeStep = 4, noise_multiple = 4,
325 multDet_timeLimit = 1, multDet_rangeLimit = 3,
329 multDet_timeLimit = 1, multDet_rangeLimit = 3,
326 phaseThresh = 20, SNRThresh = 8,
330 phaseThresh = 20, SNRThresh = 8,
327 hmin = 70, hmax=110, azimuth = 0) :
331 hmin = 70, hmax=110, azimuth = 0) :
328
332
329 '''
333 '''
330 Function DetectMeteors()
334 Function DetectMeteors()
331 Project developed with paper:
335 Project developed with paper:
332 HOLDSWORTH ET AL. 2004
336 HOLDSWORTH ET AL. 2004
333
337
334 Input:
338 Input:
335 self.dataOut.data_pre
339 self.dataOut.data_pre
336
340
337 centerReceiverIndex: From the channels, which is the center receiver
341 centerReceiverIndex: From the channels, which is the center receiver
338
342
339 hei_ref: Height reference for the Beacon signal extraction
343 hei_ref: Height reference for the Beacon signal extraction
340 tauindex:
344 tauindex:
341 predefinedPhaseShifts: Predefined phase offset for the voltge signals
345 predefinedPhaseShifts: Predefined phase offset for the voltge signals
342
346
343 cohDetection: Whether to user Coherent detection or not
347 cohDetection: Whether to user Coherent detection or not
344 cohDet_timeStep: Coherent Detection calculation time step
348 cohDet_timeStep: Coherent Detection calculation time step
345 cohDet_thresh: Coherent Detection phase threshold to correct phases
349 cohDet_thresh: Coherent Detection phase threshold to correct phases
346
350
347 noise_timeStep: Noise calculation time step
351 noise_timeStep: Noise calculation time step
348 noise_multiple: Noise multiple to define signal threshold
352 noise_multiple: Noise multiple to define signal threshold
349
353
350 multDet_timeLimit: Multiple Detection Removal time limit in seconds
354 multDet_timeLimit: Multiple Detection Removal time limit in seconds
351 multDet_rangeLimit: Multiple Detection Removal range limit in km
355 multDet_rangeLimit: Multiple Detection Removal range limit in km
352
356
353 phaseThresh: Maximum phase difference between receiver to be consider a meteor
357 phaseThresh: Maximum phase difference between receiver to be consider a meteor
354 SNRThresh: Minimum SNR threshold of the meteor signal to be consider a meteor
358 SNRThresh: Minimum SNR threshold of the meteor signal to be consider a meteor
355
359
356 hmin: Minimum Height of the meteor to use it in the further wind estimations
360 hmin: Minimum Height of the meteor to use it in the further wind estimations
357 hmax: Maximum Height of the meteor to use it in the further wind estimations
361 hmax: Maximum Height of the meteor to use it in the further wind estimations
358 azimuth: Azimuth angle correction
362 azimuth: Azimuth angle correction
359
363
360 Affected:
364 Affected:
361 self.dataOut.data_param
365 self.dataOut.data_param
362
366
363 Rejection Criteria (Errors):
367 Rejection Criteria (Errors):
364 0: No error; analysis OK
368 0: No error; analysis OK
365 1: SNR < SNR threshold
369 1: SNR < SNR threshold
366 2: angle of arrival (AOA) ambiguously determined
370 2: angle of arrival (AOA) ambiguously determined
367 3: AOA estimate not feasible
371 3: AOA estimate not feasible
368 4: Large difference in AOAs obtained from different antenna baselines
372 4: Large difference in AOAs obtained from different antenna baselines
369 5: echo at start or end of time series
373 5: echo at start or end of time series
370 6: echo less than 5 examples long; too short for analysis
374 6: echo less than 5 examples long; too short for analysis
371 7: echo rise exceeds 0.3s
375 7: echo rise exceeds 0.3s
372 8: echo decay time less than twice rise time
376 8: echo decay time less than twice rise time
373 9: large power level before echo
377 9: large power level before echo
374 10: large power level after echo
378 10: large power level after echo
375 11: poor fit to amplitude for estimation of decay time
379 11: poor fit to amplitude for estimation of decay time
376 12: poor fit to CCF phase variation for estimation of radial drift velocity
380 12: poor fit to CCF phase variation for estimation of radial drift velocity
377 13: height unresolvable echo: not valid height within 70 to 110 km
381 13: height unresolvable echo: not valid height within 70 to 110 km
378 14: height ambiguous echo: more then one possible height within 70 to 110 km
382 14: height ambiguous echo: more then one possible height within 70 to 110 km
379 15: radial drift velocity or projected horizontal velocity exceeds 200 m/s
383 15: radial drift velocity or projected horizontal velocity exceeds 200 m/s
380 16: oscilatory echo, indicating event most likely not an underdense echo
384 16: oscilatory echo, indicating event most likely not an underdense echo
381
385
382 17: phase difference in meteor Reestimation
386 17: phase difference in meteor Reestimation
383
387
384 Data Storage:
388 Data Storage:
385 Meteors for Wind Estimation (8):
389 Meteors for Wind Estimation (8):
386 Day Hour | Range Height
390 Day Hour | Range Height
387 Azimuth Zenith errorCosDir
391 Azimuth Zenith errorCosDir
388 VelRad errorVelRad
392 VelRad errorVelRad
389 TypeError
393 TypeError
390
394
391 '''
395 '''
392 #Get Beacon signal
396 #Get Beacon signal
393 newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
397 newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
394
398
395 if hei_ref != None:
399 if hei_ref != None:
396 newheis = numpy.where(self.dataOut.heightList>hei_ref)
400 newheis = numpy.where(self.dataOut.heightList>hei_ref)
397
401
398 heiRang = self.dataOut.getHeiRange()
402 heiRang = self.dataOut.getHeiRange()
399 #Pairs List
403 #Pairs List
400 pairslist = []
404 pairslist = []
401 nChannel = self.dataOut.nChannels
405 nChannel = self.dataOut.nChannels
402 for i in range(nChannel):
406 for i in range(nChannel):
403 if i != centerReceiverIndex:
407 if i != centerReceiverIndex:
404 pairslist.append((centerReceiverIndex,i))
408 pairslist.append((centerReceiverIndex,i))
405
409
406 #****************REMOVING HARDWARE PHASE DIFFERENCES***************
410 #****************REMOVING HARDWARE PHASE DIFFERENCES***************
407 # see if the user put in pre defined phase shifts
411 # see if the user put in pre defined phase shifts
408 voltsPShift = self.dataOut.data_pre.copy()
412 voltsPShift = self.dataOut.data_pre.copy()
409
413
410 if predefinedPhaseShifts != None:
414 if predefinedPhaseShifts != None:
411 hardwarePhaseShifts = numpy.array(predefinedPhaseShifts)*numpy.pi/180
415 hardwarePhaseShifts = numpy.array(predefinedPhaseShifts)*numpy.pi/180
412 else:
416 else:
413 #get hardware phase shifts using beacon signal
417 #get hardware phase shifts using beacon signal
414 hardwarePhaseShifts = self.__getHardwarePhaseDiff(self.dataOut.data_pre, pairslist, newheis, 10)
418 hardwarePhaseShifts = self.__getHardwarePhaseDiff(self.dataOut.data_pre, pairslist, newheis, 10)
415 hardwarePhaseShifts = numpy.insert(hardwarePhaseShifts,centerReceiverIndex,0)
419 hardwarePhaseShifts = numpy.insert(hardwarePhaseShifts,centerReceiverIndex,0)
416
420
417 voltsPShift = numpy.zeros((self.dataOut.data_pre.shape[0],self.dataOut.data_pre.shape[1],self.dataOut.data_pre.shape[2]), dtype = 'complex')
421 voltsPShift = numpy.zeros((self.dataOut.data_pre.shape[0],self.dataOut.data_pre.shape[1],self.dataOut.data_pre.shape[2]), dtype = 'complex')
418 for i in range(self.dataOut.data_pre.shape[0]):
422 for i in range(self.dataOut.data_pre.shape[0]):
419 voltsPShift[i,:,:] = self.__shiftPhase(self.dataOut.data_pre[i,:,:], hardwarePhaseShifts[i])
423 voltsPShift[i,:,:] = self.__shiftPhase(self.dataOut.data_pre[i,:,:], hardwarePhaseShifts[i])
420 #******************END OF REMOVING HARDWARE PHASE DIFFERENCES*********
424 #******************END OF REMOVING HARDWARE PHASE DIFFERENCES*********
421
425
422 #Remove DC
426 #Remove DC
423 voltsDC = numpy.mean(voltsPShift,1)
427 voltsDC = numpy.mean(voltsPShift,1)
424 voltsDC = numpy.mean(voltsDC,1)
428 voltsDC = numpy.mean(voltsDC,1)
425 for i in range(voltsDC.shape[0]):
429 for i in range(voltsDC.shape[0]):
426 voltsPShift[i] = voltsPShift[i] - voltsDC[i]
430 voltsPShift[i] = voltsPShift[i] - voltsDC[i]
427
431
428 #Don't considerate last heights, theyre used to calculate Hardware Phase Shift
432 #Don't considerate last heights, theyre used to calculate Hardware Phase Shift
429 voltsPShift = voltsPShift[:,:,:newheis[0][0]]
433 voltsPShift = voltsPShift[:,:,:newheis[0][0]]
430
434
431 #************ FIND POWER OF DATA W/COH OR NON COH DETECTION (3.4) **********
435 #************ FIND POWER OF DATA W/COH OR NON COH DETECTION (3.4) **********
432 #Coherent Detection
436 #Coherent Detection
433 if cohDetection:
437 if cohDetection:
434 #use coherent detection to get the net power
438 #use coherent detection to get the net power
435 cohDet_thresh = cohDet_thresh*numpy.pi/180
439 cohDet_thresh = cohDet_thresh*numpy.pi/180
436 voltsPShift = self.__coherentDetection(voltsPShift, cohDet_timeStep, self.dataOut.timeInterval, pairslist, cohDet_thresh)
440 voltsPShift = self.__coherentDetection(voltsPShift, cohDet_timeStep, self.dataOut.timeInterval, pairslist, cohDet_thresh)
437
441
438 #Non-coherent detection!
442 #Non-coherent detection!
439 powerNet = numpy.nansum(numpy.abs(voltsPShift[:,:,:])**2,0)
443 powerNet = numpy.nansum(numpy.abs(voltsPShift[:,:,:])**2,0)
440 #********** END OF COH/NON-COH POWER CALCULATION**********************
444 #********** END OF COH/NON-COH POWER CALCULATION**********************
441
445
442 #********** FIND THE NOISE LEVEL AND POSSIBLE METEORS ****************
446 #********** FIND THE NOISE LEVEL AND POSSIBLE METEORS ****************
443 #Get noise
447 #Get noise
444 noise, noise1 = self.__getNoise(powerNet, noise_timeStep, self.dataOut.timeInterval)
448 noise, noise1 = self.__getNoise(powerNet, noise_timeStep, self.dataOut.timeInterval)
445 # noise = self.getNoise1(powerNet, noise_timeStep, self.dataOut.timeInterval)
449 # noise = self.getNoise1(powerNet, noise_timeStep, self.dataOut.timeInterval)
446 #Get signal threshold
450 #Get signal threshold
447 signalThresh = noise_multiple*noise
451 signalThresh = noise_multiple*noise
448 #Meteor echoes detection
452 #Meteor echoes detection
449 listMeteors = self.__findMeteors(powerNet, signalThresh)
453 listMeteors = self.__findMeteors(powerNet, signalThresh)
450 #******* END OF NOISE LEVEL AND POSSIBLE METEORS CACULATION **********
454 #******* END OF NOISE LEVEL AND POSSIBLE METEORS CACULATION **********
451
455
452 #************** REMOVE MULTIPLE DETECTIONS (3.5) ***************************
456 #************** REMOVE MULTIPLE DETECTIONS (3.5) ***************************
453 #Parameters
457 #Parameters
454 heiRange = self.dataOut.getHeiRange()
458 heiRange = self.dataOut.getHeiRange()
455 rangeInterval = heiRange[1] - heiRange[0]
459 rangeInterval = heiRange[1] - heiRange[0]
456 rangeLimit = multDet_rangeLimit/rangeInterval
460 rangeLimit = multDet_rangeLimit/rangeInterval
457 timeLimit = multDet_timeLimit/self.dataOut.timeInterval
461 timeLimit = multDet_timeLimit/self.dataOut.timeInterval
458 #Multiple detection removals
462 #Multiple detection removals
459 listMeteors1 = self.__removeMultipleDetections(listMeteors, rangeLimit, timeLimit)
463 listMeteors1 = self.__removeMultipleDetections(listMeteors, rangeLimit, timeLimit)
460 #************ END OF REMOVE MULTIPLE DETECTIONS **********************
464 #************ END OF REMOVE MULTIPLE DETECTIONS **********************
461
465
462 #********************* METEOR REESTIMATION (3.7, 3.8, 3.9, 3.10) ********************
466 #********************* METEOR REESTIMATION (3.7, 3.8, 3.9, 3.10) ********************
463 #Parameters
467 #Parameters
464 phaseThresh = phaseThresh*numpy.pi/180
468 phaseThresh = phaseThresh*numpy.pi/180
465 thresh = [phaseThresh, noise_multiple, SNRThresh]
469 thresh = [phaseThresh, noise_multiple, SNRThresh]
466 #Meteor reestimation (Errors N 1, 6, 12, 17)
470 #Meteor reestimation (Errors N 1, 6, 12, 17)
467 listMeteors2, listMeteorsPower, listMeteorsVolts = self.__meteorReestimation(listMeteors1, voltsPShift, pairslist, thresh, noise, self.dataOut.timeInterval, self.dataOut.frequency)
471 listMeteors2, listMeteorsPower, listMeteorsVolts = self.__meteorReestimation(listMeteors1, voltsPShift, pairslist, thresh, noise, self.dataOut.timeInterval, self.dataOut.frequency)
468 # listMeteors2, listMeteorsPower, listMeteorsVolts = self.meteorReestimation3(listMeteors2, listMeteorsPower, listMeteorsVolts, voltsPShift, pairslist, thresh, noise)
472 # listMeteors2, listMeteorsPower, listMeteorsVolts = self.meteorReestimation3(listMeteors2, listMeteorsPower, listMeteorsVolts, voltsPShift, pairslist, thresh, noise)
469 #Estimation of decay times (Errors N 7, 8, 11)
473 #Estimation of decay times (Errors N 7, 8, 11)
470 listMeteors3 = self.__estimateDecayTime(listMeteors2, listMeteorsPower, self.dataOut.timeInterval, self.dataOut.frequency)
474 listMeteors3 = self.__estimateDecayTime(listMeteors2, listMeteorsPower, self.dataOut.timeInterval, self.dataOut.frequency)
471 #******************* END OF METEOR REESTIMATION *******************
475 #******************* END OF METEOR REESTIMATION *******************
472
476
473 #********************* METEOR PARAMETERS CALCULATION (3.11, 3.12, 3.13) **************************
477 #********************* METEOR PARAMETERS CALCULATION (3.11, 3.12, 3.13) **************************
474 #Calculating Radial Velocity (Error N 15)
478 #Calculating Radial Velocity (Error N 15)
475 radialStdThresh = 10
479 radialStdThresh = 10
476 listMeteors4 = self.__getRadialVelocity(listMeteors3, listMeteorsVolts, radialStdThresh, pairslist, self.dataOut.timeInterval)
480 listMeteors4 = self.__getRadialVelocity(listMeteors3, listMeteorsVolts, radialStdThresh, pairslist, self.dataOut.timeInterval)
477
481
478 if len(listMeteors4) > 0:
482 if len(listMeteors4) > 0:
479 #Setting New Array
483 #Setting New Array
480 date = repr(self.dataOut.datatime)
484 date = repr(self.dataOut.datatime)
481 arrayMeteors4, arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang)
485 arrayMeteors4, arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang)
482
486
483 #Calculate AOA (Error N 3, 4)
487 #Calculate AOA (Error N 3, 4)
484 #JONES ET AL. 1998
488 #JONES ET AL. 1998
485 AOAthresh = numpy.pi/8
489 AOAthresh = numpy.pi/8
486 error = arrayParameters[:,-1]
490 error = arrayParameters[:,-1]
487 phases = -arrayMeteors4[:,9:13]
491 phases = -arrayMeteors4[:,9:13]
488 pairsList = []
492 pairsList = []
489 pairsList.append((0,3))
493 pairsList.append((0,3))
490 pairsList.append((1,2))
494 pairsList.append((1,2))
491 arrayParameters[:,4:7], arrayParameters[:,-1] = self.__getAOA(phases, pairsList, error, AOAthresh, azimuth)
495 arrayParameters[:,4:7], arrayParameters[:,-1] = self.__getAOA(phases, pairsList, error, AOAthresh, azimuth)
492
496
493 #Calculate Heights (Error N 13 and 14)
497 #Calculate Heights (Error N 13 and 14)
494 error = arrayParameters[:,-1]
498 error = arrayParameters[:,-1]
495 Ranges = arrayParameters[:,2]
499 Ranges = arrayParameters[:,2]
496 zenith = arrayParameters[:,5]
500 zenith = arrayParameters[:,5]
497 arrayParameters[:,3], arrayParameters[:,-1] = self.__getHeights(Ranges, zenith, error, hmin, hmax)
501 arrayParameters[:,3], arrayParameters[:,-1] = self.__getHeights(Ranges, zenith, error, hmin, hmax)
498 #********************* END OF PARAMETERS CALCULATION **************************
502 #********************* END OF PARAMETERS CALCULATION **************************
499
503
500 #***************************+ SAVE DATA IN HDF5 FORMAT **********************
504 #***************************+ SAVE DATA IN HDF5 FORMAT **********************
501 self.dataOut.data_param = arrayParameters
505 self.dataOut.data_param = arrayParameters
502
506
503 return
507 return
504
508
505 def __getHardwarePhaseDiff(self, voltage0, pairslist, newheis, n):
509 def __getHardwarePhaseDiff(self, voltage0, pairslist, newheis, n):
506
510
507 minIndex = min(newheis[0])
511 minIndex = min(newheis[0])
508 maxIndex = max(newheis[0])
512 maxIndex = max(newheis[0])
509
513
510 voltage = voltage0[:,:,minIndex:maxIndex+1]
514 voltage = voltage0[:,:,minIndex:maxIndex+1]
511 nLength = voltage.shape[1]/n
515 nLength = voltage.shape[1]/n
512 nMin = 0
516 nMin = 0
513 nMax = 0
517 nMax = 0
514 phaseOffset = numpy.zeros((len(pairslist),n))
518 phaseOffset = numpy.zeros((len(pairslist),n))
515
519
516 for i in range(n):
520 for i in range(n):
517 nMax += nLength
521 nMax += nLength
518 phaseCCF = -numpy.angle(self.__calculateCCF(voltage[:,nMin:nMax,:], pairslist, [0]))
522 phaseCCF = -numpy.angle(self.__calculateCCF(voltage[:,nMin:nMax,:], pairslist, [0]))
519 phaseCCF = numpy.mean(phaseCCF, axis = 2)
523 phaseCCF = numpy.mean(phaseCCF, axis = 2)
520 phaseOffset[:,i] = phaseCCF.transpose()
524 phaseOffset[:,i] = phaseCCF.transpose()
521 nMin = nMax
525 nMin = nMax
522 # phaseDiff, phaseArrival = self.estimatePhaseDifference(voltage, pairslist)
526 # phaseDiff, phaseArrival = self.estimatePhaseDifference(voltage, pairslist)
523
527
524 #Remove Outliers
528 #Remove Outliers
525 factor = 2
529 factor = 2
526 wt = phaseOffset - signal.medfilt(phaseOffset,(1,5))
530 wt = phaseOffset - signal.medfilt(phaseOffset,(1,5))
527 dw = numpy.std(wt,axis = 1)
531 dw = numpy.std(wt,axis = 1)
528 dw = dw.reshape((dw.size,1))
532 dw = dw.reshape((dw.size,1))
529 ind = numpy.where(numpy.logical_or(wt>dw*factor,wt<-dw*factor))
533 ind = numpy.where(numpy.logical_or(wt>dw*factor,wt<-dw*factor))
530 phaseOffset[ind] = numpy.nan
534 phaseOffset[ind] = numpy.nan
531 phaseOffset = stats.nanmean(phaseOffset, axis=1)
535 phaseOffset = stats.nanmean(phaseOffset, axis=1)
532
536
533 return phaseOffset
537 return phaseOffset
534
538
535 def __shiftPhase(self, data, phaseShift):
539 def __shiftPhase(self, data, phaseShift):
536 #this will shift the phase of a complex number
540 #this will shift the phase of a complex number
537 dataShifted = numpy.abs(data) * numpy.exp((numpy.angle(data)+phaseShift)*1j)
541 dataShifted = numpy.abs(data) * numpy.exp((numpy.angle(data)+phaseShift)*1j)
538 return dataShifted
542 return dataShifted
539
543
540 def __estimatePhaseDifference(self, array, pairslist):
544 def __estimatePhaseDifference(self, array, pairslist):
541 nChannel = array.shape[0]
545 nChannel = array.shape[0]
542 nHeights = array.shape[2]
546 nHeights = array.shape[2]
543 numPairs = len(pairslist)
547 numPairs = len(pairslist)
544 # phaseCCF = numpy.zeros((nChannel, 5, nHeights))
548 # phaseCCF = numpy.zeros((nChannel, 5, nHeights))
545 phaseCCF = numpy.angle(self.__calculateCCF(array, pairslist, [-2,-1,0,1,2]))
549 phaseCCF = numpy.angle(self.__calculateCCF(array, pairslist, [-2,-1,0,1,2]))
546
550
547 #Correct phases
551 #Correct phases
548 derPhaseCCF = phaseCCF[:,1:,:] - phaseCCF[:,0:-1,:]
552 derPhaseCCF = phaseCCF[:,1:,:] - phaseCCF[:,0:-1,:]
549 indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
553 indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
550
554
551 if indDer[0].shape[0] > 0:
555 if indDer[0].shape[0] > 0:
552 for i in range(indDer[0].shape[0]):
556 for i in range(indDer[0].shape[0]):
553 signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i],indDer[2][i]])
557 signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i],indDer[2][i]])
554 phaseCCF[indDer[0][i],indDer[1][i]+1:,:] += signo*2*numpy.pi
558 phaseCCF[indDer[0][i],indDer[1][i]+1:,:] += signo*2*numpy.pi
555
559
556 # for j in range(numSides):
560 # for j in range(numSides):
557 # phaseCCFAux = self.calculateCCF(arrayCenter, arraySides[j,:,:], [-2,1,0,1,2])
561 # phaseCCFAux = self.calculateCCF(arrayCenter, arraySides[j,:,:], [-2,1,0,1,2])
558 # phaseCCF[j,:,:] = numpy.angle(phaseCCFAux)
562 # phaseCCF[j,:,:] = numpy.angle(phaseCCFAux)
559 #
563 #
560 #Linear
564 #Linear
561 phaseInt = numpy.zeros((numPairs,1))
565 phaseInt = numpy.zeros((numPairs,1))
562 angAllCCF = phaseCCF[:,[0,1,3,4],0]
566 angAllCCF = phaseCCF[:,[0,1,3,4],0]
563 for j in range(numPairs):
567 for j in range(numPairs):
564 fit = stats.linregress([-2,-1,1,2],angAllCCF[j,:])
568 fit = stats.linregress([-2,-1,1,2],angAllCCF[j,:])
565 phaseInt[j] = fit[1]
569 phaseInt[j] = fit[1]
566 #Phase Differences
570 #Phase Differences
567 phaseDiff = phaseInt - phaseCCF[:,2,:]
571 phaseDiff = phaseInt - phaseCCF[:,2,:]
568 phaseArrival = phaseInt.reshape(phaseInt.size)
572 phaseArrival = phaseInt.reshape(phaseInt.size)
569
573
570 #Dealias
574 #Dealias
571 indAlias = numpy.where(phaseArrival > numpy.pi)
575 indAlias = numpy.where(phaseArrival > numpy.pi)
572 phaseArrival[indAlias] -= 2*numpy.pi
576 phaseArrival[indAlias] -= 2*numpy.pi
573 indAlias = numpy.where(phaseArrival < -numpy.pi)
577 indAlias = numpy.where(phaseArrival < -numpy.pi)
574 phaseArrival[indAlias] += 2*numpy.pi
578 phaseArrival[indAlias] += 2*numpy.pi
575
579
576 return phaseDiff, phaseArrival
580 return phaseDiff, phaseArrival
577
581
578 def __coherentDetection(self, volts, timeSegment, timeInterval, pairslist, thresh):
582 def __coherentDetection(self, volts, timeSegment, timeInterval, pairslist, thresh):
579 #this function will run the coherent detection used in Holdworth et al. 2004 and return the net power
583 #this function will run the coherent detection used in Holdworth et al. 2004 and return the net power
580 #find the phase shifts of each channel over 1 second intervals
584 #find the phase shifts of each channel over 1 second intervals
581 #only look at ranges below the beacon signal
585 #only look at ranges below the beacon signal
582 numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
586 numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
583 numBlocks = int(volts.shape[1]/numProfPerBlock)
587 numBlocks = int(volts.shape[1]/numProfPerBlock)
584 numHeights = volts.shape[2]
588 numHeights = volts.shape[2]
585 nChannel = volts.shape[0]
589 nChannel = volts.shape[0]
586 voltsCohDet = volts.copy()
590 voltsCohDet = volts.copy()
587
591
588 pairsarray = numpy.array(pairslist)
592 pairsarray = numpy.array(pairslist)
589 indSides = pairsarray[:,1]
593 indSides = pairsarray[:,1]
590 # indSides = numpy.array(range(nChannel))
594 # indSides = numpy.array(range(nChannel))
591 # indSides = numpy.delete(indSides, indCenter)
595 # indSides = numpy.delete(indSides, indCenter)
592 #
596 #
593 # listCenter = numpy.array_split(volts[indCenter,:,:], numBlocks, 0)
597 # listCenter = numpy.array_split(volts[indCenter,:,:], numBlocks, 0)
594 listBlocks = numpy.array_split(volts, numBlocks, 1)
598 listBlocks = numpy.array_split(volts, numBlocks, 1)
595
599
596 startInd = 0
600 startInd = 0
597 endInd = 0
601 endInd = 0
598
602
599 for i in range(numBlocks):
603 for i in range(numBlocks):
600 startInd = endInd
604 startInd = endInd
601 endInd = endInd + listBlocks[i].shape[1]
605 endInd = endInd + listBlocks[i].shape[1]
602
606
603 arrayBlock = listBlocks[i]
607 arrayBlock = listBlocks[i]
604 # arrayBlockCenter = listCenter[i]
608 # arrayBlockCenter = listCenter[i]
605
609
606 #Estimate the Phase Difference
610 #Estimate the Phase Difference
607 phaseDiff, aux = self.__estimatePhaseDifference(arrayBlock, pairslist)
611 phaseDiff, aux = self.__estimatePhaseDifference(arrayBlock, pairslist)
608 #Phase Difference RMS
612 #Phase Difference RMS
609 arrayPhaseRMS = numpy.abs(phaseDiff)
613 arrayPhaseRMS = numpy.abs(phaseDiff)
610 phaseRMSaux = numpy.sum(arrayPhaseRMS < thresh,0)
614 phaseRMSaux = numpy.sum(arrayPhaseRMS < thresh,0)
611 indPhase = numpy.where(phaseRMSaux==4)
615 indPhase = numpy.where(phaseRMSaux==4)
612 #Shifting
616 #Shifting
613 if indPhase[0].shape[0] > 0:
617 if indPhase[0].shape[0] > 0:
614 for j in range(indSides.size):
618 for j in range(indSides.size):
615 arrayBlock[indSides[j],:,indPhase] = self.__shiftPhase(arrayBlock[indSides[j],:,indPhase], phaseDiff[j,indPhase].transpose())
619 arrayBlock[indSides[j],:,indPhase] = self.__shiftPhase(arrayBlock[indSides[j],:,indPhase], phaseDiff[j,indPhase].transpose())
616 voltsCohDet[:,startInd:endInd,:] = arrayBlock
620 voltsCohDet[:,startInd:endInd,:] = arrayBlock
617
621
618 return voltsCohDet
622 return voltsCohDet
619
623
620 def __calculateCCF(self, volts, pairslist ,laglist):
624 def __calculateCCF(self, volts, pairslist ,laglist):
621
625
622 nHeights = volts.shape[2]
626 nHeights = volts.shape[2]
623 nPoints = volts.shape[1]
627 nPoints = volts.shape[1]
624 voltsCCF = numpy.zeros((len(pairslist), len(laglist), nHeights),dtype = 'complex')
628 voltsCCF = numpy.zeros((len(pairslist), len(laglist), nHeights),dtype = 'complex')
625
629
626 for i in range(len(pairslist)):
630 for i in range(len(pairslist)):
627 volts1 = volts[pairslist[i][0]]
631 volts1 = volts[pairslist[i][0]]
628 volts2 = volts[pairslist[i][1]]
632 volts2 = volts[pairslist[i][1]]
629
633
630 for t in range(len(laglist)):
634 for t in range(len(laglist)):
631 idxT = laglist[t]
635 idxT = laglist[t]
632 if idxT >= 0:
636 if idxT >= 0:
633 vStacked = numpy.vstack((volts2[idxT:,:],
637 vStacked = numpy.vstack((volts2[idxT:,:],
634 numpy.zeros((idxT, nHeights),dtype='complex')))
638 numpy.zeros((idxT, nHeights),dtype='complex')))
635 else:
639 else:
636 vStacked = numpy.vstack((numpy.zeros((-idxT, nHeights),dtype='complex'),
640 vStacked = numpy.vstack((numpy.zeros((-idxT, nHeights),dtype='complex'),
637 volts2[:(nPoints + idxT),:]))
641 volts2[:(nPoints + idxT),:]))
638 voltsCCF[i,t,:] = numpy.sum((numpy.conjugate(volts1)*vStacked),axis=0)
642 voltsCCF[i,t,:] = numpy.sum((numpy.conjugate(volts1)*vStacked),axis=0)
639
643
640 vStacked = None
644 vStacked = None
641 return voltsCCF
645 return voltsCCF
642
646
643 def __getNoise(self, power, timeSegment, timeInterval):
647 def __getNoise(self, power, timeSegment, timeInterval):
644 numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
648 numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
645 numBlocks = int(power.shape[0]/numProfPerBlock)
649 numBlocks = int(power.shape[0]/numProfPerBlock)
646 numHeights = power.shape[1]
650 numHeights = power.shape[1]
647
651
648 listPower = numpy.array_split(power, numBlocks, 0)
652 listPower = numpy.array_split(power, numBlocks, 0)
649 noise = numpy.zeros((power.shape[0], power.shape[1]))
653 noise = numpy.zeros((power.shape[0], power.shape[1]))
650 noise1 = numpy.zeros((power.shape[0], power.shape[1]))
654 noise1 = numpy.zeros((power.shape[0], power.shape[1]))
651
655
652 startInd = 0
656 startInd = 0
653 endInd = 0
657 endInd = 0
654
658
655 for i in range(numBlocks): #split por canal
659 for i in range(numBlocks): #split por canal
656 startInd = endInd
660 startInd = endInd
657 endInd = endInd + listPower[i].shape[0]
661 endInd = endInd + listPower[i].shape[0]
658
662
659 arrayBlock = listPower[i]
663 arrayBlock = listPower[i]
660 noiseAux = numpy.mean(arrayBlock, 0)
664 noiseAux = numpy.mean(arrayBlock, 0)
661 # noiseAux = numpy.median(noiseAux)
665 # noiseAux = numpy.median(noiseAux)
662 # noiseAux = numpy.mean(arrayBlock)
666 # noiseAux = numpy.mean(arrayBlock)
663 noise[startInd:endInd,:] = noise[startInd:endInd,:] + noiseAux
667 noise[startInd:endInd,:] = noise[startInd:endInd,:] + noiseAux
664
668
665 noiseAux1 = numpy.mean(arrayBlock)
669 noiseAux1 = numpy.mean(arrayBlock)
666 noise1[startInd:endInd,:] = noise1[startInd:endInd,:] + noiseAux1
670 noise1[startInd:endInd,:] = noise1[startInd:endInd,:] + noiseAux1
667
671
668 return noise, noise1
672 return noise, noise1
669
673
670 def __findMeteors(self, power, thresh):
674 def __findMeteors(self, power, thresh):
671 nProf = power.shape[0]
675 nProf = power.shape[0]
672 nHeights = power.shape[1]
676 nHeights = power.shape[1]
673 listMeteors = []
677 listMeteors = []
674
678
675 for i in range(nHeights):
679 for i in range(nHeights):
676 powerAux = power[:,i]
680 powerAux = power[:,i]
677 threshAux = thresh[:,i]
681 threshAux = thresh[:,i]
678
682
679 indUPthresh = numpy.where(powerAux > threshAux)[0]
683 indUPthresh = numpy.where(powerAux > threshAux)[0]
680 indDNthresh = numpy.where(powerAux <= threshAux)[0]
684 indDNthresh = numpy.where(powerAux <= threshAux)[0]
681
685
682 j = 0
686 j = 0
683
687
684 while (j < indUPthresh.size - 2):
688 while (j < indUPthresh.size - 2):
685 if (indUPthresh[j + 2] == indUPthresh[j] + 2):
689 if (indUPthresh[j + 2] == indUPthresh[j] + 2):
686 indDNAux = numpy.where(indDNthresh > indUPthresh[j])
690 indDNAux = numpy.where(indDNthresh > indUPthresh[j])
687 indDNthresh = indDNthresh[indDNAux]
691 indDNthresh = indDNthresh[indDNAux]
688
692
689 if (indDNthresh.size > 0):
693 if (indDNthresh.size > 0):
690 indEnd = indDNthresh[0] - 1
694 indEnd = indDNthresh[0] - 1
691 indInit = indUPthresh[j]
695 indInit = indUPthresh[j]
692
696
693 meteor = powerAux[indInit:indEnd + 1]
697 meteor = powerAux[indInit:indEnd + 1]
694 indPeak = meteor.argmax() + indInit
698 indPeak = meteor.argmax() + indInit
695 FLA = sum(numpy.conj(meteor)*numpy.hstack((meteor[1:],0)))
699 FLA = sum(numpy.conj(meteor)*numpy.hstack((meteor[1:],0)))
696
700
697 listMeteors.append(numpy.array([i,indInit,indPeak,indEnd,FLA])) #CHEQUEAR!!!!!
701 listMeteors.append(numpy.array([i,indInit,indPeak,indEnd,FLA])) #CHEQUEAR!!!!!
698 j = numpy.where(indUPthresh == indEnd)[0] + 1
702 j = numpy.where(indUPthresh == indEnd)[0] + 1
699 else: j+=1
703 else: j+=1
700 else: j+=1
704 else: j+=1
701
705
702 return listMeteors
706 return listMeteors
703
707
704 def __removeMultipleDetections(self,listMeteors, rangeLimit, timeLimit):
708 def __removeMultipleDetections(self,listMeteors, rangeLimit, timeLimit):
705
709
706 arrayMeteors = numpy.asarray(listMeteors)
710 arrayMeteors = numpy.asarray(listMeteors)
707 listMeteors1 = []
711 listMeteors1 = []
708
712
709 while arrayMeteors.shape[0] > 0:
713 while arrayMeteors.shape[0] > 0:
710 FLAs = arrayMeteors[:,4]
714 FLAs = arrayMeteors[:,4]
711 maxFLA = FLAs.argmax()
715 maxFLA = FLAs.argmax()
712 listMeteors1.append(arrayMeteors[maxFLA,:])
716 listMeteors1.append(arrayMeteors[maxFLA,:])
713
717
714 MeteorInitTime = arrayMeteors[maxFLA,1]
718 MeteorInitTime = arrayMeteors[maxFLA,1]
715 MeteorEndTime = arrayMeteors[maxFLA,3]
719 MeteorEndTime = arrayMeteors[maxFLA,3]
716 MeteorHeight = arrayMeteors[maxFLA,0]
720 MeteorHeight = arrayMeteors[maxFLA,0]
717
721
718 #Check neighborhood
722 #Check neighborhood
719 maxHeightIndex = MeteorHeight + rangeLimit
723 maxHeightIndex = MeteorHeight + rangeLimit
720 minHeightIndex = MeteorHeight - rangeLimit
724 minHeightIndex = MeteorHeight - rangeLimit
721 minTimeIndex = MeteorInitTime - timeLimit
725 minTimeIndex = MeteorInitTime - timeLimit
722 maxTimeIndex = MeteorEndTime + timeLimit
726 maxTimeIndex = MeteorEndTime + timeLimit
723
727
724 #Check Heights
728 #Check Heights
725 indHeight = numpy.logical_and(arrayMeteors[:,0] >= minHeightIndex, arrayMeteors[:,0] <= maxHeightIndex)
729 indHeight = numpy.logical_and(arrayMeteors[:,0] >= minHeightIndex, arrayMeteors[:,0] <= maxHeightIndex)
726 indTime = numpy.logical_and(arrayMeteors[:,3] >= minTimeIndex, arrayMeteors[:,1] <= maxTimeIndex)
730 indTime = numpy.logical_and(arrayMeteors[:,3] >= minTimeIndex, arrayMeteors[:,1] <= maxTimeIndex)
727 indBoth = numpy.where(numpy.logical_and(indTime,indHeight))
731 indBoth = numpy.where(numpy.logical_and(indTime,indHeight))
728
732
729 arrayMeteors = numpy.delete(arrayMeteors, indBoth, axis = 0)
733 arrayMeteors = numpy.delete(arrayMeteors, indBoth, axis = 0)
730
734
731 return listMeteors1
735 return listMeteors1
732
736
733 def __meteorReestimation(self, listMeteors, volts, pairslist, thresh, noise, timeInterval,frequency):
737 def __meteorReestimation(self, listMeteors, volts, pairslist, thresh, noise, timeInterval,frequency):
734 numHeights = volts.shape[2]
738 numHeights = volts.shape[2]
735 nChannel = volts.shape[0]
739 nChannel = volts.shape[0]
736
740
737 thresholdPhase = thresh[0]
741 thresholdPhase = thresh[0]
738 thresholdNoise = thresh[1]
742 thresholdNoise = thresh[1]
739 thresholdDB = float(thresh[2])
743 thresholdDB = float(thresh[2])
740
744
741 thresholdDB1 = 10**(thresholdDB/10)
745 thresholdDB1 = 10**(thresholdDB/10)
742 pairsarray = numpy.array(pairslist)
746 pairsarray = numpy.array(pairslist)
743 indSides = pairsarray[:,1]
747 indSides = pairsarray[:,1]
744
748
745 pairslist1 = list(pairslist)
749 pairslist1 = list(pairslist)
746 pairslist1.append((0,1))
750 pairslist1.append((0,1))
747 pairslist1.append((3,4))
751 pairslist1.append((3,4))
748
752
749 listMeteors1 = []
753 listMeteors1 = []
750 listPowerSeries = []
754 listPowerSeries = []
751 listVoltageSeries = []
755 listVoltageSeries = []
752 #volts has the war data
756 #volts has the war data
753
757
754 if frequency == 30e6:
758 if frequency == 30e6:
755 timeLag = 45*10**-3
759 timeLag = 45*10**-3
756 else:
760 else:
757 timeLag = 15*10**-3
761 timeLag = 15*10**-3
758 lag = numpy.ceil(timeLag/timeInterval)
762 lag = numpy.ceil(timeLag/timeInterval)
759
763
760 for i in range(len(listMeteors)):
764 for i in range(len(listMeteors)):
761
765
762 ###################### 3.6 - 3.7 PARAMETERS REESTIMATION #########################
766 ###################### 3.6 - 3.7 PARAMETERS REESTIMATION #########################
763 meteorAux = numpy.zeros(16)
767 meteorAux = numpy.zeros(16)
764
768
765 #Loading meteor Data (mHeight, mStart, mPeak, mEnd)
769 #Loading meteor Data (mHeight, mStart, mPeak, mEnd)
766 mHeight = listMeteors[i][0]
770 mHeight = listMeteors[i][0]
767 mStart = listMeteors[i][1]
771 mStart = listMeteors[i][1]
768 mPeak = listMeteors[i][2]
772 mPeak = listMeteors[i][2]
769 mEnd = listMeteors[i][3]
773 mEnd = listMeteors[i][3]
770
774
771 #get the volt data between the start and end times of the meteor
775 #get the volt data between the start and end times of the meteor
772 meteorVolts = volts[:,mStart:mEnd+1,mHeight]
776 meteorVolts = volts[:,mStart:mEnd+1,mHeight]
773 meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
777 meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
774
778
775 #3.6. Phase Difference estimation
779 #3.6. Phase Difference estimation
776 phaseDiff, aux = self.__estimatePhaseDifference(meteorVolts, pairslist)
780 phaseDiff, aux = self.__estimatePhaseDifference(meteorVolts, pairslist)
777
781
778 #3.7. Phase difference removal & meteor start, peak and end times reestimated
782 #3.7. Phase difference removal & meteor start, peak and end times reestimated
779 #meteorVolts0.- all Channels, all Profiles
783 #meteorVolts0.- all Channels, all Profiles
780 meteorVolts0 = volts[:,:,mHeight]
784 meteorVolts0 = volts[:,:,mHeight]
781 meteorThresh = noise[:,mHeight]*thresholdNoise
785 meteorThresh = noise[:,mHeight]*thresholdNoise
782 meteorNoise = noise[:,mHeight]
786 meteorNoise = noise[:,mHeight]
783 meteorVolts0[indSides,:] = self.__shiftPhase(meteorVolts0[indSides,:], phaseDiff) #Phase Shifting
787 meteorVolts0[indSides,:] = self.__shiftPhase(meteorVolts0[indSides,:], phaseDiff) #Phase Shifting
784 powerNet0 = numpy.nansum(numpy.abs(meteorVolts0)**2, axis = 0) #Power
788 powerNet0 = numpy.nansum(numpy.abs(meteorVolts0)**2, axis = 0) #Power
785
789
786 #Times reestimation
790 #Times reestimation
787 mStart1 = numpy.where(powerNet0[:mPeak] < meteorThresh[:mPeak])[0]
791 mStart1 = numpy.where(powerNet0[:mPeak] < meteorThresh[:mPeak])[0]
788 if mStart1.size > 0:
792 if mStart1.size > 0:
789 mStart1 = mStart1[-1] + 1
793 mStart1 = mStart1[-1] + 1
790
794
791 else:
795 else:
792 mStart1 = mPeak
796 mStart1 = mPeak
793
797
794 mEnd1 = numpy.where(powerNet0[mPeak:] < meteorThresh[mPeak:])[0][0] + mPeak - 1
798 mEnd1 = numpy.where(powerNet0[mPeak:] < meteorThresh[mPeak:])[0][0] + mPeak - 1
795 mEndDecayTime1 = numpy.where(powerNet0[mPeak:] < meteorNoise[mPeak:])[0]
799 mEndDecayTime1 = numpy.where(powerNet0[mPeak:] < meteorNoise[mPeak:])[0]
796 if mEndDecayTime1.size == 0:
800 if mEndDecayTime1.size == 0:
797 mEndDecayTime1 = powerNet0.size
801 mEndDecayTime1 = powerNet0.size
798 else:
802 else:
799 mEndDecayTime1 = mEndDecayTime1[0] + mPeak - 1
803 mEndDecayTime1 = mEndDecayTime1[0] + mPeak - 1
800 # mPeak1 = meteorVolts0[mStart1:mEnd1 + 1].argmax()
804 # mPeak1 = meteorVolts0[mStart1:mEnd1 + 1].argmax()
801
805
802 #meteorVolts1.- all Channels, from start to end
806 #meteorVolts1.- all Channels, from start to end
803 meteorVolts1 = meteorVolts0[:,mStart1:mEnd1 + 1]
807 meteorVolts1 = meteorVolts0[:,mStart1:mEnd1 + 1]
804 meteorVolts2 = meteorVolts0[:,mPeak + lag:mEnd1 + 1]
808 meteorVolts2 = meteorVolts0[:,mPeak + lag:mEnd1 + 1]
805 if meteorVolts2.shape[1] == 0:
809 if meteorVolts2.shape[1] == 0:
806 meteorVolts2 = meteorVolts0[:,mPeak:mEnd1 + 1]
810 meteorVolts2 = meteorVolts0[:,mPeak:mEnd1 + 1]
807 meteorVolts1 = meteorVolts1.reshape(meteorVolts1.shape[0], meteorVolts1.shape[1], 1)
811 meteorVolts1 = meteorVolts1.reshape(meteorVolts1.shape[0], meteorVolts1.shape[1], 1)
808 meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1], 1)
812 meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1], 1)
809 ##################### END PARAMETERS REESTIMATION #########################
813 ##################### END PARAMETERS REESTIMATION #########################
810
814
811 ##################### 3.8 PHASE DIFFERENCE REESTIMATION ########################
815 ##################### 3.8 PHASE DIFFERENCE REESTIMATION ########################
812 # if mEnd1 - mStart1 > 4: #Error Number 6: echo less than 5 samples long; too short for analysis
816 # if mEnd1 - mStart1 > 4: #Error Number 6: echo less than 5 samples long; too short for analysis
813 if meteorVolts2.shape[1] > 0:
817 if meteorVolts2.shape[1] > 0:
814 #Phase Difference re-estimation
818 #Phase Difference re-estimation
815 phaseDiff1, phaseDiffint = self.__estimatePhaseDifference(meteorVolts2, pairslist1) #Phase Difference Estimation
819 phaseDiff1, phaseDiffint = self.__estimatePhaseDifference(meteorVolts2, pairslist1) #Phase Difference Estimation
816 # phaseDiff1, phaseDiffint = self.estimatePhaseDifference(meteorVolts2, pairslist)
820 # phaseDiff1, phaseDiffint = self.estimatePhaseDifference(meteorVolts2, pairslist)
817 meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1])
821 meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1])
818 phaseDiff11 = numpy.reshape(phaseDiff1, (phaseDiff1.shape[0],1))
822 phaseDiff11 = numpy.reshape(phaseDiff1, (phaseDiff1.shape[0],1))
819 meteorVolts2[indSides,:] = self.__shiftPhase(meteorVolts2[indSides,:], phaseDiff11[0:4]) #Phase Shifting
823 meteorVolts2[indSides,:] = self.__shiftPhase(meteorVolts2[indSides,:], phaseDiff11[0:4]) #Phase Shifting
820
824
821 #Phase Difference RMS
825 #Phase Difference RMS
822 phaseRMS1 = numpy.sqrt(numpy.mean(numpy.square(phaseDiff1)))
826 phaseRMS1 = numpy.sqrt(numpy.mean(numpy.square(phaseDiff1)))
823 powerNet1 = numpy.nansum(numpy.abs(meteorVolts1[:,:])**2,0)
827 powerNet1 = numpy.nansum(numpy.abs(meteorVolts1[:,:])**2,0)
824 #Data from Meteor
828 #Data from Meteor
825 mPeak1 = powerNet1.argmax() + mStart1
829 mPeak1 = powerNet1.argmax() + mStart1
826 mPeakPower1 = powerNet1.max()
830 mPeakPower1 = powerNet1.max()
827 noiseAux = sum(noise[mStart1:mEnd1 + 1,mHeight])
831 noiseAux = sum(noise[mStart1:mEnd1 + 1,mHeight])
828 mSNR1 = (sum(powerNet1)-noiseAux)/noiseAux
832 mSNR1 = (sum(powerNet1)-noiseAux)/noiseAux
829 Meteor1 = numpy.array([mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1])
833 Meteor1 = numpy.array([mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1])
830 Meteor1 = numpy.hstack((Meteor1,phaseDiffint))
834 Meteor1 = numpy.hstack((Meteor1,phaseDiffint))
831 PowerSeries = powerNet0[mStart1:mEndDecayTime1 + 1]
835 PowerSeries = powerNet0[mStart1:mEndDecayTime1 + 1]
832 #Vectorize
836 #Vectorize
833 meteorAux[0:7] = [mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1]
837 meteorAux[0:7] = [mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1]
834 meteorAux[7:11] = phaseDiffint[0:4]
838 meteorAux[7:11] = phaseDiffint[0:4]
835
839
836 #Rejection Criterions
840 #Rejection Criterions
837 if phaseRMS1 > thresholdPhase: #Error Number 17: Phase variation
841 if phaseRMS1 > thresholdPhase: #Error Number 17: Phase variation
838 meteorAux[-1] = 17
842 meteorAux[-1] = 17
839 elif mSNR1 < thresholdDB1: #Error Number 1: SNR < threshold dB
843 elif mSNR1 < thresholdDB1: #Error Number 1: SNR < threshold dB
840 meteorAux[-1] = 1
844 meteorAux[-1] = 1
841
845
842
846
843 else:
847 else:
844 meteorAux[0:4] = [mHeight, mStart, mPeak, mEnd]
848 meteorAux[0:4] = [mHeight, mStart, mPeak, mEnd]
845 meteorAux[-1] = 6 #Error Number 6: echo less than 5 samples long; too short for analysis
849 meteorAux[-1] = 6 #Error Number 6: echo less than 5 samples long; too short for analysis
846 PowerSeries = 0
850 PowerSeries = 0
847
851
848 listMeteors1.append(meteorAux)
852 listMeteors1.append(meteorAux)
849 listPowerSeries.append(PowerSeries)
853 listPowerSeries.append(PowerSeries)
850 listVoltageSeries.append(meteorVolts1)
854 listVoltageSeries.append(meteorVolts1)
851
855
852 return listMeteors1, listPowerSeries, listVoltageSeries
856 return listMeteors1, listPowerSeries, listVoltageSeries
853
857
854 def __estimateDecayTime(self, listMeteors, listPower, timeInterval, frequency):
858 def __estimateDecayTime(self, listMeteors, listPower, timeInterval, frequency):
855
859
856 threshError = 10
860 threshError = 10
857 #Depending if it is 30 or 50 MHz
861 #Depending if it is 30 or 50 MHz
858 if frequency == 30e6:
862 if frequency == 30e6:
859 timeLag = 45*10**-3
863 timeLag = 45*10**-3
860 else:
864 else:
861 timeLag = 15*10**-3
865 timeLag = 15*10**-3
862 lag = numpy.ceil(timeLag/timeInterval)
866 lag = numpy.ceil(timeLag/timeInterval)
863
867
864 listMeteors1 = []
868 listMeteors1 = []
865
869
866 for i in range(len(listMeteors)):
870 for i in range(len(listMeteors)):
867 meteorPower = listPower[i]
871 meteorPower = listPower[i]
868 meteorAux = listMeteors[i]
872 meteorAux = listMeteors[i]
869
873
870 if meteorAux[-1] == 0:
874 if meteorAux[-1] == 0:
871
875
872 try:
876 try:
873 indmax = meteorPower.argmax()
877 indmax = meteorPower.argmax()
874 indlag = indmax + lag
878 indlag = indmax + lag
875
879
876 y = meteorPower[indlag:]
880 y = meteorPower[indlag:]
877 x = numpy.arange(0, y.size)*timeLag
881 x = numpy.arange(0, y.size)*timeLag
878
882
879 #first guess
883 #first guess
880 a = y[0]
884 a = y[0]
881 tau = timeLag
885 tau = timeLag
882 #exponential fit
886 #exponential fit
883 popt, pcov = optimize.curve_fit(self.__exponential_function, x, y, p0 = [a, tau])
887 popt, pcov = optimize.curve_fit(self.__exponential_function, x, y, p0 = [a, tau])
884 y1 = self.__exponential_function(x, *popt)
888 y1 = self.__exponential_function(x, *popt)
885 #error estimation
889 #error estimation
886 error = sum((y - y1)**2)/(numpy.var(y)*(y.size - popt.size))
890 error = sum((y - y1)**2)/(numpy.var(y)*(y.size - popt.size))
887
891
888 decayTime = popt[1]
892 decayTime = popt[1]
889 riseTime = indmax*timeInterval
893 riseTime = indmax*timeInterval
890 meteorAux[11:13] = [decayTime, error]
894 meteorAux[11:13] = [decayTime, error]
891
895
892 #Table items 7, 8 and 11
896 #Table items 7, 8 and 11
893 if (riseTime > 0.3): #Number 7: Echo rise exceeds 0.3s
897 if (riseTime > 0.3): #Number 7: Echo rise exceeds 0.3s
894 meteorAux[-1] = 7
898 meteorAux[-1] = 7
895 elif (decayTime < 2*riseTime) : #Number 8: Echo decay time less than than twice rise time
899 elif (decayTime < 2*riseTime) : #Number 8: Echo decay time less than than twice rise time
896 meteorAux[-1] = 8
900 meteorAux[-1] = 8
897 if (error > threshError): #Number 11: Poor fit to amplitude for estimation of decay time
901 if (error > threshError): #Number 11: Poor fit to amplitude for estimation of decay time
898 meteorAux[-1] = 11
902 meteorAux[-1] = 11
899
903
900
904
901 except:
905 except:
902 meteorAux[-1] = 11
906 meteorAux[-1] = 11
903
907
904
908
905 listMeteors1.append(meteorAux)
909 listMeteors1.append(meteorAux)
906
910
907 return listMeteors1
911 return listMeteors1
908
912
909 #Exponential Function
913 #Exponential Function
910
914
911 def __exponential_function(self, x, a, tau):
915 def __exponential_function(self, x, a, tau):
912 y = a*numpy.exp(-x/tau)
916 y = a*numpy.exp(-x/tau)
913 return y
917 return y
914
918
915 def __getRadialVelocity(self, listMeteors, listVolts, radialStdThresh, pairslist, timeInterval):
919 def __getRadialVelocity(self, listMeteors, listVolts, radialStdThresh, pairslist, timeInterval):
916
920
917 pairslist1 = list(pairslist)
921 pairslist1 = list(pairslist)
918 pairslist1.append((0,1))
922 pairslist1.append((0,1))
919 pairslist1.append((3,4))
923 pairslist1.append((3,4))
920 numPairs = len(pairslist1)
924 numPairs = len(pairslist1)
921 #Time Lag
925 #Time Lag
922 timeLag = 45*10**-3
926 timeLag = 45*10**-3
923 c = 3e8
927 c = 3e8
924 lag = numpy.ceil(timeLag/timeInterval)
928 lag = numpy.ceil(timeLag/timeInterval)
925 freq = 30e6
929 freq = 30e6
926
930
927 listMeteors1 = []
931 listMeteors1 = []
928
932
929 for i in range(len(listMeteors)):
933 for i in range(len(listMeteors)):
930 meteor = listMeteors[i]
934 meteor = listMeteors[i]
931 meteorAux = numpy.hstack((meteor[:-1], 0, 0, meteor[-1]))
935 meteorAux = numpy.hstack((meteor[:-1], 0, 0, meteor[-1]))
932 if meteor[-1] == 0:
936 if meteor[-1] == 0:
933 mStart = listMeteors[i][1]
937 mStart = listMeteors[i][1]
934 mPeak = listMeteors[i][2]
938 mPeak = listMeteors[i][2]
935 mLag = mPeak - mStart + lag
939 mLag = mPeak - mStart + lag
936
940
937 #get the volt data between the start and end times of the meteor
941 #get the volt data between the start and end times of the meteor
938 meteorVolts = listVolts[i]
942 meteorVolts = listVolts[i]
939 meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
943 meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
940
944
941 #Get CCF
945 #Get CCF
942 allCCFs = self.__calculateCCF(meteorVolts, pairslist1, [-2,-1,0,1,2])
946 allCCFs = self.__calculateCCF(meteorVolts, pairslist1, [-2,-1,0,1,2])
943
947
944 #Method 2
948 #Method 2
945 slopes = numpy.zeros(numPairs)
949 slopes = numpy.zeros(numPairs)
946 time = numpy.array([-2,-1,1,2])*timeInterval
950 time = numpy.array([-2,-1,1,2])*timeInterval
947 angAllCCF = numpy.angle(allCCFs[:,[0,1,3,4],0])
951 angAllCCF = numpy.angle(allCCFs[:,[0,1,3,4],0])
948
952
949 #Correct phases
953 #Correct phases
950 derPhaseCCF = angAllCCF[:,1:] - angAllCCF[:,0:-1]
954 derPhaseCCF = angAllCCF[:,1:] - angAllCCF[:,0:-1]
951 indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
955 indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
952
956
953 if indDer[0].shape[0] > 0:
957 if indDer[0].shape[0] > 0:
954 for i in range(indDer[0].shape[0]):
958 for i in range(indDer[0].shape[0]):
955 signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i]])
959 signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i]])
956 angAllCCF[indDer[0][i],indDer[1][i]+1:] += signo*2*numpy.pi
960 angAllCCF[indDer[0][i],indDer[1][i]+1:] += signo*2*numpy.pi
957
961
958 # fit = scipy.stats.linregress(numpy.array([-2,-1,1,2])*timeInterval, numpy.array([phaseLagN2s[i],phaseLagN1s[i],phaseLag1s[i],phaseLag2s[i]]))
962 # fit = scipy.stats.linregress(numpy.array([-2,-1,1,2])*timeInterval, numpy.array([phaseLagN2s[i],phaseLagN1s[i],phaseLag1s[i],phaseLag2s[i]]))
959 for j in range(numPairs):
963 for j in range(numPairs):
960 fit = stats.linregress(time, angAllCCF[j,:])
964 fit = stats.linregress(time, angAllCCF[j,:])
961 slopes[j] = fit[0]
965 slopes[j] = fit[0]
962
966
963 #Remove Outlier
967 #Remove Outlier
964 # indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
968 # indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
965 # slopes = numpy.delete(slopes,indOut)
969 # slopes = numpy.delete(slopes,indOut)
966 # indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
970 # indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
967 # slopes = numpy.delete(slopes,indOut)
971 # slopes = numpy.delete(slopes,indOut)
968
972
969 radialVelocity = -numpy.mean(slopes)*(0.25/numpy.pi)*(c/freq)
973 radialVelocity = -numpy.mean(slopes)*(0.25/numpy.pi)*(c/freq)
970 radialError = numpy.std(slopes)*(0.25/numpy.pi)*(c/freq)
974 radialError = numpy.std(slopes)*(0.25/numpy.pi)*(c/freq)
971 meteorAux[-2] = radialError
975 meteorAux[-2] = radialError
972 meteorAux[-3] = radialVelocity
976 meteorAux[-3] = radialVelocity
973
977
974 #Setting Error
978 #Setting Error
975 #Number 15: Radial Drift velocity or projected horizontal velocity exceeds 200 m/s
979 #Number 15: Radial Drift velocity or projected horizontal velocity exceeds 200 m/s
976 if numpy.abs(radialVelocity) > 200:
980 if numpy.abs(radialVelocity) > 200:
977 meteorAux[-1] = 15
981 meteorAux[-1] = 15
978 #Number 12: Poor fit to CCF variation for estimation of radial drift velocity
982 #Number 12: Poor fit to CCF variation for estimation of radial drift velocity
979 elif radialError > radialStdThresh:
983 elif radialError > radialStdThresh:
980 meteorAux[-1] = 12
984 meteorAux[-1] = 12
981
985
982 listMeteors1.append(meteorAux)
986 listMeteors1.append(meteorAux)
983 return listMeteors1
987 return listMeteors1
984
988
985 def __setNewArrays(self, listMeteors, date, heiRang):
989 def __setNewArrays(self, listMeteors, date, heiRang):
986
990
987 #New arrays
991 #New arrays
988 arrayMeteors = numpy.array(listMeteors)
992 arrayMeteors = numpy.array(listMeteors)
989 arrayParameters = numpy.zeros((len(listMeteors),10))
993 arrayParameters = numpy.zeros((len(listMeteors),10))
990
994
991 #Date inclusion
995 #Date inclusion
992 date = re.findall(r'\((.*?)\)', date)
996 date = re.findall(r'\((.*?)\)', date)
993 date = date[0].split(',')
997 date = date[0].split(',')
994 date = map(int, date)
998 date = map(int, date)
995 date = [date[0]*10000 + date[1]*100 + date[2], date[3]*10000 + date[4]*100 + date[5]]
999 date = [date[0]*10000 + date[1]*100 + date[2], date[3]*10000 + date[4]*100 + date[5]]
996 arrayDate = numpy.tile(date, (len(listMeteors), 1))
1000 arrayDate = numpy.tile(date, (len(listMeteors), 1))
997
1001
998 #Meteor array
1002 #Meteor array
999 arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)]
1003 arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)]
1000 arrayMeteors = numpy.hstack((arrayDate, arrayMeteors))
1004 arrayMeteors = numpy.hstack((arrayDate, arrayMeteors))
1001
1005
1002 #Parameters Array
1006 #Parameters Array
1003 arrayParameters[:,0:3] = arrayMeteors[:,0:3]
1007 arrayParameters[:,0:3] = arrayMeteors[:,0:3]
1004 arrayParameters[:,-3:] = arrayMeteors[:,-3:]
1008 arrayParameters[:,-3:] = arrayMeteors[:,-3:]
1005
1009
1006 return arrayMeteors, arrayParameters
1010 return arrayMeteors, arrayParameters
1007
1011
1008 def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth):
1012 def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth):
1009
1013
1010 arrayAOA = numpy.zeros((phases.shape[0],3))
1014 arrayAOA = numpy.zeros((phases.shape[0],3))
1011 cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList)
1015 cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList)
1012
1016
1013 arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
1017 arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
1014 cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
1018 cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
1015 arrayAOA[:,2] = cosDirError
1019 arrayAOA[:,2] = cosDirError
1016
1020
1017 azimuthAngle = arrayAOA[:,0]
1021 azimuthAngle = arrayAOA[:,0]
1018 zenithAngle = arrayAOA[:,1]
1022 zenithAngle = arrayAOA[:,1]
1019
1023
1020 #Setting Error
1024 #Setting Error
1021 #Number 3: AOA not fesible
1025 #Number 3: AOA not fesible
1022 indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
1026 indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
1023 error[indInvalid] = 3
1027 error[indInvalid] = 3
1024 #Number 4: Large difference in AOAs obtained from different antenna baselines
1028 #Number 4: Large difference in AOAs obtained from different antenna baselines
1025 indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
1029 indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
1026 error[indInvalid] = 4
1030 error[indInvalid] = 4
1027 return arrayAOA, error
1031 return arrayAOA, error
1028
1032
1029 def __getDirectionCosines(self, arrayPhase, pairsList):
1033 def __getDirectionCosines(self, arrayPhase, pairsList):
1030
1034
1031 #Initializing some variables
1035 #Initializing some variables
1032 ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
1036 ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
1033 ang_aux = ang_aux.reshape(1,ang_aux.size)
1037 ang_aux = ang_aux.reshape(1,ang_aux.size)
1034
1038
1035 cosdir = numpy.zeros((arrayPhase.shape[0],2))
1039 cosdir = numpy.zeros((arrayPhase.shape[0],2))
1036 cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
1040 cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
1037
1041
1038
1042
1039 for i in range(2):
1043 for i in range(2):
1040 #First Estimation
1044 #First Estimation
1041 phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]]
1045 phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]]
1042 #Dealias
1046 #Dealias
1043 indcsi = numpy.where(phi0_aux > numpy.pi)
1047 indcsi = numpy.where(phi0_aux > numpy.pi)
1044 phi0_aux[indcsi] -= 2*numpy.pi
1048 phi0_aux[indcsi] -= 2*numpy.pi
1045 indcsi = numpy.where(phi0_aux < -numpy.pi)
1049 indcsi = numpy.where(phi0_aux < -numpy.pi)
1046 phi0_aux[indcsi] += 2*numpy.pi
1050 phi0_aux[indcsi] += 2*numpy.pi
1047 #Direction Cosine 0
1051 #Direction Cosine 0
1048 cosdir0[:,i] = -(phi0_aux)/(2*numpy.pi*0.5)
1052 cosdir0[:,i] = -(phi0_aux)/(2*numpy.pi*0.5)
1049
1053
1050 #Most-Accurate Second Estimation
1054 #Most-Accurate Second Estimation
1051 phi1_aux = arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]]
1055 phi1_aux = arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]]
1052 phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
1056 phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
1053 #Direction Cosine 1
1057 #Direction Cosine 1
1054 cosdir1 = -(phi1_aux + ang_aux)/(2*numpy.pi*4.5)
1058 cosdir1 = -(phi1_aux + ang_aux)/(2*numpy.pi*4.5)
1055
1059
1056 #Searching the correct Direction Cosine
1060 #Searching the correct Direction Cosine
1057 cosdir0_aux = cosdir0[:,i]
1061 cosdir0_aux = cosdir0[:,i]
1058 cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
1062 cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
1059 #Minimum Distance
1063 #Minimum Distance
1060 cosDiff = (cosdir1 - cosdir0_aux)**2
1064 cosDiff = (cosdir1 - cosdir0_aux)**2
1061 indcos = cosDiff.argmin(axis = 1)
1065 indcos = cosDiff.argmin(axis = 1)
1062 #Saving Value obtained
1066 #Saving Value obtained
1063 cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
1067 cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
1064
1068
1065 return cosdir0, cosdir
1069 return cosdir0, cosdir
1066
1070
1067 def __calculateAOA(self, cosdir, azimuth):
1071 def __calculateAOA(self, cosdir, azimuth):
1068 cosdirX = cosdir[:,0]
1072 cosdirX = cosdir[:,0]
1069 cosdirY = cosdir[:,1]
1073 cosdirY = cosdir[:,1]
1070
1074
1071 zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
1075 zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
1072 azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth #0 deg north, 90 deg east
1076 azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth #0 deg north, 90 deg east
1073 angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
1077 angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
1074
1078
1075 return angles
1079 return angles
1076
1080
1077 def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
1081 def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
1078
1082
1079 Ramb = 375 #Ramb = c/(2*PRF)
1083 Ramb = 375 #Ramb = c/(2*PRF)
1080 Re = 6371 #Earth Radius
1084 Re = 6371 #Earth Radius
1081 heights = numpy.zeros(Ranges.shape)
1085 heights = numpy.zeros(Ranges.shape)
1082
1086
1083 R_aux = numpy.array([0,1,2])*Ramb
1087 R_aux = numpy.array([0,1,2])*Ramb
1084 R_aux = R_aux.reshape(1,R_aux.size)
1088 R_aux = R_aux.reshape(1,R_aux.size)
1085
1089
1086 Ranges = Ranges.reshape(Ranges.size,1)
1090 Ranges = Ranges.reshape(Ranges.size,1)
1087
1091
1088 Ri = Ranges + R_aux
1092 Ri = Ranges + R_aux
1089 hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re
1093 hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re
1090
1094
1091 #Check if there is a height between 70 and 110 km
1095 #Check if there is a height between 70 and 110 km
1092 h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
1096 h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
1093 ind_h = numpy.where(h_bool == 1)[0]
1097 ind_h = numpy.where(h_bool == 1)[0]
1094
1098
1095 hCorr = hi[ind_h, :]
1099 hCorr = hi[ind_h, :]
1096 ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
1100 ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
1097
1101
1098 hCorr = hi[ind_hCorr]
1102 hCorr = hi[ind_hCorr]
1099 heights[ind_h] = hCorr
1103 heights[ind_h] = hCorr
1100
1104
1101 #Setting Error
1105 #Setting Error
1102 #Number 13: Height unresolvable echo: not valid height within 70 to 110 km
1106 #Number 13: Height unresolvable echo: not valid height within 70 to 110 km
1103 #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
1107 #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
1104
1108
1105 indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
1109 indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
1106 error[indInvalid2] = 14
1110 error[indInvalid2] = 14
1107 indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
1111 indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
1108 error[indInvalid1] = 13
1112 error[indInvalid1] = 13
1109
1113
1110 return heights, error
1114 return heights, error
1111
1115
1116 def SpectralFitting(self, getSNR = True, path=None, file=None, groupList=None):
1117
1118 '''
1119 Function GetMoments()
1120
1121 Input:
1122 Output:
1123 Variables modified:
1124 '''
1125 if path != None:
1126 sys.path.append(path)
1127 self.dataOut.library = importlib.import_module(file)
1128
1129 #To be inserted as a parameter
1130 groupArray = numpy.array(groupList)
1131 # groupArray = numpy.array([[0,1],[2,3]])
1132 self.dataOut.groupList = groupArray
1133
1134 nGroups = groupArray.shape[0]
1135 nChannels = self.dataIn.nChannels
1136 nHeights=self.dataIn.heightList.size
1137
1138 #Parameters Array
1139 self.dataOut.data_param = None
1140
1141 #Set constants
1142 constants = self.dataOut.library.setConstants(self.dataIn)
1143 self.dataOut.constants = constants
1144 M = self.dataIn.normFactor
1145 N = self.dataIn.nFFTPoints
1146 ippSeconds = self.dataIn.ippSeconds
1147 K = self.dataIn.nIncohInt
1148 pairsArray = numpy.array(self.dataIn.pairsList)
1149
1150 #List of possible combinations
1151 listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)
1152 indCross = numpy.zeros(len(list(listComb)), dtype = 'int')
1153
1154 if getSNR:
1155 listChannels = groupArray.reshape((groupArray.size))
1156 listChannels.sort()
1157 noise = self.dataIn.getNoise()
1158 self.dataOut.SNR = self.__getSNR(self.dataIn.data_spc[listChannels,:,:], noise[listChannels])
1159
1160 for i in range(nGroups):
1161 coord = groupArray[i,:]
1162
1163 #Input data array
1164 data = self.dataIn.data_spc[coord,:,:]/(M*N)
1165 data = data.reshape((data.shape[0]*data.shape[1],data.shape[2]))
1166
1167 #Cross Spectra data array for Covariance Matrixes
1168 ind = 0
1169 for pairs in listComb:
1170 pairsSel = numpy.array([coord[x],coord[y]])
1171 indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0])
1172 ind += 1
1173 dataCross = self.dataIn.data_cspc[indCross,:,:]/(M*N)
1174 dataCross = dataCross**2/K
1175
1176 for h in range(nHeights):
1177 # print self.dataOut.heightList[h]
1178
1179 #Input
1180 d = data[:,h]
1181
1182 #Covariance Matrix
1183 D = numpy.diag(d**2/K)
1184 ind = 0
1185 for pairs in listComb:
1186 #Coordinates in Covariance Matrix
1187 x = pairs[0]
1188 y = pairs[1]
1189 #Channel Index
1190 S12 = dataCross[ind,:,h]
1191 D12 = numpy.diag(S12)
1192 #Completing Covariance Matrix with Cross Spectras
1193 D[x*N:(x+1)*N,y*N:(y+1)*N] = D12
1194 D[y*N:(y+1)*N,x*N:(x+1)*N] = D12
1195 ind += 1
1196 Dinv=numpy.linalg.inv(D)
1197 L=numpy.linalg.cholesky(Dinv)
1198 LT=L.T
1199
1200 dp = numpy.dot(LT,d)
1201
1202 #Initial values
1203 data_spc = self.dataIn.data_spc[coord,:,h]
1204 p0 = self.dataOut.library.initialValuesFunction(data_spc, constants)
1205
1206 #Least Squares
1207 minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True)
1208 # minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))
1209 #Chi square error
1210 error0 = numpy.sum(infodict['fvec']**2)/(2*N)
1211 # error0 = 0
1212 #Error with Jacobian
1213 error1 = self.dataOut.library.errorFunction(minp,constants,LT)
1214 #Save
1215 if self.dataOut.data_param == None:
1216 self.dataOut.data_param = numpy.zeros((nGroups, minp.size, nHeights))*numpy.nan
1217 self.dataOut.error = numpy.zeros((nGroups, error1.size + 1, nHeights))*numpy.nan
1218
1219 self.dataOut.error[i,:,h] = numpy.hstack((error0,error1))
1220 self.dataOut.data_param[i,:,h] = minp
1221 return
1222
1223
1224 def __residFunction(self, p, dp, LT, constants):
1225
1226 fm = self.dataOut.library.modelFunction(p, constants)
1227 fmp=numpy.dot(LT,fm)
1228
1229 return dp-fmp
1230
1231 def __getSNR(self, z, noise):
1232
1233 avg = numpy.average(z, axis=1)
1234 SNR = (avg.T-noise)/noise
1235 SNR = SNR.T
1236 return SNR
1237
1238 def __chisq(p,chindex,hindex):
1239 #similar to Resid but calculates CHI**2
1240 [LT,d,fm]=setupLTdfm(p,chindex,hindex)
1241 dp=numpy.dot(LT,d)
1242 fmp=numpy.dot(LT,fm)
1243 chisq=numpy.dot((dp-fmp).T,(dp-fmp))
1244 return chisq
1245
1246
1112
1247
1113 class WindProfiler(Operation):
1248 class WindProfiler(Operation):
1114
1249
1115 __isConfig = False
1250 __isConfig = False
1116
1251
1117 __initime = None
1252 __initime = None
1118 __lastdatatime = None
1253 __lastdatatime = None
1119 __integrationtime = None
1254 __integrationtime = None
1120
1255
1121 __buffer = None
1256 __buffer = None
1122
1257
1123 __dataReady = False
1258 __dataReady = False
1124
1259
1125 __firstdata = None
1260 __firstdata = None
1126
1261
1127 n = None
1262 n = None
1128
1263
1129 def __init__(self):
1264 def __init__(self):
1130 Operation.__init__(self)
1265 Operation.__init__(self)
1131
1266
1132 def __calculateCosDir(self, elev, azim):
1267 def __calculateCosDir(self, elev, azim):
1133 zen = (90 - elev)*numpy.pi/180
1268 zen = (90 - elev)*numpy.pi/180
1134 azim = azim*numpy.pi/180
1269 azim = azim*numpy.pi/180
1135 cosDirX = numpy.sqrt((1-numpy.cos(zen)**2)/((1+numpy.tan(azim)**2)))
1270 cosDirX = numpy.sqrt((1-numpy.cos(zen)**2)/((1+numpy.tan(azim)**2)))
1136 cosDirY = numpy.sqrt(1-numpy.cos(zen)**2-cosDirX**2)
1271 cosDirY = numpy.sqrt(1-numpy.cos(zen)**2-cosDirX**2)
1137
1272
1138 signX = numpy.sign(numpy.cos(azim))
1273 signX = numpy.sign(numpy.cos(azim))
1139 signY = numpy.sign(numpy.sin(azim))
1274 signY = numpy.sign(numpy.sin(azim))
1140
1275
1141 cosDirX = numpy.copysign(cosDirX, signX)
1276 cosDirX = numpy.copysign(cosDirX, signX)
1142 cosDirY = numpy.copysign(cosDirY, signY)
1277 cosDirY = numpy.copysign(cosDirY, signY)
1143 return cosDirX, cosDirY
1278 return cosDirX, cosDirY
1144
1279
1145 def __calculateAngles(self, theta_x, theta_y, azimuth):
1280 def __calculateAngles(self, theta_x, theta_y, azimuth):
1146
1281
1147 dir_cosw = numpy.sqrt(1-theta_x**2-theta_y**2)
1282 dir_cosw = numpy.sqrt(1-theta_x**2-theta_y**2)
1148 zenith_arr = numpy.arccos(dir_cosw)
1283 zenith_arr = numpy.arccos(dir_cosw)
1149 azimuth_arr = numpy.arctan2(theta_x,theta_y) + azimuth*math.pi/180
1284 azimuth_arr = numpy.arctan2(theta_x,theta_y) + azimuth*math.pi/180
1150
1285
1151 dir_cosu = numpy.sin(azimuth_arr)*numpy.sin(zenith_arr)
1286 dir_cosu = numpy.sin(azimuth_arr)*numpy.sin(zenith_arr)
1152 dir_cosv = numpy.cos(azimuth_arr)*numpy.sin(zenith_arr)
1287 dir_cosv = numpy.cos(azimuth_arr)*numpy.sin(zenith_arr)
1153
1288
1154 return azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw
1289 return azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw
1155
1290
1156 def __calculateMatA(self, dir_cosu, dir_cosv, dir_cosw, horOnly):
1291 def __calculateMatA(self, dir_cosu, dir_cosv, dir_cosw, horOnly):
1157
1292
1158 #
1293 #
1159 if horOnly:
1294 if horOnly:
1160 A = numpy.c_[dir_cosu,dir_cosv]
1295 A = numpy.c_[dir_cosu,dir_cosv]
1161 else:
1296 else:
1162 A = numpy.c_[dir_cosu,dir_cosv,dir_cosw]
1297 A = numpy.c_[dir_cosu,dir_cosv,dir_cosw]
1163 A = numpy.asmatrix(A)
1298 A = numpy.asmatrix(A)
1164 A1 = numpy.linalg.inv(A.transpose()*A)*A.transpose()
1299 A1 = numpy.linalg.inv(A.transpose()*A)*A.transpose()
1165
1300
1166 return A1
1301 return A1
1167
1302
1168 def __correctValues(self, heiRang, phi, velRadial, SNR):
1303 def __correctValues(self, heiRang, phi, velRadial, SNR):
1169 listPhi = phi.tolist()
1304 listPhi = phi.tolist()
1170 maxid = listPhi.index(max(listPhi))
1305 maxid = listPhi.index(max(listPhi))
1171 minid = listPhi.index(min(listPhi))
1306 minid = listPhi.index(min(listPhi))
1172
1307
1173 rango = range(len(phi))
1308 rango = range(len(phi))
1174 # rango = numpy.delete(rango,maxid)
1309 # rango = numpy.delete(rango,maxid)
1175
1310
1176 heiRang1 = heiRang*math.cos(phi[maxid])
1311 heiRang1 = heiRang*math.cos(phi[maxid])
1177 heiRangAux = heiRang*math.cos(phi[minid])
1312 heiRangAux = heiRang*math.cos(phi[minid])
1178 indOut = (heiRang1 < heiRangAux[0]).nonzero()
1313 indOut = (heiRang1 < heiRangAux[0]).nonzero()
1179 heiRang1 = numpy.delete(heiRang1,indOut)
1314 heiRang1 = numpy.delete(heiRang1,indOut)
1180
1315
1181 velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
1316 velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
1182 SNR1 = numpy.zeros([len(phi),len(heiRang1)])
1317 SNR1 = numpy.zeros([len(phi),len(heiRang1)])
1183
1318
1184 for i in rango:
1319 for i in rango:
1185 x = heiRang*math.cos(phi[i])
1320 x = heiRang*math.cos(phi[i])
1186 y1 = velRadial[i,:]
1321 y1 = velRadial[i,:]
1187 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
1322 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
1188
1323
1189 x1 = heiRang1
1324 x1 = heiRang1
1190 y11 = f1(x1)
1325 y11 = f1(x1)
1191
1326
1192 y2 = SNR[i,:]
1327 y2 = SNR[i,:]
1193 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
1328 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
1194 y21 = f2(x1)
1329 y21 = f2(x1)
1195
1330
1196 velRadial1[i,:] = y11
1331 velRadial1[i,:] = y11
1197 SNR1[i,:] = y21
1332 SNR1[i,:] = y21
1198
1333
1199 return heiRang1, velRadial1, SNR1
1334 return heiRang1, velRadial1, SNR1
1200
1335
1201 def __calculateVelUVW(self, A, velRadial):
1336 def __calculateVelUVW(self, A, velRadial):
1202
1337
1203 #Operacion Matricial
1338 #Operacion Matricial
1204 # velUVW = numpy.zeros((velRadial.shape[1],3))
1339 # velUVW = numpy.zeros((velRadial.shape[1],3))
1205 # for ind in range(velRadial.shape[1]):
1340 # for ind in range(velRadial.shape[1]):
1206 # velUVW[ind,:] = numpy.dot(A,velRadial[:,ind])
1341 # velUVW[ind,:] = numpy.dot(A,velRadial[:,ind])
1207 # velUVW = velUVW.transpose()
1342 # velUVW = velUVW.transpose()
1208 velUVW = numpy.zeros((A.shape[0],velRadial.shape[1]))
1343 velUVW = numpy.zeros((A.shape[0],velRadial.shape[1]))
1209 velUVW[:,:] = numpy.dot(A,velRadial)
1344 velUVW[:,:] = numpy.dot(A,velRadial)
1210
1345
1211
1346
1212 return velUVW
1347 return velUVW
1213
1348
1214 def techniqueDBS(self, velRadial0, dirCosx, disrCosy, azimuth, correct, horizontalOnly, heiRang, SNR0):
1349 def techniqueDBS(self, velRadial0, dirCosx, disrCosy, azimuth, correct, horizontalOnly, heiRang, SNR0):
1215 """
1350 """
1216 Function that implements Doppler Beam Swinging (DBS) technique.
1351 Function that implements Doppler Beam Swinging (DBS) technique.
1217
1352
1218 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
1353 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
1219 Direction correction (if necessary), Ranges and SNR
1354 Direction correction (if necessary), Ranges and SNR
1220
1355
1221 Output: Winds estimation (Zonal, Meridional and Vertical)
1356 Output: Winds estimation (Zonal, Meridional and Vertical)
1222
1357
1223 Parameters affected: Winds, height range, SNR
1358 Parameters affected: Winds, height range, SNR
1224 """
1359 """
1225 azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(dirCosx, disrCosy, azimuth)
1360 azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(dirCosx, disrCosy, azimuth)
1226 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correct*velRadial0, SNR0)
1361 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correct*velRadial0, SNR0)
1227 A = self.__calculateMatA(dir_cosu, dir_cosv, dir_cosw, horizontalOnly)
1362 A = self.__calculateMatA(dir_cosu, dir_cosv, dir_cosw, horizontalOnly)
1228
1363
1229 #Calculo de Componentes de la velocidad con DBS
1364 #Calculo de Componentes de la velocidad con DBS
1230 winds = self.__calculateVelUVW(A,velRadial1)
1365 winds = self.__calculateVelUVW(A,velRadial1)
1231
1366
1232 return winds, heiRang1, SNR1
1367 return winds, heiRang1, SNR1
1233
1368
1234 def __calculateDistance(self, posx, posy, pairsCrossCorr, pairsList, pairs, azimuth = None):
1369 def __calculateDistance(self, posx, posy, pairsCrossCorr, pairsList, pairs, azimuth = None):
1235
1370
1236 posx = numpy.asarray(posx)
1371 posx = numpy.asarray(posx)
1237 posy = numpy.asarray(posy)
1372 posy = numpy.asarray(posy)
1238
1373
1239 #Rotacion Inversa para alinear con el azimuth
1374 #Rotacion Inversa para alinear con el azimuth
1240 if azimuth!= None:
1375 if azimuth!= None:
1241 azimuth = azimuth*math.pi/180
1376 azimuth = azimuth*math.pi/180
1242 posx1 = posx*math.cos(azimuth) + posy*math.sin(azimuth)
1377 posx1 = posx*math.cos(azimuth) + posy*math.sin(azimuth)
1243 posy1 = -posx*math.sin(azimuth) + posy*math.cos(azimuth)
1378 posy1 = -posx*math.sin(azimuth) + posy*math.cos(azimuth)
1244 else:
1379 else:
1245 posx1 = posx
1380 posx1 = posx
1246 posy1 = posy
1381 posy1 = posy
1247
1382
1248 #Calculo de Distancias
1383 #Calculo de Distancias
1249 distx = numpy.zeros(pairsCrossCorr.size)
1384 distx = numpy.zeros(pairsCrossCorr.size)
1250 disty = numpy.zeros(pairsCrossCorr.size)
1385 disty = numpy.zeros(pairsCrossCorr.size)
1251 dist = numpy.zeros(pairsCrossCorr.size)
1386 dist = numpy.zeros(pairsCrossCorr.size)
1252 ang = numpy.zeros(pairsCrossCorr.size)
1387 ang = numpy.zeros(pairsCrossCorr.size)
1253
1388
1254 for i in range(pairsCrossCorr.size):
1389 for i in range(pairsCrossCorr.size):
1255 distx[i] = posx1[pairsList[pairsCrossCorr[i]][1]] - posx1[pairsList[pairsCrossCorr[i]][0]]
1390 distx[i] = posx1[pairsList[pairsCrossCorr[i]][1]] - posx1[pairsList[pairsCrossCorr[i]][0]]
1256 disty[i] = posy1[pairsList[pairsCrossCorr[i]][1]] - posy1[pairsList[pairsCrossCorr[i]][0]]
1391 disty[i] = posy1[pairsList[pairsCrossCorr[i]][1]] - posy1[pairsList[pairsCrossCorr[i]][0]]
1257 dist[i] = numpy.sqrt(distx[i]**2 + disty[i]**2)
1392 dist[i] = numpy.sqrt(distx[i]**2 + disty[i]**2)
1258 ang[i] = numpy.arctan2(disty[i],distx[i])
1393 ang[i] = numpy.arctan2(disty[i],distx[i])
1259 #Calculo de Matrices
1394 #Calculo de Matrices
1260 nPairs = len(pairs)
1395 nPairs = len(pairs)
1261 ang1 = numpy.zeros((nPairs, 2, 1))
1396 ang1 = numpy.zeros((nPairs, 2, 1))
1262 dist1 = numpy.zeros((nPairs, 2, 1))
1397 dist1 = numpy.zeros((nPairs, 2, 1))
1263
1398
1264 for j in range(nPairs):
1399 for j in range(nPairs):
1265 dist1[j,0,0] = dist[pairs[j][0]]
1400 dist1[j,0,0] = dist[pairs[j][0]]
1266 dist1[j,1,0] = dist[pairs[j][1]]
1401 dist1[j,1,0] = dist[pairs[j][1]]
1267 ang1[j,0,0] = ang[pairs[j][0]]
1402 ang1[j,0,0] = ang[pairs[j][0]]
1268 ang1[j,1,0] = ang[pairs[j][1]]
1403 ang1[j,1,0] = ang[pairs[j][1]]
1269
1404
1270 return distx,disty, dist1,ang1
1405 return distx,disty, dist1,ang1
1271
1406
1272 def __calculateVelVer(self, phase, lagTRange, _lambda):
1407 def __calculateVelVer(self, phase, lagTRange, _lambda):
1273
1408
1274 Ts = lagTRange[1] - lagTRange[0]
1409 Ts = lagTRange[1] - lagTRange[0]
1275 velW = -_lambda*phase/(4*math.pi*Ts)
1410 velW = -_lambda*phase/(4*math.pi*Ts)
1276
1411
1277 return velW
1412 return velW
1278
1413
1279 def __calculateVelHorDir(self, dist, tau1, tau2, ang):
1414 def __calculateVelHorDir(self, dist, tau1, tau2, ang):
1280 nPairs = tau1.shape[0]
1415 nPairs = tau1.shape[0]
1281 vel = numpy.zeros((nPairs,3,tau1.shape[2]))
1416 vel = numpy.zeros((nPairs,3,tau1.shape[2]))
1282
1417
1283 angCos = numpy.cos(ang)
1418 angCos = numpy.cos(ang)
1284 angSin = numpy.sin(ang)
1419 angSin = numpy.sin(ang)
1285
1420
1286 vel0 = dist*tau1/(2*tau2**2)
1421 vel0 = dist*tau1/(2*tau2**2)
1287 vel[:,0,:] = (vel0*angCos).sum(axis = 1)
1422 vel[:,0,:] = (vel0*angCos).sum(axis = 1)
1288 vel[:,1,:] = (vel0*angSin).sum(axis = 1)
1423 vel[:,1,:] = (vel0*angSin).sum(axis = 1)
1289
1424
1290 ind = numpy.where(numpy.isinf(vel))
1425 ind = numpy.where(numpy.isinf(vel))
1291 vel[ind] = numpy.nan
1426 vel[ind] = numpy.nan
1292
1427
1293 return vel
1428 return vel
1294
1429
1295 def __getPairsAutoCorr(self, pairsList, nChannels):
1430 def __getPairsAutoCorr(self, pairsList, nChannels):
1296
1431
1297 pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
1432 pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
1298
1433
1299 for l in range(len(pairsList)):
1434 for l in range(len(pairsList)):
1300 firstChannel = pairsList[l][0]
1435 firstChannel = pairsList[l][0]
1301 secondChannel = pairsList[l][1]
1436 secondChannel = pairsList[l][1]
1302
1437
1303 #Obteniendo pares de Autocorrelacion
1438 #Obteniendo pares de Autocorrelacion
1304 if firstChannel == secondChannel:
1439 if firstChannel == secondChannel:
1305 pairsAutoCorr[firstChannel] = int(l)
1440 pairsAutoCorr[firstChannel] = int(l)
1306
1441
1307 pairsAutoCorr = pairsAutoCorr.astype(int)
1442 pairsAutoCorr = pairsAutoCorr.astype(int)
1308
1443
1309 pairsCrossCorr = range(len(pairsList))
1444 pairsCrossCorr = range(len(pairsList))
1310 pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
1445 pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
1311
1446
1312 return pairsAutoCorr, pairsCrossCorr
1447 return pairsAutoCorr, pairsCrossCorr
1313
1448
1314 def techniqueSA(self, pairsSelected, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, lagTRange, correctFactor):
1449 def techniqueSA(self, pairsSelected, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, lagTRange, correctFactor):
1315 """
1450 """
1316 Function that implements Spaced Antenna (SA) technique.
1451 Function that implements Spaced Antenna (SA) technique.
1317
1452
1318 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
1453 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
1319 Direction correction (if necessary), Ranges and SNR
1454 Direction correction (if necessary), Ranges and SNR
1320
1455
1321 Output: Winds estimation (Zonal, Meridional and Vertical)
1456 Output: Winds estimation (Zonal, Meridional and Vertical)
1322
1457
1323 Parameters affected: Winds
1458 Parameters affected: Winds
1324 """
1459 """
1325 #Cross Correlation pairs obtained
1460 #Cross Correlation pairs obtained
1326 pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)
1461 pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)
1327 pairsArray = numpy.array(pairsList)[pairsCrossCorr]
1462 pairsArray = numpy.array(pairsList)[pairsCrossCorr]
1328 pairsSelArray = numpy.array(pairsSelected)
1463 pairsSelArray = numpy.array(pairsSelected)
1329 pairs = []
1464 pairs = []
1330
1465
1331 #Wind estimation pairs obtained
1466 #Wind estimation pairs obtained
1332 for i in range(pairsSelArray.shape[0]/2):
1467 for i in range(pairsSelArray.shape[0]/2):
1333 ind1 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i], axis = 1))[0][0]
1468 ind1 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i], axis = 1))[0][0]
1334 ind2 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i + 1], axis = 1))[0][0]
1469 ind2 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i + 1], axis = 1))[0][0]
1335 pairs.append((ind1,ind2))
1470 pairs.append((ind1,ind2))
1336
1471
1337 indtau = tau.shape[0]/2
1472 indtau = tau.shape[0]/2
1338 tau1 = tau[:indtau,:]
1473 tau1 = tau[:indtau,:]
1339 tau2 = tau[indtau:-1,:]
1474 tau2 = tau[indtau:-1,:]
1340 tau1 = tau1[pairs,:]
1475 tau1 = tau1[pairs,:]
1341 tau2 = tau2[pairs,:]
1476 tau2 = tau2[pairs,:]
1342 phase1 = tau[-1,:]
1477 phase1 = tau[-1,:]
1343
1478
1344 #---------------------------------------------------------------------
1479 #---------------------------------------------------------------------
1345 #Metodo Directo
1480 #Metodo Directo
1346 distx, disty, dist, ang = self.__calculateDistance(position_x, position_y, pairsCrossCorr, pairsList, pairs,azimuth)
1481 distx, disty, dist, ang = self.__calculateDistance(position_x, position_y, pairsCrossCorr, pairsList, pairs,azimuth)
1347 winds = self.__calculateVelHorDir(dist, tau1, tau2, ang)
1482 winds = self.__calculateVelHorDir(dist, tau1, tau2, ang)
1348 winds = stats.nanmean(winds, axis=0)
1483 winds = stats.nanmean(winds, axis=0)
1349 #---------------------------------------------------------------------
1484 #---------------------------------------------------------------------
1350 #Metodo General
1485 #Metodo General
1351 # distx, disty, dist = self.calculateDistance(position_x,position_y,pairsCrossCorr, pairsList, azimuth)
1486 # distx, disty, dist = self.calculateDistance(position_x,position_y,pairsCrossCorr, pairsList, azimuth)
1352 # #Calculo Coeficientes de Funcion de Correlacion
1487 # #Calculo Coeficientes de Funcion de Correlacion
1353 # F,G,A,B,H = self.calculateCoef(tau1,tau2,distx,disty,n)
1488 # F,G,A,B,H = self.calculateCoef(tau1,tau2,distx,disty,n)
1354 # #Calculo de Velocidades
1489 # #Calculo de Velocidades
1355 # winds = self.calculateVelUV(F,G,A,B,H)
1490 # winds = self.calculateVelUV(F,G,A,B,H)
1356
1491
1357 #---------------------------------------------------------------------
1492 #---------------------------------------------------------------------
1358 winds[2,:] = self.__calculateVelVer(phase1, lagTRange, _lambda)
1493 winds[2,:] = self.__calculateVelVer(phase1, lagTRange, _lambda)
1359 winds = correctFactor*winds
1494 winds = correctFactor*winds
1360 return winds
1495 return winds
1361
1496
1362 def __checkTime(self, currentTime, paramInterval, windsInterval):
1497 def __checkTime(self, currentTime, paramInterval, outputInterval):
1363
1498
1364 dataTime = currentTime + paramInterval
1499 dataTime = currentTime + paramInterval
1365 deltaTime = dataTime - self.__initime
1500 deltaTime = dataTime - self.__initime
1366
1501
1367 if deltaTime >= windsInterval or deltaTime < 0:
1502 if deltaTime >= outputInterval or deltaTime < 0:
1368 self.__dataReady = True
1503 self.__dataReady = True
1369 return
1504 return
1370
1505
1371 def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax):
1506 def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax):
1372 '''
1507 '''
1373 Function that implements winds estimation technique with detected meteors.
1508 Function that implements winds estimation technique with detected meteors.
1374
1509
1375 Input: Detected meteors, Minimum meteor quantity to wind estimation
1510 Input: Detected meteors, Minimum meteor quantity to wind estimation
1376
1511
1377 Output: Winds estimation (Zonal and Meridional)
1512 Output: Winds estimation (Zonal and Meridional)
1378
1513
1379 Parameters affected: Winds
1514 Parameters affected: Winds
1380 '''
1515 '''
1381 #Settings
1516 #Settings
1382 nInt = (heightMax - heightMin)/2
1517 nInt = (heightMax - heightMin)/2
1383 winds = numpy.zeros((2,nInt))*numpy.nan
1518 winds = numpy.zeros((2,nInt))*numpy.nan
1384
1519
1385 #Filter errors
1520 #Filter errors
1386 error = numpy.where(arrayMeteor[:,-1] == 0)[0]
1521 error = numpy.where(arrayMeteor[:,-1] == 0)[0]
1387 finalMeteor = arrayMeteor[error,:]
1522 finalMeteor = arrayMeteor[error,:]
1388
1523
1389 #Meteor Histogram
1524 #Meteor Histogram
1390 finalHeights = finalMeteor[:,3]
1525 finalHeights = finalMeteor[:,3]
1391 hist = numpy.histogram(finalHeights, bins = nInt, range = (heightMin,heightMax))
1526 hist = numpy.histogram(finalHeights, bins = nInt, range = (heightMin,heightMax))
1392 nMeteorsPerI = hist[0]
1527 nMeteorsPerI = hist[0]
1393 heightPerI = hist[1]
1528 heightPerI = hist[1]
1394
1529
1395 #Sort of meteors
1530 #Sort of meteors
1396 indSort = finalHeights.argsort()
1531 indSort = finalHeights.argsort()
1397 finalMeteor2 = finalMeteor[indSort,:]
1532 finalMeteor2 = finalMeteor[indSort,:]
1398
1533
1399 # Calculating winds
1534 # Calculating winds
1400 ind1 = 0
1535 ind1 = 0
1401 ind2 = 0
1536 ind2 = 0
1402
1537
1403 for i in range(nInt):
1538 for i in range(nInt):
1404 nMet = nMeteorsPerI[i]
1539 nMet = nMeteorsPerI[i]
1405 ind1 = ind2
1540 ind1 = ind2
1406 ind2 = ind1 + nMet
1541 ind2 = ind1 + nMet
1407
1542
1408 meteorAux = finalMeteor2[ind1:ind2,:]
1543 meteorAux = finalMeteor2[ind1:ind2,:]
1409
1544
1410 if meteorAux.shape[0] >= meteorThresh:
1545 if meteorAux.shape[0] >= meteorThresh:
1411 vel = meteorAux[:, 7]
1546 vel = meteorAux[:, 7]
1412 zen = meteorAux[:, 5]*numpy.pi/180
1547 zen = meteorAux[:, 5]*numpy.pi/180
1413 azim = meteorAux[:, 4]*numpy.pi/180
1548 azim = meteorAux[:, 4]*numpy.pi/180
1414
1549
1415 n = numpy.cos(zen)
1550 n = numpy.cos(zen)
1416 # m = (1 - n**2)/(1 - numpy.tan(azim)**2)
1551 # m = (1 - n**2)/(1 - numpy.tan(azim)**2)
1417 # l = m*numpy.tan(azim)
1552 # l = m*numpy.tan(azim)
1418 l = numpy.sin(zen)*numpy.sin(azim)
1553 l = numpy.sin(zen)*numpy.sin(azim)
1419 m = numpy.sin(zen)*numpy.cos(azim)
1554 m = numpy.sin(zen)*numpy.cos(azim)
1420
1555
1421 A = numpy.vstack((l, m)).transpose()
1556 A = numpy.vstack((l, m)).transpose()
1422 A1 = numpy.dot(numpy.linalg.inv( numpy.dot(A.transpose(),A) ),A.transpose())
1557 A1 = numpy.dot(numpy.linalg.inv( numpy.dot(A.transpose(),A) ),A.transpose())
1423 windsAux = numpy.dot(A1, vel)
1558 windsAux = numpy.dot(A1, vel)
1424
1559
1425 winds[0,i] = windsAux[0]
1560 winds[0,i] = windsAux[0]
1426 winds[1,i] = windsAux[1]
1561 winds[1,i] = windsAux[1]
1427
1562
1428 return winds, heightPerI[:-1]
1563 return winds, heightPerI[:-1]
1429
1564
1430 def run(self, dataOut, technique, **kwargs):
1565 def run(self, dataOut, technique, **kwargs):
1431
1566
1432 param = dataOut.data_param
1567 param = dataOut.data_param
1433 if dataOut.abscissaRange != None:
1568 if dataOut.abscissaRange != None:
1434 absc = dataOut.abscissaRange[:-1]
1569 absc = dataOut.abscissaRange[:-1]
1435 noise = dataOut.noise
1570 noise = dataOut.noise
1436 heightRange = dataOut.getHeiRange()
1571 heightRange = dataOut.getHeiRange()
1437 SNR = dataOut.SNR
1572 SNR = dataOut.SNR
1438
1573
1439 if technique == 'DBS':
1574 if technique == 'DBS':
1440
1575
1441 if kwargs.has_key('dirCosx') and kwargs.has_key('dirCosy'):
1576 if kwargs.has_key('dirCosx') and kwargs.has_key('dirCosy'):
1442 theta_x = numpy.array(kwargs['dirCosx'])
1577 theta_x = numpy.array(kwargs['dirCosx'])
1443 theta_y = numpy.array(kwargs['dirCosy'])
1578 theta_y = numpy.array(kwargs['dirCosy'])
1444 else:
1579 else:
1445 elev = numpy.array(kwargs['elevation'])
1580 elev = numpy.array(kwargs['elevation'])
1446 azim = numpy.array(kwargs['azimuth'])
1581 azim = numpy.array(kwargs['azimuth'])
1447 theta_x, theta_y = self.__calculateCosDir(elev, azim)
1582 theta_x, theta_y = self.__calculateCosDir(elev, azim)
1448 azimuth = kwargs['correctAzimuth']
1583 azimuth = kwargs['correctAzimuth']
1449 if kwargs.has_key('horizontalOnly'):
1584 if kwargs.has_key('horizontalOnly'):
1450 horizontalOnly = kwargs['horizontalOnly']
1585 horizontalOnly = kwargs['horizontalOnly']
1451 else: horizontalOnly = False
1586 else: horizontalOnly = False
1452 if kwargs.has_key('correctFactor'):
1587 if kwargs.has_key('correctFactor'):
1453 correctFactor = kwargs['correctFactor']
1588 correctFactor = kwargs['correctFactor']
1454 else: correctFactor = 1
1589 else: correctFactor = 1
1455 if kwargs.has_key('channelList'):
1590 if kwargs.has_key('channelList'):
1456 channelList = kwargs['channelList']
1591 channelList = kwargs['channelList']
1457 if len(channelList) == 2:
1592 if len(channelList) == 2:
1458 horizontalOnly = True
1593 horizontalOnly = True
1459 arrayChannel = numpy.array(channelList)
1594 arrayChannel = numpy.array(channelList)
1460 param = param[arrayChannel,:,:]
1595 param = param[arrayChannel,:,:]
1461 theta_x = theta_x[arrayChannel]
1596 theta_x = theta_x[arrayChannel]
1462 theta_y = theta_y[arrayChannel]
1597 theta_y = theta_y[arrayChannel]
1463
1598
1464 velRadial0 = param[:,1,:] #Radial velocity
1599 velRadial0 = param[:,1,:] #Radial velocity
1465 dataOut.winds, dataOut.heightRange, dataOut.SNR = self.techniqueDBS(velRadial0, theta_x, theta_y, azimuth, correctFactor, horizontalOnly, heightRange, SNR) #DBS Function
1600 dataOut.data_output, dataOut.heightRange, dataOut.SNR = self.techniqueDBS(velRadial0, theta_x, theta_y, azimuth, correctFactor, horizontalOnly, heightRange, SNR) #DBS Function
1466
1601
1467 elif technique == 'SA':
1602 elif technique == 'SA':
1468
1603
1469 #Parameters
1604 #Parameters
1470 position_x = kwargs['positionX']
1605 position_x = kwargs['positionX']
1471 position_y = kwargs['positionY']
1606 position_y = kwargs['positionY']
1472 azimuth = kwargs['azimuth']
1607 azimuth = kwargs['azimuth']
1473
1608
1474 if kwargs.has_key('crosspairsList'):
1609 if kwargs.has_key('crosspairsList'):
1475 pairs = kwargs['crosspairsList']
1610 pairs = kwargs['crosspairsList']
1476 else:
1611 else:
1477 pairs = None
1612 pairs = None
1478
1613
1479 if kwargs.has_key('correctFactor'):
1614 if kwargs.has_key('correctFactor'):
1480 correctFactor = kwargs['correctFactor']
1615 correctFactor = kwargs['correctFactor']
1481 else:
1616 else:
1482 correctFactor = 1
1617 correctFactor = 1
1483
1618
1484 tau = dataOut.data_param
1619 tau = dataOut.data_param
1485 _lambda = dataOut.C/dataOut.frequency
1620 _lambda = dataOut.C/dataOut.frequency
1486 pairsList = dataOut.pairsList
1621 pairsList = dataOut.groupList
1487 nChannels = dataOut.nChannels
1622 nChannels = dataOut.nChannels
1488
1623
1489 dataOut.winds = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor)
1624 dataOut.data_output = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor)
1490 dataOut.initUtcTime = dataOut.ltctime
1625 dataOut.initUtcTime = dataOut.ltctime
1491 dataOut.windsInterval = dataOut.timeInterval
1626 dataOut.outputInterval = dataOut.timeInterval
1492
1627
1493 elif technique == 'Meteors':
1628 elif technique == 'Meteors':
1494 dataOut.flagNoData = True
1629 dataOut.flagNoData = True
1495 self.__dataReady = False
1630 self.__dataReady = False
1496
1631
1497 if kwargs.has_key('nHours'):
1632 if kwargs.has_key('nHours'):
1498 nHours = kwargs['nHours']
1633 nHours = kwargs['nHours']
1499 else:
1634 else:
1500 nHours = 1
1635 nHours = 1
1501
1636
1502 if kwargs.has_key('meteorsPerBin'):
1637 if kwargs.has_key('meteorsPerBin'):
1503 meteorThresh = kwargs['meteorsPerBin']
1638 meteorThresh = kwargs['meteorsPerBin']
1504 else:
1639 else:
1505 meteorThresh = 6
1640 meteorThresh = 6
1506
1641
1507 if kwargs.has_key('hmin'):
1642 if kwargs.has_key('hmin'):
1508 hmin = kwargs['hmin']
1643 hmin = kwargs['hmin']
1509 else: hmin = 70
1644 else: hmin = 70
1510 if kwargs.has_key('hmax'):
1645 if kwargs.has_key('hmax'):
1511 hmax = kwargs['hmax']
1646 hmax = kwargs['hmax']
1512 else: hmax = 110
1647 else: hmax = 110
1513
1648
1514 dataOut.windsInterval = nHours*3600
1649 dataOut.outputInterval = nHours*3600
1515
1650
1516 if self.__isConfig == False:
1651 if self.__isConfig == False:
1517 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
1652 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
1518 #Get Initial LTC time
1653 #Get Initial LTC time
1519 self.__initime = (dataOut.datatime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
1654 self.__initime = (dataOut.datatime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
1520 self.__isConfig = True
1655 self.__isConfig = True
1521
1656
1522 if self.__buffer == None:
1657 if self.__buffer == None:
1523 self.__buffer = dataOut.data_param
1658 self.__buffer = dataOut.data_param
1524 self.__firstdata = copy.copy(dataOut)
1659 self.__firstdata = copy.copy(dataOut)
1525
1660
1526 else:
1661 else:
1527 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
1662 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
1528
1663
1529 self.__checkTime(dataOut.ltctime, dataOut.paramInterval, dataOut.windsInterval) #Check if the buffer is ready
1664 self.__checkTime(dataOut.ltctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
1530
1665
1531 if self.__dataReady:
1666 if self.__dataReady:
1532 dataOut.initUtcTime = self.__initime
1667 dataOut.initUtcTime = self.__initime
1533 self.__initime = self.__initime + dataOut.windsInterval #to erase time offset
1668 self.__initime = self.__initime + dataOut.outputInterval #to erase time offset
1534
1669
1535 dataOut.winds, dataOut.heightRange = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax)
1670 dataOut.data_output, dataOut.heightRange = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax)
1536 dataOut.flagNoData = False
1671 dataOut.flagNoData = False
1537 self.__buffer = None
1672 self.__buffer = None
1538
1673
1539 return No newline at end of file
1674 return
1675
1676 class EWDriftsEstimation(Operation):
1677
1678
1679 def __init__(self):
1680 Operation.__init__(self)
1681
1682 def __correctValues(self, heiRang, phi, velRadial, SNR):
1683 listPhi = phi.tolist()
1684 maxid = listPhi.index(max(listPhi))
1685 minid = listPhi.index(min(listPhi))
1686
1687 rango = range(len(phi))
1688 # rango = numpy.delete(rango,maxid)
1689
1690 heiRang1 = heiRang*math.cos(phi[maxid])
1691 heiRangAux = heiRang*math.cos(phi[minid])
1692 indOut = (heiRang1 < heiRangAux[0]).nonzero()
1693 heiRang1 = numpy.delete(heiRang1,indOut)
1694
1695 velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
1696 SNR1 = numpy.zeros([len(phi),len(heiRang1)])
1697
1698 for i in rango:
1699 x = heiRang*math.cos(phi[i])
1700 y1 = velRadial[i,:]
1701 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
1702
1703 x1 = heiRang1
1704 y11 = f1(x1)
1705
1706 y2 = SNR[i,:]
1707 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
1708 y21 = f2(x1)
1709
1710 velRadial1[i,:] = y11
1711 SNR1[i,:] = y21
1712
1713 return heiRang1, velRadial1, SNR1
1714
1715 def run(self, dataOut, zenith, zenithCorrection):
1716 heiRang = dataOut.heightList
1717 velRadial = dataOut.data_param[:,3,:]
1718 SNR = dataOut.SNR
1719
1720 zenith = numpy.array(zenith)
1721 zenith -= zenithCorrection
1722 zenith *= numpy.pi/180
1723
1724 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR)
1725
1726 alp = zenith[0]
1727 bet = zenith[1]
1728
1729 w_w = velRadial1[0,:]
1730 w_e = velRadial1[1,:]
1731
1732 w = (w_w*numpy.sin(bet) - w_e*numpy.sin(alp))/(numpy.cos(alp)*numpy.sin(bet) - numpy.cos(bet)*numpy.sin(alp))
1733 u = (w_w*numpy.cos(bet) - w_e*numpy.cos(alp))/(numpy.sin(alp)*numpy.cos(bet) - numpy.sin(bet)*numpy.cos(alp))
1734
1735 winds = numpy.vstack((u,w))
1736
1737 dataOut.heightList = heiRang1
1738 dataOut.data_output = winds
1739 dataOut.SNR = SNR1
1740
1741 dataOut.initUtcTime = dataOut.ltctime
1742 dataOut.outputInterval = dataOut.timeInterval
1743 return
1744
1745
1746
1747
1748
1749 No newline at end of file
General Comments 0
You need to be logged in to leave comments. Login now