##// END OF EJS Templates
-Modifications to DBS module to accept azimuth and elevation angles as inputs....
Julio Valdez -
r503:ea980322a40d
parent child
Show More
@@ -1,579 +1,580
1 import os
1 import os
2 import datetime
2 import datetime
3 import numpy
3 import numpy
4
4
5 from figure import Figure, isRealtime
5 from figure import Figure, isRealtime
6
6
7 class MomentsPlot(Figure):
7 class MomentsPlot(Figure):
8
8
9 isConfig = None
9 isConfig = None
10 __nsubplots = None
10 __nsubplots = None
11
11
12 WIDTHPROF = None
12 WIDTHPROF = None
13 HEIGHTPROF = None
13 HEIGHTPROF = None
14 PREFIX = 'prm'
14 PREFIX = 'prm'
15
15
16 def __init__(self):
16 def __init__(self):
17
17
18 self.isConfig = False
18 self.isConfig = False
19 self.__nsubplots = 1
19 self.__nsubplots = 1
20
20
21 self.WIDTH = 280
21 self.WIDTH = 280
22 self.HEIGHT = 250
22 self.HEIGHT = 250
23 self.WIDTHPROF = 120
23 self.WIDTHPROF = 120
24 self.HEIGHTPROF = 0
24 self.HEIGHTPROF = 0
25 self.counter_imagwr = 0
25 self.counter_imagwr = 0
26
26
27 self.PLOT_CODE = 1
27 self.PLOT_CODE = 1
28 self.FTP_WEI = None
28 self.FTP_WEI = None
29 self.EXP_CODE = None
29 self.EXP_CODE = None
30 self.SUB_EXP_CODE = None
30 self.SUB_EXP_CODE = None
31 self.PLOT_POS = None
31 self.PLOT_POS = None
32
32
33 def getSubplots(self):
33 def getSubplots(self):
34
34
35 ncol = int(numpy.sqrt(self.nplots)+0.9)
35 ncol = int(numpy.sqrt(self.nplots)+0.9)
36 nrow = int(self.nplots*1./ncol + 0.9)
36 nrow = int(self.nplots*1./ncol + 0.9)
37
37
38 return nrow, ncol
38 return nrow, ncol
39
39
40 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
40 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
41
41
42 self.__showprofile = showprofile
42 self.__showprofile = showprofile
43 self.nplots = nplots
43 self.nplots = nplots
44
44
45 ncolspan = 1
45 ncolspan = 1
46 colspan = 1
46 colspan = 1
47 if showprofile:
47 if showprofile:
48 ncolspan = 3
48 ncolspan = 3
49 colspan = 2
49 colspan = 2
50 self.__nsubplots = 2
50 self.__nsubplots = 2
51
51
52 self.createFigure(id = id,
52 self.createFigure(id = id,
53 wintitle = wintitle,
53 wintitle = wintitle,
54 widthplot = self.WIDTH + self.WIDTHPROF,
54 widthplot = self.WIDTH + self.WIDTHPROF,
55 heightplot = self.HEIGHT + self.HEIGHTPROF,
55 heightplot = self.HEIGHT + self.HEIGHTPROF,
56 show=show)
56 show=show)
57
57
58 nrow, ncol = self.getSubplots()
58 nrow, ncol = self.getSubplots()
59
59
60 counter = 0
60 counter = 0
61 for y in range(nrow):
61 for y in range(nrow):
62 for x in range(ncol):
62 for x in range(ncol):
63
63
64 if counter >= self.nplots:
64 if counter >= self.nplots:
65 break
65 break
66
66
67 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
67 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
68
68
69 if showprofile:
69 if showprofile:
70 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
70 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
71
71
72 counter += 1
72 counter += 1
73
73
74 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=True,
74 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=True,
75 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
75 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
76 save=False, figpath='', figfile=None, show=True, ftp=False, wr_period=1,
76 save=False, figpath='', figfile=None, show=True, ftp=False, wr_period=1,
77 server=None, folder=None, username=None, password=None,
77 server=None, folder=None, username=None, password=None,
78 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False):
78 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False):
79
79
80 """
80 """
81
81
82 Input:
82 Input:
83 dataOut :
83 dataOut :
84 id :
84 id :
85 wintitle :
85 wintitle :
86 channelList :
86 channelList :
87 showProfile :
87 showProfile :
88 xmin : None,
88 xmin : None,
89 xmax : None,
89 xmax : None,
90 ymin : None,
90 ymin : None,
91 ymax : None,
91 ymax : None,
92 zmin : None,
92 zmin : None,
93 zmax : None
93 zmax : None
94 """
94 """
95
95
96 if dataOut.flagNoData:
96 if dataOut.flagNoData:
97 return None
97 return None
98
98
99 if realtime:
99 if realtime:
100 if not(isRealtime(utcdatatime = dataOut.utctime)):
100 if not(isRealtime(utcdatatime = dataOut.utctime)):
101 print 'Skipping this plot function'
101 print 'Skipping this plot function'
102 return
102 return
103
103
104 if channelList == None:
104 if channelList == None:
105 channelIndexList = dataOut.channelIndexList
105 channelIndexList = dataOut.channelIndexList
106 else:
106 else:
107 channelIndexList = []
107 channelIndexList = []
108 for channel in channelList:
108 for channel in channelList:
109 if channel not in dataOut.channelList:
109 if channel not in dataOut.channelList:
110 raise ValueError, "Channel %d is not in dataOut.channelList"
110 raise ValueError, "Channel %d is not in dataOut.channelList"
111 channelIndexList.append(dataOut.channelList.index(channel))
111 channelIndexList.append(dataOut.channelList.index(channel))
112
112
113 factor = dataOut.normFactor
113 factor = dataOut.normFactor
114 x = dataOut.abscissaRange
114 x = dataOut.abscissaRange
115 y = dataOut.heightRange
115 y = dataOut.heightRange
116
116
117 z = dataOut.data_pre[channelIndexList,:,:]/factor
117 z = dataOut.data_pre[channelIndexList,:,:]/factor
118 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
118 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
119 avg = numpy.average(z, axis=1)
119 avg = numpy.average(z, axis=1)
120 noise = dataOut.noise/factor
120 noise = dataOut.noise/factor
121
121
122 zdB = 10*numpy.log10(z)
122 zdB = 10*numpy.log10(z)
123 avgdB = 10*numpy.log10(avg)
123 avgdB = 10*numpy.log10(avg)
124 noisedB = 10*numpy.log10(noise)
124 noisedB = 10*numpy.log10(noise)
125
125
126 #thisDatetime = dataOut.datatime
126 #thisDatetime = dataOut.datatime
127 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
127 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
128 title = wintitle + " Parameters"
128 title = wintitle + " Parameters"
129 xlabel = "Velocity (m/s)"
129 xlabel = "Velocity (m/s)"
130 ylabel = "Range (Km)"
130 ylabel = "Range (Km)"
131
131
132 if not self.isConfig:
132 if not self.isConfig:
133
133
134 nplots = len(channelIndexList)
134 nplots = len(channelIndexList)
135
135
136 self.setup(id=id,
136 self.setup(id=id,
137 nplots=nplots,
137 nplots=nplots,
138 wintitle=wintitle,
138 wintitle=wintitle,
139 showprofile=showprofile,
139 showprofile=showprofile,
140 show=show)
140 show=show)
141
141
142 if xmin == None: xmin = numpy.nanmin(x)
142 if xmin == None: xmin = numpy.nanmin(x)
143 if xmax == None: xmax = numpy.nanmax(x)
143 if xmax == None: xmax = numpy.nanmax(x)
144 if ymin == None: ymin = numpy.nanmin(y)
144 if ymin == None: ymin = numpy.nanmin(y)
145 if ymax == None: ymax = numpy.nanmax(y)
145 if ymax == None: ymax = numpy.nanmax(y)
146 if zmin == None: zmin = numpy.nanmin(avgdB)*0.9
146 if zmin == None: zmin = numpy.nanmin(avgdB)*0.9
147 if zmax == None: zmax = numpy.nanmax(avgdB)*0.9
147 if zmax == None: zmax = numpy.nanmax(avgdB)*0.9
148
148
149 self.FTP_WEI = ftp_wei
149 self.FTP_WEI = ftp_wei
150 self.EXP_CODE = exp_code
150 self.EXP_CODE = exp_code
151 self.SUB_EXP_CODE = sub_exp_code
151 self.SUB_EXP_CODE = sub_exp_code
152 self.PLOT_POS = plot_pos
152 self.PLOT_POS = plot_pos
153
153
154 self.isConfig = True
154 self.isConfig = True
155
155
156 self.setWinTitle(title)
156 self.setWinTitle(title)
157
157
158 for i in range(self.nplots):
158 for i in range(self.nplots):
159 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
159 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
160 title = "Channel %d: %4.2fdB: %s" %(dataOut.channelList[i]+1, noisedB[i], str_datetime)
160 title = "Channel %d: %4.2fdB: %s" %(dataOut.channelList[i]+1, noisedB[i], str_datetime)
161 axes = self.axesList[i*self.__nsubplots]
161 axes = self.axesList[i*self.__nsubplots]
162 axes.pcolor(x, y, zdB[i,:,:],
162 axes.pcolor(x, y, zdB[i,:,:],
163 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
163 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
164 xlabel=xlabel, ylabel=ylabel, title=title,
164 xlabel=xlabel, ylabel=ylabel, title=title,
165 ticksize=9, cblabel='')
165 ticksize=9, cblabel='')
166 #Mean Line
166 #Mean Line
167 mean = dataOut.data_param[i, 1, :]
167 mean = dataOut.data_param[i, 1, :]
168 axes.addpline(mean, y, idline=0, color="black", linestyle="solid", lw=1)
168 axes.addpline(mean, y, idline=0, color="black", linestyle="solid", lw=1)
169
169
170 if self.__showprofile:
170 if self.__showprofile:
171 axes = self.axesList[i*self.__nsubplots +1]
171 axes = self.axesList[i*self.__nsubplots +1]
172 axes.pline(avgdB[i], y,
172 axes.pline(avgdB[i], y,
173 xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
173 xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
174 xlabel='dB', ylabel='', title='',
174 xlabel='dB', ylabel='', title='',
175 ytick_visible=False,
175 ytick_visible=False,
176 grid='x')
176 grid='x')
177
177
178 noiseline = numpy.repeat(noisedB[i], len(y))
178 noiseline = numpy.repeat(noisedB[i], len(y))
179 axes.addpline(noiseline, y, idline=1, color="black", linestyle="dashed", lw=2)
179 axes.addpline(noiseline, y, idline=1, color="black", linestyle="dashed", lw=2)
180
180
181 self.draw()
181 self.draw()
182
182
183 if figfile == None:
183 if figfile == None:
184 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
184 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
185 figfile = self.getFilename(name = str_datetime)
185 figfile = self.getFilename(name = str_datetime)
186
186
187 if figpath != '':
187 if figpath != '':
188 self.counter_imagwr += 1
188 self.counter_imagwr += 1
189 if (self.counter_imagwr>=wr_period):
189 if (self.counter_imagwr>=wr_period):
190 # store png plot to local folder
190 # store png plot to local folder
191 self.saveFigure(figpath, figfile)
191 self.saveFigure(figpath, figfile)
192 # store png plot to FTP server according to RT-Web format
192 # store png plot to FTP server according to RT-Web format
193 name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
193 name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
194 ftp_filename = os.path.join(figpath, name)
194 ftp_filename = os.path.join(figpath, name)
195 self.saveFigure(figpath, ftp_filename)
195 self.saveFigure(figpath, ftp_filename)
196 self.counter_imagwr = 0
196 self.counter_imagwr = 0
197
197
198 class SkyMapPlot(Figure):
198 class SkyMapPlot(Figure):
199
199
200 __isConfig = None
200 __isConfig = None
201 __nsubplots = None
201 __nsubplots = None
202
202
203 WIDTHPROF = None
203 WIDTHPROF = None
204 HEIGHTPROF = None
204 HEIGHTPROF = None
205 PREFIX = 'prm'
205 PREFIX = 'prm'
206
206
207 def __init__(self):
207 def __init__(self):
208
208
209 self.__isConfig = False
209 self.__isConfig = False
210 self.__nsubplots = 1
210 self.__nsubplots = 1
211
211
212 # self.WIDTH = 280
212 # self.WIDTH = 280
213 # self.HEIGHT = 250
213 # self.HEIGHT = 250
214 self.WIDTH = 600
214 self.WIDTH = 600
215 self.HEIGHT = 600
215 self.HEIGHT = 600
216 self.WIDTHPROF = 120
216 self.WIDTHPROF = 120
217 self.HEIGHTPROF = 0
217 self.HEIGHTPROF = 0
218 self.counter_imagwr = 0
218 self.counter_imagwr = 0
219
219
220 self.PLOT_CODE = 1
220 self.PLOT_CODE = 1
221 self.FTP_WEI = None
221 self.FTP_WEI = None
222 self.EXP_CODE = None
222 self.EXP_CODE = None
223 self.SUB_EXP_CODE = None
223 self.SUB_EXP_CODE = None
224 self.PLOT_POS = None
224 self.PLOT_POS = None
225
225
226 def getSubplots(self):
226 def getSubplots(self):
227
227
228 ncol = int(numpy.sqrt(self.nplots)+0.9)
228 ncol = int(numpy.sqrt(self.nplots)+0.9)
229 nrow = int(self.nplots*1./ncol + 0.9)
229 nrow = int(self.nplots*1./ncol + 0.9)
230
230
231 return nrow, ncol
231 return nrow, ncol
232
232
233 def setup(self, id, nplots, wintitle, showprofile=False, show=True):
233 def setup(self, id, nplots, wintitle, showprofile=False, show=True):
234
234
235 self.__showprofile = showprofile
235 self.__showprofile = showprofile
236 self.nplots = nplots
236 self.nplots = nplots
237
237
238 ncolspan = 1
238 ncolspan = 1
239 colspan = 1
239 colspan = 1
240
240
241 self.createFigure(id = id,
241 self.createFigure(id = id,
242 wintitle = wintitle,
242 wintitle = wintitle,
243 widthplot = self.WIDTH, #+ self.WIDTHPROF,
243 widthplot = self.WIDTH, #+ self.WIDTHPROF,
244 heightplot = self.HEIGHT,# + self.HEIGHTPROF,
244 heightplot = self.HEIGHT,# + self.HEIGHTPROF,
245 show=show)
245 show=show)
246
246
247 nrow, ncol = 1,1
247 nrow, ncol = 1,1
248 counter = 0
248 counter = 0
249 x = 0
249 x = 0
250 y = 0
250 y = 0
251 self.addAxes(1, 1, 0, 0, 1, 1, True)
251 self.addAxes(1, 1, 0, 0, 1, 1, True)
252
252
253 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=False,
253 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=False,
254 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
254 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
255 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
255 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
256 server=None, folder=None, username=None, password=None,
256 server=None, folder=None, username=None, password=None,
257 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False):
257 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False):
258
258
259 """
259 """
260
260
261 Input:
261 Input:
262 dataOut :
262 dataOut :
263 id :
263 id :
264 wintitle :
264 wintitle :
265 channelList :
265 channelList :
266 showProfile :
266 showProfile :
267 xmin : None,
267 xmin : None,
268 xmax : None,
268 xmax : None,
269 ymin : None,
269 ymin : None,
270 ymax : None,
270 ymax : None,
271 zmin : None,
271 zmin : None,
272 zmax : None
272 zmax : None
273 """
273 """
274
274
275 arrayParameters = dataOut.data_param
275 arrayParameters = dataOut.data_param
276 error = arrayParameters[:,-1]
276 error = arrayParameters[:,-1]
277 indValid = numpy.where(error == 0)[0]
277 indValid = numpy.where(error == 0)[0]
278 finalMeteor = arrayParameters[indValid,:]
278 finalMeteor = arrayParameters[indValid,:]
279 finalAzimuth = finalMeteor[:,4]
279 finalAzimuth = finalMeteor[:,4]
280 finalZenith = finalMeteor[:,5]
280 finalZenith = finalMeteor[:,5]
281
281
282 x = finalAzimuth*numpy.pi/180
282 x = finalAzimuth*numpy.pi/180
283 y = finalZenith
283 y = finalZenith
284
284
285
285
286 #thisDatetime = dataOut.datatime
286 #thisDatetime = dataOut.datatime
287 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
287 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
288 title = wintitle + " Parameters"
288 title = wintitle + " Parameters"
289 xlabel = "Zonal Zenith Angle (deg) "
289 xlabel = "Zonal Zenith Angle (deg) "
290 ylabel = "Meridional Zenith Angle (deg)"
290 ylabel = "Meridional Zenith Angle (deg)"
291
291
292 if not self.__isConfig:
292 if not self.__isConfig:
293
293
294 nplots = 1
294 nplots = 1
295
295
296 self.setup(id=id,
296 self.setup(id=id,
297 nplots=nplots,
297 nplots=nplots,
298 wintitle=wintitle,
298 wintitle=wintitle,
299 showprofile=showprofile,
299 showprofile=showprofile,
300 show=show)
300 show=show)
301
301
302 self.FTP_WEI = ftp_wei
302 self.FTP_WEI = ftp_wei
303 self.EXP_CODE = exp_code
303 self.EXP_CODE = exp_code
304 self.SUB_EXP_CODE = sub_exp_code
304 self.SUB_EXP_CODE = sub_exp_code
305 self.PLOT_POS = plot_pos
305 self.PLOT_POS = plot_pos
306 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
306 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
307 self.firstdate = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
307 self.firstdate = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
308 self.__isConfig = True
308 self.__isConfig = True
309
309
310 self.setWinTitle(title)
310 self.setWinTitle(title)
311
311
312 i = 0
312 i = 0
313 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
313 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
314
314
315 axes = self.axesList[i*self.__nsubplots]
315 axes = self.axesList[i*self.__nsubplots]
316 nevents = axes.x_buffer.shape[0] + x.shape[0]
316 nevents = axes.x_buffer.shape[0] + x.shape[0]
317 title = "Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n" %(self.firstdate,str_datetime,nevents)
317 title = "Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n" %(self.firstdate,str_datetime,nevents)
318 axes.polar(x, y,
318 axes.polar(x, y,
319 title=title, xlabel=xlabel, ylabel=ylabel,
319 title=title, xlabel=xlabel, ylabel=ylabel,
320 ticksize=9, cblabel='')
320 ticksize=9, cblabel='')
321
321
322 self.draw()
322 self.draw()
323
323
324 if save:
324 if save:
325
325
326 self.counter_imagwr += 1
326 self.counter_imagwr += 1
327 if (self.counter_imagwr==wr_period):
327 if (self.counter_imagwr==wr_period):
328
328
329 if figfile == None:
329 if figfile == None:
330 figfile = self.getFilename(name = self.name)
330 figfile = self.getFilename(name = self.name)
331 self.saveFigure(figpath, figfile)
331 self.saveFigure(figpath, figfile)
332
332
333 if ftp:
333 if ftp:
334 #provisionalmente envia archivos en el formato de la web en tiempo real
334 #provisionalmente envia archivos en el formato de la web en tiempo real
335 name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
335 name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
336 path = '%s%03d' %(self.PREFIX, self.id)
336 path = '%s%03d' %(self.PREFIX, self.id)
337 ftp_file = os.path.join(path,'ftp','%s.png'%name)
337 ftp_file = os.path.join(path,'ftp','%s.png'%name)
338 self.saveFigure(figpath, ftp_file)
338 self.saveFigure(figpath, ftp_file)
339 ftp_filename = os.path.join(figpath,ftp_file)
339 ftp_filename = os.path.join(figpath,ftp_file)
340
340
341
341
342 try:
342 try:
343 self.sendByFTP(ftp_filename, server, folder, username, password)
343 self.sendByFTP(ftp_filename, server, folder, username, password)
344 except:
344 except:
345 self.counter_imagwr = 0
345 self.counter_imagwr = 0
346 raise ValueError, 'Error FTP'
346 raise ValueError, 'Error FTP'
347
347
348 self.counter_imagwr = 0
348 self.counter_imagwr = 0
349
349
350
350
351 class WindProfilerPlot(Figure):
351 class WindProfilerPlot(Figure):
352
352
353 __isConfig = None
353 __isConfig = None
354 __nsubplots = None
354 __nsubplots = None
355
355
356 WIDTHPROF = None
356 WIDTHPROF = None
357 HEIGHTPROF = None
357 HEIGHTPROF = None
358 PREFIX = 'wind'
358 PREFIX = 'wind'
359
359
360 def __init__(self):
360 def __init__(self):
361
361
362 self.timerange = 2*60*60
362 self.timerange = 2*60*60
363 self.isConfig = False
363 self.isConfig = False
364 self.__nsubplots = 1
364 self.__nsubplots = 1
365
365
366 self.WIDTH = 800
366 self.WIDTH = 800
367 self.HEIGHT = 150
367 self.HEIGHT = 150
368 self.WIDTHPROF = 120
368 self.WIDTHPROF = 120
369 self.HEIGHTPROF = 0
369 self.HEIGHTPROF = 0
370 self.counter_imagwr = 0
370 self.counter_imagwr = 0
371
371
372 self.PLOT_CODE = 0
372 self.PLOT_CODE = 0
373 self.FTP_WEI = None
373 self.FTP_WEI = None
374 self.EXP_CODE = None
374 self.EXP_CODE = None
375 self.SUB_EXP_CODE = None
375 self.SUB_EXP_CODE = None
376 self.PLOT_POS = None
376 self.PLOT_POS = None
377 self.tmin = None
377 self.tmin = None
378 self.tmax = None
378 self.tmax = None
379
379
380 self.xmin = None
380 self.xmin = None
381 self.xmax = None
381 self.xmax = None
382
382
383 self.figfile = None
383 self.figfile = None
384
384
385 def getSubplots(self):
385 def getSubplots(self):
386
386
387 ncol = 1
387 ncol = 1
388 nrow = self.nplots
388 nrow = self.nplots
389
389
390 return nrow, ncol
390 return nrow, ncol
391
391
392 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
392 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
393
393
394 self.__showprofile = showprofile
394 self.__showprofile = showprofile
395 self.nplots = nplots
395 self.nplots = nplots
396
396
397 ncolspan = 1
397 ncolspan = 1
398 colspan = 1
398 colspan = 1
399
399
400 self.createFigure(id = id,
400 self.createFigure(id = id,
401 wintitle = wintitle,
401 wintitle = wintitle,
402 widthplot = self.WIDTH + self.WIDTHPROF,
402 widthplot = self.WIDTH + self.WIDTHPROF,
403 heightplot = self.HEIGHT + self.HEIGHTPROF,
403 heightplot = self.HEIGHT + self.HEIGHTPROF,
404 show=show)
404 show=show)
405
405
406 nrow, ncol = self.getSubplots()
406 nrow, ncol = self.getSubplots()
407
407
408 counter = 0
408 counter = 0
409 for y in range(nrow):
409 for y in range(nrow):
410 if counter >= self.nplots:
410 if counter >= self.nplots:
411 break
411 break
412
412
413 self.addAxes(nrow, ncol*ncolspan, y, 0, colspan, 1)
413 self.addAxes(nrow, ncol*ncolspan, y, 0, colspan, 1)
414 counter += 1
414 counter += 1
415
415
416 def run(self, dataOut, id, wintitle="", channelList=None,
416 def run(self, dataOut, id, wintitle="", channelList=None,
417 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
417 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
418 zmax_ver = None, zmin_ver = None, SNRmin = None, SNRmax = None,
418 zmax_ver = None, zmin_ver = None, SNRmin = None, SNRmax = None,
419 timerange=None, SNRthresh = None,
419 timerange=None, SNRthresh = None,
420 save=False, figpath='', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
420 save=False, figpath='', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
421 server=None, folder=None, username=None, password=None,
421 server=None, folder=None, username=None, password=None,
422 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
422 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
423 """
423 """
424
424
425 Input:
425 Input:
426 dataOut :
426 dataOut :
427 id :
427 id :
428 wintitle :
428 wintitle :
429 channelList :
429 channelList :
430 showProfile :
430 showProfile :
431 xmin : None,
431 xmin : None,
432 xmax : None,
432 xmax : None,
433 ymin : None,
433 ymin : None,
434 ymax : None,
434 ymax : None,
435 zmin : None,
435 zmin : None,
436 zmax : None
436 zmax : None
437 """
437 """
438
438
439 if channelList == None:
439 if channelList == None:
440 channelIndexList = dataOut.channelIndexList
440 channelIndexList = dataOut.channelIndexList
441 else:
441 else:
442 channelIndexList = []
442 channelIndexList = []
443 for channel in channelList:
443 for channel in channelList:
444 if channel not in dataOut.channelList:
444 if channel not in dataOut.channelList:
445 raise ValueError, "Channel %d is not in dataOut.channelList"
445 raise ValueError, "Channel %d is not in dataOut.channelList"
446 channelIndexList.append(dataOut.channelList.index(channel))
446 channelIndexList.append(dataOut.channelList.index(channel))
447
447
448 if timerange != None:
448 if timerange != None:
449 self.timerange = timerange
449 self.timerange = timerange
450
450
451 tmin = None
451 tmin = None
452 tmax = None
452 tmax = None
453
453
454 x = dataOut.getTimeRange1()
454 x = dataOut.getTimeRange1()
455 # y = dataOut.heightRange
455 y = dataOut.heightRange
456 y = dataOut.heightRange
456
457
457 z = dataOut.winds
458 z = dataOut.winds
458 nplots = z.shape[0] #Number of wind dimensions estimated
459 nplots = z.shape[0] #Number of wind dimensions estimated
459 nplotsw = nplots
460 nplotsw = nplots
460
461
461 #If there is a SNR function defined
462 #If there is a SNR function defined
462 if dataOut.SNR != None:
463 if dataOut.SNR != None:
463 nplots += 1
464 nplots += 1
464 SNR = dataOut.SNR
465 SNR = dataOut.SNR
465 SNRavg = numpy.average(SNR, axis=0)
466 SNRavg = numpy.average(SNR, axis=0)
466
467
467 SNRdB = 10*numpy.log10(SNR)
468 SNRdB = 10*numpy.log10(SNR)
468 SNRavgdB = 10*numpy.log10(SNRavg)
469 SNRavgdB = 10*numpy.log10(SNRavg)
469
470
470 if SNRthresh == None: SNRthresh = -5.0
471 if SNRthresh == None: SNRthresh = -5.0
471 ind = numpy.where(SNRavg < 10**(SNRthresh/10))[0]
472 ind = numpy.where(SNRavg < 10**(SNRthresh/10))[0]
472
473
473 for i in range(nplotsw):
474 for i in range(nplotsw):
474 z[i,ind] = numpy.nan
475 z[i,ind] = numpy.nan
475
476
476
477
477 showprofile = False
478 showprofile = False
478 # thisDatetime = dataOut.datatime
479 # thisDatetime = dataOut.datatime
479 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
480 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
480 title = wintitle + "Wind"
481 title = wintitle + "Wind"
481 xlabel = ""
482 xlabel = ""
482 ylabel = "Range (Km)"
483 ylabel = "Range (Km)"
483
484
484 if not self.__isConfig:
485 if not self.__isConfig:
485
486
486
487
487
488
488 self.setup(id=id,
489 self.setup(id=id,
489 nplots=nplots,
490 nplots=nplots,
490 wintitle=wintitle,
491 wintitle=wintitle,
491 showprofile=showprofile,
492 showprofile=showprofile,
492 show=show)
493 show=show)
493
494
494 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
495 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
495
496
496 if ymin == None: ymin = numpy.nanmin(y)
497 if ymin == None: ymin = numpy.nanmin(y)
497 if ymax == None: ymax = numpy.nanmax(y)
498 if ymax == None: ymax = numpy.nanmax(y)
498
499
499 if zmax == None: zmax = numpy.nanmax(abs(z[range(2),:]))
500 if zmax == None: zmax = numpy.nanmax(abs(z[range(2),:]))
500 #if numpy.isnan(zmax): zmax = 50
501 #if numpy.isnan(zmax): zmax = 50
501 if zmin == None: zmin = -zmax
502 if zmin == None: zmin = -zmax
502
503
503 if nplotsw == 3:
504 if nplotsw == 3:
504 if zmax_ver == None: zmax_ver = numpy.nanmax(abs(z[2,:]))
505 if zmax_ver == None: zmax_ver = numpy.nanmax(abs(z[2,:]))
505 if zmin_ver == None: zmin_ver = -zmax_ver
506 if zmin_ver == None: zmin_ver = -zmax_ver
506
507
507 if dataOut.SNR != None:
508 if dataOut.SNR != None:
508 if SNRmin == None: SNRmin = numpy.nanmin(SNRavgdB)
509 if SNRmin == None: SNRmin = numpy.nanmin(SNRavgdB)
509 if SNRmax == None: SNRmax = numpy.nanmax(SNRavgdB)
510 if SNRmax == None: SNRmax = numpy.nanmax(SNRavgdB)
510
511
511 self.FTP_WEI = ftp_wei
512 self.FTP_WEI = ftp_wei
512 self.EXP_CODE = exp_code
513 self.EXP_CODE = exp_code
513 self.SUB_EXP_CODE = sub_exp_code
514 self.SUB_EXP_CODE = sub_exp_code
514 self.PLOT_POS = plot_pos
515 self.PLOT_POS = plot_pos
515
516
516 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
517 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
517 self.__isConfig = True
518 self.__isConfig = True
518
519
519
520
520 self.setWinTitle(title)
521 self.setWinTitle(title)
521
522
522 if ((self.xmax - x[1]) < (x[1]-x[0])):
523 if ((self.xmax - x[1]) < (x[1]-x[0])):
523 x[1] = self.xmax
524 x[1] = self.xmax
524
525
525 strWind = ['Zonal', 'Meridional', 'Vertical']
526 strWind = ['Zonal', 'Meridional', 'Vertical']
526 strCb = ['Velocity (m/s)','Velocity (m/s)','Velocity (cm/s)']
527 strCb = ['Velocity (m/s)','Velocity (m/s)','Velocity (cm/s)']
527 zmaxVector = [zmax, zmax, zmax_ver]
528 zmaxVector = [zmax, zmax, zmax_ver]
528 zminVector = [zmin, zmin, zmin_ver]
529 zminVector = [zmin, zmin, zmin_ver]
529 windFactor = [1,1,100]
530 windFactor = [1,1,100]
530
531
531 for i in range(nplotsw):
532 for i in range(nplotsw):
532
533
533 title = "%s Wind: %s" %(strWind[i], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
534 title = "%s Wind: %s" %(strWind[i], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
534 axes = self.axesList[i*self.__nsubplots]
535 axes = self.axesList[i*self.__nsubplots]
535
536
536 z1 = z[i,:].reshape((1,-1))*windFactor[i]
537 z1 = z[i,:].reshape((1,-1))*windFactor[i]
537
538
538 axes.pcolorbuffer(x, y, z1,
539 axes.pcolorbuffer(x, y, z1,
539 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zminVector[i], zmax=zmaxVector[i],
540 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zminVector[i], zmax=zmaxVector[i],
540 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
541 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
541 ticksize=9, cblabel=strCb[i], cbsize="1%", colormap="RdBu_r" )
542 ticksize=9, cblabel=strCb[i], cbsize="1%", colormap="RdBu_r" )
542
543
543 if dataOut.SNR != None:
544 if dataOut.SNR != None:
544 i += 1
545 i += 1
545 title = "Signal Noise Ratio (SNR): %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
546 title = "Signal Noise Ratio (SNR): %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
546 axes = self.axesList[i*self.__nsubplots]
547 axes = self.axesList[i*self.__nsubplots]
547
548
548 SNRavgdB = SNRavgdB.reshape((1,-1))
549 SNRavgdB = SNRavgdB.reshape((1,-1))
549
550
550 axes.pcolorbuffer(x, y, SNRavgdB,
551 axes.pcolorbuffer(x, y, SNRavgdB,
551 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
552 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
552 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
553 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
553 ticksize=9, cblabel='', cbsize="1%", colormap="jet")
554 ticksize=9, cblabel='', cbsize="1%", colormap="jet")
554
555
555 self.draw()
556 self.draw()
556
557
557 if x[1] >= self.axesList[0].xmax:
558 if x[1] >= self.axesList[0].xmax:
558 self.counter_imagwr = wr_period
559 self.counter_imagwr = wr_period
559 self.__isConfig = False
560 self.__isConfig = False
560
561
561
562
562 if self.figfile == None:
563 if self.figfile == None:
563 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
564 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
564 self.figfile = self.getFilename(name = str_datetime)
565 self.figfile = self.getFilename(name = str_datetime)
565
566
566 if figpath != '':
567 if figpath != '':
567
568
568 self.counter_imagwr += 1
569 self.counter_imagwr += 1
569 if (self.counter_imagwr>=wr_period):
570 if (self.counter_imagwr>=wr_period):
570 # store png plot to local folder
571 # store png plot to local folder
571 self.saveFigure(figpath, self.figfile)
572 self.saveFigure(figpath, self.figfile)
572 # store png plot to FTP server according to RT-Web format
573 # store png plot to FTP server according to RT-Web format
573 name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
574 name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
574 ftp_filename = os.path.join(figpath, name)
575 ftp_filename = os.path.join(figpath, name)
575 self.saveFigure(figpath, ftp_filename)
576 self.saveFigure(figpath, ftp_filename)
576
577
577 self.counter_imagwr = 0
578 self.counter_imagwr = 0
578
579
579 No newline at end of file
580
@@ -1,1521 +1,1539
1 import numpy
1 import numpy
2 import math
2 import math
3 from scipy import optimize
3 from scipy import optimize
4 from scipy import interpolate
4 from scipy import interpolate
5 from scipy import signal
5 from scipy import signal
6 from scipy import stats
6 from scipy import stats
7 import re
7 import re
8 import datetime
8 import datetime
9 import copy
9 import copy
10
10
11
11
12 from jroproc_base import ProcessingUnit, Operation
12 from jroproc_base import ProcessingUnit, Operation
13 from model.data.jrodata import Parameters
13 from model.data.jrodata import Parameters
14
14
15
15
16 class ParametersProc(ProcessingUnit):
16 class ParametersProc(ProcessingUnit):
17
17
18 nSeconds = None
18 nSeconds = None
19
19
20 def __init__(self):
20 def __init__(self):
21 ProcessingUnit.__init__(self)
21 ProcessingUnit.__init__(self)
22
22
23 self.objectDict = {}
23 self.objectDict = {}
24 self.buffer = None
24 self.buffer = None
25 self.firstdatatime = None
25 self.firstdatatime = None
26 self.profIndex = 0
26 self.profIndex = 0
27 self.dataOut = Parameters()
27 self.dataOut = Parameters()
28
28
29 def __updateObjFromInput(self):
29 def __updateObjFromInput(self):
30
30
31 self.dataOut.inputUnit = self.dataIn.type
31 self.dataOut.inputUnit = self.dataIn.type
32
32
33 self.dataOut.timeZone = self.dataIn.timeZone
33 self.dataOut.timeZone = self.dataIn.timeZone
34 self.dataOut.dstFlag = self.dataIn.dstFlag
34 self.dataOut.dstFlag = self.dataIn.dstFlag
35 self.dataOut.errorCount = self.dataIn.errorCount
35 self.dataOut.errorCount = self.dataIn.errorCount
36 self.dataOut.useLocalTime = self.dataIn.useLocalTime
36 self.dataOut.useLocalTime = self.dataIn.useLocalTime
37
37
38 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
38 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
39 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
39 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
40 self.dataOut.channelList = self.dataIn.channelList
40 self.dataOut.channelList = self.dataIn.channelList
41 self.dataOut.heightList = self.dataIn.heightList
41 self.dataOut.heightList = self.dataIn.heightList
42 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
42 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
43 # self.dataOut.nHeights = self.dataIn.nHeights
43 # self.dataOut.nHeights = self.dataIn.nHeights
44 # self.dataOut.nChannels = self.dataIn.nChannels
44 # self.dataOut.nChannels = self.dataIn.nChannels
45 self.dataOut.nBaud = self.dataIn.nBaud
45 self.dataOut.nBaud = self.dataIn.nBaud
46 self.dataOut.nCode = self.dataIn.nCode
46 self.dataOut.nCode = self.dataIn.nCode
47 self.dataOut.code = self.dataIn.code
47 self.dataOut.code = self.dataIn.code
48 # self.dataOut.nProfiles = self.dataOut.nFFTPoints
48 # self.dataOut.nProfiles = self.dataOut.nFFTPoints
49 self.dataOut.flagTimeBlock = self.dataIn.flagTimeBlock
49 self.dataOut.flagTimeBlock = self.dataIn.flagTimeBlock
50 self.dataOut.utctime = self.firstdatatime
50 self.dataOut.utctime = self.firstdatatime
51 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
51 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
52 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
52 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
53 # self.dataOut.nCohInt = self.dataIn.nCohInt
53 # self.dataOut.nCohInt = self.dataIn.nCohInt
54 # self.dataOut.nIncohInt = 1
54 # self.dataOut.nIncohInt = 1
55 self.dataOut.ippSeconds = self.dataIn.ippSeconds
55 self.dataOut.ippSeconds = self.dataIn.ippSeconds
56 # self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
56 # self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
57 self.dataOut.timeInterval = self.dataIn.timeInterval
57 self.dataOut.timeInterval = self.dataIn.timeInterval
58 self.dataOut.heightRange = self.dataIn.getHeiRange()
58 self.dataOut.heightRange = self.dataIn.getHeiRange()
59 self.dataOut.frequency = self.dataIn.frequency
59 self.dataOut.frequency = self.dataIn.frequency
60
60
61 def run(self, nSeconds = None, nProfiles = None):
61 def run(self, nSeconds = None, nProfiles = None):
62
62
63 self.dataOut.flagNoData = True
63 self.dataOut.flagNoData = True
64
64
65 if self.firstdatatime == None:
65 if self.firstdatatime == None:
66 self.firstdatatime = self.dataIn.utctime
66 self.firstdatatime = self.dataIn.utctime
67
67
68 #---------------------- Voltage Data ---------------------------
68 #---------------------- Voltage Data ---------------------------
69
69
70 if self.dataIn.type == "Voltage":
70 if self.dataIn.type == "Voltage":
71 if nSeconds != None:
71 if nSeconds != None:
72 self.nSeconds = nSeconds
72 self.nSeconds = nSeconds
73 self.nProfiles= int(numpy.floor(nSeconds/(self.dataIn.ippSeconds*self.dataIn.nCohInt)))
73 self.nProfiles= int(numpy.floor(nSeconds/(self.dataIn.ippSeconds*self.dataIn.nCohInt)))
74
74
75 if self.buffer == None:
75 if self.buffer == None:
76 self.buffer = numpy.zeros((self.dataIn.nChannels,
76 self.buffer = numpy.zeros((self.dataIn.nChannels,
77 self.nProfiles,
77 self.nProfiles,
78 self.dataIn.nHeights),
78 self.dataIn.nHeights),
79 dtype='complex')
79 dtype='complex')
80
80
81 self.buffer[:,self.profIndex,:] = self.dataIn.data.copy()
81 self.buffer[:,self.profIndex,:] = self.dataIn.data.copy()
82 self.profIndex += 1
82 self.profIndex += 1
83
83
84 if self.profIndex == self.nProfiles:
84 if self.profIndex == self.nProfiles:
85
85
86 self.__updateObjFromInput()
86 self.__updateObjFromInput()
87 self.dataOut.data_pre = self.buffer.copy()
87 self.dataOut.data_pre = self.buffer.copy()
88 self.dataOut.paramInterval = nSeconds
88 self.dataOut.paramInterval = nSeconds
89 self.dataOut.flagNoData = False
89 self.dataOut.flagNoData = False
90
90
91 self.buffer = None
91 self.buffer = None
92 self.firstdatatime = None
92 self.firstdatatime = None
93 self.profIndex = 0
93 self.profIndex = 0
94
94
95 #---------------------- Spectra Data ---------------------------
95 #---------------------- Spectra Data ---------------------------
96
96
97 if self.dataIn.type == "Spectra":
97 if self.dataIn.type == "Spectra":
98 self.dataOut.data_pre = self.dataIn.data_spc.copy()
98 self.dataOut.data_pre = self.dataIn.data_spc.copy()
99 self.dataOut.abscissaRange = self.dataIn.getVelRange(1)
99 self.dataOut.abscissaRange = self.dataIn.getVelRange(1)
100 self.dataOut.noise = self.dataIn.getNoise()
100 self.dataOut.noise = self.dataIn.getNoise()
101 self.dataOut.normFactor = self.dataIn.normFactor
101 self.dataOut.normFactor = self.dataIn.normFactor
102
102
103 self.__updateObjFromInput()
103 self.__updateObjFromInput()
104 self.dataOut.flagNoData = False
104 self.dataOut.flagNoData = False
105 self.firstdatatime = None
105 self.firstdatatime = None
106
106
107 #---------------------- Correlation Data ---------------------------
107 #---------------------- Correlation Data ---------------------------
108
108
109 if self.dataIn.type == "Correlation":
109 if self.dataIn.type == "Correlation":
110 lagRRange = self.dataIn.lagR
110 lagRRange = self.dataIn.lagR
111 indR = numpy.where(lagRRange == 0)[0][0]
111 indR = numpy.where(lagRRange == 0)[0][0]
112
112
113 self.dataOut.data_pre = self.dataIn.data_corr.copy()[:,:,indR,:]
113 self.dataOut.data_pre = self.dataIn.data_corr.copy()[:,:,indR,:]
114 self.dataOut.abscissaRange = self.dataIn.getLagTRange(1)
114 self.dataOut.abscissaRange = self.dataIn.getLagTRange(1)
115 self.dataOut.noise = self.dataIn.noise
115 self.dataOut.noise = self.dataIn.noise
116 self.dataOut.normFactor = self.dataIn.normFactor
116 self.dataOut.normFactor = self.dataIn.normFactor
117 self.dataOut.SNR = self.dataIn.SNR
117 self.dataOut.SNR = self.dataIn.SNR
118 self.dataOut.pairsList = self.dataIn.pairsList
118 self.dataOut.pairsList = self.dataIn.pairsList
119
119
120 self.__updateObjFromInput()
120 self.__updateObjFromInput()
121 self.dataOut.flagNoData = False
121 self.dataOut.flagNoData = False
122 self.firstdatatime = None
122 self.firstdatatime = None
123
123
124 #------------------- Get Moments ----------------------------------
124 #------------------- Get Moments ----------------------------------
125 def GetMoments(self, channelList = None):
125 def GetMoments(self, channelList = None):
126 '''
126 '''
127 Function GetMoments()
127 Function GetMoments()
128
128
129 Input:
129 Input:
130 channelList : simple channel list to select e.g. [2,3,7]
130 channelList : simple channel list to select e.g. [2,3,7]
131 self.dataOut.data_pre
131 self.dataOut.data_pre
132 self.dataOut.abscissaRange
132 self.dataOut.abscissaRange
133 self.dataOut.noise
133 self.dataOut.noise
134
134
135 Affected:
135 Affected:
136 self.dataOut.data_param
136 self.dataOut.data_param
137 self.dataOut.SNR
137 self.dataOut.SNR
138
138
139 '''
139 '''
140 data = self.dataOut.data_pre
140 data = self.dataOut.data_pre
141 absc = self.dataOut.abscissaRange[:-1]
141 absc = self.dataOut.abscissaRange[:-1]
142 noise = self.dataOut.noise
142 noise = self.dataOut.noise
143
143
144 data_param = numpy.zeros((data.shape[0], 4, data.shape[2]))
144 data_param = numpy.zeros((data.shape[0], 4, data.shape[2]))
145
145
146 if channelList== None: channelList = self.dataOut.channelList
146 if channelList== None: channelList = self.dataOut.channelList
147
147
148 for ind in channelList:
148 for ind in channelList:
149 data_param[ind,:,:] = self.__calculateMoments(data[ind,:,:], absc, noise[ind])
149 data_param[ind,:,:] = self.__calculateMoments(data[ind,:,:], absc, noise[ind])
150
150
151 self.dataOut.data_param = data_param[:,1:]
151 self.dataOut.data_param = data_param[:,1:]
152 self.dataOut.SNR = data_param[:,0]
152 self.dataOut.SNR = data_param[:,0]
153 return
153 return
154
154
155 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):
155 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):
156
156
157 if (nicoh == None): nicoh = 1
157 if (nicoh == None): nicoh = 1
158 if (graph == None): graph = 0
158 if (graph == None): graph = 0
159 if (smooth == None): smooth = 0
159 if (smooth == None): smooth = 0
160 elif (self.smooth < 3): smooth = 0
160 elif (self.smooth < 3): smooth = 0
161
161
162 if (type1 == None): type1 = 0
162 if (type1 == None): type1 = 0
163 if (fwindow == None): fwindow = numpy.zeros(oldfreq.size) + 1
163 if (fwindow == None): fwindow = numpy.zeros(oldfreq.size) + 1
164 if (snrth == None): snrth = -3
164 if (snrth == None): snrth = -3
165 if (dc == None): dc = 0
165 if (dc == None): dc = 0
166 if (aliasing == None): aliasing = 0
166 if (aliasing == None): aliasing = 0
167 if (oldfd == None): oldfd = 0
167 if (oldfd == None): oldfd = 0
168 if (wwauto == None): wwauto = 0
168 if (wwauto == None): wwauto = 0
169
169
170 if (n0 < 1.e-20): n0 = 1.e-20
170 if (n0 < 1.e-20): n0 = 1.e-20
171
171
172 freq = oldfreq
172 freq = oldfreq
173 vec_power = numpy.zeros(oldspec.shape[1])
173 vec_power = numpy.zeros(oldspec.shape[1])
174 vec_fd = numpy.zeros(oldspec.shape[1])
174 vec_fd = numpy.zeros(oldspec.shape[1])
175 vec_w = numpy.zeros(oldspec.shape[1])
175 vec_w = numpy.zeros(oldspec.shape[1])
176 vec_snr = numpy.zeros(oldspec.shape[1])
176 vec_snr = numpy.zeros(oldspec.shape[1])
177
177
178 for ind in range(oldspec.shape[1]):
178 for ind in range(oldspec.shape[1]):
179
179
180 spec = oldspec[:,ind]
180 spec = oldspec[:,ind]
181 aux = spec*fwindow
181 aux = spec*fwindow
182 max_spec = aux.max()
182 max_spec = aux.max()
183 m = list(aux).index(max_spec)
183 m = list(aux).index(max_spec)
184
184
185 #Smooth
185 #Smooth
186 if (smooth == 0): spec2 = spec
186 if (smooth == 0): spec2 = spec
187 else: spec2 = scipy.ndimage.filters.uniform_filter1d(spec,size=smooth)
187 else: spec2 = scipy.ndimage.filters.uniform_filter1d(spec,size=smooth)
188
188
189 # Calculo de Momentos
189 # Calculo de Momentos
190 bb = spec2[range(m,spec2.size)]
190 bb = spec2[range(m,spec2.size)]
191 bb = (bb<n0).nonzero()
191 bb = (bb<n0).nonzero()
192 bb = bb[0]
192 bb = bb[0]
193
193
194 ss = spec2[range(0,m + 1)]
194 ss = spec2[range(0,m + 1)]
195 ss = (ss<n0).nonzero()
195 ss = (ss<n0).nonzero()
196 ss = ss[0]
196 ss = ss[0]
197
197
198 if (bb.size == 0):
198 if (bb.size == 0):
199 bb0 = spec.size - 1 - m
199 bb0 = spec.size - 1 - m
200 else:
200 else:
201 bb0 = bb[0] - 1
201 bb0 = bb[0] - 1
202 if (bb0 < 0):
202 if (bb0 < 0):
203 bb0 = 0
203 bb0 = 0
204
204
205 if (ss.size == 0): ss1 = 1
205 if (ss.size == 0): ss1 = 1
206 else: ss1 = max(ss) + 1
206 else: ss1 = max(ss) + 1
207
207
208 if (ss1 > m): ss1 = m
208 if (ss1 > m): ss1 = m
209
209
210 valid = numpy.asarray(range(int(m + bb0 - ss1 + 1))) + ss1
210 valid = numpy.asarray(range(int(m + bb0 - ss1 + 1))) + ss1
211 power = ((spec2[valid] - n0)*fwindow[valid]).sum()
211 power = ((spec2[valid] - n0)*fwindow[valid]).sum()
212 fd = ((spec2[valid]- n0)*freq[valid]*fwindow[valid]).sum()/power
212 fd = ((spec2[valid]- n0)*freq[valid]*fwindow[valid]).sum()/power
213 w = math.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum()/power)
213 w = math.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum()/power)
214 snr = (spec2.mean()-n0)/n0
214 snr = (spec2.mean()-n0)/n0
215
215
216 if (snr < 1.e-20) :
216 if (snr < 1.e-20) :
217 snr = 1.e-20
217 snr = 1.e-20
218
218
219 vec_power[ind] = power
219 vec_power[ind] = power
220 vec_fd[ind] = fd
220 vec_fd[ind] = fd
221 vec_w[ind] = w
221 vec_w[ind] = w
222 vec_snr[ind] = snr
222 vec_snr[ind] = snr
223
223
224 moments = numpy.vstack((vec_snr, vec_power, vec_fd, vec_w))
224 moments = numpy.vstack((vec_snr, vec_power, vec_fd, vec_w))
225 return moments
225 return moments
226
226
227 #------------------- Get Lags ----------------------------------
227 #------------------- Get Lags ----------------------------------
228
228
229 def GetLags(self):
229 def GetLags(self):
230 '''
230 '''
231 Function GetMoments()
231 Function GetMoments()
232
232
233 Input:
233 Input:
234 self.dataOut.data_pre
234 self.dataOut.data_pre
235 self.dataOut.abscissaRange
235 self.dataOut.abscissaRange
236 self.dataOut.noise
236 self.dataOut.noise
237 self.dataOut.normFactor
237 self.dataOut.normFactor
238 self.dataOut.SNR
238 self.dataOut.SNR
239 self.dataOut.pairsList
239 self.dataOut.pairsList
240 self.dataOut.nChannels
240 self.dataOut.nChannels
241
241
242 Affected:
242 Affected:
243 self.dataOut.data_param
243 self.dataOut.data_param
244
244
245 '''
245 '''
246 data = self.dataOut.data_pre
246 data = self.dataOut.data_pre
247 normFactor = self.dataOut.normFactor
247 normFactor = self.dataOut.normFactor
248 nHeights = self.dataOut.nHeights
248 nHeights = self.dataOut.nHeights
249 absc = self.dataOut.abscissaRange[:-1]
249 absc = self.dataOut.abscissaRange[:-1]
250 noise = self.dataOut.noise
250 noise = self.dataOut.noise
251 SNR = self.dataOut.SNR
251 SNR = self.dataOut.SNR
252 pairsList = self.dataOut.pairsList
252 pairsList = self.dataOut.pairsList
253 nChannels = self.dataOut.nChannels
253 nChannels = self.dataOut.nChannels
254 pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)
254 pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)
255 self.dataOut.data_param = numpy.zeros((len(pairsCrossCorr)*2 + 1, nHeights))
255 self.dataOut.data_param = numpy.zeros((len(pairsCrossCorr)*2 + 1, nHeights))
256
256
257 dataNorm = numpy.abs(data)
257 dataNorm = numpy.abs(data)
258 for l in range(len(pairsList)):
258 for l in range(len(pairsList)):
259 dataNorm[l,:,:] = dataNorm[l,:,:]/normFactor[l,:]
259 dataNorm[l,:,:] = dataNorm[l,:,:]/normFactor[l,:]
260
260
261 self.dataOut.data_param[:-1,:] = self.__calculateTaus(dataNorm, pairsCrossCorr, pairsAutoCorr, absc)
261 self.dataOut.data_param[:-1,:] = self.__calculateTaus(dataNorm, pairsCrossCorr, pairsAutoCorr, absc)
262 self.dataOut.data_param[-1,:] = self.__calculateLag1Phase(data, pairsAutoCorr, absc)
262 self.dataOut.data_param[-1,:] = self.__calculateLag1Phase(data, pairsAutoCorr, absc)
263 return
263 return
264
264
265 def __getPairsAutoCorr(self, pairsList, nChannels):
265 def __getPairsAutoCorr(self, pairsList, nChannels):
266
266
267 pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
267 pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
268
268
269 for l in range(len(pairsList)):
269 for l in range(len(pairsList)):
270 firstChannel = pairsList[l][0]
270 firstChannel = pairsList[l][0]
271 secondChannel = pairsList[l][1]
271 secondChannel = pairsList[l][1]
272
272
273 #Obteniendo pares de Autocorrelacion
273 #Obteniendo pares de Autocorrelacion
274 if firstChannel == secondChannel:
274 if firstChannel == secondChannel:
275 pairsAutoCorr[firstChannel] = int(l)
275 pairsAutoCorr[firstChannel] = int(l)
276
276
277 pairsAutoCorr = pairsAutoCorr.astype(int)
277 pairsAutoCorr = pairsAutoCorr.astype(int)
278
278
279 pairsCrossCorr = range(len(pairsList))
279 pairsCrossCorr = range(len(pairsList))
280 pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
280 pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
281
281
282 return pairsAutoCorr, pairsCrossCorr
282 return pairsAutoCorr, pairsCrossCorr
283
283
284 def __calculateTaus(self, data, pairsCrossCorr, pairsAutoCorr, lagTRange):
284 def __calculateTaus(self, data, pairsCrossCorr, pairsAutoCorr, lagTRange):
285
285
286 Pt0 = data.shape[1]/2
286 Pt0 = data.shape[1]/2
287 #Funcion de Autocorrelacion
287 #Funcion de Autocorrelacion
288 dataAutoCorr = stats.nanmean(data[pairsAutoCorr,:,:], axis = 0)
288 dataAutoCorr = stats.nanmean(data[pairsAutoCorr,:,:], axis = 0)
289
289
290 #Obtencion Indice de TauCross
290 #Obtencion Indice de TauCross
291 indCross = data[pairsCrossCorr,:,:].argmax(axis = 1)
291 indCross = data[pairsCrossCorr,:,:].argmax(axis = 1)
292 #Obtencion Indice de TauAuto
292 #Obtencion Indice de TauAuto
293 indAuto = numpy.zeros(indCross.shape,dtype = 'int')
293 indAuto = numpy.zeros(indCross.shape,dtype = 'int')
294 CCValue = data[pairsCrossCorr,Pt0,:]
294 CCValue = data[pairsCrossCorr,Pt0,:]
295 for i in range(pairsCrossCorr.size):
295 for i in range(pairsCrossCorr.size):
296 indAuto[i,:] = numpy.abs(dataAutoCorr - CCValue[i,:]).argmin(axis = 0)
296 indAuto[i,:] = numpy.abs(dataAutoCorr - CCValue[i,:]).argmin(axis = 0)
297
297
298 #Obtencion de TauCross y TauAuto
298 #Obtencion de TauCross y TauAuto
299 tauCross = lagTRange[indCross]
299 tauCross = lagTRange[indCross]
300 tauAuto = lagTRange[indAuto]
300 tauAuto = lagTRange[indAuto]
301
301
302 Nan1, Nan2 = numpy.where(tauCross == lagTRange[0])
302 Nan1, Nan2 = numpy.where(tauCross == lagTRange[0])
303
303
304 tauCross[Nan1,Nan2] = numpy.nan
304 tauCross[Nan1,Nan2] = numpy.nan
305 tauAuto[Nan1,Nan2] = numpy.nan
305 tauAuto[Nan1,Nan2] = numpy.nan
306 tau = numpy.vstack((tauCross,tauAuto))
306 tau = numpy.vstack((tauCross,tauAuto))
307
307
308 return tau
308 return tau
309
309
310 def __calculateLag1Phase(self, data, pairs, lagTRange):
310 def __calculateLag1Phase(self, data, pairs, lagTRange):
311 data1 = stats.nanmean(data[pairs,:,:], axis = 0)
311 data1 = stats.nanmean(data[pairs,:,:], axis = 0)
312 lag1 = numpy.where(lagTRange == 0)[0][0] + 1
312 lag1 = numpy.where(lagTRange == 0)[0][0] + 1
313
313
314 phase = numpy.angle(data1[lag1,:])
314 phase = numpy.angle(data1[lag1,:])
315
315
316 return phase
316 return phase
317 #------------------- Detect Meteors ------------------------------
317 #------------------- Detect Meteors ------------------------------
318
318
319 def DetectMeteors(self, hei_ref = None, tauindex = 0,
319 def DetectMeteors(self, hei_ref = None, tauindex = 0,
320 predefinedPhaseShifts = None, centerReceiverIndex = 2,
320 predefinedPhaseShifts = None, centerReceiverIndex = 2,
321 cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25,
321 cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25,
322 noise_timeStep = 4, noise_multiple = 4,
322 noise_timeStep = 4, noise_multiple = 4,
323 multDet_timeLimit = 1, multDet_rangeLimit = 3,
323 multDet_timeLimit = 1, multDet_rangeLimit = 3,
324 phaseThresh = 20, SNRThresh = 8,
324 phaseThresh = 20, SNRThresh = 8,
325 hmin = 70, hmax=110, azimuth = 0) :
325 hmin = 70, hmax=110, azimuth = 0) :
326
326
327 '''
327 '''
328 Function DetectMeteors()
328 Function DetectMeteors()
329 Project developed with paper:
329 Project developed with paper:
330 HOLDSWORTH ET AL. 2004
330 HOLDSWORTH ET AL. 2004
331
331
332 Input:
332 Input:
333 self.dataOut.data_pre
333 self.dataOut.data_pre
334
334
335 centerReceiverIndex: From the channels, which is the center receiver
335 centerReceiverIndex: From the channels, which is the center receiver
336
336
337 hei_ref: Height reference for the Beacon signal extraction
337 hei_ref: Height reference for the Beacon signal extraction
338 tauindex:
338 tauindex:
339 predefinedPhaseShifts: Predefined phase offset for the voltge signals
339 predefinedPhaseShifts: Predefined phase offset for the voltge signals
340
340
341 cohDetection: Whether to user Coherent detection or not
341 cohDetection: Whether to user Coherent detection or not
342 cohDet_timeStep: Coherent Detection calculation time step
342 cohDet_timeStep: Coherent Detection calculation time step
343 cohDet_thresh: Coherent Detection phase threshold to correct phases
343 cohDet_thresh: Coherent Detection phase threshold to correct phases
344
344
345 noise_timeStep: Noise calculation time step
345 noise_timeStep: Noise calculation time step
346 noise_multiple: Noise multiple to define signal threshold
346 noise_multiple: Noise multiple to define signal threshold
347
347
348 multDet_timeLimit: Multiple Detection Removal time limit in seconds
348 multDet_timeLimit: Multiple Detection Removal time limit in seconds
349 multDet_rangeLimit: Multiple Detection Removal range limit in km
349 multDet_rangeLimit: Multiple Detection Removal range limit in km
350
350
351 phaseThresh: Maximum phase difference between receiver to be consider a meteor
351 phaseThresh: Maximum phase difference between receiver to be consider a meteor
352 SNRThresh: Minimum SNR threshold of the meteor signal to be consider a meteor
352 SNRThresh: Minimum SNR threshold of the meteor signal to be consider a meteor
353
353
354 hmin: Minimum Height of the meteor to use it in the further wind estimations
354 hmin: Minimum Height of the meteor to use it in the further wind estimations
355 hmax: Maximum Height of the meteor to use it in the further wind estimations
355 hmax: Maximum Height of the meteor to use it in the further wind estimations
356 azimuth: Azimuth angle correction
356 azimuth: Azimuth angle correction
357
357
358 Affected:
358 Affected:
359 self.dataOut.data_param
359 self.dataOut.data_param
360
360
361 Rejection Criteria (Errors):
361 Rejection Criteria (Errors):
362 0: No error; analysis OK
362 0: No error; analysis OK
363 1: SNR < SNR threshold
363 1: SNR < SNR threshold
364 2: angle of arrival (AOA) ambiguously determined
364 2: angle of arrival (AOA) ambiguously determined
365 3: AOA estimate not feasible
365 3: AOA estimate not feasible
366 4: Large difference in AOAs obtained from different antenna baselines
366 4: Large difference in AOAs obtained from different antenna baselines
367 5: echo at start or end of time series
367 5: echo at start or end of time series
368 6: echo less than 5 examples long; too short for analysis
368 6: echo less than 5 examples long; too short for analysis
369 7: echo rise exceeds 0.3s
369 7: echo rise exceeds 0.3s
370 8: echo decay time less than twice rise time
370 8: echo decay time less than twice rise time
371 9: large power level before echo
371 9: large power level before echo
372 10: large power level after echo
372 10: large power level after echo
373 11: poor fit to amplitude for estimation of decay time
373 11: poor fit to amplitude for estimation of decay time
374 12: poor fit to CCF phase variation for estimation of radial drift velocity
374 12: poor fit to CCF phase variation for estimation of radial drift velocity
375 13: height unresolvable echo: not valid height within 70 to 110 km
375 13: height unresolvable echo: not valid height within 70 to 110 km
376 14: height ambiguous echo: more then one possible height within 70 to 110 km
376 14: height ambiguous echo: more then one possible height within 70 to 110 km
377 15: radial drift velocity or projected horizontal velocity exceeds 200 m/s
377 15: radial drift velocity or projected horizontal velocity exceeds 200 m/s
378 16: oscilatory echo, indicating event most likely not an underdense echo
378 16: oscilatory echo, indicating event most likely not an underdense echo
379
379
380 17: phase difference in meteor Reestimation
380 17: phase difference in meteor Reestimation
381
381
382 Data Storage:
382 Data Storage:
383 Meteors for Wind Estimation (8):
383 Meteors for Wind Estimation (8):
384 Day Hour | Range Height
384 Day Hour | Range Height
385 Azimuth Zenith errorCosDir
385 Azimuth Zenith errorCosDir
386 VelRad errorVelRad
386 VelRad errorVelRad
387 TypeError
387 TypeError
388
388
389 '''
389 '''
390 #Get Beacon signal
390 #Get Beacon signal
391 newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
391 newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
392
392
393 if hei_ref != None:
393 if hei_ref != None:
394 newheis = numpy.where(self.dataOut.heightList>hei_ref)
394 newheis = numpy.where(self.dataOut.heightList>hei_ref)
395
395
396 heiRang = self.dataOut.getHeiRange()
396 heiRang = self.dataOut.getHeiRange()
397 #Pairs List
397 #Pairs List
398 pairslist = []
398 pairslist = []
399 nChannel = self.dataOut.nChannels
399 nChannel = self.dataOut.nChannels
400 for i in range(nChannel):
400 for i in range(nChannel):
401 if i != centerReceiverIndex:
401 if i != centerReceiverIndex:
402 pairslist.append((centerReceiverIndex,i))
402 pairslist.append((centerReceiverIndex,i))
403
403
404 #****************REMOVING HARDWARE PHASE DIFFERENCES***************
404 #****************REMOVING HARDWARE PHASE DIFFERENCES***************
405 # see if the user put in pre defined phase shifts
405 # see if the user put in pre defined phase shifts
406 voltsPShift = self.dataOut.data_pre.copy()
406 voltsPShift = self.dataOut.data_pre.copy()
407
407
408 if predefinedPhaseShifts != None:
408 if predefinedPhaseShifts != None:
409 hardwarePhaseShifts = numpy.array(predefinedPhaseShifts)*numpy.pi/180
409 hardwarePhaseShifts = numpy.array(predefinedPhaseShifts)*numpy.pi/180
410 else:
410 else:
411 #get hardware phase shifts using beacon signal
411 #get hardware phase shifts using beacon signal
412 hardwarePhaseShifts = self.__getHardwarePhaseDiff(self.dataOut.data_pre, pairslist, newheis, 10)
412 hardwarePhaseShifts = self.__getHardwarePhaseDiff(self.dataOut.data_pre, pairslist, newheis, 10)
413 hardwarePhaseShifts = numpy.insert(hardwarePhaseShifts,centerReceiverIndex,0)
413 hardwarePhaseShifts = numpy.insert(hardwarePhaseShifts,centerReceiverIndex,0)
414
414
415 voltsPShift = numpy.zeros((self.dataOut.data_pre.shape[0],self.dataOut.data_pre.shape[1],self.dataOut.data_pre.shape[2]), dtype = 'complex')
415 voltsPShift = numpy.zeros((self.dataOut.data_pre.shape[0],self.dataOut.data_pre.shape[1],self.dataOut.data_pre.shape[2]), dtype = 'complex')
416 for i in range(self.dataOut.data_pre.shape[0]):
416 for i in range(self.dataOut.data_pre.shape[0]):
417 voltsPShift[i,:,:] = self.__shiftPhase(self.dataOut.data_pre[i,:,:], hardwarePhaseShifts[i])
417 voltsPShift[i,:,:] = self.__shiftPhase(self.dataOut.data_pre[i,:,:], hardwarePhaseShifts[i])
418 #******************END OF REMOVING HARDWARE PHASE DIFFERENCES*********
418 #******************END OF REMOVING HARDWARE PHASE DIFFERENCES*********
419
419
420 #Remove DC
420 #Remove DC
421 voltsDC = numpy.mean(voltsPShift,1)
421 voltsDC = numpy.mean(voltsPShift,1)
422 voltsDC = numpy.mean(voltsDC,1)
422 voltsDC = numpy.mean(voltsDC,1)
423 for i in range(voltsDC.shape[0]):
423 for i in range(voltsDC.shape[0]):
424 voltsPShift[i] = voltsPShift[i] - voltsDC[i]
424 voltsPShift[i] = voltsPShift[i] - voltsDC[i]
425
425
426 #Don't considerate last heights, theyre used to calculate Hardware Phase Shift
426 #Don't considerate last heights, theyre used to calculate Hardware Phase Shift
427 voltsPShift = voltsPShift[:,:,:newheis[0][0]]
427 voltsPShift = voltsPShift[:,:,:newheis[0][0]]
428
428
429 #************ FIND POWER OF DATA W/COH OR NON COH DETECTION (3.4) **********
429 #************ FIND POWER OF DATA W/COH OR NON COH DETECTION (3.4) **********
430 #Coherent Detection
430 #Coherent Detection
431 if cohDetection:
431 if cohDetection:
432 #use coherent detection to get the net power
432 #use coherent detection to get the net power
433 cohDet_thresh = cohDet_thresh*numpy.pi/180
433 cohDet_thresh = cohDet_thresh*numpy.pi/180
434 voltsPShift = self.__coherentDetection(voltsPShift, cohDet_timeStep, self.dataOut.timeInterval, pairslist, cohDet_thresh)
434 voltsPShift = self.__coherentDetection(voltsPShift, cohDet_timeStep, self.dataOut.timeInterval, pairslist, cohDet_thresh)
435
435
436 #Non-coherent detection!
436 #Non-coherent detection!
437 powerNet = numpy.nansum(numpy.abs(voltsPShift[:,:,:])**2,0)
437 powerNet = numpy.nansum(numpy.abs(voltsPShift[:,:,:])**2,0)
438 #********** END OF COH/NON-COH POWER CALCULATION**********************
438 #********** END OF COH/NON-COH POWER CALCULATION**********************
439
439
440 #********** FIND THE NOISE LEVEL AND POSSIBLE METEORS ****************
440 #********** FIND THE NOISE LEVEL AND POSSIBLE METEORS ****************
441 #Get noise
441 #Get noise
442 noise, noise1 = self.__getNoise(powerNet, noise_timeStep, self.dataOut.timeInterval)
442 noise, noise1 = self.__getNoise(powerNet, noise_timeStep, self.dataOut.timeInterval)
443 # noise = self.getNoise1(powerNet, noise_timeStep, self.dataOut.timeInterval)
443 # noise = self.getNoise1(powerNet, noise_timeStep, self.dataOut.timeInterval)
444 #Get signal threshold
444 #Get signal threshold
445 signalThresh = noise_multiple*noise
445 signalThresh = noise_multiple*noise
446 #Meteor echoes detection
446 #Meteor echoes detection
447 listMeteors = self.__findMeteors(powerNet, signalThresh)
447 listMeteors = self.__findMeteors(powerNet, signalThresh)
448 #******* END OF NOISE LEVEL AND POSSIBLE METEORS CACULATION **********
448 #******* END OF NOISE LEVEL AND POSSIBLE METEORS CACULATION **********
449
449
450 #************** REMOVE MULTIPLE DETECTIONS (3.5) ***************************
450 #************** REMOVE MULTIPLE DETECTIONS (3.5) ***************************
451 #Parameters
451 #Parameters
452 heiRange = self.dataOut.getHeiRange()
452 heiRange = self.dataOut.getHeiRange()
453 rangeInterval = heiRange[1] - heiRange[0]
453 rangeInterval = heiRange[1] - heiRange[0]
454 rangeLimit = multDet_rangeLimit/rangeInterval
454 rangeLimit = multDet_rangeLimit/rangeInterval
455 timeLimit = multDet_timeLimit/self.dataOut.timeInterval
455 timeLimit = multDet_timeLimit/self.dataOut.timeInterval
456 #Multiple detection removals
456 #Multiple detection removals
457 listMeteors1 = self.__removeMultipleDetections(listMeteors, rangeLimit, timeLimit)
457 listMeteors1 = self.__removeMultipleDetections(listMeteors, rangeLimit, timeLimit)
458 #************ END OF REMOVE MULTIPLE DETECTIONS **********************
458 #************ END OF REMOVE MULTIPLE DETECTIONS **********************
459
459
460 #********************* METEOR REESTIMATION (3.7, 3.8, 3.9, 3.10) ********************
460 #********************* METEOR REESTIMATION (3.7, 3.8, 3.9, 3.10) ********************
461 #Parameters
461 #Parameters
462 phaseThresh = phaseThresh*numpy.pi/180
462 phaseThresh = phaseThresh*numpy.pi/180
463 thresh = [phaseThresh, noise_multiple, SNRThresh]
463 thresh = [phaseThresh, noise_multiple, SNRThresh]
464 #Meteor reestimation (Errors N 1, 6, 12, 17)
464 #Meteor reestimation (Errors N 1, 6, 12, 17)
465 listMeteors2, listMeteorsPower, listMeteorsVolts = self.__meteorReestimation(listMeteors1, voltsPShift, pairslist, thresh, noise, self.dataOut.timeInterval, self.dataOut.frequency)
465 listMeteors2, listMeteorsPower, listMeteorsVolts = self.__meteorReestimation(listMeteors1, voltsPShift, pairslist, thresh, noise, self.dataOut.timeInterval, self.dataOut.frequency)
466 # listMeteors2, listMeteorsPower, listMeteorsVolts = self.meteorReestimation3(listMeteors2, listMeteorsPower, listMeteorsVolts, voltsPShift, pairslist, thresh, noise)
466 # listMeteors2, listMeteorsPower, listMeteorsVolts = self.meteorReestimation3(listMeteors2, listMeteorsPower, listMeteorsVolts, voltsPShift, pairslist, thresh, noise)
467 #Estimation of decay times (Errors N 7, 8, 11)
467 #Estimation of decay times (Errors N 7, 8, 11)
468 listMeteors3 = self.__estimateDecayTime(listMeteors2, listMeteorsPower, self.dataOut.timeInterval, self.dataOut.frequency)
468 listMeteors3 = self.__estimateDecayTime(listMeteors2, listMeteorsPower, self.dataOut.timeInterval, self.dataOut.frequency)
469 #******************* END OF METEOR REESTIMATION *******************
469 #******************* END OF METEOR REESTIMATION *******************
470
470
471 #********************* METEOR PARAMETERS CALCULATION (3.11, 3.12, 3.13) **************************
471 #********************* METEOR PARAMETERS CALCULATION (3.11, 3.12, 3.13) **************************
472 #Calculating Radial Velocity (Error N 15)
472 #Calculating Radial Velocity (Error N 15)
473 radialStdThresh = 10
473 radialStdThresh = 10
474 listMeteors4 = self.__getRadialVelocity(listMeteors3, listMeteorsVolts, radialStdThresh, pairslist, self.dataOut.timeInterval)
474 listMeteors4 = self.__getRadialVelocity(listMeteors3, listMeteorsVolts, radialStdThresh, pairslist, self.dataOut.timeInterval)
475
475
476 if len(listMeteors4) > 0:
476 if len(listMeteors4) > 0:
477 #Setting New Array
477 #Setting New Array
478 date = repr(self.dataOut.datatime)
478 date = repr(self.dataOut.datatime)
479 arrayMeteors4, arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang)
479 arrayMeteors4, arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang)
480
480
481 #Calculate AOA (Error N 3, 4)
481 #Calculate AOA (Error N 3, 4)
482 #JONES ET AL. 1998
482 #JONES ET AL. 1998
483 AOAthresh = numpy.pi/8
483 AOAthresh = numpy.pi/8
484 error = arrayParameters[:,-1]
484 error = arrayParameters[:,-1]
485 phases = -arrayMeteors4[:,9:13]
485 phases = -arrayMeteors4[:,9:13]
486 pairsList = []
486 pairsList = []
487 pairsList.append((0,3))
487 pairsList.append((0,3))
488 pairsList.append((1,2))
488 pairsList.append((1,2))
489 arrayParameters[:,4:7], arrayParameters[:,-1] = self.__getAOA(phases, pairsList, error, AOAthresh, azimuth)
489 arrayParameters[:,4:7], arrayParameters[:,-1] = self.__getAOA(phases, pairsList, error, AOAthresh, azimuth)
490
490
491 #Calculate Heights (Error N 13 and 14)
491 #Calculate Heights (Error N 13 and 14)
492 error = arrayParameters[:,-1]
492 error = arrayParameters[:,-1]
493 Ranges = arrayParameters[:,2]
493 Ranges = arrayParameters[:,2]
494 zenith = arrayParameters[:,5]
494 zenith = arrayParameters[:,5]
495 arrayParameters[:,3], arrayParameters[:,-1] = self.__getHeights(Ranges, zenith, error, hmin, hmax)
495 arrayParameters[:,3], arrayParameters[:,-1] = self.__getHeights(Ranges, zenith, error, hmin, hmax)
496 #********************* END OF PARAMETERS CALCULATION **************************
496 #********************* END OF PARAMETERS CALCULATION **************************
497
497
498 #***************************+ SAVE DATA IN HDF5 FORMAT **********************
498 #***************************+ SAVE DATA IN HDF5 FORMAT **********************
499 self.dataOut.data_param = arrayParameters
499 self.dataOut.data_param = arrayParameters
500
500
501 return
501 return
502
502
503 def __getHardwarePhaseDiff(self, voltage0, pairslist, newheis, n):
503 def __getHardwarePhaseDiff(self, voltage0, pairslist, newheis, n):
504
504
505 minIndex = min(newheis[0])
505 minIndex = min(newheis[0])
506 maxIndex = max(newheis[0])
506 maxIndex = max(newheis[0])
507
507
508 voltage = voltage0[:,:,minIndex:maxIndex+1]
508 voltage = voltage0[:,:,minIndex:maxIndex+1]
509 nLength = voltage.shape[1]/n
509 nLength = voltage.shape[1]/n
510 nMin = 0
510 nMin = 0
511 nMax = 0
511 nMax = 0
512 phaseOffset = numpy.zeros((len(pairslist),n))
512 phaseOffset = numpy.zeros((len(pairslist),n))
513
513
514 for i in range(n):
514 for i in range(n):
515 nMax += nLength
515 nMax += nLength
516 phaseCCF = -numpy.angle(self.__calculateCCF(voltage[:,nMin:nMax,:], pairslist, [0]))
516 phaseCCF = -numpy.angle(self.__calculateCCF(voltage[:,nMin:nMax,:], pairslist, [0]))
517 phaseCCF = numpy.mean(phaseCCF, axis = 2)
517 phaseCCF = numpy.mean(phaseCCF, axis = 2)
518 phaseOffset[:,i] = phaseCCF.transpose()
518 phaseOffset[:,i] = phaseCCF.transpose()
519 nMin = nMax
519 nMin = nMax
520 # phaseDiff, phaseArrival = self.estimatePhaseDifference(voltage, pairslist)
520 # phaseDiff, phaseArrival = self.estimatePhaseDifference(voltage, pairslist)
521
521
522 #Remove Outliers
522 #Remove Outliers
523 factor = 2
523 factor = 2
524 wt = phaseOffset - signal.medfilt(phaseOffset,(1,5))
524 wt = phaseOffset - signal.medfilt(phaseOffset,(1,5))
525 dw = numpy.std(wt,axis = 1)
525 dw = numpy.std(wt,axis = 1)
526 dw = dw.reshape((dw.size,1))
526 dw = dw.reshape((dw.size,1))
527 ind = numpy.where(numpy.logical_or(wt>dw*factor,wt<-dw*factor))
527 ind = numpy.where(numpy.logical_or(wt>dw*factor,wt<-dw*factor))
528 phaseOffset[ind] = numpy.nan
528 phaseOffset[ind] = numpy.nan
529 phaseOffset = stats.nanmean(phaseOffset, axis=1)
529 phaseOffset = stats.nanmean(phaseOffset, axis=1)
530
530
531 return phaseOffset
531 return phaseOffset
532
532
533 def __shiftPhase(self, data, phaseShift):
533 def __shiftPhase(self, data, phaseShift):
534 #this will shift the phase of a complex number
534 #this will shift the phase of a complex number
535 dataShifted = numpy.abs(data) * numpy.exp((numpy.angle(data)+phaseShift)*1j)
535 dataShifted = numpy.abs(data) * numpy.exp((numpy.angle(data)+phaseShift)*1j)
536 return dataShifted
536 return dataShifted
537
537
538 def __estimatePhaseDifference(self, array, pairslist):
538 def __estimatePhaseDifference(self, array, pairslist):
539 nChannel = array.shape[0]
539 nChannel = array.shape[0]
540 nHeights = array.shape[2]
540 nHeights = array.shape[2]
541 numPairs = len(pairslist)
541 numPairs = len(pairslist)
542 # phaseCCF = numpy.zeros((nChannel, 5, nHeights))
542 # phaseCCF = numpy.zeros((nChannel, 5, nHeights))
543 phaseCCF = numpy.angle(self.__calculateCCF(array, pairslist, [-2,-1,0,1,2]))
543 phaseCCF = numpy.angle(self.__calculateCCF(array, pairslist, [-2,-1,0,1,2]))
544
544
545 #Correct phases
545 #Correct phases
546 derPhaseCCF = phaseCCF[:,1:,:] - phaseCCF[:,0:-1,:]
546 derPhaseCCF = phaseCCF[:,1:,:] - phaseCCF[:,0:-1,:]
547 indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
547 indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
548
548
549 if indDer[0].shape[0] > 0:
549 if indDer[0].shape[0] > 0:
550 for i in range(indDer[0].shape[0]):
550 for i in range(indDer[0].shape[0]):
551 signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i],indDer[2][i]])
551 signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i],indDer[2][i]])
552 phaseCCF[indDer[0][i],indDer[1][i]+1:,:] += signo*2*numpy.pi
552 phaseCCF[indDer[0][i],indDer[1][i]+1:,:] += signo*2*numpy.pi
553
553
554 # for j in range(numSides):
554 # for j in range(numSides):
555 # phaseCCFAux = self.calculateCCF(arrayCenter, arraySides[j,:,:], [-2,1,0,1,2])
555 # phaseCCFAux = self.calculateCCF(arrayCenter, arraySides[j,:,:], [-2,1,0,1,2])
556 # phaseCCF[j,:,:] = numpy.angle(phaseCCFAux)
556 # phaseCCF[j,:,:] = numpy.angle(phaseCCFAux)
557 #
557 #
558 #Linear
558 #Linear
559 phaseInt = numpy.zeros((numPairs,1))
559 phaseInt = numpy.zeros((numPairs,1))
560 angAllCCF = phaseCCF[:,[0,1,3,4],0]
560 angAllCCF = phaseCCF[:,[0,1,3,4],0]
561 for j in range(numPairs):
561 for j in range(numPairs):
562 fit = stats.linregress([-2,-1,1,2],angAllCCF[j,:])
562 fit = stats.linregress([-2,-1,1,2],angAllCCF[j,:])
563 phaseInt[j] = fit[1]
563 phaseInt[j] = fit[1]
564 #Phase Differences
564 #Phase Differences
565 phaseDiff = phaseInt - phaseCCF[:,2,:]
565 phaseDiff = phaseInt - phaseCCF[:,2,:]
566 phaseArrival = phaseInt.reshape(phaseInt.size)
566 phaseArrival = phaseInt.reshape(phaseInt.size)
567
567
568 #Dealias
568 #Dealias
569 indAlias = numpy.where(phaseArrival > numpy.pi)
569 indAlias = numpy.where(phaseArrival > numpy.pi)
570 phaseArrival[indAlias] -= 2*numpy.pi
570 phaseArrival[indAlias] -= 2*numpy.pi
571 indAlias = numpy.where(phaseArrival < -numpy.pi)
571 indAlias = numpy.where(phaseArrival < -numpy.pi)
572 phaseArrival[indAlias] += 2*numpy.pi
572 phaseArrival[indAlias] += 2*numpy.pi
573
573
574 return phaseDiff, phaseArrival
574 return phaseDiff, phaseArrival
575
575
576 def __coherentDetection(self, volts, timeSegment, timeInterval, pairslist, thresh):
576 def __coherentDetection(self, volts, timeSegment, timeInterval, pairslist, thresh):
577 #this function will run the coherent detection used in Holdworth et al. 2004 and return the net power
577 #this function will run the coherent detection used in Holdworth et al. 2004 and return the net power
578 #find the phase shifts of each channel over 1 second intervals
578 #find the phase shifts of each channel over 1 second intervals
579 #only look at ranges below the beacon signal
579 #only look at ranges below the beacon signal
580 numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
580 numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
581 numBlocks = int(volts.shape[1]/numProfPerBlock)
581 numBlocks = int(volts.shape[1]/numProfPerBlock)
582 numHeights = volts.shape[2]
582 numHeights = volts.shape[2]
583 nChannel = volts.shape[0]
583 nChannel = volts.shape[0]
584 voltsCohDet = volts.copy()
584 voltsCohDet = volts.copy()
585
585
586 pairsarray = numpy.array(pairslist)
586 pairsarray = numpy.array(pairslist)
587 indSides = pairsarray[:,1]
587 indSides = pairsarray[:,1]
588 # indSides = numpy.array(range(nChannel))
588 # indSides = numpy.array(range(nChannel))
589 # indSides = numpy.delete(indSides, indCenter)
589 # indSides = numpy.delete(indSides, indCenter)
590 #
590 #
591 # listCenter = numpy.array_split(volts[indCenter,:,:], numBlocks, 0)
591 # listCenter = numpy.array_split(volts[indCenter,:,:], numBlocks, 0)
592 listBlocks = numpy.array_split(volts, numBlocks, 1)
592 listBlocks = numpy.array_split(volts, numBlocks, 1)
593
593
594 startInd = 0
594 startInd = 0
595 endInd = 0
595 endInd = 0
596
596
597 for i in range(numBlocks):
597 for i in range(numBlocks):
598 startInd = endInd
598 startInd = endInd
599 endInd = endInd + listBlocks[i].shape[1]
599 endInd = endInd + listBlocks[i].shape[1]
600
600
601 arrayBlock = listBlocks[i]
601 arrayBlock = listBlocks[i]
602 # arrayBlockCenter = listCenter[i]
602 # arrayBlockCenter = listCenter[i]
603
603
604 #Estimate the Phase Difference
604 #Estimate the Phase Difference
605 phaseDiff, aux = self.__estimatePhaseDifference(arrayBlock, pairslist)
605 phaseDiff, aux = self.__estimatePhaseDifference(arrayBlock, pairslist)
606 #Phase Difference RMS
606 #Phase Difference RMS
607 arrayPhaseRMS = numpy.abs(phaseDiff)
607 arrayPhaseRMS = numpy.abs(phaseDiff)
608 phaseRMSaux = numpy.sum(arrayPhaseRMS < thresh,0)
608 phaseRMSaux = numpy.sum(arrayPhaseRMS < thresh,0)
609 indPhase = numpy.where(phaseRMSaux==4)
609 indPhase = numpy.where(phaseRMSaux==4)
610 #Shifting
610 #Shifting
611 if indPhase[0].shape[0] > 0:
611 if indPhase[0].shape[0] > 0:
612 for j in range(indSides.size):
612 for j in range(indSides.size):
613 arrayBlock[indSides[j],:,indPhase] = self.__shiftPhase(arrayBlock[indSides[j],:,indPhase], phaseDiff[j,indPhase].transpose())
613 arrayBlock[indSides[j],:,indPhase] = self.__shiftPhase(arrayBlock[indSides[j],:,indPhase], phaseDiff[j,indPhase].transpose())
614 voltsCohDet[:,startInd:endInd,:] = arrayBlock
614 voltsCohDet[:,startInd:endInd,:] = arrayBlock
615
615
616 return voltsCohDet
616 return voltsCohDet
617
617
618 def __calculateCCF(self, volts, pairslist ,laglist):
618 def __calculateCCF(self, volts, pairslist ,laglist):
619
619
620 nHeights = volts.shape[2]
620 nHeights = volts.shape[2]
621 nPoints = volts.shape[1]
621 nPoints = volts.shape[1]
622 voltsCCF = numpy.zeros((len(pairslist), len(laglist), nHeights),dtype = 'complex')
622 voltsCCF = numpy.zeros((len(pairslist), len(laglist), nHeights),dtype = 'complex')
623
623
624 for i in range(len(pairslist)):
624 for i in range(len(pairslist)):
625 volts1 = volts[pairslist[i][0]]
625 volts1 = volts[pairslist[i][0]]
626 volts2 = volts[pairslist[i][1]]
626 volts2 = volts[pairslist[i][1]]
627
627
628 for t in range(len(laglist)):
628 for t in range(len(laglist)):
629 idxT = laglist[t]
629 idxT = laglist[t]
630 if idxT >= 0:
630 if idxT >= 0:
631 vStacked = numpy.vstack((volts2[idxT:,:],
631 vStacked = numpy.vstack((volts2[idxT:,:],
632 numpy.zeros((idxT, nHeights),dtype='complex')))
632 numpy.zeros((idxT, nHeights),dtype='complex')))
633 else:
633 else:
634 vStacked = numpy.vstack((numpy.zeros((-idxT, nHeights),dtype='complex'),
634 vStacked = numpy.vstack((numpy.zeros((-idxT, nHeights),dtype='complex'),
635 volts2[:(nPoints + idxT),:]))
635 volts2[:(nPoints + idxT),:]))
636 voltsCCF[i,t,:] = numpy.sum((numpy.conjugate(volts1)*vStacked),axis=0)
636 voltsCCF[i,t,:] = numpy.sum((numpy.conjugate(volts1)*vStacked),axis=0)
637
637
638 vStacked = None
638 vStacked = None
639 return voltsCCF
639 return voltsCCF
640
640
641 def __getNoise(self, power, timeSegment, timeInterval):
641 def __getNoise(self, power, timeSegment, timeInterval):
642 numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
642 numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
643 numBlocks = int(power.shape[0]/numProfPerBlock)
643 numBlocks = int(power.shape[0]/numProfPerBlock)
644 numHeights = power.shape[1]
644 numHeights = power.shape[1]
645
645
646 listPower = numpy.array_split(power, numBlocks, 0)
646 listPower = numpy.array_split(power, numBlocks, 0)
647 noise = numpy.zeros((power.shape[0], power.shape[1]))
647 noise = numpy.zeros((power.shape[0], power.shape[1]))
648 noise1 = numpy.zeros((power.shape[0], power.shape[1]))
648 noise1 = numpy.zeros((power.shape[0], power.shape[1]))
649
649
650 startInd = 0
650 startInd = 0
651 endInd = 0
651 endInd = 0
652
652
653 for i in range(numBlocks): #split por canal
653 for i in range(numBlocks): #split por canal
654 startInd = endInd
654 startInd = endInd
655 endInd = endInd + listPower[i].shape[0]
655 endInd = endInd + listPower[i].shape[0]
656
656
657 arrayBlock = listPower[i]
657 arrayBlock = listPower[i]
658 noiseAux = numpy.mean(arrayBlock, 0)
658 noiseAux = numpy.mean(arrayBlock, 0)
659 # noiseAux = numpy.median(noiseAux)
659 # noiseAux = numpy.median(noiseAux)
660 # noiseAux = numpy.mean(arrayBlock)
660 # noiseAux = numpy.mean(arrayBlock)
661 noise[startInd:endInd,:] = noise[startInd:endInd,:] + noiseAux
661 noise[startInd:endInd,:] = noise[startInd:endInd,:] + noiseAux
662
662
663 noiseAux1 = numpy.mean(arrayBlock)
663 noiseAux1 = numpy.mean(arrayBlock)
664 noise1[startInd:endInd,:] = noise1[startInd:endInd,:] + noiseAux1
664 noise1[startInd:endInd,:] = noise1[startInd:endInd,:] + noiseAux1
665
665
666 return noise, noise1
666 return noise, noise1
667
667
668 def __findMeteors(self, power, thresh):
668 def __findMeteors(self, power, thresh):
669 nProf = power.shape[0]
669 nProf = power.shape[0]
670 nHeights = power.shape[1]
670 nHeights = power.shape[1]
671 listMeteors = []
671 listMeteors = []
672
672
673 for i in range(nHeights):
673 for i in range(nHeights):
674 powerAux = power[:,i]
674 powerAux = power[:,i]
675 threshAux = thresh[:,i]
675 threshAux = thresh[:,i]
676
676
677 indUPthresh = numpy.where(powerAux > threshAux)[0]
677 indUPthresh = numpy.where(powerAux > threshAux)[0]
678 indDNthresh = numpy.where(powerAux <= threshAux)[0]
678 indDNthresh = numpy.where(powerAux <= threshAux)[0]
679
679
680 j = 0
680 j = 0
681
681
682 while (j < indUPthresh.size - 2):
682 while (j < indUPthresh.size - 2):
683 if (indUPthresh[j + 2] == indUPthresh[j] + 2):
683 if (indUPthresh[j + 2] == indUPthresh[j] + 2):
684 indDNAux = numpy.where(indDNthresh > indUPthresh[j])
684 indDNAux = numpy.where(indDNthresh > indUPthresh[j])
685 indDNthresh = indDNthresh[indDNAux]
685 indDNthresh = indDNthresh[indDNAux]
686
686
687 if (indDNthresh.size > 0):
687 if (indDNthresh.size > 0):
688 indEnd = indDNthresh[0] - 1
688 indEnd = indDNthresh[0] - 1
689 indInit = indUPthresh[j]
689 indInit = indUPthresh[j]
690
690
691 meteor = powerAux[indInit:indEnd + 1]
691 meteor = powerAux[indInit:indEnd + 1]
692 indPeak = meteor.argmax() + indInit
692 indPeak = meteor.argmax() + indInit
693 FLA = sum(numpy.conj(meteor)*numpy.hstack((meteor[1:],0)))
693 FLA = sum(numpy.conj(meteor)*numpy.hstack((meteor[1:],0)))
694
694
695 listMeteors.append(numpy.array([i,indInit,indPeak,indEnd,FLA])) #CHEQUEAR!!!!!
695 listMeteors.append(numpy.array([i,indInit,indPeak,indEnd,FLA])) #CHEQUEAR!!!!!
696 j = numpy.where(indUPthresh == indEnd)[0] + 1
696 j = numpy.where(indUPthresh == indEnd)[0] + 1
697 else: j+=1
697 else: j+=1
698 else: j+=1
698 else: j+=1
699
699
700 return listMeteors
700 return listMeteors
701
701
702 def __removeMultipleDetections(self,listMeteors, rangeLimit, timeLimit):
702 def __removeMultipleDetections(self,listMeteors, rangeLimit, timeLimit):
703
703
704 arrayMeteors = numpy.asarray(listMeteors)
704 arrayMeteors = numpy.asarray(listMeteors)
705 listMeteors1 = []
705 listMeteors1 = []
706
706
707 while arrayMeteors.shape[0] > 0:
707 while arrayMeteors.shape[0] > 0:
708 FLAs = arrayMeteors[:,4]
708 FLAs = arrayMeteors[:,4]
709 maxFLA = FLAs.argmax()
709 maxFLA = FLAs.argmax()
710 listMeteors1.append(arrayMeteors[maxFLA,:])
710 listMeteors1.append(arrayMeteors[maxFLA,:])
711
711
712 MeteorInitTime = arrayMeteors[maxFLA,1]
712 MeteorInitTime = arrayMeteors[maxFLA,1]
713 MeteorEndTime = arrayMeteors[maxFLA,3]
713 MeteorEndTime = arrayMeteors[maxFLA,3]
714 MeteorHeight = arrayMeteors[maxFLA,0]
714 MeteorHeight = arrayMeteors[maxFLA,0]
715
715
716 #Check neighborhood
716 #Check neighborhood
717 maxHeightIndex = MeteorHeight + rangeLimit
717 maxHeightIndex = MeteorHeight + rangeLimit
718 minHeightIndex = MeteorHeight - rangeLimit
718 minHeightIndex = MeteorHeight - rangeLimit
719 minTimeIndex = MeteorInitTime - timeLimit
719 minTimeIndex = MeteorInitTime - timeLimit
720 maxTimeIndex = MeteorEndTime + timeLimit
720 maxTimeIndex = MeteorEndTime + timeLimit
721
721
722 #Check Heights
722 #Check Heights
723 indHeight = numpy.logical_and(arrayMeteors[:,0] >= minHeightIndex, arrayMeteors[:,0] <= maxHeightIndex)
723 indHeight = numpy.logical_and(arrayMeteors[:,0] >= minHeightIndex, arrayMeteors[:,0] <= maxHeightIndex)
724 indTime = numpy.logical_and(arrayMeteors[:,3] >= minTimeIndex, arrayMeteors[:,1] <= maxTimeIndex)
724 indTime = numpy.logical_and(arrayMeteors[:,3] >= minTimeIndex, arrayMeteors[:,1] <= maxTimeIndex)
725 indBoth = numpy.where(numpy.logical_and(indTime,indHeight))
725 indBoth = numpy.where(numpy.logical_and(indTime,indHeight))
726
726
727 arrayMeteors = numpy.delete(arrayMeteors, indBoth, axis = 0)
727 arrayMeteors = numpy.delete(arrayMeteors, indBoth, axis = 0)
728
728
729 return listMeteors1
729 return listMeteors1
730
730
731 def __meteorReestimation(self, listMeteors, volts, pairslist, thresh, noise, timeInterval,frequency):
731 def __meteorReestimation(self, listMeteors, volts, pairslist, thresh, noise, timeInterval,frequency):
732 numHeights = volts.shape[2]
732 numHeights = volts.shape[2]
733 nChannel = volts.shape[0]
733 nChannel = volts.shape[0]
734
734
735 thresholdPhase = thresh[0]
735 thresholdPhase = thresh[0]
736 thresholdNoise = thresh[1]
736 thresholdNoise = thresh[1]
737 thresholdDB = float(thresh[2])
737 thresholdDB = float(thresh[2])
738
738
739 thresholdDB1 = 10**(thresholdDB/10)
739 thresholdDB1 = 10**(thresholdDB/10)
740 pairsarray = numpy.array(pairslist)
740 pairsarray = numpy.array(pairslist)
741 indSides = pairsarray[:,1]
741 indSides = pairsarray[:,1]
742
742
743 pairslist1 = list(pairslist)
743 pairslist1 = list(pairslist)
744 pairslist1.append((0,1))
744 pairslist1.append((0,1))
745 pairslist1.append((3,4))
745 pairslist1.append((3,4))
746
746
747 listMeteors1 = []
747 listMeteors1 = []
748 listPowerSeries = []
748 listPowerSeries = []
749 listVoltageSeries = []
749 listVoltageSeries = []
750 #volts has the war data
750 #volts has the war data
751
751
752 if frequency == 30e6:
752 if frequency == 30e6:
753 timeLag = 45*10**-3
753 timeLag = 45*10**-3
754 else:
754 else:
755 timeLag = 15*10**-3
755 timeLag = 15*10**-3
756 lag = numpy.ceil(timeLag/timeInterval)
756 lag = numpy.ceil(timeLag/timeInterval)
757
757
758 for i in range(len(listMeteors)):
758 for i in range(len(listMeteors)):
759
759
760 ###################### 3.6 - 3.7 PARAMETERS REESTIMATION #########################
760 ###################### 3.6 - 3.7 PARAMETERS REESTIMATION #########################
761 meteorAux = numpy.zeros(16)
761 meteorAux = numpy.zeros(16)
762
762
763 #Loading meteor Data (mHeight, mStart, mPeak, mEnd)
763 #Loading meteor Data (mHeight, mStart, mPeak, mEnd)
764 mHeight = listMeteors[i][0]
764 mHeight = listMeteors[i][0]
765 mStart = listMeteors[i][1]
765 mStart = listMeteors[i][1]
766 mPeak = listMeteors[i][2]
766 mPeak = listMeteors[i][2]
767 mEnd = listMeteors[i][3]
767 mEnd = listMeteors[i][3]
768
768
769 #get the volt data between the start and end times of the meteor
769 #get the volt data between the start and end times of the meteor
770 meteorVolts = volts[:,mStart:mEnd+1,mHeight]
770 meteorVolts = volts[:,mStart:mEnd+1,mHeight]
771 meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
771 meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
772
772
773 #3.6. Phase Difference estimation
773 #3.6. Phase Difference estimation
774 phaseDiff, aux = self.__estimatePhaseDifference(meteorVolts, pairslist)
774 phaseDiff, aux = self.__estimatePhaseDifference(meteorVolts, pairslist)
775
775
776 #3.7. Phase difference removal & meteor start, peak and end times reestimated
776 #3.7. Phase difference removal & meteor start, peak and end times reestimated
777 #meteorVolts0.- all Channels, all Profiles
777 #meteorVolts0.- all Channels, all Profiles
778 meteorVolts0 = volts[:,:,mHeight]
778 meteorVolts0 = volts[:,:,mHeight]
779 meteorThresh = noise[:,mHeight]*thresholdNoise
779 meteorThresh = noise[:,mHeight]*thresholdNoise
780 meteorNoise = noise[:,mHeight]
780 meteorNoise = noise[:,mHeight]
781 meteorVolts0[indSides,:] = self.__shiftPhase(meteorVolts0[indSides,:], phaseDiff) #Phase Shifting
781 meteorVolts0[indSides,:] = self.__shiftPhase(meteorVolts0[indSides,:], phaseDiff) #Phase Shifting
782 powerNet0 = numpy.nansum(numpy.abs(meteorVolts0)**2, axis = 0) #Power
782 powerNet0 = numpy.nansum(numpy.abs(meteorVolts0)**2, axis = 0) #Power
783
783
784 #Times reestimation
784 #Times reestimation
785 mStart1 = numpy.where(powerNet0[:mPeak] < meteorThresh[:mPeak])[0]
785 mStart1 = numpy.where(powerNet0[:mPeak] < meteorThresh[:mPeak])[0]
786 if mStart1.size > 0:
786 if mStart1.size > 0:
787 mStart1 = mStart1[-1] + 1
787 mStart1 = mStart1[-1] + 1
788
788
789 else:
789 else:
790 mStart1 = mPeak
790 mStart1 = mPeak
791
791
792 mEnd1 = numpy.where(powerNet0[mPeak:] < meteorThresh[mPeak:])[0][0] + mPeak - 1
792 mEnd1 = numpy.where(powerNet0[mPeak:] < meteorThresh[mPeak:])[0][0] + mPeak - 1
793 mEndDecayTime1 = numpy.where(powerNet0[mPeak:] < meteorNoise[mPeak:])[0]
793 mEndDecayTime1 = numpy.where(powerNet0[mPeak:] < meteorNoise[mPeak:])[0]
794 if mEndDecayTime1.size == 0:
794 if mEndDecayTime1.size == 0:
795 mEndDecayTime1 = powerNet0.size
795 mEndDecayTime1 = powerNet0.size
796 else:
796 else:
797 mEndDecayTime1 = mEndDecayTime1[0] + mPeak - 1
797 mEndDecayTime1 = mEndDecayTime1[0] + mPeak - 1
798 # mPeak1 = meteorVolts0[mStart1:mEnd1 + 1].argmax()
798 # mPeak1 = meteorVolts0[mStart1:mEnd1 + 1].argmax()
799
799
800 #meteorVolts1.- all Channels, from start to end
800 #meteorVolts1.- all Channels, from start to end
801 meteorVolts1 = meteorVolts0[:,mStart1:mEnd1 + 1]
801 meteorVolts1 = meteorVolts0[:,mStart1:mEnd1 + 1]
802 meteorVolts2 = meteorVolts0[:,mPeak + lag:mEnd1 + 1]
802 meteorVolts2 = meteorVolts0[:,mPeak + lag:mEnd1 + 1]
803 if meteorVolts2.shape[1] == 0:
803 if meteorVolts2.shape[1] == 0:
804 meteorVolts2 = meteorVolts0[:,mPeak:mEnd1 + 1]
804 meteorVolts2 = meteorVolts0[:,mPeak:mEnd1 + 1]
805 meteorVolts1 = meteorVolts1.reshape(meteorVolts1.shape[0], meteorVolts1.shape[1], 1)
805 meteorVolts1 = meteorVolts1.reshape(meteorVolts1.shape[0], meteorVolts1.shape[1], 1)
806 meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1], 1)
806 meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1], 1)
807 ##################### END PARAMETERS REESTIMATION #########################
807 ##################### END PARAMETERS REESTIMATION #########################
808
808
809 ##################### 3.8 PHASE DIFFERENCE REESTIMATION ########################
809 ##################### 3.8 PHASE DIFFERENCE REESTIMATION ########################
810 # if mEnd1 - mStart1 > 4: #Error Number 6: echo less than 5 samples long; too short for analysis
810 # if mEnd1 - mStart1 > 4: #Error Number 6: echo less than 5 samples long; too short for analysis
811 if meteorVolts2.shape[1] > 0:
811 if meteorVolts2.shape[1] > 0:
812 #Phase Difference re-estimation
812 #Phase Difference re-estimation
813 phaseDiff1, phaseDiffint = self.__estimatePhaseDifference(meteorVolts2, pairslist1) #Phase Difference Estimation
813 phaseDiff1, phaseDiffint = self.__estimatePhaseDifference(meteorVolts2, pairslist1) #Phase Difference Estimation
814 # phaseDiff1, phaseDiffint = self.estimatePhaseDifference(meteorVolts2, pairslist)
814 # phaseDiff1, phaseDiffint = self.estimatePhaseDifference(meteorVolts2, pairslist)
815 meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1])
815 meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1])
816 phaseDiff11 = numpy.reshape(phaseDiff1, (phaseDiff1.shape[0],1))
816 phaseDiff11 = numpy.reshape(phaseDiff1, (phaseDiff1.shape[0],1))
817 meteorVolts2[indSides,:] = self.__shiftPhase(meteorVolts2[indSides,:], phaseDiff11[0:4]) #Phase Shifting
817 meteorVolts2[indSides,:] = self.__shiftPhase(meteorVolts2[indSides,:], phaseDiff11[0:4]) #Phase Shifting
818
818
819 #Phase Difference RMS
819 #Phase Difference RMS
820 phaseRMS1 = numpy.sqrt(numpy.mean(numpy.square(phaseDiff1)))
820 phaseRMS1 = numpy.sqrt(numpy.mean(numpy.square(phaseDiff1)))
821 powerNet1 = numpy.nansum(numpy.abs(meteorVolts1[:,:])**2,0)
821 powerNet1 = numpy.nansum(numpy.abs(meteorVolts1[:,:])**2,0)
822 #Data from Meteor
822 #Data from Meteor
823 mPeak1 = powerNet1.argmax() + mStart1
823 mPeak1 = powerNet1.argmax() + mStart1
824 mPeakPower1 = powerNet1.max()
824 mPeakPower1 = powerNet1.max()
825 noiseAux = sum(noise[mStart1:mEnd1 + 1,mHeight])
825 noiseAux = sum(noise[mStart1:mEnd1 + 1,mHeight])
826 mSNR1 = (sum(powerNet1)-noiseAux)/noiseAux
826 mSNR1 = (sum(powerNet1)-noiseAux)/noiseAux
827 Meteor1 = numpy.array([mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1])
827 Meteor1 = numpy.array([mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1])
828 Meteor1 = numpy.hstack((Meteor1,phaseDiffint))
828 Meteor1 = numpy.hstack((Meteor1,phaseDiffint))
829 PowerSeries = powerNet0[mStart1:mEndDecayTime1 + 1]
829 PowerSeries = powerNet0[mStart1:mEndDecayTime1 + 1]
830 #Vectorize
830 #Vectorize
831 meteorAux[0:7] = [mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1]
831 meteorAux[0:7] = [mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1]
832 meteorAux[7:11] = phaseDiffint[0:4]
832 meteorAux[7:11] = phaseDiffint[0:4]
833
833
834 #Rejection Criterions
834 #Rejection Criterions
835 if phaseRMS1 > thresholdPhase: #Error Number 17: Phase variation
835 if phaseRMS1 > thresholdPhase: #Error Number 17: Phase variation
836 meteorAux[-1] = 17
836 meteorAux[-1] = 17
837 elif mSNR1 < thresholdDB1: #Error Number 1: SNR < threshold dB
837 elif mSNR1 < thresholdDB1: #Error Number 1: SNR < threshold dB
838 meteorAux[-1] = 1
838 meteorAux[-1] = 1
839
839
840
840
841 else:
841 else:
842 meteorAux[0:4] = [mHeight, mStart, mPeak, mEnd]
842 meteorAux[0:4] = [mHeight, mStart, mPeak, mEnd]
843 meteorAux[-1] = 6 #Error Number 6: echo less than 5 samples long; too short for analysis
843 meteorAux[-1] = 6 #Error Number 6: echo less than 5 samples long; too short for analysis
844 PowerSeries = 0
844 PowerSeries = 0
845
845
846 listMeteors1.append(meteorAux)
846 listMeteors1.append(meteorAux)
847 listPowerSeries.append(PowerSeries)
847 listPowerSeries.append(PowerSeries)
848 listVoltageSeries.append(meteorVolts1)
848 listVoltageSeries.append(meteorVolts1)
849
849
850 return listMeteors1, listPowerSeries, listVoltageSeries
850 return listMeteors1, listPowerSeries, listVoltageSeries
851
851
852 def __estimateDecayTime(self, listMeteors, listPower, timeInterval, frequency):
852 def __estimateDecayTime(self, listMeteors, listPower, timeInterval, frequency):
853
853
854 threshError = 10
854 threshError = 10
855 #Depending if it is 30 or 50 MHz
855 #Depending if it is 30 or 50 MHz
856 if frequency == 30e6:
856 if frequency == 30e6:
857 timeLag = 45*10**-3
857 timeLag = 45*10**-3
858 else:
858 else:
859 timeLag = 15*10**-3
859 timeLag = 15*10**-3
860 lag = numpy.ceil(timeLag/timeInterval)
860 lag = numpy.ceil(timeLag/timeInterval)
861
861
862 listMeteors1 = []
862 listMeteors1 = []
863
863
864 for i in range(len(listMeteors)):
864 for i in range(len(listMeteors)):
865 meteorPower = listPower[i]
865 meteorPower = listPower[i]
866 meteorAux = listMeteors[i]
866 meteorAux = listMeteors[i]
867
867
868 if meteorAux[-1] == 0:
868 if meteorAux[-1] == 0:
869
869
870 try:
870 try:
871 indmax = meteorPower.argmax()
871 indmax = meteorPower.argmax()
872 indlag = indmax + lag
872 indlag = indmax + lag
873
873
874 y = meteorPower[indlag:]
874 y = meteorPower[indlag:]
875 x = numpy.arange(0, y.size)*timeLag
875 x = numpy.arange(0, y.size)*timeLag
876
876
877 #first guess
877 #first guess
878 a = y[0]
878 a = y[0]
879 tau = timeLag
879 tau = timeLag
880 #exponential fit
880 #exponential fit
881 popt, pcov = optimize.curve_fit(self.__exponential_function, x, y, p0 = [a, tau])
881 popt, pcov = optimize.curve_fit(self.__exponential_function, x, y, p0 = [a, tau])
882 y1 = self.__exponential_function(x, *popt)
882 y1 = self.__exponential_function(x, *popt)
883 #error estimation
883 #error estimation
884 error = sum((y - y1)**2)/(numpy.var(y)*(y.size - popt.size))
884 error = sum((y - y1)**2)/(numpy.var(y)*(y.size - popt.size))
885
885
886 decayTime = popt[1]
886 decayTime = popt[1]
887 riseTime = indmax*timeInterval
887 riseTime = indmax*timeInterval
888 meteorAux[11:13] = [decayTime, error]
888 meteorAux[11:13] = [decayTime, error]
889
889
890 #Table items 7, 8 and 11
890 #Table items 7, 8 and 11
891 if (riseTime > 0.3): #Number 7: Echo rise exceeds 0.3s
891 if (riseTime > 0.3): #Number 7: Echo rise exceeds 0.3s
892 meteorAux[-1] = 7
892 meteorAux[-1] = 7
893 elif (decayTime < 2*riseTime) : #Number 8: Echo decay time less than than twice rise time
893 elif (decayTime < 2*riseTime) : #Number 8: Echo decay time less than than twice rise time
894 meteorAux[-1] = 8
894 meteorAux[-1] = 8
895 if (error > threshError): #Number 11: Poor fit to amplitude for estimation of decay time
895 if (error > threshError): #Number 11: Poor fit to amplitude for estimation of decay time
896 meteorAux[-1] = 11
896 meteorAux[-1] = 11
897
897
898
898
899 except:
899 except:
900 meteorAux[-1] = 11
900 meteorAux[-1] = 11
901
901
902
902
903 listMeteors1.append(meteorAux)
903 listMeteors1.append(meteorAux)
904
904
905 return listMeteors1
905 return listMeteors1
906
906
907 #Exponential Function
907 #Exponential Function
908
908
909 def __exponential_function(self, x, a, tau):
909 def __exponential_function(self, x, a, tau):
910 y = a*numpy.exp(-x/tau)
910 y = a*numpy.exp(-x/tau)
911 return y
911 return y
912
912
913 def __getRadialVelocity(self, listMeteors, listVolts, radialStdThresh, pairslist, timeInterval):
913 def __getRadialVelocity(self, listMeteors, listVolts, radialStdThresh, pairslist, timeInterval):
914
914
915 pairslist1 = list(pairslist)
915 pairslist1 = list(pairslist)
916 pairslist1.append((0,1))
916 pairslist1.append((0,1))
917 pairslist1.append((3,4))
917 pairslist1.append((3,4))
918 numPairs = len(pairslist1)
918 numPairs = len(pairslist1)
919 #Time Lag
919 #Time Lag
920 timeLag = 45*10**-3
920 timeLag = 45*10**-3
921 c = 3e8
921 c = 3e8
922 lag = numpy.ceil(timeLag/timeInterval)
922 lag = numpy.ceil(timeLag/timeInterval)
923 freq = 30e6
923 freq = 30e6
924
924
925 listMeteors1 = []
925 listMeteors1 = []
926
926
927 for i in range(len(listMeteors)):
927 for i in range(len(listMeteors)):
928 meteor = listMeteors[i]
928 meteor = listMeteors[i]
929 meteorAux = numpy.hstack((meteor[:-1], 0, 0, meteor[-1]))
929 meteorAux = numpy.hstack((meteor[:-1], 0, 0, meteor[-1]))
930 if meteor[-1] == 0:
930 if meteor[-1] == 0:
931 mStart = listMeteors[i][1]
931 mStart = listMeteors[i][1]
932 mPeak = listMeteors[i][2]
932 mPeak = listMeteors[i][2]
933 mLag = mPeak - mStart + lag
933 mLag = mPeak - mStart + lag
934
934
935 #get the volt data between the start and end times of the meteor
935 #get the volt data between the start and end times of the meteor
936 meteorVolts = listVolts[i]
936 meteorVolts = listVolts[i]
937 meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
937 meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
938
938
939 #Get CCF
939 #Get CCF
940 allCCFs = self.__calculateCCF(meteorVolts, pairslist1, [-2,-1,0,1,2])
940 allCCFs = self.__calculateCCF(meteorVolts, pairslist1, [-2,-1,0,1,2])
941
941
942 #Method 2
942 #Method 2
943 slopes = numpy.zeros(numPairs)
943 slopes = numpy.zeros(numPairs)
944 time = numpy.array([-2,-1,1,2])*timeInterval
944 time = numpy.array([-2,-1,1,2])*timeInterval
945 angAllCCF = numpy.angle(allCCFs[:,[0,1,3,4],0])
945 angAllCCF = numpy.angle(allCCFs[:,[0,1,3,4],0])
946
946
947 #Correct phases
947 #Correct phases
948 derPhaseCCF = angAllCCF[:,1:] - angAllCCF[:,0:-1]
948 derPhaseCCF = angAllCCF[:,1:] - angAllCCF[:,0:-1]
949 indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
949 indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
950
950
951 if indDer[0].shape[0] > 0:
951 if indDer[0].shape[0] > 0:
952 for i in range(indDer[0].shape[0]):
952 for i in range(indDer[0].shape[0]):
953 signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i]])
953 signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i]])
954 angAllCCF[indDer[0][i],indDer[1][i]+1:] += signo*2*numpy.pi
954 angAllCCF[indDer[0][i],indDer[1][i]+1:] += signo*2*numpy.pi
955
955
956 # fit = scipy.stats.linregress(numpy.array([-2,-1,1,2])*timeInterval, numpy.array([phaseLagN2s[i],phaseLagN1s[i],phaseLag1s[i],phaseLag2s[i]]))
956 # fit = scipy.stats.linregress(numpy.array([-2,-1,1,2])*timeInterval, numpy.array([phaseLagN2s[i],phaseLagN1s[i],phaseLag1s[i],phaseLag2s[i]]))
957 for j in range(numPairs):
957 for j in range(numPairs):
958 fit = stats.linregress(time, angAllCCF[j,:])
958 fit = stats.linregress(time, angAllCCF[j,:])
959 slopes[j] = fit[0]
959 slopes[j] = fit[0]
960
960
961 #Remove Outlier
961 #Remove Outlier
962 # indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
962 # indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
963 # slopes = numpy.delete(slopes,indOut)
963 # slopes = numpy.delete(slopes,indOut)
964 # indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
964 # indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
965 # slopes = numpy.delete(slopes,indOut)
965 # slopes = numpy.delete(slopes,indOut)
966
966
967 radialVelocity = -numpy.mean(slopes)*(0.25/numpy.pi)*(c/freq)
967 radialVelocity = -numpy.mean(slopes)*(0.25/numpy.pi)*(c/freq)
968 radialError = numpy.std(slopes)*(0.25/numpy.pi)*(c/freq)
968 radialError = numpy.std(slopes)*(0.25/numpy.pi)*(c/freq)
969 meteorAux[-2] = radialError
969 meteorAux[-2] = radialError
970 meteorAux[-3] = radialVelocity
970 meteorAux[-3] = radialVelocity
971
971
972 #Setting Error
972 #Setting Error
973 #Number 15: Radial Drift velocity or projected horizontal velocity exceeds 200 m/s
973 #Number 15: Radial Drift velocity or projected horizontal velocity exceeds 200 m/s
974 if numpy.abs(radialVelocity) > 200:
974 if numpy.abs(radialVelocity) > 200:
975 meteorAux[-1] = 15
975 meteorAux[-1] = 15
976 #Number 12: Poor fit to CCF variation for estimation of radial drift velocity
976 #Number 12: Poor fit to CCF variation for estimation of radial drift velocity
977 elif radialError > radialStdThresh:
977 elif radialError > radialStdThresh:
978 meteorAux[-1] = 12
978 meteorAux[-1] = 12
979
979
980 listMeteors1.append(meteorAux)
980 listMeteors1.append(meteorAux)
981 return listMeteors1
981 return listMeteors1
982
982
983 def __setNewArrays(self, listMeteors, date, heiRang):
983 def __setNewArrays(self, listMeteors, date, heiRang):
984
984
985 #New arrays
985 #New arrays
986 arrayMeteors = numpy.array(listMeteors)
986 arrayMeteors = numpy.array(listMeteors)
987 arrayParameters = numpy.zeros((len(listMeteors),10))
987 arrayParameters = numpy.zeros((len(listMeteors),10))
988
988
989 #Date inclusion
989 #Date inclusion
990 date = re.findall(r'\((.*?)\)', date)
990 date = re.findall(r'\((.*?)\)', date)
991 date = date[0].split(',')
991 date = date[0].split(',')
992 date = map(int, date)
992 date = map(int, date)
993 date = [date[0]*10000 + date[1]*100 + date[2], date[3]*10000 + date[4]*100 + date[5]]
993 date = [date[0]*10000 + date[1]*100 + date[2], date[3]*10000 + date[4]*100 + date[5]]
994 arrayDate = numpy.tile(date, (len(listMeteors), 1))
994 arrayDate = numpy.tile(date, (len(listMeteors), 1))
995
995
996 #Meteor array
996 #Meteor array
997 arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)]
997 arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)]
998 arrayMeteors = numpy.hstack((arrayDate, arrayMeteors))
998 arrayMeteors = numpy.hstack((arrayDate, arrayMeteors))
999
999
1000 #Parameters Array
1000 #Parameters Array
1001 arrayParameters[:,0:3] = arrayMeteors[:,0:3]
1001 arrayParameters[:,0:3] = arrayMeteors[:,0:3]
1002 arrayParameters[:,-3:] = arrayMeteors[:,-3:]
1002 arrayParameters[:,-3:] = arrayMeteors[:,-3:]
1003
1003
1004 return arrayMeteors, arrayParameters
1004 return arrayMeteors, arrayParameters
1005
1005
1006 def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth):
1006 def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth):
1007
1007
1008 arrayAOA = numpy.zeros((phases.shape[0],3))
1008 arrayAOA = numpy.zeros((phases.shape[0],3))
1009 cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList)
1009 cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList)
1010
1010
1011 arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
1011 arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
1012 cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
1012 cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
1013 arrayAOA[:,2] = cosDirError
1013 arrayAOA[:,2] = cosDirError
1014
1014
1015 azimuthAngle = arrayAOA[:,0]
1015 azimuthAngle = arrayAOA[:,0]
1016 zenithAngle = arrayAOA[:,1]
1016 zenithAngle = arrayAOA[:,1]
1017
1017
1018 #Setting Error
1018 #Setting Error
1019 #Number 3: AOA not fesible
1019 #Number 3: AOA not fesible
1020 indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
1020 indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
1021 error[indInvalid] = 3
1021 error[indInvalid] = 3
1022 #Number 4: Large difference in AOAs obtained from different antenna baselines
1022 #Number 4: Large difference in AOAs obtained from different antenna baselines
1023 indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
1023 indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
1024 error[indInvalid] = 4
1024 error[indInvalid] = 4
1025 return arrayAOA, error
1025 return arrayAOA, error
1026
1026
1027 def __getDirectionCosines(self, arrayPhase, pairsList):
1027 def __getDirectionCosines(self, arrayPhase, pairsList):
1028
1028
1029 #Initializing some variables
1029 #Initializing some variables
1030 ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
1030 ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
1031 ang_aux = ang_aux.reshape(1,ang_aux.size)
1031 ang_aux = ang_aux.reshape(1,ang_aux.size)
1032
1032
1033 cosdir = numpy.zeros((arrayPhase.shape[0],2))
1033 cosdir = numpy.zeros((arrayPhase.shape[0],2))
1034 cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
1034 cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
1035
1035
1036
1036
1037 for i in range(2):
1037 for i in range(2):
1038 #First Estimation
1038 #First Estimation
1039 phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]]
1039 phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]]
1040 #Dealias
1040 #Dealias
1041 indcsi = numpy.where(phi0_aux > numpy.pi)
1041 indcsi = numpy.where(phi0_aux > numpy.pi)
1042 phi0_aux[indcsi] -= 2*numpy.pi
1042 phi0_aux[indcsi] -= 2*numpy.pi
1043 indcsi = numpy.where(phi0_aux < -numpy.pi)
1043 indcsi = numpy.where(phi0_aux < -numpy.pi)
1044 phi0_aux[indcsi] += 2*numpy.pi
1044 phi0_aux[indcsi] += 2*numpy.pi
1045 #Direction Cosine 0
1045 #Direction Cosine 0
1046 cosdir0[:,i] = -(phi0_aux)/(2*numpy.pi*0.5)
1046 cosdir0[:,i] = -(phi0_aux)/(2*numpy.pi*0.5)
1047
1047
1048 #Most-Accurate Second Estimation
1048 #Most-Accurate Second Estimation
1049 phi1_aux = arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]]
1049 phi1_aux = arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]]
1050 phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
1050 phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
1051 #Direction Cosine 1
1051 #Direction Cosine 1
1052 cosdir1 = -(phi1_aux + ang_aux)/(2*numpy.pi*4.5)
1052 cosdir1 = -(phi1_aux + ang_aux)/(2*numpy.pi*4.5)
1053
1053
1054 #Searching the correct Direction Cosine
1054 #Searching the correct Direction Cosine
1055 cosdir0_aux = cosdir0[:,i]
1055 cosdir0_aux = cosdir0[:,i]
1056 cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
1056 cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
1057 #Minimum Distance
1057 #Minimum Distance
1058 cosDiff = (cosdir1 - cosdir0_aux)**2
1058 cosDiff = (cosdir1 - cosdir0_aux)**2
1059 indcos = cosDiff.argmin(axis = 1)
1059 indcos = cosDiff.argmin(axis = 1)
1060 #Saving Value obtained
1060 #Saving Value obtained
1061 cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
1061 cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
1062
1062
1063 return cosdir0, cosdir
1063 return cosdir0, cosdir
1064
1064
1065 def __calculateAOA(self, cosdir, azimuth):
1065 def __calculateAOA(self, cosdir, azimuth):
1066 cosdirX = cosdir[:,0]
1066 cosdirX = cosdir[:,0]
1067 cosdirY = cosdir[:,1]
1067 cosdirY = cosdir[:,1]
1068
1068
1069 zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
1069 zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
1070 azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth #0 deg north, 90 deg east
1070 azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth #0 deg north, 90 deg east
1071 angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
1071 angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
1072
1072
1073 return angles
1073 return angles
1074
1074
1075 def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
1075 def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
1076
1076
1077 Ramb = 375 #Ramb = c/(2*PRF)
1077 Ramb = 375 #Ramb = c/(2*PRF)
1078 Re = 6371 #Earth Radius
1078 Re = 6371 #Earth Radius
1079 heights = numpy.zeros(Ranges.shape)
1079 heights = numpy.zeros(Ranges.shape)
1080
1080
1081 R_aux = numpy.array([0,1,2])*Ramb
1081 R_aux = numpy.array([0,1,2])*Ramb
1082 R_aux = R_aux.reshape(1,R_aux.size)
1082 R_aux = R_aux.reshape(1,R_aux.size)
1083
1083
1084 Ranges = Ranges.reshape(Ranges.size,1)
1084 Ranges = Ranges.reshape(Ranges.size,1)
1085
1085
1086 Ri = Ranges + R_aux
1086 Ri = Ranges + R_aux
1087 hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re
1087 hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re
1088
1088
1089 #Check if there is a height between 70 and 110 km
1089 #Check if there is a height between 70 and 110 km
1090 h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
1090 h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
1091 ind_h = numpy.where(h_bool == 1)[0]
1091 ind_h = numpy.where(h_bool == 1)[0]
1092
1092
1093 hCorr = hi[ind_h, :]
1093 hCorr = hi[ind_h, :]
1094 ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
1094 ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
1095
1095
1096 hCorr = hi[ind_hCorr]
1096 hCorr = hi[ind_hCorr]
1097 heights[ind_h] = hCorr
1097 heights[ind_h] = hCorr
1098
1098
1099 #Setting Error
1099 #Setting Error
1100 #Number 13: Height unresolvable echo: not valid height within 70 to 110 km
1100 #Number 13: Height unresolvable echo: not valid height within 70 to 110 km
1101 #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
1101 #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
1102
1102
1103 indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
1103 indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
1104 error[indInvalid2] = 14
1104 error[indInvalid2] = 14
1105 indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
1105 indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
1106 error[indInvalid1] = 13
1106 error[indInvalid1] = 13
1107
1107
1108 return heights, error
1108 return heights, error
1109
1109
1110
1110
1111 class WindProfiler(Operation):
1111 class WindProfiler(Operation):
1112
1112
1113 __isConfig = False
1113 __isConfig = False
1114
1114
1115 __initime = None
1115 __initime = None
1116 __lastdatatime = None
1116 __lastdatatime = None
1117 __integrationtime = None
1117 __integrationtime = None
1118
1118
1119 __buffer = None
1119 __buffer = None
1120
1120
1121 __dataReady = False
1121 __dataReady = False
1122
1122
1123 __firstdata = None
1123 __firstdata = None
1124
1124
1125 n = None
1125 n = None
1126
1126
1127 def __init__(self):
1127 def __init__(self):
1128 Operation.__init__(self)
1128 Operation.__init__(self)
1129
1129
1130 def __calculateCosDir(self, elev, azim):
1131 zen = (90 - elev)*numpy.pi/180
1132 azim = azim*numpy.pi/180
1133 cosDirX = numpy.sqrt((1-numpy.cos(zen)**2)/((1+numpy.tan(azim)**2)))
1134 cosDirY = numpy.sqrt(1-numpy.cos(zen)**2-cosDirX**2)
1135
1136 signX = numpy.sign(numpy.cos(azim))
1137 signY = numpy.sign(numpy.sin(azim))
1138
1139 cosDirX = numpy.copysign(cosDirX, signX)
1140 cosDirY = numpy.copysign(cosDirY, signY)
1141 return cosDirX, cosDirY
1142
1130 def __calculateAngles(self, theta_x, theta_y, azimuth):
1143 def __calculateAngles(self, theta_x, theta_y, azimuth):
1131
1144
1132 dir_cosw = numpy.sqrt(1-theta_x**2-theta_y**2)
1145 dir_cosw = numpy.sqrt(1-theta_x**2-theta_y**2)
1133 zenith_arr = numpy.arccos(dir_cosw)
1146 zenith_arr = numpy.arccos(dir_cosw)
1134 azimuth_arr = numpy.arctan2(theta_x,theta_y) + azimuth*math.pi/180
1147 azimuth_arr = numpy.arctan2(theta_x,theta_y) + azimuth*math.pi/180
1135
1148
1136 dir_cosu = numpy.sin(azimuth_arr)*numpy.sin(zenith_arr)
1149 dir_cosu = numpy.sin(azimuth_arr)*numpy.sin(zenith_arr)
1137 dir_cosv = numpy.cos(azimuth_arr)*numpy.sin(zenith_arr)
1150 dir_cosv = numpy.cos(azimuth_arr)*numpy.sin(zenith_arr)
1138
1151
1139 return azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw
1152 return azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw
1140
1153
1141 def __calculateMatA(self, dir_cosu, dir_cosv, dir_cosw, horOnly):
1154 def __calculateMatA(self, dir_cosu, dir_cosv, dir_cosw, horOnly):
1142
1155
1143 #
1156 #
1144 if horOnly:
1157 if horOnly:
1145 A = numpy.c_[dir_cosu,dir_cosv]
1158 A = numpy.c_[dir_cosu,dir_cosv]
1146 else:
1159 else:
1147 A = numpy.c_[dir_cosu,dir_cosv,dir_cosw]
1160 A = numpy.c_[dir_cosu,dir_cosv,dir_cosw]
1148 A = numpy.asmatrix(A)
1161 A = numpy.asmatrix(A)
1149 A1 = numpy.linalg.inv(A.transpose()*A)*A.transpose()
1162 A1 = numpy.linalg.inv(A.transpose()*A)*A.transpose()
1150
1163
1151 return A1
1164 return A1
1152
1165
1153 def __correctValues(self, heiRang, phi, velRadial, SNR):
1166 def __correctValues(self, heiRang, phi, velRadial, SNR):
1154 listPhi = phi.tolist()
1167 listPhi = phi.tolist()
1155 maxid = listPhi.index(max(listPhi))
1168 maxid = listPhi.index(max(listPhi))
1156 minid = listPhi.index(min(listPhi))
1169 minid = listPhi.index(min(listPhi))
1157
1170
1158 rango = range(len(phi))
1171 rango = range(len(phi))
1159 # rango = numpy.delete(rango,maxid)
1172 # rango = numpy.delete(rango,maxid)
1160
1173
1161 heiRang1 = heiRang*math.cos(phi[maxid])
1174 heiRang1 = heiRang*math.cos(phi[maxid])
1162 heiRangAux = heiRang*math.cos(phi[minid])
1175 heiRangAux = heiRang*math.cos(phi[minid])
1163 indOut = (heiRang1 < heiRangAux[0]).nonzero()
1176 indOut = (heiRang1 < heiRangAux[0]).nonzero()
1164 heiRang1 = numpy.delete(heiRang1,indOut)
1177 heiRang1 = numpy.delete(heiRang1,indOut)
1165
1178
1166 velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
1179 velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
1167 SNR1 = numpy.zeros([len(phi),len(heiRang1)])
1180 SNR1 = numpy.zeros([len(phi),len(heiRang1)])
1168
1181
1169 for i in rango:
1182 for i in rango:
1170 x = heiRang*math.cos(phi[i])
1183 x = heiRang*math.cos(phi[i])
1171 y1 = velRadial[i,:]
1184 y1 = velRadial[i,:]
1172 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
1185 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
1173
1186
1174 x1 = heiRang1
1187 x1 = heiRang1
1175 y11 = f1(x1)
1188 y11 = f1(x1)
1176
1189
1177 y2 = SNR[i,:]
1190 y2 = SNR[i,:]
1178 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
1191 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
1179 y21 = f2(x1)
1192 y21 = f2(x1)
1180
1193
1181 velRadial1[i,:] = y11
1194 velRadial1[i,:] = y11
1182 SNR1[i,:] = y21
1195 SNR1[i,:] = y21
1183
1196
1184 return heiRang1, velRadial1, SNR1
1197 return heiRang1, velRadial1, SNR1
1185
1198
1186 def __calculateVelUVW(self, A, velRadial):
1199 def __calculateVelUVW(self, A, velRadial):
1187
1200
1188 #Operacion Matricial
1201 #Operacion Matricial
1189 # velUVW = numpy.zeros((velRadial.shape[1],3))
1202 # velUVW = numpy.zeros((velRadial.shape[1],3))
1190 # for ind in range(velRadial.shape[1]):
1203 # for ind in range(velRadial.shape[1]):
1191 # velUVW[ind,:] = numpy.dot(A,velRadial[:,ind])
1204 # velUVW[ind,:] = numpy.dot(A,velRadial[:,ind])
1192 # velUVW = velUVW.transpose()
1205 # velUVW = velUVW.transpose()
1193 velUVW = numpy.zeros((A.shape[0],velRadial.shape[1]))
1206 velUVW = numpy.zeros((A.shape[0],velRadial.shape[1]))
1194 velUVW[:,:] = numpy.dot(A,velRadial)
1207 velUVW[:,:] = numpy.dot(A,velRadial)
1195
1208
1196
1209
1197 return velUVW
1210 return velUVW
1198
1211
1199 def techniqueDBS(self, velRadial0, dirCosx, disrCosy, azimuth, correct, horizontalOnly, heiRang, SNR0):
1212 def techniqueDBS(self, velRadial0, dirCosx, disrCosy, azimuth, correct, horizontalOnly, heiRang, SNR0):
1200 """
1213 """
1201 Function that implements Doppler Beam Swinging (DBS) technique.
1214 Function that implements Doppler Beam Swinging (DBS) technique.
1202
1215
1203 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
1216 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
1204 Direction correction (if necessary), Ranges and SNR
1217 Direction correction (if necessary), Ranges and SNR
1205
1218
1206 Output: Winds estimation (Zonal, Meridional and Vertical)
1219 Output: Winds estimation (Zonal, Meridional and Vertical)
1207
1220
1208 Parameters affected: Winds, height range, SNR
1221 Parameters affected: Winds, height range, SNR
1209 """
1222 """
1210 azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(dirCosx, disrCosy, azimuth)
1223 azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(dirCosx, disrCosy, azimuth)
1211 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correct*velRadial0, SNR0)
1224 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correct*velRadial0, SNR0)
1212 A = self.__calculateMatA(dir_cosu, dir_cosv, dir_cosw, horizontalOnly)
1225 A = self.__calculateMatA(dir_cosu, dir_cosv, dir_cosw, horizontalOnly)
1213
1226
1214 #Calculo de Componentes de la velocidad con DBS
1227 #Calculo de Componentes de la velocidad con DBS
1215 winds = self.__calculateVelUVW(A,velRadial1)
1228 winds = self.__calculateVelUVW(A,velRadial1)
1216
1229
1217 return winds, heiRang1, SNR1
1230 return winds, heiRang1, SNR1
1218
1231
1219 def __calculateDistance(self, posx, posy, pairsCrossCorr, pairsList, pairs, azimuth = None):
1232 def __calculateDistance(self, posx, posy, pairsCrossCorr, pairsList, pairs, azimuth = None):
1220
1233
1221 posx = numpy.asarray(posx)
1234 posx = numpy.asarray(posx)
1222 posy = numpy.asarray(posy)
1235 posy = numpy.asarray(posy)
1223
1236
1224 #Rotacion Inversa para alinear con el azimuth
1237 #Rotacion Inversa para alinear con el azimuth
1225 if azimuth!= None:
1238 if azimuth!= None:
1226 azimuth = azimuth*math.pi/180
1239 azimuth = azimuth*math.pi/180
1227 posx1 = posx*math.cos(azimuth) + posy*math.sin(azimuth)
1240 posx1 = posx*math.cos(azimuth) + posy*math.sin(azimuth)
1228 posy1 = -posx*math.sin(azimuth) + posy*math.cos(azimuth)
1241 posy1 = -posx*math.sin(azimuth) + posy*math.cos(azimuth)
1229 else:
1242 else:
1230 posx1 = posx
1243 posx1 = posx
1231 posy1 = posy
1244 posy1 = posy
1232
1245
1233 #Calculo de Distancias
1246 #Calculo de Distancias
1234 distx = numpy.zeros(pairsCrossCorr.size)
1247 distx = numpy.zeros(pairsCrossCorr.size)
1235 disty = numpy.zeros(pairsCrossCorr.size)
1248 disty = numpy.zeros(pairsCrossCorr.size)
1236 dist = numpy.zeros(pairsCrossCorr.size)
1249 dist = numpy.zeros(pairsCrossCorr.size)
1237 ang = numpy.zeros(pairsCrossCorr.size)
1250 ang = numpy.zeros(pairsCrossCorr.size)
1238
1251
1239 for i in range(pairsCrossCorr.size):
1252 for i in range(pairsCrossCorr.size):
1240 distx[i] = posx1[pairsList[pairsCrossCorr[i]][1]] - posx1[pairsList[pairsCrossCorr[i]][0]]
1253 distx[i] = posx1[pairsList[pairsCrossCorr[i]][1]] - posx1[pairsList[pairsCrossCorr[i]][0]]
1241 disty[i] = posy1[pairsList[pairsCrossCorr[i]][1]] - posy1[pairsList[pairsCrossCorr[i]][0]]
1254 disty[i] = posy1[pairsList[pairsCrossCorr[i]][1]] - posy1[pairsList[pairsCrossCorr[i]][0]]
1242 dist[i] = numpy.sqrt(distx[i]**2 + disty[i]**2)
1255 dist[i] = numpy.sqrt(distx[i]**2 + disty[i]**2)
1243 ang[i] = numpy.arctan2(disty[i],distx[i])
1256 ang[i] = numpy.arctan2(disty[i],distx[i])
1244 #Calculo de Matrices
1257 #Calculo de Matrices
1245 nPairs = len(pairs)
1258 nPairs = len(pairs)
1246 ang1 = numpy.zeros((nPairs, 2, 1))
1259 ang1 = numpy.zeros((nPairs, 2, 1))
1247 dist1 = numpy.zeros((nPairs, 2, 1))
1260 dist1 = numpy.zeros((nPairs, 2, 1))
1248
1261
1249 for j in range(nPairs):
1262 for j in range(nPairs):
1250 dist1[j,0,0] = dist[pairs[j][0]]
1263 dist1[j,0,0] = dist[pairs[j][0]]
1251 dist1[j,1,0] = dist[pairs[j][1]]
1264 dist1[j,1,0] = dist[pairs[j][1]]
1252 ang1[j,0,0] = ang[pairs[j][0]]
1265 ang1[j,0,0] = ang[pairs[j][0]]
1253 ang1[j,1,0] = ang[pairs[j][1]]
1266 ang1[j,1,0] = ang[pairs[j][1]]
1254
1267
1255 return distx,disty, dist1,ang1
1268 return distx,disty, dist1,ang1
1256
1269
1257 def __calculateVelVer(self, phase, lagTRange, _lambda):
1270 def __calculateVelVer(self, phase, lagTRange, _lambda):
1258
1271
1259 Ts = lagTRange[1] - lagTRange[0]
1272 Ts = lagTRange[1] - lagTRange[0]
1260 velW = -_lambda*phase/(4*math.pi*Ts)
1273 velW = -_lambda*phase/(4*math.pi*Ts)
1261
1274
1262 return velW
1275 return velW
1263
1276
1264 def __calculateVelHorDir(self, dist, tau1, tau2, ang):
1277 def __calculateVelHorDir(self, dist, tau1, tau2, ang):
1265 nPairs = tau1.shape[0]
1278 nPairs = tau1.shape[0]
1266 vel = numpy.zeros((nPairs,3,tau1.shape[2]))
1279 vel = numpy.zeros((nPairs,3,tau1.shape[2]))
1267
1280
1268 angCos = numpy.cos(ang)
1281 angCos = numpy.cos(ang)
1269 angSin = numpy.sin(ang)
1282 angSin = numpy.sin(ang)
1270
1283
1271 vel0 = dist*tau1/(2*tau2**2)
1284 vel0 = dist*tau1/(2*tau2**2)
1272 vel[:,0,:] = (vel0*angCos).sum(axis = 1)
1285 vel[:,0,:] = (vel0*angCos).sum(axis = 1)
1273 vel[:,1,:] = (vel0*angSin).sum(axis = 1)
1286 vel[:,1,:] = (vel0*angSin).sum(axis = 1)
1274
1287
1275 ind = numpy.where(numpy.isinf(vel))
1288 ind = numpy.where(numpy.isinf(vel))
1276 vel[ind] = numpy.nan
1289 vel[ind] = numpy.nan
1277
1290
1278 return vel
1291 return vel
1279
1292
1280 def __getPairsAutoCorr(self, pairsList, nChannels):
1293 def __getPairsAutoCorr(self, pairsList, nChannels):
1281
1294
1282 pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
1295 pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
1283
1296
1284 for l in range(len(pairsList)):
1297 for l in range(len(pairsList)):
1285 firstChannel = pairsList[l][0]
1298 firstChannel = pairsList[l][0]
1286 secondChannel = pairsList[l][1]
1299 secondChannel = pairsList[l][1]
1287
1300
1288 #Obteniendo pares de Autocorrelacion
1301 #Obteniendo pares de Autocorrelacion
1289 if firstChannel == secondChannel:
1302 if firstChannel == secondChannel:
1290 pairsAutoCorr[firstChannel] = int(l)
1303 pairsAutoCorr[firstChannel] = int(l)
1291
1304
1292 pairsAutoCorr = pairsAutoCorr.astype(int)
1305 pairsAutoCorr = pairsAutoCorr.astype(int)
1293
1306
1294 pairsCrossCorr = range(len(pairsList))
1307 pairsCrossCorr = range(len(pairsList))
1295 pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
1308 pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
1296
1309
1297 return pairsAutoCorr, pairsCrossCorr
1310 return pairsAutoCorr, pairsCrossCorr
1298
1311
1299 def techniqueSA(self, pairsSelected, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, lagTRange, correctFactor):
1312 def techniqueSA(self, pairsSelected, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, lagTRange, correctFactor):
1300 """
1313 """
1301 Function that implements Spaced Antenna (SA) technique.
1314 Function that implements Spaced Antenna (SA) technique.
1302
1315
1303 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
1316 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
1304 Direction correction (if necessary), Ranges and SNR
1317 Direction correction (if necessary), Ranges and SNR
1305
1318
1306 Output: Winds estimation (Zonal, Meridional and Vertical)
1319 Output: Winds estimation (Zonal, Meridional and Vertical)
1307
1320
1308 Parameters affected: Winds
1321 Parameters affected: Winds
1309 """
1322 """
1310 #Cross Correlation pairs obtained
1323 #Cross Correlation pairs obtained
1311 pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)
1324 pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)
1312 pairsArray = numpy.array(pairsList)[pairsCrossCorr]
1325 pairsArray = numpy.array(pairsList)[pairsCrossCorr]
1313 pairsSelArray = numpy.array(pairsSelected)
1326 pairsSelArray = numpy.array(pairsSelected)
1314 pairs = []
1327 pairs = []
1315
1328
1316 #Wind estimation pairs obtained
1329 #Wind estimation pairs obtained
1317 for i in range(pairsSelArray.shape[0]/2):
1330 for i in range(pairsSelArray.shape[0]/2):
1318 ind1 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i], axis = 1))[0][0]
1331 ind1 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i], axis = 1))[0][0]
1319 ind2 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i + 1], axis = 1))[0][0]
1332 ind2 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i + 1], axis = 1))[0][0]
1320 pairs.append((ind1,ind2))
1333 pairs.append((ind1,ind2))
1321
1334
1322 indtau = tau.shape[0]/2
1335 indtau = tau.shape[0]/2
1323 tau1 = tau[:indtau,:]
1336 tau1 = tau[:indtau,:]
1324 tau2 = tau[indtau:-1,:]
1337 tau2 = tau[indtau:-1,:]
1325 tau1 = tau1[pairs,:]
1338 tau1 = tau1[pairs,:]
1326 tau2 = tau2[pairs,:]
1339 tau2 = tau2[pairs,:]
1327 phase1 = tau[-1,:]
1340 phase1 = tau[-1,:]
1328
1341
1329 #---------------------------------------------------------------------
1342 #---------------------------------------------------------------------
1330 #Metodo Directo
1343 #Metodo Directo
1331 distx, disty, dist, ang = self.__calculateDistance(position_x, position_y, pairsCrossCorr, pairsList, pairs,azimuth)
1344 distx, disty, dist, ang = self.__calculateDistance(position_x, position_y, pairsCrossCorr, pairsList, pairs,azimuth)
1332 winds = self.__calculateVelHorDir(dist, tau1, tau2, ang)
1345 winds = self.__calculateVelHorDir(dist, tau1, tau2, ang)
1333 winds = stats.nanmean(winds, axis=0)
1346 winds = stats.nanmean(winds, axis=0)
1334 #---------------------------------------------------------------------
1347 #---------------------------------------------------------------------
1335 #Metodo General
1348 #Metodo General
1336 # distx, disty, dist = self.calculateDistance(position_x,position_y,pairsCrossCorr, pairsList, azimuth)
1349 # distx, disty, dist = self.calculateDistance(position_x,position_y,pairsCrossCorr, pairsList, azimuth)
1337 # #Calculo Coeficientes de Funcion de Correlacion
1350 # #Calculo Coeficientes de Funcion de Correlacion
1338 # F,G,A,B,H = self.calculateCoef(tau1,tau2,distx,disty,n)
1351 # F,G,A,B,H = self.calculateCoef(tau1,tau2,distx,disty,n)
1339 # #Calculo de Velocidades
1352 # #Calculo de Velocidades
1340 # winds = self.calculateVelUV(F,G,A,B,H)
1353 # winds = self.calculateVelUV(F,G,A,B,H)
1341
1354
1342 #---------------------------------------------------------------------
1355 #---------------------------------------------------------------------
1343 winds[2,:] = self.__calculateVelVer(phase1, lagTRange, _lambda)
1356 winds[2,:] = self.__calculateVelVer(phase1, lagTRange, _lambda)
1344 winds = correctFactor*winds
1357 winds = correctFactor*winds
1345 return winds
1358 return winds
1346
1359
1347 def __checkTime(self, currentTime, paramInterval, windsInterval):
1360 def __checkTime(self, currentTime, paramInterval, windsInterval):
1348
1361
1349 dataTime = currentTime + paramInterval
1362 dataTime = currentTime + paramInterval
1350 deltaTime = dataTime - self.__initime
1363 deltaTime = dataTime - self.__initime
1351
1364
1352 if deltaTime >= windsInterval or deltaTime < 0:
1365 if deltaTime >= windsInterval or deltaTime < 0:
1353 self.__dataReady = True
1366 self.__dataReady = True
1354 return
1367 return
1355
1368
1356 def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax):
1369 def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax):
1357 '''
1370 '''
1358 Function that implements winds estimation technique with detected meteors.
1371 Function that implements winds estimation technique with detected meteors.
1359
1372
1360 Input: Detected meteors, Minimum meteor quantity to wind estimation
1373 Input: Detected meteors, Minimum meteor quantity to wind estimation
1361
1374
1362 Output: Winds estimation (Zonal and Meridional)
1375 Output: Winds estimation (Zonal and Meridional)
1363
1376
1364 Parameters affected: Winds
1377 Parameters affected: Winds
1365 '''
1378 '''
1366 #Settings
1379 #Settings
1367 nInt = (heightMax - heightMin)/2
1380 nInt = (heightMax - heightMin)/2
1368 winds = numpy.zeros((2,nInt))*numpy.nan
1381 winds = numpy.zeros((2,nInt))*numpy.nan
1369
1382
1370 #Filter errors
1383 #Filter errors
1371 error = numpy.where(arrayMeteor[:,-1] == 0)[0]
1384 error = numpy.where(arrayMeteor[:,-1] == 0)[0]
1372 finalMeteor = arrayMeteor[error,:]
1385 finalMeteor = arrayMeteor[error,:]
1373
1386
1374 #Meteor Histogram
1387 #Meteor Histogram
1375 finalHeights = finalMeteor[:,3]
1388 finalHeights = finalMeteor[:,3]
1376 hist = numpy.histogram(finalHeights, bins = nInt, range = (heightMin,heightMax))
1389 hist = numpy.histogram(finalHeights, bins = nInt, range = (heightMin,heightMax))
1377 nMeteorsPerI = hist[0]
1390 nMeteorsPerI = hist[0]
1378 heightPerI = hist[1]
1391 heightPerI = hist[1]
1379
1392
1380 #Sort of meteors
1393 #Sort of meteors
1381 indSort = finalHeights.argsort()
1394 indSort = finalHeights.argsort()
1382 finalMeteor2 = finalMeteor[indSort,:]
1395 finalMeteor2 = finalMeteor[indSort,:]
1383
1396
1384 # Calculating winds
1397 # Calculating winds
1385 ind1 = 0
1398 ind1 = 0
1386 ind2 = 0
1399 ind2 = 0
1387
1400
1388 for i in range(nInt):
1401 for i in range(nInt):
1389 nMet = nMeteorsPerI[i]
1402 nMet = nMeteorsPerI[i]
1390 ind1 = ind2
1403 ind1 = ind2
1391 ind2 = ind1 + nMet
1404 ind2 = ind1 + nMet
1392
1405
1393 meteorAux = finalMeteor2[ind1:ind2,:]
1406 meteorAux = finalMeteor2[ind1:ind2,:]
1394
1407
1395 if meteorAux.shape[0] >= meteorThresh:
1408 if meteorAux.shape[0] >= meteorThresh:
1396 vel = meteorAux[:, 7]
1409 vel = meteorAux[:, 7]
1397 zen = meteorAux[:, 5]*numpy.pi/180
1410 zen = meteorAux[:, 5]*numpy.pi/180
1398 azim = meteorAux[:, 4]*numpy.pi/180
1411 azim = meteorAux[:, 4]*numpy.pi/180
1399
1412
1400 n = numpy.cos(zen)
1413 n = numpy.cos(zen)
1401 # m = (1 - n**2)/(1 - numpy.tan(azim)**2)
1414 # m = (1 - n**2)/(1 - numpy.tan(azim)**2)
1402 # l = m*numpy.tan(azim)
1415 # l = m*numpy.tan(azim)
1403 l = numpy.sin(zen)*numpy.sin(azim)
1416 l = numpy.sin(zen)*numpy.sin(azim)
1404 m = numpy.sin(zen)*numpy.cos(azim)
1417 m = numpy.sin(zen)*numpy.cos(azim)
1405
1418
1406 A = numpy.vstack((l, m)).transpose()
1419 A = numpy.vstack((l, m)).transpose()
1407 A1 = numpy.dot(numpy.linalg.inv( numpy.dot(A.transpose(),A) ),A.transpose())
1420 A1 = numpy.dot(numpy.linalg.inv( numpy.dot(A.transpose(),A) ),A.transpose())
1408 windsAux = numpy.dot(A1, vel)
1421 windsAux = numpy.dot(A1, vel)
1409
1422
1410 winds[0,i] = windsAux[0]
1423 winds[0,i] = windsAux[0]
1411 winds[1,i] = windsAux[1]
1424 winds[1,i] = windsAux[1]
1412
1425
1413 return winds, heightPerI[:-1]
1426 return winds, heightPerI[:-1]
1414
1427
1415 def run(self, dataOut, technique, **kwargs):
1428 def run(self, dataOut, technique, **kwargs):
1416
1429
1417 param = dataOut.data_param
1430 param = dataOut.data_param
1418 if dataOut.abscissaRange != None:
1431 if dataOut.abscissaRange != None:
1419 absc = dataOut.abscissaRange[:-1]
1432 absc = dataOut.abscissaRange[:-1]
1420 noise = dataOut.noise
1433 noise = dataOut.noise
1421 heightRange = dataOut.getHeiRange()
1434 heightRange = dataOut.getHeiRange()
1422 SNR = dataOut.SNR
1435 SNR = dataOut.SNR
1423
1436
1424 if technique == 'DBS':
1437 if technique == 'DBS':
1425
1438
1426 theta_x = numpy.array(kwargs['dirCosx'])
1439 if kwargs.has_key('dirCosx') and kwargs.has_key('dirCosy'):
1427 theta_y = numpy.array(kwargs['dirCosy'])
1440 theta_x = numpy.array(kwargs['dirCosx'])
1428 azimuth = kwargs['azimuth']
1441 theta_y = numpy.array(kwargs['dirCosy'])
1442 else:
1443 elev = numpy.array(kwargs['elevation'])
1444 azim = numpy.array(kwargs['azimuth'])
1445 theta_x, theta_y = self.__calculateCosDir(elev, azim)
1446 azimuth = kwargs['correctAzimuth']
1429 if kwargs.has_key('horizontalOnly'):
1447 if kwargs.has_key('horizontalOnly'):
1430 horizontalOnly = kwargs['horizontalOnly']
1448 horizontalOnly = kwargs['horizontalOnly']
1431 else: horizontalOnly = False
1449 else: horizontalOnly = False
1432 if kwargs.has_key('correctFactor'):
1450 if kwargs.has_key('correctFactor'):
1433 correctFactor = kwargs['correctFactor']
1451 correctFactor = kwargs['correctFactor']
1434 else: correctFactor = 1
1452 else: correctFactor = 1
1435 if kwargs.has_key('channelList'):
1453 if kwargs.has_key('channelList'):
1436 channelList = kwargs['channelList']
1454 channelList = kwargs['channelList']
1437 if len(channelList) == 2:
1455 if len(channelList) == 2:
1438 horizontalOnly = True
1456 horizontalOnly = True
1439 arrayChannel = numpy.array(channelList)
1457 arrayChannel = numpy.array(channelList)
1440 param = param[arrayChannel,:,:]
1458 param = param[arrayChannel,:,:]
1441 theta_x = theta_x[arrayChannel]
1459 theta_x = theta_x[arrayChannel]
1442 theta_y = theta_y[arrayChannel]
1460 theta_y = theta_y[arrayChannel]
1443
1461
1444 velRadial0 = param[:,1,:] #Radial velocity
1462 velRadial0 = param[:,1,:] #Radial velocity
1445 dataOut.winds, dataOut.heightRange, dataOut.SNR = self.techniqueDBS(velRadial0, theta_x, theta_y, azimuth, correctFactor, horizontalOnly, heightRange, SNR) #DBS Function
1463 dataOut.winds, dataOut.heightRange, dataOut.SNR = self.techniqueDBS(velRadial0, theta_x, theta_y, azimuth, correctFactor, horizontalOnly, heightRange, SNR) #DBS Function
1446 dataOut.initUtcTime = dataOut.ltctime
1464 dataOut.initUtcTime = dataOut.ltctime
1447 dataOut.windsInterval = dataOut.timeInterval
1465 dataOut.windsInterval = dataOut.timeInterval
1448
1466
1449 elif technique == 'SA':
1467 elif technique == 'SA':
1450
1468
1451 #Parameters
1469 #Parameters
1452 position_x = kwargs['positionX']
1470 position_x = kwargs['positionX']
1453 position_y = kwargs['positionY']
1471 position_y = kwargs['positionY']
1454 azimuth = kwargs['azimuth']
1472 azimuth = kwargs['azimuth']
1455
1473
1456 if kwargs.has_key('crosspairsList'):
1474 if kwargs.has_key('crosspairsList'):
1457 pairs = kwargs['crosspairsList']
1475 pairs = kwargs['crosspairsList']
1458 else:
1476 else:
1459 pairs = None
1477 pairs = None
1460
1478
1461 if kwargs.has_key('correctFactor'):
1479 if kwargs.has_key('correctFactor'):
1462 correctFactor = kwargs['correctFactor']
1480 correctFactor = kwargs['correctFactor']
1463 else:
1481 else:
1464 correctFactor = 1
1482 correctFactor = 1
1465
1483
1466 tau = dataOut.data_param
1484 tau = dataOut.data_param
1467 _lambda = dataOut.C/dataOut.frequency
1485 _lambda = dataOut.C/dataOut.frequency
1468 pairsList = dataOut.pairsList
1486 pairsList = dataOut.pairsList
1469 nChannels = dataOut.nChannels
1487 nChannels = dataOut.nChannels
1470
1488
1471 dataOut.winds = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor)
1489 dataOut.winds = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor)
1472 dataOut.initUtcTime = dataOut.ltctime
1490 dataOut.initUtcTime = dataOut.ltctime
1473 dataOut.windsInterval = dataOut.timeInterval
1491 dataOut.windsInterval = dataOut.timeInterval
1474
1492
1475 elif technique == 'Meteors':
1493 elif technique == 'Meteors':
1476 dataOut.flagNoData = True
1494 dataOut.flagNoData = True
1477 self.__dataReady = False
1495 self.__dataReady = False
1478
1496
1479 if kwargs.has_key('nHours'):
1497 if kwargs.has_key('nHours'):
1480 nHours = kwargs['nHours']
1498 nHours = kwargs['nHours']
1481 else:
1499 else:
1482 nHours = 1
1500 nHours = 1
1483
1501
1484 if kwargs.has_key('meteorsPerBin'):
1502 if kwargs.has_key('meteorsPerBin'):
1485 meteorThresh = kwargs['meteorsPerBin']
1503 meteorThresh = kwargs['meteorsPerBin']
1486 else:
1504 else:
1487 meteorThresh = 6
1505 meteorThresh = 6
1488
1506
1489 if kwargs.has_key('hmin'):
1507 if kwargs.has_key('hmin'):
1490 hmin = kwargs['hmin']
1508 hmin = kwargs['hmin']
1491 else: hmin = 70
1509 else: hmin = 70
1492 if kwargs.has_key('hmax'):
1510 if kwargs.has_key('hmax'):
1493 hmax = kwargs['hmax']
1511 hmax = kwargs['hmax']
1494 else: hmax = 110
1512 else: hmax = 110
1495
1513
1496 dataOut.windsInterval = nHours*3600
1514 dataOut.windsInterval = nHours*3600
1497
1515
1498 if self.__isConfig == False:
1516 if self.__isConfig == False:
1499 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
1517 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
1500 #Get Initial LTC time
1518 #Get Initial LTC time
1501 self.__initime = (dataOut.datatime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
1519 self.__initime = (dataOut.datatime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
1502 self.__isConfig = True
1520 self.__isConfig = True
1503
1521
1504 if self.__buffer == None:
1522 if self.__buffer == None:
1505 self.__buffer = dataOut.data_param
1523 self.__buffer = dataOut.data_param
1506 self.__firstdata = copy.copy(dataOut)
1524 self.__firstdata = copy.copy(dataOut)
1507
1525
1508 else:
1526 else:
1509 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
1527 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
1510
1528
1511 self.__checkTime(dataOut.ltctime, dataOut.paramInterval, dataOut.windsInterval) #Check if the buffer is ready
1529 self.__checkTime(dataOut.ltctime, dataOut.paramInterval, dataOut.windsInterval) #Check if the buffer is ready
1512
1530
1513 if self.__dataReady:
1531 if self.__dataReady:
1514 dataOut.initUtcTime = self.__initime
1532 dataOut.initUtcTime = self.__initime
1515 self.__initime = self.__initime + dataOut.windsInterval #to erase time offset
1533 self.__initime = self.__initime + dataOut.windsInterval #to erase time offset
1516
1534
1517 dataOut.winds, dataOut.heightRange = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax)
1535 dataOut.winds, dataOut.heightRange = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax)
1518 dataOut.flagNoData = False
1536 dataOut.flagNoData = False
1519 self.__buffer = None
1537 self.__buffer = None
1520
1538
1521 return No newline at end of file
1539 return
@@ -1,144 +1,144
1 # DIAS 19 Y 20 FEB 2014
1 # DIAS 19 Y 20 FEB 2014
2 # Comprobacion de Resultados DBS con SA
2 # Comprobacion de Resultados DBS con SA
3
3
4 import os, sys
4 import os, sys
5
5
6 path = os.path.split(os.getcwd())[0]
6 path = os.path.split(os.getcwd())[0]
7 sys.path.append(path)
7 sys.path.append(path)
8
8
9 from controller import *
9 from controller import *
10
10
11 desc = "DBS Experiment Test"
11 desc = "DBS Experiment Test"
12 filename = "DBStest.xml"
12 filename = "DBStest.xml"
13
13
14 controllerObj = Project()
14 controllerObj = Project()
15
15
16 controllerObj.setup(id = '191', name='test01', description=desc)
16 controllerObj.setup(id = '191', name='test01', description=desc)
17
17
18 #Experimentos
18 #Experimentos
19
19
20 #2014050 19 Feb 2014
20 #2014050 19 Feb 2014
21 # path = '/home/soporte/Documents/MST_Data/DBS/d2014050'
21 # path = '/home/soporte/Documents/MST_Data/DBS/d2014050'
22 # pathFigure = '/home/soporte/workspace/Graficos/DBS/d2014050p/'
22 # pathFigure = '/home/soporte/workspace/Graficos/DBS/d2014050p/'
23 # xmin = '15.5'
23 # xmin = '15.5'
24 # xmax = '23.99999999'
24 # xmax = '23.99999999'
25 # startTime = '17:25:00'
25 # startTime = '17:25:00'
26 # filehdf5 = "DBS_2014050.hdf5"
26 # filehdf5 = "DBS_2014050.hdf5"
27
27
28 #2014051 20 Feb 2014
28 #2014051 20 Feb 2014
29 path = '/home/soporte/Data/MST/DBS/d2014051'
29 path = '/home/soporte/Data/MST/DBS/d2014051'
30 pathFigure = '/home/soporte/workspace/Graficos/DBS/prueba1/'
30 pathFigure = '/home/soporte/workspace/Graficos/DBS/prueba1/'
31 xmin = '0.0'
31 xmin = '0.0'
32 xmax = '8.0'
32 xmax = '8.0'
33 startTime = '00:00:00'
33 startTime = '00:00:00'
34 filehdf5 = "DBS_2014051.hdf5"
34 filehdf5 = "DBS_2014051.hdf5"
35
35
36
36
37
37
38 #------------------------------------------------------------------------------------------------
38 #------------------------------------------------------------------------------------------------
39 readUnitConfObj = controllerObj.addReadUnit(datatype='VoltageReader',
39 readUnitConfObj = controllerObj.addReadUnit(datatype='VoltageReader',
40 path=path,
40 path=path,
41 startDate='2014/01/31',
41 startDate='2014/01/31',
42 endDate='2014/03/31',
42 endDate='2014/03/31',
43 startTime=startTime,
43 startTime=startTime,
44 endTime='23:59:59',
44 endTime='23:59:59',
45 online=0,
45 online=0,
46 delay=5,
46 delay=5,
47 walk=0)
47 walk=0)
48
48
49 opObj11 = readUnitConfObj.addOperation(name='printNumberOfBlock')
49 opObj11 = readUnitConfObj.addOperation(name='printNumberOfBlock')
50
50
51
51
52 #--------------------------------------------------------------------------------------------------
52 #--------------------------------------------------------------------------------------------------
53
53
54 procUnitConfObj0 = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())
54 procUnitConfObj0 = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())
55
55
56 opObj11 = procUnitConfObj0.addOperation(name='Decoder', optype='other')
56 opObj11 = procUnitConfObj0.addOperation(name='Decoder', optype='other')
57
57
58 opObj11 = procUnitConfObj0.addOperation(name='CohInt', optype='other')
58 opObj11 = procUnitConfObj0.addOperation(name='CohInt', optype='other')
59 opObj11.addParameter(name='n', value='256', format='int')
59 opObj11.addParameter(name='n', value='256', format='int')
60 # opObj11.addParameter(name='n', value='16', format='int')
60 # opObj11.addParameter(name='n', value='16', format='int')
61
61
62 opObj11 = procUnitConfObj0.addOperation(name='selectHeightsByIndex')
62 opObj11 = procUnitConfObj0.addOperation(name='selectHeightsByIndex')
63 opObj11.addParameter(name='minIndex', value='10', format='float')
63 opObj11.addParameter(name='minIndex', value='10', format='float')
64 opObj11.addParameter(name='maxIndex', value='60', format='float')
64 opObj11.addParameter(name='maxIndex', value='60', format='float')
65
65
66 #---------------------------------------------------------------------------------------------------
66 #---------------------------------------------------------------------------------------------------
67
67
68 procUnitConfObj1 = controllerObj.addProcUnit(datatype='SpectraProc', inputId=procUnitConfObj0.getId())
68 procUnitConfObj1 = controllerObj.addProcUnit(datatype='SpectraProc', inputId=procUnitConfObj0.getId())
69 procUnitConfObj1.addParameter(name='nFFTPoints', value='64', format='int')
69 procUnitConfObj1.addParameter(name='nFFTPoints', value='64', format='int')
70 procUnitConfObj1.addParameter(name='nProfiles', value='64', format='int')
70 procUnitConfObj1.addParameter(name='nProfiles', value='64', format='int')
71 # procUnitConfObj1.addParameter(name='ippFactor', value='2', format='int')
71 # procUnitConfObj1.addParameter(name='ippFactor', value='2', format='int')
72 procUnitConfObj1.addParameter(name='pairsList', value='(0,0),(0,1),(2,1)', format='pairsList')
72 procUnitConfObj1.addParameter(name='pairsList', value='(0,0),(0,1),(2,1)', format='pairsList')
73
73
74 opObj11 = procUnitConfObj1.addOperation(name='IncohInt', optype='other')
74 opObj11 = procUnitConfObj1.addOperation(name='IncohInt', optype='other')
75 opObj11.addParameter(name='n', value='5', format='int')
75 opObj11.addParameter(name='n', value='5', format='int')
76
76
77 opObj14 = procUnitConfObj1.addOperation(name='SpectraPlot', optype='other')
77 opObj14 = procUnitConfObj1.addOperation(name='SpectraPlot', optype='other')
78 opObj14.addParameter(name='id', value='1', format='int')
78 opObj14.addParameter(name='id', value='1', format='int')
79 opObj14.addParameter(name='wintitle', value='Con interf', format='str')
79 opObj14.addParameter(name='wintitle', value='Con interf', format='str')
80 opObj14.addParameter(name='save', value='1', format='bool')
80 opObj14.addParameter(name='save', value='1', format='bool')
81 opObj14.addParameter(name='figpath', value=pathFigure, format='str')
81 opObj14.addParameter(name='figpath', value=pathFigure, format='str')
82 opObj14.addParameter(name='zmin', value='5', format='int')
82 opObj14.addParameter(name='zmin', value='5', format='int')
83 opObj14.addParameter(name='zmax', value='90', format='int')
83 opObj14.addParameter(name='zmax', value='90', format='int')
84
84
85 opObj12 = procUnitConfObj1.addOperation(name='removeInterference')
85 opObj12 = procUnitConfObj1.addOperation(name='removeInterference')
86 opObj13 = procUnitConfObj1.addOperation(name='removeDC')
86 opObj13 = procUnitConfObj1.addOperation(name='removeDC')
87 opObj13.addParameter(name='mode', value='1', format='int')
87 opObj13.addParameter(name='mode', value='1', format='int')
88
88
89 opObj12 = procUnitConfObj1.addOperation(name='RTIPlot', optype='other')
89 opObj12 = procUnitConfObj1.addOperation(name='RTIPlot', optype='other')
90 opObj12.addParameter(name='id', value='2', format='int')
90 opObj12.addParameter(name='id', value='2', format='int')
91 opObj12.addParameter(name='wintitle', value='RTI Plot', format='str')
91 opObj12.addParameter(name='wintitle', value='RTI Plot', format='str')
92 opObj12.addParameter(name='save', value='1', format='bool')
92 opObj12.addParameter(name='save', value='1', format='bool')
93 opObj12.addParameter(name='figpath', value = pathFigure, format='str')
93 opObj12.addParameter(name='figpath', value = pathFigure, format='str')
94 opObj12.addParameter(name='xmin', value=xmin, format='float')
94 opObj12.addParameter(name='xmin', value=xmin, format='float')
95 opObj12.addParameter(name='xmax', value=xmax, format='float')
95 opObj12.addParameter(name='xmax', value=xmax, format='float')
96 opObj12.addParameter(name='zmin', value='5', format='int')
96 opObj12.addParameter(name='zmin', value='5', format='int')
97 opObj12.addParameter(name='zmax', value='90', format='int')
97 opObj12.addParameter(name='zmax', value='90', format='int')
98
98
99 #--------------------------------------------------------------------------------------------------
99 #--------------------------------------------------------------------------------------------------
100
100
101 procUnitConfObj2 = controllerObj.addProcUnit(datatype='ParametersProc', inputId=procUnitConfObj1.getId())
101 procUnitConfObj2 = controllerObj.addProcUnit(datatype='ParametersProc', inputId=procUnitConfObj1.getId())
102 opObj20 = procUnitConfObj2.addOperation(name='GetMoments')
102 opObj20 = procUnitConfObj2.addOperation(name='GetMoments')
103
103
104 opObj21 = procUnitConfObj2.addOperation(name='MomentsPlot', optype='other')
104 opObj21 = procUnitConfObj2.addOperation(name='MomentsPlot', optype='other')
105 opObj21.addParameter(name='id', value='3', format='int')
105 opObj21.addParameter(name='id', value='3', format='int')
106 opObj21.addParameter(name='wintitle', value='Moments Plot', format='str')
106 opObj21.addParameter(name='wintitle', value='Moments Plot', format='str')
107 opObj21.addParameter(name='save', value='1', format='bool')
107 opObj21.addParameter(name='save', value='1', format='bool')
108 opObj21.addParameter(name='figpath', value=pathFigure, format='str')
108 opObj21.addParameter(name='figpath', value=pathFigure, format='str')
109 opObj21.addParameter(name='zmin', value='5', format='int')
109 opObj21.addParameter(name='zmin', value='5', format='int')
110 opObj21.addParameter(name='zmax', value='90', format='int')
110 opObj21.addParameter(name='zmax', value='90', format='int')
111
111
112 opObj22 = procUnitConfObj2.addOperation(name='WindProfiler', optype='other')
112 opObj22 = procUnitConfObj2.addOperation(name='WindProfiler', optype='other')
113 opObj22.addParameter(name='technique', value='DBS', format='str')
113 opObj22.addParameter(name='technique', value='DBS', format='str')
114 opObj22.addParameter(name='azimuth', value='51.06', format='float')
114 opObj22.addParameter(name='correctAzimuth', value='51.06', format='float')
115 opObj22.addParameter(name='correctFactor', value='-1', format='float')
115 opObj22.addParameter(name='correctFactor', value='-1', format='float')
116 opObj22.addParameter(name='dirCosx', value='0.041016, 0, -0.054688', format='floatlist')
116 opObj22.addParameter(name='dirCosx', value='0.041016, 0, -0.054688', format='floatlist')
117 opObj22.addParameter(name='dirCosy', value='-0.041016, 0.025391, -0.023438', format='floatlist')
117 opObj22.addParameter(name='dirCosy', value='-0.041016, 0.025391, -0.023438', format='floatlist')
118 # opObj22.addParameter(name='horizontalOnly', value='1', format='bool')
118 # opObj22.addParameter(name='horizontalOnly', value='1', format='bool')
119 # opObj22.addParameter(name='channelList', value='1,2', format='intlist')
119 # opObj22.addParameter(name='channelList', value='1,2', format='intlist')
120
120
121 opObj23 = procUnitConfObj2.addOperation(name='WindProfilerPlot', optype='other')
121 opObj23 = procUnitConfObj2.addOperation(name='WindProfilerPlot', optype='other')
122 opObj23.addParameter(name='id', value='4', format='int')
122 opObj23.addParameter(name='id', value='4', format='int')
123 opObj23.addParameter(name='wintitle', value='Wind Profiler', format='str')
123 opObj23.addParameter(name='wintitle', value='Wind Profiler', format='str')
124 opObj23.addParameter(name='save', value='1', format='bool')
124 opObj23.addParameter(name='save', value='1', format='bool')
125 opObj23.addParameter(name='figpath', value = pathFigure, format='str')
125 opObj23.addParameter(name='figpath', value = pathFigure, format='str')
126 opObj23.addParameter(name='zmin', value='-10', format='int')
126 opObj23.addParameter(name='zmin', value='-10', format='int')
127 opObj23.addParameter(name='zmax', value='10', format='int')
127 opObj23.addParameter(name='zmax', value='10', format='int')
128 opObj23.addParameter(name='zmin_ver', value='-80', format='float')
128 opObj23.addParameter(name='zmin_ver', value='-80', format='float')
129 opObj23.addParameter(name='zmax_ver', value='80', format='float')
129 opObj23.addParameter(name='zmax_ver', value='80', format='float')
130 opObj23.addParameter(name='SNRmin', value='-10', format='int')
130 opObj23.addParameter(name='SNRmin', value='-10', format='int')
131 opObj23.addParameter(name='SNRmax', value='60', format='int')
131 opObj23.addParameter(name='SNRmax', value='60', format='int')
132 opObj23.addParameter(name='SNRthresh', value='0', format='float')
132 opObj23.addParameter(name='SNRthresh', value='0', format='float')
133 opObj23.addParameter(name='xmin', value=xmin, format='float')
133 opObj23.addParameter(name='xmin', value=xmin, format='float')
134 opObj23.addParameter(name='xmax', value=xmax, format='float')
134 opObj23.addParameter(name='xmax', value=xmax, format='float')
135
135
136 #--------------------------------------------------------------------------------------------------
136 #--------------------------------------------------------------------------------------------------
137 print "Escribiendo el archivo XML"
137 print "Escribiendo el archivo XML"
138 controllerObj.writeXml(filename)
138 controllerObj.writeXml(filename)
139 print "Leyendo el archivo XML"
139 print "Leyendo el archivo XML"
140 controllerObj.readXml(filename)
140 controllerObj.readXml(filename)
141
141
142 controllerObj.createObjects()
142 controllerObj.createObjects()
143 controllerObj.connectObjects()
143 controllerObj.connectObjects()
144 controllerObj.run() No newline at end of file
144 controllerObj.run()
@@ -1,131 +1,146
1 import os, sys
1 import os, sys
2
2
3 path = os.path.split(os.getcwd())[0]
3 path = os.path.split(os.getcwd())[0]
4 sys.path.append(path)
4 sys.path.append(path)
5
5
6 from controller import *
6 from controller import *
7
7
8 desc = "AMISR Experiment"
8 desc = "AMISR Experiment"
9
9
10 filename = "amisr_reader.xml"
10 filename = "amisr_reader.xml"
11
11
12 controllerObj = Project()
12 controllerObj = Project()
13
13
14 controllerObj.setup(id = '191', name='test01', description=desc)
14 controllerObj.setup(id = '191', name='test01', description=desc)
15
15
16
16
17 path = os.path.join(os.environ['HOME'],'Development/amisr/data')
17 path = os.path.join(os.environ['HOME'],'Development/amisr/data')
18 path = '/home/soporte/Data/AMISR'
18 path = '/home/soporte/Data/AMISR'
19 figpath = os.path.join(os.environ['HOME'],'Pictures/amisr')
19 figpath = os.path.join(os.environ['HOME'],'Pictures/amisr')
20
20
21 pathFigure = '/home/soporte/workspace/Graficos/DBS/amisr/'
21 pathFigure = '/home/soporte/workspace/Graficos/DBS/amisr/'
22 xmin = '12.0'
22 xmin = '16.0'
23 xmax = '14.0'
23 xmax = '17.0'
24
24
25 readUnitConfObj = controllerObj.addReadUnit(datatype='AMISRReader',
25 readUnitConfObj = controllerObj.addReadUnit(datatype='AMISRReader',
26 path=path,
26 path=path,
27 startDate='2014/10/21',
27 startDate='2014/10/21',
28 endDate='2014/10/21',
28 endDate='2014/10/21',
29 startTime='00:00:00',
29 startTime='00:00:00',
30 endTime='23:59:59',
30 endTime='23:59:59',
31 walk=1,
31 walk=1,
32 timezone='lt')
32 timezone='lt')
33
33
34 #AMISR Processing Unit
34 #AMISR Processing Unit
35 procUnitAMISRBeam0 = controllerObj.addProcUnit(datatype='AMISRProc', inputId=readUnitConfObj.getId())
35 procUnitAMISRBeam0 = controllerObj.addProcUnit(datatype='AMISRProc', inputId=readUnitConfObj.getId())
36
36
37 opObj11 = procUnitAMISRBeam0.addOperation(name='PrintInfo', optype='other')
37 opObj11 = procUnitAMISRBeam0.addOperation(name='PrintInfo', optype='other')
38
38
39 #Reshaper
39 #Reshaper
40 opObj11 = procUnitAMISRBeam0.addOperation(name='ProfileToChannels', optype='other')
40 opObj11 = procUnitAMISRBeam0.addOperation(name='ProfileToChannels', optype='other')
41
41
42
42
43 #Beam Selector
43 #Beam Selector
44 #opObj11 = procUnitAMISRBeam0.addOperation(name='BeamSelector', optype='other')
44 #opObj11 = procUnitAMISRBeam0.addOperation(name='BeamSelector', optype='other')
45 #opObj11.addParameter(name='beam', value='0', format='int')
45 #opObj11.addParameter(name='beam', value='0', format='int')
46
46
47 #Voltage Processing Unit
47 #Voltage Processing Unit
48 procUnitConfObjBeam0 = controllerObj.addProcUnit(datatype='VoltageProc', inputId=procUnitAMISRBeam0.getId())
48 procUnitConfObjBeam0 = controllerObj.addProcUnit(datatype='VoltageProc', inputId=procUnitAMISRBeam0.getId())
49 opObj10 = procUnitConfObjBeam0.addOperation(name='setRadarFrequency')
50 opObj10.addParameter(name='frequency', value='445e6', format='float')
51
52 opObj12 = procUnitConfObjBeam0.addOperation(name='selectHeights')
53 opObj12.addParameter(name='minHei', value='0', format='float')
54 opObj12.addParameter(name='maxHei', value='10', format='float')
49 #Coherent Integration
55 #Coherent Integration
50 opObj11 = procUnitConfObjBeam0.addOperation(name='CohInt', optype='other')
56 opObj11 = procUnitConfObjBeam0.addOperation(name='CohInt', optype='other')
51 opObj11.addParameter(name='n', value='8', format='int')
57 opObj11.addParameter(name='n', value='8', format='int')
52 #Spectra Unit Processing, getting spectras with nProfiles and nFFTPoints
58 #Spectra Unit Processing, getting spectras with nProfiles and nFFTPoints
53 procUnitConfObjSpectraBeam0 = controllerObj.addProcUnit(datatype='SpectraProc', inputId=procUnitConfObjBeam0.getId())
59 procUnitConfObjSpectraBeam0 = controllerObj.addProcUnit(datatype='SpectraProc', inputId=procUnitConfObjBeam0.getId())
54 procUnitConfObjSpectraBeam0.addParameter(name='nFFTPoints', value=32, format='int')
60 procUnitConfObjSpectraBeam0.addParameter(name='nFFTPoints', value=32, format='int')
55 procUnitConfObjSpectraBeam0.addParameter(name='nProfiles', value=32, format='int')
61 procUnitConfObjSpectraBeam0.addParameter(name='nProfiles', value=32, format='int')
62
63 opObj11 = procUnitConfObjSpectraBeam0.addOperation(name='IncohInt', optype='other')
64 opObj11.addParameter(name='n', value='16', format='int')
65
56 #RemoveDc
66 #RemoveDc
57 opObj11 = procUnitConfObjSpectraBeam0.addOperation(name='removeDC')
67 opObj11 = procUnitConfObjSpectraBeam0.addOperation(name='removeDC')
58
68
59 #Noise Estimation
69 #Noise Estimation
60 opObj11 = procUnitConfObjSpectraBeam0.addOperation(name='getNoise')
70 opObj11 = procUnitConfObjSpectraBeam0.addOperation(name='getNoise')
61 opObj11.addParameter(name='minHei', value='5', format='float')
71 opObj11.addParameter(name='minHei', value='5', format='float')
62 opObj11.addParameter(name='maxHei', value='20', format='float')
72 opObj11.addParameter(name='maxHei', value='20', format='float')
63
73
64 #SpectraPlot
74 #SpectraPlot
65 opObj11 = procUnitConfObjSpectraBeam0.addOperation(name='SpectraPlot', optype='other')
75 opObj11 = procUnitConfObjSpectraBeam0.addOperation(name='SpectraPlot', optype='other')
66 opObj11.addParameter(name='id', value='100', format='int')
76 opObj11.addParameter(name='id', value='100', format='int')
67 opObj11.addParameter(name='wintitle', value='AMISR Beam 0', format='str')
77 opObj11.addParameter(name='wintitle', value='AMISR Beam 0', format='str')
68 opObj11.addParameter(name='zmin', value='30', format='int')
78 opObj11.addParameter(name='zmin', value='30', format='int')
69 opObj11.addParameter(name='zmax', value='80', format='int')
79 opObj11.addParameter(name='zmax', value='80', format='int')
70
80 opObj11.addParameter(name='save', value='1', format='bool')
81 opObj11.addParameter(name='figpath', value = pathFigure, format='str')
71 #RTIPlot
82 #RTIPlot
72 #title0 = 'RTI AMISR Beam 0'
83 #title0 = 'RTI AMISR Beam 0'
73 #opObj11 = procUnitConfObjSpectraBeam0.addOperation(name='RTIPlot', optype='other')
84 opObj11 = procUnitConfObjSpectraBeam0.addOperation(name='RTIPlot', optype='other')
74 #opObj11.addParameter(name='id', value='200', format='int')
85 opObj11.addParameter(name='id', value='200', format='int')
75 #opObj11.addParameter(name='wintitle', value=title0, format='str')
86 # opObj11.addParameter(name='wintitle', value=title0, format='str')
76 #opObj11.addParameter(name='showprofile', value='0', format='int')
87 opObj11.addParameter(name='showprofile', value='0', format='int')
77 ##Setting RTI time using xmin,xmax
88 #Setting RTI time using xmin,xmax
78 #opObj11.addParameter(name='xmin', value='15', format='int')
89 opObj11.addParameter(name='xmin', value=xmin, format='float')
79 #opObj11.addParameter(name='xmax', value='23', format='int')
90 opObj11.addParameter(name='xmax', value=xmax, format='float')
80 #Setting dB range with zmin, zmax
91 # Setting dB range with zmin, zmax
81 #opObj11.addParameter(name='zmin', value='45', format='int')
92 opObj11.addParameter(name='zmin', value='45', format='int')
82 #opObj11.addParameter(name='zmax', value='70', format='int')
93 opObj11.addParameter(name='zmax', value='70', format='int')
83 #Save RTI
94 # Save RTI
84 #figfile0 = 'amisr_rti_beam0.png'
95 # figfile0 = 'amisr_rti_beam0.png'
85 #opObj11.addParameter(name='figpath', value=figpath, format='str')
96 # opObj11.addParameter(name='figpath', value=figpath, format='str')
86 #opObj11.addParameter(name='figfile', value=figfile0, format='str')
97 # opObj11.addParameter(name='figfile', value=figfile0, format='str')
87
98
88 #-----------------------------------------------------------------------------------------------
99 #-----------------------------------------------------------------------------------------------
89 procUnitConfObj2 = controllerObj.addProcUnit(datatype='ParametersProc', inputId=procUnitConfObjSpectraBeam0 .getId())
100 procUnitConfObj2 = controllerObj.addProcUnit(datatype='ParametersProc', inputId=procUnitConfObjSpectraBeam0 .getId())
90 opObj20 = procUnitConfObj2.addOperation(name='GetMoments')
101 opObj20 = procUnitConfObj2.addOperation(name='GetMoments')
91
102
92 # opObj21 = procUnitConfObj2.addOperation(name='MomentsPlot', optype='other')
103 opObj21 = procUnitConfObj2.addOperation(name='MomentsPlot', optype='other')
93 # opObj21.addParameter(name='id', value='3', format='int')
104 opObj21.addParameter(name='id', value='3', format='int')
94 # opObj21.addParameter(name='wintitle', value='Moments Plot', format='str')
105 opObj21.addParameter(name='wintitle', value='Moments Plot', format='str')
95 # opObj21.addParameter(name='save', value='1', format='bool')
106 opObj21.addParameter(name='save', value='1', format='bool')
96 # opObj21.addParameter(name='figpath', value=pathFigure, format='str')
107 opObj21.addParameter(name='figpath', value=pathFigure, format='str')
97 # opObj21.addParameter(name='zmin', value='5', format='int')
108 opObj21.addParameter(name='zmin', value='30', format='int')
98 # opObj21.addParameter(name='zmax', value='90', format='int')
109 opObj21.addParameter(name='zmax', value='80', format='int')
99
110
100 opObj22 = procUnitConfObj2.addOperation(name='WindProfiler', optype='other')
111 opObj22 = procUnitConfObj2.addOperation(name='WindProfiler', optype='other')
101 opObj22.addParameter(name='technique', value='DBS', format='str')
112 opObj22.addParameter(name='technique', value='DBS', format='str')
102 opObj22.addParameter(name='azimuth', value='51.06', format='float')
113 opObj22.addParameter(name='correctAzimuth', value='51.06', format='float')
103 opObj22.addParameter(name='correctFactor', value='-1', format='float')
114 opObj22.addParameter(name='correctFactor', value='-1', format='float')
104 opObj22.addParameter(name='dirCosx', value='9.63770247e-01, 5.93066547e-17, 1, 5.93066547e-17,-9.68583161e-01', format='floatlist')
115 opObj22.addParameter(name='azimuth', value='0,180,90,0,-90', format='floatlist')
105 opObj22.addParameter(name='dirCosy', value=' 0,-0.96858316,0,0.96858316,0', format='floatlist')
116 opObj22.addParameter(name='elevation', value='74.53,75.60,75.60,90.0,75.60', format='floatlist')
106 # opObj22.addParameter(name='horizontalOnly', value='1', format='bool')
117 # opObj22.addParameter(name='horizontalOnly', value='1', format='bool')
107 # opObj22.addParameter(name='channelList', value='1,2', format='intlist')
118 # opObj22.addParameter(name='channelList', value='1,2', format='intlist')
108
119
109 opObj23 = procUnitConfObj2.addOperation(name='WindProfilerPlot', optype='other')
120 opObj23 = procUnitConfObj2.addOperation(name='WindProfilerPlot', optype='other')
110 opObj23.addParameter(name='id', value='4', format='int')
121 opObj23.addParameter(name='id', value='4', format='int')
111 opObj23.addParameter(name='wintitle', value='Wind Profiler', format='str')
122 opObj23.addParameter(name='wintitle', value='Wind Profiler', format='str')
112 opObj23.addParameter(name='save', value='1', format='bool')
123 opObj23.addParameter(name='save', value='1', format='bool')
113 opObj23.addParameter(name='figpath', value = pathFigure, format='str')
124 opObj23.addParameter(name='figpath', value = pathFigure, format='str')
114 opObj23.addParameter(name='zmin', value='-10', format='int')
125 opObj23.addParameter(name='zmin', value='-20', format='int')
115 opObj23.addParameter(name='zmax', value='10', format='int')
126 opObj23.addParameter(name='zmax', value='20', format='int')
116 opObj23.addParameter(name='zmin_ver', value='-80', format='float')
127 opObj23.addParameter(name='zmin_ver', value='-100', format='float')
117 opObj23.addParameter(name='zmax_ver', value='80', format='float')
128 opObj23.addParameter(name='zmax_ver', value='100', format='float')
118 opObj23.addParameter(name='SNRmin', value='-10', format='int')
129 opObj23.addParameter(name='SNRmin', value='-10', format='int')
119 opObj23.addParameter(name='SNRmax', value='60', format='int')
130 opObj23.addParameter(name='SNRmax', value='60', format='int')
120 opObj23.addParameter(name='SNRthresh', value='0', format='float')
131 opObj23.addParameter(name='SNRthresh', value='0', format='float')
121 opObj23.addParameter(name='xmin', value=xmin, format='float')
132 opObj23.addParameter(name='xmin', value=xmin, format='float')
122 opObj23.addParameter(name='xmax', value=xmax, format='float')
133 opObj23.addParameter(name='xmax', value=xmax, format='float')
123
134
124 print "Escribiendo el archivo XML"
135 print "Escribiendo el archivo XML"
125 controllerObj.writeXml(filename)
136 controllerObj.writeXml(filename)
126 print "Leyendo el archivo XML"
137 print "Leyendo el archivo XML"
127 controllerObj.readXml(filename)
138 controllerObj.readXml(filename)
128
139
129 controllerObj.createObjects()
140 controllerObj.createObjects()
130 controllerObj.connectObjects()
141 controllerObj.connectObjects()
131 controllerObj.run()
142 controllerObj.run()
143
144 #21 3 pm
145
146
General Comments 0
You need to be logged in to leave comments. Login now