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