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