##// END OF EJS Templates
schain
RSS
fix pulse pair spectra
avaldez -
r1283:47c7b4917f5b
parent child
Show More
Side by Side Unified
Context file:
r1283:47c7b4917f5b Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,1809 +1,1809
1 1 '''
2 2 Created on Jul 9, 2014
3 3
4 4 @author: roj-idl71
5 5 '''
6 6 import os
7 7 import datetime
8 8 import numpy
9 9
10 10 from .figure import Figure, isRealtime, isTimeInHourRange
11 11 from .plotting_codes import *
12 12 from schainpy.model.proc.jroproc_base import MPDecorator
13 13
14 14 from schainpy.utils import log
15 15
16 16 @MPDecorator
17 17 class SpectraPlot_(Figure):
18 18
19 19 isConfig = None
20 20 __nsubplots = None
21 21
22 22 WIDTHPROF = None
23 23 HEIGHTPROF = None
24 24 PREFIX = 'spc'
25 25
26 26 def __init__(self):
27 27 Figure.__init__(self)
28 28 self.isConfig = False
29 29 self.__nsubplots = 1
30 30 self.WIDTH = 250
31 31 self.HEIGHT = 250
32 32 self.WIDTHPROF = 120
33 33 self.HEIGHTPROF = 0
34 34 self.counter_imagwr = 0
35 35
36 36 self.PLOT_CODE = SPEC_CODE
37 37
38 38 self.FTP_WEI = None
39 39 self.EXP_CODE = None
40 40 self.SUB_EXP_CODE = None
41 41 self.PLOT_POS = None
42 42
43 43 self.__xfilter_ena = False
44 44 self.__yfilter_ena = False
45 45
46 46 self.indice=1
47 47
48 48 def getSubplots(self):
49 49
50 50 ncol = int(numpy.sqrt(self.nplots)+0.9)
51 51 nrow = int(self.nplots*1./ncol + 0.9)
52 52
53 53 return nrow, ncol
54 54
55 55 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
56 56
57 57 self.__showprofile = showprofile
58 58 self.nplots = nplots
59 59
60 60 ncolspan = 1
61 61 colspan = 1
62 62 if showprofile:
63 63 ncolspan = 3
64 64 colspan = 2
65 65 self.__nsubplots = 2
66 66
67 67 self.createFigure(id = id,
68 68 wintitle = wintitle,
69 69 widthplot = self.WIDTH + self.WIDTHPROF,
70 70 heightplot = self.HEIGHT + self.HEIGHTPROF,
71 71 show=show)
72 72
73 73 nrow, ncol = self.getSubplots()
74 74
75 75 counter = 0
76 76 for y in range(nrow):
77 77 for x in range(ncol):
78 78
79 79 if counter >= self.nplots:
80 80 break
81 81
82 82 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
83 83
84 84 if showprofile:
85 85 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
86 86
87 87 counter += 1
88 88
89 89 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=True,
90 90 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
91 91 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
92 92 server=None, folder=None, username=None, password=None,
93 93 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False,
94 94 xaxis="frequency", colormap='jet', normFactor=None):
95 95
96 96 """
97 97
98 98 Input:
99 99 dataOut :
100 100 id :
101 101 wintitle :
102 102 channelList :
103 103 showProfile :
104 104 xmin : None,
105 105 xmax : None,
106 106 ymin : None,
107 107 ymax : None,
108 108 zmin : None,
109 109 zmax : None
110 110 """
111 111 if dataOut.flagNoData:
112 112 return dataOut
113 113
114 114 if realtime:
115 115 if not(isRealtime(utcdatatime = dataOut.utctime)):
116 116 print('Skipping this plot function')
117 117 return
118 118
119 119 if channelList == None:
120 120 channelIndexList = dataOut.channelIndexList
121 121 else:
122 122 channelIndexList = []
123 123 for channel in channelList:
124 124 if channel not in dataOut.channelList:
125 125 raise ValueError("Channel %d is not in dataOut.channelList" %channel)
126 126 channelIndexList.append(dataOut.channelList.index(channel))
127 127
128 128 if normFactor is None:
129 129 factor = dataOut.normFactor
130 130 else:
131 131 factor = normFactor
132 132 if xaxis == "frequency":
133 133 x = dataOut.getFreqRange(1)/1000.
134 134 xlabel = "Frequency (kHz)"
135 135
136 136 elif xaxis == "time":
137 137 x = dataOut.getAcfRange(1)
138 138 xlabel = "Time (ms)"
139 139
140 140 else:
141 141 x = dataOut.getVelRange(1)
142 142 xlabel = "Velocity (m/s)"
143 143
144 144 ylabel = "Range (km)"
145 145
146 146 y = dataOut.getHeiRange()
147 147 z = dataOut.data_spc/factor
148 148 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
149 149 zdB = 10*numpy.log10(z)
150 150
151 151 avg = numpy.average(z, axis=1)
152 152 avgdB = 10*numpy.log10(avg)
153 153
154 154 noise = dataOut.getNoise()/factor
155 155 noisedB = 10*numpy.log10(noise)
156 156
157 157 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
158 158 title = wintitle + " Spectra"
159 159
160 160 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
161 161 title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith)
162 162
163 163 if not self.isConfig:
164 164
165 165 nplots = len(channelIndexList)
166 166
167 167 self.setup(id=id,
168 168 nplots=nplots,
169 169 wintitle=wintitle,
170 170 showprofile=showprofile,
171 171 show=show)
172 172
173 173 if xmin == None: xmin = numpy.nanmin(x)
174 174 if xmax == None: xmax = numpy.nanmax(x)
175 175 if ymin == None: ymin = numpy.nanmin(y)
176 176 if ymax == None: ymax = numpy.nanmax(y)
177 177 if zmin == None: zmin = numpy.floor(numpy.nanmin(noisedB)) - 3
178 178 if zmax == None: zmax = numpy.ceil(numpy.nanmax(avgdB)) + 3
179 179
180 180 self.FTP_WEI = ftp_wei
181 181 self.EXP_CODE = exp_code
182 182 self.SUB_EXP_CODE = sub_exp_code
183 183 self.PLOT_POS = plot_pos
184 184
185 185 self.isConfig = True
186 186
187 187 self.setWinTitle(title)
188 188
189 189 for i in range(self.nplots):
190 190 index = channelIndexList[i]
191 191 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
192 192 title = "Channel %d: %4.2fdB: %s" %(dataOut.channelList[index], noisedB[index], str_datetime)
193 193 if len(dataOut.beam.codeList) != 0:
194 194 title = "Ch%d:%4.2fdB,%2.2f,%2.2f:%s" %(dataOut.channelList[index], noisedB[index], dataOut.beam.azimuthList[index], dataOut.beam.zenithList[index], str_datetime)
195 195
196 196 axes = self.axesList[i*self.__nsubplots]
197 197 axes.pcolor(x, y, zdB[index,:,:],
198 198 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
199 199 xlabel=xlabel, ylabel=ylabel, title=title, colormap=colormap,
200 200 ticksize=9, cblabel='')
201 201
202 202 if self.__showprofile:
203 203 axes = self.axesList[i*self.__nsubplots +1]
204 204 axes.pline(avgdB[index,:], y,
205 205 xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
206 206 xlabel='dB', ylabel='', title='',
207 207 ytick_visible=False,
208 208 grid='x')
209 209
210 210 noiseline = numpy.repeat(noisedB[index], len(y))
211 211 axes.addpline(noiseline, y, idline=1, color="black", linestyle="dashed", lw=2)
212 212
213 213 self.draw()
214 214
215 215 if figfile == None:
216 216 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
217 217 name = str_datetime
218 218 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
219 219 name = name + '_az' + '_%2.2f'%(dataOut.azimuth) + '_zn' + '_%2.2f'%(dataOut.zenith)
220 220 figfile = self.getFilename(name)
221 221
222 222 self.save(figpath=figpath,
223 223 figfile=figfile,
224 224 save=save,
225 225 ftp=ftp,
226 226 wr_period=wr_period,
227 227 thisDatetime=thisDatetime)
228 228
229 229
230 230 return dataOut
231 231
232 232 @MPDecorator
233 233 class WpowerPlot_(Figure):
234 234
235 235 isConfig = None
236 236 __nsubplots = None
237 237
238 238 WIDTHPROF = None
239 239 HEIGHTPROF = None
240 240 PREFIX = 'wpo'
241 241
242 242 def __init__(self):
243 243 Figure.__init__(self)
244 244 self.isConfig = False
245 245 self.__nsubplots = 1
246 246 self.WIDTH = 250
247 247 self.HEIGHT = 250
248 248 self.WIDTHPROF = 120
249 249 self.HEIGHTPROF = 0
250 250 self.counter_imagwr = 0
251 251
252 252 self.PLOT_CODE = WPO_CODE
253 253
254 254 self.FTP_WEI = None
255 255 self.EXP_CODE = None
256 256 self.SUB_EXP_CODE = None
257 257 self.PLOT_POS = None
258 258
259 259 self.__xfilter_ena = False
260 260 self.__yfilter_ena = False
261 261
262 262 self.indice=1
263 263
264 264 def getSubplots(self):
265 265
266 266 ncol = int(numpy.sqrt(self.nplots)+0.9)
267 267 nrow = int(self.nplots*1./ncol + 0.9)
268 268
269 269 return nrow, ncol
270 270
271 271 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
272 272
273 273 self.__showprofile = showprofile
274 274 self.nplots = nplots
275 275
276 276 ncolspan = 1
277 277 colspan = 1
278 278 if showprofile:
279 279 ncolspan = 3
280 280 colspan = 2
281 281 self.__nsubplots = 2
282 282
283 283 self.createFigure(id = id,
284 284 wintitle = wintitle,
285 285 widthplot = self.WIDTH + self.WIDTHPROF,
286 286 heightplot = self.HEIGHT + self.HEIGHTPROF,
287 287 show=show)
288 288
289 289 nrow, ncol = self.getSubplots()
290 290
291 291 counter = 0
292 292 for y in range(nrow):
293 293 for x in range(ncol):
294 294
295 295 if counter >= self.nplots:
296 296 break
297 297
298 298 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
299 299
300 300 if showprofile:
301 301 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
302 302
303 303 counter += 1
304 304
305 305 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=True,
306 306 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
307 307 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
308 308 server=None, folder=None, username=None, password=None,
309 309 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False,
310 310 xaxis="frequency", colormap='jet', normFactor=None):
311 311
312 312 """
313 313
314 314 Input:
315 315 dataOut :
316 316 id :
317 317 wintitle :
318 318 channelList :
319 319 showProfile :
320 320 xmin : None,
321 321 xmax : None,
322 322 ymin : None,
323 323 ymax : None,
324 324 zmin : None,
325 325 zmax : None
326 326 """
327 print("***************PLOTEO******************")
328 print("DATAOUT SHAPE : ",dataOut.data.shape)
327 #print("***************PLOTEO******************")
328 #print("DATAOUT SHAPE : ",dataOut.data.shape)
329 329 if dataOut.flagNoData:
330 330 return dataOut
331 331
332 332 if realtime:
333 333 if not(isRealtime(utcdatatime = dataOut.utctime)):
334 334 print('Skipping this plot function')
335 335 return
336 336
337 337 if channelList == None:
338 338 channelIndexList = dataOut.channelIndexList
339 339 else:
340 340 channelIndexList = []
341 341 for channel in channelList:
342 342 if channel not in dataOut.channelList:
343 343 raise ValueError("Channel %d is not in dataOut.channelList" %channel)
344 344 channelIndexList.append(dataOut.channelList.index(channel))
345 345
346 346
347 print("channelIndexList",channelIndexList)
347 #print("channelIndexList",channelIndexList)
348 348 if normFactor is None:
349 349 factor = dataOut.normFactor
350 350 else:
351 351 factor = normFactor
352 352 if xaxis == "frequency":
353 353 x = dataOut.getFreqRange(1)/1000.
354 354 xlabel = "Frequency (kHz)"
355 355
356 356 elif xaxis == "time":
357 357 x = dataOut.getAcfRange(1)
358 358 xlabel = "Time (ms)"
359 359
360 360 else:
361 361 x = dataOut.getVelRange(1)
362 362 xlabel = "Velocity (m/s)"
363 363
364 364 ylabel = "Range (km)"
365 365
366 366 y = dataOut.getHeiRange()
367 print("factor",factor)
367 #print("factor",factor)
368 368
369 369 z = dataOut.data/factor # dividido /factor
370 370 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
371 371 zdB = 10*numpy.log10(z)
372 372
373 373 avg = numpy.average(z, axis=1)
374 374 avgdB = 10*numpy.log10(avg)
375 375
376 376 noise = dataOut.getNoise()/factor
377 377 noisedB = 10*numpy.log10(noise)
378 378
379 379 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
380 380 title = wintitle + "Weather Power"
381 381
382 382 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
383 383 title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith)
384 384
385 385 if not self.isConfig:
386 386
387 387 nplots = len(channelIndexList)
388 388
389 389 self.setup(id=id,
390 390 nplots=nplots,
391 391 wintitle=wintitle,
392 392 showprofile=showprofile,
393 393 show=show)
394 394
395 395 if xmin == None: xmin = numpy.nanmin(x)
396 396 if xmax == None: xmax = numpy.nanmax(x)
397 397 if ymin == None: ymin = numpy.nanmin(y)
398 398 if ymax == None: ymax = numpy.nanmax(y)
399 399 if zmin == None: zmin = numpy.floor(numpy.nanmin(noisedB)) - 3
400 400 if zmax == None: zmax = numpy.ceil(numpy.nanmax(avgdB)) + 3
401 401
402 402 self.FTP_WEI = ftp_wei
403 403 self.EXP_CODE = exp_code
404 404 self.SUB_EXP_CODE = sub_exp_code
405 405 self.PLOT_POS = plot_pos
406 406
407 407 self.isConfig = True
408 408
409 409 self.setWinTitle(title)
410 410
411 411 for i in range(self.nplots):
412 412 index = channelIndexList[i]
413 413 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
414 414 title = "Channel %d: %4.2fdB: %s" %(dataOut.channelList[index], noisedB[index], str_datetime)
415 415 if len(dataOut.beam.codeList) != 0:
416 416 title = "Ch%d:%4.2fdB,%2.2f,%2.2f:%s" %(dataOut.channelList[index], noisedB[index], dataOut.beam.azimuthList[index], dataOut.beam.zenithList[index], str_datetime)
417 417
418 418 axes = self.axesList[i*self.__nsubplots]
419 419 axes.pcolor(x, y, zdB[index,:,:],
420 420 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
421 421 xlabel=xlabel, ylabel=ylabel, title=title, colormap=colormap,
422 422 ticksize=9, cblabel='')
423 423
424 424 if self.__showprofile:
425 425 axes = self.axesList[i*self.__nsubplots +1]
426 426 axes.pline(avgdB[index,:], y,
427 427 xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
428 428 xlabel='dB', ylabel='', title='',
429 429 ytick_visible=False,
430 430 grid='x')
431 431
432 432 noiseline = numpy.repeat(noisedB[index], len(y))
433 433 axes.addpline(noiseline, y, idline=1, color="black", linestyle="dashed", lw=2)
434 434
435 435 self.draw()
436 436
437 437 if figfile == None:
438 438 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
439 439 name = str_datetime
440 440 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
441 441 name = name + '_az' + '_%2.2f'%(dataOut.azimuth) + '_zn' + '_%2.2f'%(dataOut.zenith)
442 442 figfile = self.getFilename(name)
443 443
444 444 self.save(figpath=figpath,
445 445 figfile=figfile,
446 446 save=save,
447 447 ftp=ftp,
448 448 wr_period=wr_period,
449 449 thisDatetime=thisDatetime)
450 450 return dataOut
451 451
452 452 @MPDecorator
453 453 class CrossSpectraPlot_(Figure):
454 454
455 455 isConfig = None
456 456 __nsubplots = None
457 457
458 458 WIDTH = None
459 459 HEIGHT = None
460 460 WIDTHPROF = None
461 461 HEIGHTPROF = None
462 462 PREFIX = 'cspc'
463 463
464 464 def __init__(self):
465 465 Figure.__init__(self)
466 466 self.isConfig = False
467 467 self.__nsubplots = 4
468 468 self.counter_imagwr = 0
469 469 self.WIDTH = 250
470 470 self.HEIGHT = 250
471 471 self.WIDTHPROF = 0
472 472 self.HEIGHTPROF = 0
473 473
474 474 self.PLOT_CODE = CROSS_CODE
475 475 self.FTP_WEI = None
476 476 self.EXP_CODE = None
477 477 self.SUB_EXP_CODE = None
478 478 self.PLOT_POS = None
479 479
480 480 self.indice=0
481 481
482 482 def getSubplots(self):
483 483
484 484 ncol = 4
485 485 nrow = self.nplots
486 486
487 487 return nrow, ncol
488 488
489 489 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
490 490
491 491 self.__showprofile = showprofile
492 492 self.nplots = nplots
493 493
494 494 ncolspan = 1
495 495 colspan = 1
496 496
497 497 self.createFigure(id = id,
498 498 wintitle = wintitle,
499 499 widthplot = self.WIDTH + self.WIDTHPROF,
500 500 heightplot = self.HEIGHT + self.HEIGHTPROF,
501 501 show=True)
502 502
503 503 nrow, ncol = self.getSubplots()
504 504
505 505 counter = 0
506 506 for y in range(nrow):
507 507 for x in range(ncol):
508 508 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
509 509
510 510 counter += 1
511 511
512 512 def run(self, dataOut, id, wintitle="", pairsList=None,
513 513 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
514 514 coh_min=None, coh_max=None, phase_min=None, phase_max=None,
515 515 save=False, figpath='./', figfile=None, ftp=False, wr_period=1,
516 516 power_cmap='jet', coherence_cmap='jet', phase_cmap='RdBu_r', show=True,
517 517 server=None, folder=None, username=None, password=None,
518 518 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, normFactor=None,
519 519 xaxis='frequency'):
520 520
521 521 """
522 522
523 523 Input:
524 524 dataOut :
525 525 id :
526 526 wintitle :
527 527 channelList :
528 528 showProfile :
529 529 xmin : None,
530 530 xmax : None,
531 531 ymin : None,
532 532 ymax : None,
533 533 zmin : None,
534 534 zmax : None
535 535 """
536 536
537 537 if dataOut.flagNoData:
538 538 return dataOut
539 539
540 540 if pairsList == None:
541 541 pairsIndexList = dataOut.pairsIndexList
542 542 else:
543 543 pairsIndexList = []
544 544 for pair in pairsList:
545 545 if pair not in dataOut.pairsList:
546 546 raise ValueError("Pair %s is not in dataOut.pairsList" %str(pair))
547 547 pairsIndexList.append(dataOut.pairsList.index(pair))
548 548
549 549 if not pairsIndexList:
550 550 return
551 551
552 552 if len(pairsIndexList) > 4:
553 553 pairsIndexList = pairsIndexList[0:4]
554 554
555 555 if normFactor is None:
556 556 factor = dataOut.normFactor
557 557 else:
558 558 factor = normFactor
559 559 x = dataOut.getVelRange(1)
560 560 y = dataOut.getHeiRange()
561 561 z = dataOut.data_spc[:,:,:]/factor
562 562 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
563 563
564 564 noise = dataOut.noise/factor
565 565
566 566 zdB = 10*numpy.log10(z)
567 567 noisedB = 10*numpy.log10(noise)
568 568
569 569 if coh_min == None:
570 570 coh_min = 0.0
571 571 if coh_max == None:
572 572 coh_max = 1.0
573 573
574 574 if phase_min == None:
575 575 phase_min = -180
576 576 if phase_max == None:
577 577 phase_max = 180
578 578
579 579 #thisDatetime = dataOut.datatime
580 580 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
581 581 title = wintitle + " Cross-Spectra: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
582 582 # xlabel = "Velocity (m/s)"
583 583 ylabel = "Range (Km)"
584 584
585 585 if xaxis == "frequency":
586 586 x = dataOut.getFreqRange(1)/1000.
587 587 xlabel = "Frequency (kHz)"
588 588
589 589 elif xaxis == "time":
590 590 x = dataOut.getAcfRange(1)
591 591 xlabel = "Time (ms)"
592 592
593 593 else:
594 594 x = dataOut.getVelRange(1)
595 595 xlabel = "Velocity (m/s)"
596 596
597 597 if not self.isConfig:
598 598
599 599 nplots = len(pairsIndexList)
600 600
601 601 self.setup(id=id,
602 602 nplots=nplots,
603 603 wintitle=wintitle,
604 604 showprofile=False,
605 605 show=show)
606 606
607 607 avg = numpy.abs(numpy.average(z, axis=1))
608 608 avgdB = 10*numpy.log10(avg)
609 609
610 610 if xmin == None: xmin = numpy.nanmin(x)
611 611 if xmax == None: xmax = numpy.nanmax(x)
612 612 if ymin == None: ymin = numpy.nanmin(y)
613 613 if ymax == None: ymax = numpy.nanmax(y)
614 614 if zmin == None: zmin = numpy.floor(numpy.nanmin(noisedB)) - 3
615 615 if zmax == None: zmax = numpy.ceil(numpy.nanmax(avgdB)) + 3
616 616
617 617 self.FTP_WEI = ftp_wei
618 618 self.EXP_CODE = exp_code
619 619 self.SUB_EXP_CODE = sub_exp_code
620 620 self.PLOT_POS = plot_pos
621 621
622 622 self.isConfig = True
623 623
624 624 self.setWinTitle(title)
625 625
626 626
627 627 for i in range(self.nplots):
628 628 pair = dataOut.pairsList[pairsIndexList[i]]
629 629
630 630 chan_index0 = dataOut.channelList.index(pair[0])
631 631 chan_index1 = dataOut.channelList.index(pair[1])
632 632
633 633 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
634 634 title = "Ch%d: %4.2fdB: %s" %(pair[0], noisedB[chan_index0], str_datetime)
635 635 zdB = 10.*numpy.log10(dataOut.data_spc[chan_index0,:,:]/factor)
636 636 axes0 = self.axesList[i*self.__nsubplots]
637 637 axes0.pcolor(x, y, zdB,
638 638 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
639 639 xlabel=xlabel, ylabel=ylabel, title=title,
640 640 ticksize=9, colormap=power_cmap, cblabel='')
641 641
642 642 title = "Ch%d: %4.2fdB: %s" %(pair[1], noisedB[chan_index1], str_datetime)
643 643 zdB = 10.*numpy.log10(dataOut.data_spc[chan_index1,:,:]/factor)
644 644 axes0 = self.axesList[i*self.__nsubplots+1]
645 645 axes0.pcolor(x, y, zdB,
646 646 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
647 647 xlabel=xlabel, ylabel=ylabel, title=title,
648 648 ticksize=9, colormap=power_cmap, cblabel='')
649 649
650 650 coherenceComplex = dataOut.data_cspc[pairsIndexList[i],:,:] / numpy.sqrt( dataOut.data_spc[chan_index0,:,:]*dataOut.data_spc[chan_index1,:,:] )
651 651 coherence = numpy.abs(coherenceComplex)
652 652 # phase = numpy.arctan(-1*coherenceComplex.imag/coherenceComplex.real)*180/numpy.pi
653 653 phase = numpy.arctan2(coherenceComplex.imag, coherenceComplex.real)*180/numpy.pi
654 654
655 655 title = "Coherence Ch%d * Ch%d" %(pair[0], pair[1])
656 656 axes0 = self.axesList[i*self.__nsubplots+2]
657 657 axes0.pcolor(x, y, coherence,
658 658 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=coh_min, zmax=coh_max,
659 659 xlabel=xlabel, ylabel=ylabel, title=title,
660 660 ticksize=9, colormap=coherence_cmap, cblabel='')
661 661
662 662 title = "Phase Ch%d * Ch%d" %(pair[0], pair[1])
663 663 axes0 = self.axesList[i*self.__nsubplots+3]
664 664 axes0.pcolor(x, y, phase,
665 665 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=phase_min, zmax=phase_max,
666 666 xlabel=xlabel, ylabel=ylabel, title=title,
667 667 ticksize=9, colormap=phase_cmap, cblabel='')
668 668
669 669 self.draw()
670 670
671 671 self.save(figpath=figpath,
672 672 figfile=figfile,
673 673 save=save,
674 674 ftp=ftp,
675 675 wr_period=wr_period,
676 676 thisDatetime=thisDatetime)
677 677
678 678 return dataOut
679 679
680 680 @MPDecorator
681 681 class RTIPlot_(Figure):
682 682
683 683 __isConfig = None
684 684 __nsubplots = None
685 685
686 686 WIDTHPROF = None
687 687 HEIGHTPROF = None
688 688 PREFIX = 'rti'
689 689
690 690 def __init__(self):
691 691
692 692 Figure.__init__(self)
693 693 self.timerange = None
694 694 self.isConfig = False
695 695 self.__nsubplots = 1
696 696
697 697 self.WIDTH = 800
698 698 self.HEIGHT = 250
699 699 self.WIDTHPROF = 120
700 700 self.HEIGHTPROF = 0
701 701 self.counter_imagwr = 0
702 702
703 703 self.PLOT_CODE = RTI_CODE
704 704
705 705 self.FTP_WEI = None
706 706 self.EXP_CODE = None
707 707 self.SUB_EXP_CODE = None
708 708 self.PLOT_POS = None
709 709 self.tmin = None
710 710 self.tmax = None
711 711
712 712 self.xmin = None
713 713 self.xmax = None
714 714
715 715 self.figfile = None
716 716
717 717 def getSubplots(self):
718 718
719 719 ncol = 1
720 720 nrow = self.nplots
721 721
722 722 return nrow, ncol
723 723
724 724 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
725 725
726 726 self.__showprofile = showprofile
727 727 self.nplots = nplots
728 728
729 729 ncolspan = 1
730 730 colspan = 1
731 731 if showprofile:
732 732 ncolspan = 7
733 733 colspan = 6
734 734 self.__nsubplots = 2
735 735
736 736 self.createFigure(id = id,
737 737 wintitle = wintitle,
738 738 widthplot = self.WIDTH + self.WIDTHPROF,
739 739 heightplot = self.HEIGHT + self.HEIGHTPROF,
740 740 show=show)
741 741
742 742 nrow, ncol = self.getSubplots()
743 743
744 744 counter = 0
745 745 for y in range(nrow):
746 746 for x in range(ncol):
747 747
748 748 if counter >= self.nplots:
749 749 break
750 750
751 751 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
752 752
753 753 if showprofile:
754 754 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
755 755
756 756 counter += 1
757 757
758 758 def run(self, dataOut, id, wintitle="", channelList=None, showprofile='True',
759 759 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
760 760 timerange=None, colormap='jet',
761 761 save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
762 762 server=None, folder=None, username=None, password=None,
763 763 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, normFactor=None, HEIGHT=None):
764 764
765 765 """
766 766
767 767 Input:
768 768 dataOut :
769 769 id :
770 770 wintitle :
771 771 channelList :
772 772 showProfile :
773 773 xmin : None,
774 774 xmax : None,
775 775 ymin : None,
776 776 ymax : None,
777 777 zmin : None,
778 778 zmax : None
779 779 """
780 780 if dataOut.flagNoData:
781 781 return dataOut
782 782
783 783 #colormap = kwargs.get('colormap', 'jet')
784 784 if HEIGHT is not None:
785 785 self.HEIGHT = HEIGHT
786 786
787 787 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
788 788 return
789 789
790 790 if channelList == None:
791 791 channelIndexList = dataOut.channelIndexList
792 792 else:
793 793 channelIndexList = []
794 794 for channel in channelList:
795 795 if channel not in dataOut.channelList:
796 796 raise ValueError("Channel %d is not in dataOut.channelList")
797 797 channelIndexList.append(dataOut.channelList.index(channel))
798 798
799 799 if normFactor is None:
800 800 factor = dataOut.normFactor
801 801 else:
802 802 factor = normFactor
803 803
804 804 #factor = dataOut.normFactor
805 805 x = dataOut.getTimeRange()
806 806 y = dataOut.getHeiRange()
807 807
808 808 z = dataOut.data_spc/factor
809 809 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
810 810 avg = numpy.average(z, axis=1)
811 811 avgdB = 10.*numpy.log10(avg)
812 812 # avgdB = dataOut.getPower()
813 813
814 814
815 815 thisDatetime = dataOut.datatime
816 816 #thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
817 817 title = wintitle + " RTI" #: %s" %(thisDatetime.strftime("%d-%b-%Y"))
818 818 xlabel = ""
819 819 ylabel = "Range (Km)"
820 820
821 821 update_figfile = False
822 822
823 823 if self.xmax is not None and dataOut.ltctime >= self.xmax: #yong
824 824 self.counter_imagwr = wr_period
825 825 self.isConfig = False
826 826 update_figfile = True
827 827
828 828 if not self.isConfig:
829 829
830 830 nplots = len(channelIndexList)
831 831
832 832 self.setup(id=id,
833 833 nplots=nplots,
834 834 wintitle=wintitle,
835 835 showprofile=showprofile,
836 836 show=show)
837 837
838 838 if timerange != None:
839 839 self.timerange = timerange
840 840
841 841 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
842 842
843 843 noise = dataOut.noise/factor
844 844 noisedB = 10*numpy.log10(noise)
845 845
846 846 if ymin == None: ymin = numpy.nanmin(y)
847 847 if ymax == None: ymax = numpy.nanmax(y)
848 848 if zmin == None: zmin = numpy.floor(numpy.nanmin(noisedB)) - 3
849 849 if zmax == None: zmax = numpy.ceil(numpy.nanmax(avgdB)) + 3
850 850
851 851 self.FTP_WEI = ftp_wei
852 852 self.EXP_CODE = exp_code
853 853 self.SUB_EXP_CODE = sub_exp_code
854 854 self.PLOT_POS = plot_pos
855 855
856 856 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
857 857 self.isConfig = True
858 858 self.figfile = figfile
859 859 update_figfile = True
860 860
861 861 self.setWinTitle(title)
862 862
863 863 for i in range(self.nplots):
864 864 index = channelIndexList[i]
865 865 title = "Channel %d: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
866 866 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
867 867 title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith)
868 868 axes = self.axesList[i*self.__nsubplots]
869 869 zdB = avgdB[index].reshape((1,-1))
870 870 axes.pcolorbuffer(x, y, zdB,
871 871 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
872 872 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
873 873 ticksize=9, cblabel='', cbsize="1%", colormap=colormap)
874 874
875 875 if self.__showprofile:
876 876 axes = self.axesList[i*self.__nsubplots +1]
877 877 axes.pline(avgdB[index], y,
878 878 xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
879 879 xlabel='dB', ylabel='', title='',
880 880 ytick_visible=False,
881 881 grid='x')
882 882
883 883 self.draw()
884 884
885 885 self.save(figpath=figpath,
886 886 figfile=figfile,
887 887 save=save,
888 888 ftp=ftp,
889 889 wr_period=wr_period,
890 890 thisDatetime=thisDatetime,
891 891 update_figfile=update_figfile)
892 892 return dataOut
893 893
894 894 @MPDecorator
895 895 class CoherenceMap_(Figure):
896 896 isConfig = None
897 897 __nsubplots = None
898 898
899 899 WIDTHPROF = None
900 900 HEIGHTPROF = None
901 901 PREFIX = 'cmap'
902 902
903 903 def __init__(self):
904 904 Figure.__init__(self)
905 905 self.timerange = 2*60*60
906 906 self.isConfig = False
907 907 self.__nsubplots = 1
908 908
909 909 self.WIDTH = 800
910 910 self.HEIGHT = 180
911 911 self.WIDTHPROF = 120
912 912 self.HEIGHTPROF = 0
913 913 self.counter_imagwr = 0
914 914
915 915 self.PLOT_CODE = COH_CODE
916 916
917 917 self.FTP_WEI = None
918 918 self.EXP_CODE = None
919 919 self.SUB_EXP_CODE = None
920 920 self.PLOT_POS = None
921 921 self.counter_imagwr = 0
922 922
923 923 self.xmin = None
924 924 self.xmax = None
925 925
926 926 def getSubplots(self):
927 927 ncol = 1
928 928 nrow = self.nplots*2
929 929
930 930 return nrow, ncol
931 931
932 932 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
933 933 self.__showprofile = showprofile
934 934 self.nplots = nplots
935 935
936 936 ncolspan = 1
937 937 colspan = 1
938 938 if showprofile:
939 939 ncolspan = 7
940 940 colspan = 6
941 941 self.__nsubplots = 2
942 942
943 943 self.createFigure(id = id,
944 944 wintitle = wintitle,
945 945 widthplot = self.WIDTH + self.WIDTHPROF,
946 946 heightplot = self.HEIGHT + self.HEIGHTPROF,
947 947 show=True)
948 948
949 949 nrow, ncol = self.getSubplots()
950 950
951 951 for y in range(nrow):
952 952 for x in range(ncol):
953 953
954 954 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
955 955
956 956 if showprofile:
957 957 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
958 958
959 959 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
960 960 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
961 961 timerange=None, phase_min=None, phase_max=None,
962 962 save=False, figpath='./', figfile=None, ftp=False, wr_period=1,
963 963 coherence_cmap='jet', phase_cmap='RdBu_r', show=True,
964 964 server=None, folder=None, username=None, password=None,
965 965 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
966 966
967 967
968 968 if dataOut.flagNoData:
969 969 return dataOut
970 970
971 971 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
972 972 return
973 973
974 974 if pairsList == None:
975 975 pairsIndexList = dataOut.pairsIndexList
976 976 else:
977 977 pairsIndexList = []
978 978 for pair in pairsList:
979 979 if pair not in dataOut.pairsList:
980 980 raise ValueError("Pair %s is not in dataOut.pairsList" %(pair))
981 981 pairsIndexList.append(dataOut.pairsList.index(pair))
982 982
983 983 if pairsIndexList == []:
984 984 return
985 985
986 986 if len(pairsIndexList) > 4:
987 987 pairsIndexList = pairsIndexList[0:4]
988 988
989 989 if phase_min == None:
990 990 phase_min = -180
991 991 if phase_max == None:
992 992 phase_max = 180
993 993
994 994 x = dataOut.getTimeRange()
995 995 y = dataOut.getHeiRange()
996 996
997 997 thisDatetime = dataOut.datatime
998 998
999 999 title = wintitle + " CoherenceMap" #: %s" %(thisDatetime.strftime("%d-%b-%Y"))
1000 1000 xlabel = ""
1001 1001 ylabel = "Range (Km)"
1002 1002 update_figfile = False
1003 1003
1004 1004 if not self.isConfig:
1005 1005 nplots = len(pairsIndexList)
1006 1006 self.setup(id=id,
1007 1007 nplots=nplots,
1008 1008 wintitle=wintitle,
1009 1009 showprofile=showprofile,
1010 1010 show=show)
1011 1011
1012 1012 if timerange != None:
1013 1013 self.timerange = timerange
1014 1014
1015 1015 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
1016 1016
1017 1017 if ymin == None: ymin = numpy.nanmin(y)
1018 1018 if ymax == None: ymax = numpy.nanmax(y)
1019 1019 if zmin == None: zmin = 0.
1020 1020 if zmax == None: zmax = 1.
1021 1021
1022 1022 self.FTP_WEI = ftp_wei
1023 1023 self.EXP_CODE = exp_code
1024 1024 self.SUB_EXP_CODE = sub_exp_code
1025 1025 self.PLOT_POS = plot_pos
1026 1026
1027 1027 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1028 1028
1029 1029 self.isConfig = True
1030 1030 update_figfile = True
1031 1031
1032 1032 self.setWinTitle(title)
1033 1033
1034 1034 for i in range(self.nplots):
1035 1035
1036 1036 pair = dataOut.pairsList[pairsIndexList[i]]
1037 1037
1038 1038 ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i],:,:],axis=0)
1039 1039 powa = numpy.average(dataOut.data_spc[pair[0],:,:],axis=0)
1040 1040 powb = numpy.average(dataOut.data_spc[pair[1],:,:],axis=0)
1041 1041
1042 1042
1043 1043 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
1044 1044 coherence = numpy.abs(avgcoherenceComplex)
1045 1045
1046 1046 z = coherence.reshape((1,-1))
1047 1047
1048 1048 counter = 0
1049 1049
1050 1050 title = "Coherence Ch%d * Ch%d: %s" %(pair[0], pair[1], thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
1051 1051 axes = self.axesList[i*self.__nsubplots*2]
1052 1052 axes.pcolorbuffer(x, y, z,
1053 1053 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
1054 1054 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1055 1055 ticksize=9, cblabel='', colormap=coherence_cmap, cbsize="1%")
1056 1056
1057 1057 if self.__showprofile:
1058 1058 counter += 1
1059 1059 axes = self.axesList[i*self.__nsubplots*2 + counter]
1060 1060 axes.pline(coherence, y,
1061 1061 xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
1062 1062 xlabel='', ylabel='', title='', ticksize=7,
1063 1063 ytick_visible=False, nxticks=5,
1064 1064 grid='x')
1065 1065
1066 1066 counter += 1
1067 1067
1068 1068 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
1069 1069
1070 1070 z = phase.reshape((1,-1))
1071 1071
1072 1072 title = "Phase Ch%d * Ch%d: %s" %(pair[0], pair[1], thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
1073 1073 axes = self.axesList[i*self.__nsubplots*2 + counter]
1074 1074 axes.pcolorbuffer(x, y, z,
1075 1075 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=phase_min, zmax=phase_max,
1076 1076 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1077 1077 ticksize=9, cblabel='', colormap=phase_cmap, cbsize="1%")
1078 1078
1079 1079 if self.__showprofile:
1080 1080 counter += 1
1081 1081 axes = self.axesList[i*self.__nsubplots*2 + counter]
1082 1082 axes.pline(phase, y,
1083 1083 xmin=phase_min, xmax=phase_max, ymin=ymin, ymax=ymax,
1084 1084 xlabel='', ylabel='', title='', ticksize=7,
1085 1085 ytick_visible=False, nxticks=4,
1086 1086 grid='x')
1087 1087
1088 1088 self.draw()
1089 1089
1090 1090 if dataOut.ltctime >= self.xmax:
1091 1091 self.counter_imagwr = wr_period
1092 1092 self.isConfig = False
1093 1093 update_figfile = True
1094 1094
1095 1095 self.save(figpath=figpath,
1096 1096 figfile=figfile,
1097 1097 save=save,
1098 1098 ftp=ftp,
1099 1099 wr_period=wr_period,
1100 1100 thisDatetime=thisDatetime,
1101 1101 update_figfile=update_figfile)
1102 1102
1103 1103 return dataOut
1104 1104
1105 1105 @MPDecorator
1106 1106 class PowerProfilePlot_(Figure):
1107 1107
1108 1108 isConfig = None
1109 1109 __nsubplots = None
1110 1110
1111 1111 WIDTHPROF = None
1112 1112 HEIGHTPROF = None
1113 1113 PREFIX = 'spcprofile'
1114 1114
1115 1115 def __init__(self):
1116 1116 Figure.__init__(self)
1117 1117 self.isConfig = False
1118 1118 self.__nsubplots = 1
1119 1119
1120 1120 self.PLOT_CODE = POWER_CODE
1121 1121
1122 1122 self.WIDTH = 300
1123 1123 self.HEIGHT = 500
1124 1124 self.counter_imagwr = 0
1125 1125
1126 1126 def getSubplots(self):
1127 1127 ncol = 1
1128 1128 nrow = 1
1129 1129
1130 1130 return nrow, ncol
1131 1131
1132 1132 def setup(self, id, nplots, wintitle, show):
1133 1133
1134 1134 self.nplots = nplots
1135 1135
1136 1136 ncolspan = 1
1137 1137 colspan = 1
1138 1138
1139 1139 self.createFigure(id = id,
1140 1140 wintitle = wintitle,
1141 1141 widthplot = self.WIDTH,
1142 1142 heightplot = self.HEIGHT,
1143 1143 show=show)
1144 1144
1145 1145 nrow, ncol = self.getSubplots()
1146 1146
1147 1147 counter = 0
1148 1148 for y in range(nrow):
1149 1149 for x in range(ncol):
1150 1150 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
1151 1151
1152 1152 def run(self, dataOut, id, wintitle="", channelList=None,
1153 1153 xmin=None, xmax=None, ymin=None, ymax=None,
1154 1154 save=False, figpath='./', figfile=None, show=True,
1155 1155 ftp=False, wr_period=1, server=None,
1156 1156 folder=None, username=None, password=None):
1157 1157
1158 1158 if dataOut.flagNoData:
1159 1159 return dataOut
1160 1160
1161 1161
1162 1162 if channelList == None:
1163 1163 channelIndexList = dataOut.channelIndexList
1164 1164 channelList = dataOut.channelList
1165 1165 else:
1166 1166 channelIndexList = []
1167 1167 for channel in channelList:
1168 1168 if channel not in dataOut.channelList:
1169 1169 raise ValueError("Channel %d is not in dataOut.channelList")
1170 1170 channelIndexList.append(dataOut.channelList.index(channel))
1171 1171
1172 1172 factor = dataOut.normFactor
1173 1173
1174 1174 y = dataOut.getHeiRange()
1175 1175
1176 1176 #for voltage
1177 1177 if dataOut.type == 'Voltage':
1178 1178 x = dataOut.data[channelIndexList,:] * numpy.conjugate(dataOut.data[channelIndexList,:])
1179 1179 x = x.real
1180 1180 x = numpy.where(numpy.isfinite(x), x, numpy.NAN)
1181 1181
1182 1182 #for spectra
1183 1183 if dataOut.type == 'Spectra':
1184 1184 x = dataOut.data_spc[channelIndexList,:,:]/factor
1185 1185 x = numpy.where(numpy.isfinite(x), x, numpy.NAN)
1186 1186 x = numpy.average(x, axis=1)
1187 1187
1188 1188
1189 1189 xdB = 10*numpy.log10(x)
1190 1190
1191 1191 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
1192 1192 title = wintitle + " Power Profile %s" %(thisDatetime.strftime("%d-%b-%Y"))
1193 1193 xlabel = "dB"
1194 1194 ylabel = "Range (Km)"
1195 1195
1196 1196 if not self.isConfig:
1197 1197
1198 1198 nplots = 1
1199 1199
1200 1200 self.setup(id=id,
1201 1201 nplots=nplots,
1202 1202 wintitle=wintitle,
1203 1203 show=show)
1204 1204
1205 1205 if ymin == None: ymin = numpy.nanmin(y)
1206 1206 if ymax == None: ymax = numpy.nanmax(y)
1207 1207 if xmin == None: xmin = numpy.nanmin(xdB)*0.9
1208 1208 if xmax == None: xmax = numpy.nanmax(xdB)*1.1
1209 1209
1210 1210 self.isConfig = True
1211 1211
1212 1212 self.setWinTitle(title)
1213 1213
1214 1214 title = "Power Profile: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
1215 1215 axes = self.axesList[0]
1216 1216
1217 1217 legendlabels = ["channel %d"%x for x in channelList]
1218 1218 axes.pmultiline(xdB, y,
1219 1219 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
1220 1220 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels,
1221 1221 ytick_visible=True, nxticks=5,
1222 1222 grid='x')
1223 1223
1224 1224 self.draw()
1225 1225
1226 1226 self.save(figpath=figpath,
1227 1227 figfile=figfile,
1228 1228 save=save,
1229 1229 ftp=ftp,
1230 1230 wr_period=wr_period,
1231 1231 thisDatetime=thisDatetime)
1232 1232
1233 1233 return dataOut
1234 1234
1235 1235 @MPDecorator
1236 1236 class SpectraCutPlot_(Figure):
1237 1237
1238 1238 isConfig = None
1239 1239 __nsubplots = None
1240 1240
1241 1241 WIDTHPROF = None
1242 1242 HEIGHTPROF = None
1243 1243 PREFIX = 'spc_cut'
1244 1244
1245 1245 def __init__(self):
1246 1246 Figure.__init__(self)
1247 1247 self.isConfig = False
1248 1248 self.__nsubplots = 1
1249 1249
1250 1250 self.PLOT_CODE = POWER_CODE
1251 1251
1252 1252 self.WIDTH = 700
1253 1253 self.HEIGHT = 500
1254 1254 self.counter_imagwr = 0
1255 1255
1256 1256 def getSubplots(self):
1257 1257 ncol = 1
1258 1258 nrow = 1
1259 1259
1260 1260 return nrow, ncol
1261 1261
1262 1262 def setup(self, id, nplots, wintitle, show):
1263 1263
1264 1264 self.nplots = nplots
1265 1265
1266 1266 ncolspan = 1
1267 1267 colspan = 1
1268 1268
1269 1269 self.createFigure(id = id,
1270 1270 wintitle = wintitle,
1271 1271 widthplot = self.WIDTH,
1272 1272 heightplot = self.HEIGHT,
1273 1273 show=show)
1274 1274
1275 1275 nrow, ncol = self.getSubplots()
1276 1276
1277 1277 counter = 0
1278 1278 for y in range(nrow):
1279 1279 for x in range(ncol):
1280 1280 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
1281 1281
1282 1282 def run(self, dataOut, id, wintitle="", channelList=None,
1283 1283 xmin=None, xmax=None, ymin=None, ymax=None,
1284 1284 save=False, figpath='./', figfile=None, show=True,
1285 1285 ftp=False, wr_period=1, server=None,
1286 1286 folder=None, username=None, password=None,
1287 1287 xaxis="frequency"):
1288 1288
1289 1289 if dataOut.flagNoData:
1290 1290 return dataOut
1291 1291
1292 1292 if channelList == None:
1293 1293 channelIndexList = dataOut.channelIndexList
1294 1294 channelList = dataOut.channelList
1295 1295 else:
1296 1296 channelIndexList = []
1297 1297 for channel in channelList:
1298 1298 if channel not in dataOut.channelList:
1299 1299 raise ValueError("Channel %d is not in dataOut.channelList")
1300 1300 channelIndexList.append(dataOut.channelList.index(channel))
1301 1301
1302 1302 factor = dataOut.normFactor
1303 1303
1304 1304 y = dataOut.getHeiRange()
1305 1305
1306 1306 z = dataOut.data_spc/factor
1307 1307 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
1308 1308
1309 1309 hei_index = numpy.arange(25)*3 + 20
1310 1310
1311 1311 if xaxis == "frequency":
1312 1312 x = dataOut.getFreqRange()/1000.
1313 1313 zdB = 10*numpy.log10(z[0,:,hei_index])
1314 1314 xlabel = "Frequency (kHz)"
1315 1315 ylabel = "Power (dB)"
1316 1316
1317 1317 elif xaxis == "time":
1318 1318 x = dataOut.getAcfRange()
1319 1319 zdB = z[0,:,hei_index]
1320 1320 xlabel = "Time (ms)"
1321 1321 ylabel = "ACF"
1322 1322
1323 1323 else:
1324 1324 x = dataOut.getVelRange()
1325 1325 zdB = 10*numpy.log10(z[0,:,hei_index])
1326 1326 xlabel = "Velocity (m/s)"
1327 1327 ylabel = "Power (dB)"
1328 1328
1329 1329 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
1330 1330 title = wintitle + " Range Cuts %s" %(thisDatetime.strftime("%d-%b-%Y"))
1331 1331
1332 1332 if not self.isConfig:
1333 1333
1334 1334 nplots = 1
1335 1335
1336 1336 self.setup(id=id,
1337 1337 nplots=nplots,
1338 1338 wintitle=wintitle,
1339 1339 show=show)
1340 1340
1341 1341 if xmin == None: xmin = numpy.nanmin(x)*0.9
1342 1342 if xmax == None: xmax = numpy.nanmax(x)*1.1
1343 1343 if ymin == None: ymin = numpy.nanmin(zdB)
1344 1344 if ymax == None: ymax = numpy.nanmax(zdB)
1345 1345
1346 1346 self.isConfig = True
1347 1347
1348 1348 self.setWinTitle(title)
1349 1349
1350 1350 title = "Spectra Cuts: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
1351 1351 axes = self.axesList[0]
1352 1352
1353 1353 legendlabels = ["Range = %dKm" %y[i] for i in hei_index]
1354 1354
1355 1355 axes.pmultilineyaxis( x, zdB,
1356 1356 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
1357 1357 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels,
1358 1358 ytick_visible=True, nxticks=5,
1359 1359 grid='x')
1360 1360
1361 1361 self.draw()
1362 1362
1363 1363 self.save(figpath=figpath,
1364 1364 figfile=figfile,
1365 1365 save=save,
1366 1366 ftp=ftp,
1367 1367 wr_period=wr_period,
1368 1368 thisDatetime=thisDatetime)
1369 1369
1370 1370 return dataOut
1371 1371
1372 1372 @MPDecorator
1373 1373 class Noise_(Figure):
1374 1374
1375 1375 isConfig = None
1376 1376 __nsubplots = None
1377 1377
1378 1378 PREFIX = 'noise'
1379 1379
1380 1380
1381 1381 def __init__(self):
1382 1382 Figure.__init__(self)
1383 1383 self.timerange = 24*60*60
1384 1384 self.isConfig = False
1385 1385 self.__nsubplots = 1
1386 1386 self.counter_imagwr = 0
1387 1387 self.WIDTH = 800
1388 1388 self.HEIGHT = 400
1389 1389 self.WIDTHPROF = 120
1390 1390 self.HEIGHTPROF = 0
1391 1391 self.xdata = None
1392 1392 self.ydata = None
1393 1393
1394 1394 self.PLOT_CODE = NOISE_CODE
1395 1395
1396 1396 self.FTP_WEI = None
1397 1397 self.EXP_CODE = None
1398 1398 self.SUB_EXP_CODE = None
1399 1399 self.PLOT_POS = None
1400 1400 self.figfile = None
1401 1401
1402 1402 self.xmin = None
1403 1403 self.xmax = None
1404 1404
1405 1405 def getSubplots(self):
1406 1406
1407 1407 ncol = 1
1408 1408 nrow = 1
1409 1409
1410 1410 return nrow, ncol
1411 1411
1412 1412 def openfile(self, filename):
1413 1413 dirname = os.path.dirname(filename)
1414 1414
1415 1415 if not os.path.exists(dirname):
1416 1416 os.mkdir(dirname)
1417 1417
1418 1418 f = open(filename,'w+')
1419 1419 f.write('\n\n')
1420 1420 f.write('JICAMARCA RADIO OBSERVATORY - Noise \n')
1421 1421 f.write('DD MM YYYY HH MM SS Channel0 Channel1 Channel2 Channel3\n\n' )
1422 1422 f.close()
1423 1423
1424 1424 def save_data(self, filename_phase, data, data_datetime):
1425 1425
1426 1426 f=open(filename_phase,'a')
1427 1427
1428 1428 timetuple_data = data_datetime.timetuple()
1429 1429 day = str(timetuple_data.tm_mday)
1430 1430 month = str(timetuple_data.tm_mon)
1431 1431 year = str(timetuple_data.tm_year)
1432 1432 hour = str(timetuple_data.tm_hour)
1433 1433 minute = str(timetuple_data.tm_min)
1434 1434 second = str(timetuple_data.tm_sec)
1435 1435
1436 1436 data_msg = ''
1437 1437 for i in range(len(data)):
1438 1438 data_msg += str(data[i]) + ' '
1439 1439
1440 1440 f.write(day+' '+month+' '+year+' '+hour+' '+minute+' '+second+' ' + data_msg + '\n')
1441 1441 f.close()
1442 1442
1443 1443
1444 1444 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
1445 1445
1446 1446 self.__showprofile = showprofile
1447 1447 self.nplots = nplots
1448 1448
1449 1449 ncolspan = 7
1450 1450 colspan = 6
1451 1451 self.__nsubplots = 2
1452 1452
1453 1453 self.createFigure(id = id,
1454 1454 wintitle = wintitle,
1455 1455 widthplot = self.WIDTH+self.WIDTHPROF,
1456 1456 heightplot = self.HEIGHT+self.HEIGHTPROF,
1457 1457 show=show)
1458 1458
1459 1459 nrow, ncol = self.getSubplots()
1460 1460
1461 1461 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
1462 1462
1463 1463
1464 1464 def run(self, dataOut, id, wintitle="", channelList=None, showprofile='True',
1465 1465 xmin=None, xmax=None, ymin=None, ymax=None,
1466 1466 timerange=None,
1467 1467 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
1468 1468 server=None, folder=None, username=None, password=None,
1469 1469 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1470 1470
1471 1471 if dataOut.flagNoData:
1472 1472 return dataOut
1473 1473
1474 1474 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
1475 1475 return
1476 1476
1477 1477 if channelList == None:
1478 1478 channelIndexList = dataOut.channelIndexList
1479 1479 channelList = dataOut.channelList
1480 1480 else:
1481 1481 channelIndexList = []
1482 1482 for channel in channelList:
1483 1483 if channel not in dataOut.channelList:
1484 1484 raise ValueError("Channel %d is not in dataOut.channelList")
1485 1485 channelIndexList.append(dataOut.channelList.index(channel))
1486 1486
1487 1487 x = dataOut.getTimeRange()
1488 1488 #y = dataOut.getHeiRange()
1489 1489 factor = dataOut.normFactor
1490 1490 noise = dataOut.noise[channelIndexList]/factor
1491 1491 noisedB = 10*numpy.log10(noise)
1492 1492
1493 1493 thisDatetime = dataOut.datatime
1494 1494
1495 1495 title = wintitle + " Noise" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
1496 1496 xlabel = ""
1497 1497 ylabel = "Intensity (dB)"
1498 1498 update_figfile = False
1499 1499
1500 1500 if not self.isConfig:
1501 1501
1502 1502 nplots = 1
1503 1503
1504 1504 self.setup(id=id,
1505 1505 nplots=nplots,
1506 1506 wintitle=wintitle,
1507 1507 showprofile=showprofile,
1508 1508 show=show)
1509 1509
1510 1510 if timerange != None:
1511 1511 self.timerange = timerange
1512 1512
1513 1513 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
1514 1514
1515 1515 if ymin == None: ymin = numpy.floor(numpy.nanmin(noisedB)) - 10.0
1516 1516 if ymax == None: ymax = numpy.nanmax(noisedB) + 10.0
1517 1517
1518 1518 self.FTP_WEI = ftp_wei
1519 1519 self.EXP_CODE = exp_code
1520 1520 self.SUB_EXP_CODE = sub_exp_code
1521 1521 self.PLOT_POS = plot_pos
1522 1522
1523 1523
1524 1524 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1525 1525 self.isConfig = True
1526 1526 self.figfile = figfile
1527 1527 self.xdata = numpy.array([])
1528 1528 self.ydata = numpy.array([])
1529 1529
1530 1530 update_figfile = True
1531 1531
1532 1532 #open file beacon phase
1533 1533 path = '%s%03d' %(self.PREFIX, self.id)
1534 1534 noise_file = os.path.join(path,'%s.txt'%self.name)
1535 1535 self.filename_noise = os.path.join(figpath,noise_file)
1536 1536
1537 1537 self.setWinTitle(title)
1538 1538
1539 1539 title = "Noise %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1540 1540
1541 1541 legendlabels = ["channel %d"%(idchannel) for idchannel in channelList]
1542 1542 axes = self.axesList[0]
1543 1543
1544 1544 self.xdata = numpy.hstack((self.xdata, x[0:1]))
1545 1545
1546 1546 if len(self.ydata)==0:
1547 1547 self.ydata = noisedB.reshape(-1,1)
1548 1548 else:
1549 1549 self.ydata = numpy.hstack((self.ydata, noisedB.reshape(-1,1)))
1550 1550
1551 1551
1552 1552 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
1553 1553 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
1554 1554 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
1555 1555 XAxisAsTime=True, grid='both'
1556 1556 )
1557 1557
1558 1558 self.draw()
1559 1559
1560 1560 if dataOut.ltctime >= self.xmax:
1561 1561 self.counter_imagwr = wr_period
1562 1562 self.isConfig = False
1563 1563 update_figfile = True
1564 1564
1565 1565 self.save(figpath=figpath,
1566 1566 figfile=figfile,
1567 1567 save=save,
1568 1568 ftp=ftp,
1569 1569 wr_period=wr_period,
1570 1570 thisDatetime=thisDatetime,
1571 1571 update_figfile=update_figfile)
1572 1572
1573 1573 #store data beacon phase
1574 1574 if save:
1575 1575 self.save_data(self.filename_noise, noisedB, thisDatetime)
1576 1576
1577 1577 return dataOut
1578 1578
1579 1579 @MPDecorator
1580 1580 class BeaconPhase_(Figure):
1581 1581
1582 1582 __isConfig = None
1583 1583 __nsubplots = None
1584 1584
1585 1585 PREFIX = 'beacon_phase'
1586 1586
1587 1587 def __init__(self):
1588 1588 Figure.__init__(self)
1589 1589 self.timerange = 24*60*60
1590 1590 self.isConfig = False
1591 1591 self.__nsubplots = 1
1592 1592 self.counter_imagwr = 0
1593 1593 self.WIDTH = 800
1594 1594 self.HEIGHT = 400
1595 1595 self.WIDTHPROF = 120
1596 1596 self.HEIGHTPROF = 0
1597 1597 self.xdata = None
1598 1598 self.ydata = None
1599 1599
1600 1600 self.PLOT_CODE = BEACON_CODE
1601 1601
1602 1602 self.FTP_WEI = None
1603 1603 self.EXP_CODE = None
1604 1604 self.SUB_EXP_CODE = None
1605 1605 self.PLOT_POS = None
1606 1606
1607 1607 self.filename_phase = None
1608 1608
1609 1609 self.figfile = None
1610 1610
1611 1611 self.xmin = None
1612 1612 self.xmax = None
1613 1613
1614 1614 def getSubplots(self):
1615 1615
1616 1616 ncol = 1
1617 1617 nrow = 1
1618 1618
1619 1619 return nrow, ncol
1620 1620
1621 1621 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
1622 1622
1623 1623 self.__showprofile = showprofile
1624 1624 self.nplots = nplots
1625 1625
1626 1626 ncolspan = 7
1627 1627 colspan = 6
1628 1628 self.__nsubplots = 2
1629 1629
1630 1630 self.createFigure(id = id,
1631 1631 wintitle = wintitle,
1632 1632 widthplot = self.WIDTH+self.WIDTHPROF,
1633 1633 heightplot = self.HEIGHT+self.HEIGHTPROF,
1634 1634 show=show)
1635 1635
1636 1636 nrow, ncol = self.getSubplots()
1637 1637
1638 1638 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
1639 1639
1640 1640 def save_phase(self, filename_phase):
1641 1641 f = open(filename_phase,'w+')
1642 1642 f.write('\n\n')
1643 1643 f.write('JICAMARCA RADIO OBSERVATORY - Beacon Phase \n')
1644 1644 f.write('DD MM YYYY HH MM SS pair(2,0) pair(2,1) pair(2,3) pair(2,4)\n\n' )
1645 1645 f.close()
1646 1646
1647 1647 def save_data(self, filename_phase, data, data_datetime):
1648 1648 f=open(filename_phase,'a')
1649 1649 timetuple_data = data_datetime.timetuple()
1650 1650 day = str(timetuple_data.tm_mday)
1651 1651 month = str(timetuple_data.tm_mon)
1652 1652 year = str(timetuple_data.tm_year)
1653 1653 hour = str(timetuple_data.tm_hour)
1654 1654 minute = str(timetuple_data.tm_min)
1655 1655 second = str(timetuple_data.tm_sec)
1656 1656 f.write(day+' '+month+' '+year+' '+hour+' '+minute+' '+second+' '+str(data[0])+' '+str(data[1])+' '+str(data[2])+' '+str(data[3])+'\n')
1657 1657 f.close()
1658 1658
1659 1659
1660 1660 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
1661 1661 xmin=None, xmax=None, ymin=None, ymax=None, hmin=None, hmax=None,
1662 1662 timerange=None,
1663 1663 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
1664 1664 server=None, folder=None, username=None, password=None,
1665 1665 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1666 1666
1667 1667 if dataOut.flagNoData:
1668 1668 return dataOut
1669 1669
1670 1670 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
1671 1671 return
1672 1672
1673 1673 if pairsList == None:
1674 1674 pairsIndexList = dataOut.pairsIndexList[:10]
1675 1675 else:
1676 1676 pairsIndexList = []
1677 1677 for pair in pairsList:
1678 1678 if pair not in dataOut.pairsList:
1679 1679 raise ValueError("Pair %s is not in dataOut.pairsList" %(pair))
1680 1680 pairsIndexList.append(dataOut.pairsList.index(pair))
1681 1681
1682 1682 if pairsIndexList == []:
1683 1683 return
1684 1684
1685 1685 # if len(pairsIndexList) > 4:
1686 1686 # pairsIndexList = pairsIndexList[0:4]
1687 1687
1688 1688 hmin_index = None
1689 1689 hmax_index = None
1690 1690
1691 1691 if hmin != None and hmax != None:
1692 1692 indexes = numpy.arange(dataOut.nHeights)
1693 1693 hmin_list = indexes[dataOut.heightList >= hmin]
1694 1694 hmax_list = indexes[dataOut.heightList <= hmax]
1695 1695
1696 1696 if hmin_list.any():
1697 1697 hmin_index = hmin_list[0]
1698 1698
1699 1699 if hmax_list.any():
1700 1700 hmax_index = hmax_list[-1]+1
1701 1701
1702 1702 x = dataOut.getTimeRange()
1703 1703 #y = dataOut.getHeiRange()
1704 1704
1705 1705
1706 1706 thisDatetime = dataOut.datatime
1707 1707
1708 1708 title = wintitle + " Signal Phase" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
1709 1709 xlabel = "Local Time"
1710 1710 ylabel = "Phase (degrees)"
1711 1711
1712 1712 update_figfile = False
1713 1713
1714 1714 nplots = len(pairsIndexList)
1715 1715 #phase = numpy.zeros((len(pairsIndexList),len(dataOut.beacon_heiIndexList)))
1716 1716 phase_beacon = numpy.zeros(len(pairsIndexList))
1717 1717 for i in range(nplots):
1718 1718 pair = dataOut.pairsList[pairsIndexList[i]]
1719 1719 ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i], :, hmin_index:hmax_index], axis=0)
1720 1720 powa = numpy.average(dataOut.data_spc[pair[0], :, hmin_index:hmax_index], axis=0)
1721 1721 powb = numpy.average(dataOut.data_spc[pair[1], :, hmin_index:hmax_index], axis=0)
1722 1722 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
1723 1723 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
1724 1724
1725 1725 if dataOut.beacon_heiIndexList:
1726 1726 phase_beacon[i] = numpy.average(phase[dataOut.beacon_heiIndexList])
1727 1727 else:
1728 1728 phase_beacon[i] = numpy.average(phase)
1729 1729
1730 1730 if not self.isConfig:
1731 1731
1732 1732 nplots = len(pairsIndexList)
1733 1733
1734 1734 self.setup(id=id,
1735 1735 nplots=nplots,
1736 1736 wintitle=wintitle,
1737 1737 showprofile=showprofile,
1738 1738 show=show)
1739 1739
1740 1740 if timerange != None:
1741 1741 self.timerange = timerange
1742 1742
1743 1743 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
1744 1744
1745 1745 if ymin == None: ymin = 0
1746 1746 if ymax == None: ymax = 360
1747 1747
1748 1748 self.FTP_WEI = ftp_wei
1749 1749 self.EXP_CODE = exp_code
1750 1750 self.SUB_EXP_CODE = sub_exp_code
1751 1751 self.PLOT_POS = plot_pos
1752 1752
1753 1753 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1754 1754 self.isConfig = True
1755 1755 self.figfile = figfile
1756 1756 self.xdata = numpy.array([])
1757 1757 self.ydata = numpy.array([])
1758 1758
1759 1759 update_figfile = True
1760 1760
1761 1761 #open file beacon phase
1762 1762 path = '%s%03d' %(self.PREFIX, self.id)
1763 1763 beacon_file = os.path.join(path,'%s.txt'%self.name)
1764 1764 self.filename_phase = os.path.join(figpath,beacon_file)
1765 1765 #self.save_phase(self.filename_phase)
1766 1766
1767 1767
1768 1768 #store data beacon phase
1769 1769 #self.save_data(self.filename_phase, phase_beacon, thisDatetime)
1770 1770
1771 1771 self.setWinTitle(title)
1772 1772
1773 1773
1774 1774 title = "Phase Plot %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1775 1775
1776 1776 legendlabels = ["Pair (%d,%d)"%(pair[0], pair[1]) for pair in dataOut.pairsList]
1777 1777
1778 1778 axes = self.axesList[0]
1779 1779
1780 1780 self.xdata = numpy.hstack((self.xdata, x[0:1]))
1781 1781
1782 1782 if len(self.ydata)==0:
1783 1783 self.ydata = phase_beacon.reshape(-1,1)
1784 1784 else:
1785 1785 self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1)))
1786 1786
1787 1787
1788 1788 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
1789 1789 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
1790 1790 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
1791 1791 XAxisAsTime=True, grid='both'
1792 1792 )
1793 1793
1794 1794 self.draw()
1795 1795
1796 1796 if dataOut.ltctime >= self.xmax:
1797 1797 self.counter_imagwr = wr_period
1798 1798 self.isConfig = False
1799 1799 update_figfile = True
1800 1800
1801 1801 self.save(figpath=figpath,
1802 1802 figfile=figfile,
1803 1803 save=save,
1804 1804 ftp=ftp,
1805 1805 wr_period=wr_period,
1806 1806 thisDatetime=thisDatetime,
1807 1807 update_figfile=update_figfile)
1808 1808
1809 1809 return dataOut
Show file before | Show file after | Hide comments
@@ -1,1260 +1,1266
1 1 import itertools
2 2
3 3 import numpy
4 4
5 5 from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator, Operation
6 6 from schainpy.model.data.jrodata import Spectra
7 7 from schainpy.model.data.jrodata import hildebrand_sekhon
8 8 from schainpy.utils import log
9 9
10 10 @MPDecorator
11 11 class SpectraProc(ProcessingUnit):
12 12
13 13
14 14 def __init__(self):
15 15
16 16 ProcessingUnit.__init__(self)
17 17
18 18 self.buffer = None
19 19 self.firstdatatime = None
20 20 self.profIndex = 0
21 21 self.dataOut = Spectra()
22 22 self.id_min = None
23 23 self.id_max = None
24 24 self.setupReq = False #Agregar a todas las unidades de proc
25 25
26 26 def __updateSpecFromVoltage(self):
27 27
28 28 self.dataOut.timeZone = self.dataIn.timeZone
29 29 self.dataOut.dstFlag = self.dataIn.dstFlag
30 30 self.dataOut.errorCount = self.dataIn.errorCount
31 31 self.dataOut.useLocalTime = self.dataIn.useLocalTime
32 32 try:
33 33 self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
34 34 except:
35 35 pass
36 36 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
37 37 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
38 38 self.dataOut.channelList = self.dataIn.channelList
39 39 self.dataOut.heightList = self.dataIn.heightList
40 40 self.dataOut.dtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
41 41
42 42 self.dataOut.nBaud = self.dataIn.nBaud
43 43 self.dataOut.nCode = self.dataIn.nCode
44 44 self.dataOut.code = self.dataIn.code
45 45 self.dataOut.nProfiles = self.dataOut.nFFTPoints
46 46
47 47 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
48 48 self.dataOut.utctime = self.firstdatatime
49 49 # asumo q la data esta decodificada
50 50 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData
51 51 # asumo q la data esta sin flip
52 52 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData
53 53 self.dataOut.flagShiftFFT = False
54 54
55 55 self.dataOut.nCohInt = self.dataIn.nCohInt
56 56 self.dataOut.nIncohInt = 1
57 57
58 58 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
59 59
60 60 self.dataOut.frequency = self.dataIn.frequency
61 61 self.dataOut.realtime = self.dataIn.realtime
62 62
63 63 self.dataOut.azimuth = self.dataIn.azimuth
64 64 self.dataOut.zenith = self.dataIn.zenith
65 65
66 66 self.dataOut.beam.codeList = self.dataIn.beam.codeList
67 67 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
68 68 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
69 69
70 70 def __getFft(self):
71 71 """
72 72 Convierte valores de Voltaje a Spectra
73 73
74 74 Affected:
75 75 self.dataOut.data_spc
76 76 self.dataOut.data_cspc
77 77 self.dataOut.data_dc
78 78 self.dataOut.heightList
79 79 self.profIndex
80 80 self.buffer
81 81 self.dataOut.flagNoData
82 82 """
83 83 fft_volt = numpy.fft.fft(
84 84 self.buffer, n=self.dataOut.nFFTPoints, axis=1)
85 85 fft_volt = fft_volt.astype(numpy.dtype('complex'))
86 86 dc = fft_volt[:, 0, :]
87 87
88 88 # calculo de self-spectra
89 89 fft_volt = numpy.fft.fftshift(fft_volt, axes=(1,))
90 90 spc = fft_volt * numpy.conjugate(fft_volt)
91 #print("spcch0",spc[0])
91 92 spc = spc.real
92 93
93 94 blocksize = 0
94 95 blocksize += dc.size
95 96 blocksize += spc.size
96 97
97 98 #print("spc :",spc.shape)
98 99 data_wr = None
99 100 if self.dataOut.flagWR:
100 101 data_wr = fft_volt
101 102 blocksize = fft_volt.size
102 103
103 104 cspc = None
104 105 pairIndex = 0
105 106 if self.dataOut.pairsList != None:
106 107 # calculo de cross-spectra
107 108 cspc = numpy.zeros(
108 109 (self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
109 110 for pair in self.dataOut.pairsList:
110 111 if pair[0] not in self.dataOut.channelList:
111 112 raise ValueError("Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" % (
112 113 str(pair), str(self.dataOut.channelList)))
113 114 if pair[1] not in self.dataOut.channelList:
114 115 raise ValueError("Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" % (
115 116 str(pair), str(self.dataOut.channelList)))
116 117
117 118 cspc[pairIndex, :, :] = fft_volt[pair[0], :, :] * \
118 119 numpy.conjugate(fft_volt[pair[1], :, :])
119 120 pairIndex += 1
120 121 blocksize += cspc.size
121 122
122 123 self.dataOut.data_spc = spc
123 124 self.dataOut.data_cspc = cspc
124 125 self.dataOut.data_wr = data_wr
125 126 self.dataOut.data_dc = dc
126 127 self.dataOut.blockSize = blocksize
127 128 self.dataOut.flagShiftFFT = False
128 129
129 130 def run(self, nProfiles=None, nFFTPoints=None, pairsList=[], ippFactor=None, shift_fft=False,flagWR= 0):
130 131
131 132 self.dataOut.flagWR = flagWR
132 133
133 134 if self.dataIn.type == "Spectra":
134 135 self.dataOut.copy(self.dataIn)
135 136
136 137 if shift_fft:
137 138 #desplaza a la derecha en el eje 2 determinadas posiciones
138 139 shift = int(self.dataOut.nFFTPoints/2)
139 140 self.dataOut.data_spc = numpy.roll(self.dataOut.data_spc, shift , axis=1)
140 141
141 142 if self.dataOut.data_cspc is not None:
142 143 #desplaza a la derecha en el eje 2 determinadas posiciones
143 144 self.dataOut.data_cspc = numpy.roll(self.dataOut.data_cspc, shift, axis=1)
144 145
145 146 return True
146 147
147 148 if self.dataIn.type == "Voltage":
148 149 #print("VOLTAGE INPUT SPECTRA")
149 150 self.dataOut.flagNoData = True
150 151
151 152 if nFFTPoints == None:
152 153 raise ValueError("This SpectraProc.run() need nFFTPoints input variable")
153 154
154 155 if nProfiles == None:
155 156 nProfiles = nFFTPoints
156 157
157 158 if ippFactor == None:
158 159 ippFactor = 1
159 160
160 161 self.dataOut.ippFactor = ippFactor
161 162
162 163 self.dataOut.nFFTPoints = nFFTPoints
163 164 self.dataOut.pairsList = pairsList
164 165
165 166 if self.buffer is None:
166 167 self.buffer = numpy.zeros((self.dataIn.nChannels,
167 168 nProfiles,
168 169 self.dataIn.nHeights),
169 170 dtype='complex')
170 171 #print("buffer :",self.buffer.shape)
171 172
172 173 if self.dataIn.flagDataAsBlock:
173 174 nVoltProfiles = self.dataIn.data.shape[1]
174 175
175 176 if nVoltProfiles == nProfiles:
176 177 self.buffer = self.dataIn.data.copy()
177 178 self.profIndex = nVoltProfiles
178 179
179 180 elif nVoltProfiles < nProfiles:
180 181
181 182 if self.profIndex == 0:
182 183 self.id_min = 0
183 184 self.id_max = nVoltProfiles
184 185
185 186 self.buffer[:, self.id_min:self.id_max,
186 187 :] = self.dataIn.data
187 188 self.profIndex += nVoltProfiles
188 189 self.id_min += nVoltProfiles
189 190 self.id_max += nVoltProfiles
190 191 else:
191 192 raise ValueError("The type object %s has %d profiles, it should just has %d profiles" % (
192 193 self.dataIn.type, self.dataIn.data.shape[1], nProfiles))
193 194 self.dataOut.flagNoData = True
194 195 return 0
195 196 else:
196 197 #print("Spectra ",self.profIndex)
197 198 self.buffer[:, self.profIndex, :] = self.dataIn.data.copy()
198 199 self.profIndex += 1
199 200
200 201 if self.firstdatatime == None:
201 202 self.firstdatatime = self.dataIn.utctime
202 203
203 204 if self.profIndex == nProfiles:
204 205 self.__updateSpecFromVoltage()
205 206 self.__getFft()
206 207 #print(" DATAOUT SHAPE SPEC",self.dataOut.data_spc.shape)
207 208
208 209 self.dataOut.flagNoData = False
209 210 self.firstdatatime = None
210 211 self.profIndex = 0
211 212
212 213 return True
213 214
214 215 raise ValueError("The type of input object '%s' is not valid" % (
215 216 self.dataIn.type))
216 217
217 218 def __selectPairs(self, pairsList):
218 219
219 220 if not pairsList:
220 221 return
221 222
222 223 pairs = []
223 224 pairsIndex = []
224 225
225 226 for pair in pairsList:
226 227 if pair[0] not in self.dataOut.channelList or pair[1] not in self.dataOut.channelList:
227 228 continue
228 229 pairs.append(pair)
229 230 pairsIndex.append(pairs.index(pair))
230 231
231 232 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndex]
232 233 self.dataOut.pairsList = pairs
233 234
234 235 return
235 236
236 237 def __selectPairsByChannel(self, channelList=None):
237 238
238 239 if channelList == None:
239 240 return
240 241
241 242 pairsIndexListSelected = []
242 243 for pairIndex in self.dataOut.pairsIndexList:
243 244 # First pair
244 245 if self.dataOut.pairsList[pairIndex][0] not in channelList:
245 246 continue
246 247 # Second pair
247 248 if self.dataOut.pairsList[pairIndex][1] not in channelList:
248 249 continue
249 250
250 251 pairsIndexListSelected.append(pairIndex)
251 252
252 253 if not pairsIndexListSelected:
253 254 self.dataOut.data_cspc = None
254 255 self.dataOut.pairsList = []
255 256 return
256 257
257 258 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
258 259 self.dataOut.pairsList = [self.dataOut.pairsList[i]
259 260 for i in pairsIndexListSelected]
260 261
261 262 return
262 263
263 264 def selectChannels(self, channelList):
264 265
265 266 channelIndexList = []
266 267
267 268 for channel in channelList:
268 269 if channel not in self.dataOut.channelList:
269 270 raise ValueError("Error selecting channels, Channel %d is not valid.\nAvailable channels = %s" % (
270 271 channel, str(self.dataOut.channelList)))
271 272
272 273 index = self.dataOut.channelList.index(channel)
273 274 channelIndexList.append(index)
274 275
275 276 self.selectChannelsByIndex(channelIndexList)
276 277
277 278 def selectChannelsByIndex(self, channelIndexList):
278 279 """
279 280 Selecciona un bloque de datos en base a canales segun el channelIndexList
280 281
281 282 Input:
282 283 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
283 284
284 285 Affected:
285 286 self.dataOut.data_spc
286 287 self.dataOut.channelIndexList
287 288 self.dataOut.nChannels
288 289
289 290 Return:
290 291 None
291 292 """
292 293
293 294 for channelIndex in channelIndexList:
294 295 if channelIndex not in self.dataOut.channelIndexList:
295 296 raise ValueError("Error selecting channels: The value %d in channelIndexList is not valid.\nAvailable channel indexes = " % (
296 297 channelIndex, self.dataOut.channelIndexList))
297 298
298 299 data_spc = self.dataOut.data_spc[channelIndexList, :]
299 300 data_dc = self.dataOut.data_dc[channelIndexList, :]
300 301
301 302 self.dataOut.data_spc = data_spc
302 303 self.dataOut.data_dc = data_dc
303 304
304 305 # self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
305 306 self.dataOut.channelList = range(len(channelIndexList))
306 307 self.__selectPairsByChannel(channelIndexList)
307 308
308 309 return 1
309 310
310 311
311 312 def selectFFTs(self, minFFT, maxFFT ):
312 313 """
313 314 Selecciona un bloque de datos en base a un grupo de valores de puntos FFTs segun el rango
314 315 minFFT<= FFT <= maxFFT
315 316 """
316 317
317 318 if (minFFT > maxFFT):
318 319 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minFFT, maxFFT))
319 320
320 321 if (minFFT < self.dataOut.getFreqRange()[0]):
321 322 minFFT = self.dataOut.getFreqRange()[0]
322 323
323 324 if (maxFFT > self.dataOut.getFreqRange()[-1]):
324 325 maxFFT = self.dataOut.getFreqRange()[-1]
325 326
326 327 minIndex = 0
327 328 maxIndex = 0
328 329 FFTs = self.dataOut.getFreqRange()
329 330
330 331 inda = numpy.where(FFTs >= minFFT)
331 332 indb = numpy.where(FFTs <= maxFFT)
332 333
333 334 try:
334 335 minIndex = inda[0][0]
335 336 except:
336 337 minIndex = 0
337 338
338 339 try:
339 340 maxIndex = indb[0][-1]
340 341 except:
341 342 maxIndex = len(FFTs)
342 343
343 344 self.selectFFTsByIndex(minIndex, maxIndex)
344 345
345 346 return 1
346 347
347 348
348 349 def setH0(self, h0, deltaHeight = None):
349 350
350 351 if not deltaHeight:
351 352 deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0]
352 353
353 354 nHeights = self.dataOut.nHeights
354 355
355 356 newHeiRange = h0 + numpy.arange(nHeights)*deltaHeight
356 357
357 358 self.dataOut.heightList = newHeiRange
358 359
359 360
360 361 def selectHeights(self, minHei, maxHei):
361 362 """
362 363 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
363 364 minHei <= height <= maxHei
364 365
365 366 Input:
366 367 minHei : valor minimo de altura a considerar
367 368 maxHei : valor maximo de altura a considerar
368 369
369 370 Affected:
370 371 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
371 372
372 373 Return:
373 374 1 si el metodo se ejecuto con exito caso contrario devuelve 0
374 375 """
375 376
376 377
377 378 if (minHei > maxHei):
378 379 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minHei, maxHei))
379 380
380 381 if (minHei < self.dataOut.heightList[0]):
381 382 minHei = self.dataOut.heightList[0]
382 383
383 384 if (maxHei > self.dataOut.heightList[-1]):
384 385 maxHei = self.dataOut.heightList[-1]
385 386
386 387 minIndex = 0
387 388 maxIndex = 0
388 389 heights = self.dataOut.heightList
389 390
390 391 inda = numpy.where(heights >= minHei)
391 392 indb = numpy.where(heights <= maxHei)
392 393
393 394 try:
394 395 minIndex = inda[0][0]
395 396 except:
396 397 minIndex = 0
397 398
398 399 try:
399 400 maxIndex = indb[0][-1]
400 401 except:
401 402 maxIndex = len(heights)
402 403
403 404 self.selectHeightsByIndex(minIndex, maxIndex)
404 405
405 406
406 407 return 1
407 408
408 409 def getBeaconSignal(self, tauindex=0, channelindex=0, hei_ref=None):
409 410 newheis = numpy.where(
410 411 self.dataOut.heightList > self.dataOut.radarControllerHeaderObj.Taus[tauindex])
411 412
412 413 if hei_ref != None:
413 414 newheis = numpy.where(self.dataOut.heightList > hei_ref)
414 415
415 416 minIndex = min(newheis[0])
416 417 maxIndex = max(newheis[0])
417 418 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
418 419 heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
419 420
420 421 # determina indices
421 422 nheis = int(self.dataOut.radarControllerHeaderObj.txB /
422 423 (self.dataOut.heightList[1] - self.dataOut.heightList[0]))
423 424 avg_dB = 10 * \
424 425 numpy.log10(numpy.sum(data_spc[channelindex, :, :], axis=0))
425 426 beacon_dB = numpy.sort(avg_dB)[-nheis:]
426 427 beacon_heiIndexList = []
427 428 for val in avg_dB.tolist():
428 429 if val >= beacon_dB[0]:
429 430 beacon_heiIndexList.append(avg_dB.tolist().index(val))
430 431
431 432 #data_spc = data_spc[:,:,beacon_heiIndexList]
432 433 data_cspc = None
433 434 if self.dataOut.data_cspc is not None:
434 435 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
435 436 #data_cspc = data_cspc[:,:,beacon_heiIndexList]
436 437
437 438 data_dc = None
438 439 if self.dataOut.data_dc is not None:
439 440 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
440 441 #data_dc = data_dc[:,beacon_heiIndexList]
441 442
442 443 self.dataOut.data_spc = data_spc
443 444 self.dataOut.data_cspc = data_cspc
444 445 self.dataOut.data_dc = data_dc
445 446 self.dataOut.heightList = heightList
446 447 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
447 448
448 449 return 1
449 450
450 451 def selectFFTsByIndex(self, minIndex, maxIndex):
451 452 """
452 453
453 454 """
454 455
455 456 if (minIndex < 0) or (minIndex > maxIndex):
456 457 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex))
457 458
458 459 if (maxIndex >= self.dataOut.nProfiles):
459 460 maxIndex = self.dataOut.nProfiles-1
460 461
461 462 #Spectra
462 463 data_spc = self.dataOut.data_spc[:,minIndex:maxIndex+1,:]
463 464
464 465 data_cspc = None
465 466 if self.dataOut.data_cspc is not None:
466 467 data_cspc = self.dataOut.data_cspc[:,minIndex:maxIndex+1,:]
467 468
468 469 data_dc = None
469 470 if self.dataOut.data_dc is not None:
470 471 data_dc = self.dataOut.data_dc[minIndex:maxIndex+1,:]
471 472
472 473 self.dataOut.data_spc = data_spc
473 474 self.dataOut.data_cspc = data_cspc
474 475 self.dataOut.data_dc = data_dc
475 476
476 477 self.dataOut.ippSeconds = self.dataOut.ippSeconds*(self.dataOut.nFFTPoints / numpy.shape(data_cspc)[1])
477 478 self.dataOut.nFFTPoints = numpy.shape(data_cspc)[1]
478 479 self.dataOut.profilesPerBlock = numpy.shape(data_cspc)[1]
479 480
480 481 return 1
481 482
482 483
483 484
484 485 def selectHeightsByIndex(self, minIndex, maxIndex):
485 486 """
486 487 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
487 488 minIndex <= index <= maxIndex
488 489
489 490 Input:
490 491 minIndex : valor de indice minimo de altura a considerar
491 492 maxIndex : valor de indice maximo de altura a considerar
492 493
493 494 Affected:
494 495 self.dataOut.data_spc
495 496 self.dataOut.data_cspc
496 497 self.dataOut.data_dc
497 498 self.dataOut.heightList
498 499
499 500 Return:
500 501 1 si el metodo se ejecuto con exito caso contrario devuelve 0
501 502 """
502 503
503 504 if (minIndex < 0) or (minIndex > maxIndex):
504 505 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (
505 506 minIndex, maxIndex))
506 507
507 508 if (maxIndex >= self.dataOut.nHeights):
508 509 maxIndex = self.dataOut.nHeights - 1
509 510
510 511 # Spectra
511 512 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
512 513
513 514 data_cspc = None
514 515 if self.dataOut.data_cspc is not None:
515 516 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
516 517
517 518 data_dc = None
518 519 if self.dataOut.data_dc is not None:
519 520 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
520 521
521 522 self.dataOut.data_spc = data_spc
522 523 self.dataOut.data_cspc = data_cspc
523 524 self.dataOut.data_dc = data_dc
524 525
525 526 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
526 527
527 528 return 1
528 529
529 530 def removeDC(self, mode=2):
530 531 jspectra = self.dataOut.data_spc
531 532 jcspectra = self.dataOut.data_cspc
532 533
533 534 num_chan = jspectra.shape[0]
534 535 num_hei = jspectra.shape[2]
535 536
536 537 if jcspectra is not None:
537 538 jcspectraExist = True
538 539 num_pairs = jcspectra.shape[0]
539 540 else:
540 541 jcspectraExist = False
541 542
542 543 freq_dc = int(jspectra.shape[1] / 2)
543 544 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
544 545 ind_vel = ind_vel.astype(int)
545 546
546 547 if ind_vel[0] < 0:
547 548 ind_vel[list(range(0, 1))] = ind_vel[list(range(0, 1))] + self.num_prof
548 549
549 550 if mode == 1:
550 551 jspectra[:, freq_dc, :] = (
551 552 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
552 553
553 554 if jcspectraExist:
554 555 jcspectra[:, freq_dc, :] = (
555 556 jcspectra[:, ind_vel[1], :] + jcspectra[:, ind_vel[2], :]) / 2
556 557
557 558 if mode == 2:
558 559
559 560 vel = numpy.array([-2, -1, 1, 2])
560 561 xx = numpy.zeros([4, 4])
561 562
562 563 for fil in range(4):
563 564 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
564 565
565 566 xx_inv = numpy.linalg.inv(xx)
566 567 xx_aux = xx_inv[0, :]
567 568
568 569 for ich in range(num_chan):
569 570 yy = jspectra[ich, ind_vel, :]
570 571 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
571 572
572 573 junkid = jspectra[ich, freq_dc, :] <= 0
573 574 cjunkid = sum(junkid)
574 575
575 576 if cjunkid.any():
576 577 jspectra[ich, freq_dc, junkid.nonzero()] = (
577 578 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
578 579
579 580 if jcspectraExist:
580 581 for ip in range(num_pairs):
581 582 yy = jcspectra[ip, ind_vel, :]
582 583 jcspectra[ip, freq_dc, :] = numpy.dot(xx_aux, yy)
583 584
584 585 self.dataOut.data_spc = jspectra
585 586 self.dataOut.data_cspc = jcspectra
586 587
587 588 return 1
588 589
589 590 def removeInterference2(self):
590 591
591 592 cspc = self.dataOut.data_cspc
592 593 spc = self.dataOut.data_spc
593 594 Heights = numpy.arange(cspc.shape[2])
594 595 realCspc = numpy.abs(cspc)
595 596
596 597 for i in range(cspc.shape[0]):
597 598 LinePower= numpy.sum(realCspc[i], axis=0)
598 599 Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)]
599 600 SelectedHeights = Heights[ numpy.where( LinePower < Threshold ) ]
600 601 InterferenceSum = numpy.sum( realCspc[i,:,SelectedHeights], axis=0 )
601 602 InterferenceThresholdMin = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)]
602 603 InterferenceThresholdMax = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.99)]
603 604
604 605
605 606 InterferenceRange = numpy.where( ([InterferenceSum > InterferenceThresholdMin]))# , InterferenceSum < InterferenceThresholdMax]) )
606 607 #InterferenceRange = numpy.where( ([InterferenceRange < InterferenceThresholdMax]))
607 608 if len(InterferenceRange)<int(cspc.shape[1]*0.3):
608 609 cspc[i,InterferenceRange,:] = numpy.NaN
609 610
610 611
611 612
612 613 self.dataOut.data_cspc = cspc
613 614
614 615 def removeInterference(self, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None):
615 616
616 617 jspectra = self.dataOut.data_spc
617 618 jcspectra = self.dataOut.data_cspc
618 619 jnoise = self.dataOut.getNoise()
619 620 num_incoh = self.dataOut.nIncohInt
620 621
621 622 num_channel = jspectra.shape[0]
622 623 num_prof = jspectra.shape[1]
623 624 num_hei = jspectra.shape[2]
624 625
625 626 # hei_interf
626 627 if hei_interf is None:
627 628 count_hei = int(num_hei / 2)
628 629 hei_interf = numpy.asmatrix(list(range(count_hei))) + num_hei - count_hei
629 630 hei_interf = numpy.asarray(hei_interf)[0]
630 631 # nhei_interf
631 632 if (nhei_interf == None):
632 633 nhei_interf = 5
633 634 if (nhei_interf < 1):
634 635 nhei_interf = 1
635 636 if (nhei_interf > count_hei):
636 637 nhei_interf = count_hei
637 638 if (offhei_interf == None):
638 639 offhei_interf = 0
639 640
640 641 ind_hei = list(range(num_hei))
641 642 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
642 643 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
643 644 mask_prof = numpy.asarray(list(range(num_prof)))
644 645 num_mask_prof = mask_prof.size
645 646 comp_mask_prof = [0, num_prof / 2]
646 647
647 648 # noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
648 649 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
649 650 jnoise = numpy.nan
650 651 noise_exist = jnoise[0] < numpy.Inf
651 652
652 653 # Subrutina de Remocion de la Interferencia
653 654 for ich in range(num_channel):
654 655 # Se ordena los espectros segun su potencia (menor a mayor)
655 656 power = jspectra[ich, mask_prof, :]
656 657 power = power[:, hei_interf]
657 658 power = power.sum(axis=0)
658 659 psort = power.ravel().argsort()
659 660
660 661 # Se estima la interferencia promedio en los Espectros de Potencia empleando
661 662 junkspc_interf = jspectra[ich, :, hei_interf[psort[list(range(
662 663 offhei_interf, nhei_interf + offhei_interf))]]]
663 664
664 665 if noise_exist:
665 666 # tmp_noise = jnoise[ich] / num_prof
666 667 tmp_noise = jnoise[ich]
667 668 junkspc_interf = junkspc_interf - tmp_noise
668 669 #junkspc_interf[:,comp_mask_prof] = 0
669 670
670 671 jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf
671 672 jspc_interf = jspc_interf.transpose()
672 673 # Calculando el espectro de interferencia promedio
673 674 noiseid = numpy.where(
674 675 jspc_interf <= tmp_noise / numpy.sqrt(num_incoh))
675 676 noiseid = noiseid[0]
676 677 cnoiseid = noiseid.size
677 678 interfid = numpy.where(
678 679 jspc_interf > tmp_noise / numpy.sqrt(num_incoh))
679 680 interfid = interfid[0]
680 681 cinterfid = interfid.size
681 682
682 683 if (cnoiseid > 0):
683 684 jspc_interf[noiseid] = 0
684 685
685 686 # Expandiendo los perfiles a limpiar
686 687 if (cinterfid > 0):
687 688 new_interfid = (
688 689 numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof) % num_prof
689 690 new_interfid = numpy.asarray(new_interfid)
690 691 new_interfid = {x for x in new_interfid}
691 692 new_interfid = numpy.array(list(new_interfid))
692 693 new_cinterfid = new_interfid.size
693 694 else:
694 695 new_cinterfid = 0
695 696
696 697 for ip in range(new_cinterfid):
697 698 ind = junkspc_interf[:, new_interfid[ip]].ravel().argsort()
698 699 jspc_interf[new_interfid[ip]
699 700 ] = junkspc_interf[ind[nhei_interf // 2], new_interfid[ip]]
700 701
701 702 jspectra[ich, :, ind_hei] = jspectra[ich, :,
702 703 ind_hei] - jspc_interf # Corregir indices
703 704
704 705 # Removiendo la interferencia del punto de mayor interferencia
705 706 ListAux = jspc_interf[mask_prof].tolist()
706 707 maxid = ListAux.index(max(ListAux))
707 708
708 709 if cinterfid > 0:
709 710 for ip in range(cinterfid * (interf == 2) - 1):
710 711 ind = (jspectra[ich, interfid[ip], :] < tmp_noise *
711 712 (1 + 1 / numpy.sqrt(num_incoh))).nonzero()
712 713 cind = len(ind)
713 714
714 715 if (cind > 0):
715 716 jspectra[ich, interfid[ip], ind] = tmp_noise * \
716 717 (1 + (numpy.random.uniform(cind) - 0.5) /
717 718 numpy.sqrt(num_incoh))
718 719
719 720 ind = numpy.array([-2, -1, 1, 2])
720 721 xx = numpy.zeros([4, 4])
721 722
722 723 for id1 in range(4):
723 724 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
724 725
725 726 xx_inv = numpy.linalg.inv(xx)
726 727 xx = xx_inv[:, 0]
727 728 ind = (ind + maxid + num_mask_prof) % num_mask_prof
728 729 yy = jspectra[ich, mask_prof[ind], :]
729 730 jspectra[ich, mask_prof[maxid], :] = numpy.dot(
730 731 yy.transpose(), xx)
731 732
732 733 indAux = (jspectra[ich, :, :] < tmp_noise *
733 734 (1 - 1 / numpy.sqrt(num_incoh))).nonzero()
734 735 jspectra[ich, indAux[0], indAux[1]] = tmp_noise * \
735 736 (1 - 1 / numpy.sqrt(num_incoh))
736 737
737 738 # Remocion de Interferencia en el Cross Spectra
738 739 if jcspectra is None:
739 740 return jspectra, jcspectra
740 741 num_pairs = int(jcspectra.size / (num_prof * num_hei))
741 742 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
742 743
743 744 for ip in range(num_pairs):
744 745
745 746 #-------------------------------------------
746 747
747 748 cspower = numpy.abs(jcspectra[ip, mask_prof, :])
748 749 cspower = cspower[:, hei_interf]
749 750 cspower = cspower.sum(axis=0)
750 751
751 752 cspsort = cspower.ravel().argsort()
752 753 junkcspc_interf = jcspectra[ip, :, hei_interf[cspsort[list(range(
753 754 offhei_interf, nhei_interf + offhei_interf))]]]
754 755 junkcspc_interf = junkcspc_interf.transpose()
755 756 jcspc_interf = junkcspc_interf.sum(axis=1) / nhei_interf
756 757
757 758 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
758 759
759 760 median_real = int(numpy.median(numpy.real(
760 761 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
761 762 median_imag = int(numpy.median(numpy.imag(
762 763 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
763 764 comp_mask_prof = [int(e) for e in comp_mask_prof]
764 765 junkcspc_interf[comp_mask_prof, :] = numpy.complex(
765 766 median_real, median_imag)
766 767
767 768 for iprof in range(num_prof):
768 769 ind = numpy.abs(junkcspc_interf[iprof, :]).ravel().argsort()
769 770 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf // 2]]
770 771
771 772 # Removiendo la Interferencia
772 773 jcspectra[ip, :, ind_hei] = jcspectra[ip,
773 774 :, ind_hei] - jcspc_interf
774 775
775 776 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
776 777 maxid = ListAux.index(max(ListAux))
777 778
778 779 ind = numpy.array([-2, -1, 1, 2])
779 780 xx = numpy.zeros([4, 4])
780 781
781 782 for id1 in range(4):
782 783 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
783 784
784 785 xx_inv = numpy.linalg.inv(xx)
785 786 xx = xx_inv[:, 0]
786 787
787 788 ind = (ind + maxid + num_mask_prof) % num_mask_prof
788 789 yy = jcspectra[ip, mask_prof[ind], :]
789 790 jcspectra[ip, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
790 791
791 792 # Guardar Resultados
792 793 self.dataOut.data_spc = jspectra
793 794 self.dataOut.data_cspc = jcspectra
794 795
795 796 return 1
796 797
797 798 def setRadarFrequency(self, frequency=None):
798 799
799 800 if frequency != None:
800 801 self.dataOut.frequency = frequency
801 802
802 803 return 1
803 804
804 805 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
805 806 # validacion de rango
806 807 if minHei == None:
807 808 minHei = self.dataOut.heightList[0]
808 809
809 810 if maxHei == None:
810 811 maxHei = self.dataOut.heightList[-1]
811 812
812 813 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
813 814 print('minHei: %.2f is out of the heights range' % (minHei))
814 815 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
815 816 minHei = self.dataOut.heightList[0]
816 817
817 818 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
818 819 print('maxHei: %.2f is out of the heights range' % (maxHei))
819 820 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
820 821 maxHei = self.dataOut.heightList[-1]
821 822
822 823 # validacion de velocidades
823 824 velrange = self.dataOut.getVelRange(1)
824 825
825 826 if minVel == None:
826 827 minVel = velrange[0]
827 828
828 829 if maxVel == None:
829 830 maxVel = velrange[-1]
830 831
831 832 if (minVel < velrange[0]) or (minVel > maxVel):
832 833 print('minVel: %.2f is out of the velocity range' % (minVel))
833 834 print('minVel is setting to %.2f' % (velrange[0]))
834 835 minVel = velrange[0]
835 836
836 837 if (maxVel > velrange[-1]) or (maxVel < minVel):
837 838 print('maxVel: %.2f is out of the velocity range' % (maxVel))
838 839 print('maxVel is setting to %.2f' % (velrange[-1]))
839 840 maxVel = velrange[-1]
840 841
841 842 # seleccion de indices para rango
842 843 minIndex = 0
843 844 maxIndex = 0
844 845 heights = self.dataOut.heightList
845 846
846 847 inda = numpy.where(heights >= minHei)
847 848 indb = numpy.where(heights <= maxHei)
848 849
849 850 try:
850 851 minIndex = inda[0][0]
851 852 except:
852 853 minIndex = 0
853 854
854 855 try:
855 856 maxIndex = indb[0][-1]
856 857 except:
857 858 maxIndex = len(heights)
858 859
859 860 if (minIndex < 0) or (minIndex > maxIndex):
860 861 raise ValueError("some value in (%d,%d) is not valid" % (
861 862 minIndex, maxIndex))
862 863
863 864 if (maxIndex >= self.dataOut.nHeights):
864 865 maxIndex = self.dataOut.nHeights - 1
865 866
866 867 # seleccion de indices para velocidades
867 868 indminvel = numpy.where(velrange >= minVel)
868 869 indmaxvel = numpy.where(velrange <= maxVel)
869 870 try:
870 871 minIndexVel = indminvel[0][0]
871 872 except:
872 873 minIndexVel = 0
873 874
874 875 try:
875 876 maxIndexVel = indmaxvel[0][-1]
876 877 except:
877 878 maxIndexVel = len(velrange)
878 879
879 880 # seleccion del espectro
880 881 data_spc = self.dataOut.data_spc[:,
881 882 minIndexVel:maxIndexVel + 1, minIndex:maxIndex + 1]
882 883 # estimacion de ruido
883 884 noise = numpy.zeros(self.dataOut.nChannels)
884 885
885 886 for channel in range(self.dataOut.nChannels):
886 887 daux = data_spc[channel, :, :]
887 888 noise[channel] = hildebrand_sekhon(daux, self.dataOut.nIncohInt)
888 889
889 890 self.dataOut.noise_estimation = noise.copy()
890 891
891 892 return 1
892 893
893 894
894 895 class IncohInt(Operation):
895 896
896 897 __profIndex = 0
897 898 __withOverapping = False
898 899
899 900 __byTime = False
900 901 __initime = None
901 902 __lastdatatime = None
902 903 __integrationtime = None
903 904
904 905 __buffer_spc = None
905 906 __buffer_cspc = None
906 907 __buffer_dc = None
907 908
908 909 __dataReady = False
909 910
910 911 __timeInterval = None
911 912
912 913 n = None
913 914
914 915 def __init__(self):
915 916
916 917 Operation.__init__(self)
917 918
918 919 def setup(self, n=None, timeInterval=None, overlapping=False):
919 920 """
920 921 Set the parameters of the integration class.
921 922
922 923 Inputs:
923 924
924 925 n : Number of coherent integrations
925 926 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
926 927 overlapping :
927 928
928 929 """
929 930
930 931 self.__initime = None
931 932 self.__lastdatatime = 0
932 933
933 934 self.__buffer_spc = 0
934 935 self.__buffer_cspc = 0
935 936 self.__buffer_dc = 0
936 937
937 938 self.__profIndex = 0
938 939 self.__dataReady = False
939 940 self.__byTime = False
940 941
941 942 if n is None and timeInterval is None:
942 943 raise ValueError("n or timeInterval should be specified ...")
943 944
944 945 if n is not None:
945 946 self.n = int(n)
946 947 else:
947 948
948 949 self.__integrationtime = int(timeInterval)
949 950 self.n = None
950 951 self.__byTime = True
951 952
952 953 def putData(self, data_spc, data_cspc, data_dc):
953 954 """
954 955 Add a profile to the __buffer_spc and increase in one the __profileIndex
955 956
956 957 """
957 958 print("profIndex: ",self.__profIndex)
958 959 print("data_spc.shape: ",data_spc.shape)
959 960 print("data_spc.shape: ",data_spc[0,0,:])
960 961
961 962 self.__buffer_spc += data_spc
962 963
963 964 if data_cspc is None:
964 965 self.__buffer_cspc = None
965 966 else:
966 967 self.__buffer_cspc += data_cspc
967 968
968 969 if data_dc is None:
969 970 self.__buffer_dc = None
970 971 else:
971 972 self.__buffer_dc += data_dc
972 973
973 974 self.__profIndex += 1
974 975
975 976 return
976 977
977 978 def pushData(self):
978 979 """
979 980 Return the sum of the last profiles and the profiles used in the sum.
980 981
981 982 Affected:
982 983
983 984 self.__profileIndex
984 985
985 986 """
986 987
987 988 data_spc = self.__buffer_spc
988 989 data_cspc = self.__buffer_cspc
989 990 data_dc = self.__buffer_dc
990 991 n = self.__profIndex
991 992
992 993 self.__buffer_spc = 0
993 994 self.__buffer_cspc = 0
994 995 self.__buffer_dc = 0
995 996 self.__profIndex = 0
996 997
997 998 return data_spc, data_cspc, data_dc, n
998 999
999 1000 def byProfiles(self, *args):
1000 1001
1001 1002 self.__dataReady = False
1002 1003 avgdata_spc = None
1003 1004 avgdata_cspc = None
1004 1005 avgdata_dc = None
1005 1006
1006 1007 self.putData(*args)
1007 1008
1008 1009 if self.__profIndex == self.n:
1009 1010
1010 1011 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1011 1012 self.n = n
1012 1013 self.__dataReady = True
1013 1014
1014 1015 return avgdata_spc, avgdata_cspc, avgdata_dc
1015 1016
1016 1017 def byTime(self, datatime, *args):
1017 1018
1018 1019 self.__dataReady = False
1019 1020 avgdata_spc = None
1020 1021 avgdata_cspc = None
1021 1022 avgdata_dc = None
1022 1023
1023 1024 self.putData(*args)
1024 1025
1025 1026 if (datatime - self.__initime) >= self.__integrationtime:
1026 1027 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1027 1028 self.n = n
1028 1029 self.__dataReady = True
1029 1030
1030 1031 return avgdata_spc, avgdata_cspc, avgdata_dc
1031 1032
1032 1033 def integrate(self, datatime, *args):
1033 1034
1034 1035 if self.__profIndex == 0:
1035 1036 self.__initime = datatime
1036 1037
1037 1038 if self.__byTime:
1038 1039 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
1039 1040 datatime, *args)
1040 1041 else:
1041 1042 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
1042 1043
1043 1044 if not self.__dataReady:
1044 1045 return None, None, None, None
1045 1046
1046 1047 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
1047 1048
1048 1049 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
1049 1050 if n == 1:
1050 1051 return
1051 1052
1052 1053 dataOut.flagNoData = True
1053 1054
1054 1055 if not self.isConfig:
1055 1056 self.setup(n, timeInterval, overlapping)
1056 1057 self.isConfig = True
1057 1058
1058 1059 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
1059 1060 dataOut.data_spc,
1060 1061 dataOut.data_cspc,
1061 1062 dataOut.data_dc)
1062 1063
1063 1064 if self.__dataReady:
1064 1065
1065 1066 dataOut.data_spc = avgdata_spc
1066 1067 dataOut.data_cspc = avgdata_cspc
1067 1068 dataOut.data_dc = avgdata_dc
1068 1069 dataOut.nIncohInt *= self.n
1069 1070 dataOut.utctime = avgdatatime
1070 1071 dataOut.flagNoData = False
1071 1072
1072 1073 return dataOut
1073 1074
1074 1075
1075 1076 class PulsePair(Operation):
1076 1077 isConfig = False
1077 1078 __profIndex = 0
1078 1079 __profIndex2 = 0
1079 1080 __initime = None
1080 1081 __lastdatatime = None
1081 1082 __buffer = None
1082 1083 __buffer2 = []
1083 1084 __buffer3 = None
1084 1085 __dataReady = False
1085 1086 n = None
1086 1087
1087 1088 __nch =0
1088 1089 __nProf =0
1089 1090 __nHeis =0
1090 1091
1091 1092 def __init__(self,**kwargs):
1092 1093 Operation.__init__(self,**kwargs)
1093 1094
1094 1095 def setup(self,dataOut,n =None, m = None):
1095 1096
1096 1097 self.__initime = None
1097 1098 self.__lastdatatime = 0
1098 1099 self.__buffer = 0
1099 1100 self.__bufferV = 0
1100 1101 #self.__buffer2 = []
1101 1102 self.__buffer3 = 0
1102 1103 self.__dataReady = False
1103 1104 self.__profIndex = 0
1104 1105 self.__profIndex2 = 0
1105 1106 self.count = 0
1106 1107
1107 1108
1108 1109 self.__nch = dataOut.nChannels
1109 1110 self.__nHeis = dataOut.nHeights
1110 1111 self.__nProf = dataOut.nProfiles
1111 1112 self.__nFFT = dataOut.nFFTPoints
1112 1113 #print("Valores de Ch,Samples,Perfiles,nFFT",self.__nch,self.__nHeis,self.__nProf, self.__nFFT)
1113 1114 #print("EL VALOR DE n es:",n)
1114 1115 if n == None:
1115 1116 raise ValueError("n Should be specified.")
1116 1117
1117 1118 if n != None:
1118 1119 if n<2:
1119 1120 raise ValueError("n Should be greather than 2 ")
1120 1121 self.n = n
1121 1122 if m == None:
1122 1123 m = n
1123 1124 if m != None:
1124 1125 if m<2:
1125 1126 raise ValueError("n Should be greather than 2 ")
1126 1127
1127 1128 self.m = m
1128 1129 self.__buffer2 = numpy.zeros((self.__nch,self.m,self.__nHeis))
1129 1130 self.__bufferV2 = numpy.zeros((self.__nch,self.m,self.__nHeis))
1130 1131
1131 1132
1132 1133
1133 1134 def putData(self,data):
1134 1135 #print("###################################################")
1135 1136 '''
1136 1137 data_tmp = numpy.zeros(self.__nch,self.n,self.__nHeis, dtype= complex)
1137 1138 if self.count < self.__nProf:
1138 1139
1139 1140 for i in range(self.n):
1140 1141 data_tmp[:,i,:] = data[:,i+self.count,:]
1141 1142
1142 1143 self.__buffer = data_tmp*numpy.conjugate(data_tmp)
1143 1144
1144 1145
1145 1146 #####self.__buffer = data*numpy.conjugate(data)
1146 1147 #####self.__bufferV = data[:,(self.__nProf-1):,:]*numpy.conjugate(data[:,1:,:])
1147 1148
1148 1149 #self.__buffer2.append(numpy.conjugate(data))
1149 1150
1150 1151 #####self.__profIndex = data.shape[1]
1151 1152 self.count = self.count + self.n -1
1152 1153 self.__profIndex = self.n
1153 1154 '''
1154 1155 self.__buffer = data*numpy.conjugate(data)
1155 1156 self.__bufferV = data[:,(self.__nProf-1):,:]*numpy.conjugate(data[:,1:,:])
1156 1157 self.__profIndex = self.n
1158 #print("spcch0",self.__buffer)
1157 1159 return
1158 1160
1159 1161 def pushData(self):
1160 1162
1161 1163 data_I = numpy.zeros((self.__nch,self.__nHeis))
1162 1164 data_IV = numpy.zeros((self.__nch,self.__nHeis))
1163 1165
1164 1166 for i in range(self.__nch):
1165 data_I[i,:] = numpy.sum(numpy.sum(self.__buffer[i],axis=0),axis=0)/self.n
1166 data_IV[i,:] = numpy.sum(numpy.sum(self.__bufferV[i],axis=0),axis=0)/(self.n-1)
1167
1167 data_I[i,:] = numpy.sum(self.__buffer[i],axis=0)/self.n
1168 data_IV[i,:] = numpy.sum(self.__bufferV[i],axis=0)/(self.n-1)
1169 ##print("******")
1170 #print("data_I",data_I[0])
1171 #print(self.__buffer.shape)
1172 #a=numpy.average(self.__buffer,axis=1)
1173 #print("average", a)
1168 1174 n = self.__profIndex
1169 1175 ####data_intensity = numpy.sum(numpy.sum(self.__buffer,axis=0),axis=0)/self.n
1170 1176 #print("data_intensity push data",data_intensity.shape)
1171 1177 #data_velocity = self.__buffer3/(self.n-1)
1172 1178 ####n = self.__profIndex
1173 1179
1174 1180 self.__buffer = 0
1175 1181 self.__buffer3 = 0
1176 1182 self.__profIndex = 0
1177 1183
1178 1184 #return data_intensity,data_velocity,n
1179 1185 return data_I,data_IV,n
1180 1186
1181 1187 def pulsePairbyProfiles(self,data):
1182 1188 self.__dataReady = False
1183 1189 data_intensity = None
1184 1190 data_velocity = None
1185 1191
1186 1192 self.putData(data)
1187 1193
1188 1194 if self.__profIndex == self.n:
1189 1195 #data_intensity,data_velocity,n = self.pushData()
1190 1196 data_intensity,data_velocity,n = self.pushData()
1191 1197 #print(data_intensity.shape)
1192 1198 #print("self.__profIndex2", self.__profIndex2)
1193 1199 if self.__profIndex2 == 0:
1194 1200 #print("PRIMERA VEZ")
1195 1201 #print("self.__buffer2",self.__buffer2)
1196 1202 for i in range(self.__nch):
1197 1203 self.__buffer2[i][self.__profIndex2] = data_intensity[i]
1198 1204 self.__bufferV2[i][self.__profIndex2] = data_velocity[i]
1199 1205 self.__profIndex2 += 1
1200 1206 return None,None
1201 1207
1202 1208 if self.__profIndex2 > 0:
1203 1209 for i in range(self.__nch):
1204 1210 self.__buffer2[i][self.__profIndex2] = data_intensity[i]
1205 1211 self.__bufferV2[i][self.__profIndex2] = data_velocity[i]
1206 1212 #print("Dentro del bucle",self.__buffer2)
1207 1213 self.__profIndex2 += 1
1208 1214 if self.__profIndex2 == self.m :
1209 1215 data_i = self.__buffer2
1210 1216 data_v = self.__bufferV2
1211 1217 #print(data_i.shape)
1212 1218 self.__dataReady = True
1213 1219 self.__profIndex2 = 0
1214 1220 self.__buffer2 = numpy.zeros((self.__nch,self.m,self.__nHeis))
1215 1221 self.__bufferV2 = numpy.zeros((self.__nch,self.m,self.__nHeis))
1216 1222 return data_i,data_v
1217 1223 return None,None
1218 1224
1219 1225 def pulsePairOp(self,data,datatime=None):
1220 1226 if self.__initime == None:
1221 1227 self.__initime = datatime
1222 1228
1223 1229 data_intensity,data_velocity = self.pulsePairbyProfiles(data)
1224 1230 self.__lastdatatime = datatime
1225 1231
1226 1232 if data_intensity is None:
1227 1233 return None,None,None
1228 1234
1229 1235 avgdatatime = self.__initime
1230 1236 self.__initime = datatime
1231 1237
1232 1238 return data_intensity,data_velocity,avgdatatime
1233 1239
1234 1240 def run(self,dataOut,n =None,m=None):
1235 1241
1236 1242 if not self.isConfig:
1237 1243 self.setup(dataOut = dataOut, n = n, m = m)
1238 1244 self.isConfig = True
1239 1245
1240 1246 data_intensity,data_velocity,avgdatatime = self.pulsePairOp(dataOut.data_wr,dataOut.utctime)
1241 1247 dataOut.flagNoData = True
1242 1248
1243 1249 if self.__dataReady:
1244 1250 #print(" DATA " , data_intensity.shape)
1245 1251 #dataOut.data = numpy.array([data_intensity])#aqui amigo revisa
1246 1252 #tmp = numpy.zeros([1,data_intensity.shape[0],data_intensity.shape[1]])
1247 1253 #tmp[0] = data_intensity
1248 1254 dataOut.data = data_intensity
1249 1255 dataOut.data_velocity = data_velocity
1250 1256 #dataOut.data = tmp
1251 1257 #print(" DATA " , dataOut.data.shape)
1252 1258 dataOut.nIncohInt *= self.n
1253 1259 dataOut.nProfiles = self.m
1254 1260 dataOut.nFFTPoints = self.m
1255 1261 #dataOut.data_intensity = data_intensity
1256 1262 dataOut.PRFbyAngle = self.n
1257 1263 dataOut.utctime = avgdatatime
1258 1264 dataOut.flagNoData = False
1259 1265 #####print("TIEMPO: ",dataOut.utctime)
1260 1266 return dataOut
General Comments 0
You need to be logged in to leave comments. Login now