##// END OF EJS Templates
Fix SpectralMoments Plot
Juan C. Espinoza -
r1207:15feaa6e0b57
parent child
Show More
@@ -1,1353 +1,1358
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 import json
11 11
12 12 from schainpy.utils import log
13 13 from .jroheaderIO import SystemHeader, RadarControllerHeader
14 14
15 15
16 16 def getNumpyDtype(dataTypeCode):
17 17
18 18 if dataTypeCode == 0:
19 19 numpyDtype = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
20 20 elif dataTypeCode == 1:
21 21 numpyDtype = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
22 22 elif dataTypeCode == 2:
23 23 numpyDtype = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
24 24 elif dataTypeCode == 3:
25 25 numpyDtype = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
26 26 elif dataTypeCode == 4:
27 27 numpyDtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
28 28 elif dataTypeCode == 5:
29 29 numpyDtype = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
30 30 else:
31 31 raise ValueError('dataTypeCode was not defined')
32 32
33 33 return numpyDtype
34 34
35 35
36 36 def getDataTypeCode(numpyDtype):
37 37
38 38 if numpyDtype == numpy.dtype([('real', '<i1'), ('imag', '<i1')]):
39 39 datatype = 0
40 40 elif numpyDtype == numpy.dtype([('real', '<i2'), ('imag', '<i2')]):
41 41 datatype = 1
42 42 elif numpyDtype == numpy.dtype([('real', '<i4'), ('imag', '<i4')]):
43 43 datatype = 2
44 44 elif numpyDtype == numpy.dtype([('real', '<i8'), ('imag', '<i8')]):
45 45 datatype = 3
46 46 elif numpyDtype == numpy.dtype([('real', '<f4'), ('imag', '<f4')]):
47 47 datatype = 4
48 48 elif numpyDtype == numpy.dtype([('real', '<f8'), ('imag', '<f8')]):
49 49 datatype = 5
50 50 else:
51 51 datatype = None
52 52
53 53 return datatype
54 54
55 55
56 56 def hildebrand_sekhon(data, navg):
57 57 """
58 58 This method is for the objective determination of the noise level in Doppler spectra. This
59 59 implementation technique is based on the fact that the standard deviation of the spectral
60 60 densities is equal to the mean spectral density for white Gaussian noise
61 61
62 62 Inputs:
63 63 Data : heights
64 64 navg : numbers of averages
65 65
66 66 Return:
67 67 mean : noise's level
68 68 """
69 69
70 70 sortdata = numpy.sort(data, axis=None)
71 71 lenOfData = len(sortdata)
72 72 nums_min = lenOfData*0.2
73 73
74 74 if nums_min <= 5:
75 75
76 76 nums_min = 5
77 77
78 78 sump = 0.
79 79 sumq = 0.
80 80
81 81 j = 0
82 82 cont = 1
83 83
84 84 while((cont == 1)and(j < lenOfData)):
85 85
86 86 sump += sortdata[j]
87 87 sumq += sortdata[j]**2
88 88
89 89 if j > nums_min:
90 90 rtest = float(j)/(j-1) + 1.0/navg
91 91 if ((sumq*j) > (rtest*sump**2)):
92 92 j = j - 1
93 93 sump = sump - sortdata[j]
94 94 sumq = sumq - sortdata[j]**2
95 95 cont = 0
96 96
97 97 j += 1
98 98
99 99 lnoise = sump / j
100 100
101 101 return lnoise
102 102
103 103
104 104 class Beam:
105 105
106 106 def __init__(self):
107 107 self.codeList = []
108 108 self.azimuthList = []
109 109 self.zenithList = []
110 110
111 111
112 112 class GenericData(object):
113 113
114 114 flagNoData = True
115 115
116 116 def copy(self, inputObj=None):
117 117
118 118 if inputObj == None:
119 119 return copy.deepcopy(self)
120 120
121 121 for key in list(inputObj.__dict__.keys()):
122 122
123 123 attribute = inputObj.__dict__[key]
124 124
125 125 # If this attribute is a tuple or list
126 126 if type(inputObj.__dict__[key]) in (tuple, list):
127 127 self.__dict__[key] = attribute[:]
128 128 continue
129 129
130 130 # If this attribute is another object or instance
131 131 if hasattr(attribute, '__dict__'):
132 132 self.__dict__[key] = attribute.copy()
133 133 continue
134 134
135 135 self.__dict__[key] = inputObj.__dict__[key]
136 136
137 137 def deepcopy(self):
138 138
139 139 return copy.deepcopy(self)
140 140
141 141 def isEmpty(self):
142 142
143 143 return self.flagNoData
144 144
145 145
146 146 class JROData(GenericData):
147 147
148 148 # m_BasicHeader = BasicHeader()
149 149 # m_ProcessingHeader = ProcessingHeader()
150 150
151 151 systemHeaderObj = SystemHeader()
152 152 radarControllerHeaderObj = RadarControllerHeader()
153 153 # data = None
154 154 type = None
155 155 datatype = None # dtype but in string
156 156 # dtype = None
157 157 # nChannels = None
158 158 # nHeights = None
159 159 nProfiles = None
160 160 heightList = None
161 161 channelList = None
162 162 flagDiscontinuousBlock = False
163 163 useLocalTime = False
164 164 utctime = None
165 165 timeZone = None
166 166 dstFlag = None
167 167 errorCount = None
168 168 blocksize = None
169 169 # nCode = None
170 170 # nBaud = None
171 171 # code = None
172 172 flagDecodeData = False # asumo q la data no esta decodificada
173 173 flagDeflipData = False # asumo q la data no esta sin flip
174 174 flagShiftFFT = False
175 175 # ippSeconds = None
176 176 # timeInterval = None
177 177 nCohInt = None
178 178 # noise = None
179 179 windowOfFilter = 1
180 180 # Speed of ligth
181 181 C = 3e8
182 182 frequency = 49.92e6
183 183 realtime = False
184 184 beacon_heiIndexList = None
185 185 last_block = None
186 186 blocknow = None
187 187 azimuth = None
188 188 zenith = None
189 189 beam = Beam()
190 190 profileIndex = None
191 191 error = None
192 192 data = None
193 193 nmodes = None
194 194
195 195 def __str__(self):
196 196
197 197 return '{} - {}'.format(self.type, self.getDatatime())
198 198
199 199 def getNoise(self):
200 200
201 201 raise NotImplementedError
202 202
203 203 def getNChannels(self):
204 204
205 205 return len(self.channelList)
206 206
207 207 def getChannelIndexList(self):
208 208
209 209 return list(range(self.nChannels))
210 210
211 211 def getNHeights(self):
212 212
213 213 return len(self.heightList)
214 214
215 215 def getHeiRange(self, extrapoints=0):
216 216
217 217 heis = self.heightList
218 218 # deltah = self.heightList[1] - self.heightList[0]
219 219 #
220 220 # heis.append(self.heightList[-1])
221 221
222 222 return heis
223 223
224 224 def getDeltaH(self):
225 225
226 226 delta = self.heightList[1] - self.heightList[0]
227 227
228 228 return delta
229 229
230 230 def getltctime(self):
231 231
232 232 if self.useLocalTime:
233 233 return self.utctime - self.timeZone * 60
234 234
235 235 return self.utctime
236 236
237 237 def getDatatime(self):
238 238
239 239 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
240 240 return datatimeValue
241 241
242 242 def getTimeRange(self):
243 243
244 244 datatime = []
245 245
246 246 datatime.append(self.ltctime)
247 247 datatime.append(self.ltctime + self.timeInterval + 1)
248 248
249 249 datatime = numpy.array(datatime)
250 250
251 251 return datatime
252 252
253 253 def getFmaxTimeResponse(self):
254 254
255 255 period = (10**-6) * self.getDeltaH() / (0.15)
256 256
257 257 PRF = 1. / (period * self.nCohInt)
258 258
259 259 fmax = PRF
260 260
261 261 return fmax
262 262
263 263 def getFmax(self):
264 264 PRF = 1. / (self.ippSeconds * self.nCohInt)
265 265
266 266 fmax = PRF
267 267 return fmax
268 268
269 269 def getVmax(self):
270 270
271 271 _lambda = self.C / self.frequency
272 272
273 273 vmax = self.getFmax() * _lambda / 2
274 274
275 275 return vmax
276 276
277 277 def get_ippSeconds(self):
278 278 '''
279 279 '''
280 280 return self.radarControllerHeaderObj.ippSeconds
281 281
282 282 def set_ippSeconds(self, ippSeconds):
283 283 '''
284 284 '''
285 285
286 286 self.radarControllerHeaderObj.ippSeconds = ippSeconds
287 287
288 288 return
289 289
290 290 def get_dtype(self):
291 291 '''
292 292 '''
293 293 return getNumpyDtype(self.datatype)
294 294
295 295 def set_dtype(self, numpyDtype):
296 296 '''
297 297 '''
298 298
299 299 self.datatype = getDataTypeCode(numpyDtype)
300 300
301 301 def get_code(self):
302 302 '''
303 303 '''
304 304 return self.radarControllerHeaderObj.code
305 305
306 306 def set_code(self, code):
307 307 '''
308 308 '''
309 309 self.radarControllerHeaderObj.code = code
310 310
311 311 return
312 312
313 313 def get_ncode(self):
314 314 '''
315 315 '''
316 316 return self.radarControllerHeaderObj.nCode
317 317
318 318 def set_ncode(self, nCode):
319 319 '''
320 320 '''
321 321 self.radarControllerHeaderObj.nCode = nCode
322 322
323 323 return
324 324
325 325 def get_nbaud(self):
326 326 '''
327 327 '''
328 328 return self.radarControllerHeaderObj.nBaud
329 329
330 330 def set_nbaud(self, nBaud):
331 331 '''
332 332 '''
333 333 self.radarControllerHeaderObj.nBaud = nBaud
334 334
335 335 return
336 336
337 337 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
338 338 channelIndexList = property(
339 339 getChannelIndexList, "I'm the 'channelIndexList' property.")
340 340 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
341 341 #noise = property(getNoise, "I'm the 'nHeights' property.")
342 342 datatime = property(getDatatime, "I'm the 'datatime' property")
343 343 ltctime = property(getltctime, "I'm the 'ltctime' property")
344 344 ippSeconds = property(get_ippSeconds, set_ippSeconds)
345 345 dtype = property(get_dtype, set_dtype)
346 346 # timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
347 347 code = property(get_code, set_code)
348 348 nCode = property(get_ncode, set_ncode)
349 349 nBaud = property(get_nbaud, set_nbaud)
350 350
351 351
352 352 class Voltage(JROData):
353 353
354 354 # data es un numpy array de 2 dmensiones (canales, alturas)
355 355 data = None
356 356
357 357 def __init__(self):
358 358 '''
359 359 Constructor
360 360 '''
361 361
362 362 self.useLocalTime = True
363 363 self.radarControllerHeaderObj = RadarControllerHeader()
364 364 self.systemHeaderObj = SystemHeader()
365 365 self.type = "Voltage"
366 366 self.data = None
367 367 # self.dtype = None
368 368 # self.nChannels = 0
369 369 # self.nHeights = 0
370 370 self.nProfiles = None
371 371 self.heightList = None
372 372 self.channelList = None
373 373 # self.channelIndexList = None
374 374 self.flagNoData = True
375 375 self.flagDiscontinuousBlock = False
376 376 self.utctime = None
377 377 self.timeZone = None
378 378 self.dstFlag = None
379 379 self.errorCount = None
380 380 self.nCohInt = None
381 381 self.blocksize = None
382 382 self.flagDecodeData = False # asumo q la data no esta decodificada
383 383 self.flagDeflipData = False # asumo q la data no esta sin flip
384 384 self.flagShiftFFT = False
385 385 self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
386 386 self.profileIndex = 0
387 387
388 388 def getNoisebyHildebrand(self, channel=None):
389 389 """
390 390 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
391 391
392 392 Return:
393 393 noiselevel
394 394 """
395 395
396 396 if channel != None:
397 397 data = self.data[channel]
398 398 nChannels = 1
399 399 else:
400 400 data = self.data
401 401 nChannels = self.nChannels
402 402
403 403 noise = numpy.zeros(nChannels)
404 404 power = data * numpy.conjugate(data)
405 405
406 406 for thisChannel in range(nChannels):
407 407 if nChannels == 1:
408 408 daux = power[:].real
409 409 else:
410 410 daux = power[thisChannel, :].real
411 411 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
412 412
413 413 return noise
414 414
415 415 def getNoise(self, type=1, channel=None):
416 416
417 417 if type == 1:
418 418 noise = self.getNoisebyHildebrand(channel)
419 419
420 420 return noise
421 421
422 422 def getPower(self, channel=None):
423 423
424 424 if channel != None:
425 425 data = self.data[channel]
426 426 else:
427 427 data = self.data
428 428
429 429 power = data * numpy.conjugate(data)
430 430 powerdB = 10 * numpy.log10(power.real)
431 431 powerdB = numpy.squeeze(powerdB)
432 432
433 433 return powerdB
434 434
435 435 def getTimeInterval(self):
436 436
437 437 timeInterval = self.ippSeconds * self.nCohInt
438 438
439 439 return timeInterval
440 440
441 441 noise = property(getNoise, "I'm the 'nHeights' property.")
442 442 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
443 443
444 444
445 445 class Spectra(JROData):
446 446
447 447 # data spc es un numpy array de 2 dmensiones (canales, perfiles, alturas)
448 448 data_spc = None
449 449 # data cspc es un numpy array de 2 dmensiones (canales, pares, alturas)
450 450 data_cspc = None
451 451 # data dc es un numpy array de 2 dmensiones (canales, alturas)
452 452 data_dc = None
453 453 # data power
454 454 data_pwr = None
455 455 nFFTPoints = None
456 456 # nPairs = None
457 457 pairsList = None
458 458 nIncohInt = None
459 459 wavelength = None # Necesario para cacular el rango de velocidad desde la frecuencia
460 460 nCohInt = None # se requiere para determinar el valor de timeInterval
461 461 ippFactor = None
462 462 profileIndex = 0
463 463 plotting = "spectra"
464 464
465 465 def __init__(self):
466 466 '''
467 467 Constructor
468 468 '''
469 469
470 470 self.useLocalTime = True
471 471 self.radarControllerHeaderObj = RadarControllerHeader()
472 472 self.systemHeaderObj = SystemHeader()
473 473 self.type = "Spectra"
474 474 # self.data = None
475 475 # self.dtype = None
476 476 # self.nChannels = 0
477 477 # self.nHeights = 0
478 478 self.nProfiles = None
479 479 self.heightList = None
480 480 self.channelList = None
481 481 # self.channelIndexList = None
482 482 self.pairsList = None
483 483 self.flagNoData = True
484 484 self.flagDiscontinuousBlock = False
485 485 self.utctime = None
486 486 self.nCohInt = None
487 487 self.nIncohInt = None
488 488 self.blocksize = None
489 489 self.nFFTPoints = None
490 490 self.wavelength = None
491 491 self.flagDecodeData = False # asumo q la data no esta decodificada
492 492 self.flagDeflipData = False # asumo q la data no esta sin flip
493 493 self.flagShiftFFT = False
494 494 self.ippFactor = 1
495 495 #self.noise = None
496 496 self.beacon_heiIndexList = []
497 497 self.noise_estimation = None
498 498
499 499 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
500 500 """
501 501 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
502 502
503 503 Return:
504 504 noiselevel
505 505 """
506 506
507 507 noise = numpy.zeros(self.nChannels)
508 508
509 509 for channel in range(self.nChannels):
510 510 daux = self.data_spc[channel,
511 511 xmin_index:xmax_index, ymin_index:ymax_index]
512 512 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
513 513
514 514 return noise
515 515
516 516 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
517 517
518 518 if self.noise_estimation is not None:
519 519 # this was estimated by getNoise Operation defined in jroproc_spectra.py
520 520 return self.noise_estimation
521 521 else:
522 522 noise = self.getNoisebyHildebrand(
523 523 xmin_index, xmax_index, ymin_index, ymax_index)
524 524 return noise
525 525
526 526 def getFreqRangeTimeResponse(self, extrapoints=0):
527 527
528 528 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints * self.ippFactor)
529 529 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
530 530
531 531 return freqrange
532 532
533 533 def getAcfRange(self, extrapoints=0):
534 534
535 535 deltafreq = 10. / (self.getFmax() / (self.nFFTPoints * self.ippFactor))
536 536 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
537 537
538 538 return freqrange
539 539
540 540 def getFreqRange(self, extrapoints=0):
541 541
542 542 deltafreq = self.getFmax() / (self.nFFTPoints * self.ippFactor)
543 543 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
544 544
545 545 return freqrange
546 546
547 547 def getVelRange(self, extrapoints=0):
548 548
549 549 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
550 550 velrange = deltav * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.)
551 551
552 552 if self.nmodes:
553 553 return velrange/self.nmodes
554 554 else:
555 555 return velrange
556 556
557 557 def getNPairs(self):
558 558
559 559 return len(self.pairsList)
560 560
561 561 def getPairsIndexList(self):
562 562
563 563 return list(range(self.nPairs))
564 564
565 565 def getNormFactor(self):
566 566
567 567 pwcode = 1
568 568
569 569 if self.flagDecodeData:
570 570 pwcode = numpy.sum(self.code[0]**2)
571 571 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
572 572 normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
573 573
574 574 return normFactor
575 575
576 576 def getFlagCspc(self):
577 577
578 578 if self.data_cspc is None:
579 579 return True
580 580
581 581 return False
582 582
583 583 def getFlagDc(self):
584 584
585 585 if self.data_dc is None:
586 586 return True
587 587
588 588 return False
589 589
590 590 def getTimeInterval(self):
591 591
592 592 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles * self.ippFactor
593 593
594 594 return timeInterval
595 595
596 596 def getPower(self):
597 597
598 598 factor = self.normFactor
599 599 z = self.data_spc / factor
600 600 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
601 601 avg = numpy.average(z, axis=1)
602 602
603 603 return 10 * numpy.log10(avg)
604 604
605 605 def getCoherence(self, pairsList=None, phase=False):
606 606
607 607 z = []
608 608 if pairsList is None:
609 609 pairsIndexList = self.pairsIndexList
610 610 else:
611 611 pairsIndexList = []
612 612 for pair in pairsList:
613 613 if pair not in self.pairsList:
614 614 raise ValueError("Pair %s is not in dataOut.pairsList" % (
615 615 pair))
616 616 pairsIndexList.append(self.pairsList.index(pair))
617 617 for i in range(len(pairsIndexList)):
618 618 pair = self.pairsList[pairsIndexList[i]]
619 619 ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
620 620 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
621 621 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
622 622 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
623 623 if phase:
624 624 data = numpy.arctan2(avgcoherenceComplex.imag,
625 625 avgcoherenceComplex.real) * 180 / numpy.pi
626 626 else:
627 627 data = numpy.abs(avgcoherenceComplex)
628 628
629 629 z.append(data)
630 630
631 631 return numpy.array(z)
632 632
633 633 def setValue(self, value):
634 634
635 635 print("This property should not be initialized")
636 636
637 637 return
638 638
639 639 nPairs = property(getNPairs, setValue, "I'm the 'nPairs' property.")
640 640 pairsIndexList = property(
641 641 getPairsIndexList, setValue, "I'm the 'pairsIndexList' property.")
642 642 normFactor = property(getNormFactor, setValue,
643 643 "I'm the 'getNormFactor' property.")
644 644 flag_cspc = property(getFlagCspc, setValue)
645 645 flag_dc = property(getFlagDc, setValue)
646 646 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
647 647 timeInterval = property(getTimeInterval, setValue,
648 648 "I'm the 'timeInterval' property")
649 649
650 650
651 651 class SpectraHeis(Spectra):
652 652
653 653 data_spc = None
654 654 data_cspc = None
655 655 data_dc = None
656 656 nFFTPoints = None
657 657 # nPairs = None
658 658 pairsList = None
659 659 nCohInt = None
660 660 nIncohInt = None
661 661
662 662 def __init__(self):
663 663
664 664 self.radarControllerHeaderObj = RadarControllerHeader()
665 665
666 666 self.systemHeaderObj = SystemHeader()
667 667
668 668 self.type = "SpectraHeis"
669 669
670 670 # self.dtype = None
671 671
672 672 # self.nChannels = 0
673 673
674 674 # self.nHeights = 0
675 675
676 676 self.nProfiles = None
677 677
678 678 self.heightList = None
679 679
680 680 self.channelList = None
681 681
682 682 # self.channelIndexList = None
683 683
684 684 self.flagNoData = True
685 685
686 686 self.flagDiscontinuousBlock = False
687 687
688 688 # self.nPairs = 0
689 689
690 690 self.utctime = None
691 691
692 692 self.blocksize = None
693 693
694 694 self.profileIndex = 0
695 695
696 696 self.nCohInt = 1
697 697
698 698 self.nIncohInt = 1
699 699
700 700 def getNormFactor(self):
701 701 pwcode = 1
702 702 if self.flagDecodeData:
703 703 pwcode = numpy.sum(self.code[0]**2)
704 704
705 705 normFactor = self.nIncohInt * self.nCohInt * pwcode
706 706
707 707 return normFactor
708 708
709 709 def getTimeInterval(self):
710 710
711 711 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
712 712
713 713 return timeInterval
714 714
715 715 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
716 716 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
717 717
718 718
719 719 class Fits(JROData):
720 720
721 721 heightList = None
722 722 channelList = None
723 723 flagNoData = True
724 724 flagDiscontinuousBlock = False
725 725 useLocalTime = False
726 726 utctime = None
727 727 timeZone = None
728 728 # ippSeconds = None
729 729 # timeInterval = None
730 730 nCohInt = None
731 731 nIncohInt = None
732 732 noise = None
733 733 windowOfFilter = 1
734 734 # Speed of ligth
735 735 C = 3e8
736 736 frequency = 49.92e6
737 737 realtime = False
738 738
739 739 def __init__(self):
740 740
741 741 self.type = "Fits"
742 742
743 743 self.nProfiles = None
744 744
745 745 self.heightList = None
746 746
747 747 self.channelList = None
748 748
749 749 # self.channelIndexList = None
750 750
751 751 self.flagNoData = True
752 752
753 753 self.utctime = None
754 754
755 755 self.nCohInt = 1
756 756
757 757 self.nIncohInt = 1
758 758
759 759 self.useLocalTime = True
760 760
761 761 self.profileIndex = 0
762 762
763 763 # self.utctime = None
764 764 # self.timeZone = None
765 765 # self.ltctime = None
766 766 # self.timeInterval = None
767 767 # self.header = None
768 768 # self.data_header = None
769 769 # self.data = None
770 770 # self.datatime = None
771 771 # self.flagNoData = False
772 772 # self.expName = ''
773 773 # self.nChannels = None
774 774 # self.nSamples = None
775 775 # self.dataBlocksPerFile = None
776 776 # self.comments = ''
777 777 #
778 778
779 779 def getltctime(self):
780 780
781 781 if self.useLocalTime:
782 782 return self.utctime - self.timeZone * 60
783 783
784 784 return self.utctime
785 785
786 786 def getDatatime(self):
787 787
788 788 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
789 789 return datatime
790 790
791 791 def getTimeRange(self):
792 792
793 793 datatime = []
794 794
795 795 datatime.append(self.ltctime)
796 796 datatime.append(self.ltctime + self.timeInterval)
797 797
798 798 datatime = numpy.array(datatime)
799 799
800 800 return datatime
801 801
802 802 def getHeiRange(self):
803 803
804 804 heis = self.heightList
805 805
806 806 return heis
807 807
808 808 def getNHeights(self):
809 809
810 810 return len(self.heightList)
811 811
812 812 def getNChannels(self):
813 813
814 814 return len(self.channelList)
815 815
816 816 def getChannelIndexList(self):
817 817
818 818 return list(range(self.nChannels))
819 819
820 820 def getNoise(self, type=1):
821 821
822 822 #noise = numpy.zeros(self.nChannels)
823 823
824 824 if type == 1:
825 825 noise = self.getNoisebyHildebrand()
826 826
827 827 if type == 2:
828 828 noise = self.getNoisebySort()
829 829
830 830 if type == 3:
831 831 noise = self.getNoisebyWindow()
832 832
833 833 return noise
834 834
835 835 def getTimeInterval(self):
836 836
837 837 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
838 838
839 839 return timeInterval
840 840
841 841 def get_ippSeconds(self):
842 842 '''
843 843 '''
844 844 return self.ipp_sec
845 845
846 846
847 847 datatime = property(getDatatime, "I'm the 'datatime' property")
848 848 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
849 849 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
850 850 channelIndexList = property(
851 851 getChannelIndexList, "I'm the 'channelIndexList' property.")
852 852 noise = property(getNoise, "I'm the 'nHeights' property.")
853 853
854 854 ltctime = property(getltctime, "I'm the 'ltctime' property")
855 855 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
856 856 ippSeconds = property(get_ippSeconds, '')
857 857
858 858 class Correlation(JROData):
859 859
860 860 noise = None
861 861 SNR = None
862 862 #--------------------------------------------------
863 863 mode = None
864 864 split = False
865 865 data_cf = None
866 866 lags = None
867 867 lagRange = None
868 868 pairsList = None
869 869 normFactor = None
870 870 #--------------------------------------------------
871 871 # calculateVelocity = None
872 872 nLags = None
873 873 nPairs = None
874 874 nAvg = None
875 875
876 876 def __init__(self):
877 877 '''
878 878 Constructor
879 879 '''
880 880 self.radarControllerHeaderObj = RadarControllerHeader()
881 881
882 882 self.systemHeaderObj = SystemHeader()
883 883
884 884 self.type = "Correlation"
885 885
886 886 self.data = None
887 887
888 888 self.dtype = None
889 889
890 890 self.nProfiles = None
891 891
892 892 self.heightList = None
893 893
894 894 self.channelList = None
895 895
896 896 self.flagNoData = True
897 897
898 898 self.flagDiscontinuousBlock = False
899 899
900 900 self.utctime = None
901 901
902 902 self.timeZone = None
903 903
904 904 self.dstFlag = None
905 905
906 906 self.errorCount = None
907 907
908 908 self.blocksize = None
909 909
910 910 self.flagDecodeData = False # asumo q la data no esta decodificada
911 911
912 912 self.flagDeflipData = False # asumo q la data no esta sin flip
913 913
914 914 self.pairsList = None
915 915
916 916 self.nPoints = None
917 917
918 918 def getPairsList(self):
919 919
920 920 return self.pairsList
921 921
922 922 def getNoise(self, mode=2):
923 923
924 924 indR = numpy.where(self.lagR == 0)[0][0]
925 925 indT = numpy.where(self.lagT == 0)[0][0]
926 926
927 927 jspectra0 = self.data_corr[:, :, indR, :]
928 928 jspectra = copy.copy(jspectra0)
929 929
930 930 num_chan = jspectra.shape[0]
931 931 num_hei = jspectra.shape[2]
932 932
933 933 freq_dc = jspectra.shape[1] / 2
934 934 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
935 935
936 936 if ind_vel[0] < 0:
937 937 ind_vel[list(range(0, 1))] = ind_vel[list(
938 938 range(0, 1))] + self.num_prof
939 939
940 940 if mode == 1:
941 941 jspectra[:, freq_dc, :] = (
942 942 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
943 943
944 944 if mode == 2:
945 945
946 946 vel = numpy.array([-2, -1, 1, 2])
947 947 xx = numpy.zeros([4, 4])
948 948
949 949 for fil in range(4):
950 950 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
951 951
952 952 xx_inv = numpy.linalg.inv(xx)
953 953 xx_aux = xx_inv[0, :]
954 954
955 955 for ich in range(num_chan):
956 956 yy = jspectra[ich, ind_vel, :]
957 957 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
958 958
959 959 junkid = jspectra[ich, freq_dc, :] <= 0
960 960 cjunkid = sum(junkid)
961 961
962 962 if cjunkid.any():
963 963 jspectra[ich, freq_dc, junkid.nonzero()] = (
964 964 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
965 965
966 966 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
967 967
968 968 return noise
969 969
970 970 def getTimeInterval(self):
971 971
972 972 timeInterval = self.ippSeconds * self.nCohInt * self.nProfiles
973 973
974 974 return timeInterval
975 975
976 976 def splitFunctions(self):
977 977
978 978 pairsList = self.pairsList
979 979 ccf_pairs = []
980 980 acf_pairs = []
981 981 ccf_ind = []
982 982 acf_ind = []
983 983 for l in range(len(pairsList)):
984 984 chan0 = pairsList[l][0]
985 985 chan1 = pairsList[l][1]
986 986
987 987 # Obteniendo pares de Autocorrelacion
988 988 if chan0 == chan1:
989 989 acf_pairs.append(chan0)
990 990 acf_ind.append(l)
991 991 else:
992 992 ccf_pairs.append(pairsList[l])
993 993 ccf_ind.append(l)
994 994
995 995 data_acf = self.data_cf[acf_ind]
996 996 data_ccf = self.data_cf[ccf_ind]
997 997
998 998 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
999 999
1000 1000 def getNormFactor(self):
1001 1001 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
1002 1002 acf_pairs = numpy.array(acf_pairs)
1003 1003 normFactor = numpy.zeros((self.nPairs, self.nHeights))
1004 1004
1005 1005 for p in range(self.nPairs):
1006 1006 pair = self.pairsList[p]
1007 1007
1008 1008 ch0 = pair[0]
1009 1009 ch1 = pair[1]
1010 1010
1011 1011 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
1012 1012 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
1013 1013 normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
1014 1014
1015 1015 return normFactor
1016 1016
1017 1017 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
1018 1018 normFactor = property(getNormFactor, "I'm the 'normFactor property'")
1019 1019
1020 1020
1021 1021 class Parameters(Spectra):
1022 1022
1023 1023 experimentInfo = None # Information about the experiment
1024 1024 # Information from previous data
1025 1025 inputUnit = None # Type of data to be processed
1026 1026 operation = None # Type of operation to parametrize
1027 1027 # normFactor = None #Normalization Factor
1028 1028 groupList = None # List of Pairs, Groups, etc
1029 1029 # Parameters
1030 1030 data_param = None # Parameters obtained
1031 1031 data_pre = None # Data Pre Parametrization
1032 1032 data_SNR = None # Signal to Noise Ratio
1033 1033 # heightRange = None #Heights
1034 1034 abscissaList = None # Abscissa, can be velocities, lags or time
1035 1035 # noise = None #Noise Potency
1036 1036 utctimeInit = None # Initial UTC time
1037 1037 paramInterval = None # Time interval to calculate Parameters in seconds
1038 1038 useLocalTime = True
1039 1039 # Fitting
1040 1040 data_error = None # Error of the estimation
1041 1041 constants = None
1042 1042 library = None
1043 1043 # Output signal
1044 1044 outputInterval = None # Time interval to calculate output signal in seconds
1045 1045 data_output = None # Out signal
1046 1046 nAvg = None
1047 1047 noise_estimation = None
1048 1048 GauSPC = None # Fit gaussian SPC
1049 1049
1050 1050 def __init__(self):
1051 1051 '''
1052 1052 Constructor
1053 1053 '''
1054 1054 self.radarControllerHeaderObj = RadarControllerHeader()
1055 1055
1056 1056 self.systemHeaderObj = SystemHeader()
1057 1057
1058 1058 self.type = "Parameters"
1059 1059
1060 1060 def getTimeRange1(self, interval):
1061 1061
1062 1062 datatime = []
1063 1063
1064 1064 if self.useLocalTime:
1065 1065 time1 = self.utctimeInit - self.timeZone * 60
1066 1066 else:
1067 1067 time1 = self.utctimeInit
1068 1068
1069 1069 datatime.append(time1)
1070 1070 datatime.append(time1 + interval)
1071 1071 datatime = numpy.array(datatime)
1072 1072
1073 1073 return datatime
1074 1074
1075 1075 def getTimeInterval(self):
1076 1076
1077 1077 if hasattr(self, 'timeInterval1'):
1078 1078 return self.timeInterval1
1079 1079 else:
1080 1080 return self.paramInterval
1081 1081
1082 1082 def setValue(self, value):
1083 1083
1084 1084 print("This property should not be initialized")
1085 1085
1086 1086 return
1087 1087
1088 1088 def getNoise(self):
1089 1089
1090 1090 return self.spc_noise
1091 1091
1092 1092 timeInterval = property(getTimeInterval)
1093 1093 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
1094 1094
1095 1095
1096 1096 class PlotterData(object):
1097 1097 '''
1098 1098 Object to hold data to be plotted
1099 1099 '''
1100 1100
1101 1101 MAXNUMX = 100
1102 1102 MAXNUMY = 100
1103 1103
1104 1104 def __init__(self, code, throttle_value, exp_code, buffering=True):
1105 1105
1106 1106 self.throttle = throttle_value
1107 1107 self.exp_code = exp_code
1108 1108 self.buffering = buffering
1109 1109 self.ready = False
1110 1110 self.localtime = False
1111 1111 self.data = {}
1112 1112 self.meta = {}
1113 1113 self.__times = []
1114 1114 self.__heights = []
1115 1115
1116 1116 if 'snr' in code:
1117 1117 self.plottypes = ['snr']
1118 1118 elif code == 'spc':
1119 1119 self.plottypes = ['spc', 'noise', 'rti']
1120 1120 elif code == 'rti':
1121 1121 self.plottypes = ['noise', 'rti']
1122 1122 else:
1123 1123 self.plottypes = [code]
1124 1124
1125 1125 for plot in self.plottypes:
1126 1126 self.data[plot] = {}
1127 1127
1128 1128 def __str__(self):
1129 1129 dum = ['{}{}'.format(key, self.shape(key)) for key in self.data]
1130 1130 return 'Data[{}][{}]'.format(';'.join(dum), len(self.__times))
1131 1131
1132 1132 def __len__(self):
1133 1133 return len(self.__times)
1134 1134
1135 1135 def __getitem__(self, key):
1136 1136
1137 1137 if key not in self.data:
1138 1138 raise KeyError(log.error('Missing key: {}'.format(key)))
1139 1139 if 'spc' in key or not self.buffering:
1140 1140 ret = self.data[key]
1141 1141 elif 'scope' in key:
1142 1142 ret = numpy.array(self.data[key][float(self.tm)])
1143 1143 else:
1144 1144 ret = numpy.array([self.data[key][x] for x in self.times])
1145 1145 if ret.ndim > 1:
1146 1146 ret = numpy.swapaxes(ret, 0, 1)
1147 1147 return ret
1148 1148
1149 1149 def __contains__(self, key):
1150 1150 return key in self.data
1151 1151
1152 1152 def setup(self):
1153 1153 '''
1154 1154 Configure object
1155 1155 '''
1156 1156
1157 1157 self.type = ''
1158 1158 self.ready = False
1159 1159 self.data = {}
1160 1160 self.__times = []
1161 1161 self.__heights = []
1162 1162 self.__all_heights = set()
1163 1163 for plot in self.plottypes:
1164 1164 if 'snr' in plot:
1165 1165 plot = 'snr'
1166 elif 'spc_moments' == plot:
1167 plot = 'moments'
1166 1168 self.data[plot] = {}
1167 1169
1168 if 'spc' in self.data or 'rti' in self.data or 'cspc' in self.data:
1170 if 'spc' in self.data or 'rti' in self.data or 'cspc' in self.data or 'moments' in self.data:
1169 1171 self.data['noise'] = {}
1170 1172 if 'noise' not in self.plottypes:
1171 1173 self.plottypes.append('noise')
1172 1174
1173 1175 def shape(self, key):
1174 1176 '''
1175 1177 Get the shape of the one-element data for the given key
1176 1178 '''
1177 1179
1178 1180 if len(self.data[key]):
1179 1181 if 'spc' in key or not self.buffering:
1180 1182 return self.data[key].shape
1181 1183 return self.data[key][self.__times[0]].shape
1182 1184 return (0,)
1183 1185
1184 1186 def update(self, dataOut, tm):
1185 1187 '''
1186 1188 Update data object with new dataOut
1187 1189 '''
1188 1190
1189 1191 if tm in self.__times:
1190 1192 return
1191 1193 self.profileIndex = dataOut.profileIndex
1192 1194 self.tm = tm
1193 1195 self.type = dataOut.type
1194 1196 self.parameters = getattr(dataOut, 'parameters', [])
1195 1197 if hasattr(dataOut, 'pairsList'):
1196 1198 self.pairs = dataOut.pairsList
1197 1199 if hasattr(dataOut, 'meta'):
1198 1200 self.meta = dataOut.meta
1199 1201 self.channels = dataOut.channelList
1200 1202 self.interval = dataOut.getTimeInterval()
1201 1203 self.localtime = dataOut.useLocalTime
1202 if 'spc' in self.plottypes or 'cspc' in self.plottypes:
1204 if 'spc' in self.plottypes or 'cspc' in self.plottypes or 'spc_moments' in self.plottypes:
1203 1205 self.xrange = (dataOut.getFreqRange(1)/1000.,
1204 1206 dataOut.getAcfRange(1), dataOut.getVelRange(1))
1205 1207 self.factor = dataOut.normFactor
1206 1208 self.__heights.append(dataOut.heightList)
1207 1209 self.__all_heights.update(dataOut.heightList)
1208 1210 self.__times.append(tm)
1209 1211
1210 1212 for plot in self.plottypes:
1211 if plot == 'spc':
1213 if plot in ('spc', 'spc_moments'):
1212 1214 z = dataOut.data_spc/dataOut.normFactor
1213 1215 buffer = 10*numpy.log10(z)
1214 1216 if plot == 'cspc':
1215 1217 z = dataOut.data_spc/dataOut.normFactor
1216 1218 buffer = (dataOut.data_spc, dataOut.data_cspc)
1217 1219 if plot == 'noise':
1218 1220 buffer = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
1219 1221 if plot == 'rti':
1220 1222 buffer = dataOut.getPower()
1221 1223 if plot == 'snr_db':
1222 1224 buffer = dataOut.data_SNR
1223 1225 if plot == 'snr':
1224 1226 buffer = 10*numpy.log10(dataOut.data_SNR)
1225 1227 if plot == 'dop':
1226 1228 buffer = 10*numpy.log10(dataOut.data_DOP)
1227 1229 if plot == 'mean':
1228 1230 buffer = dataOut.data_MEAN
1229 1231 if plot == 'std':
1230 1232 buffer = dataOut.data_STD
1231 1233 if plot == 'coh':
1232 1234 buffer = dataOut.getCoherence()
1233 1235 if plot == 'phase':
1234 1236 buffer = dataOut.getCoherence(phase=True)
1235 1237 if plot == 'output':
1236 1238 buffer = dataOut.data_output
1237 1239 if plot == 'param':
1238 1240 buffer = dataOut.data_param
1239 1241 if plot == 'scope':
1240 1242 buffer = dataOut.data
1241 1243 self.flagDataAsBlock = dataOut.flagDataAsBlock
1242 1244 self.nProfiles = dataOut.nProfiles
1243 1245
1244 1246 if plot == 'spc':
1245 1247 self.data[plot] = buffer
1246 1248 elif plot == 'cspc':
1247 1249 self.data['spc'] = buffer[0]
1248 1250 self.data['cspc'] = buffer[1]
1251 elif plot == 'spc_moments':
1252 self.data['spc'] = buffer
1253 self.data['moments'][tm] = dataOut.moments
1249 1254 else:
1250 1255 if self.buffering:
1251 1256 self.data[plot][tm] = buffer
1252 1257 else:
1253 1258 self.data[plot] = buffer
1254 1259
1255 1260 def normalize_heights(self):
1256 1261 '''
1257 1262 Ensure same-dimension of the data for different heighList
1258 1263 '''
1259 1264
1260 1265 H = numpy.array(list(self.__all_heights))
1261 1266 H.sort()
1262 1267 for key in self.data:
1263 1268 shape = self.shape(key)[:-1] + H.shape
1264 1269 for tm, obj in list(self.data[key].items()):
1265 1270 h = self.__heights[self.__times.index(tm)]
1266 1271 if H.size == h.size:
1267 1272 continue
1268 1273 index = numpy.where(numpy.in1d(H, h))[0]
1269 1274 dummy = numpy.zeros(shape) + numpy.nan
1270 1275 if len(shape) == 2:
1271 1276 dummy[:, index] = obj
1272 1277 else:
1273 1278 dummy[index] = obj
1274 1279 self.data[key][tm] = dummy
1275 1280
1276 1281 self.__heights = [H for tm in self.__times]
1277 1282
1278 1283 def jsonify(self, decimate=False):
1279 1284 '''
1280 1285 Convert data to json
1281 1286 '''
1282 1287
1283 1288 data = {}
1284 1289 tm = self.times[-1]
1285 1290 dy = int(self.heights.size/self.MAXNUMY) + 1
1286 1291 for key in self.data:
1287 1292 if key in ('spc', 'cspc') or not self.buffering:
1288 1293 dx = int(self.data[key].shape[1]/self.MAXNUMX) + 1
1289 1294 data[key] = self.roundFloats(
1290 1295 self.data[key][::, ::dx, ::dy].tolist())
1291 1296 else:
1292 1297 data[key] = self.roundFloats(self.data[key][tm].tolist())
1293 1298
1294 1299 ret = {'data': data}
1295 1300 ret['exp_code'] = self.exp_code
1296 1301 ret['time'] = float(tm)
1297 1302 ret['interval'] = float(self.interval)
1298 1303 ret['localtime'] = self.localtime
1299 1304 ret['yrange'] = self.roundFloats(self.heights[::dy].tolist())
1300 1305 if 'spc' in self.data or 'cspc' in self.data:
1301 1306 ret['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
1302 1307 else:
1303 1308 ret['xrange'] = []
1304 1309 if hasattr(self, 'pairs'):
1305 1310 ret['pairs'] = [(int(p[0]), int(p[1])) for p in self.pairs]
1306 1311 else:
1307 1312 ret['pairs'] = []
1308 1313
1309 1314 for key, value in list(self.meta.items()):
1310 1315 ret[key] = value
1311 1316
1312 1317 return json.dumps(ret)
1313 1318
1314 1319 @property
1315 1320 def times(self):
1316 1321 '''
1317 1322 Return the list of times of the current data
1318 1323 '''
1319 1324
1320 1325 ret = numpy.array(self.__times)
1321 1326 ret.sort()
1322 1327 return ret
1323 1328
1324 1329 @property
1325 1330 def min_time(self):
1326 1331 '''
1327 1332 Return the minimun time value
1328 1333 '''
1329 1334
1330 1335 return self.times[0]
1331 1336
1332 1337 @property
1333 1338 def max_time(self):
1334 1339 '''
1335 1340 Return the maximun time value
1336 1341 '''
1337 1342
1338 1343 return self.times[-1]
1339 1344
1340 1345 @property
1341 1346 def heights(self):
1342 1347 '''
1343 1348 Return the list of heights of the current data
1344 1349 '''
1345 1350
1346 1351 return numpy.array(self.__heights[-1])
1347 1352
1348 1353 @staticmethod
1349 1354 def roundFloats(obj):
1350 1355 if isinstance(obj, list):
1351 1356 return list(map(PlotterData.roundFloats, obj))
1352 1357 elif isinstance(obj, float):
1353 1358 return round(obj, 2)
@@ -1,748 +1,747
1 1 '''
2 2 New Plots Operations
3 3
4 4 @author: juan.espinoza@jro.igp.gob.pe
5 5 '''
6 6
7 7
8 8 import time
9 9 import datetime
10 10 import numpy
11 11
12 12 from schainpy.model.graphics.jroplot_base import Plot, plt
13 13 from schainpy.utils import log
14 14
15 15 EARTH_RADIUS = 6.3710e3
16 16
17 17
18 18 def ll2xy(lat1, lon1, lat2, lon2):
19 19
20 20 p = 0.017453292519943295
21 21 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
22 22 numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
23 23 r = 12742 * numpy.arcsin(numpy.sqrt(a))
24 24 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
25 25 * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
26 26 theta = -theta + numpy.pi/2
27 27 return r*numpy.cos(theta), r*numpy.sin(theta)
28 28
29 29
30 30 def km2deg(km):
31 31 '''
32 32 Convert distance in km to degrees
33 33 '''
34 34
35 35 return numpy.rad2deg(km/EARTH_RADIUS)
36 36
37 37
38 38 class SpectraPlot(Plot):
39 39 '''
40 40 Plot for Spectra data
41 41 '''
42 42
43 43 CODE = 'spc'
44 44 colormap = 'jro'
45 45
46 46 def setup(self):
47 47 self.nplots = len(self.data.channels)
48 48 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
49 49 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
50 50 self.width = 3.4 * self.ncols
51 51 self.height = 3 * self.nrows
52 52 self.cb_label = 'dB'
53 53 if self.showprofile:
54 54 self.width += 0.8 * self.ncols
55 55
56 56 self.ylabel = 'Range [km]'
57 57
58 58 def plot(self):
59 59 if self.xaxis == "frequency":
60 60 x = self.data.xrange[0]
61 61 self.xlabel = "Frequency (kHz)"
62 62 elif self.xaxis == "time":
63 63 x = self.data.xrange[1]
64 64 self.xlabel = "Time (ms)"
65 65 else:
66 66 x = self.data.xrange[2]
67 67 self.xlabel = "Velocity (m/s)"
68 68
69 if self.CODE == 'spc_mean':
69 if self.CODE == 'spc_moments':
70 70 x = self.data.xrange[2]
71 71 self.xlabel = "Velocity (m/s)"
72 72
73 73 self.titles = []
74 74
75 75 y = self.data.heights
76 76 self.y = y
77 77 z = self.data['spc']
78 78
79 79 for n, ax in enumerate(self.axes):
80 80 noise = self.data['noise'][n][-1]
81 if self.CODE == 'spc_mean':
82 mean = self.data['mean'][n][-1]
81 if self.CODE == 'spc_moments':
82 mean = self.data['moments'][n, :, 1, :][-1]
83 83 if ax.firsttime:
84 84 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
85 85 self.xmin = self.xmin if self.xmin else -self.xmax
86 86 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
87 87 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
88 88 ax.plt = ax.pcolormesh(x, y, z[n].T,
89 89 vmin=self.zmin,
90 90 vmax=self.zmax,
91 91 cmap=plt.get_cmap(self.colormap)
92 92 )
93 93
94 94 if self.showprofile:
95 95 ax.plt_profile = self.pf_axes[n].plot(
96 96 self.data['rti'][n][-1], y)[0]
97 97 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
98 98 color="k", linestyle="dashed", lw=1)[0]
99 if self.CODE == 'spc_mean':
99 if self.CODE == 'spc_moments':
100 100 ax.plt_mean = ax.plot(mean, y, color='k')[0]
101 101 else:
102 102 ax.plt.set_array(z[n].T.ravel())
103 103 if self.showprofile:
104 104 ax.plt_profile.set_data(self.data['rti'][n][-1], y)
105 105 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
106 if self.CODE == 'spc_mean':
106 if self.CODE == 'spc_moments':
107 107 ax.plt_mean.set_data(mean, y)
108
109 108 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
110 109
111 110
112 111 class CrossSpectraPlot(Plot):
113 112
114 113 CODE = 'cspc'
115 114 colormap = 'jet'
116 115 zmin_coh = None
117 116 zmax_coh = None
118 117 zmin_phase = None
119 118 zmax_phase = None
120 119
121 120 def setup(self):
122 121
123 122 self.ncols = 4
124 123 self.nrows = len(self.data.pairs)
125 124 self.nplots = self.nrows * 4
126 125 self.width = 3.4 * self.ncols
127 126 self.height = 3 * self.nrows
128 127 self.ylabel = 'Range [km]'
129 128 self.showprofile = False
130 129
131 130 def plot(self):
132 131
133 132 if self.xaxis == "frequency":
134 133 x = self.data.xrange[0]
135 134 self.xlabel = "Frequency (kHz)"
136 135 elif self.xaxis == "time":
137 136 x = self.data.xrange[1]
138 137 self.xlabel = "Time (ms)"
139 138 else:
140 139 x = self.data.xrange[2]
141 140 self.xlabel = "Velocity (m/s)"
142 141
143 142 self.titles = []
144 143
145 144 y = self.data.heights
146 145 self.y = y
147 146 spc = self.data['spc']
148 147 cspc = self.data['cspc']
149 148
150 149 for n in range(self.nrows):
151 150 noise = self.data['noise'][n][-1]
152 151 pair = self.data.pairs[n]
153 152 ax = self.axes[4 * n]
154 153 spc0 = 10.*numpy.log10(spc[pair[0]]/self.data.factor)
155 154 if ax.firsttime:
156 155 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
157 156 self.xmin = self.xmin if self.xmin else -self.xmax
158 157 self.zmin = self.zmin if self.zmin else numpy.nanmin(spc)
159 158 self.zmax = self.zmax if self.zmax else numpy.nanmax(spc)
160 159 ax.plt = ax.pcolormesh(x , y , spc0.T,
161 160 vmin=self.zmin,
162 161 vmax=self.zmax,
163 162 cmap=plt.get_cmap(self.colormap)
164 163 )
165 164 else:
166 165 ax.plt.set_array(spc0.T.ravel())
167 166 self.titles.append('CH {}: {:3.2f}dB'.format(pair[0], noise))
168 167
169 168 ax = self.axes[4 * n + 1]
170 169 spc1 = 10.*numpy.log10(spc[pair[1]]/self.data.factor)
171 170 if ax.firsttime:
172 171 ax.plt = ax.pcolormesh(x , y, spc1.T,
173 172 vmin=self.zmin,
174 173 vmax=self.zmax,
175 174 cmap=plt.get_cmap(self.colormap)
176 175 )
177 176 else:
178 177 ax.plt.set_array(spc1.T.ravel())
179 178 self.titles.append('CH {}: {:3.2f}dB'.format(pair[1], noise))
180 179
181 180 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
182 181 coh = numpy.abs(out)
183 182 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
184 183
185 184 ax = self.axes[4 * n + 2]
186 185 if ax.firsttime:
187 186 ax.plt = ax.pcolormesh(x, y, coh.T,
188 187 vmin=0,
189 188 vmax=1,
190 189 cmap=plt.get_cmap(self.colormap_coh)
191 190 )
192 191 else:
193 192 ax.plt.set_array(coh.T.ravel())
194 193 self.titles.append(
195 194 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
196 195
197 196 ax = self.axes[4 * n + 3]
198 197 if ax.firsttime:
199 198 ax.plt = ax.pcolormesh(x, y, phase.T,
200 199 vmin=-180,
201 200 vmax=180,
202 201 cmap=plt.get_cmap(self.colormap_phase)
203 202 )
204 203 else:
205 204 ax.plt.set_array(phase.T.ravel())
206 205 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
207 206
208 207
209 class SpectraMeanPlot(SpectraPlot):
208 class SpectralMomentsPlot(SpectraPlot):
210 209 '''
211 Plot for Spectra and Mean
210 Plot for Spectral Moments
212 211 '''
213 CODE = 'spc_mean'
212 CODE = 'spc_moments'
214 213 colormap = 'jro'
215 214
216 215
217 216 class RTIPlot(Plot):
218 217 '''
219 218 Plot for RTI data
220 219 '''
221 220
222 221 CODE = 'rti'
223 222 colormap = 'jro'
224 223
225 224 def setup(self):
226 225 self.xaxis = 'time'
227 226 self.ncols = 1
228 227 self.nrows = len(self.data.channels)
229 228 self.nplots = len(self.data.channels)
230 229 self.ylabel = 'Range [km]'
231 230 self.cb_label = 'dB'
232 231 self.titles = ['{} Channel {}'.format(
233 232 self.CODE.upper(), x) for x in range(self.nrows)]
234 233
235 234 def plot(self):
236 235 self.x = self.data.times
237 236 self.y = self.data.heights
238 237 self.z = self.data[self.CODE]
239 238 self.z = numpy.ma.masked_invalid(self.z)
240 239
241 240 if self.decimation is None:
242 241 x, y, z = self.fill_gaps(self.x, self.y, self.z)
243 242 else:
244 243 x, y, z = self.fill_gaps(*self.decimate())
245 244
246 245 for n, ax in enumerate(self.axes):
247 246 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
248 247 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
249 248 if ax.firsttime:
250 249 ax.plt = ax.pcolormesh(x, y, z[n].T,
251 250 vmin=self.zmin,
252 251 vmax=self.zmax,
253 252 cmap=plt.get_cmap(self.colormap)
254 253 )
255 254 if self.showprofile:
256 255 ax.plot_profile = self.pf_axes[n].plot(
257 256 self.data['rti'][n][-1], self.y)[0]
258 257 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(self.data['noise'][n][-1], len(self.y)), self.y,
259 258 color="k", linestyle="dashed", lw=1)[0]
260 259 else:
261 260 ax.collections.remove(ax.collections[0])
262 261 ax.plt = ax.pcolormesh(x, y, z[n].T,
263 262 vmin=self.zmin,
264 263 vmax=self.zmax,
265 264 cmap=plt.get_cmap(self.colormap)
266 265 )
267 266 if self.showprofile:
268 267 ax.plot_profile.set_data(self.data['rti'][n][-1], self.y)
269 268 ax.plot_noise.set_data(numpy.repeat(
270 269 self.data['noise'][n][-1], len(self.y)), self.y)
271 270
272 271
273 272 class CoherencePlot(RTIPlot):
274 273 '''
275 274 Plot for Coherence data
276 275 '''
277 276
278 277 CODE = 'coh'
279 278
280 279 def setup(self):
281 280 self.xaxis = 'time'
282 281 self.ncols = 1
283 282 self.nrows = len(self.data.pairs)
284 283 self.nplots = len(self.data.pairs)
285 284 self.ylabel = 'Range [km]'
286 285 if self.CODE == 'coh':
287 286 self.cb_label = ''
288 287 self.titles = [
289 288 'Coherence Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
290 289 else:
291 290 self.cb_label = 'Degrees'
292 291 self.titles = [
293 292 'Phase Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
294 293
295 294
296 295 class PhasePlot(CoherencePlot):
297 296 '''
298 297 Plot for Phase map data
299 298 '''
300 299
301 300 CODE = 'phase'
302 301 colormap = 'seismic'
303 302
304 303
305 304 class NoisePlot(Plot):
306 305 '''
307 306 Plot for noise
308 307 '''
309 308
310 309 CODE = 'noise'
311 310
312 311 def setup(self):
313 312 self.xaxis = 'time'
314 313 self.ncols = 1
315 314 self.nrows = 1
316 315 self.nplots = 1
317 316 self.ylabel = 'Intensity [dB]'
318 317 self.titles = ['Noise']
319 318 self.colorbar = False
320 319
321 320 def plot(self):
322 321
323 322 x = self.data.times
324 323 xmin = self.data.min_time
325 324 xmax = xmin + self.xrange * 60 * 60
326 325 Y = self.data[self.CODE]
327 326
328 327 if self.axes[0].firsttime:
329 328 for ch in self.data.channels:
330 329 y = Y[ch]
331 330 self.axes[0].plot(x, y, lw=1, label='Ch{}'.format(ch))
332 331 plt.legend()
333 332 else:
334 333 for ch in self.data.channels:
335 334 y = Y[ch]
336 335 self.axes[0].lines[ch].set_data(x, y)
337 336
338 337 self.ymin = numpy.nanmin(Y) - 5
339 338 self.ymax = numpy.nanmax(Y) + 5
340 339
341 340
342 341 class SnrPlot(RTIPlot):
343 342 '''
344 343 Plot for SNR Data
345 344 '''
346 345
347 346 CODE = 'snr'
348 347 colormap = 'jet'
349 348
350 349
351 350 class DopplerPlot(RTIPlot):
352 351 '''
353 352 Plot for DOPPLER Data
354 353 '''
355 354
356 355 CODE = 'dop'
357 356 colormap = 'jet'
358 357
359 358
360 359 class SkyMapPlot(Plot):
361 360 '''
362 361 Plot for meteors detection data
363 362 '''
364 363
365 364 CODE = 'param'
366 365
367 366 def setup(self):
368 367
369 368 self.ncols = 1
370 369 self.nrows = 1
371 370 self.width = 7.2
372 371 self.height = 7.2
373 372 self.nplots = 1
374 373 self.xlabel = 'Zonal Zenith Angle (deg)'
375 374 self.ylabel = 'Meridional Zenith Angle (deg)'
376 375 self.polar = True
377 376 self.ymin = -180
378 377 self.ymax = 180
379 378 self.colorbar = False
380 379
381 380 def plot(self):
382 381
383 382 arrayParameters = numpy.concatenate(self.data['param'])
384 383 error = arrayParameters[:, -1]
385 384 indValid = numpy.where(error == 0)[0]
386 385 finalMeteor = arrayParameters[indValid, :]
387 386 finalAzimuth = finalMeteor[:, 3]
388 387 finalZenith = finalMeteor[:, 4]
389 388
390 389 x = finalAzimuth * numpy.pi / 180
391 390 y = finalZenith
392 391
393 392 ax = self.axes[0]
394 393
395 394 if ax.firsttime:
396 395 ax.plot = ax.plot(x, y, 'bo', markersize=5)[0]
397 396 else:
398 397 ax.plot.set_data(x, y)
399 398
400 399 dt1 = self.getDateTime(self.data.min_time).strftime('%y/%m/%d %H:%M:%S')
401 400 dt2 = self.getDateTime(self.data.max_time).strftime('%y/%m/%d %H:%M:%S')
402 401 title = 'Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n' % (dt1,
403 402 dt2,
404 403 len(x))
405 404 self.titles[0] = title
406 405
407 406
408 407 class ParametersPlot(RTIPlot):
409 408 '''
410 409 Plot for data_param object
411 410 '''
412 411
413 412 CODE = 'param'
414 413 colormap = 'seismic'
415 414
416 415 def setup(self):
417 416 self.xaxis = 'time'
418 417 self.ncols = 1
419 418 self.nrows = self.data.shape(self.CODE)[0]
420 419 self.nplots = self.nrows
421 420 if self.showSNR:
422 421 self.nrows += 1
423 422 self.nplots += 1
424 423
425 424 self.ylabel = 'Height [km]'
426 425 if not self.titles:
427 426 self.titles = self.data.parameters \
428 427 if self.data.parameters else ['Param {}'.format(x) for x in range(self.nrows)]
429 428 if self.showSNR:
430 429 self.titles.append('SNR')
431 430
432 431 def plot(self):
433 432 self.data.normalize_heights()
434 433 self.x = self.data.times
435 434 self.y = self.data.heights
436 435 if self.showSNR:
437 436 self.z = numpy.concatenate(
438 437 (self.data[self.CODE], self.data['snr'])
439 438 )
440 439 else:
441 440 self.z = self.data[self.CODE]
442 441
443 442 self.z = numpy.ma.masked_invalid(self.z)
444 443
445 444 if self.decimation is None:
446 445 x, y, z = self.fill_gaps(self.x, self.y, self.z)
447 446 else:
448 447 x, y, z = self.fill_gaps(*self.decimate())
449 448
450 449 for n, ax in enumerate(self.axes):
451 450
452 451 self.zmax = self.zmax if self.zmax is not None else numpy.max(
453 452 self.z[n])
454 453 self.zmin = self.zmin if self.zmin is not None else numpy.min(
455 454 self.z[n])
456 455
457 456 if ax.firsttime:
458 457 if self.zlimits is not None:
459 458 self.zmin, self.zmax = self.zlimits[n]
460 459
461 460 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
462 461 vmin=self.zmin,
463 462 vmax=self.zmax,
464 463 cmap=self.cmaps[n]
465 464 )
466 465 else:
467 466 if self.zlimits is not None:
468 467 self.zmin, self.zmax = self.zlimits[n]
469 468 ax.collections.remove(ax.collections[0])
470 469 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
471 470 vmin=self.zmin,
472 471 vmax=self.zmax,
473 472 cmap=self.cmaps[n]
474 473 )
475 474
476 475
477 476 class OutputPlot(ParametersPlot):
478 477 '''
479 478 Plot data_output object
480 479 '''
481 480
482 481 CODE = 'output'
483 482 colormap = 'seismic'
484 483
485 484
486 485 class PolarMapPlot(Plot):
487 486 '''
488 487 Plot for weather radar
489 488 '''
490 489
491 490 CODE = 'param'
492 491 colormap = 'seismic'
493 492
494 493 def setup(self):
495 494 self.ncols = 1
496 495 self.nrows = 1
497 496 self.width = 9
498 497 self.height = 8
499 498 self.mode = self.data.meta['mode']
500 499 if self.channels is not None:
501 500 self.nplots = len(self.channels)
502 501 self.nrows = len(self.channels)
503 502 else:
504 503 self.nplots = self.data.shape(self.CODE)[0]
505 504 self.nrows = self.nplots
506 505 self.channels = list(range(self.nplots))
507 506 if self.mode == 'E':
508 507 self.xlabel = 'Longitude'
509 508 self.ylabel = 'Latitude'
510 509 else:
511 510 self.xlabel = 'Range (km)'
512 511 self.ylabel = 'Height (km)'
513 512 self.bgcolor = 'white'
514 513 self.cb_labels = self.data.meta['units']
515 514 self.lat = self.data.meta['latitude']
516 515 self.lon = self.data.meta['longitude']
517 516 self.xmin, self.xmax = float(
518 517 km2deg(self.xmin) + self.lon), float(km2deg(self.xmax) + self.lon)
519 518 self.ymin, self.ymax = float(
520 519 km2deg(self.ymin) + self.lat), float(km2deg(self.ymax) + self.lat)
521 520 # self.polar = True
522 521
523 522 def plot(self):
524 523
525 524 for n, ax in enumerate(self.axes):
526 525 data = self.data['param'][self.channels[n]]
527 526
528 527 zeniths = numpy.linspace(
529 528 0, self.data.meta['max_range'], data.shape[1])
530 529 if self.mode == 'E':
531 530 azimuths = -numpy.radians(self.data.heights)+numpy.pi/2
532 531 r, theta = numpy.meshgrid(zeniths, azimuths)
533 532 x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])), r*numpy.sin(
534 533 theta)*numpy.cos(numpy.radians(self.data.meta['elevation']))
535 534 x = km2deg(x) + self.lon
536 535 y = km2deg(y) + self.lat
537 536 else:
538 537 azimuths = numpy.radians(self.data.heights)
539 538 r, theta = numpy.meshgrid(zeniths, azimuths)
540 539 x, y = r*numpy.cos(theta), r*numpy.sin(theta)
541 540 self.y = zeniths
542 541
543 542 if ax.firsttime:
544 543 if self.zlimits is not None:
545 544 self.zmin, self.zmax = self.zlimits[n]
546 545 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
547 546 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
548 547 vmin=self.zmin,
549 548 vmax=self.zmax,
550 549 cmap=self.cmaps[n])
551 550 else:
552 551 if self.zlimits is not None:
553 552 self.zmin, self.zmax = self.zlimits[n]
554 553 ax.collections.remove(ax.collections[0])
555 554 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
556 555 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
557 556 vmin=self.zmin,
558 557 vmax=self.zmax,
559 558 cmap=self.cmaps[n])
560 559
561 560 if self.mode == 'A':
562 561 continue
563 562
564 563 # plot district names
565 564 f = open('/data/workspace/schain_scripts/distrito.csv')
566 565 for line in f:
567 566 label, lon, lat = [s.strip() for s in line.split(',') if s]
568 567 lat = float(lat)
569 568 lon = float(lon)
570 569 # ax.plot(lon, lat, '.b', ms=2)
571 570 ax.text(lon, lat, label.decode('utf8'), ha='center',
572 571 va='bottom', size='8', color='black')
573 572
574 573 # plot limites
575 574 limites = []
576 575 tmp = []
577 576 for line in open('/data/workspace/schain_scripts/lima.csv'):
578 577 if '#' in line:
579 578 if tmp:
580 579 limites.append(tmp)
581 580 tmp = []
582 581 continue
583 582 values = line.strip().split(',')
584 583 tmp.append((float(values[0]), float(values[1])))
585 584 for points in limites:
586 585 ax.add_patch(
587 586 Polygon(points, ec='k', fc='none', ls='--', lw=0.5))
588 587
589 588 # plot Cuencas
590 589 for cuenca in ('rimac', 'lurin', 'mala', 'chillon', 'chilca', 'chancay-huaral'):
591 590 f = open('/data/workspace/schain_scripts/{}.csv'.format(cuenca))
592 591 values = [line.strip().split(',') for line in f]
593 592 points = [(float(s[0]), float(s[1])) for s in values]
594 593 ax.add_patch(Polygon(points, ec='b', fc='none'))
595 594
596 595 # plot grid
597 596 for r in (15, 30, 45, 60):
598 597 ax.add_artist(plt.Circle((self.lon, self.lat),
599 598 km2deg(r), color='0.6', fill=False, lw=0.2))
600 599 ax.text(
601 600 self.lon + (km2deg(r))*numpy.cos(60*numpy.pi/180),
602 601 self.lat + (km2deg(r))*numpy.sin(60*numpy.pi/180),
603 602 '{}km'.format(r),
604 603 ha='center', va='bottom', size='8', color='0.6', weight='heavy')
605 604
606 605 if self.mode == 'E':
607 606 title = 'El={}$^\circ$'.format(self.data.meta['elevation'])
608 607 label = 'E{:02d}'.format(int(self.data.meta['elevation']))
609 608 else:
610 609 title = 'Az={}$^\circ$'.format(self.data.meta['azimuth'])
611 610 label = 'A{:02d}'.format(int(self.data.meta['azimuth']))
612 611
613 612 self.save_labels = ['{}-{}'.format(lbl, label) for lbl in self.labels]
614 613 self.titles = ['{} {}'.format(
615 614 self.data.parameters[x], title) for x in self.channels]
616 615
617 616 class ScopePlot(Plot):
618 617
619 618 '''
620 619 Plot for Scope
621 620 '''
622 621
623 622 CODE = 'scope'
624 623
625 624 def setup(self):
626 625
627 626 self.xaxis = 'Range (Km)'
628 627 self.ncols = 1
629 628 self.nrows = 1
630 629 self.nplots = 1
631 630 self.ylabel = 'Intensity [dB]'
632 631 self.titles = ['Scope']
633 632 self.colorbar = False
634 633 colspan = 3
635 634 rowspan = 1
636 635
637 636 def plot_iq(self, x, y, channelIndexList, thisDatetime, wintitle):
638 637
639 638 yreal = y[channelIndexList,:].real
640 639 yimag = y[channelIndexList,:].imag
641 640 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y"))
642 641 self.xlabel = "Range (Km)"
643 642 self.ylabel = "Intensity - IQ"
644 643
645 644 self.y = yreal
646 645 self.x = x
647 646 self.xmin = min(x)
648 647 self.xmax = max(x)
649 648
650 649
651 650 self.titles[0] = title
652 651
653 652 for i,ax in enumerate(self.axes):
654 653 title = "Channel %d" %(i)
655 654 if ax.firsttime:
656 655 ax.plt_r = ax.plot(x, yreal[i,:], color='b')[0]
657 656 ax.plt_i = ax.plot(x, yimag[i,:], color='r')[0]
658 657 else:
659 658 #pass
660 659 ax.plt_r.set_data(x, yreal[i,:])
661 660 ax.plt_i.set_data(x, yimag[i,:])
662 661
663 662 def plot_power(self, x, y, channelIndexList, thisDatetime, wintitle):
664 663 y = y[channelIndexList,:] * numpy.conjugate(y[channelIndexList,:])
665 664 yreal = y.real
666 665 self.y = yreal
667 666 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y"))
668 667 self.xlabel = "Range (Km)"
669 668 self.ylabel = "Intensity"
670 669 self.xmin = min(x)
671 670 self.xmax = max(x)
672 671
673 672
674 673 self.titles[0] = title
675 674
676 675 for i,ax in enumerate(self.axes):
677 676 title = "Channel %d" %(i)
678 677
679 678 ychannel = yreal[i,:]
680 679
681 680 if ax.firsttime:
682 681 ax.plt_r = ax.plot(x, ychannel)[0]
683 682 else:
684 683 #pass
685 684 ax.plt_r.set_data(x, ychannel)
686 685
687 686
688 687 def plot(self):
689 688
690 689 if self.channels:
691 690 channels = self.channels
692 691 else:
693 692 channels = self.data.channels
694 693
695 694
696 695
697 696 thisDatetime = datetime.datetime.utcfromtimestamp(self.data.times[-1])
698 697
699 698 scope = self.data['scope']
700 699
701 700
702 701 if self.data.flagDataAsBlock:
703 702
704 703 for i in range(self.data.nProfiles):
705 704
706 705 wintitle1 = " [Profile = %d] " %i
707 706
708 707 if self.type == "power":
709 708 self.plot_power(self.data.heights,
710 709 scope[:,i,:],
711 710 channels,
712 711 thisDatetime,
713 712 wintitle1
714 713 )
715 714
716 715 if self.type == "iq":
717 716 self.plot_iq(self.data.heights,
718 717 scope[:,i,:],
719 718 channels,
720 719 thisDatetime,
721 720 wintitle1
722 721 )
723 722
724 723
725 724
726 725
727 726
728 727 else:
729 728 wintitle = " [Profile = %d] " %self.data.profileIndex
730 729
731 730 if self.type == "power":
732 731 self.plot_power(self.data.heights,
733 732 scope,
734 733 channels,
735 734 thisDatetime,
736 735 wintitle
737 736 )
738 737
739 738 if self.type == "iq":
740 739 self.plot_iq(self.data.heights,
741 740 scope,
742 741 channels,
743 742 thisDatetime,
744 743 wintitle
745 744 )
746 745
747 746
748 747 No newline at end of file
1 NO CONTENT: modified file
@@ -1,1588 +1,1589
1 1 '''
2 2 Created on Jul 9, 2014
3 3
4 4 @author: roj-idl71
5 5 '''
6 6 import os
7 7 import datetime
8 8 import numpy
9 9
10 10 from .figure import Figure, isRealtime, isTimeInHourRange
11 11 from .plotting_codes import *
12 12 from schainpy.model.proc.jroproc_base import MPDecorator
13 13
14 14 from schainpy.utils import log
15 15
16 16 @MPDecorator
17 17 class SpectraPlot_(Figure):
18 18
19 19 isConfig = None
20 20 __nsubplots = None
21 21
22 22 WIDTHPROF = None
23 23 HEIGHTPROF = None
24 24 PREFIX = 'spc'
25 25
26 26 def __init__(self):
27 27 Figure.__init__(self)
28 28 self.isConfig = False
29 29 self.__nsubplots = 1
30 30 self.WIDTH = 250
31 31 self.HEIGHT = 250
32 32 self.WIDTHPROF = 120
33 33 self.HEIGHTPROF = 0
34 34 self.counter_imagwr = 0
35 35
36 36 self.PLOT_CODE = SPEC_CODE
37 37
38 38 self.FTP_WEI = None
39 39 self.EXP_CODE = None
40 40 self.SUB_EXP_CODE = None
41 41 self.PLOT_POS = None
42 42
43 43 self.__xfilter_ena = False
44 44 self.__yfilter_ena = False
45 45
46 46 self.indice=1
47 47
48 48 def getSubplots(self):
49 49
50 50 ncol = int(numpy.sqrt(self.nplots)+0.9)
51 51 nrow = int(self.nplots*1./ncol + 0.9)
52 52
53 53 return nrow, ncol
54 54
55 55 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
56 56
57 57 self.__showprofile = showprofile
58 58 self.nplots = nplots
59 59
60 60 ncolspan = 1
61 61 colspan = 1
62 62 if showprofile:
63 63 ncolspan = 3
64 64 colspan = 2
65 65 self.__nsubplots = 2
66 66
67 67 self.createFigure(id = id,
68 68 wintitle = wintitle,
69 69 widthplot = self.WIDTH + self.WIDTHPROF,
70 70 heightplot = self.HEIGHT + self.HEIGHTPROF,
71 71 show=show)
72 72
73 73 nrow, ncol = self.getSubplots()
74 74
75 75 counter = 0
76 76 for y in range(nrow):
77 77 for x in range(ncol):
78 78
79 79 if counter >= self.nplots:
80 80 break
81 81
82 82 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
83 83
84 84 if showprofile:
85 85 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
86 86
87 87 counter += 1
88 88
89 89 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=True,
90 90 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
91 91 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
92 92 server=None, folder=None, username=None, password=None,
93 93 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False,
94 94 xaxis="frequency", colormap='jet', normFactor=None):
95 95
96 96 """
97 97
98 98 Input:
99 99 dataOut :
100 100 id :
101 101 wintitle :
102 102 channelList :
103 103 showProfile :
104 104 xmin : None,
105 105 xmax : None,
106 106 ymin : None,
107 107 ymax : None,
108 108 zmin : None,
109 109 zmax : None
110 110 """
111 111 if dataOut.flagNoData:
112 112 return dataOut
113 113
114 114 if realtime:
115 115 if not(isRealtime(utcdatatime = dataOut.utctime)):
116 116 print('Skipping this plot function')
117 117 return
118 118
119 119 if channelList == None:
120 120 channelIndexList = dataOut.channelIndexList
121 121 else:
122 122 channelIndexList = []
123 123 for channel in channelList:
124 124 if channel not in dataOut.channelList:
125 125 raise ValueError("Channel %d is not in dataOut.channelList" %channel)
126 126 channelIndexList.append(dataOut.channelList.index(channel))
127 127
128 128 if normFactor is None:
129 129 factor = dataOut.normFactor
130 130 else:
131 131 factor = normFactor
132 132 if xaxis == "frequency":
133 133 x = dataOut.getFreqRange(1)/1000.
134 134 xlabel = "Frequency (kHz)"
135 135
136 136 elif xaxis == "time":
137 137 x = dataOut.getAcfRange(1)
138 138 xlabel = "Time (ms)"
139 139
140 140 else:
141 141 x = dataOut.getVelRange(1)
142 142 xlabel = "Velocity (m/s)"
143 143
144 144 ylabel = "Range (km)"
145 145
146 146 y = dataOut.getHeiRange()
147 147 z = dataOut.data_spc/factor
148 148 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
149 149 zdB = 10*numpy.log10(z)
150 150
151 151 avg = numpy.average(z, axis=1)
152 152 avgdB = 10*numpy.log10(avg)
153 153
154 154 noise = dataOut.getNoise()/factor
155 155 noisedB = 10*numpy.log10(noise)
156 156
157 157 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
158 158 title = wintitle + " Spectra"
159 159
160 160 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
161 161 title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith)
162 162
163 163 if not self.isConfig:
164 164
165 165 nplots = len(channelIndexList)
166 166
167 167 self.setup(id=id,
168 168 nplots=nplots,
169 169 wintitle=wintitle,
170 170 showprofile=showprofile,
171 171 show=show)
172 172
173 173 if xmin == None: xmin = numpy.nanmin(x)
174 174 if xmax == None: xmax = numpy.nanmax(x)
175 175 if ymin == None: ymin = numpy.nanmin(y)
176 176 if ymax == None: ymax = numpy.nanmax(y)
177 177 if zmin == None: zmin = numpy.floor(numpy.nanmin(noisedB)) - 3
178 178 if zmax == None: zmax = numpy.ceil(numpy.nanmax(avgdB)) + 3
179 179
180 180 self.FTP_WEI = ftp_wei
181 181 self.EXP_CODE = exp_code
182 182 self.SUB_EXP_CODE = sub_exp_code
183 183 self.PLOT_POS = plot_pos
184 184
185 185 self.isConfig = True
186 186
187 187 self.setWinTitle(title)
188 188
189 189 for i in range(self.nplots):
190 190 index = channelIndexList[i]
191 191 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
192 192 title = "Channel %d: %4.2fdB: %s" %(dataOut.channelList[index], noisedB[index], str_datetime)
193 193 if len(dataOut.beam.codeList) != 0:
194 194 title = "Ch%d:%4.2fdB,%2.2f,%2.2f:%s" %(dataOut.channelList[index], noisedB[index], dataOut.beam.azimuthList[index], dataOut.beam.zenithList[index], str_datetime)
195 195
196 196 axes = self.axesList[i*self.__nsubplots]
197 197 axes.pcolor(x, y, zdB[index,:,:],
198 198 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
199 199 xlabel=xlabel, ylabel=ylabel, title=title, colormap=colormap,
200 200 ticksize=9, cblabel='')
201 201
202 202 if self.__showprofile:
203 203 axes = self.axesList[i*self.__nsubplots +1]
204 204 axes.pline(avgdB[index,:], y,
205 205 xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
206 206 xlabel='dB', ylabel='', title='',
207 207 ytick_visible=False,
208 208 grid='x')
209 209
210 210 noiseline = numpy.repeat(noisedB[index], len(y))
211 211 axes.addpline(noiseline, y, idline=1, color="black", linestyle="dashed", lw=2)
212 212
213 213 self.draw()
214 214
215 215 if figfile == None:
216 216 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
217 217 name = str_datetime
218 218 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
219 219 name = name + '_az' + '_%2.2f'%(dataOut.azimuth) + '_zn' + '_%2.2f'%(dataOut.zenith)
220 220 figfile = self.getFilename(name)
221 221
222 222 self.save(figpath=figpath,
223 223 figfile=figfile,
224 224 save=save,
225 225 ftp=ftp,
226 226 wr_period=wr_period,
227 227 thisDatetime=thisDatetime)
228 228
229 229
230 230 return dataOut
231
231 232 @MPDecorator
232 233 class CrossSpectraPlot_(Figure):
233 234
234 235 isConfig = None
235 236 __nsubplots = None
236 237
237 238 WIDTH = None
238 239 HEIGHT = None
239 240 WIDTHPROF = None
240 241 HEIGHTPROF = None
241 242 PREFIX = 'cspc'
242 243
243 244 def __init__(self):
244 245 Figure.__init__(self)
245 246 self.isConfig = False
246 247 self.__nsubplots = 4
247 248 self.counter_imagwr = 0
248 249 self.WIDTH = 250
249 250 self.HEIGHT = 250
250 251 self.WIDTHPROF = 0
251 252 self.HEIGHTPROF = 0
252 253
253 254 self.PLOT_CODE = CROSS_CODE
254 255 self.FTP_WEI = None
255 256 self.EXP_CODE = None
256 257 self.SUB_EXP_CODE = None
257 258 self.PLOT_POS = None
258 259
259 260 self.indice=0
260 261
261 262 def getSubplots(self):
262 263
263 264 ncol = 4
264 265 nrow = self.nplots
265 266
266 267 return nrow, ncol
267 268
268 269 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
269 270
270 271 self.__showprofile = showprofile
271 272 self.nplots = nplots
272 273
273 274 ncolspan = 1
274 275 colspan = 1
275 276
276 277 self.createFigure(id = id,
277 278 wintitle = wintitle,
278 279 widthplot = self.WIDTH + self.WIDTHPROF,
279 280 heightplot = self.HEIGHT + self.HEIGHTPROF,
280 281 show=True)
281 282
282 283 nrow, ncol = self.getSubplots()
283 284
284 285 counter = 0
285 286 for y in range(nrow):
286 287 for x in range(ncol):
287 288 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
288 289
289 290 counter += 1
290 291
291 292 def run(self, dataOut, id, wintitle="", pairsList=None,
292 293 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
293 294 coh_min=None, coh_max=None, phase_min=None, phase_max=None,
294 295 save=False, figpath='./', figfile=None, ftp=False, wr_period=1,
295 296 power_cmap='jet', coherence_cmap='jet', phase_cmap='RdBu_r', show=True,
296 297 server=None, folder=None, username=None, password=None,
297 298 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, normFactor=None,
298 299 xaxis='frequency'):
299 300
300 301 """
301 302
302 303 Input:
303 304 dataOut :
304 305 id :
305 306 wintitle :
306 307 channelList :
307 308 showProfile :
308 309 xmin : None,
309 310 xmax : None,
310 311 ymin : None,
311 312 ymax : None,
312 313 zmin : None,
313 314 zmax : None
314 315 """
315 316
316 317 if dataOut.flagNoData:
317 318 return dataOut
318 319
319 320 if pairsList == None:
320 321 pairsIndexList = dataOut.pairsIndexList
321 322 else:
322 323 pairsIndexList = []
323 324 for pair in pairsList:
324 325 if pair not in dataOut.pairsList:
325 326 raise ValueError("Pair %s is not in dataOut.pairsList" %str(pair))
326 327 pairsIndexList.append(dataOut.pairsList.index(pair))
327 328
328 329 if not pairsIndexList:
329 330 return
330 331
331 332 if len(pairsIndexList) > 4:
332 333 pairsIndexList = pairsIndexList[0:4]
333 334
334 335 if normFactor is None:
335 336 factor = dataOut.normFactor
336 337 else:
337 338 factor = normFactor
338 339 x = dataOut.getVelRange(1)
339 340 y = dataOut.getHeiRange()
340 341 z = dataOut.data_spc[:,:,:]/factor
341 342 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
342 343
343 344 noise = dataOut.noise/factor
344 345
345 346 zdB = 10*numpy.log10(z)
346 347 noisedB = 10*numpy.log10(noise)
347 348
348 349 if coh_min == None:
349 350 coh_min = 0.0
350 351 if coh_max == None:
351 352 coh_max = 1.0
352 353
353 354 if phase_min == None:
354 355 phase_min = -180
355 356 if phase_max == None:
356 357 phase_max = 180
357 358
358 359 #thisDatetime = dataOut.datatime
359 360 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
360 361 title = wintitle + " Cross-Spectra: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
361 362 # xlabel = "Velocity (m/s)"
362 363 ylabel = "Range (Km)"
363 364
364 365 if xaxis == "frequency":
365 366 x = dataOut.getFreqRange(1)/1000.
366 367 xlabel = "Frequency (kHz)"
367 368
368 369 elif xaxis == "time":
369 370 x = dataOut.getAcfRange(1)
370 371 xlabel = "Time (ms)"
371 372
372 373 else:
373 374 x = dataOut.getVelRange(1)
374 375 xlabel = "Velocity (m/s)"
375 376
376 377 if not self.isConfig:
377 378
378 379 nplots = len(pairsIndexList)
379 380
380 381 self.setup(id=id,
381 382 nplots=nplots,
382 383 wintitle=wintitle,
383 384 showprofile=False,
384 385 show=show)
385 386
386 387 avg = numpy.abs(numpy.average(z, axis=1))
387 388 avgdB = 10*numpy.log10(avg)
388 389
389 390 if xmin == None: xmin = numpy.nanmin(x)
390 391 if xmax == None: xmax = numpy.nanmax(x)
391 392 if ymin == None: ymin = numpy.nanmin(y)
392 393 if ymax == None: ymax = numpy.nanmax(y)
393 394 if zmin == None: zmin = numpy.floor(numpy.nanmin(noisedB)) - 3
394 395 if zmax == None: zmax = numpy.ceil(numpy.nanmax(avgdB)) + 3
395 396
396 397 self.FTP_WEI = ftp_wei
397 398 self.EXP_CODE = exp_code
398 399 self.SUB_EXP_CODE = sub_exp_code
399 400 self.PLOT_POS = plot_pos
400 401
401 402 self.isConfig = True
402 403
403 404 self.setWinTitle(title)
404 405
405 406
406 407 for i in range(self.nplots):
407 408 pair = dataOut.pairsList[pairsIndexList[i]]
408 409
409 410 chan_index0 = dataOut.channelList.index(pair[0])
410 411 chan_index1 = dataOut.channelList.index(pair[1])
411 412
412 413 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
413 414 title = "Ch%d: %4.2fdB: %s" %(pair[0], noisedB[chan_index0], str_datetime)
414 415 zdB = 10.*numpy.log10(dataOut.data_spc[chan_index0,:,:]/factor)
415 416 axes0 = self.axesList[i*self.__nsubplots]
416 417 axes0.pcolor(x, y, zdB,
417 418 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
418 419 xlabel=xlabel, ylabel=ylabel, title=title,
419 420 ticksize=9, colormap=power_cmap, cblabel='')
420 421
421 422 title = "Ch%d: %4.2fdB: %s" %(pair[1], noisedB[chan_index1], str_datetime)
422 423 zdB = 10.*numpy.log10(dataOut.data_spc[chan_index1,:,:]/factor)
423 424 axes0 = self.axesList[i*self.__nsubplots+1]
424 425 axes0.pcolor(x, y, zdB,
425 426 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
426 427 xlabel=xlabel, ylabel=ylabel, title=title,
427 428 ticksize=9, colormap=power_cmap, cblabel='')
428 429
429 430 coherenceComplex = dataOut.data_cspc[pairsIndexList[i],:,:] / numpy.sqrt( dataOut.data_spc[chan_index0,:,:]*dataOut.data_spc[chan_index1,:,:] )
430 431 coherence = numpy.abs(coherenceComplex)
431 432 # phase = numpy.arctan(-1*coherenceComplex.imag/coherenceComplex.real)*180/numpy.pi
432 433 phase = numpy.arctan2(coherenceComplex.imag, coherenceComplex.real)*180/numpy.pi
433 434
434 435 title = "Coherence Ch%d * Ch%d" %(pair[0], pair[1])
435 436 axes0 = self.axesList[i*self.__nsubplots+2]
436 437 axes0.pcolor(x, y, coherence,
437 438 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=coh_min, zmax=coh_max,
438 439 xlabel=xlabel, ylabel=ylabel, title=title,
439 440 ticksize=9, colormap=coherence_cmap, cblabel='')
440 441
441 442 title = "Phase Ch%d * Ch%d" %(pair[0], pair[1])
442 443 axes0 = self.axesList[i*self.__nsubplots+3]
443 444 axes0.pcolor(x, y, phase,
444 445 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=phase_min, zmax=phase_max,
445 446 xlabel=xlabel, ylabel=ylabel, title=title,
446 447 ticksize=9, colormap=phase_cmap, cblabel='')
447 448
448 449 self.draw()
449 450
450 451 self.save(figpath=figpath,
451 452 figfile=figfile,
452 453 save=save,
453 454 ftp=ftp,
454 455 wr_period=wr_period,
455 456 thisDatetime=thisDatetime)
456 457
457 458 return dataOut
458 459
459 460 @MPDecorator
460 461 class RTIPlot_(Figure):
461 462
462 463 __isConfig = None
463 464 __nsubplots = None
464 465
465 466 WIDTHPROF = None
466 467 HEIGHTPROF = None
467 468 PREFIX = 'rti'
468 469
469 470 def __init__(self):
470 471
471 472 Figure.__init__(self)
472 473 self.timerange = None
473 474 self.isConfig = False
474 475 self.__nsubplots = 1
475 476
476 477 self.WIDTH = 800
477 478 self.HEIGHT = 250
478 479 self.WIDTHPROF = 120
479 480 self.HEIGHTPROF = 0
480 481 self.counter_imagwr = 0
481 482
482 483 self.PLOT_CODE = RTI_CODE
483 484
484 485 self.FTP_WEI = None
485 486 self.EXP_CODE = None
486 487 self.SUB_EXP_CODE = None
487 488 self.PLOT_POS = None
488 489 self.tmin = None
489 490 self.tmax = None
490 491
491 492 self.xmin = None
492 493 self.xmax = None
493 494
494 495 self.figfile = None
495 496
496 497 def getSubplots(self):
497 498
498 499 ncol = 1
499 500 nrow = self.nplots
500 501
501 502 return nrow, ncol
502 503
503 504 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
504 505
505 506 self.__showprofile = showprofile
506 507 self.nplots = nplots
507 508
508 509 ncolspan = 1
509 510 colspan = 1
510 511 if showprofile:
511 512 ncolspan = 7
512 513 colspan = 6
513 514 self.__nsubplots = 2
514 515
515 516 self.createFigure(id = id,
516 517 wintitle = wintitle,
517 518 widthplot = self.WIDTH + self.WIDTHPROF,
518 519 heightplot = self.HEIGHT + self.HEIGHTPROF,
519 520 show=show)
520 521
521 522 nrow, ncol = self.getSubplots()
522 523
523 524 counter = 0
524 525 for y in range(nrow):
525 526 for x in range(ncol):
526 527
527 528 if counter >= self.nplots:
528 529 break
529 530
530 531 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
531 532
532 533 if showprofile:
533 534 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
534 535
535 536 counter += 1
536 537
537 538 def run(self, dataOut, id, wintitle="", channelList=None, showprofile='True',
538 539 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
539 540 timerange=None, colormap='jet',
540 541 save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
541 542 server=None, folder=None, username=None, password=None,
542 543 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, normFactor=None, HEIGHT=None):
543 544
544 545 """
545 546
546 547 Input:
547 548 dataOut :
548 549 id :
549 550 wintitle :
550 551 channelList :
551 552 showProfile :
552 553 xmin : None,
553 554 xmax : None,
554 555 ymin : None,
555 556 ymax : None,
556 557 zmin : None,
557 558 zmax : None
558 559 """
559 560 if dataOut.flagNoData:
560 561 return dataOut
561 562
562 563 #colormap = kwargs.get('colormap', 'jet')
563 564 if HEIGHT is not None:
564 565 self.HEIGHT = HEIGHT
565 566
566 567 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
567 568 return
568 569
569 570 if channelList == None:
570 571 channelIndexList = dataOut.channelIndexList
571 572 else:
572 573 channelIndexList = []
573 574 for channel in channelList:
574 575 if channel not in dataOut.channelList:
575 576 raise ValueError("Channel %d is not in dataOut.channelList")
576 577 channelIndexList.append(dataOut.channelList.index(channel))
577 578
578 579 if normFactor is None:
579 580 factor = dataOut.normFactor
580 581 else:
581 582 factor = normFactor
582 583
583 584 #factor = dataOut.normFactor
584 585 x = dataOut.getTimeRange()
585 586 y = dataOut.getHeiRange()
586 587
587 588 z = dataOut.data_spc/factor
588 589 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
589 590 avg = numpy.average(z, axis=1)
590 591 avgdB = 10.*numpy.log10(avg)
591 592 # avgdB = dataOut.getPower()
592 593
593 594
594 595 thisDatetime = dataOut.datatime
595 596 #thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
596 597 title = wintitle + " RTI" #: %s" %(thisDatetime.strftime("%d-%b-%Y"))
597 598 xlabel = ""
598 599 ylabel = "Range (Km)"
599 600
600 601 update_figfile = False
601 602
602 603 if self.xmax is not None and dataOut.ltctime >= self.xmax: #yong
603 604 self.counter_imagwr = wr_period
604 605 self.isConfig = False
605 606 update_figfile = True
606 607
607 608 if not self.isConfig:
608 609
609 610 nplots = len(channelIndexList)
610 611
611 612 self.setup(id=id,
612 613 nplots=nplots,
613 614 wintitle=wintitle,
614 615 showprofile=showprofile,
615 616 show=show)
616 617
617 618 if timerange != None:
618 619 self.timerange = timerange
619 620
620 621 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
621 622
622 623 noise = dataOut.noise/factor
623 624 noisedB = 10*numpy.log10(noise)
624 625
625 626 if ymin == None: ymin = numpy.nanmin(y)
626 627 if ymax == None: ymax = numpy.nanmax(y)
627 628 if zmin == None: zmin = numpy.floor(numpy.nanmin(noisedB)) - 3
628 629 if zmax == None: zmax = numpy.ceil(numpy.nanmax(avgdB)) + 3
629 630
630 631 self.FTP_WEI = ftp_wei
631 632 self.EXP_CODE = exp_code
632 633 self.SUB_EXP_CODE = sub_exp_code
633 634 self.PLOT_POS = plot_pos
634 635
635 636 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
636 637 self.isConfig = True
637 638 self.figfile = figfile
638 639 update_figfile = True
639 640
640 641 self.setWinTitle(title)
641 642
642 643 for i in range(self.nplots):
643 644 index = channelIndexList[i]
644 645 title = "Channel %d: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
645 646 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
646 647 title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith)
647 648 axes = self.axesList[i*self.__nsubplots]
648 649 zdB = avgdB[index].reshape((1,-1))
649 650 axes.pcolorbuffer(x, y, zdB,
650 651 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
651 652 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
652 653 ticksize=9, cblabel='', cbsize="1%", colormap=colormap)
653 654
654 655 if self.__showprofile:
655 656 axes = self.axesList[i*self.__nsubplots +1]
656 657 axes.pline(avgdB[index], y,
657 658 xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
658 659 xlabel='dB', ylabel='', title='',
659 660 ytick_visible=False,
660 661 grid='x')
661 662
662 663 self.draw()
663 664
664 665 self.save(figpath=figpath,
665 666 figfile=figfile,
666 667 save=save,
667 668 ftp=ftp,
668 669 wr_period=wr_period,
669 670 thisDatetime=thisDatetime,
670 671 update_figfile=update_figfile)
671 672 return dataOut
672 673
673 674 @MPDecorator
674 675 class CoherenceMap_(Figure):
675 676 isConfig = None
676 677 __nsubplots = None
677 678
678 679 WIDTHPROF = None
679 680 HEIGHTPROF = None
680 681 PREFIX = 'cmap'
681 682
682 683 def __init__(self):
683 684 Figure.__init__(self)
684 685 self.timerange = 2*60*60
685 686 self.isConfig = False
686 687 self.__nsubplots = 1
687 688
688 689 self.WIDTH = 800
689 690 self.HEIGHT = 180
690 691 self.WIDTHPROF = 120
691 692 self.HEIGHTPROF = 0
692 693 self.counter_imagwr = 0
693 694
694 695 self.PLOT_CODE = COH_CODE
695 696
696 697 self.FTP_WEI = None
697 698 self.EXP_CODE = None
698 699 self.SUB_EXP_CODE = None
699 700 self.PLOT_POS = None
700 701 self.counter_imagwr = 0
701 702
702 703 self.xmin = None
703 704 self.xmax = None
704 705
705 706 def getSubplots(self):
706 707 ncol = 1
707 708 nrow = self.nplots*2
708 709
709 710 return nrow, ncol
710 711
711 712 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
712 713 self.__showprofile = showprofile
713 714 self.nplots = nplots
714 715
715 716 ncolspan = 1
716 717 colspan = 1
717 718 if showprofile:
718 719 ncolspan = 7
719 720 colspan = 6
720 721 self.__nsubplots = 2
721 722
722 723 self.createFigure(id = id,
723 724 wintitle = wintitle,
724 725 widthplot = self.WIDTH + self.WIDTHPROF,
725 726 heightplot = self.HEIGHT + self.HEIGHTPROF,
726 727 show=True)
727 728
728 729 nrow, ncol = self.getSubplots()
729 730
730 731 for y in range(nrow):
731 732 for x in range(ncol):
732 733
733 734 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
734 735
735 736 if showprofile:
736 737 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
737 738
738 739 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
739 740 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
740 741 timerange=None, phase_min=None, phase_max=None,
741 742 save=False, figpath='./', figfile=None, ftp=False, wr_period=1,
742 743 coherence_cmap='jet', phase_cmap='RdBu_r', show=True,
743 744 server=None, folder=None, username=None, password=None,
744 745 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
745 746
746 747
747 748 if dataOut.flagNoData:
748 749 return dataOut
749 750
750 751 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
751 752 return
752 753
753 754 if pairsList == None:
754 755 pairsIndexList = dataOut.pairsIndexList
755 756 else:
756 757 pairsIndexList = []
757 758 for pair in pairsList:
758 759 if pair not in dataOut.pairsList:
759 760 raise ValueError("Pair %s is not in dataOut.pairsList" %(pair))
760 761 pairsIndexList.append(dataOut.pairsList.index(pair))
761 762
762 763 if pairsIndexList == []:
763 764 return
764 765
765 766 if len(pairsIndexList) > 4:
766 767 pairsIndexList = pairsIndexList[0:4]
767 768
768 769 if phase_min == None:
769 770 phase_min = -180
770 771 if phase_max == None:
771 772 phase_max = 180
772 773
773 774 x = dataOut.getTimeRange()
774 775 y = dataOut.getHeiRange()
775 776
776 777 thisDatetime = dataOut.datatime
777 778
778 779 title = wintitle + " CoherenceMap" #: %s" %(thisDatetime.strftime("%d-%b-%Y"))
779 780 xlabel = ""
780 781 ylabel = "Range (Km)"
781 782 update_figfile = False
782 783
783 784 if not self.isConfig:
784 785 nplots = len(pairsIndexList)
785 786 self.setup(id=id,
786 787 nplots=nplots,
787 788 wintitle=wintitle,
788 789 showprofile=showprofile,
789 790 show=show)
790 791
791 792 if timerange != None:
792 793 self.timerange = timerange
793 794
794 795 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
795 796
796 797 if ymin == None: ymin = numpy.nanmin(y)
797 798 if ymax == None: ymax = numpy.nanmax(y)
798 799 if zmin == None: zmin = 0.
799 800 if zmax == None: zmax = 1.
800 801
801 802 self.FTP_WEI = ftp_wei
802 803 self.EXP_CODE = exp_code
803 804 self.SUB_EXP_CODE = sub_exp_code
804 805 self.PLOT_POS = plot_pos
805 806
806 807 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
807 808
808 809 self.isConfig = True
809 810 update_figfile = True
810 811
811 812 self.setWinTitle(title)
812 813
813 814 for i in range(self.nplots):
814 815
815 816 pair = dataOut.pairsList[pairsIndexList[i]]
816 817
817 818 ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i],:,:],axis=0)
818 819 powa = numpy.average(dataOut.data_spc[pair[0],:,:],axis=0)
819 820 powb = numpy.average(dataOut.data_spc[pair[1],:,:],axis=0)
820 821
821 822
822 823 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
823 824 coherence = numpy.abs(avgcoherenceComplex)
824 825
825 826 z = coherence.reshape((1,-1))
826 827
827 828 counter = 0
828 829
829 830 title = "Coherence Ch%d * Ch%d: %s" %(pair[0], pair[1], thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
830 831 axes = self.axesList[i*self.__nsubplots*2]
831 832 axes.pcolorbuffer(x, y, z,
832 833 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
833 834 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
834 835 ticksize=9, cblabel='', colormap=coherence_cmap, cbsize="1%")
835 836
836 837 if self.__showprofile:
837 838 counter += 1
838 839 axes = self.axesList[i*self.__nsubplots*2 + counter]
839 840 axes.pline(coherence, y,
840 841 xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
841 842 xlabel='', ylabel='', title='', ticksize=7,
842 843 ytick_visible=False, nxticks=5,
843 844 grid='x')
844 845
845 846 counter += 1
846 847
847 848 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
848 849
849 850 z = phase.reshape((1,-1))
850 851
851 852 title = "Phase Ch%d * Ch%d: %s" %(pair[0], pair[1], thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
852 853 axes = self.axesList[i*self.__nsubplots*2 + counter]
853 854 axes.pcolorbuffer(x, y, z,
854 855 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=phase_min, zmax=phase_max,
855 856 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
856 857 ticksize=9, cblabel='', colormap=phase_cmap, cbsize="1%")
857 858
858 859 if self.__showprofile:
859 860 counter += 1
860 861 axes = self.axesList[i*self.__nsubplots*2 + counter]
861 862 axes.pline(phase, y,
862 863 xmin=phase_min, xmax=phase_max, ymin=ymin, ymax=ymax,
863 864 xlabel='', ylabel='', title='', ticksize=7,
864 865 ytick_visible=False, nxticks=4,
865 866 grid='x')
866 867
867 868 self.draw()
868 869
869 870 if dataOut.ltctime >= self.xmax:
870 871 self.counter_imagwr = wr_period
871 872 self.isConfig = False
872 873 update_figfile = True
873 874
874 875 self.save(figpath=figpath,
875 876 figfile=figfile,
876 877 save=save,
877 878 ftp=ftp,
878 879 wr_period=wr_period,
879 880 thisDatetime=thisDatetime,
880 881 update_figfile=update_figfile)
881 882
882 883 return dataOut
883 884
884 885 @MPDecorator
885 886 class PowerProfilePlot_(Figure):
886 887
887 888 isConfig = None
888 889 __nsubplots = None
889 890
890 891 WIDTHPROF = None
891 892 HEIGHTPROF = None
892 893 PREFIX = 'spcprofile'
893 894
894 895 def __init__(self):
895 896 Figure.__init__(self)
896 897 self.isConfig = False
897 898 self.__nsubplots = 1
898 899
899 900 self.PLOT_CODE = POWER_CODE
900 901
901 902 self.WIDTH = 300
902 903 self.HEIGHT = 500
903 904 self.counter_imagwr = 0
904 905
905 906 def getSubplots(self):
906 907 ncol = 1
907 908 nrow = 1
908 909
909 910 return nrow, ncol
910 911
911 912 def setup(self, id, nplots, wintitle, show):
912 913
913 914 self.nplots = nplots
914 915
915 916 ncolspan = 1
916 917 colspan = 1
917 918
918 919 self.createFigure(id = id,
919 920 wintitle = wintitle,
920 921 widthplot = self.WIDTH,
921 922 heightplot = self.HEIGHT,
922 923 show=show)
923 924
924 925 nrow, ncol = self.getSubplots()
925 926
926 927 counter = 0
927 928 for y in range(nrow):
928 929 for x in range(ncol):
929 930 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
930 931
931 932 def run(self, dataOut, id, wintitle="", channelList=None,
932 933 xmin=None, xmax=None, ymin=None, ymax=None,
933 934 save=False, figpath='./', figfile=None, show=True,
934 935 ftp=False, wr_period=1, server=None,
935 936 folder=None, username=None, password=None):
936 937
937 938 if dataOut.flagNoData:
938 939 return dataOut
939 940
940 941
941 942 if channelList == None:
942 943 channelIndexList = dataOut.channelIndexList
943 944 channelList = dataOut.channelList
944 945 else:
945 946 channelIndexList = []
946 947 for channel in channelList:
947 948 if channel not in dataOut.channelList:
948 949 raise ValueError("Channel %d is not in dataOut.channelList")
949 950 channelIndexList.append(dataOut.channelList.index(channel))
950 951
951 952 factor = dataOut.normFactor
952 953
953 954 y = dataOut.getHeiRange()
954 955
955 956 #for voltage
956 957 if dataOut.type == 'Voltage':
957 958 x = dataOut.data[channelIndexList,:] * numpy.conjugate(dataOut.data[channelIndexList,:])
958 959 x = x.real
959 960 x = numpy.where(numpy.isfinite(x), x, numpy.NAN)
960 961
961 962 #for spectra
962 963 if dataOut.type == 'Spectra':
963 964 x = dataOut.data_spc[channelIndexList,:,:]/factor
964 965 x = numpy.where(numpy.isfinite(x), x, numpy.NAN)
965 966 x = numpy.average(x, axis=1)
966 967
967 968
968 969 xdB = 10*numpy.log10(x)
969 970
970 971 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
971 972 title = wintitle + " Power Profile %s" %(thisDatetime.strftime("%d-%b-%Y"))
972 973 xlabel = "dB"
973 974 ylabel = "Range (Km)"
974 975
975 976 if not self.isConfig:
976 977
977 978 nplots = 1
978 979
979 980 self.setup(id=id,
980 981 nplots=nplots,
981 982 wintitle=wintitle,
982 983 show=show)
983 984
984 985 if ymin == None: ymin = numpy.nanmin(y)
985 986 if ymax == None: ymax = numpy.nanmax(y)
986 987 if xmin == None: xmin = numpy.nanmin(xdB)*0.9
987 988 if xmax == None: xmax = numpy.nanmax(xdB)*1.1
988 989
989 990 self.isConfig = True
990 991
991 992 self.setWinTitle(title)
992 993
993 994 title = "Power Profile: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
994 995 axes = self.axesList[0]
995 996
996 997 legendlabels = ["channel %d"%x for x in channelList]
997 998 axes.pmultiline(xdB, y,
998 999 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
999 1000 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels,
1000 1001 ytick_visible=True, nxticks=5,
1001 1002 grid='x')
1002 1003
1003 1004 self.draw()
1004 1005
1005 1006 self.save(figpath=figpath,
1006 1007 figfile=figfile,
1007 1008 save=save,
1008 1009 ftp=ftp,
1009 1010 wr_period=wr_period,
1010 1011 thisDatetime=thisDatetime)
1011 1012
1012 1013 return dataOut
1013 1014
1014 1015 @MPDecorator
1015 1016 class SpectraCutPlot_(Figure):
1016 1017
1017 1018 isConfig = None
1018 1019 __nsubplots = None
1019 1020
1020 1021 WIDTHPROF = None
1021 1022 HEIGHTPROF = None
1022 1023 PREFIX = 'spc_cut'
1023 1024
1024 1025 def __init__(self):
1025 1026 Figure.__init__(self)
1026 1027 self.isConfig = False
1027 1028 self.__nsubplots = 1
1028 1029
1029 1030 self.PLOT_CODE = POWER_CODE
1030 1031
1031 1032 self.WIDTH = 700
1032 1033 self.HEIGHT = 500
1033 1034 self.counter_imagwr = 0
1034 1035
1035 1036 def getSubplots(self):
1036 1037 ncol = 1
1037 1038 nrow = 1
1038 1039
1039 1040 return nrow, ncol
1040 1041
1041 1042 def setup(self, id, nplots, wintitle, show):
1042 1043
1043 1044 self.nplots = nplots
1044 1045
1045 1046 ncolspan = 1
1046 1047 colspan = 1
1047 1048
1048 1049 self.createFigure(id = id,
1049 1050 wintitle = wintitle,
1050 1051 widthplot = self.WIDTH,
1051 1052 heightplot = self.HEIGHT,
1052 1053 show=show)
1053 1054
1054 1055 nrow, ncol = self.getSubplots()
1055 1056
1056 1057 counter = 0
1057 1058 for y in range(nrow):
1058 1059 for x in range(ncol):
1059 1060 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
1060 1061
1061 1062 def run(self, dataOut, id, wintitle="", channelList=None,
1062 1063 xmin=None, xmax=None, ymin=None, ymax=None,
1063 1064 save=False, figpath='./', figfile=None, show=True,
1064 1065 ftp=False, wr_period=1, server=None,
1065 1066 folder=None, username=None, password=None,
1066 1067 xaxis="frequency"):
1067 1068
1068 1069 if dataOut.flagNoData:
1069 1070 return dataOut
1070 1071
1071 1072 if channelList == None:
1072 1073 channelIndexList = dataOut.channelIndexList
1073 1074 channelList = dataOut.channelList
1074 1075 else:
1075 1076 channelIndexList = []
1076 1077 for channel in channelList:
1077 1078 if channel not in dataOut.channelList:
1078 1079 raise ValueError("Channel %d is not in dataOut.channelList")
1079 1080 channelIndexList.append(dataOut.channelList.index(channel))
1080 1081
1081 1082 factor = dataOut.normFactor
1082 1083
1083 1084 y = dataOut.getHeiRange()
1084 1085
1085 1086 z = dataOut.data_spc/factor
1086 1087 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
1087 1088
1088 1089 hei_index = numpy.arange(25)*3 + 20
1089 1090
1090 1091 if xaxis == "frequency":
1091 1092 x = dataOut.getFreqRange()/1000.
1092 1093 zdB = 10*numpy.log10(z[0,:,hei_index])
1093 1094 xlabel = "Frequency (kHz)"
1094 1095 ylabel = "Power (dB)"
1095 1096
1096 1097 elif xaxis == "time":
1097 1098 x = dataOut.getAcfRange()
1098 1099 zdB = z[0,:,hei_index]
1099 1100 xlabel = "Time (ms)"
1100 1101 ylabel = "ACF"
1101 1102
1102 1103 else:
1103 1104 x = dataOut.getVelRange()
1104 1105 zdB = 10*numpy.log10(z[0,:,hei_index])
1105 1106 xlabel = "Velocity (m/s)"
1106 1107 ylabel = "Power (dB)"
1107 1108
1108 1109 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
1109 1110 title = wintitle + " Range Cuts %s" %(thisDatetime.strftime("%d-%b-%Y"))
1110 1111
1111 1112 if not self.isConfig:
1112 1113
1113 1114 nplots = 1
1114 1115
1115 1116 self.setup(id=id,
1116 1117 nplots=nplots,
1117 1118 wintitle=wintitle,
1118 1119 show=show)
1119 1120
1120 1121 if xmin == None: xmin = numpy.nanmin(x)*0.9
1121 1122 if xmax == None: xmax = numpy.nanmax(x)*1.1
1122 1123 if ymin == None: ymin = numpy.nanmin(zdB)
1123 1124 if ymax == None: ymax = numpy.nanmax(zdB)
1124 1125
1125 1126 self.isConfig = True
1126 1127
1127 1128 self.setWinTitle(title)
1128 1129
1129 1130 title = "Spectra Cuts: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
1130 1131 axes = self.axesList[0]
1131 1132
1132 1133 legendlabels = ["Range = %dKm" %y[i] for i in hei_index]
1133 1134
1134 1135 axes.pmultilineyaxis( x, zdB,
1135 1136 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
1136 1137 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels,
1137 1138 ytick_visible=True, nxticks=5,
1138 1139 grid='x')
1139 1140
1140 1141 self.draw()
1141 1142
1142 1143 self.save(figpath=figpath,
1143 1144 figfile=figfile,
1144 1145 save=save,
1145 1146 ftp=ftp,
1146 1147 wr_period=wr_period,
1147 1148 thisDatetime=thisDatetime)
1148 1149
1149 1150 return dataOut
1150 1151
1151 1152 @MPDecorator
1152 1153 class Noise_(Figure):
1153 1154
1154 1155 isConfig = None
1155 1156 __nsubplots = None
1156 1157
1157 1158 PREFIX = 'noise'
1158 1159
1159 1160
1160 1161 def __init__(self):
1161 1162 Figure.__init__(self)
1162 1163 self.timerange = 24*60*60
1163 1164 self.isConfig = False
1164 1165 self.__nsubplots = 1
1165 1166 self.counter_imagwr = 0
1166 1167 self.WIDTH = 800
1167 1168 self.HEIGHT = 400
1168 1169 self.WIDTHPROF = 120
1169 1170 self.HEIGHTPROF = 0
1170 1171 self.xdata = None
1171 1172 self.ydata = None
1172 1173
1173 1174 self.PLOT_CODE = NOISE_CODE
1174 1175
1175 1176 self.FTP_WEI = None
1176 1177 self.EXP_CODE = None
1177 1178 self.SUB_EXP_CODE = None
1178 1179 self.PLOT_POS = None
1179 1180 self.figfile = None
1180 1181
1181 1182 self.xmin = None
1182 1183 self.xmax = None
1183 1184
1184 1185 def getSubplots(self):
1185 1186
1186 1187 ncol = 1
1187 1188 nrow = 1
1188 1189
1189 1190 return nrow, ncol
1190 1191
1191 1192 def openfile(self, filename):
1192 1193 dirname = os.path.dirname(filename)
1193 1194
1194 1195 if not os.path.exists(dirname):
1195 1196 os.mkdir(dirname)
1196 1197
1197 1198 f = open(filename,'w+')
1198 1199 f.write('\n\n')
1199 1200 f.write('JICAMARCA RADIO OBSERVATORY - Noise \n')
1200 1201 f.write('DD MM YYYY HH MM SS Channel0 Channel1 Channel2 Channel3\n\n' )
1201 1202 f.close()
1202 1203
1203 1204 def save_data(self, filename_phase, data, data_datetime):
1204 1205
1205 1206 f=open(filename_phase,'a')
1206 1207
1207 1208 timetuple_data = data_datetime.timetuple()
1208 1209 day = str(timetuple_data.tm_mday)
1209 1210 month = str(timetuple_data.tm_mon)
1210 1211 year = str(timetuple_data.tm_year)
1211 1212 hour = str(timetuple_data.tm_hour)
1212 1213 minute = str(timetuple_data.tm_min)
1213 1214 second = str(timetuple_data.tm_sec)
1214 1215
1215 1216 data_msg = ''
1216 1217 for i in range(len(data)):
1217 1218 data_msg += str(data[i]) + ' '
1218 1219
1219 1220 f.write(day+' '+month+' '+year+' '+hour+' '+minute+' '+second+' ' + data_msg + '\n')
1220 1221 f.close()
1221 1222
1222 1223
1223 1224 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
1224 1225
1225 1226 self.__showprofile = showprofile
1226 1227 self.nplots = nplots
1227 1228
1228 1229 ncolspan = 7
1229 1230 colspan = 6
1230 1231 self.__nsubplots = 2
1231 1232
1232 1233 self.createFigure(id = id,
1233 1234 wintitle = wintitle,
1234 1235 widthplot = self.WIDTH+self.WIDTHPROF,
1235 1236 heightplot = self.HEIGHT+self.HEIGHTPROF,
1236 1237 show=show)
1237 1238
1238 1239 nrow, ncol = self.getSubplots()
1239 1240
1240 1241 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
1241 1242
1242 1243
1243 1244 def run(self, dataOut, id, wintitle="", channelList=None, showprofile='True',
1244 1245 xmin=None, xmax=None, ymin=None, ymax=None,
1245 1246 timerange=None,
1246 1247 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
1247 1248 server=None, folder=None, username=None, password=None,
1248 1249 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1249 1250
1250 1251 if dataOut.flagNoData:
1251 1252 return dataOut
1252 1253
1253 1254 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
1254 1255 return
1255 1256
1256 1257 if channelList == None:
1257 1258 channelIndexList = dataOut.channelIndexList
1258 1259 channelList = dataOut.channelList
1259 1260 else:
1260 1261 channelIndexList = []
1261 1262 for channel in channelList:
1262 1263 if channel not in dataOut.channelList:
1263 1264 raise ValueError("Channel %d is not in dataOut.channelList")
1264 1265 channelIndexList.append(dataOut.channelList.index(channel))
1265 1266
1266 1267 x = dataOut.getTimeRange()
1267 1268 #y = dataOut.getHeiRange()
1268 1269 factor = dataOut.normFactor
1269 1270 noise = dataOut.noise[channelIndexList]/factor
1270 1271 noisedB = 10*numpy.log10(noise)
1271 1272
1272 1273 thisDatetime = dataOut.datatime
1273 1274
1274 1275 title = wintitle + " Noise" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
1275 1276 xlabel = ""
1276 1277 ylabel = "Intensity (dB)"
1277 1278 update_figfile = False
1278 1279
1279 1280 if not self.isConfig:
1280 1281
1281 1282 nplots = 1
1282 1283
1283 1284 self.setup(id=id,
1284 1285 nplots=nplots,
1285 1286 wintitle=wintitle,
1286 1287 showprofile=showprofile,
1287 1288 show=show)
1288 1289
1289 1290 if timerange != None:
1290 1291 self.timerange = timerange
1291 1292
1292 1293 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
1293 1294
1294 1295 if ymin == None: ymin = numpy.floor(numpy.nanmin(noisedB)) - 10.0
1295 1296 if ymax == None: ymax = numpy.nanmax(noisedB) + 10.0
1296 1297
1297 1298 self.FTP_WEI = ftp_wei
1298 1299 self.EXP_CODE = exp_code
1299 1300 self.SUB_EXP_CODE = sub_exp_code
1300 1301 self.PLOT_POS = plot_pos
1301 1302
1302 1303
1303 1304 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1304 1305 self.isConfig = True
1305 1306 self.figfile = figfile
1306 1307 self.xdata = numpy.array([])
1307 1308 self.ydata = numpy.array([])
1308 1309
1309 1310 update_figfile = True
1310 1311
1311 1312 #open file beacon phase
1312 1313 path = '%s%03d' %(self.PREFIX, self.id)
1313 1314 noise_file = os.path.join(path,'%s.txt'%self.name)
1314 1315 self.filename_noise = os.path.join(figpath,noise_file)
1315 1316
1316 1317 self.setWinTitle(title)
1317 1318
1318 1319 title = "Noise %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1319 1320
1320 1321 legendlabels = ["channel %d"%(idchannel) for idchannel in channelList]
1321 1322 axes = self.axesList[0]
1322 1323
1323 1324 self.xdata = numpy.hstack((self.xdata, x[0:1]))
1324 1325
1325 1326 if len(self.ydata)==0:
1326 1327 self.ydata = noisedB.reshape(-1,1)
1327 1328 else:
1328 1329 self.ydata = numpy.hstack((self.ydata, noisedB.reshape(-1,1)))
1329 1330
1330 1331
1331 1332 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
1332 1333 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
1333 1334 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
1334 1335 XAxisAsTime=True, grid='both'
1335 1336 )
1336 1337
1337 1338 self.draw()
1338 1339
1339 1340 if dataOut.ltctime >= self.xmax:
1340 1341 self.counter_imagwr = wr_period
1341 1342 self.isConfig = False
1342 1343 update_figfile = True
1343 1344
1344 1345 self.save(figpath=figpath,
1345 1346 figfile=figfile,
1346 1347 save=save,
1347 1348 ftp=ftp,
1348 1349 wr_period=wr_period,
1349 1350 thisDatetime=thisDatetime,
1350 1351 update_figfile=update_figfile)
1351 1352
1352 1353 #store data beacon phase
1353 1354 if save:
1354 1355 self.save_data(self.filename_noise, noisedB, thisDatetime)
1355 1356
1356 1357 return dataOut
1357 1358
1358 1359 @MPDecorator
1359 1360 class BeaconPhase_(Figure):
1360 1361
1361 1362 __isConfig = None
1362 1363 __nsubplots = None
1363 1364
1364 1365 PREFIX = 'beacon_phase'
1365 1366
1366 1367 def __init__(self):
1367 1368 Figure.__init__(self)
1368 1369 self.timerange = 24*60*60
1369 1370 self.isConfig = False
1370 1371 self.__nsubplots = 1
1371 1372 self.counter_imagwr = 0
1372 1373 self.WIDTH = 800
1373 1374 self.HEIGHT = 400
1374 1375 self.WIDTHPROF = 120
1375 1376 self.HEIGHTPROF = 0
1376 1377 self.xdata = None
1377 1378 self.ydata = None
1378 1379
1379 1380 self.PLOT_CODE = BEACON_CODE
1380 1381
1381 1382 self.FTP_WEI = None
1382 1383 self.EXP_CODE = None
1383 1384 self.SUB_EXP_CODE = None
1384 1385 self.PLOT_POS = None
1385 1386
1386 1387 self.filename_phase = None
1387 1388
1388 1389 self.figfile = None
1389 1390
1390 1391 self.xmin = None
1391 1392 self.xmax = None
1392 1393
1393 1394 def getSubplots(self):
1394 1395
1395 1396 ncol = 1
1396 1397 nrow = 1
1397 1398
1398 1399 return nrow, ncol
1399 1400
1400 1401 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
1401 1402
1402 1403 self.__showprofile = showprofile
1403 1404 self.nplots = nplots
1404 1405
1405 1406 ncolspan = 7
1406 1407 colspan = 6
1407 1408 self.__nsubplots = 2
1408 1409
1409 1410 self.createFigure(id = id,
1410 1411 wintitle = wintitle,
1411 1412 widthplot = self.WIDTH+self.WIDTHPROF,
1412 1413 heightplot = self.HEIGHT+self.HEIGHTPROF,
1413 1414 show=show)
1414 1415
1415 1416 nrow, ncol = self.getSubplots()
1416 1417
1417 1418 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
1418 1419
1419 1420 def save_phase(self, filename_phase):
1420 1421 f = open(filename_phase,'w+')
1421 1422 f.write('\n\n')
1422 1423 f.write('JICAMARCA RADIO OBSERVATORY - Beacon Phase \n')
1423 1424 f.write('DD MM YYYY HH MM SS pair(2,0) pair(2,1) pair(2,3) pair(2,4)\n\n' )
1424 1425 f.close()
1425 1426
1426 1427 def save_data(self, filename_phase, data, data_datetime):
1427 1428 f=open(filename_phase,'a')
1428 1429 timetuple_data = data_datetime.timetuple()
1429 1430 day = str(timetuple_data.tm_mday)
1430 1431 month = str(timetuple_data.tm_mon)
1431 1432 year = str(timetuple_data.tm_year)
1432 1433 hour = str(timetuple_data.tm_hour)
1433 1434 minute = str(timetuple_data.tm_min)
1434 1435 second = str(timetuple_data.tm_sec)
1435 1436 f.write(day+' '+month+' '+year+' '+hour+' '+minute+' '+second+' '+str(data[0])+' '+str(data[1])+' '+str(data[2])+' '+str(data[3])+'\n')
1436 1437 f.close()
1437 1438
1438 1439
1439 1440 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
1440 1441 xmin=None, xmax=None, ymin=None, ymax=None, hmin=None, hmax=None,
1441 1442 timerange=None,
1442 1443 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
1443 1444 server=None, folder=None, username=None, password=None,
1444 1445 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1445 1446
1446 1447 if dataOut.flagNoData:
1447 1448 return dataOut
1448 1449
1449 1450 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
1450 1451 return
1451 1452
1452 1453 if pairsList == None:
1453 1454 pairsIndexList = dataOut.pairsIndexList[:10]
1454 1455 else:
1455 1456 pairsIndexList = []
1456 1457 for pair in pairsList:
1457 1458 if pair not in dataOut.pairsList:
1458 1459 raise ValueError("Pair %s is not in dataOut.pairsList" %(pair))
1459 1460 pairsIndexList.append(dataOut.pairsList.index(pair))
1460 1461
1461 1462 if pairsIndexList == []:
1462 1463 return
1463 1464
1464 1465 # if len(pairsIndexList) > 4:
1465 1466 # pairsIndexList = pairsIndexList[0:4]
1466 1467
1467 1468 hmin_index = None
1468 1469 hmax_index = None
1469 1470
1470 1471 if hmin != None and hmax != None:
1471 1472 indexes = numpy.arange(dataOut.nHeights)
1472 1473 hmin_list = indexes[dataOut.heightList >= hmin]
1473 1474 hmax_list = indexes[dataOut.heightList <= hmax]
1474 1475
1475 1476 if hmin_list.any():
1476 1477 hmin_index = hmin_list[0]
1477 1478
1478 1479 if hmax_list.any():
1479 1480 hmax_index = hmax_list[-1]+1
1480 1481
1481 1482 x = dataOut.getTimeRange()
1482 1483 #y = dataOut.getHeiRange()
1483 1484
1484 1485
1485 1486 thisDatetime = dataOut.datatime
1486 1487
1487 1488 title = wintitle + " Signal Phase" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
1488 1489 xlabel = "Local Time"
1489 1490 ylabel = "Phase (degrees)"
1490 1491
1491 1492 update_figfile = False
1492 1493
1493 1494 nplots = len(pairsIndexList)
1494 1495 #phase = numpy.zeros((len(pairsIndexList),len(dataOut.beacon_heiIndexList)))
1495 1496 phase_beacon = numpy.zeros(len(pairsIndexList))
1496 1497 for i in range(nplots):
1497 1498 pair = dataOut.pairsList[pairsIndexList[i]]
1498 1499 ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i], :, hmin_index:hmax_index], axis=0)
1499 1500 powa = numpy.average(dataOut.data_spc[pair[0], :, hmin_index:hmax_index], axis=0)
1500 1501 powb = numpy.average(dataOut.data_spc[pair[1], :, hmin_index:hmax_index], axis=0)
1501 1502 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
1502 1503 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
1503 1504
1504 1505 if dataOut.beacon_heiIndexList:
1505 1506 phase_beacon[i] = numpy.average(phase[dataOut.beacon_heiIndexList])
1506 1507 else:
1507 1508 phase_beacon[i] = numpy.average(phase)
1508 1509
1509 1510 if not self.isConfig:
1510 1511
1511 1512 nplots = len(pairsIndexList)
1512 1513
1513 1514 self.setup(id=id,
1514 1515 nplots=nplots,
1515 1516 wintitle=wintitle,
1516 1517 showprofile=showprofile,
1517 1518 show=show)
1518 1519
1519 1520 if timerange != None:
1520 1521 self.timerange = timerange
1521 1522
1522 1523 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
1523 1524
1524 1525 if ymin == None: ymin = 0
1525 1526 if ymax == None: ymax = 360
1526 1527
1527 1528 self.FTP_WEI = ftp_wei
1528 1529 self.EXP_CODE = exp_code
1529 1530 self.SUB_EXP_CODE = sub_exp_code
1530 1531 self.PLOT_POS = plot_pos
1531 1532
1532 1533 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1533 1534 self.isConfig = True
1534 1535 self.figfile = figfile
1535 1536 self.xdata = numpy.array([])
1536 1537 self.ydata = numpy.array([])
1537 1538
1538 1539 update_figfile = True
1539 1540
1540 1541 #open file beacon phase
1541 1542 path = '%s%03d' %(self.PREFIX, self.id)
1542 1543 beacon_file = os.path.join(path,'%s.txt'%self.name)
1543 1544 self.filename_phase = os.path.join(figpath,beacon_file)
1544 1545 #self.save_phase(self.filename_phase)
1545 1546
1546 1547
1547 1548 #store data beacon phase
1548 1549 #self.save_data(self.filename_phase, phase_beacon, thisDatetime)
1549 1550
1550 1551 self.setWinTitle(title)
1551 1552
1552 1553
1553 1554 title = "Phase Plot %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1554 1555
1555 1556 legendlabels = ["Pair (%d,%d)"%(pair[0], pair[1]) for pair in dataOut.pairsList]
1556 1557
1557 1558 axes = self.axesList[0]
1558 1559
1559 1560 self.xdata = numpy.hstack((self.xdata, x[0:1]))
1560 1561
1561 1562 if len(self.ydata)==0:
1562 1563 self.ydata = phase_beacon.reshape(-1,1)
1563 1564 else:
1564 1565 self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1)))
1565 1566
1566 1567
1567 1568 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
1568 1569 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
1569 1570 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
1570 1571 XAxisAsTime=True, grid='both'
1571 1572 )
1572 1573
1573 1574 self.draw()
1574 1575
1575 1576 if dataOut.ltctime >= self.xmax:
1576 1577 self.counter_imagwr = wr_period
1577 1578 self.isConfig = False
1578 1579 update_figfile = True
1579 1580
1580 1581 self.save(figpath=figpath,
1581 1582 figfile=figfile,
1582 1583 save=save,
1583 1584 ftp=ftp,
1584 1585 wr_period=wr_period,
1585 1586 thisDatetime=thisDatetime,
1586 1587 update_figfile=update_figfile)
1587 1588
1588 1589 return dataOut No newline at end of file
@@ -1,500 +1,470
1 1 import os
2 2 import sys
3 3 import datetime
4 4 import numpy
5 import matplotlib
6
7 if 'BACKEND' in os.environ:
8 matplotlib.use(os.environ['BACKEND'])
9 elif 'linux' in sys.platform:
10 matplotlib.use("TkAgg")
11 elif 'darwin' in sys.platform:
12 matplotlib.use('TkAgg')
13 else:
14 from schainpy.utils import log
15 log.warning('Using default Backend="Agg"', 'INFO')
16 matplotlib.use('Agg')
17 # Qt4Agg', 'GTK', 'GTKAgg', 'ps', 'agg', 'cairo', 'MacOSX', 'GTKCairo', 'WXAgg', 'template', 'TkAgg', 'GTK3Cairo', 'GTK3Agg', 'svg', 'WebAgg', 'CocoaAgg', 'emf', 'gdk', 'WX'
18 import matplotlib.pyplot
19
20 from mpl_toolkits.axes_grid1 import make_axes_locatable
21 from matplotlib.ticker import FuncFormatter, LinearLocator
22
23 ###########################################
24 # Actualizacion de las funciones del driver
25 ###########################################
26
27 # create jro colormap
28
29 jet_values = matplotlib.pyplot.get_cmap("jet", 100)(numpy.arange(100))[10:90]
30 blu_values = matplotlib.pyplot.get_cmap(
31 "seismic_r", 20)(numpy.arange(20))[10:15]
32 ncmap = matplotlib.colors.LinearSegmentedColormap.from_list(
33 "jro", numpy.vstack((blu_values, jet_values)))
34 matplotlib.pyplot.register_cmap(cmap=ncmap)
35
5 from .jroplot_base import matplotlib, make_axes_locatable, FuncFormatter, LinearLocator
36 6
37 7 def createFigure(id, wintitle, width, height, facecolor="w", show=True, dpi=80):
38 8
39 9 matplotlib.pyplot.ioff()
40 10
41 11 fig = matplotlib.pyplot.figure(num=id, facecolor=facecolor, figsize=(
42 12 1.0 * width / dpi, 1.0 * height / dpi))
43 13 fig.canvas.manager.set_window_title(wintitle)
44 14 # fig.canvas.manager.resize(width, height)
45 15 matplotlib.pyplot.ion()
46 16
47 17 if show:
48 18 matplotlib.pyplot.show()
49 19
50 20 return fig
51 21
52 22
53 23 def closeFigure(show=False, fig=None):
54 24
55 25 # matplotlib.pyplot.ioff()
56 26 # matplotlib.pyplot.pause(0)
57 27
58 28 if show:
59 29 matplotlib.pyplot.show()
60 30
61 31 if fig != None:
62 32 matplotlib.pyplot.close(fig)
63 33 # matplotlib.pyplot.pause(0)
64 34 # matplotlib.pyplot.ion()
65 35
66 36 return
67 37
68 38 matplotlib.pyplot.close("all")
69 39 # matplotlib.pyplot.pause(0)
70 40 # matplotlib.pyplot.ion()
71 41
72 42 return
73 43
74 44
75 45 def saveFigure(fig, filename):
76 46
77 47 # matplotlib.pyplot.ioff()
78 48 fig.savefig(filename, dpi=matplotlib.pyplot.gcf().dpi)
79 49 # matplotlib.pyplot.ion()
80 50
81 51
82 52 def clearFigure(fig):
83 53
84 54 fig.clf()
85 55
86 56
87 57 def setWinTitle(fig, title):
88 58
89 59 fig.canvas.manager.set_window_title(title)
90 60
91 61
92 62 def setTitle(fig, title):
93 63
94 64 fig.suptitle(title)
95 65
96 66
97 67 def createAxes(fig, nrow, ncol, xpos, ypos, colspan, rowspan, polar=False):
98 68
99 69 matplotlib.pyplot.ioff()
100 70 matplotlib.pyplot.figure(fig.number)
101 71 axes = matplotlib.pyplot.subplot2grid((nrow, ncol),
102 72 (xpos, ypos),
103 73 colspan=colspan,
104 74 rowspan=rowspan,
105 75 polar=polar)
106 76
107 77 matplotlib.pyplot.ion()
108 78 return axes
109 79
110 80
111 81 def setAxesText(ax, text):
112 82
113 83 ax.annotate(text,
114 84 xy=(.1, .99),
115 85 xycoords='figure fraction',
116 86 horizontalalignment='left',
117 87 verticalalignment='top',
118 88 fontsize=10)
119 89
120 90
121 91 def printLabels(ax, xlabel, ylabel, title):
122 92
123 93 ax.set_xlabel(xlabel, size=11)
124 94 ax.set_ylabel(ylabel, size=11)
125 95 ax.set_title(title, size=8)
126 96
127 97
128 98 def createPline(ax, x, y, xmin, xmax, ymin, ymax, xlabel='', ylabel='', title='',
129 99 ticksize=9, xtick_visible=True, ytick_visible=True,
130 100 nxticks=4, nyticks=10,
131 101 grid=None, color='blue'):
132 102 """
133 103
134 104 Input:
135 105 grid : None, 'both', 'x', 'y'
136 106 """
137 107
138 108 matplotlib.pyplot.ioff()
139 109
140 110 ax.set_xlim([xmin, xmax])
141 111 ax.set_ylim([ymin, ymax])
142 112
143 113 printLabels(ax, xlabel, ylabel, title)
144 114
145 115 ######################################################
146 116 if (xmax - xmin) <= 1:
147 117 xtickspos = numpy.linspace(xmin, xmax, nxticks)
148 118 xtickspos = numpy.array([float("%.1f" % i) for i in xtickspos])
149 119 ax.set_xticks(xtickspos)
150 120 else:
151 121 xtickspos = numpy.arange(nxticks) * \
152 122 int((xmax - xmin) / (nxticks)) + int(xmin)
153 123 # xtickspos = numpy.arange(nxticks)*float(xmax-xmin)/float(nxticks) + int(xmin)
154 124 ax.set_xticks(xtickspos)
155 125
156 126 for tick in ax.get_xticklabels():
157 127 tick.set_visible(xtick_visible)
158 128
159 129 for tick in ax.xaxis.get_major_ticks():
160 130 tick.label.set_fontsize(ticksize)
161 131
162 132 ######################################################
163 133 for tick in ax.get_yticklabels():
164 134 tick.set_visible(ytick_visible)
165 135
166 136 for tick in ax.yaxis.get_major_ticks():
167 137 tick.label.set_fontsize(ticksize)
168 138
169 139 ax.plot(x, y, color=color)
170 140 iplot = ax.lines[-1]
171 141
172 142 ######################################################
173 143 if '0.' in matplotlib.__version__[0:2]:
174 144 print("The matplotlib version has to be updated to 1.1 or newer")
175 145 return iplot
176 146
177 147 if '1.0.' in matplotlib.__version__[0:4]:
178 148 print("The matplotlib version has to be updated to 1.1 or newer")
179 149 return iplot
180 150
181 151 if grid != None:
182 152 ax.grid(b=True, which='major', axis=grid)
183 153
184 154 matplotlib.pyplot.tight_layout()
185 155
186 156 matplotlib.pyplot.ion()
187 157
188 158 return iplot
189 159
190 160
191 161 def set_linedata(ax, x, y, idline):
192 162
193 163 ax.lines[idline].set_data(x, y)
194 164
195 165
196 166 def pline(iplot, x, y, xlabel='', ylabel='', title=''):
197 167
198 168 ax = iplot.axes
199 169
200 170 printLabels(ax, xlabel, ylabel, title)
201 171
202 172 set_linedata(ax, x, y, idline=0)
203 173
204 174
205 175 def addpline(ax, x, y, color, linestyle, lw):
206 176
207 177 ax.plot(x, y, color=color, linestyle=linestyle, lw=lw)
208 178
209 179
210 180 def createPcolor(ax, x, y, z, xmin, xmax, ymin, ymax, zmin, zmax,
211 181 xlabel='', ylabel='', title='', ticksize=9,
212 182 colormap='jet', cblabel='', cbsize="5%",
213 183 XAxisAsTime=False):
214 184
215 185 matplotlib.pyplot.ioff()
216 186
217 187 divider = make_axes_locatable(ax)
218 188 ax_cb = divider.new_horizontal(size=cbsize, pad=0.05)
219 189 fig = ax.get_figure()
220 190 fig.add_axes(ax_cb)
221 191
222 192 ax.set_xlim([xmin, xmax])
223 193 ax.set_ylim([ymin, ymax])
224 194
225 195 printLabels(ax, xlabel, ylabel, title)
226 196
227 197 z = numpy.ma.masked_invalid(z)
228 198 cmap = matplotlib.pyplot.get_cmap(colormap)
229 199 cmap.set_bad('white', 1.)
230 200 imesh = ax.pcolormesh(x, y, z.T, vmin=zmin, vmax=zmax, cmap=cmap)
231 201 cb = matplotlib.pyplot.colorbar(imesh, cax=ax_cb)
232 202 cb.set_label(cblabel)
233 203
234 204 # for tl in ax_cb.get_yticklabels():
235 205 # tl.set_visible(True)
236 206
237 207 for tick in ax.yaxis.get_major_ticks():
238 208 tick.label.set_fontsize(ticksize)
239 209
240 210 for tick in ax.xaxis.get_major_ticks():
241 211 tick.label.set_fontsize(ticksize)
242 212
243 213 for tick in cb.ax.get_yticklabels():
244 214 tick.set_fontsize(ticksize)
245 215
246 216 ax_cb.yaxis.tick_right()
247 217
248 218 if '0.' in matplotlib.__version__[0:2]:
249 219 print("The matplotlib version has to be updated to 1.1 or newer")
250 220 return imesh
251 221
252 222 if '1.0.' in matplotlib.__version__[0:4]:
253 223 print("The matplotlib version has to be updated to 1.1 or newer")
254 224 return imesh
255 225
256 226 matplotlib.pyplot.tight_layout()
257 227
258 228 if XAxisAsTime:
259 229
260 230 def func(x, pos): return ('%s') % (
261 231 datetime.datetime.utcfromtimestamp(x).strftime("%H:%M:%S"))
262 232 ax.xaxis.set_major_formatter(FuncFormatter(func))
263 233 ax.xaxis.set_major_locator(LinearLocator(7))
264 234
265 235 matplotlib.pyplot.ion()
266 236 return imesh
267 237
268 238
269 239 def pcolor(imesh, z, xlabel='', ylabel='', title=''):
270 240
271 241 z = z.T
272 242 ax = imesh.axes
273 243 printLabels(ax, xlabel, ylabel, title)
274 244 imesh.set_array(z.ravel())
275 245
276 246
277 247 def addpcolor(ax, x, y, z, zmin, zmax, xlabel='', ylabel='', title='', colormap='jet'):
278 248
279 249 printLabels(ax, xlabel, ylabel, title)
280 250
281 251 ax.pcolormesh(x, y, z.T, vmin=zmin, vmax=zmax,
282 252 cmap=matplotlib.pyplot.get_cmap(colormap))
283 253
284 254
285 255 def addpcolorbuffer(ax, x, y, z, zmin, zmax, xlabel='', ylabel='', title='', colormap='jet'):
286 256
287 257 printLabels(ax, xlabel, ylabel, title)
288 258
289 259 ax.collections.remove(ax.collections[0])
290 260
291 261 z = numpy.ma.masked_invalid(z)
292 262
293 263 cmap = matplotlib.pyplot.get_cmap(colormap)
294 264 cmap.set_bad('white', 1.)
295 265
296 266 ax.pcolormesh(x, y, z.T, vmin=zmin, vmax=zmax, cmap=cmap)
297 267
298 268
299 269 def createPmultiline(ax, x, y, xmin, xmax, ymin, ymax, xlabel='', ylabel='', title='', legendlabels=None,
300 270 ticksize=9, xtick_visible=True, ytick_visible=True,
301 271 nxticks=4, nyticks=10,
302 272 grid=None):
303 273 """
304 274
305 275 Input:
306 276 grid : None, 'both', 'x', 'y'
307 277 """
308 278
309 279 matplotlib.pyplot.ioff()
310 280
311 281 lines = ax.plot(x.T, y)
312 282 leg = ax.legend(lines, legendlabels, loc='upper right')
313 283 leg.get_frame().set_alpha(0.5)
314 284 ax.set_xlim([xmin, xmax])
315 285 ax.set_ylim([ymin, ymax])
316 286 printLabels(ax, xlabel, ylabel, title)
317 287
318 288 xtickspos = numpy.arange(nxticks) * \
319 289 int((xmax - xmin) / (nxticks)) + int(xmin)
320 290 ax.set_xticks(xtickspos)
321 291
322 292 for tick in ax.get_xticklabels():
323 293 tick.set_visible(xtick_visible)
324 294
325 295 for tick in ax.xaxis.get_major_ticks():
326 296 tick.label.set_fontsize(ticksize)
327 297
328 298 for tick in ax.get_yticklabels():
329 299 tick.set_visible(ytick_visible)
330 300
331 301 for tick in ax.yaxis.get_major_ticks():
332 302 tick.label.set_fontsize(ticksize)
333 303
334 304 iplot = ax.lines[-1]
335 305
336 306 if '0.' in matplotlib.__version__[0:2]:
337 307 print("The matplotlib version has to be updated to 1.1 or newer")
338 308 return iplot
339 309
340 310 if '1.0.' in matplotlib.__version__[0:4]:
341 311 print("The matplotlib version has to be updated to 1.1 or newer")
342 312 return iplot
343 313
344 314 if grid != None:
345 315 ax.grid(b=True, which='major', axis=grid)
346 316
347 317 matplotlib.pyplot.tight_layout()
348 318
349 319 matplotlib.pyplot.ion()
350 320
351 321 return iplot
352 322
353 323
354 324 def pmultiline(iplot, x, y, xlabel='', ylabel='', title=''):
355 325
356 326 ax = iplot.axes
357 327
358 328 printLabels(ax, xlabel, ylabel, title)
359 329
360 330 for i in range(len(ax.lines)):
361 331 line = ax.lines[i]
362 332 line.set_data(x[i, :], y)
363 333
364 334
365 335 def createPmultilineYAxis(ax, x, y, xmin, xmax, ymin, ymax, xlabel='', ylabel='', title='', legendlabels=None,
366 336 ticksize=9, xtick_visible=True, ytick_visible=True,
367 337 nxticks=4, nyticks=10, marker='.', markersize=10, linestyle="None",
368 338 grid=None, XAxisAsTime=False):
369 339 """
370 340
371 341 Input:
372 342 grid : None, 'both', 'x', 'y'
373 343 """
374 344
375 345 matplotlib.pyplot.ioff()
376 346
377 347 # lines = ax.plot(x, y.T, marker=marker,markersize=markersize,linestyle=linestyle)
378 348 lines = ax.plot(x, y.T)
379 349 # leg = ax.legend(lines, legendlabels, loc=2, bbox_to_anchor=(1.01, 1.00), numpoints=1, handlelength=1.5, \
380 350 # handletextpad=0.5, borderpad=0.5, labelspacing=0.5, borderaxespad=0.)
381 351
382 352 leg = ax.legend(lines, legendlabels,
383 353 loc='upper right', bbox_to_anchor=(1.16, 1), borderaxespad=0)
384 354
385 355 for label in leg.get_texts():
386 356 label.set_fontsize(9)
387 357
388 358 ax.set_xlim([xmin, xmax])
389 359 ax.set_ylim([ymin, ymax])
390 360 printLabels(ax, xlabel, ylabel, title)
391 361
392 362 # xtickspos = numpy.arange(nxticks)*int((xmax-xmin)/(nxticks)) + int(xmin)
393 363 # ax.set_xticks(xtickspos)
394 364
395 365 for tick in ax.get_xticklabels():
396 366 tick.set_visible(xtick_visible)
397 367
398 368 for tick in ax.xaxis.get_major_ticks():
399 369 tick.label.set_fontsize(ticksize)
400 370
401 371 for tick in ax.get_yticklabels():
402 372 tick.set_visible(ytick_visible)
403 373
404 374 for tick in ax.yaxis.get_major_ticks():
405 375 tick.label.set_fontsize(ticksize)
406 376
407 377 iplot = ax.lines[-1]
408 378
409 379 if '0.' in matplotlib.__version__[0:2]:
410 380 print("The matplotlib version has to be updated to 1.1 or newer")
411 381 return iplot
412 382
413 383 if '1.0.' in matplotlib.__version__[0:4]:
414 384 print("The matplotlib version has to be updated to 1.1 or newer")
415 385 return iplot
416 386
417 387 if grid != None:
418 388 ax.grid(b=True, which='major', axis=grid)
419 389
420 390 matplotlib.pyplot.tight_layout()
421 391
422 392 if XAxisAsTime:
423 393
424 394 def func(x, pos): return ('%s') % (
425 395 datetime.datetime.utcfromtimestamp(x).strftime("%H:%M:%S"))
426 396 ax.xaxis.set_major_formatter(FuncFormatter(func))
427 397 ax.xaxis.set_major_locator(LinearLocator(7))
428 398
429 399 matplotlib.pyplot.ion()
430 400
431 401 return iplot
432 402
433 403
434 404 def pmultilineyaxis(iplot, x, y, xlabel='', ylabel='', title=''):
435 405
436 406 ax = iplot.axes
437 407 printLabels(ax, xlabel, ylabel, title)
438 408
439 409 for i in range(len(ax.lines)):
440 410 line = ax.lines[i]
441 411 line.set_data(x, y[i, :])
442 412
443 413
444 414 def createPolar(ax, x, y,
445 415 xlabel='', ylabel='', title='', ticksize=9,
446 416 colormap='jet', cblabel='', cbsize="5%",
447 417 XAxisAsTime=False):
448 418
449 419 matplotlib.pyplot.ioff()
450 420
451 421 ax.plot(x, y, 'bo', markersize=5)
452 422 # ax.set_rmax(90)
453 423 ax.set_ylim(0, 90)
454 424 ax.set_yticks(numpy.arange(0, 90, 20))
455 425 # ax.text(0, -110, ylabel, rotation='vertical', va ='center', ha = 'center' ,size='11')
456 426 # ax.text(0, 50, ylabel, rotation='vertical', va ='center', ha = 'left' ,size='11')
457 427 # ax.text(100, 100, 'example', ha='left', va='center', rotation='vertical')
458 428 ax.yaxis.labelpad = 40
459 429 printLabels(ax, xlabel, ylabel, title)
460 430 iplot = ax.lines[-1]
461 431
462 432 if '0.' in matplotlib.__version__[0:2]:
463 433 print("The matplotlib version has to be updated to 1.1 or newer")
464 434 return iplot
465 435
466 436 if '1.0.' in matplotlib.__version__[0:4]:
467 437 print("The matplotlib version has to be updated to 1.1 or newer")
468 438 return iplot
469 439
470 440 # if grid != None:
471 441 # ax.grid(b=True, which='major', axis=grid)
472 442
473 443 matplotlib.pyplot.tight_layout()
474 444
475 445 matplotlib.pyplot.ion()
476 446
477 447 return iplot
478 448
479 449
480 450 def polar(iplot, x, y, xlabel='', ylabel='', title=''):
481 451
482 452 ax = iplot.axes
483 453
484 454 # ax.text(0, -110, ylabel, rotation='vertical', va ='center', ha = 'center',size='11')
485 455 printLabels(ax, xlabel, ylabel, title)
486 456
487 457 set_linedata(ax, x, y, idline=0)
488 458
489 459
490 460 def draw(fig):
491 461
492 462 if type(fig) == 'int':
493 463 raise ValueError("Error drawing: Fig parameter should be a matplotlib figure object figure")
494 464
495 465 fig.canvas.draw()
496 466
497 467
498 468 def pause(interval=0.000001):
499 469
500 470 matplotlib.pyplot.pause(interval) No newline at end of file
1 NO CONTENT: modified file
General Comments 0
You need to be logged in to leave comments. Login now