##// END OF EJS Templates
schain
RSS
-Parameters Plot corrected...
Julio Valdez -
r832:2ef6a65d9dbc
parent child
Show More
Show file before | Show file after | Hide comments
@@ -1,1169 +1,1183
1 1 '''
2 2
3 3 $Author: murco $
4 4 $Id: JROData.py 173 2012-11-20 15:06:21Z murco $
5 5 '''
6 6
7 7 import copy
8 8 import numpy
9 9 import datetime
10 10
11 11 from jroheaderIO import SystemHeader, RadarControllerHeader
12 12
13 13 def getNumpyDtype(dataTypeCode):
14 14
15 15 if dataTypeCode == 0:
16 16 numpyDtype = numpy.dtype([('real','<i1'),('imag','<i1')])
17 17 elif dataTypeCode == 1:
18 18 numpyDtype = numpy.dtype([('real','<i2'),('imag','<i2')])
19 19 elif dataTypeCode == 2:
20 20 numpyDtype = numpy.dtype([('real','<i4'),('imag','<i4')])
21 21 elif dataTypeCode == 3:
22 22 numpyDtype = numpy.dtype([('real','<i8'),('imag','<i8')])
23 23 elif dataTypeCode == 4:
24 24 numpyDtype = numpy.dtype([('real','<f4'),('imag','<f4')])
25 25 elif dataTypeCode == 5:
26 26 numpyDtype = numpy.dtype([('real','<f8'),('imag','<f8')])
27 27 else:
28 28 raise ValueError, 'dataTypeCode was not defined'
29 29
30 30 return numpyDtype
31 31
32 32 def getDataTypeCode(numpyDtype):
33 33
34 34 if numpyDtype == numpy.dtype([('real','<i1'),('imag','<i1')]):
35 35 datatype = 0
36 36 elif numpyDtype == numpy.dtype([('real','<i2'),('imag','<i2')]):
37 37 datatype = 1
38 38 elif numpyDtype == numpy.dtype([('real','<i4'),('imag','<i4')]):
39 39 datatype = 2
40 40 elif numpyDtype == numpy.dtype([('real','<i8'),('imag','<i8')]):
41 41 datatype = 3
42 42 elif numpyDtype == numpy.dtype([('real','<f4'),('imag','<f4')]):
43 43 datatype = 4
44 44 elif numpyDtype == numpy.dtype([('real','<f8'),('imag','<f8')]):
45 45 datatype = 5
46 46 else:
47 47 datatype = None
48 48
49 49 return datatype
50 50
51 51 def hildebrand_sekhon(data, navg):
52 52 """
53 53 This method is for the objective determination of the noise level in Doppler spectra. This
54 54 implementation technique is based on the fact that the standard deviation of the spectral
55 55 densities is equal to the mean spectral density for white Gaussian noise
56 56
57 57 Inputs:
58 58 Data : heights
59 59 navg : numbers of averages
60 60
61 61 Return:
62 62 -1 : any error
63 63 anoise : noise's level
64 64 """
65 65
66 66 sortdata = numpy.sort(data,axis=None)
67 67 lenOfData = len(sortdata)
68 68 nums_min = lenOfData*0.2
69 69
70 70 if nums_min <= 5:
71 71 nums_min = 5
72 72
73 73 sump = 0.
74 74
75 75 sumq = 0.
76 76
77 77 j = 0
78 78
79 79 cont = 1
80 80
81 81 while((cont==1)and(j<lenOfData)):
82 82
83 83 sump += sortdata[j]
84 84
85 85 sumq += sortdata[j]**2
86 86
87 87 if j > nums_min:
88 88 rtest = float(j)/(j-1) + 1.0/navg
89 89 if ((sumq*j) > (rtest*sump**2)):
90 90 j = j - 1
91 91 sump = sump - sortdata[j]
92 92 sumq = sumq - sortdata[j]**2
93 93 cont = 0
94 94
95 95 j += 1
96 96
97 97 lnoise = sump /j
98 98 # stdv = numpy.sqrt((sumq - lnoise**2)/(j - 1))
99 99 return lnoise
100 100
101 101 class Beam:
102 102 def __init__(self):
103 103 self.codeList = []
104 104 self.azimuthList = []
105 105 self.zenithList = []
106 106
107 107 class GenericData(object):
108 108
109 109 flagNoData = True
110 110
111 111 def __init__(self):
112 112
113 113 raise NotImplementedError
114 114
115 115 def copy(self, inputObj=None):
116 116
117 117 if inputObj == None:
118 118 return copy.deepcopy(self)
119 119
120 120 for key in inputObj.__dict__.keys():
121 121
122 122 attribute = inputObj.__dict__[key]
123 123
124 124 #If this attribute is a tuple or list
125 125 if type(inputObj.__dict__[key]) in (tuple, list):
126 126 self.__dict__[key] = attribute[:]
127 127 continue
128 128
129 129 #If this attribute is another object or instance
130 130 if hasattr(attribute, '__dict__'):
131 131 self.__dict__[key] = attribute.copy()
132 132 continue
133 133
134 134 self.__dict__[key] = inputObj.__dict__[key]
135 135
136 136 def deepcopy(self):
137 137
138 138 return copy.deepcopy(self)
139 139
140 140 def isEmpty(self):
141 141
142 142 return self.flagNoData
143 143
144 144 class JROData(GenericData):
145 145
146 146 # m_BasicHeader = BasicHeader()
147 147 # m_ProcessingHeader = ProcessingHeader()
148 148
149 149 systemHeaderObj = SystemHeader()
150 150
151 151 radarControllerHeaderObj = RadarControllerHeader()
152 152
153 153 # data = None
154 154
155 155 type = None
156 156
157 157 datatype = None #dtype but in string
158 158
159 159 # dtype = None
160 160
161 161 # nChannels = None
162 162
163 163 # nHeights = None
164 164
165 165 nProfiles = None
166 166
167 167 heightList = None
168 168
169 169 channelList = None
170 170
171 171 flagDiscontinuousBlock = False
172 172
173 173 useLocalTime = False
174 174
175 175 utctime = None
176 176
177 177 timeZone = None
178 178
179 179 dstFlag = None
180 180
181 181 errorCount = None
182 182
183 183 blocksize = None
184 184
185 185 # nCode = None
186 186 #
187 187 # nBaud = None
188 188 #
189 189 # code = None
190 190
191 191 flagDecodeData = False #asumo q la data no esta decodificada
192 192
193 193 flagDeflipData = False #asumo q la data no esta sin flip
194 194
195 195 flagShiftFFT = False
196 196
197 197 # ippSeconds = None
198 198
199 199 # timeInterval = None
200 200
201 201 nCohInt = None
202 202
203 203 # noise = None
204 204
205 205 windowOfFilter = 1
206 206
207 207 #Speed of ligth
208 208 C = 3e8
209 209
210 210 frequency = 49.92e6
211 211
212 212 realtime = False
213 213
214 214 beacon_heiIndexList = None
215 215
216 216 last_block = None
217 217
218 218 blocknow = None
219 219
220 220 azimuth = None
221 221
222 222 zenith = None
223 223
224 224 beam = Beam()
225 225
226 226 profileIndex = None
227 227
228 228 def __init__(self):
229 229
230 230 raise NotImplementedError
231 231
232 232 def getNoise(self):
233 233
234 234 raise NotImplementedError
235 235
236 236 def getNChannels(self):
237 237
238 238 return len(self.channelList)
239 239
240 240 def getChannelIndexList(self):
241 241
242 242 return range(self.nChannels)
243 243
244 244 def getNHeights(self):
245 245
246 246 return len(self.heightList)
247 247
248 248 def getHeiRange(self, extrapoints=0):
249 249
250 250 heis = self.heightList
251 251 # deltah = self.heightList[1] - self.heightList[0]
252 252 #
253 253 # heis.append(self.heightList[-1])
254 254
255 255 return heis
256 256
257 257 def getDeltaH(self):
258 258
259 259 delta = self.heightList[1] - self.heightList[0]
260 260
261 261 return delta
262 262
263 263 def getltctime(self):
264 264
265 265 if self.useLocalTime:
266 266 return self.utctime - self.timeZone*60
267 267
268 268 return self.utctime
269 269
270 270 def getDatatime(self):
271 271
272 272 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
273 273 return datatimeValue
274 274
275 275 def getTimeRange(self):
276 276
277 277 datatime = []
278 278
279 279 datatime.append(self.ltctime)
280 280 datatime.append(self.ltctime + self.timeInterval+1)
281 281
282 282 datatime = numpy.array(datatime)
283 283
284 284 return datatime
285 285
286 286 def getFmaxTimeResponse(self):
287 287
288 288 period = (10**-6)*self.getDeltaH()/(0.15)
289 289
290 290 PRF = 1./(period * self.nCohInt)
291 291
292 292 fmax = PRF
293 293
294 294 return fmax
295 295
296 296 def getFmax(self):
297 297
298 298 PRF = 1./(self.ippSeconds * self.nCohInt)
299 299
300 300 fmax = PRF
301 301
302 302 return fmax
303 303
304 304 def getVmax(self):
305 305
306 306 _lambda = self.C/self.frequency
307 307
308 308 vmax = self.getFmax() * _lambda
309 309
310 310 return vmax
311 311
312 312 def get_ippSeconds(self):
313 313 '''
314 314 '''
315 315 return self.radarControllerHeaderObj.ippSeconds
316 316
317 317 def set_ippSeconds(self, ippSeconds):
318 318 '''
319 319 '''
320 320
321 321 self.radarControllerHeaderObj.ippSeconds = ippSeconds
322 322
323 323 return
324 324
325 325 def get_dtype(self):
326 326 '''
327 327 '''
328 328 return getNumpyDtype(self.datatype)
329 329
330 330 def set_dtype(self, numpyDtype):
331 331 '''
332 332 '''
333 333
334 334 self.datatype = getDataTypeCode(numpyDtype)
335 335
336 336 def get_code(self):
337 337 '''
338 338 '''
339 339 return self.radarControllerHeaderObj.code
340 340
341 341 def set_code(self, code):
342 342 '''
343 343 '''
344 344 self.radarControllerHeaderObj.code = code
345 345
346 346 return
347 347
348 348 def get_ncode(self):
349 349 '''
350 350 '''
351 351 return self.radarControllerHeaderObj.nCode
352 352
353 353 def set_ncode(self, nCode):
354 354 '''
355 355 '''
356 356 self.radarControllerHeaderObj.nCode = nCode
357 357
358 358 return
359 359
360 360 def get_nbaud(self):
361 361 '''
362 362 '''
363 363 return self.radarControllerHeaderObj.nBaud
364 364
365 365 def set_nbaud(self, nBaud):
366 366 '''
367 367 '''
368 368 self.radarControllerHeaderObj.nBaud = nBaud
369 369
370 370 return
371 371
372 372 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
373 373 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
374 374 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
375 375 #noise = property(getNoise, "I'm the 'nHeights' property.")
376 376 datatime = property(getDatatime, "I'm the 'datatime' property")
377 377 ltctime = property(getltctime, "I'm the 'ltctime' property")
378 378 ippSeconds = property(get_ippSeconds, set_ippSeconds)
379 379 dtype = property(get_dtype, set_dtype)
380 380 # timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
381 381 code = property(get_code, set_code)
382 382 nCode = property(get_ncode, set_ncode)
383 383 nBaud = property(get_nbaud, set_nbaud)
384 384
385 385 class Voltage(JROData):
386 386
387 387 #data es un numpy array de 2 dmensiones (canales, alturas)
388 388 data = None
389 389
390 390 def __init__(self):
391 391 '''
392 392 Constructor
393 393 '''
394 394
395 395 self.useLocalTime = True
396 396
397 397 self.radarControllerHeaderObj = RadarControllerHeader()
398 398
399 399 self.systemHeaderObj = SystemHeader()
400 400
401 401 self.type = "Voltage"
402 402
403 403 self.data = None
404 404
405 405 # self.dtype = None
406 406
407 407 # self.nChannels = 0
408 408
409 409 # self.nHeights = 0
410 410
411 411 self.nProfiles = None
412 412
413 413 self.heightList = None
414 414
415 415 self.channelList = None
416 416
417 417 # self.channelIndexList = None
418 418
419 419 self.flagNoData = True
420 420
421 421 self.flagDiscontinuousBlock = False
422 422
423 423 self.utctime = None
424 424
425 425 self.timeZone = None
426 426
427 427 self.dstFlag = None
428 428
429 429 self.errorCount = None
430 430
431 431 self.nCohInt = None
432 432
433 433 self.blocksize = None
434 434
435 435 self.flagDecodeData = False #asumo q la data no esta decodificada
436 436
437 437 self.flagDeflipData = False #asumo q la data no esta sin flip
438 438
439 439 self.flagShiftFFT = False
440 440
441 441 self.flagDataAsBlock = False #Asumo que la data es leida perfil a perfil
442 442
443 443 self.profileIndex = 0
444 444
445 445 def getNoisebyHildebrand(self, channel = None):
446 446 """
447 447 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
448 448
449 449 Return:
450 450 noiselevel
451 451 """
452 452
453 453 if channel != None:
454 454 data = self.data[channel]
455 455 nChannels = 1
456 456 else:
457 457 data = self.data
458 458 nChannels = self.nChannels
459 459
460 460 noise = numpy.zeros(nChannels)
461 461 power = data * numpy.conjugate(data)
462 462
463 463 for thisChannel in range(nChannels):
464 464 if nChannels == 1:
465 465 daux = power[:].real
466 466 else:
467 467 daux = power[thisChannel,:].real
468 468 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
469 469
470 470 return noise
471 471
472 472 def getNoise(self, type = 1, channel = None):
473 473
474 474 if type == 1:
475 475 noise = self.getNoisebyHildebrand(channel)
476 476
477 477 return noise
478 478
479 479 def getPower(self, channel = None):
480 480
481 481 if channel != None:
482 482 data = self.data[channel]
483 483 else:
484 484 data = self.data
485 485
486 486 power = data * numpy.conjugate(data)
487 487
488 488 return 10*numpy.log10(power.real)
489 489
490 490 def getTimeInterval(self):
491 491
492 492 timeInterval = self.ippSeconds * self.nCohInt
493 493
494 494 return timeInterval
495 495
496 496 noise = property(getNoise, "I'm the 'nHeights' property.")
497 497 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
498 498
499 499 class Spectra(JROData):
500 500
501 #data es un numpy array de 2 dmensiones (canales, perfiles, alturas)
501 #data spc es un numpy array de 2 dmensiones (canales, perfiles, alturas)
502 502 data_spc = None
503 503
504 #data es un numpy array de 2 dmensiones (canales, pares, alturas)
504 #data cspc es un numpy array de 2 dmensiones (canales, pares, alturas)
505 505 data_cspc = None
506 506
507 #data es un numpy array de 2 dmensiones (canales, alturas)
507 #data dc es un numpy array de 2 dmensiones (canales, alturas)
508 508 data_dc = None
509 509
510 #data power
511 data_pwr = None
512
510 513 nFFTPoints = None
511 514
512 515 # nPairs = None
513 516
514 517 pairsList = None
515 518
516 519 nIncohInt = None
517 520
518 521 wavelength = None #Necesario para cacular el rango de velocidad desde la frecuencia
519 522
520 523 nCohInt = None #se requiere para determinar el valor de timeInterval
521 524
522 525 ippFactor = None
523 526
524 527 profileIndex = 0
525 528
526 529 plotting = "spectra"
527 530
528 531 def __init__(self):
529 532 '''
530 533 Constructor
531 534 '''
532 535
533 536 self.useLocalTime = True
534 537
535 538 self.radarControllerHeaderObj = RadarControllerHeader()
536 539
537 540 self.systemHeaderObj = SystemHeader()
538 541
539 542 self.type = "Spectra"
540 543
541 544 # self.data = None
542 545
543 546 # self.dtype = None
544 547
545 548 # self.nChannels = 0
546 549
547 550 # self.nHeights = 0
548 551
549 552 self.nProfiles = None
550 553
551 554 self.heightList = None
552 555
553 556 self.channelList = None
554 557
555 558 # self.channelIndexList = None
556 559
557 560 self.pairsList = None
558 561
559 562 self.flagNoData = True
560 563
561 564 self.flagDiscontinuousBlock = False
562 565
563 566 self.utctime = None
564 567
565 568 self.nCohInt = None
566 569
567 570 self.nIncohInt = None
568 571
569 572 self.blocksize = None
570 573
571 574 self.nFFTPoints = None
572 575
573 576 self.wavelength = None
574 577
575 578 self.flagDecodeData = False #asumo q la data no esta decodificada
576 579
577 580 self.flagDeflipData = False #asumo q la data no esta sin flip
578 581
579 582 self.flagShiftFFT = False
580 583
581 584 self.ippFactor = 1
582 585
583 586 #self.noise = None
584 587
585 588 self.beacon_heiIndexList = []
586 589
587 590 self.noise_estimation = None
588 591
589 592
590 593 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
591 594 """
592 595 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
593 596
594 597 Return:
595 598 noiselevel
596 599 """
597 600
598 601 noise = numpy.zeros(self.nChannels)
599 602
600 603 for channel in range(self.nChannels):
601 604 daux = self.data_spc[channel,xmin_index:xmax_index,ymin_index:ymax_index]
602 605 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
603 606
604 607 return noise
605 608
606 609 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
607 610
608 611 if self.noise_estimation is not None:
609 612 return self.noise_estimation #this was estimated by getNoise Operation defined in jroproc_spectra.py
610 613 else:
611 614 noise = self.getNoisebyHildebrand(xmin_index, xmax_index, ymin_index, ymax_index)
612 615 return noise
613 616
614 617 def getFreqRangeTimeResponse(self, extrapoints=0):
615 618
616 619 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints*self.ippFactor)
617 620 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
618 621
619 622 return freqrange
620 623
621 624 def getAcfRange(self, extrapoints=0):
622 625
623 626 deltafreq = 10./(self.getFmax() / (self.nFFTPoints*self.ippFactor))
624 627 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
625 628
626 629 return freqrange
627 630
628 631 def getFreqRange(self, extrapoints=0):
629 632
630 633 deltafreq = self.getFmax() / (self.nFFTPoints*self.ippFactor)
631 634 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
632 635
633 636 return freqrange
634 637
635 638 def getVelRange(self, extrapoints=0):
636 639
637 640 deltav = self.getVmax() / (self.nFFTPoints*self.ippFactor)
638 641 velrange = deltav*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltav/2
639 642
640 643 return velrange
641 644
642 645 def getNPairs(self):
643 646
644 647 return len(self.pairsList)
645 648
646 649 def getPairsIndexList(self):
647 650
648 651 return range(self.nPairs)
649 652
650 653 def getNormFactor(self):
651 654
652 655 pwcode = 1
653 656
654 657 if self.flagDecodeData:
655 658 pwcode = numpy.sum(self.code[0]**2)
656 659 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
657 660 normFactor = self.nProfiles*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
658 661
659 662 return normFactor
660 663
661 664 def getFlagCspc(self):
662 665
663 666 if self.data_cspc is None:
664 667 return True
665 668
666 669 return False
667 670
668 671 def getFlagDc(self):
669 672
670 673 if self.data_dc is None:
671 674 return True
672 675
673 676 return False
674 677
675 678 def getTimeInterval(self):
676 679
677 680 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles
678 681
679 682 return timeInterval
680 683
684 def getPower(self):
685
686 factor = self.normFactor
687 z = self.data_spc/factor
688 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
689 avg = numpy.average(z, axis=1)
690
691 return 10*numpy.log10(avg)
692
681 693 def setValue(self, value):
682 694
683 695 print "This property should not be initialized"
684 696
685 697 return
686 698
687 699 nPairs = property(getNPairs, setValue, "I'm the 'nPairs' property.")
688 700 pairsIndexList = property(getPairsIndexList, setValue, "I'm the 'pairsIndexList' property.")
689 701 normFactor = property(getNormFactor, setValue, "I'm the 'getNormFactor' property.")
690 702 flag_cspc = property(getFlagCspc, setValue)
691 703 flag_dc = property(getFlagDc, setValue)
692 704 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
693 705 timeInterval = property(getTimeInterval, setValue, "I'm the 'timeInterval' property")
694 706
695 707 class SpectraHeis(Spectra):
696 708
697 709 data_spc = None
698 710
699 711 data_cspc = None
700 712
701 713 data_dc = None
702 714
703 715 nFFTPoints = None
704 716
705 717 # nPairs = None
706 718
707 719 pairsList = None
708 720
709 721 nCohInt = None
710 722
711 723 nIncohInt = None
712 724
713 725 def __init__(self):
714 726
715 727 self.radarControllerHeaderObj = RadarControllerHeader()
716 728
717 729 self.systemHeaderObj = SystemHeader()
718 730
719 731 self.type = "SpectraHeis"
720 732
721 733 # self.dtype = None
722 734
723 735 # self.nChannels = 0
724 736
725 737 # self.nHeights = 0
726 738
727 739 self.nProfiles = None
728 740
729 741 self.heightList = None
730 742
731 743 self.channelList = None
732 744
733 745 # self.channelIndexList = None
734 746
735 747 self.flagNoData = True
736 748
737 749 self.flagDiscontinuousBlock = False
738 750
739 751 # self.nPairs = 0
740 752
741 753 self.utctime = None
742 754
743 755 self.blocksize = None
744 756
745 757 self.profileIndex = 0
746 758
747 759 self.nCohInt = 1
748 760
749 761 self.nIncohInt = 1
750 762
751 763 def getNormFactor(self):
752 764 pwcode = 1
753 765 if self.flagDecodeData:
754 766 pwcode = numpy.sum(self.code[0]**2)
755 767
756 768 normFactor = self.nIncohInt*self.nCohInt*pwcode
757 769
758 770 return normFactor
759 771
760 772 def getTimeInterval(self):
761 773
762 774 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
763 775
764 776 return timeInterval
765 777
766 778 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
767 779 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
768 780
769 781 class Fits(JROData):
770 782
771 783 heightList = None
772 784
773 785 channelList = None
774 786
775 787 flagNoData = True
776 788
777 789 flagDiscontinuousBlock = False
778 790
779 791 useLocalTime = False
780 792
781 793 utctime = None
782 794
783 795 timeZone = None
784 796
785 797 # ippSeconds = None
786 798
787 799 # timeInterval = None
788 800
789 801 nCohInt = None
790 802
791 803 nIncohInt = None
792 804
793 805 noise = None
794 806
795 807 windowOfFilter = 1
796 808
797 809 #Speed of ligth
798 810 C = 3e8
799 811
800 812 frequency = 49.92e6
801 813
802 814 realtime = False
803 815
804 816
805 817 def __init__(self):
806 818
807 819 self.type = "Fits"
808 820
809 821 self.nProfiles = None
810 822
811 823 self.heightList = None
812 824
813 825 self.channelList = None
814 826
815 827 # self.channelIndexList = None
816 828
817 829 self.flagNoData = True
818 830
819 831 self.utctime = None
820 832
821 833 self.nCohInt = 1
822 834
823 835 self.nIncohInt = 1
824 836
825 837 self.useLocalTime = True
826 838
827 839 self.profileIndex = 0
828 840
829 841 # self.utctime = None
830 842 # self.timeZone = None
831 843 # self.ltctime = None
832 844 # self.timeInterval = None
833 845 # self.header = None
834 846 # self.data_header = None
835 847 # self.data = None
836 848 # self.datatime = None
837 849 # self.flagNoData = False
838 850 # self.expName = ''
839 851 # self.nChannels = None
840 852 # self.nSamples = None
841 853 # self.dataBlocksPerFile = None
842 854 # self.comments = ''
843 855 #
844 856
845 857
846 858 def getltctime(self):
847 859
848 860 if self.useLocalTime:
849 861 return self.utctime - self.timeZone*60
850 862
851 863 return self.utctime
852 864
853 865 def getDatatime(self):
854 866
855 867 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
856 868 return datatime
857 869
858 870 def getTimeRange(self):
859 871
860 872 datatime = []
861 873
862 874 datatime.append(self.ltctime)
863 875 datatime.append(self.ltctime + self.timeInterval)
864 876
865 877 datatime = numpy.array(datatime)
866 878
867 879 return datatime
868 880
869 881 def getHeiRange(self):
870 882
871 883 heis = self.heightList
872 884
873 885 return heis
874 886
875 887 def getNHeights(self):
876 888
877 889 return len(self.heightList)
878 890
879 891 def getNChannels(self):
880 892
881 893 return len(self.channelList)
882 894
883 895 def getChannelIndexList(self):
884 896
885 897 return range(self.nChannels)
886 898
887 899 def getNoise(self, type = 1):
888 900
889 901 #noise = numpy.zeros(self.nChannels)
890 902
891 903 if type == 1:
892 904 noise = self.getNoisebyHildebrand()
893 905
894 906 if type == 2:
895 907 noise = self.getNoisebySort()
896 908
897 909 if type == 3:
898 910 noise = self.getNoisebyWindow()
899 911
900 912 return noise
901 913
902 914 def getTimeInterval(self):
903 915
904 916 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
905 917
906 918 return timeInterval
907 919
908 920 datatime = property(getDatatime, "I'm the 'datatime' property")
909 921 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
910 922 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
911 923 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
912 924 noise = property(getNoise, "I'm the 'nHeights' property.")
913 925
914 926 ltctime = property(getltctime, "I'm the 'ltctime' property")
915 927 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
916 928
917 929 class Correlation(JROData):
918 930
919 931 noise = None
920 932
921 933 SNR = None
922 934
923 935 pairsAutoCorr = None #Pairs of Autocorrelation
924 936
925 937 #--------------------------------------------------
926 938
927 939 data_corr = None
928 940
929 941 data_volt = None
930 942
931 943 lagT = None # each element value is a profileIndex
932 944
933 945 lagR = None # each element value is in km
934 946
935 947 pairsList = None
936 948
937 949 calculateVelocity = None
938 950
939 951 nPoints = None
940 952
941 953 nAvg = None
942 954
943 955 bufferSize = None
944 956
945 957 def __init__(self):
946 958 '''
947 959 Constructor
948 960 '''
949 961 self.radarControllerHeaderObj = RadarControllerHeader()
950 962
951 963 self.systemHeaderObj = SystemHeader()
952 964
953 965 self.type = "Correlation"
954 966
955 967 self.data = None
956 968
957 969 self.dtype = None
958 970
959 971 self.nProfiles = None
960 972
961 973 self.heightList = None
962 974
963 975 self.channelList = None
964 976
965 977 self.flagNoData = True
966 978
967 979 self.flagDiscontinuousBlock = False
968 980
969 981 self.utctime = None
970 982
971 983 self.timeZone = None
972 984
973 985 self.dstFlag = None
974 986
975 987 self.errorCount = None
976 988
977 989 self.blocksize = None
978 990
979 991 self.flagDecodeData = False #asumo q la data no esta decodificada
980 992
981 993 self.flagDeflipData = False #asumo q la data no esta sin flip
982 994
983 995 self.pairsList = None
984 996
985 997 self.nPoints = None
986 998
987 999 def getLagTRange(self, extrapoints=0):
988 1000
989 1001 lagTRange = self.lagT
990 1002 diff = lagTRange[1] - lagTRange[0]
991 1003 extra = numpy.arange(1,extrapoints + 1)*diff + lagTRange[-1]
992 1004 lagTRange = numpy.hstack((lagTRange, extra))
993 1005
994 1006 return lagTRange
995 1007
996 1008 def getLagRRange(self, extrapoints=0):
997 1009
998 1010 return self.lagR
999 1011
1000 1012 def getPairsList(self):
1001 1013
1002 1014 return self.pairsList
1003 1015
1004 1016 def getCalculateVelocity(self):
1005 1017
1006 1018 return self.calculateVelocity
1007 1019
1008 1020 def getNPoints(self):
1009 1021
1010 1022 return self.nPoints
1011 1023
1012 1024 def getNAvg(self):
1013 1025
1014 1026 return self.nAvg
1015 1027
1016 1028 def getBufferSize(self):
1017 1029
1018 1030 return self.bufferSize
1019 1031
1020 1032 def getPairsAutoCorr(self):
1021 1033 pairsList = self.pairsList
1022 1034 pairsAutoCorr = numpy.zeros(self.nChannels, dtype = 'int')*numpy.nan
1023 1035
1024 1036 for l in range(len(pairsList)):
1025 1037 firstChannel = pairsList[l][0]
1026 1038 secondChannel = pairsList[l][1]
1027 1039
1028 1040 #Obteniendo pares de Autocorrelacion
1029 1041 if firstChannel == secondChannel:
1030 1042 pairsAutoCorr[firstChannel] = int(l)
1031 1043
1032 1044 pairsAutoCorr = pairsAutoCorr.astype(int)
1033 1045
1034 1046 return pairsAutoCorr
1035 1047
1036 1048 def getNoise(self, mode = 2):
1037 1049
1038 1050 indR = numpy.where(self.lagR == 0)[0][0]
1039 1051 indT = numpy.where(self.lagT == 0)[0][0]
1040 1052
1041 1053 jspectra0 = self.data_corr[:,:,indR,:]
1042 1054 jspectra = copy.copy(jspectra0)
1043 1055
1044 1056 num_chan = jspectra.shape[0]
1045 1057 num_hei = jspectra.shape[2]
1046 1058
1047 1059 freq_dc = jspectra.shape[1]/2
1048 1060 ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
1049 1061
1050 1062 if ind_vel[0]<0:
1051 1063 ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
1052 1064
1053 1065 if mode == 1:
1054 1066 jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
1055 1067
1056 1068 if mode == 2:
1057 1069
1058 1070 vel = numpy.array([-2,-1,1,2])
1059 1071 xx = numpy.zeros([4,4])
1060 1072
1061 1073 for fil in range(4):
1062 1074 xx[fil,:] = vel[fil]**numpy.asarray(range(4))
1063 1075
1064 1076 xx_inv = numpy.linalg.inv(xx)
1065 1077 xx_aux = xx_inv[0,:]
1066 1078
1067 1079 for ich in range(num_chan):
1068 1080 yy = jspectra[ich,ind_vel,:]
1069 1081 jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy)
1070 1082
1071 1083 junkid = jspectra[ich,freq_dc,:]<=0
1072 1084 cjunkid = sum(junkid)
1073 1085
1074 1086 if cjunkid.any():
1075 1087 jspectra[ich,freq_dc,junkid.nonzero()] = (jspectra[ich,ind_vel[1],junkid] + jspectra[ich,ind_vel[2],junkid])/2
1076 1088
1077 1089 noise = jspectra0[:,freq_dc,:] - jspectra[:,freq_dc,:]
1078 1090
1079 1091 return noise
1080 1092
1081 1093 def getTimeInterval(self):
1082 1094
1083 1095 timeInterval = self.ippSeconds * self.nCohInt * self.nPoints
1084 1096
1085 1097 return timeInterval
1086 1098
1087 1099 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
1088 1100 # pairsList = property(getPairsList, "I'm the 'pairsList' property.")
1089 1101 # nPoints = property(getNPoints, "I'm the 'nPoints' property.")
1090 1102 calculateVelocity = property(getCalculateVelocity, "I'm the 'calculateVelocity' property.")
1091 1103 nAvg = property(getNAvg, "I'm the 'nAvg' property.")
1092 1104 bufferSize = property(getBufferSize, "I'm the 'bufferSize' property.")
1093 1105
1094 1106
1095 1107 class Parameters(JROData):
1096 1108
1109 experimentInfo = None #Information about the experiment
1110
1097 1111 #Information from previous data
1098 1112
1099 1113 inputUnit = None #Type of data to be processed
1100 1114
1101 1115 operation = None #Type of operation to parametrize
1102 1116
1103 1117 normFactor = None #Normalization Factor
1104 1118
1105 1119 groupList = None #List of Pairs, Groups, etc
1106 1120
1107 1121 #Parameters
1108 1122
1109 1123 data_param = None #Parameters obtained
1110 1124
1111 1125 data_pre = None #Data Pre Parametrization
1112 1126
1113 1127 data_SNR = None #Signal to Noise Ratio
1114 1128
1115 1129 # heightRange = None #Heights
1116 1130
1117 1131 abscissaList = None #Abscissa, can be velocities, lags or time
1118 1132
1119 1133 noise = None #Noise Potency
1120 1134
1121 1135 utctimeInit = None #Initial UTC time
1122 1136
1123 1137 paramInterval = None #Time interval to calculate Parameters in seconds
1124 1138
1125 1139 useLocalTime = True
1126 1140
1127 1141 #Fitting
1128 1142
1129 1143 data_error = None #Error of the estimation
1130 1144
1131 1145 constants = None
1132 1146
1133 1147 library = None
1134 1148
1135 1149 #Output signal
1136 1150
1137 1151 outputInterval = None #Time interval to calculate output signal in seconds
1138 1152
1139 1153 data_output = None #Out signal
1140 1154
1141 1155
1142 1156
1143 1157 def __init__(self):
1144 1158 '''
1145 1159 Constructor
1146 1160 '''
1147 1161 self.radarControllerHeaderObj = RadarControllerHeader()
1148 1162
1149 1163 self.systemHeaderObj = SystemHeader()
1150 1164
1151 1165 self.type = "Parameters"
1152 1166
1153 1167 def getTimeRange1(self, interval):
1154 1168
1155 1169 datatime = []
1156 1170
1157 1171 if self.useLocalTime:
1158 1172 time1 = self.utctimeInit - self.timeZone*60
1159 1173 else:
1160 1174 time1 = self.utctimeInit
1161 1175
1162 1176 # datatime.append(self.utctimeInit)
1163 1177 # datatime.append(self.utctimeInit + self.outputInterval)
1164 1178 datatime.append(time1)
1165 1179 datatime.append(time1 + interval)
1166 1180
1167 1181 datatime = numpy.array(datatime)
1168 1182
1169 1183 return datatime
Show file before | Show file after | Hide comments
@@ -1,1364 +1,1569
1 1 import os
2 2 import datetime
3 3 import numpy
4 4
5 from figure import Figure, isRealtime
5 from figure import Figure, isRealtime, isTimeInHourRange
6 6 from plotting_codes import *
7 7
8 8 class MomentsPlot(Figure):
9 9
10 10 isConfig = None
11 11 __nsubplots = None
12 12
13 13 WIDTHPROF = None
14 14 HEIGHTPROF = None
15 15 PREFIX = 'prm'
16 16
17 17 def __init__(self):
18 18
19 19 self.isConfig = False
20 20 self.__nsubplots = 1
21 21
22 22 self.WIDTH = 280
23 23 self.HEIGHT = 250
24 24 self.WIDTHPROF = 120
25 25 self.HEIGHTPROF = 0
26 26 self.counter_imagwr = 0
27 27
28 28 self.PLOT_CODE = MOMENTS_CODE
29 29
30 30 self.FTP_WEI = None
31 31 self.EXP_CODE = None
32 32 self.SUB_EXP_CODE = None
33 33 self.PLOT_POS = None
34 34
35 35 def getSubplots(self):
36 36
37 37 ncol = int(numpy.sqrt(self.nplots)+0.9)
38 38 nrow = int(self.nplots*1./ncol + 0.9)
39 39
40 40 return nrow, ncol
41 41
42 42 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
43 43
44 44 self.__showprofile = showprofile
45 45 self.nplots = nplots
46 46
47 47 ncolspan = 1
48 48 colspan = 1
49 49 if showprofile:
50 50 ncolspan = 3
51 51 colspan = 2
52 52 self.__nsubplots = 2
53 53
54 54 self.createFigure(id = id,
55 55 wintitle = wintitle,
56 56 widthplot = self.WIDTH + self.WIDTHPROF,
57 57 heightplot = self.HEIGHT + self.HEIGHTPROF,
58 58 show=show)
59 59
60 60 nrow, ncol = self.getSubplots()
61 61
62 62 counter = 0
63 63 for y in range(nrow):
64 64 for x in range(ncol):
65 65
66 66 if counter >= self.nplots:
67 67 break
68 68
69 69 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
70 70
71 71 if showprofile:
72 72 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
73 73
74 74 counter += 1
75 75
76 76 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=True,
77 77 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
78 78 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
79 79 server=None, folder=None, username=None, password=None,
80 80 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False):
81 81
82 82 """
83 83
84 84 Input:
85 85 dataOut :
86 86 id :
87 87 wintitle :
88 88 channelList :
89 89 showProfile :
90 90 xmin : None,
91 91 xmax : None,
92 92 ymin : None,
93 93 ymax : None,
94 94 zmin : None,
95 95 zmax : None
96 96 """
97 97
98 98 if dataOut.flagNoData:
99 99 return None
100 100
101 101 if realtime:
102 102 if not(isRealtime(utcdatatime = dataOut.utctime)):
103 103 print 'Skipping this plot function'
104 104 return
105 105
106 106 if channelList == None:
107 107 channelIndexList = dataOut.channelIndexList
108 108 else:
109 109 channelIndexList = []
110 110 for channel in channelList:
111 111 if channel not in dataOut.channelList:
112 112 raise ValueError, "Channel %d is not in dataOut.channelList"
113 113 channelIndexList.append(dataOut.channelList.index(channel))
114 114
115 115 factor = dataOut.normFactor
116 116 x = dataOut.abscissaList
117 117 y = dataOut.heightList
118 118
119 119 z = dataOut.data_pre[channelIndexList,:,:]/factor
120 120 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
121 121 avg = numpy.average(z, axis=1)
122 122 noise = dataOut.noise/factor
123 123
124 124 zdB = 10*numpy.log10(z)
125 125 avgdB = 10*numpy.log10(avg)
126 126 noisedB = 10*numpy.log10(noise)
127 127
128 128 #thisDatetime = dataOut.datatime
129 129 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
130 130 title = wintitle + " Parameters"
131 131 xlabel = "Velocity (m/s)"
132 132 ylabel = "Range (Km)"
133 133
134 134 update_figfile = False
135 135
136 136 if not self.isConfig:
137 137
138 138 nplots = len(channelIndexList)
139 139
140 140 self.setup(id=id,
141 141 nplots=nplots,
142 142 wintitle=wintitle,
143 143 showprofile=showprofile,
144 144 show=show)
145 145
146 146 if xmin == None: xmin = numpy.nanmin(x)
147 147 if xmax == None: xmax = numpy.nanmax(x)
148 148 if ymin == None: ymin = numpy.nanmin(y)
149 149 if ymax == None: ymax = numpy.nanmax(y)
150 150 if zmin == None: zmin = numpy.nanmin(avgdB)*0.9
151 151 if zmax == None: zmax = numpy.nanmax(avgdB)*0.9
152 152
153 153 self.FTP_WEI = ftp_wei
154 154 self.EXP_CODE = exp_code
155 155 self.SUB_EXP_CODE = sub_exp_code
156 156 self.PLOT_POS = plot_pos
157 157
158 158 self.isConfig = True
159 159 update_figfile = True
160 160
161 161 self.setWinTitle(title)
162 162
163 163 for i in range(self.nplots):
164 164 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
165 165 title = "Channel %d: %4.2fdB: %s" %(dataOut.channelList[i], noisedB[i], str_datetime)
166 166 axes = self.axesList[i*self.__nsubplots]
167 167 axes.pcolor(x, y, zdB[i,:,:],
168 168 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
169 169 xlabel=xlabel, ylabel=ylabel, title=title,
170 170 ticksize=9, cblabel='')
171 171 #Mean Line
172 172 mean = dataOut.data_param[i, 1, :]
173 173 axes.addpline(mean, y, idline=0, color="black", linestyle="solid", lw=1)
174 174
175 175 if self.__showprofile:
176 176 axes = self.axesList[i*self.__nsubplots +1]
177 177 axes.pline(avgdB[i], y,
178 178 xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
179 179 xlabel='dB', ylabel='', title='',
180 180 ytick_visible=False,
181 181 grid='x')
182 182
183 183 noiseline = numpy.repeat(noisedB[i], len(y))
184 184 axes.addpline(noiseline, y, idline=1, color="black", linestyle="dashed", lw=2)
185 185
186 186 self.draw()
187 187
188 188 self.save(figpath=figpath,
189 189 figfile=figfile,
190 190 save=save,
191 191 ftp=ftp,
192 192 wr_period=wr_period,
193 193 thisDatetime=thisDatetime)
194 194
195 195
196 196
197 197 class SkyMapPlot(Figure):
198 198
199 199 __isConfig = None
200 200 __nsubplots = None
201 201
202 202 WIDTHPROF = None
203 203 HEIGHTPROF = None
204 204 PREFIX = 'mmap'
205 205
206 206 def __init__(self):
207 207
208 208 self.isConfig = False
209 209 self.__nsubplots = 1
210 210
211 211 # self.WIDTH = 280
212 212 # self.HEIGHT = 250
213 213 self.WIDTH = 600
214 214 self.HEIGHT = 600
215 215 self.WIDTHPROF = 120
216 216 self.HEIGHTPROF = 0
217 217 self.counter_imagwr = 0
218 218
219 219 self.PLOT_CODE = MSKYMAP_CODE
220 220
221 221 self.FTP_WEI = None
222 222 self.EXP_CODE = None
223 223 self.SUB_EXP_CODE = None
224 224 self.PLOT_POS = None
225 225
226 226 def getSubplots(self):
227 227
228 228 ncol = int(numpy.sqrt(self.nplots)+0.9)
229 229 nrow = int(self.nplots*1./ncol + 0.9)
230 230
231 231 return nrow, ncol
232 232
233 233 def setup(self, id, nplots, wintitle, showprofile=False, show=True):
234 234
235 235 self.__showprofile = showprofile
236 236 self.nplots = nplots
237 237
238 238 ncolspan = 1
239 239 colspan = 1
240 240
241 241 self.createFigure(id = id,
242 242 wintitle = wintitle,
243 243 widthplot = self.WIDTH, #+ self.WIDTHPROF,
244 244 heightplot = self.HEIGHT,# + self.HEIGHTPROF,
245 245 show=show)
246 246
247 247 nrow, ncol = 1,1
248 248 counter = 0
249 249 x = 0
250 250 y = 0
251 251 self.addAxes(1, 1, 0, 0, 1, 1, True)
252 252
253 253 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=False,
254 254 tmin=0, tmax=24, timerange=None,
255 255 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
256 256 server=None, folder=None, username=None, password=None,
257 257 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False):
258 258
259 259 """
260 260
261 261 Input:
262 262 dataOut :
263 263 id :
264 264 wintitle :
265 265 channelList :
266 266 showProfile :
267 267 xmin : None,
268 268 xmax : None,
269 269 ymin : None,
270 270 ymax : None,
271 271 zmin : None,
272 272 zmax : None
273 273 """
274 274
275 275 arrayParameters = dataOut.data_param
276 276 error = arrayParameters[:,-1]
277 277 indValid = numpy.where(error == 0)[0]
278 278 finalMeteor = arrayParameters[indValid,:]
279 279 finalAzimuth = finalMeteor[:,3]
280 280 finalZenith = finalMeteor[:,4]
281 281
282 282 x = finalAzimuth*numpy.pi/180
283 283 y = finalZenith
284 284 x1 = [dataOut.ltctime, dataOut.ltctime]
285 285
286 286 #thisDatetime = dataOut.datatime
287 287 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.ltctime)
288 288 title = wintitle + " Parameters"
289 289 xlabel = "Zonal Zenith Angle (deg) "
290 290 ylabel = "Meridional Zenith Angle (deg)"
291 291 update_figfile = False
292 292
293 293 if not self.isConfig:
294 294
295 295 nplots = 1
296 296
297 297 self.setup(id=id,
298 298 nplots=nplots,
299 299 wintitle=wintitle,
300 300 showprofile=showprofile,
301 301 show=show)
302 302
303 303 if self.xmin is None and self.xmax is None:
304 304 self.xmin, self.xmax = self.getTimeLim(x1, tmin, tmax, timerange)
305 305
306 306 if timerange != None:
307 307 self.timerange = timerange
308 308 else:
309 309 self.timerange = self.xmax - self.xmin
310 310
311 311 self.FTP_WEI = ftp_wei
312 312 self.EXP_CODE = exp_code
313 313 self.SUB_EXP_CODE = sub_exp_code
314 314 self.PLOT_POS = plot_pos
315 315 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
316 316 self.firstdate = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
317 317 self.isConfig = True
318 318 update_figfile = True
319 319
320 320 self.setWinTitle(title)
321 321
322 322 i = 0
323 323 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
324 324
325 325 axes = self.axesList[i*self.__nsubplots]
326 326 nevents = axes.x_buffer.shape[0] + x.shape[0]
327 327 title = "Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n" %(self.firstdate,str_datetime,nevents)
328 328 axes.polar(x, y,
329 329 title=title, xlabel=xlabel, ylabel=ylabel,
330 330 ticksize=9, cblabel='')
331 331
332 332 self.draw()
333 333
334 334 self.save(figpath=figpath,
335 335 figfile=figfile,
336 336 save=save,
337 337 ftp=ftp,
338 338 wr_period=wr_period,
339 339 thisDatetime=thisDatetime,
340 340 update_figfile=update_figfile)
341 341
342 342 if dataOut.ltctime >= self.xmax:
343 343 self.isConfigmagwr = wr_period
344 344 self.isConfig = False
345 345 update_figfile = True
346 346 axes.__firsttime = True
347 347 self.xmin += self.timerange
348 348 self.xmax += self.timerange
349 349
350 350
351 351
352 352
353 353 class WindProfilerPlot(Figure):
354 354
355 355 __isConfig = None
356 356 __nsubplots = None
357 357
358 358 WIDTHPROF = None
359 359 HEIGHTPROF = None
360 360 PREFIX = 'wind'
361 361
362 362 def __init__(self):
363 363
364 364 self.timerange = None
365 365 self.isConfig = False
366 366 self.__nsubplots = 1
367 367
368 368 self.WIDTH = 800
369 369 self.HEIGHT = 150
370 370 self.WIDTHPROF = 120
371 371 self.HEIGHTPROF = 0
372 372 self.counter_imagwr = 0
373 373
374 374 self.PLOT_CODE = WIND_CODE
375 375
376 376 self.FTP_WEI = None
377 377 self.EXP_CODE = None
378 378 self.SUB_EXP_CODE = None
379 379 self.PLOT_POS = None
380 380 self.tmin = None
381 381 self.tmax = None
382 382
383 383 self.xmin = None
384 384 self.xmax = None
385 385
386 386 self.figfile = None
387 387
388 388 def getSubplots(self):
389 389
390 390 ncol = 1
391 391 nrow = self.nplots
392 392
393 393 return nrow, ncol
394 394
395 395 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
396 396
397 397 self.__showprofile = showprofile
398 398 self.nplots = nplots
399 399
400 400 ncolspan = 1
401 401 colspan = 1
402 402
403 403 self.createFigure(id = id,
404 404 wintitle = wintitle,
405 405 widthplot = self.WIDTH + self.WIDTHPROF,
406 406 heightplot = self.HEIGHT + self.HEIGHTPROF,
407 407 show=show)
408 408
409 409 nrow, ncol = self.getSubplots()
410 410
411 411 counter = 0
412 412 for y in range(nrow):
413 413 if counter >= self.nplots:
414 414 break
415 415
416 416 self.addAxes(nrow, ncol*ncolspan, y, 0, colspan, 1)
417 417 counter += 1
418 418
419 419 def run(self, dataOut, id, wintitle="", channelList=None, showprofile='False',
420 420 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
421 421 zmax_ver = None, zmin_ver = None, SNRmin = None, SNRmax = None,
422 422 timerange=None, SNRthresh = None,
423 423 save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
424 424 server=None, folder=None, username=None, password=None,
425 425 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
426 426 """
427 427
428 428 Input:
429 429 dataOut :
430 430 id :
431 431 wintitle :
432 432 channelList :
433 433 showProfile :
434 434 xmin : None,
435 435 xmax : None,
436 436 ymin : None,
437 437 ymax : None,
438 438 zmin : None,
439 439 zmax : None
440 440 """
441 441
442 442 # if timerange is not None:
443 443 # self.timerange = timerange
444 444 #
445 445 # tmin = None
446 446 # tmax = None
447 447
448 448 x = dataOut.getTimeRange1(dataOut.outputInterval)
449 449 # y = dataOut.heightList
450 450 y = dataOut.heightList
451 451
452 452 z = dataOut.data_output.copy()
453 453 nplots = z.shape[0] #Number of wind dimensions estimated
454 454 nplotsw = nplots
455 455
456 456 #If there is a SNR function defined
457 457 if dataOut.data_SNR is not None:
458 458 nplots += 1
459 459 SNR = dataOut.data_SNR
460 460 SNRavg = numpy.average(SNR, axis=0)
461 461
462 462 SNRdB = 10*numpy.log10(SNR)
463 463 SNRavgdB = 10*numpy.log10(SNRavg)
464 464
465 465 if SNRthresh == None: SNRthresh = -5.0
466 466 ind = numpy.where(SNRavg < 10**(SNRthresh/10))[0]
467 467
468 468 for i in range(nplotsw):
469 469 z[i,ind] = numpy.nan
470 470
471 471
472 472 # showprofile = False
473 473 # thisDatetime = dataOut.datatime
474 474 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.ltctime)
475 475 title = wintitle + "Wind"
476 476 xlabel = ""
477 477 ylabel = "Range (Km)"
478 478 update_figfile = False
479 479
480 480 if not self.isConfig:
481 481
482 482 self.setup(id=id,
483 483 nplots=nplots,
484 484 wintitle=wintitle,
485 485 showprofile=showprofile,
486 486 show=show)
487 487
488 488 if timerange is not None:
489 489 self.timerange = timerange
490 490
491 491 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
492 492
493 493 if ymin == None: ymin = numpy.nanmin(y)
494 494 if ymax == None: ymax = numpy.nanmax(y)
495 495
496 496 if zmax == None: zmax = numpy.nanmax(abs(z[range(2),:]))
497 497 #if numpy.isnan(zmax): zmax = 50
498 498 if zmin == None: zmin = -zmax
499 499
500 500 if nplotsw == 3:
501 501 if zmax_ver == None: zmax_ver = numpy.nanmax(abs(z[2,:]))
502 502 if zmin_ver == None: zmin_ver = -zmax_ver
503 503
504 504 if dataOut.data_SNR is not None:
505 505 if SNRmin == None: SNRmin = numpy.nanmin(SNRavgdB)
506 506 if SNRmax == None: SNRmax = numpy.nanmax(SNRavgdB)
507 507
508 508
509 509 self.FTP_WEI = ftp_wei
510 510 self.EXP_CODE = exp_code
511 511 self.SUB_EXP_CODE = sub_exp_code
512 512 self.PLOT_POS = plot_pos
513 513
514 514 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
515 515 self.isConfig = True
516 516 self.figfile = figfile
517 517 update_figfile = True
518 518
519 519 self.setWinTitle(title)
520 520
521 521 if ((self.xmax - x[1]) < (x[1]-x[0])):
522 522 x[1] = self.xmax
523 523
524 524 strWind = ['Zonal', 'Meridional', 'Vertical']
525 525 strCb = ['Velocity (m/s)','Velocity (m/s)','Velocity (cm/s)']
526 526 zmaxVector = [zmax, zmax, zmax_ver]
527 527 zminVector = [zmin, zmin, zmin_ver]
528 528 windFactor = [1,1,100]
529 529
530 530 for i in range(nplotsw):
531 531
532 532 title = "%s Wind: %s" %(strWind[i], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
533 533 axes = self.axesList[i*self.__nsubplots]
534 534
535 535 z1 = z[i,:].reshape((1,-1))*windFactor[i]
536 536
537 537 axes.pcolorbuffer(x, y, z1,
538 538 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zminVector[i], zmax=zmaxVector[i],
539 539 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
540 540 ticksize=9, cblabel=strCb[i], cbsize="1%", colormap="RdBu_r" )
541 541
542 542 if dataOut.data_SNR is not None:
543 543 i += 1
544 544 title = "Signal Noise Ratio (SNR): %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
545 545 axes = self.axesList[i*self.__nsubplots]
546 546
547 547 SNRavgdB = SNRavgdB.reshape((1,-1))
548 548
549 549 axes.pcolorbuffer(x, y, SNRavgdB,
550 550 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
551 551 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
552 552 ticksize=9, cblabel='', cbsize="1%", colormap="jet")
553 553
554 554 self.draw()
555 555
556 556 if dataOut.ltctime >= self.xmax:
557 557 self.counter_imagwr = wr_period
558 558 self.isConfig = False
559 559 update_figfile = True
560 560
561 561 self.save(figpath=figpath,
562 562 figfile=figfile,
563 563 save=save,
564 564 ftp=ftp,
565 565 wr_period=wr_period,
566 566 thisDatetime=thisDatetime,
567 567 update_figfile=update_figfile)
568 568
569 class ParametersPlot(Figure):
569 570
571 __isConfig = None
572 __nsubplots = None
570 573
571 class ParametersPlot(Figure):
574 WIDTHPROF = None
575 HEIGHTPROF = None
576 PREFIX = 'param'
577
578 nplots = None
579 nchan = None
580
581 def __init__(self):
582
583 self.timerange = None
584 self.isConfig = False
585 self.__nsubplots = 1
586
587 self.WIDTH = 800
588 self.HEIGHT = 180
589 self.WIDTHPROF = 120
590 self.HEIGHTPROF = 0
591 self.counter_imagwr = 0
592
593 self.PLOT_CODE = RTI_CODE
594
595 self.FTP_WEI = None
596 self.EXP_CODE = None
597 self.SUB_EXP_CODE = None
598 self.PLOT_POS = None
599 self.tmin = None
600 self.tmax = None
601
602 self.xmin = None
603 self.xmax = None
604
605 self.figfile = None
606
607 def getSubplots(self):
608
609 ncol = 1
610 nrow = self.nplots
611
612 return nrow, ncol
613
614 def setup(self, id, nplots, wintitle, show=True):
615
616 self.nplots = nplots
617
618 ncolspan = 1
619 colspan = 1
620
621 self.createFigure(id = id,
622 wintitle = wintitle,
623 widthplot = self.WIDTH + self.WIDTHPROF,
624 heightplot = self.HEIGHT + self.HEIGHTPROF,
625 show=show)
626
627 nrow, ncol = self.getSubplots()
628
629 counter = 0
630 for y in range(nrow):
631 for x in range(ncol):
632
633 if counter >= self.nplots:
634 break
635
636 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
637
638 counter += 1
639
640 def run(self, dataOut, id, wintitle="", channelList=None, paramIndex = 0, colormap=True,
641 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None, timerange=None,
642 showSNR=False, SNRthresh = -numpy.inf, SNRmin=None, SNRmax=None,
643 save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
644 server=None, folder=None, username=None, password=None,
645 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
646
647 """
648
649 Input:
650 dataOut :
651 id :
652 wintitle :
653 channelList :
654 showProfile :
655 xmin : None,
656 xmax : None,
657 ymin : None,
658 ymax : None,
659 zmin : None,
660 zmax : None
661 """
662
663 if colormap:
664 colormap="jet"
665 else:
666 colormap="RdBu_r"
667
668 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
669 return
670
671 if channelList == None:
672 channelIndexList = dataOut.channelIndexList
673 else:
674 channelIndexList = []
675 for channel in channelList:
676 if channel not in dataOut.channelList:
677 raise ValueError, "Channel %d is not in dataOut.channelList"
678 channelIndexList.append(dataOut.channelList.index(channel))
679
680 x = dataOut.getTimeRange1(dataOut.paramInterval)
681 y = dataOut.getHeiRange()
682 z = dataOut.data_param[channelIndexList,paramIndex,:]
683
684 if showSNR:
685 #SNR data
686 SNRarray = dataOut.data_SNR[channelIndexList,:]
687 SNRdB = 10*numpy.log10(SNRarray)
688 ind = numpy.where(SNRdB < SNRthresh)
689 z[ind] = numpy.nan
690
691 thisDatetime = dataOut.datatime
692 # thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
693 title = wintitle + " Parameters Plot" #: %s" %(thisDatetime.strftime("%d-%b-%Y"))
694 xlabel = ""
695 ylabel = "Range (Km)"
696
697 update_figfile = False
698
699 if not self.isConfig:
700
701 nchan = len(channelIndexList)
702 self.nchan = nchan
703 self.plotFact = 1
704 nplots = nchan
705
706 if showSNR:
707 nplots = nchan*2
708 self.plotFact = 2
709 if SNRmin == None: SNRmin = numpy.nanmin(SNRdB)
710 if SNRmax == None: SNRmax = numpy.nanmax(SNRdB)
711
712 self.setup(id=id,
713 nplots=nplots,
714 wintitle=wintitle,
715 show=show)
716
717 if timerange != None:
718 self.timerange = timerange
719
720 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
721
722 if ymin == None: ymin = numpy.nanmin(y)
723 if ymax == None: ymax = numpy.nanmax(y)
724 if zmin == None: zmin = dataOut.abscissaList[0]
725 if zmax == None: zmax = dataOut.abscissaList[-1]
726
727 self.FTP_WEI = ftp_wei
728 self.EXP_CODE = exp_code
729 self.SUB_EXP_CODE = sub_exp_code
730 self.PLOT_POS = plot_pos
731
732 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
733 self.isConfig = True
734 self.figfile = figfile
735 update_figfile = True
736
737 self.setWinTitle(title)
738
739 for i in range(self.nchan):
740 index = channelIndexList[i]
741 title = "Channel %d: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
742 axes = self.axesList[i*self.plotFact]
743 z1 = z[i,:].reshape((1,-1))
744 axes.pcolorbuffer(x, y, z1,
745 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
746 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
747 ticksize=9, cblabel='', cbsize="1%",colormap=colormap)
748
749 if showSNR:
750 title = "Channel %d SNR: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
751 axes = self.axesList[i*self.plotFact + 1]
752 SNRdB1 = SNRdB[i,:].reshape((1,-1))
753 axes.pcolorbuffer(x, y, SNRdB1,
754 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
755 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
756 ticksize=9, cblabel='', cbsize="1%",colormap='jet')
757
758
759 self.draw()
760
761 if dataOut.ltctime >= self.xmax:
762 self.counter_imagwr = wr_period
763 self.isConfig = False
764 update_figfile = True
765
766 self.save(figpath=figpath,
767 figfile=figfile,
768 save=save,
769 ftp=ftp,
770 wr_period=wr_period,
771 thisDatetime=thisDatetime,
772 update_figfile=update_figfile)
773
774
775
776 class Parameters1Plot(Figure):
572 777
573 778 __isConfig = None
574 779 __nsubplots = None
575 780
576 781 WIDTHPROF = None
577 782 HEIGHTPROF = None
578 783 PREFIX = 'prm'
579 784
580 785 def __init__(self):
581 786
582 787 self.timerange = 2*60*60
583 788 self.isConfig = False
584 789 self.__nsubplots = 1
585 790
586 791 self.WIDTH = 800
587 792 self.HEIGHT = 150
588 793 self.WIDTHPROF = 120
589 794 self.HEIGHTPROF = 0
590 795 self.counter_imagwr = 0
591 796
592 797 self.PLOT_CODE = PARMS_CODE
593 798
594 799 self.FTP_WEI = None
595 800 self.EXP_CODE = None
596 801 self.SUB_EXP_CODE = None
597 802 self.PLOT_POS = None
598 803 self.tmin = None
599 804 self.tmax = None
600 805
601 806 self.xmin = None
602 807 self.xmax = None
603 808
604 809 self.figfile = None
605 810
606 811 def getSubplots(self):
607 812
608 813 ncol = 1
609 814 nrow = self.nplots
610 815
611 816 return nrow, ncol
612 817
613 818 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
614 819
615 820 self.__showprofile = showprofile
616 821 self.nplots = nplots
617 822
618 823 ncolspan = 1
619 824 colspan = 1
620 825
621 826 self.createFigure(id = id,
622 827 wintitle = wintitle,
623 828 widthplot = self.WIDTH + self.WIDTHPROF,
624 829 heightplot = self.HEIGHT + self.HEIGHTPROF,
625 830 show=show)
626 831
627 832 nrow, ncol = self.getSubplots()
628 833
629 834 counter = 0
630 835 for y in range(nrow):
631 836 for x in range(ncol):
632 837
633 838 if counter >= self.nplots:
634 839 break
635 840
636 841 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
637 842
638 843 if showprofile:
639 844 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
640 845
641 846 counter += 1
642 847
643 848 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=False,
644 849 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,timerange=None,
645 850 parameterIndex = None, onlyPositive = False,
646 851 SNRthresh = -numpy.inf, SNR = True, SNRmin = None, SNRmax = None, onlySNR = False,
647 852 DOP = True,
648 853 zlabel = "", parameterName = "", parameterObject = "data_param",
649 854 save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
650 855 server=None, folder=None, username=None, password=None,
651 856 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
652 857
653 858 """
654 859
655 860 Input:
656 861 dataOut :
657 862 id :
658 863 wintitle :
659 864 channelList :
660 865 showProfile :
661 866 xmin : None,
662 867 xmax : None,
663 868 ymin : None,
664 869 ymax : None,
665 870 zmin : None,
666 871 zmax : None
667 872 """
668 873
669 874 data_param = getattr(dataOut, parameterObject)
670 875
671 876 if channelList == None:
672 877 channelIndexList = numpy.arange(data_param.shape[0])
673 878 else:
674 879 channelIndexList = numpy.array(channelList)
675 880
676 881 nchan = len(channelIndexList) #Number of channels being plotted
677 882
678 883 if nchan < 1:
679 884 return
680 885
681 886 nGraphsByChannel = 0
682 887
683 888 if SNR:
684 889 nGraphsByChannel += 1
685 890 if DOP:
686 891 nGraphsByChannel += 1
687 892
688 893 if nGraphsByChannel < 1:
689 894 return
690 895
691 896 nplots = nGraphsByChannel*nchan
692 897
693 898 if timerange is not None:
694 899 self.timerange = timerange
695 900
696 901 #tmin = None
697 902 #tmax = None
698 903 if parameterIndex == None:
699 904 parameterIndex = 1
700 905
701 906 x = dataOut.getTimeRange1(dataOut.paramInterval)
702 907 y = dataOut.heightList
703 908 z = data_param[channelIndexList,parameterIndex,:].copy()
704 909
705 910 zRange = dataOut.abscissaList
706 911 # nChannels = z.shape[0] #Number of wind dimensions estimated
707 912 # thisDatetime = dataOut.datatime
708 913
709 914 if dataOut.data_SNR is not None:
710 915 SNRarray = dataOut.data_SNR[channelIndexList,:]
711 916 SNRdB = 10*numpy.log10(SNRarray)
712 917 # SNRavgdB = 10*numpy.log10(SNRavg)
713 918 ind = numpy.where(SNRdB < 10**(SNRthresh/10))
714 919 z[ind] = numpy.nan
715 920
716 921 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
717 922 title = wintitle + " Parameters Plot" #: %s" %(thisDatetime.strftime("%d-%b-%Y"))
718 923 xlabel = ""
719 924 ylabel = "Range (Km)"
720 925
721 926 if (SNR and not onlySNR): nplots = 2*nplots
722 927
723 928 if onlyPositive:
724 929 colormap = "jet"
725 930 zmin = 0
726 931 else: colormap = "RdBu_r"
727 932
728 933 if not self.isConfig:
729 934
730 935 self.setup(id=id,
731 936 nplots=nplots,
732 937 wintitle=wintitle,
733 938 showprofile=showprofile,
734 939 show=show)
735 940
736 941 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
737 942
738 943 if ymin == None: ymin = numpy.nanmin(y)
739 944 if ymax == None: ymax = numpy.nanmax(y)
740 945 if zmin == None: zmin = numpy.nanmin(zRange)
741 946 if zmax == None: zmax = numpy.nanmax(zRange)
742 947
743 948 if SNR:
744 949 if SNRmin == None: SNRmin = numpy.nanmin(SNRdB)
745 950 if SNRmax == None: SNRmax = numpy.nanmax(SNRdB)
746 951
747 952 self.FTP_WEI = ftp_wei
748 953 self.EXP_CODE = exp_code
749 954 self.SUB_EXP_CODE = sub_exp_code
750 955 self.PLOT_POS = plot_pos
751 956
752 957 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
753 958 self.isConfig = True
754 959 self.figfile = figfile
755 960
756 961 self.setWinTitle(title)
757 962
758 963 if ((self.xmax - x[1]) < (x[1]-x[0])):
759 964 x[1] = self.xmax
760 965
761 966 for i in range(nchan):
762 967
763 968 if (SNR and not onlySNR): j = 2*i
764 969 else: j = i
765 970
766 971 j = nGraphsByChannel*i
767 972
768 973 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
769 974 title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith)
770 975
771 976 if not onlySNR:
772 977 axes = self.axesList[j*self.__nsubplots]
773 978 z1 = z[i,:].reshape((1,-1))
774 979 axes.pcolorbuffer(x, y, z1,
775 980 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
776 981 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,colormap=colormap,
777 982 ticksize=9, cblabel=zlabel, cbsize="1%")
778 983
779 984 if DOP:
780 985 title = "%s Channel %d: %s" %(parameterName, channelIndexList[i], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
781 986
782 987 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
783 988 title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith)
784 989 axes = self.axesList[j]
785 990 z1 = z[i,:].reshape((1,-1))
786 991 axes.pcolorbuffer(x, y, z1,
787 992 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
788 993 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,colormap=colormap,
789 994 ticksize=9, cblabel=zlabel, cbsize="1%")
790 995
791 996 if SNR:
792 997 title = "Channel %d Signal Noise Ratio (SNR): %s" %(channelIndexList[i], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
793 998 axes = self.axesList[(j)*self.__nsubplots]
794 999 if not onlySNR:
795 1000 axes = self.axesList[(j + 1)*self.__nsubplots]
796 1001
797 1002 axes = self.axesList[(j + nGraphsByChannel-1)]
798 1003
799 1004 z1 = SNRdB[i,:].reshape((1,-1))
800 1005 axes.pcolorbuffer(x, y, z1,
801 1006 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
802 1007 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,colormap="jet",
803 1008 ticksize=9, cblabel=zlabel, cbsize="1%")
804 1009
805 1010
806 1011
807 1012 self.draw()
808 1013
809 1014 if x[1] >= self.axesList[0].xmax:
810 1015 self.counter_imagwr = wr_period
811 1016 self.isConfig = False
812 1017 self.figfile = None
813 1018
814 1019 self.save(figpath=figpath,
815 1020 figfile=figfile,
816 1021 save=save,
817 1022 ftp=ftp,
818 1023 wr_period=wr_period,
819 1024 thisDatetime=thisDatetime,
820 1025 update_figfile=False)
821 1026
822 1027 class SpectralFittingPlot(Figure):
823 1028
824 1029 __isConfig = None
825 1030 __nsubplots = None
826 1031
827 1032 WIDTHPROF = None
828 1033 HEIGHTPROF = None
829 1034 PREFIX = 'prm'
830 1035
831 1036
832 1037 N = None
833 1038 ippSeconds = None
834 1039
835 1040 def __init__(self):
836 1041 self.isConfig = False
837 1042 self.__nsubplots = 1
838 1043
839 1044 self.PLOT_CODE = SPECFIT_CODE
840 1045
841 1046 self.WIDTH = 450
842 1047 self.HEIGHT = 250
843 1048 self.WIDTHPROF = 0
844 1049 self.HEIGHTPROF = 0
845 1050
846 1051 def getSubplots(self):
847 1052
848 1053 ncol = int(numpy.sqrt(self.nplots)+0.9)
849 1054 nrow = int(self.nplots*1./ncol + 0.9)
850 1055
851 1056 return nrow, ncol
852 1057
853 1058 def setup(self, id, nplots, wintitle, showprofile=False, show=True):
854 1059
855 1060 showprofile = False
856 1061 self.__showprofile = showprofile
857 1062 self.nplots = nplots
858 1063
859 1064 ncolspan = 5
860 1065 colspan = 4
861 1066 if showprofile:
862 1067 ncolspan = 5
863 1068 colspan = 4
864 1069 self.__nsubplots = 2
865 1070
866 1071 self.createFigure(id = id,
867 1072 wintitle = wintitle,
868 1073 widthplot = self.WIDTH + self.WIDTHPROF,
869 1074 heightplot = self.HEIGHT + self.HEIGHTPROF,
870 1075 show=show)
871 1076
872 1077 nrow, ncol = self.getSubplots()
873 1078
874 1079 counter = 0
875 1080 for y in range(nrow):
876 1081 for x in range(ncol):
877 1082
878 1083 if counter >= self.nplots:
879 1084 break
880 1085
881 1086 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
882 1087
883 1088 if showprofile:
884 1089 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
885 1090
886 1091 counter += 1
887 1092
888 1093 def run(self, dataOut, id, cutHeight=None, fit=False, wintitle="", channelList=None, showprofile=True,
889 1094 xmin=None, xmax=None, ymin=None, ymax=None,
890 1095 save=False, figpath='./', figfile=None, show=True):
891 1096
892 1097 """
893 1098
894 1099 Input:
895 1100 dataOut :
896 1101 id :
897 1102 wintitle :
898 1103 channelList :
899 1104 showProfile :
900 1105 xmin : None,
901 1106 xmax : None,
902 1107 zmin : None,
903 1108 zmax : None
904 1109 """
905 1110
906 1111 if cutHeight==None:
907 1112 h=270
908 1113 heightindex = numpy.abs(cutHeight - dataOut.heightList).argmin()
909 1114 cutHeight = dataOut.heightList[heightindex]
910 1115
911 1116 factor = dataOut.normFactor
912 1117 x = dataOut.abscissaList[:-1]
913 1118 #y = dataOut.getHeiRange()
914 1119
915 1120 z = dataOut.data_pre[:,:,heightindex]/factor
916 1121 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
917 1122 avg = numpy.average(z, axis=1)
918 1123 listChannels = z.shape[0]
919 1124
920 1125 #Reconstruct Function
921 1126 if fit==True:
922 1127 groupArray = dataOut.groupList
923 1128 listChannels = groupArray.reshape((groupArray.size))
924 1129 listChannels.sort()
925 1130 spcFitLine = numpy.zeros(z.shape)
926 1131 constants = dataOut.constants
927 1132
928 1133 nGroups = groupArray.shape[0]
929 1134 nChannels = groupArray.shape[1]
930 1135 nProfiles = z.shape[1]
931 1136
932 1137 for f in range(nGroups):
933 1138 groupChann = groupArray[f,:]
934 1139 p = dataOut.data_param[f,:,heightindex]
935 1140 # p = numpy.array([ 89.343967,0.14036615,0.17086219,18.89835291,1.58388365,1.55099167])
936 1141 fitLineAux = dataOut.library.modelFunction(p, constants)*nProfiles
937 1142 fitLineAux = fitLineAux.reshape((nChannels,nProfiles))
938 1143 spcFitLine[groupChann,:] = fitLineAux
939 1144 # spcFitLine = spcFitLine/factor
940 1145
941 1146 z = z[listChannels,:]
942 1147 spcFitLine = spcFitLine[listChannels,:]
943 1148 spcFitLinedB = 10*numpy.log10(spcFitLine)
944 1149
945 1150 zdB = 10*numpy.log10(z)
946 1151 #thisDatetime = dataOut.datatime
947 1152 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
948 1153 title = wintitle + " Doppler Spectra: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
949 1154 xlabel = "Velocity (m/s)"
950 1155 ylabel = "Spectrum"
951 1156
952 1157 if not self.isConfig:
953 1158
954 1159 nplots = listChannels.size
955 1160
956 1161 self.setup(id=id,
957 1162 nplots=nplots,
958 1163 wintitle=wintitle,
959 1164 showprofile=showprofile,
960 1165 show=show)
961 1166
962 1167 if xmin == None: xmin = numpy.nanmin(x)
963 1168 if xmax == None: xmax = numpy.nanmax(x)
964 1169 if ymin == None: ymin = numpy.nanmin(zdB)
965 1170 if ymax == None: ymax = numpy.nanmax(zdB)+2
966 1171
967 1172 self.isConfig = True
968 1173
969 1174 self.setWinTitle(title)
970 1175 for i in range(self.nplots):
971 1176 # title = "Channel %d: %4.2fdB" %(dataOut.channelList[i]+1, noisedB[i])
972 1177 title = "Height %4.1f km\nChannel %d:" %(cutHeight, listChannels[i])
973 1178 axes = self.axesList[i*self.__nsubplots]
974 1179 if fit == False:
975 1180 axes.pline(x, zdB[i,:],
976 1181 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
977 1182 xlabel=xlabel, ylabel=ylabel, title=title
978 1183 )
979 1184 if fit == True:
980 1185 fitline=spcFitLinedB[i,:]
981 1186 y=numpy.vstack([zdB[i,:],fitline] )
982 1187 legendlabels=['Data','Fitting']
983 1188 axes.pmultilineyaxis(x, y,
984 1189 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
985 1190 xlabel=xlabel, ylabel=ylabel, title=title,
986 1191 legendlabels=legendlabels, marker=None,
987 1192 linestyle='solid', grid='both')
988 1193
989 1194 self.draw()
990 1195
991 1196 self.save(figpath=figpath,
992 1197 figfile=figfile,
993 1198 save=save,
994 1199 ftp=ftp,
995 1200 wr_period=wr_period,
996 1201 thisDatetime=thisDatetime)
997 1202
998 1203
999 1204 class EWDriftsPlot(Figure):
1000 1205
1001 1206 __isConfig = None
1002 1207 __nsubplots = None
1003 1208
1004 1209 WIDTHPROF = None
1005 1210 HEIGHTPROF = None
1006 1211 PREFIX = 'drift'
1007 1212
1008 1213 def __init__(self):
1009 1214
1010 1215 self.timerange = 2*60*60
1011 1216 self.isConfig = False
1012 1217 self.__nsubplots = 1
1013 1218
1014 1219 self.WIDTH = 800
1015 1220 self.HEIGHT = 150
1016 1221 self.WIDTHPROF = 120
1017 1222 self.HEIGHTPROF = 0
1018 1223 self.counter_imagwr = 0
1019 1224
1020 1225 self.PLOT_CODE = EWDRIFT_CODE
1021 1226
1022 1227 self.FTP_WEI = None
1023 1228 self.EXP_CODE = None
1024 1229 self.SUB_EXP_CODE = None
1025 1230 self.PLOT_POS = None
1026 1231 self.tmin = None
1027 1232 self.tmax = None
1028 1233
1029 1234 self.xmin = None
1030 1235 self.xmax = None
1031 1236
1032 1237 self.figfile = None
1033 1238
1034 1239 def getSubplots(self):
1035 1240
1036 1241 ncol = 1
1037 1242 nrow = self.nplots
1038 1243
1039 1244 return nrow, ncol
1040 1245
1041 1246 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
1042 1247
1043 1248 self.__showprofile = showprofile
1044 1249 self.nplots = nplots
1045 1250
1046 1251 ncolspan = 1
1047 1252 colspan = 1
1048 1253
1049 1254 self.createFigure(id = id,
1050 1255 wintitle = wintitle,
1051 1256 widthplot = self.WIDTH + self.WIDTHPROF,
1052 1257 heightplot = self.HEIGHT + self.HEIGHTPROF,
1053 1258 show=show)
1054 1259
1055 1260 nrow, ncol = self.getSubplots()
1056 1261
1057 1262 counter = 0
1058 1263 for y in range(nrow):
1059 1264 if counter >= self.nplots:
1060 1265 break
1061 1266
1062 1267 self.addAxes(nrow, ncol*ncolspan, y, 0, colspan, 1)
1063 1268 counter += 1
1064 1269
1065 1270 def run(self, dataOut, id, wintitle="", channelList=None,
1066 1271 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
1067 1272 zmaxVertical = None, zminVertical = None, zmaxZonal = None, zminZonal = None,
1068 1273 timerange=None, SNRthresh = -numpy.inf, SNRmin = None, SNRmax = None, SNR_1 = False,
1069 1274 save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
1070 1275 server=None, folder=None, username=None, password=None,
1071 1276 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1072 1277 """
1073 1278
1074 1279 Input:
1075 1280 dataOut :
1076 1281 id :
1077 1282 wintitle :
1078 1283 channelList :
1079 1284 showProfile :
1080 1285 xmin : None,
1081 1286 xmax : None,
1082 1287 ymin : None,
1083 1288 ymax : None,
1084 1289 zmin : None,
1085 1290 zmax : None
1086 1291 """
1087 1292
1088 1293 if timerange is not None:
1089 1294 self.timerange = timerange
1090 1295
1091 1296 tmin = None
1092 1297 tmax = None
1093 1298
1094 1299 x = dataOut.getTimeRange1(dataOut.outputInterval)
1095 1300 # y = dataOut.heightList
1096 1301 y = dataOut.heightList
1097 1302
1098 1303 z = dataOut.data_output
1099 1304 nplots = z.shape[0] #Number of wind dimensions estimated
1100 1305 nplotsw = nplots
1101 1306
1102 1307 #If there is a SNR function defined
1103 1308 if dataOut.data_SNR is not None:
1104 1309 nplots += 1
1105 1310 SNR = dataOut.data_SNR
1106 1311
1107 1312 if SNR_1:
1108 1313 SNR += 1
1109 1314
1110 1315 SNRavg = numpy.average(SNR, axis=0)
1111 1316
1112 1317 SNRdB = 10*numpy.log10(SNR)
1113 1318 SNRavgdB = 10*numpy.log10(SNRavg)
1114 1319
1115 1320 ind = numpy.where(SNRavg < 10**(SNRthresh/10))[0]
1116 1321
1117 1322 for i in range(nplotsw):
1118 1323 z[i,ind] = numpy.nan
1119 1324
1120 1325
1121 1326 showprofile = False
1122 1327 # thisDatetime = dataOut.datatime
1123 1328 thisDatetime = datetime.datetime.utcfromtimestamp(x[1])
1124 1329 title = wintitle + " EW Drifts"
1125 1330 xlabel = ""
1126 1331 ylabel = "Height (Km)"
1127 1332
1128 1333 if not self.isConfig:
1129 1334
1130 1335 self.setup(id=id,
1131 1336 nplots=nplots,
1132 1337 wintitle=wintitle,
1133 1338 showprofile=showprofile,
1134 1339 show=show)
1135 1340
1136 1341 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
1137 1342
1138 1343 if ymin == None: ymin = numpy.nanmin(y)
1139 1344 if ymax == None: ymax = numpy.nanmax(y)
1140 1345
1141 1346 if zmaxZonal == None: zmaxZonal = numpy.nanmax(abs(z[0,:]))
1142 1347 if zminZonal == None: zminZonal = -zmaxZonal
1143 1348 if zmaxVertical == None: zmaxVertical = numpy.nanmax(abs(z[1,:]))
1144 1349 if zminVertical == None: zminVertical = -zmaxVertical
1145 1350
1146 1351 if dataOut.data_SNR is not None:
1147 1352 if SNRmin == None: SNRmin = numpy.nanmin(SNRavgdB)
1148 1353 if SNRmax == None: SNRmax = numpy.nanmax(SNRavgdB)
1149 1354
1150 1355 self.FTP_WEI = ftp_wei
1151 1356 self.EXP_CODE = exp_code
1152 1357 self.SUB_EXP_CODE = sub_exp_code
1153 1358 self.PLOT_POS = plot_pos
1154 1359
1155 1360 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1156 1361 self.isConfig = True
1157 1362
1158 1363
1159 1364 self.setWinTitle(title)
1160 1365
1161 1366 if ((self.xmax - x[1]) < (x[1]-x[0])):
1162 1367 x[1] = self.xmax
1163 1368
1164 1369 strWind = ['Zonal','Vertical']
1165 1370 strCb = 'Velocity (m/s)'
1166 1371 zmaxVector = [zmaxZonal, zmaxVertical]
1167 1372 zminVector = [zminZonal, zminVertical]
1168 1373
1169 1374 for i in range(nplotsw):
1170 1375
1171 1376 title = "%s Drifts: %s" %(strWind[i], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1172 1377 axes = self.axesList[i*self.__nsubplots]
1173 1378
1174 1379 z1 = z[i,:].reshape((1,-1))
1175 1380
1176 1381 axes.pcolorbuffer(x, y, z1,
1177 1382 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zminVector[i], zmax=zmaxVector[i],
1178 1383 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1179 1384 ticksize=9, cblabel=strCb, cbsize="1%", colormap="RdBu_r")
1180 1385
1181 1386 if dataOut.data_SNR is not None:
1182 1387 i += 1
1183 1388 if SNR_1:
1184 1389 title = "Signal Noise Ratio + 1 (SNR+1): %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1185 1390 else:
1186 1391 title = "Signal Noise Ratio (SNR): %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1187 1392 axes = self.axesList[i*self.__nsubplots]
1188 1393 SNRavgdB = SNRavgdB.reshape((1,-1))
1189 1394
1190 1395 axes.pcolorbuffer(x, y, SNRavgdB,
1191 1396 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
1192 1397 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1193 1398 ticksize=9, cblabel='', cbsize="1%", colormap="jet")
1194 1399
1195 1400 self.draw()
1196 1401
1197 1402 if x[1] >= self.axesList[0].xmax:
1198 1403 self.counter_imagwr = wr_period
1199 1404 self.isConfig = False
1200 1405 self.figfile = None
1201 1406
1202 1407
1203 1408
1204 1409
1205 1410 class PhasePlot(Figure):
1206 1411
1207 1412 __isConfig = None
1208 1413 __nsubplots = None
1209 1414
1210 1415 PREFIX = 'mphase'
1211 1416
1212 1417 def __init__(self):
1213 1418
1214 1419 self.timerange = 24*60*60
1215 1420 self.isConfig = False
1216 1421 self.__nsubplots = 1
1217 1422 self.counter_imagwr = 0
1218 1423 self.WIDTH = 600
1219 1424 self.HEIGHT = 300
1220 1425 self.WIDTHPROF = 120
1221 1426 self.HEIGHTPROF = 0
1222 1427 self.xdata = None
1223 1428 self.ydata = None
1224 1429
1225 1430 self.PLOT_CODE = MPHASE_CODE
1226 1431
1227 1432 self.FTP_WEI = None
1228 1433 self.EXP_CODE = None
1229 1434 self.SUB_EXP_CODE = None
1230 1435 self.PLOT_POS = None
1231 1436
1232 1437
1233 1438 self.filename_phase = None
1234 1439
1235 1440 self.figfile = None
1236 1441
1237 1442 def getSubplots(self):
1238 1443
1239 1444 ncol = 1
1240 1445 nrow = 1
1241 1446
1242 1447 return nrow, ncol
1243 1448
1244 1449 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
1245 1450
1246 1451 self.__showprofile = showprofile
1247 1452 self.nplots = nplots
1248 1453
1249 1454 ncolspan = 7
1250 1455 colspan = 6
1251 1456 self.__nsubplots = 2
1252 1457
1253 1458 self.createFigure(id = id,
1254 1459 wintitle = wintitle,
1255 1460 widthplot = self.WIDTH+self.WIDTHPROF,
1256 1461 heightplot = self.HEIGHT+self.HEIGHTPROF,
1257 1462 show=show)
1258 1463
1259 1464 nrow, ncol = self.getSubplots()
1260 1465
1261 1466 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
1262 1467
1263 1468
1264 1469 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
1265 1470 xmin=None, xmax=None, ymin=None, ymax=None,
1266 1471 timerange=None,
1267 1472 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
1268 1473 server=None, folder=None, username=None, password=None,
1269 1474 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1270 1475
1271 1476
1272 1477 tmin = None
1273 1478 tmax = None
1274 1479 x = dataOut.getTimeRange1(dataOut.outputInterval)
1275 1480 y = dataOut.getHeiRange()
1276 1481
1277 1482
1278 1483 #thisDatetime = dataOut.datatime
1279 1484 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
1280 1485 title = wintitle + " Phase of Beacon Signal" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
1281 1486 xlabel = "Local Time"
1282 1487 ylabel = "Phase"
1283 1488
1284 1489
1285 1490 #phase = numpy.zeros((len(pairsIndexList),len(dataOut.beacon_heiIndexList)))
1286 1491 phase_beacon = dataOut.data_output
1287 1492 update_figfile = False
1288 1493
1289 1494 if not self.isConfig:
1290 1495
1291 1496 self.nplots = phase_beacon.size
1292 1497
1293 1498 self.setup(id=id,
1294 1499 nplots=self.nplots,
1295 1500 wintitle=wintitle,
1296 1501 showprofile=showprofile,
1297 1502 show=show)
1298 1503
1299 1504 if timerange is not None:
1300 1505 self.timerange = timerange
1301 1506
1302 1507 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
1303 1508
1304 1509 if ymin == None: ymin = numpy.nanmin(phase_beacon) - 10.0
1305 1510 if ymax == None: ymax = numpy.nanmax(phase_beacon) + 10.0
1306 1511
1307 1512 self.FTP_WEI = ftp_wei
1308 1513 self.EXP_CODE = exp_code
1309 1514 self.SUB_EXP_CODE = sub_exp_code
1310 1515 self.PLOT_POS = plot_pos
1311 1516
1312 1517 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1313 1518 self.isConfig = True
1314 1519 self.figfile = figfile
1315 1520 self.xdata = numpy.array([])
1316 1521 self.ydata = numpy.array([])
1317 1522
1318 1523 #open file beacon phase
1319 1524 path = '%s%03d' %(self.PREFIX, self.id)
1320 1525 beacon_file = os.path.join(path,'%s.txt'%self.name)
1321 1526 self.filename_phase = os.path.join(figpath,beacon_file)
1322 1527 update_figfile = True
1323 1528
1324 1529
1325 1530 #store data beacon phase
1326 1531 #self.save_data(self.filename_phase, phase_beacon, thisDatetime)
1327 1532
1328 1533 self.setWinTitle(title)
1329 1534
1330 1535
1331 1536 title = "Phase Offset %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1332 1537
1333 1538 legendlabels = ["phase %d"%(chan) for chan in numpy.arange(self.nplots)]
1334 1539
1335 1540 axes = self.axesList[0]
1336 1541
1337 1542 self.xdata = numpy.hstack((self.xdata, x[0:1]))
1338 1543
1339 1544 if len(self.ydata)==0:
1340 1545 self.ydata = phase_beacon.reshape(-1,1)
1341 1546 else:
1342 1547 self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1)))
1343 1548
1344 1549
1345 1550 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
1346 1551 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
1347 1552 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
1348 1553 XAxisAsTime=True, grid='both'
1349 1554 )
1350 1555
1351 1556 self.draw()
1352 1557
1353 1558 if dataOut.ltctime >= self.xmax:
1354 1559 self.counter_imagwr = wr_period
1355 1560 self.isConfig = False
1356 1561 update_figfile = True
1357 1562
1358 1563 self.save(figpath=figpath,
1359 1564 figfile=figfile,
1360 1565 save=save,
1361 1566 ftp=ftp,
1362 1567 wr_period=wr_period,
1363 1568 thisDatetime=thisDatetime,
1364 1569 update_figfile=update_figfile)
Show file before | Show file after | Hide comments
@@ -1,1522 +1,1527
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
13 13 class SpectraPlot(Figure):
14 14
15 15 isConfig = None
16 16 __nsubplots = None
17 17
18 18 WIDTHPROF = None
19 19 HEIGHTPROF = None
20 20 PREFIX = 'spc'
21 21
22 22 def __init__(self):
23 23
24 24 self.isConfig = False
25 25 self.__nsubplots = 1
26 26
27 27 self.WIDTH = 250
28 28 self.HEIGHT = 250
29 29 self.WIDTHPROF = 120
30 30 self.HEIGHTPROF = 0
31 31 self.counter_imagwr = 0
32 32
33 33 self.PLOT_CODE = SPEC_CODE
34 34
35 35 self.FTP_WEI = None
36 36 self.EXP_CODE = None
37 37 self.SUB_EXP_CODE = None
38 38 self.PLOT_POS = None
39 39
40 40 self.__xfilter_ena = False
41 41 self.__yfilter_ena = False
42 42
43 43 def getSubplots(self):
44 44
45 45 ncol = int(numpy.sqrt(self.nplots)+0.9)
46 46 nrow = int(self.nplots*1./ncol + 0.9)
47 47
48 48 return nrow, ncol
49 49
50 50 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
51 51
52 52 self.__showprofile = showprofile
53 53 self.nplots = nplots
54 54
55 55 ncolspan = 1
56 56 colspan = 1
57 57 if showprofile:
58 58 ncolspan = 3
59 59 colspan = 2
60 60 self.__nsubplots = 2
61 61
62 62 self.createFigure(id = id,
63 63 wintitle = wintitle,
64 64 widthplot = self.WIDTH + self.WIDTHPROF,
65 65 heightplot = self.HEIGHT + self.HEIGHTPROF,
66 66 show=show)
67 67
68 68 nrow, ncol = self.getSubplots()
69 69
70 70 counter = 0
71 71 for y in range(nrow):
72 72 for x in range(ncol):
73 73
74 74 if counter >= self.nplots:
75 75 break
76 76
77 77 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
78 78
79 79 if showprofile:
80 80 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
81 81
82 82 counter += 1
83 83
84 84 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=True,
85 85 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
86 86 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
87 87 server=None, folder=None, username=None, password=None,
88 88 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False,
89 89 xaxis="frequency"):
90 90
91 91 """
92 92
93 93 Input:
94 94 dataOut :
95 95 id :
96 96 wintitle :
97 97 channelList :
98 98 showProfile :
99 99 xmin : None,
100 100 xmax : None,
101 101 ymin : None,
102 102 ymax : None,
103 103 zmin : None,
104 104 zmax : None
105 105 """
106 106
107 107 if realtime:
108 108 if not(isRealtime(utcdatatime = dataOut.utctime)):
109 109 print 'Skipping this plot function'
110 110 return
111 111
112 112 if channelList == None:
113 113 channelIndexList = dataOut.channelIndexList
114 114 else:
115 115 channelIndexList = []
116 116 for channel in channelList:
117 117 if channel not in dataOut.channelList:
118 118 raise ValueError, "Channel %d is not in dataOut.channelList" %channel
119 119 channelIndexList.append(dataOut.channelList.index(channel))
120 120
121 121 factor = dataOut.normFactor
122 122
123 123 if xaxis == "frequency":
124 124 x = dataOut.getFreqRange(1)/1000.
125 125 xlabel = "Frequency (kHz)"
126 126
127 127 elif xaxis == "time":
128 128 x = dataOut.getAcfRange(1)
129 129 xlabel = "Time (ms)"
130 130
131 131 else:
132 132 x = dataOut.getVelRange(1)
133 133 xlabel = "Velocity (m/s)"
134 134
135 135 ylabel = "Range (Km)"
136 136
137 137 y = dataOut.getHeiRange()
138 138
139 139 z = dataOut.data_spc/factor
140 140 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
141 141 zdB = 10*numpy.log10(z)
142 142
143 143 avg = numpy.average(z, axis=1)
144 144 avgdB = 10*numpy.log10(avg)
145 145
146 146 noise = dataOut.getNoise()/factor
147 147 noisedB = 10*numpy.log10(noise)
148 148
149 149 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
150 150 title = wintitle + " Spectra"
151 151 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
152 152 title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith)
153 153
154 154 if not self.isConfig:
155 155
156 156 nplots = len(channelIndexList)
157 157
158 158 self.setup(id=id,
159 159 nplots=nplots,
160 160 wintitle=wintitle,
161 161 showprofile=showprofile,
162 162 show=show)
163 163
164 164 if xmin == None: xmin = numpy.nanmin(x)
165 165 if xmax == None: xmax = numpy.nanmax(x)
166 166 if ymin == None: ymin = numpy.nanmin(y)
167 167 if ymax == None: ymax = numpy.nanmax(y)
168 168 if zmin == None: zmin = numpy.floor(numpy.nanmin(noisedB)) - 3
169 169 if zmax == None: zmax = numpy.ceil(numpy.nanmax(avgdB)) + 3
170 170
171 171 self.FTP_WEI = ftp_wei
172 172 self.EXP_CODE = exp_code
173 173 self.SUB_EXP_CODE = sub_exp_code
174 174 self.PLOT_POS = plot_pos
175 175
176 176 self.isConfig = True
177 177
178 178 self.setWinTitle(title)
179 179
180 180 for i in range(self.nplots):
181 181 index = channelIndexList[i]
182 182 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
183 183 title = "Channel %d: %4.2fdB: %s" %(dataOut.channelList[index], noisedB[index], str_datetime)
184 184 if len(dataOut.beam.codeList) != 0:
185 185 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)
186 186
187 187 axes = self.axesList[i*self.__nsubplots]
188 188 axes.pcolor(x, y, zdB[index,:,:],
189 189 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
190 190 xlabel=xlabel, ylabel=ylabel, title=title,
191 191 ticksize=9, cblabel='')
192 192
193 193 if self.__showprofile:
194 194 axes = self.axesList[i*self.__nsubplots +1]
195 195 axes.pline(avgdB[index,:], y,
196 196 xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
197 197 xlabel='dB', ylabel='', title='',
198 198 ytick_visible=False,
199 199 grid='x')
200 200
201 201 noiseline = numpy.repeat(noisedB[index], len(y))
202 202 axes.addpline(noiseline, y, idline=1, color="black", linestyle="dashed", lw=2)
203 203
204 204 self.draw()
205 205
206 206 if figfile == None:
207 207 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
208 208 name = str_datetime
209 209 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
210 210 name = name + '_az' + '_%2.2f'%(dataOut.azimuth) + '_zn' + '_%2.2f'%(dataOut.zenith)
211 211 figfile = self.getFilename(name)
212 212
213 213 self.save(figpath=figpath,
214 214 figfile=figfile,
215 215 save=save,
216 216 ftp=ftp,
217 217 wr_period=wr_period,
218 218 thisDatetime=thisDatetime)
219 219
220 220 class CrossSpectraPlot(Figure):
221 221
222 222 isConfig = None
223 223 __nsubplots = None
224 224
225 225 WIDTH = None
226 226 HEIGHT = None
227 227 WIDTHPROF = None
228 228 HEIGHTPROF = None
229 229 PREFIX = 'cspc'
230 230
231 231 def __init__(self):
232 232
233 233 self.isConfig = False
234 234 self.__nsubplots = 4
235 235 self.counter_imagwr = 0
236 236 self.WIDTH = 250
237 237 self.HEIGHT = 250
238 238 self.WIDTHPROF = 0
239 239 self.HEIGHTPROF = 0
240 240
241 241 self.PLOT_CODE = CROSS_CODE
242 242 self.FTP_WEI = None
243 243 self.EXP_CODE = None
244 244 self.SUB_EXP_CODE = None
245 245 self.PLOT_POS = None
246 246
247 247 def getSubplots(self):
248 248
249 249 ncol = 4
250 250 nrow = self.nplots
251 251
252 252 return nrow, ncol
253 253
254 254 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
255 255
256 256 self.__showprofile = showprofile
257 257 self.nplots = nplots
258 258
259 259 ncolspan = 1
260 260 colspan = 1
261 261
262 262 self.createFigure(id = id,
263 263 wintitle = wintitle,
264 264 widthplot = self.WIDTH + self.WIDTHPROF,
265 265 heightplot = self.HEIGHT + self.HEIGHTPROF,
266 266 show=True)
267 267
268 268 nrow, ncol = self.getSubplots()
269 269
270 270 counter = 0
271 271 for y in range(nrow):
272 272 for x in range(ncol):
273 273 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
274 274
275 275 counter += 1
276 276
277 277 def run(self, dataOut, id, wintitle="", pairsList=None,
278 278 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
279 279 coh_min=None, coh_max=None, phase_min=None, phase_max=None,
280 280 save=False, figpath='./', figfile=None, ftp=False, wr_period=1,
281 281 power_cmap='jet', coherence_cmap='jet', phase_cmap='RdBu_r', show=True,
282 282 server=None, folder=None, username=None, password=None,
283 283 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0,
284 284 xaxis='frequency'):
285 285
286 286 """
287 287
288 288 Input:
289 289 dataOut :
290 290 id :
291 291 wintitle :
292 292 channelList :
293 293 showProfile :
294 294 xmin : None,
295 295 xmax : None,
296 296 ymin : None,
297 297 ymax : None,
298 298 zmin : None,
299 299 zmax : None
300 300 """
301 301
302 302 if pairsList == None:
303 303 pairsIndexList = dataOut.pairsIndexList
304 304 else:
305 305 pairsIndexList = []
306 306 for pair in pairsList:
307 307 if pair not in dataOut.pairsList:
308 308 raise ValueError, "Pair %s is not in dataOut.pairsList" %str(pair)
309 309 pairsIndexList.append(dataOut.pairsList.index(pair))
310 310
311 311 if not pairsIndexList:
312 312 return
313 313
314 314 if len(pairsIndexList) > 4:
315 315 pairsIndexList = pairsIndexList[0:4]
316 316
317 317 factor = dataOut.normFactor
318 318 x = dataOut.getVelRange(1)
319 319 y = dataOut.getHeiRange()
320 320 z = dataOut.data_spc[:,:,:]/factor
321 321 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
322 322
323 323 noise = dataOut.noise/factor
324 324
325 325 zdB = 10*numpy.log10(z)
326 326 noisedB = 10*numpy.log10(noise)
327 327
328 328 if coh_min == None:
329 329 coh_min = 0.0
330 330 if coh_max == None:
331 331 coh_max = 1.0
332 332
333 333 if phase_min == None:
334 334 phase_min = -180
335 335 if phase_max == None:
336 336 phase_max = 180
337 337
338 338 #thisDatetime = dataOut.datatime
339 339 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
340 340 title = wintitle + " Cross-Spectra: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
341 341 # xlabel = "Velocity (m/s)"
342 342 ylabel = "Range (Km)"
343 343
344 344 if xaxis == "frequency":
345 345 x = dataOut.getFreqRange(1)/1000.
346 346 xlabel = "Frequency (kHz)"
347 347
348 348 elif xaxis == "time":
349 349 x = dataOut.getAcfRange(1)
350 350 xlabel = "Time (ms)"
351 351
352 352 else:
353 353 x = dataOut.getVelRange(1)
354 354 xlabel = "Velocity (m/s)"
355 355
356 356 if not self.isConfig:
357 357
358 358 nplots = len(pairsIndexList)
359 359
360 360 self.setup(id=id,
361 361 nplots=nplots,
362 362 wintitle=wintitle,
363 363 showprofile=False,
364 364 show=show)
365 365
366 366 avg = numpy.abs(numpy.average(z, axis=1))
367 367 avgdB = 10*numpy.log10(avg)
368 368
369 369 if xmin == None: xmin = numpy.nanmin(x)
370 370 if xmax == None: xmax = numpy.nanmax(x)
371 371 if ymin == None: ymin = numpy.nanmin(y)
372 372 if ymax == None: ymax = numpy.nanmax(y)
373 373 if zmin == None: zmin = numpy.floor(numpy.nanmin(noisedB)) - 3
374 374 if zmax == None: zmax = numpy.ceil(numpy.nanmax(avgdB)) + 3
375 375
376 376 self.FTP_WEI = ftp_wei
377 377 self.EXP_CODE = exp_code
378 378 self.SUB_EXP_CODE = sub_exp_code
379 379 self.PLOT_POS = plot_pos
380 380
381 381 self.isConfig = True
382 382
383 383 self.setWinTitle(title)
384 384
385 385 for i in range(self.nplots):
386 386 pair = dataOut.pairsList[pairsIndexList[i]]
387 387
388 388 chan_index0 = dataOut.channelList.index(pair[0])
389 389 chan_index1 = dataOut.channelList.index(pair[1])
390 390
391 391 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
392 392 title = "Ch%d: %4.2fdB: %s" %(pair[0], noisedB[chan_index0], str_datetime)
393 393 zdB = 10.*numpy.log10(dataOut.data_spc[chan_index0,:,:]/factor)
394 394 axes0 = self.axesList[i*self.__nsubplots]
395 395 axes0.pcolor(x, y, zdB,
396 396 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
397 397 xlabel=xlabel, ylabel=ylabel, title=title,
398 398 ticksize=9, colormap=power_cmap, cblabel='')
399 399
400 400 title = "Ch%d: %4.2fdB: %s" %(pair[1], noisedB[chan_index1], str_datetime)
401 401 zdB = 10.*numpy.log10(dataOut.data_spc[chan_index1,:,:]/factor)
402 402 axes0 = self.axesList[i*self.__nsubplots+1]
403 403 axes0.pcolor(x, y, zdB,
404 404 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
405 405 xlabel=xlabel, ylabel=ylabel, title=title,
406 406 ticksize=9, colormap=power_cmap, cblabel='')
407 407
408 408 coherenceComplex = dataOut.data_cspc[pairsIndexList[i],:,:]/numpy.sqrt(dataOut.data_spc[chan_index0,:,:]*dataOut.data_spc[chan_index1,:,:])
409 409 coherence = numpy.abs(coherenceComplex)
410 410 # phase = numpy.arctan(-1*coherenceComplex.imag/coherenceComplex.real)*180/numpy.pi
411 411 phase = numpy.arctan2(coherenceComplex.imag, coherenceComplex.real)*180/numpy.pi
412 412
413 413 title = "Coherence Ch%d * Ch%d" %(pair[0], pair[1])
414 414 axes0 = self.axesList[i*self.__nsubplots+2]
415 415 axes0.pcolor(x, y, coherence,
416 416 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=coh_min, zmax=coh_max,
417 417 xlabel=xlabel, ylabel=ylabel, title=title,
418 418 ticksize=9, colormap=coherence_cmap, cblabel='')
419 419
420 420 title = "Phase Ch%d * Ch%d" %(pair[0], pair[1])
421 421 axes0 = self.axesList[i*self.__nsubplots+3]
422 422 axes0.pcolor(x, y, phase,
423 423 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=phase_min, zmax=phase_max,
424 424 xlabel=xlabel, ylabel=ylabel, title=title,
425 425 ticksize=9, colormap=phase_cmap, cblabel='')
426 426
427 427
428 428
429 429 self.draw()
430 430
431 431 self.save(figpath=figpath,
432 432 figfile=figfile,
433 433 save=save,
434 434 ftp=ftp,
435 435 wr_period=wr_period,
436 436 thisDatetime=thisDatetime)
437 437
438 438
439 439 class RTIPlot(Figure):
440 440
441 441 __isConfig = None
442 442 __nsubplots = None
443 443
444 444 WIDTHPROF = None
445 445 HEIGHTPROF = None
446 446 PREFIX = 'rti'
447 447
448 448 def __init__(self):
449 449
450 450 self.timerange = None
451 451 self.isConfig = False
452 452 self.__nsubplots = 1
453 453
454 454 self.WIDTH = 800
455 455 self.HEIGHT = 180
456 456 self.WIDTHPROF = 120
457 457 self.HEIGHTPROF = 0
458 458 self.counter_imagwr = 0
459 459
460 460 self.PLOT_CODE = RTI_CODE
461 461
462 462 self.FTP_WEI = None
463 463 self.EXP_CODE = None
464 464 self.SUB_EXP_CODE = None
465 465 self.PLOT_POS = None
466 466 self.tmin = None
467 467 self.tmax = None
468 468
469 469 self.xmin = None
470 470 self.xmax = None
471 471
472 472 self.figfile = None
473 473
474 474 def getSubplots(self):
475 475
476 476 ncol = 1
477 477 nrow = self.nplots
478 478
479 479 return nrow, ncol
480 480
481 481 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
482 482
483 483 self.__showprofile = showprofile
484 484 self.nplots = nplots
485 485
486 486 ncolspan = 1
487 487 colspan = 1
488 488 if showprofile:
489 489 ncolspan = 7
490 490 colspan = 6
491 491 self.__nsubplots = 2
492 492
493 493 self.createFigure(id = id,
494 494 wintitle = wintitle,
495 495 widthplot = self.WIDTH + self.WIDTHPROF,
496 496 heightplot = self.HEIGHT + self.HEIGHTPROF,
497 497 show=show)
498 498
499 499 nrow, ncol = self.getSubplots()
500 500
501 501 counter = 0
502 502 for y in range(nrow):
503 503 for x in range(ncol):
504 504
505 505 if counter >= self.nplots:
506 506 break
507 507
508 508 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
509 509
510 510 if showprofile:
511 511 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
512 512
513 513 counter += 1
514 514
515 515 def run(self, dataOut, id, wintitle="", channelList=None, showprofile='True',
516 516 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
517 517 timerange=None,
518 518 save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
519 519 server=None, folder=None, username=None, password=None,
520 520 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
521 521
522 522 """
523 523
524 524 Input:
525 525 dataOut :
526 526 id :
527 527 wintitle :
528 528 channelList :
529 529 showProfile :
530 530 xmin : None,
531 531 xmax : None,
532 532 ymin : None,
533 533 ymax : None,
534 534 zmin : None,
535 535 zmax : None
536 536 """
537 537
538 538 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
539 539 return
540 540
541 541 if channelList == None:
542 542 channelIndexList = dataOut.channelIndexList
543 543 else:
544 544 channelIndexList = []
545 545 for channel in channelList:
546 546 if channel not in dataOut.channelList:
547 547 raise ValueError, "Channel %d is not in dataOut.channelList"
548 548 channelIndexList.append(dataOut.channelList.index(channel))
549 549
550 if hasattr(dataOut, 'normFactor'):
550 551 factor = dataOut.normFactor
552 else:
553 factor = 1
554
555 # factor = dataOut.normFactor
551 556 x = dataOut.getTimeRange()
552 557 y = dataOut.getHeiRange()
553 558
554 z = dataOut.data_spc/factor
555 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
556 avg = numpy.average(z, axis=1)
557
558 avgdB = 10.*numpy.log10(avg)
559 # z = dataOut.data_spc/factor
560 # z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
561 # avg = numpy.average(z, axis=1)
562 # avgdB = 10.*numpy.log10(avg)
563 avgdB = dataOut.getPower()
559 564
560 565 thisDatetime = dataOut.datatime
561 566 # thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
562 567 title = wintitle + " RTI" #: %s" %(thisDatetime.strftime("%d-%b-%Y"))
563 568 xlabel = ""
564 569 ylabel = "Range (Km)"
565 570
566 571 update_figfile = False
567 572
568 573 if dataOut.ltctime >= self.xmax:
569 574 self.counter_imagwr = wr_period
570 575 self.isConfig = False
571 576 update_figfile = True
572 577
573 578 if not self.isConfig:
574 579
575 580 nplots = len(channelIndexList)
576 581
577 582 self.setup(id=id,
578 583 nplots=nplots,
579 584 wintitle=wintitle,
580 585 showprofile=showprofile,
581 586 show=show)
582 587
583 588 if timerange != None:
584 589 self.timerange = timerange
585 590
586 591 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
587 592
588 593 noise = dataOut.noise/factor
589 594 noisedB = 10*numpy.log10(noise)
590 595
591 596 if ymin == None: ymin = numpy.nanmin(y)
592 597 if ymax == None: ymax = numpy.nanmax(y)
593 598 if zmin == None: zmin = numpy.floor(numpy.nanmin(noisedB)) - 3
594 599 if zmax == None: zmax = numpy.ceil(numpy.nanmax(avgdB)) + 3
595 600
596 601 self.FTP_WEI = ftp_wei
597 602 self.EXP_CODE = exp_code
598 603 self.SUB_EXP_CODE = sub_exp_code
599 604 self.PLOT_POS = plot_pos
600 605
601 606 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
602 607 self.isConfig = True
603 608 self.figfile = figfile
604 609 update_figfile = True
605 610
606 611 self.setWinTitle(title)
607 612
608 613 for i in range(self.nplots):
609 614 index = channelIndexList[i]
610 615 title = "Channel %d: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
611 616 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
612 617 title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith)
613 618 axes = self.axesList[i*self.__nsubplots]
614 619 zdB = avgdB[index].reshape((1,-1))
615 620 axes.pcolorbuffer(x, y, zdB,
616 621 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
617 622 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
618 623 ticksize=9, cblabel='', cbsize="1%")
619 624
620 625 if self.__showprofile:
621 626 axes = self.axesList[i*self.__nsubplots +1]
622 627 axes.pline(avgdB[index], y,
623 628 xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
624 629 xlabel='dB', ylabel='', title='',
625 630 ytick_visible=False,
626 631 grid='x')
627 632
628 633 self.draw()
629 634
630 635 self.save(figpath=figpath,
631 636 figfile=figfile,
632 637 save=save,
633 638 ftp=ftp,
634 639 wr_period=wr_period,
635 640 thisDatetime=thisDatetime,
636 641 update_figfile=update_figfile)
637 642
638 643 class CoherenceMap(Figure):
639 644 isConfig = None
640 645 __nsubplots = None
641 646
642 647 WIDTHPROF = None
643 648 HEIGHTPROF = None
644 649 PREFIX = 'cmap'
645 650
646 651 def __init__(self):
647 652 self.timerange = 2*60*60
648 653 self.isConfig = False
649 654 self.__nsubplots = 1
650 655
651 656 self.WIDTH = 800
652 657 self.HEIGHT = 180
653 658 self.WIDTHPROF = 120
654 659 self.HEIGHTPROF = 0
655 660 self.counter_imagwr = 0
656 661
657 662 self.PLOT_CODE = COH_CODE
658 663
659 664 self.FTP_WEI = None
660 665 self.EXP_CODE = None
661 666 self.SUB_EXP_CODE = None
662 667 self.PLOT_POS = None
663 668 self.counter_imagwr = 0
664 669
665 670 self.xmin = None
666 671 self.xmax = None
667 672
668 673 def getSubplots(self):
669 674 ncol = 1
670 675 nrow = self.nplots*2
671 676
672 677 return nrow, ncol
673 678
674 679 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
675 680 self.__showprofile = showprofile
676 681 self.nplots = nplots
677 682
678 683 ncolspan = 1
679 684 colspan = 1
680 685 if showprofile:
681 686 ncolspan = 7
682 687 colspan = 6
683 688 self.__nsubplots = 2
684 689
685 690 self.createFigure(id = id,
686 691 wintitle = wintitle,
687 692 widthplot = self.WIDTH + self.WIDTHPROF,
688 693 heightplot = self.HEIGHT + self.HEIGHTPROF,
689 694 show=True)
690 695
691 696 nrow, ncol = self.getSubplots()
692 697
693 698 for y in range(nrow):
694 699 for x in range(ncol):
695 700
696 701 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
697 702
698 703 if showprofile:
699 704 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
700 705
701 706 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
702 707 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
703 708 timerange=None, phase_min=None, phase_max=None,
704 709 save=False, figpath='./', figfile=None, ftp=False, wr_period=1,
705 710 coherence_cmap='jet', phase_cmap='RdBu_r', show=True,
706 711 server=None, folder=None, username=None, password=None,
707 712 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
708 713
709 714 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
710 715 return
711 716
712 717 if pairsList == None:
713 718 pairsIndexList = dataOut.pairsIndexList
714 719 else:
715 720 pairsIndexList = []
716 721 for pair in pairsList:
717 722 if pair not in dataOut.pairsList:
718 723 raise ValueError, "Pair %s is not in dataOut.pairsList" %(pair)
719 724 pairsIndexList.append(dataOut.pairsList.index(pair))
720 725
721 726 if pairsIndexList == []:
722 727 return
723 728
724 729 if len(pairsIndexList) > 4:
725 730 pairsIndexList = pairsIndexList[0:4]
726 731
727 732 if phase_min == None:
728 733 phase_min = -180
729 734 if phase_max == None:
730 735 phase_max = 180
731 736
732 737 x = dataOut.getTimeRange()
733 738 y = dataOut.getHeiRange()
734 739
735 740 thisDatetime = dataOut.datatime
736 741
737 742 title = wintitle + " CoherenceMap" #: %s" %(thisDatetime.strftime("%d-%b-%Y"))
738 743 xlabel = ""
739 744 ylabel = "Range (Km)"
740 745 update_figfile = False
741 746
742 747 if not self.isConfig:
743 748 nplots = len(pairsIndexList)
744 749 self.setup(id=id,
745 750 nplots=nplots,
746 751 wintitle=wintitle,
747 752 showprofile=showprofile,
748 753 show=show)
749 754
750 755 if timerange != None:
751 756 self.timerange = timerange
752 757
753 758 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
754 759
755 760 if ymin == None: ymin = numpy.nanmin(y)
756 761 if ymax == None: ymax = numpy.nanmax(y)
757 762 if zmin == None: zmin = 0.
758 763 if zmax == None: zmax = 1.
759 764
760 765 self.FTP_WEI = ftp_wei
761 766 self.EXP_CODE = exp_code
762 767 self.SUB_EXP_CODE = sub_exp_code
763 768 self.PLOT_POS = plot_pos
764 769
765 770 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
766 771
767 772 self.isConfig = True
768 773 update_figfile = True
769 774
770 775 self.setWinTitle(title)
771 776
772 777 for i in range(self.nplots):
773 778
774 779 pair = dataOut.pairsList[pairsIndexList[i]]
775 780
776 781 ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i],:,:],axis=0)
777 782 powa = numpy.average(dataOut.data_spc[pair[0],:,:],axis=0)
778 783 powb = numpy.average(dataOut.data_spc[pair[1],:,:],axis=0)
779 784
780 785
781 786 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
782 787 coherence = numpy.abs(avgcoherenceComplex)
783 788
784 789 z = coherence.reshape((1,-1))
785 790
786 791 counter = 0
787 792
788 793 title = "Coherence Ch%d * Ch%d: %s" %(pair[0], pair[1], thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
789 794 axes = self.axesList[i*self.__nsubplots*2]
790 795 axes.pcolorbuffer(x, y, z,
791 796 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
792 797 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
793 798 ticksize=9, cblabel='', colormap=coherence_cmap, cbsize="1%")
794 799
795 800 if self.__showprofile:
796 801 counter += 1
797 802 axes = self.axesList[i*self.__nsubplots*2 + counter]
798 803 axes.pline(coherence, y,
799 804 xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
800 805 xlabel='', ylabel='', title='', ticksize=7,
801 806 ytick_visible=False, nxticks=5,
802 807 grid='x')
803 808
804 809 counter += 1
805 810
806 811 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
807 812
808 813 z = phase.reshape((1,-1))
809 814
810 815 title = "Phase Ch%d * Ch%d: %s" %(pair[0], pair[1], thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
811 816 axes = self.axesList[i*self.__nsubplots*2 + counter]
812 817 axes.pcolorbuffer(x, y, z,
813 818 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=phase_min, zmax=phase_max,
814 819 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
815 820 ticksize=9, cblabel='', colormap=phase_cmap, cbsize="1%")
816 821
817 822 if self.__showprofile:
818 823 counter += 1
819 824 axes = self.axesList[i*self.__nsubplots*2 + counter]
820 825 axes.pline(phase, y,
821 826 xmin=phase_min, xmax=phase_max, ymin=ymin, ymax=ymax,
822 827 xlabel='', ylabel='', title='', ticksize=7,
823 828 ytick_visible=False, nxticks=4,
824 829 grid='x')
825 830
826 831 self.draw()
827 832
828 833 if dataOut.ltctime >= self.xmax:
829 834 self.counter_imagwr = wr_period
830 835 self.isConfig = False
831 836 update_figfile = True
832 837
833 838 self.save(figpath=figpath,
834 839 figfile=figfile,
835 840 save=save,
836 841 ftp=ftp,
837 842 wr_period=wr_period,
838 843 thisDatetime=thisDatetime,
839 844 update_figfile=update_figfile)
840 845
841 846 class PowerProfilePlot(Figure):
842 847
843 848 isConfig = None
844 849 __nsubplots = None
845 850
846 851 WIDTHPROF = None
847 852 HEIGHTPROF = None
848 853 PREFIX = 'spcprofile'
849 854
850 855 def __init__(self):
851 856 self.isConfig = False
852 857 self.__nsubplots = 1
853 858
854 859 self.PLOT_CODE = POWER_CODE
855 860
856 861 self.WIDTH = 300
857 862 self.HEIGHT = 500
858 863 self.counter_imagwr = 0
859 864
860 865 def getSubplots(self):
861 866 ncol = 1
862 867 nrow = 1
863 868
864 869 return nrow, ncol
865 870
866 871 def setup(self, id, nplots, wintitle, show):
867 872
868 873 self.nplots = nplots
869 874
870 875 ncolspan = 1
871 876 colspan = 1
872 877
873 878 self.createFigure(id = id,
874 879 wintitle = wintitle,
875 880 widthplot = self.WIDTH,
876 881 heightplot = self.HEIGHT,
877 882 show=show)
878 883
879 884 nrow, ncol = self.getSubplots()
880 885
881 886 counter = 0
882 887 for y in range(nrow):
883 888 for x in range(ncol):
884 889 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
885 890
886 891 def run(self, dataOut, id, wintitle="", channelList=None,
887 892 xmin=None, xmax=None, ymin=None, ymax=None,
888 893 save=False, figpath='./', figfile=None, show=True,
889 894 ftp=False, wr_period=1, server=None,
890 895 folder=None, username=None, password=None):
891 896
892 897
893 898 if channelList == None:
894 899 channelIndexList = dataOut.channelIndexList
895 900 channelList = dataOut.channelList
896 901 else:
897 902 channelIndexList = []
898 903 for channel in channelList:
899 904 if channel not in dataOut.channelList:
900 905 raise ValueError, "Channel %d is not in dataOut.channelList"
901 906 channelIndexList.append(dataOut.channelList.index(channel))
902 907
903 908 factor = dataOut.normFactor
904 909
905 910 y = dataOut.getHeiRange()
906 911
907 912 #for voltage
908 913 if dataOut.type == 'Voltage':
909 914 x = dataOut.data[channelIndexList,:] * numpy.conjugate(dataOut.data[channelIndexList,:])
910 915 x = x.real
911 916 x = numpy.where(numpy.isfinite(x), x, numpy.NAN)
912 917
913 918 #for spectra
914 919 if dataOut.type == 'Spectra':
915 920 x = dataOut.data_spc[channelIndexList,:,:]/factor
916 921 x = numpy.where(numpy.isfinite(x), x, numpy.NAN)
917 922 x = numpy.average(x, axis=1)
918 923
919 924
920 925 xdB = 10*numpy.log10(x)
921 926
922 927 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
923 928 title = wintitle + " Power Profile %s" %(thisDatetime.strftime("%d-%b-%Y"))
924 929 xlabel = "dB"
925 930 ylabel = "Range (Km)"
926 931
927 932 if not self.isConfig:
928 933
929 934 nplots = 1
930 935
931 936 self.setup(id=id,
932 937 nplots=nplots,
933 938 wintitle=wintitle,
934 939 show=show)
935 940
936 941 if ymin == None: ymin = numpy.nanmin(y)
937 942 if ymax == None: ymax = numpy.nanmax(y)
938 943 if xmin == None: xmin = numpy.nanmin(xdB)*0.9
939 944 if xmax == None: xmax = numpy.nanmax(xdB)*1.1
940 945
941 946 self.isConfig = True
942 947
943 948 self.setWinTitle(title)
944 949
945 950 title = "Power Profile: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
946 951 axes = self.axesList[0]
947 952
948 953 legendlabels = ["channel %d"%x for x in channelList]
949 954 axes.pmultiline(xdB, y,
950 955 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
951 956 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels,
952 957 ytick_visible=True, nxticks=5,
953 958 grid='x')
954 959
955 960 self.draw()
956 961
957 962 self.save(figpath=figpath,
958 963 figfile=figfile,
959 964 save=save,
960 965 ftp=ftp,
961 966 wr_period=wr_period,
962 967 thisDatetime=thisDatetime)
963 968
964 969 class SpectraCutPlot(Figure):
965 970
966 971 isConfig = None
967 972 __nsubplots = None
968 973
969 974 WIDTHPROF = None
970 975 HEIGHTPROF = None
971 976 PREFIX = 'spc_cut'
972 977
973 978 def __init__(self):
974 979 self.isConfig = False
975 980 self.__nsubplots = 1
976 981
977 982 self.PLOT_CODE = POWER_CODE
978 983
979 984 self.WIDTH = 700
980 985 self.HEIGHT = 500
981 986 self.counter_imagwr = 0
982 987
983 988 def getSubplots(self):
984 989 ncol = 1
985 990 nrow = 1
986 991
987 992 return nrow, ncol
988 993
989 994 def setup(self, id, nplots, wintitle, show):
990 995
991 996 self.nplots = nplots
992 997
993 998 ncolspan = 1
994 999 colspan = 1
995 1000
996 1001 self.createFigure(id = id,
997 1002 wintitle = wintitle,
998 1003 widthplot = self.WIDTH,
999 1004 heightplot = self.HEIGHT,
1000 1005 show=show)
1001 1006
1002 1007 nrow, ncol = self.getSubplots()
1003 1008
1004 1009 counter = 0
1005 1010 for y in range(nrow):
1006 1011 for x in range(ncol):
1007 1012 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
1008 1013
1009 1014 def run(self, dataOut, id, wintitle="", channelList=None,
1010 1015 xmin=None, xmax=None, ymin=None, ymax=None,
1011 1016 save=False, figpath='./', figfile=None, show=True,
1012 1017 ftp=False, wr_period=1, server=None,
1013 1018 folder=None, username=None, password=None,
1014 1019 xaxis="frequency"):
1015 1020
1016 1021
1017 1022 if channelList == None:
1018 1023 channelIndexList = dataOut.channelIndexList
1019 1024 channelList = dataOut.channelList
1020 1025 else:
1021 1026 channelIndexList = []
1022 1027 for channel in channelList:
1023 1028 if channel not in dataOut.channelList:
1024 1029 raise ValueError, "Channel %d is not in dataOut.channelList"
1025 1030 channelIndexList.append(dataOut.channelList.index(channel))
1026 1031
1027 1032 factor = dataOut.normFactor
1028 1033
1029 1034 y = dataOut.getHeiRange()
1030 1035
1031 1036 z = dataOut.data_spc/factor
1032 1037 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
1033 1038
1034 1039 hei_index = numpy.arange(25)*3 + 20
1035 1040
1036 1041 if xaxis == "frequency":
1037 1042 x = dataOut.getFreqRange()/1000.
1038 1043 zdB = 10*numpy.log10(z[0,:,hei_index])
1039 1044 xlabel = "Frequency (kHz)"
1040 1045 ylabel = "Power (dB)"
1041 1046
1042 1047 elif xaxis == "time":
1043 1048 x = dataOut.getAcfRange()
1044 1049 zdB = z[0,:,hei_index]
1045 1050 xlabel = "Time (ms)"
1046 1051 ylabel = "ACF"
1047 1052
1048 1053 else:
1049 1054 x = dataOut.getVelRange()
1050 1055 zdB = 10*numpy.log10(z[0,:,hei_index])
1051 1056 xlabel = "Velocity (m/s)"
1052 1057 ylabel = "Power (dB)"
1053 1058
1054 1059 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
1055 1060 title = wintitle + " Range Cuts %s" %(thisDatetime.strftime("%d-%b-%Y"))
1056 1061
1057 1062 if not self.isConfig:
1058 1063
1059 1064 nplots = 1
1060 1065
1061 1066 self.setup(id=id,
1062 1067 nplots=nplots,
1063 1068 wintitle=wintitle,
1064 1069 show=show)
1065 1070
1066 1071 if xmin == None: xmin = numpy.nanmin(x)*0.9
1067 1072 if xmax == None: xmax = numpy.nanmax(x)*1.1
1068 1073 if ymin == None: ymin = numpy.nanmin(zdB)
1069 1074 if ymax == None: ymax = numpy.nanmax(zdB)
1070 1075
1071 1076 self.isConfig = True
1072 1077
1073 1078 self.setWinTitle(title)
1074 1079
1075 1080 title = "Spectra Cuts: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
1076 1081 axes = self.axesList[0]
1077 1082
1078 1083 legendlabels = ["Range = %dKm" %y[i] for i in hei_index]
1079 1084
1080 1085 axes.pmultilineyaxis( x, zdB,
1081 1086 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
1082 1087 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels,
1083 1088 ytick_visible=True, nxticks=5,
1084 1089 grid='x')
1085 1090
1086 1091 self.draw()
1087 1092
1088 1093 self.save(figpath=figpath,
1089 1094 figfile=figfile,
1090 1095 save=save,
1091 1096 ftp=ftp,
1092 1097 wr_period=wr_period,
1093 1098 thisDatetime=thisDatetime)
1094 1099
1095 1100 class Noise(Figure):
1096 1101
1097 1102 isConfig = None
1098 1103 __nsubplots = None
1099 1104
1100 1105 PREFIX = 'noise'
1101 1106
1102 1107 def __init__(self):
1103 1108
1104 1109 self.timerange = 24*60*60
1105 1110 self.isConfig = False
1106 1111 self.__nsubplots = 1
1107 1112 self.counter_imagwr = 0
1108 1113 self.WIDTH = 800
1109 1114 self.HEIGHT = 400
1110 1115 self.WIDTHPROF = 120
1111 1116 self.HEIGHTPROF = 0
1112 1117 self.xdata = None
1113 1118 self.ydata = None
1114 1119
1115 1120 self.PLOT_CODE = NOISE_CODE
1116 1121
1117 1122 self.FTP_WEI = None
1118 1123 self.EXP_CODE = None
1119 1124 self.SUB_EXP_CODE = None
1120 1125 self.PLOT_POS = None
1121 1126 self.figfile = None
1122 1127
1123 1128 self.xmin = None
1124 1129 self.xmax = None
1125 1130
1126 1131 def getSubplots(self):
1127 1132
1128 1133 ncol = 1
1129 1134 nrow = 1
1130 1135
1131 1136 return nrow, ncol
1132 1137
1133 1138 def openfile(self, filename):
1134 1139 dirname = os.path.dirname(filename)
1135 1140
1136 1141 if not os.path.exists(dirname):
1137 1142 os.mkdir(dirname)
1138 1143
1139 1144 f = open(filename,'w+')
1140 1145 f.write('\n\n')
1141 1146 f.write('JICAMARCA RADIO OBSERVATORY - Noise \n')
1142 1147 f.write('DD MM YYYY HH MM SS Channel0 Channel1 Channel2 Channel3\n\n' )
1143 1148 f.close()
1144 1149
1145 1150 def save_data(self, filename_phase, data, data_datetime):
1146 1151
1147 1152 f=open(filename_phase,'a')
1148 1153
1149 1154 timetuple_data = data_datetime.timetuple()
1150 1155 day = str(timetuple_data.tm_mday)
1151 1156 month = str(timetuple_data.tm_mon)
1152 1157 year = str(timetuple_data.tm_year)
1153 1158 hour = str(timetuple_data.tm_hour)
1154 1159 minute = str(timetuple_data.tm_min)
1155 1160 second = str(timetuple_data.tm_sec)
1156 1161
1157 1162 data_msg = ''
1158 1163 for i in range(len(data)):
1159 1164 data_msg += str(data[i]) + ' '
1160 1165
1161 1166 f.write(day+' '+month+' '+year+' '+hour+' '+minute+' '+second+' ' + data_msg + '\n')
1162 1167 f.close()
1163 1168
1164 1169
1165 1170 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
1166 1171
1167 1172 self.__showprofile = showprofile
1168 1173 self.nplots = nplots
1169 1174
1170 1175 ncolspan = 7
1171 1176 colspan = 6
1172 1177 self.__nsubplots = 2
1173 1178
1174 1179 self.createFigure(id = id,
1175 1180 wintitle = wintitle,
1176 1181 widthplot = self.WIDTH+self.WIDTHPROF,
1177 1182 heightplot = self.HEIGHT+self.HEIGHTPROF,
1178 1183 show=show)
1179 1184
1180 1185 nrow, ncol = self.getSubplots()
1181 1186
1182 1187 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
1183 1188
1184 1189
1185 1190 def run(self, dataOut, id, wintitle="", channelList=None, showprofile='True',
1186 1191 xmin=None, xmax=None, ymin=None, ymax=None,
1187 1192 timerange=None,
1188 1193 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
1189 1194 server=None, folder=None, username=None, password=None,
1190 1195 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1191 1196
1192 1197 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
1193 1198 return
1194 1199
1195 1200 if channelList == None:
1196 1201 channelIndexList = dataOut.channelIndexList
1197 1202 channelList = dataOut.channelList
1198 1203 else:
1199 1204 channelIndexList = []
1200 1205 for channel in channelList:
1201 1206 if channel not in dataOut.channelList:
1202 1207 raise ValueError, "Channel %d is not in dataOut.channelList"
1203 1208 channelIndexList.append(dataOut.channelList.index(channel))
1204 1209
1205 1210 x = dataOut.getTimeRange()
1206 1211 #y = dataOut.getHeiRange()
1207 1212 factor = dataOut.normFactor
1208 1213 noise = dataOut.noise[channelIndexList]/factor
1209 1214 noisedB = 10*numpy.log10(noise)
1210 1215
1211 1216 thisDatetime = dataOut.datatime
1212 1217
1213 1218 title = wintitle + " Noise" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
1214 1219 xlabel = ""
1215 1220 ylabel = "Intensity (dB)"
1216 1221 update_figfile = False
1217 1222
1218 1223 if not self.isConfig:
1219 1224
1220 1225 nplots = 1
1221 1226
1222 1227 self.setup(id=id,
1223 1228 nplots=nplots,
1224 1229 wintitle=wintitle,
1225 1230 showprofile=showprofile,
1226 1231 show=show)
1227 1232
1228 1233 if timerange != None:
1229 1234 self.timerange = timerange
1230 1235
1231 1236 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
1232 1237
1233 1238 if ymin == None: ymin = numpy.floor(numpy.nanmin(noisedB)) - 10.0
1234 1239 if ymax == None: ymax = numpy.nanmax(noisedB) + 10.0
1235 1240
1236 1241 self.FTP_WEI = ftp_wei
1237 1242 self.EXP_CODE = exp_code
1238 1243 self.SUB_EXP_CODE = sub_exp_code
1239 1244 self.PLOT_POS = plot_pos
1240 1245
1241 1246
1242 1247 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1243 1248 self.isConfig = True
1244 1249 self.figfile = figfile
1245 1250 self.xdata = numpy.array([])
1246 1251 self.ydata = numpy.array([])
1247 1252
1248 1253 update_figfile = True
1249 1254
1250 1255 #open file beacon phase
1251 1256 path = '%s%03d' %(self.PREFIX, self.id)
1252 1257 noise_file = os.path.join(path,'%s.txt'%self.name)
1253 1258 self.filename_noise = os.path.join(figpath,noise_file)
1254 1259
1255 1260 self.setWinTitle(title)
1256 1261
1257 1262 title = "Noise %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1258 1263
1259 1264 legendlabels = ["channel %d"%(idchannel) for idchannel in channelList]
1260 1265 axes = self.axesList[0]
1261 1266
1262 1267 self.xdata = numpy.hstack((self.xdata, x[0:1]))
1263 1268
1264 1269 if len(self.ydata)==0:
1265 1270 self.ydata = noisedB.reshape(-1,1)
1266 1271 else:
1267 1272 self.ydata = numpy.hstack((self.ydata, noisedB.reshape(-1,1)))
1268 1273
1269 1274
1270 1275 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
1271 1276 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
1272 1277 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
1273 1278 XAxisAsTime=True, grid='both'
1274 1279 )
1275 1280
1276 1281 self.draw()
1277 1282
1278 1283 if dataOut.ltctime >= self.xmax:
1279 1284 self.counter_imagwr = wr_period
1280 1285 self.isConfig = False
1281 1286 update_figfile = True
1282 1287
1283 1288 self.save(figpath=figpath,
1284 1289 figfile=figfile,
1285 1290 save=save,
1286 1291 ftp=ftp,
1287 1292 wr_period=wr_period,
1288 1293 thisDatetime=thisDatetime,
1289 1294 update_figfile=update_figfile)
1290 1295
1291 1296 #store data beacon phase
1292 1297 if save:
1293 1298 self.save_data(self.filename_noise, noisedB, thisDatetime)
1294 1299
1295 1300 class BeaconPhase(Figure):
1296 1301
1297 1302 __isConfig = None
1298 1303 __nsubplots = None
1299 1304
1300 1305 PREFIX = 'beacon_phase'
1301 1306
1302 1307 def __init__(self):
1303 1308
1304 1309 self.timerange = 24*60*60
1305 1310 self.isConfig = False
1306 1311 self.__nsubplots = 1
1307 1312 self.counter_imagwr = 0
1308 1313 self.WIDTH = 800
1309 1314 self.HEIGHT = 400
1310 1315 self.WIDTHPROF = 120
1311 1316 self.HEIGHTPROF = 0
1312 1317 self.xdata = None
1313 1318 self.ydata = None
1314 1319
1315 1320 self.PLOT_CODE = BEACON_CODE
1316 1321
1317 1322 self.FTP_WEI = None
1318 1323 self.EXP_CODE = None
1319 1324 self.SUB_EXP_CODE = None
1320 1325 self.PLOT_POS = None
1321 1326
1322 1327 self.filename_phase = None
1323 1328
1324 1329 self.figfile = None
1325 1330
1326 1331 self.xmin = None
1327 1332 self.xmax = None
1328 1333
1329 1334 def getSubplots(self):
1330 1335
1331 1336 ncol = 1
1332 1337 nrow = 1
1333 1338
1334 1339 return nrow, ncol
1335 1340
1336 1341 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
1337 1342
1338 1343 self.__showprofile = showprofile
1339 1344 self.nplots = nplots
1340 1345
1341 1346 ncolspan = 7
1342 1347 colspan = 6
1343 1348 self.__nsubplots = 2
1344 1349
1345 1350 self.createFigure(id = id,
1346 1351 wintitle = wintitle,
1347 1352 widthplot = self.WIDTH+self.WIDTHPROF,
1348 1353 heightplot = self.HEIGHT+self.HEIGHTPROF,
1349 1354 show=show)
1350 1355
1351 1356 nrow, ncol = self.getSubplots()
1352 1357
1353 1358 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
1354 1359
1355 1360 def save_phase(self, filename_phase):
1356 1361 f = open(filename_phase,'w+')
1357 1362 f.write('\n\n')
1358 1363 f.write('JICAMARCA RADIO OBSERVATORY - Beacon Phase \n')
1359 1364 f.write('DD MM YYYY HH MM SS pair(2,0) pair(2,1) pair(2,3) pair(2,4)\n\n' )
1360 1365 f.close()
1361 1366
1362 1367 def save_data(self, filename_phase, data, data_datetime):
1363 1368 f=open(filename_phase,'a')
1364 1369 timetuple_data = data_datetime.timetuple()
1365 1370 day = str(timetuple_data.tm_mday)
1366 1371 month = str(timetuple_data.tm_mon)
1367 1372 year = str(timetuple_data.tm_year)
1368 1373 hour = str(timetuple_data.tm_hour)
1369 1374 minute = str(timetuple_data.tm_min)
1370 1375 second = str(timetuple_data.tm_sec)
1371 1376 f.write(day+' '+month+' '+year+' '+hour+' '+minute+' '+second+' '+str(data[0])+' '+str(data[1])+' '+str(data[2])+' '+str(data[3])+'\n')
1372 1377 f.close()
1373 1378
1374 1379
1375 1380 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
1376 1381 xmin=None, xmax=None, ymin=None, ymax=None, hmin=None, hmax=None,
1377 1382 timerange=None,
1378 1383 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
1379 1384 server=None, folder=None, username=None, password=None,
1380 1385 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1381 1386
1382 1387 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
1383 1388 return
1384 1389
1385 1390 if pairsList == None:
1386 1391 pairsIndexList = dataOut.pairsIndexList[:10]
1387 1392 else:
1388 1393 pairsIndexList = []
1389 1394 for pair in pairsList:
1390 1395 if pair not in dataOut.pairsList:
1391 1396 raise ValueError, "Pair %s is not in dataOut.pairsList" %(pair)
1392 1397 pairsIndexList.append(dataOut.pairsList.index(pair))
1393 1398
1394 1399 if pairsIndexList == []:
1395 1400 return
1396 1401
1397 1402 # if len(pairsIndexList) > 4:
1398 1403 # pairsIndexList = pairsIndexList[0:4]
1399 1404
1400 1405 hmin_index = None
1401 1406 hmax_index = None
1402 1407
1403 1408 if hmin != None and hmax != None:
1404 1409 indexes = numpy.arange(dataOut.nHeights)
1405 1410 hmin_list = indexes[dataOut.heightList >= hmin]
1406 1411 hmax_list = indexes[dataOut.heightList <= hmax]
1407 1412
1408 1413 if hmin_list.any():
1409 1414 hmin_index = hmin_list[0]
1410 1415
1411 1416 if hmax_list.any():
1412 1417 hmax_index = hmax_list[-1]+1
1413 1418
1414 1419 x = dataOut.getTimeRange()
1415 1420 #y = dataOut.getHeiRange()
1416 1421
1417 1422
1418 1423 thisDatetime = dataOut.datatime
1419 1424
1420 1425 title = wintitle + " Signal Phase" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
1421 1426 xlabel = "Local Time"
1422 1427 ylabel = "Phase (degrees)"
1423 1428
1424 1429 update_figfile = False
1425 1430
1426 1431 nplots = len(pairsIndexList)
1427 1432 #phase = numpy.zeros((len(pairsIndexList),len(dataOut.beacon_heiIndexList)))
1428 1433 phase_beacon = numpy.zeros(len(pairsIndexList))
1429 1434 for i in range(nplots):
1430 1435 pair = dataOut.pairsList[pairsIndexList[i]]
1431 1436 ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i], :, hmin_index:hmax_index], axis=0)
1432 1437 powa = numpy.average(dataOut.data_spc[pair[0], :, hmin_index:hmax_index], axis=0)
1433 1438 powb = numpy.average(dataOut.data_spc[pair[1], :, hmin_index:hmax_index], axis=0)
1434 1439 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
1435 1440 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
1436 1441
1437 1442 #print "Phase %d%d" %(pair[0], pair[1])
1438 1443 #print phase[dataOut.beacon_heiIndexList]
1439 1444
1440 1445 if dataOut.beacon_heiIndexList:
1441 1446 phase_beacon[i] = numpy.average(phase[dataOut.beacon_heiIndexList])
1442 1447 else:
1443 1448 phase_beacon[i] = numpy.average(phase)
1444 1449
1445 1450 if not self.isConfig:
1446 1451
1447 1452 nplots = len(pairsIndexList)
1448 1453
1449 1454 self.setup(id=id,
1450 1455 nplots=nplots,
1451 1456 wintitle=wintitle,
1452 1457 showprofile=showprofile,
1453 1458 show=show)
1454 1459
1455 1460 if timerange != None:
1456 1461 self.timerange = timerange
1457 1462
1458 1463 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
1459 1464
1460 1465 if ymin == None: ymin = 0
1461 1466 if ymax == None: ymax = 360
1462 1467
1463 1468 self.FTP_WEI = ftp_wei
1464 1469 self.EXP_CODE = exp_code
1465 1470 self.SUB_EXP_CODE = sub_exp_code
1466 1471 self.PLOT_POS = plot_pos
1467 1472
1468 1473 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1469 1474 self.isConfig = True
1470 1475 self.figfile = figfile
1471 1476 self.xdata = numpy.array([])
1472 1477 self.ydata = numpy.array([])
1473 1478
1474 1479 update_figfile = True
1475 1480
1476 1481 #open file beacon phase
1477 1482 path = '%s%03d' %(self.PREFIX, self.id)
1478 1483 beacon_file = os.path.join(path,'%s.txt'%self.name)
1479 1484 self.filename_phase = os.path.join(figpath,beacon_file)
1480 1485 #self.save_phase(self.filename_phase)
1481 1486
1482 1487
1483 1488 #store data beacon phase
1484 1489 #self.save_data(self.filename_phase, phase_beacon, thisDatetime)
1485 1490
1486 1491 self.setWinTitle(title)
1487 1492
1488 1493
1489 1494 title = "Phase Plot %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1490 1495
1491 1496 legendlabels = ["Pair (%d,%d)"%(pair[0], pair[1]) for pair in dataOut.pairsList]
1492 1497
1493 1498 axes = self.axesList[0]
1494 1499
1495 1500 self.xdata = numpy.hstack((self.xdata, x[0:1]))
1496 1501
1497 1502 if len(self.ydata)==0:
1498 1503 self.ydata = phase_beacon.reshape(-1,1)
1499 1504 else:
1500 1505 self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1)))
1501 1506
1502 1507
1503 1508 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
1504 1509 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
1505 1510 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
1506 1511 XAxisAsTime=True, grid='both'
1507 1512 )
1508 1513
1509 1514 self.draw()
1510 1515
1511 1516 if dataOut.ltctime >= self.xmax:
1512 1517 self.counter_imagwr = wr_period
1513 1518 self.isConfig = False
1514 1519 update_figfile = True
1515 1520
1516 1521 self.save(figpath=figpath,
1517 1522 figfile=figfile,
1518 1523 save=save,
1519 1524 ftp=ftp,
1520 1525 wr_period=wr_period,
1521 1526 thisDatetime=thisDatetime,
1522 1527 update_figfile=update_figfile)
General Comments 0
You need to be logged in to leave comments. Login now