##// END OF EJS Templates
schain-jc
RSS
Update plots for realtime app
Juan C. Espinoza -
r1221:94b941a5310a
parent child
Show More
Side by Side Unified
Context file:
r1221:94b941a5310a Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,1366 +1,1365
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 if self.nmodes:
594 594 return self.nmodes*timeInterval
595 595 else:
596 596 return timeInterval
597 597
598 598 def getPower(self):
599 599
600 600 factor = self.normFactor
601 601 z = self.data_spc / factor
602 602 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
603 603 avg = numpy.average(z, axis=1)
604 604
605 605 return 10 * numpy.log10(avg)
606 606
607 607 def getCoherence(self, pairsList=None, phase=False):
608 608
609 609 z = []
610 610 if pairsList is None:
611 611 pairsIndexList = self.pairsIndexList
612 612 else:
613 613 pairsIndexList = []
614 614 for pair in pairsList:
615 615 if pair not in self.pairsList:
616 616 raise ValueError("Pair %s is not in dataOut.pairsList" % (
617 617 pair))
618 618 pairsIndexList.append(self.pairsList.index(pair))
619 619 for i in range(len(pairsIndexList)):
620 620 pair = self.pairsList[pairsIndexList[i]]
621 621 ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
622 622 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
623 623 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
624 624 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
625 625 if phase:
626 626 data = numpy.arctan2(avgcoherenceComplex.imag,
627 627 avgcoherenceComplex.real) * 180 / numpy.pi
628 628 else:
629 629 data = numpy.abs(avgcoherenceComplex)
630 630
631 631 z.append(data)
632 632
633 633 return numpy.array(z)
634 634
635 635 def setValue(self, value):
636 636
637 637 print("This property should not be initialized")
638 638
639 639 return
640 640
641 641 nPairs = property(getNPairs, setValue, "I'm the 'nPairs' property.")
642 642 pairsIndexList = property(
643 643 getPairsIndexList, setValue, "I'm the 'pairsIndexList' property.")
644 644 normFactor = property(getNormFactor, setValue,
645 645 "I'm the 'getNormFactor' property.")
646 646 flag_cspc = property(getFlagCspc, setValue)
647 647 flag_dc = property(getFlagDc, setValue)
648 648 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
649 649 timeInterval = property(getTimeInterval, setValue,
650 650 "I'm the 'timeInterval' property")
651 651
652 652
653 653 class SpectraHeis(Spectra):
654 654
655 655 data_spc = None
656 656 data_cspc = None
657 657 data_dc = None
658 658 nFFTPoints = None
659 659 # nPairs = None
660 660 pairsList = None
661 661 nCohInt = None
662 662 nIncohInt = None
663 663
664 664 def __init__(self):
665 665
666 666 self.radarControllerHeaderObj = RadarControllerHeader()
667 667
668 668 self.systemHeaderObj = SystemHeader()
669 669
670 670 self.type = "SpectraHeis"
671 671
672 672 # self.dtype = None
673 673
674 674 # self.nChannels = 0
675 675
676 676 # self.nHeights = 0
677 677
678 678 self.nProfiles = None
679 679
680 680 self.heightList = None
681 681
682 682 self.channelList = None
683 683
684 684 # self.channelIndexList = None
685 685
686 686 self.flagNoData = True
687 687
688 688 self.flagDiscontinuousBlock = False
689 689
690 690 # self.nPairs = 0
691 691
692 692 self.utctime = None
693 693
694 694 self.blocksize = None
695 695
696 696 self.profileIndex = 0
697 697
698 698 self.nCohInt = 1
699 699
700 700 self.nIncohInt = 1
701 701
702 702 def getNormFactor(self):
703 703 pwcode = 1
704 704 if self.flagDecodeData:
705 705 pwcode = numpy.sum(self.code[0]**2)
706 706
707 707 normFactor = self.nIncohInt * self.nCohInt * pwcode
708 708
709 709 return normFactor
710 710
711 711 def getTimeInterval(self):
712 712
713 713 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
714 714
715 715 return timeInterval
716 716
717 717 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
718 718 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
719 719
720 720
721 721 class Fits(JROData):
722 722
723 723 heightList = None
724 724 channelList = None
725 725 flagNoData = True
726 726 flagDiscontinuousBlock = False
727 727 useLocalTime = False
728 728 utctime = None
729 729 timeZone = None
730 730 # ippSeconds = None
731 731 # timeInterval = None
732 732 nCohInt = None
733 733 nIncohInt = None
734 734 noise = None
735 735 windowOfFilter = 1
736 736 # Speed of ligth
737 737 C = 3e8
738 738 frequency = 49.92e6
739 739 realtime = False
740 740
741 741 def __init__(self):
742 742
743 743 self.type = "Fits"
744 744
745 745 self.nProfiles = None
746 746
747 747 self.heightList = None
748 748
749 749 self.channelList = None
750 750
751 751 # self.channelIndexList = None
752 752
753 753 self.flagNoData = True
754 754
755 755 self.utctime = None
756 756
757 757 self.nCohInt = 1
758 758
759 759 self.nIncohInt = 1
760 760
761 761 self.useLocalTime = True
762 762
763 763 self.profileIndex = 0
764 764
765 765 # self.utctime = None
766 766 # self.timeZone = None
767 767 # self.ltctime = None
768 768 # self.timeInterval = None
769 769 # self.header = None
770 770 # self.data_header = None
771 771 # self.data = None
772 772 # self.datatime = None
773 773 # self.flagNoData = False
774 774 # self.expName = ''
775 775 # self.nChannels = None
776 776 # self.nSamples = None
777 777 # self.dataBlocksPerFile = None
778 778 # self.comments = ''
779 779 #
780 780
781 781 def getltctime(self):
782 782
783 783 if self.useLocalTime:
784 784 return self.utctime - self.timeZone * 60
785 785
786 786 return self.utctime
787 787
788 788 def getDatatime(self):
789 789
790 790 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
791 791 return datatime
792 792
793 793 def getTimeRange(self):
794 794
795 795 datatime = []
796 796
797 797 datatime.append(self.ltctime)
798 798 datatime.append(self.ltctime + self.timeInterval)
799 799
800 800 datatime = numpy.array(datatime)
801 801
802 802 return datatime
803 803
804 804 def getHeiRange(self):
805 805
806 806 heis = self.heightList
807 807
808 808 return heis
809 809
810 810 def getNHeights(self):
811 811
812 812 return len(self.heightList)
813 813
814 814 def getNChannels(self):
815 815
816 816 return len(self.channelList)
817 817
818 818 def getChannelIndexList(self):
819 819
820 820 return list(range(self.nChannels))
821 821
822 822 def getNoise(self, type=1):
823 823
824 824 #noise = numpy.zeros(self.nChannels)
825 825
826 826 if type == 1:
827 827 noise = self.getNoisebyHildebrand()
828 828
829 829 if type == 2:
830 830 noise = self.getNoisebySort()
831 831
832 832 if type == 3:
833 833 noise = self.getNoisebyWindow()
834 834
835 835 return noise
836 836
837 837 def getTimeInterval(self):
838 838
839 839 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
840 840
841 841 return timeInterval
842 842
843 843 def get_ippSeconds(self):
844 844 '''
845 845 '''
846 846 return self.ipp_sec
847 847
848 848
849 849 datatime = property(getDatatime, "I'm the 'datatime' property")
850 850 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
851 851 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
852 852 channelIndexList = property(
853 853 getChannelIndexList, "I'm the 'channelIndexList' property.")
854 854 noise = property(getNoise, "I'm the 'nHeights' property.")
855 855
856 856 ltctime = property(getltctime, "I'm the 'ltctime' property")
857 857 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
858 858 ippSeconds = property(get_ippSeconds, '')
859 859
860 860 class Correlation(JROData):
861 861
862 862 noise = None
863 863 SNR = None
864 864 #--------------------------------------------------
865 865 mode = None
866 866 split = False
867 867 data_cf = None
868 868 lags = None
869 869 lagRange = None
870 870 pairsList = None
871 871 normFactor = None
872 872 #--------------------------------------------------
873 873 # calculateVelocity = None
874 874 nLags = None
875 875 nPairs = None
876 876 nAvg = None
877 877
878 878 def __init__(self):
879 879 '''
880 880 Constructor
881 881 '''
882 882 self.radarControllerHeaderObj = RadarControllerHeader()
883 883
884 884 self.systemHeaderObj = SystemHeader()
885 885
886 886 self.type = "Correlation"
887 887
888 888 self.data = None
889 889
890 890 self.dtype = None
891 891
892 892 self.nProfiles = None
893 893
894 894 self.heightList = None
895 895
896 896 self.channelList = None
897 897
898 898 self.flagNoData = True
899 899
900 900 self.flagDiscontinuousBlock = False
901 901
902 902 self.utctime = None
903 903
904 904 self.timeZone = None
905 905
906 906 self.dstFlag = None
907 907
908 908 self.errorCount = None
909 909
910 910 self.blocksize = None
911 911
912 912 self.flagDecodeData = False # asumo q la data no esta decodificada
913 913
914 914 self.flagDeflipData = False # asumo q la data no esta sin flip
915 915
916 916 self.pairsList = None
917 917
918 918 self.nPoints = None
919 919
920 920 def getPairsList(self):
921 921
922 922 return self.pairsList
923 923
924 924 def getNoise(self, mode=2):
925 925
926 926 indR = numpy.where(self.lagR == 0)[0][0]
927 927 indT = numpy.where(self.lagT == 0)[0][0]
928 928
929 929 jspectra0 = self.data_corr[:, :, indR, :]
930 930 jspectra = copy.copy(jspectra0)
931 931
932 932 num_chan = jspectra.shape[0]
933 933 num_hei = jspectra.shape[2]
934 934
935 935 freq_dc = jspectra.shape[1] / 2
936 936 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
937 937
938 938 if ind_vel[0] < 0:
939 939 ind_vel[list(range(0, 1))] = ind_vel[list(
940 940 range(0, 1))] + self.num_prof
941 941
942 942 if mode == 1:
943 943 jspectra[:, freq_dc, :] = (
944 944 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
945 945
946 946 if mode == 2:
947 947
948 948 vel = numpy.array([-2, -1, 1, 2])
949 949 xx = numpy.zeros([4, 4])
950 950
951 951 for fil in range(4):
952 952 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
953 953
954 954 xx_inv = numpy.linalg.inv(xx)
955 955 xx_aux = xx_inv[0, :]
956 956
957 957 for ich in range(num_chan):
958 958 yy = jspectra[ich, ind_vel, :]
959 959 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
960 960
961 961 junkid = jspectra[ich, freq_dc, :] <= 0
962 962 cjunkid = sum(junkid)
963 963
964 964 if cjunkid.any():
965 965 jspectra[ich, freq_dc, junkid.nonzero()] = (
966 966 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
967 967
968 968 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
969 969
970 970 return noise
971 971
972 972 def getTimeInterval(self):
973 973
974 974 timeInterval = self.ippSeconds * self.nCohInt * self.nProfiles
975 975
976 976 return timeInterval
977 977
978 978 def splitFunctions(self):
979 979
980 980 pairsList = self.pairsList
981 981 ccf_pairs = []
982 982 acf_pairs = []
983 983 ccf_ind = []
984 984 acf_ind = []
985 985 for l in range(len(pairsList)):
986 986 chan0 = pairsList[l][0]
987 987 chan1 = pairsList[l][1]
988 988
989 989 # Obteniendo pares de Autocorrelacion
990 990 if chan0 == chan1:
991 991 acf_pairs.append(chan0)
992 992 acf_ind.append(l)
993 993 else:
994 994 ccf_pairs.append(pairsList[l])
995 995 ccf_ind.append(l)
996 996
997 997 data_acf = self.data_cf[acf_ind]
998 998 data_ccf = self.data_cf[ccf_ind]
999 999
1000 1000 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
1001 1001
1002 1002 def getNormFactor(self):
1003 1003 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
1004 1004 acf_pairs = numpy.array(acf_pairs)
1005 1005 normFactor = numpy.zeros((self.nPairs, self.nHeights))
1006 1006
1007 1007 for p in range(self.nPairs):
1008 1008 pair = self.pairsList[p]
1009 1009
1010 1010 ch0 = pair[0]
1011 1011 ch1 = pair[1]
1012 1012
1013 1013 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
1014 1014 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
1015 1015 normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
1016 1016
1017 1017 return normFactor
1018 1018
1019 1019 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
1020 1020 normFactor = property(getNormFactor, "I'm the 'normFactor property'")
1021 1021
1022 1022
1023 1023 class Parameters(Spectra):
1024 1024
1025 1025 experimentInfo = None # Information about the experiment
1026 1026 # Information from previous data
1027 1027 inputUnit = None # Type of data to be processed
1028 1028 operation = None # Type of operation to parametrize
1029 1029 # normFactor = None #Normalization Factor
1030 1030 groupList = None # List of Pairs, Groups, etc
1031 1031 # Parameters
1032 1032 data_param = None # Parameters obtained
1033 1033 data_pre = None # Data Pre Parametrization
1034 1034 data_SNR = None # Signal to Noise Ratio
1035 1035 # heightRange = None #Heights
1036 1036 abscissaList = None # Abscissa, can be velocities, lags or time
1037 1037 # noise = None #Noise Potency
1038 1038 utctimeInit = None # Initial UTC time
1039 1039 paramInterval = None # Time interval to calculate Parameters in seconds
1040 1040 useLocalTime = True
1041 1041 # Fitting
1042 1042 data_error = None # Error of the estimation
1043 1043 constants = None
1044 1044 library = None
1045 1045 # Output signal
1046 1046 outputInterval = None # Time interval to calculate output signal in seconds
1047 1047 data_output = None # Out signal
1048 1048 nAvg = None
1049 1049 noise_estimation = None
1050 1050 GauSPC = None # Fit gaussian SPC
1051 1051
1052 1052 def __init__(self):
1053 1053 '''
1054 1054 Constructor
1055 1055 '''
1056 1056 self.radarControllerHeaderObj = RadarControllerHeader()
1057 1057
1058 1058 self.systemHeaderObj = SystemHeader()
1059 1059
1060 1060 self.type = "Parameters"
1061 1061
1062 1062 def getTimeRange1(self, interval):
1063 1063
1064 1064 datatime = []
1065 1065
1066 1066 if self.useLocalTime:
1067 1067 time1 = self.utctimeInit - self.timeZone * 60
1068 1068 else:
1069 1069 time1 = self.utctimeInit
1070 1070
1071 1071 datatime.append(time1)
1072 1072 datatime.append(time1 + interval)
1073 1073 datatime = numpy.array(datatime)
1074 1074
1075 1075 return datatime
1076 1076
1077 1077 def getTimeInterval(self):
1078 1078
1079 1079 if hasattr(self, 'timeInterval1'):
1080 1080 return self.timeInterval1
1081 1081 else:
1082 1082 return self.paramInterval
1083 1083
1084 1084 def setValue(self, value):
1085 1085
1086 1086 print("This property should not be initialized")
1087 1087
1088 1088 return
1089 1089
1090 1090 def getNoise(self):
1091 1091
1092 1092 return self.spc_noise
1093 1093
1094 1094 timeInterval = property(getTimeInterval)
1095 1095 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
1096 1096
1097 1097
1098 1098 class PlotterData(object):
1099 1099 '''
1100 1100 Object to hold data to be plotted
1101 1101 '''
1102 1102
1103 1103 MAXNUMX = 100
1104 1104 MAXNUMY = 100
1105 1105
1106 1106 def __init__(self, code, throttle_value, exp_code, buffering=True, snr=False):
1107 1107
1108 self.key = code
1108 1109 self.throttle = throttle_value
1109 1110 self.exp_code = exp_code
1110 1111 self.buffering = buffering
1111 1112 self.ready = False
1112 1113 self.localtime = False
1113 1114 self.data = {}
1114 1115 self.meta = {}
1115 1116 self.__times = []
1116 1117 self.__heights = []
1117 1118
1118 1119 if 'snr' in code:
1119 1120 self.plottypes = ['snr']
1120 1121 elif code == 'spc':
1121 1122 self.plottypes = ['spc', 'noise', 'rti']
1122 1123 elif code == 'rti':
1123 1124 self.plottypes = ['noise', 'rti']
1124 1125 else:
1125 1126 self.plottypes = [code]
1126 1127
1127 1128 if 'snr' not in self.plottypes and snr:
1128 1129 self.plottypes.append('snr')
1129 1130
1130 1131 for plot in self.plottypes:
1131 1132 self.data[plot] = {}
1132 1133
1133 1134 def __str__(self):
1134 1135 dum = ['{}{}'.format(key, self.shape(key)) for key in self.data]
1135 1136 return 'Data[{}][{}]'.format(';'.join(dum), len(self.__times))
1136 1137
1137 1138 def __len__(self):
1138 1139 return len(self.__times)
1139 1140
1140 1141 def __getitem__(self, key):
1141 1142
1142 1143 if key not in self.data:
1143 1144 raise KeyError(log.error('Missing key: {}'.format(key)))
1144 1145 if 'spc' in key or not self.buffering:
1145 1146 ret = self.data[key]
1146 1147 elif 'scope' in key:
1147 1148 ret = numpy.array(self.data[key][float(self.tm)])
1148 1149 else:
1149 1150 ret = numpy.array([self.data[key][x] for x in self.times])
1150 1151 if ret.ndim > 1:
1151 1152 ret = numpy.swapaxes(ret, 0, 1)
1152 1153 return ret
1153 1154
1154 1155 def __contains__(self, key):
1155 1156 return key in self.data
1156 1157
1157 1158 def setup(self):
1158 1159 '''
1159 1160 Configure object
1160 1161 '''
1161 1162
1162 1163 self.type = ''
1163 1164 self.ready = False
1164 1165 self.data = {}
1165 1166 self.__times = []
1166 1167 self.__heights = []
1167 1168 self.__all_heights = set()
1168 1169 for plot in self.plottypes:
1169 1170 if 'snr' in plot:
1170 1171 plot = 'snr'
1171 1172 elif 'spc_moments' == plot:
1172 1173 plot = 'moments'
1173 1174 self.data[plot] = {}
1174 1175
1175 1176 if 'spc' in self.data or 'rti' in self.data or 'cspc' in self.data or 'moments' in self.data:
1176 1177 self.data['noise'] = {}
1177 1178 self.data['rti'] = {}
1178 1179 if 'noise' not in self.plottypes:
1179 1180 self.plottypes.append('noise')
1180 1181 if 'rti' not in self.plottypes:
1181 1182 self.plottypes.append('rti')
1182 1183
1183 1184 def shape(self, key):
1184 1185 '''
1185 1186 Get the shape of the one-element data for the given key
1186 1187 '''
1187 1188
1188 1189 if len(self.data[key]):
1189 1190 if 'spc' in key or not self.buffering:
1190 1191 return self.data[key].shape
1191 1192 return self.data[key][self.__times[0]].shape
1192 1193 return (0,)
1193 1194
1194 1195 def update(self, dataOut, tm):
1195 1196 '''
1196 1197 Update data object with new dataOut
1197 1198 '''
1198 1199
1199 1200 if tm in self.__times:
1200 1201 return
1201 1202 self.profileIndex = dataOut.profileIndex
1202 1203 self.tm = tm
1203 1204 self.type = dataOut.type
1204 1205 self.parameters = getattr(dataOut, 'parameters', [])
1205 if hasattr(dataOut, 'pairsList'):
1206 self.pairs = dataOut.pairsList
1207 1206 if hasattr(dataOut, 'meta'):
1208 self.meta = dataOut.meta
1207 self.meta.update(dataOut.meta)
1209 1208 self.channels = dataOut.channelList
1210 1209 self.interval = dataOut.getTimeInterval()
1211 1210 self.localtime = dataOut.useLocalTime
1212 1211 if 'spc' in self.plottypes or 'cspc' in self.plottypes or 'spc_moments' in self.plottypes:
1213 1212 self.xrange = (dataOut.getFreqRange(1)/1000.,
1214 1213 dataOut.getAcfRange(1), dataOut.getVelRange(1))
1215 1214 self.factor = dataOut.normFactor
1216 1215 self.__heights.append(dataOut.heightList)
1217 1216 self.__all_heights.update(dataOut.heightList)
1218 1217 self.__times.append(tm)
1219 1218
1220 1219 for plot in self.plottypes:
1221 1220 if plot in ('spc', 'spc_moments'):
1222 1221 z = dataOut.data_spc/dataOut.normFactor
1223 1222 buffer = 10*numpy.log10(z)
1224 1223 if plot == 'cspc':
1225 1224 z = dataOut.data_spc/dataOut.normFactor
1226 1225 buffer = (dataOut.data_spc, dataOut.data_cspc)
1227 1226 if plot == 'noise':
1228 1227 buffer = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
1229 1228 if plot == 'rti':
1230 1229 buffer = dataOut.getPower()
1231 1230 if plot == 'snr_db':
1232 1231 buffer = dataOut.data_SNR
1233 1232 if plot == 'snr':
1234 1233 buffer = 10*numpy.log10(dataOut.data_SNR)
1235 1234 if plot == 'dop':
1236 1235 buffer = 10*numpy.log10(dataOut.data_DOP)
1237 1236 if plot == 'mean':
1238 1237 buffer = dataOut.data_MEAN
1239 1238 if plot == 'std':
1240 1239 buffer = dataOut.data_STD
1241 1240 if plot == 'coh':
1242 1241 buffer = dataOut.getCoherence()
1243 1242 if plot == 'phase':
1244 1243 buffer = dataOut.getCoherence(phase=True)
1245 1244 if plot == 'output':
1246 1245 buffer = dataOut.data_output
1247 1246 if plot == 'param':
1248 1247 buffer = dataOut.data_param
1249 1248 if plot == 'scope':
1250 1249 buffer = dataOut.data
1251 1250 self.flagDataAsBlock = dataOut.flagDataAsBlock
1252 1251 self.nProfiles = dataOut.nProfiles
1253 1252
1254 1253 if plot == 'spc':
1255 1254 self.data['spc'] = buffer
1256 1255 elif plot == 'cspc':
1257 1256 self.data['spc'] = buffer[0]
1258 1257 self.data['cspc'] = buffer[1]
1259 1258 elif plot == 'spc_moments':
1260 1259 self.data['spc'] = buffer
1261 1260 self.data['moments'][tm] = dataOut.moments
1262 1261 else:
1263 1262 if self.buffering:
1264 1263 self.data[plot][tm] = buffer
1265 1264 else:
1266 1265 self.data[plot] = buffer
1267 1266
1268 1267 def normalize_heights(self):
1269 1268 '''
1270 1269 Ensure same-dimension of the data for different heighList
1271 1270 '''
1272 1271
1273 1272 H = numpy.array(list(self.__all_heights))
1274 1273 H.sort()
1275 1274 for key in self.data:
1276 1275 shape = self.shape(key)[:-1] + H.shape
1277 1276 for tm, obj in list(self.data[key].items()):
1278 1277 h = self.__heights[self.__times.index(tm)]
1279 1278 if H.size == h.size:
1280 1279 continue
1281 1280 index = numpy.where(numpy.in1d(H, h))[0]
1282 1281 dummy = numpy.zeros(shape) + numpy.nan
1283 1282 if len(shape) == 2:
1284 1283 dummy[:, index] = obj
1285 1284 else:
1286 1285 dummy[index] = obj
1287 1286 self.data[key][tm] = dummy
1288 1287
1289 1288 self.__heights = [H for tm in self.__times]
1290 1289
1291 def jsonify(self, decimate=False):
1290 def jsonify(self, plot_name, plot_type, decimate=False):
1292 1291 '''
1293 1292 Convert data to json
1294 1293 '''
1295 1294
1296 data = {}
1297 1295 tm = self.times[-1]
1298 1296 dy = int(self.heights.size/self.MAXNUMY) + 1
1299 for key in self.data:
1300 if key in ('spc', 'cspc') or not self.buffering:
1301 dx = int(self.data[key].shape[1]/self.MAXNUMX) + 1
1302 data[key] = self.roundFloats(
1303 self.data[key][::, ::dx, ::dy].tolist())
1304 else:
1305 data[key] = self.roundFloats(self.data[key][tm].tolist())
1306
1307 ret = {'data': data}
1308 ret['exp_code'] = self.exp_code
1309 ret['time'] = float(tm)
1310 ret['interval'] = float(self.interval)
1311 ret['localtime'] = self.localtime
1312 ret['yrange'] = self.roundFloats(self.heights[::dy].tolist())
1313 if 'spc' in self.data or 'cspc' in self.data:
1314 ret['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
1297 if self.key in ('spc', 'cspc') or not self.buffering:
1298 dx = int(self.data[self.key].shape[1]/self.MAXNUMX) + 1
1299 data = self.roundFloats(
1300 self.data[self.key][::, ::dx, ::dy].tolist())
1315 1301 else:
1316 ret['xrange'] = []
1317 if hasattr(self, 'pairs'):
1318 ret['pairs'] = [(int(p[0]), int(p[1])) for p in self.pairs]
1302 data = self.roundFloats(self.data[self.key][tm].tolist())
1303 if self.key is 'noise':
1304 data = [[x] for x in data]
1305
1306 meta = {}
1307 ret = {
1308 'plot': plot_name,
1309 'code': self.exp_code,
1310 'time': float(tm),
1311 'data': data,
1312 }
1313 meta['type'] = plot_type
1314 meta['interval'] = float(self.interval)
1315 meta['localtime'] = self.localtime
1316 meta['yrange'] = self.roundFloats(self.heights[::dy].tolist())
1317 if 'spc' in self.data or 'cspc' in self.data:
1318 meta['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
1319 1319 else:
1320 ret['pairs'] = []
1321
1322 for key, value in list(self.meta.items()):
1323 ret[key] = value
1320 meta['xrange'] = []
1324 1321
1322 meta.update(self.meta)
1323 ret['metadata'] = meta
1325 1324 return json.dumps(ret)
1326 1325
1327 1326 @property
1328 1327 def times(self):
1329 1328 '''
1330 1329 Return the list of times of the current data
1331 1330 '''
1332 1331
1333 1332 ret = numpy.array(self.__times)
1334 1333 ret.sort()
1335 1334 return ret
1336 1335
1337 1336 @property
1338 1337 def min_time(self):
1339 1338 '''
1340 1339 Return the minimun time value
1341 1340 '''
1342 1341
1343 1342 return self.times[0]
1344 1343
1345 1344 @property
1346 1345 def max_time(self):
1347 1346 '''
1348 1347 Return the maximun time value
1349 1348 '''
1350 1349
1351 1350 return self.times[-1]
1352 1351
1353 1352 @property
1354 1353 def heights(self):
1355 1354 '''
1356 1355 Return the list of heights of the current data
1357 1356 '''
1358 1357
1359 1358 return numpy.array(self.__heights[-1])
1360 1359
1361 1360 @staticmethod
1362 1361 def roundFloats(obj):
1363 1362 if isinstance(obj, list):
1364 1363 return list(map(PlotterData.roundFloats, obj))
1365 1364 elif isinstance(obj, float):
1366 1365 return round(obj, 2)
Show file before | Show file after | Hide comments
@@ -1,806 +1,808
1 1
2 2 import os
3 3 import sys
4 4 import zmq
5 5 import time
6 6 import numpy
7 7 import datetime
8 8 from functools import wraps
9 9 from threading import Thread
10 10 import matplotlib
11 11
12 12 if 'BACKEND' in os.environ:
13 13 matplotlib.use(os.environ['BACKEND'])
14 14 elif 'linux' in sys.platform:
15 15 matplotlib.use("TkAgg")
16 16 elif 'darwin' in sys.platform:
17 17 matplotlib.use('WxAgg')
18 18 else:
19 19 from schainpy.utils import log
20 20 log.warning('Using default Backend="Agg"', 'INFO')
21 21 matplotlib.use('Agg')
22 22
23 23 import matplotlib.pyplot as plt
24 24 from matplotlib.patches import Polygon
25 25 from mpl_toolkits.axes_grid1 import make_axes_locatable
26 26 from matplotlib.ticker import FuncFormatter, LinearLocator, MultipleLocator
27 27
28 28 from schainpy.model.data.jrodata import PlotterData
29 29 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
30 30 from schainpy.utils import log
31 31
32 32 jet_values = matplotlib.pyplot.get_cmap('jet', 100)(numpy.arange(100))[10:90]
33 33 blu_values = matplotlib.pyplot.get_cmap(
34 34 'seismic_r', 20)(numpy.arange(20))[10:15]
35 35 ncmap = matplotlib.colors.LinearSegmentedColormap.from_list(
36 36 'jro', numpy.vstack((blu_values, jet_values)))
37 37 matplotlib.pyplot.register_cmap(cmap=ncmap)
38 38
39 39 CMAPS = [plt.get_cmap(s) for s in ('jro', 'jet', 'viridis',
40 40 'plasma', 'inferno', 'Greys', 'seismic', 'bwr', 'coolwarm')]
41 41
42 42 EARTH_RADIUS = 6.3710e3
43 43
44 44 def ll2xy(lat1, lon1, lat2, lon2):
45 45
46 46 p = 0.017453292519943295
47 47 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
48 48 numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
49 49 r = 12742 * numpy.arcsin(numpy.sqrt(a))
50 50 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
51 51 * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
52 52 theta = -theta + numpy.pi/2
53 53 return r*numpy.cos(theta), r*numpy.sin(theta)
54 54
55 55
56 56 def km2deg(km):
57 57 '''
58 58 Convert distance in km to degrees
59 59 '''
60 60
61 61 return numpy.rad2deg(km/EARTH_RADIUS)
62 62
63 63
64 64 def figpause(interval):
65 65 backend = plt.rcParams['backend']
66 66 if backend in matplotlib.rcsetup.interactive_bk:
67 67 figManager = matplotlib._pylab_helpers.Gcf.get_active()
68 68 if figManager is not None:
69 69 canvas = figManager.canvas
70 70 if canvas.figure.stale:
71 71 canvas.draw()
72 72 try:
73 73 canvas.start_event_loop(interval)
74 74 except:
75 75 pass
76 76 return
77 77
78 78
79 79 def popup(message):
80 80 '''
81 81 '''
82 82
83 83 fig = plt.figure(figsize=(12, 8), facecolor='r')
84 84 text = '\n'.join([s.strip() for s in message.split(':')])
85 85 fig.text(0.01, 0.5, text, ha='left', va='center',
86 86 size='20', weight='heavy', color='w')
87 87 fig.show()
88 88 figpause(1000)
89 89
90 90
91 91 class Throttle(object):
92 92 '''
93 93 Decorator that prevents a function from being called more than once every
94 94 time period.
95 95 To create a function that cannot be called more than once a minute, but
96 96 will sleep until it can be called:
97 97 @Throttle(minutes=1)
98 98 def foo():
99 99 pass
100 100
101 101 for i in range(10):
102 102 foo()
103 103 print "This function has run %s times." % i
104 104 '''
105 105
106 106 def __init__(self, seconds=0, minutes=0, hours=0):
107 107 self.throttle_period = datetime.timedelta(
108 108 seconds=seconds, minutes=minutes, hours=hours
109 109 )
110 110
111 111 self.time_of_last_call = datetime.datetime.min
112 112
113 113 def __call__(self, fn):
114 114 @wraps(fn)
115 115 def wrapper(*args, **kwargs):
116 116 coerce = kwargs.pop('coerce', None)
117 117 if coerce:
118 118 self.time_of_last_call = datetime.datetime.now()
119 119 return fn(*args, **kwargs)
120 120 else:
121 121 now = datetime.datetime.now()
122 122 time_since_last_call = now - self.time_of_last_call
123 123 time_left = self.throttle_period - time_since_last_call
124 124
125 125 if time_left > datetime.timedelta(seconds=0):
126 126 return
127 127
128 128 self.time_of_last_call = datetime.datetime.now()
129 129 return fn(*args, **kwargs)
130 130
131 131 return wrapper
132 132
133 133 def apply_throttle(value):
134 134
135 135 @Throttle(seconds=value)
136 136 def fnThrottled(fn):
137 137 fn()
138 138
139 139 return fnThrottled
140 140
141 141
142 142 @MPDecorator
143 143 class Plot(Operation):
144 144 '''
145 145 Base class for Schain plotting operations
146 146 '''
147 147
148 148 CODE = 'Figure'
149 149 colormap = 'jet'
150 150 bgcolor = 'white'
151 151 __missing = 1E30
152 152
153 153 __attrs__ = ['show', 'save', 'xmin', 'xmax', 'ymin', 'ymax', 'zmin', 'zmax',
154 154 'zlimits', 'xlabel', 'ylabel', 'xaxis', 'cb_label', 'title',
155 155 'colorbar', 'bgcolor', 'width', 'height', 'localtime', 'oneFigure',
156 156 'showprofile', 'decimation', 'pause']
157 157
158 158 def __init__(self):
159 159
160 160 Operation.__init__(self)
161 161 self.isConfig = False
162 162 self.isPlotConfig = False
163 163 self.save_counter = 1
164 164 self.sender_counter = 1
165 165 self.data = None
166 166
167 167 def __fmtTime(self, x, pos):
168 168 '''
169 169 '''
170 170
171 171 return '{}'.format(self.getDateTime(x).strftime('%H:%M'))
172 172
173 173 def __setup(self, **kwargs):
174 174 '''
175 175 Initialize variables
176 176 '''
177 177
178 178 self.figures = []
179 179 self.axes = []
180 180 self.cb_axes = []
181 181 self.localtime = kwargs.pop('localtime', True)
182 182 self.show = kwargs.get('show', True)
183 183 self.save = kwargs.get('save', False)
184 184 self.save_period = kwargs.get('save_period', 2)
185 185 self.ftp = kwargs.get('ftp', False)
186 186 self.colormap = kwargs.get('colormap', self.colormap)
187 187 self.colormap_coh = kwargs.get('colormap_coh', 'jet')
188 188 self.colormap_phase = kwargs.get('colormap_phase', 'RdBu_r')
189 189 self.colormaps = kwargs.get('colormaps', None)
190 190 self.bgcolor = kwargs.get('bgcolor', self.bgcolor)
191 191 self.showprofile = kwargs.get('showprofile', False)
192 192 self.title = kwargs.get('wintitle', self.CODE.upper())
193 193 self.cb_label = kwargs.get('cb_label', None)
194 194 self.cb_labels = kwargs.get('cb_labels', None)
195 195 self.labels = kwargs.get('labels', None)
196 196 self.xaxis = kwargs.get('xaxis', 'frequency')
197 197 self.zmin = kwargs.get('zmin', None)
198 198 self.zmax = kwargs.get('zmax', None)
199 199 self.zlimits = kwargs.get('zlimits', None)
200 200 self.xmin = kwargs.get('xmin', None)
201 201 self.xmax = kwargs.get('xmax', None)
202 202 self.xrange = kwargs.get('xrange', 24)
203 203 self.xscale = kwargs.get('xscale', None)
204 204 self.ymin = kwargs.get('ymin', None)
205 205 self.ymax = kwargs.get('ymax', None)
206 206 self.yscale = kwargs.get('yscale', None)
207 207 self.xlabel = kwargs.get('xlabel', None)
208 208 self.decimation = kwargs.get('decimation', None)
209 209 self.showSNR = kwargs.get('showSNR', False)
210 210 self.oneFigure = kwargs.get('oneFigure', True)
211 211 self.width = kwargs.get('width', None)
212 212 self.height = kwargs.get('height', None)
213 213 self.colorbar = kwargs.get('colorbar', True)
214 214 self.factors = kwargs.get('factors', [1, 1, 1, 1, 1, 1, 1, 1])
215 215 self.channels = kwargs.get('channels', None)
216 216 self.titles = kwargs.get('titles', [])
217 217 self.polar = False
218 218 self.type = kwargs.get('type', 'iq')
219 219 self.grid = kwargs.get('grid', False)
220 220 self.pause = kwargs.get('pause', False)
221 221 self.save_labels = kwargs.get('save_labels', None)
222 222 self.realtime = kwargs.get('realtime', True)
223 223 self.buffering = kwargs.get('buffering', True)
224 224 self.throttle = kwargs.get('throttle', 2)
225 225 self.exp_code = kwargs.get('exp_code', None)
226 226 self.plot_server = kwargs.get('plot_server', False)
227 227 self.sender_period = kwargs.get('sender_period', 2)
228 228 self.__throttle_plot = apply_throttle(self.throttle)
229 229 self.data = PlotterData(
230 230 self.CODE, self.throttle, self.exp_code, self.buffering, snr=self.showSNR)
231 231
232 232 if self.plot_server:
233 233 if not self.plot_server.startswith('tcp://'):
234 234 self.plot_server = 'tcp://{}'.format(self.plot_server)
235 235 log.success(
236 236 'Sending to server: {}'.format(self.plot_server),
237 237 self.name
238 238 )
239 if 'plot_name' in kwargs:
240 self.plot_name = kwargs['plot_name']
239 241
240 242 def __setup_plot(self):
241 243 '''
242 244 Common setup for all figures, here figures and axes are created
243 245 '''
244 246
245 247 self.setup()
246 248
247 249 self.time_label = 'LT' if self.localtime else 'UTC'
248 250
249 251 if self.width is None:
250 252 self.width = 8
251 253
252 254 self.figures = []
253 255 self.axes = []
254 256 self.cb_axes = []
255 257 self.pf_axes = []
256 258 self.cmaps = []
257 259
258 260 size = '15%' if self.ncols == 1 else '30%'
259 261 pad = '4%' if self.ncols == 1 else '8%'
260 262
261 263 if self.oneFigure:
262 264 if self.height is None:
263 265 self.height = 1.4 * self.nrows + 1
264 266 fig = plt.figure(figsize=(self.width, self.height),
265 267 edgecolor='k',
266 268 facecolor='w')
267 269 self.figures.append(fig)
268 270 for n in range(self.nplots):
269 271 ax = fig.add_subplot(self.nrows, self.ncols,
270 272 n + 1, polar=self.polar)
271 273 ax.tick_params(labelsize=8)
272 274 ax.firsttime = True
273 275 ax.index = 0
274 276 ax.press = None
275 277 self.axes.append(ax)
276 278 if self.showprofile:
277 279 cax = self.__add_axes(ax, size=size, pad=pad)
278 280 cax.tick_params(labelsize=8)
279 281 self.pf_axes.append(cax)
280 282 else:
281 283 if self.height is None:
282 284 self.height = 3
283 285 for n in range(self.nplots):
284 286 fig = plt.figure(figsize=(self.width, self.height),
285 287 edgecolor='k',
286 288 facecolor='w')
287 289 ax = fig.add_subplot(1, 1, 1, polar=self.polar)
288 290 ax.tick_params(labelsize=8)
289 291 ax.firsttime = True
290 292 ax.index = 0
291 293 ax.press = None
292 294 self.figures.append(fig)
293 295 self.axes.append(ax)
294 296 if self.showprofile:
295 297 cax = self.__add_axes(ax, size=size, pad=pad)
296 298 cax.tick_params(labelsize=8)
297 299 self.pf_axes.append(cax)
298 300
299 301 for n in range(self.nrows):
300 302 if self.colormaps is not None:
301 303 cmap = plt.get_cmap(self.colormaps[n])
302 304 else:
303 305 cmap = plt.get_cmap(self.colormap)
304 306 cmap.set_bad(self.bgcolor, 1.)
305 307 self.cmaps.append(cmap)
306 308
307 309 for fig in self.figures:
308 310 fig.canvas.mpl_connect('key_press_event', self.OnKeyPress)
309 311 fig.canvas.mpl_connect('scroll_event', self.OnBtnScroll)
310 312 fig.canvas.mpl_connect('button_press_event', self.onBtnPress)
311 313 fig.canvas.mpl_connect('motion_notify_event', self.onMotion)
312 314 fig.canvas.mpl_connect('button_release_event', self.onBtnRelease)
313 315
314 316 def OnKeyPress(self, event):
315 317 '''
316 318 Event for pressing keys (up, down) change colormap
317 319 '''
318 320 ax = event.inaxes
319 321 if ax in self.axes:
320 322 if event.key == 'down':
321 323 ax.index += 1
322 324 elif event.key == 'up':
323 325 ax.index -= 1
324 326 if ax.index < 0:
325 327 ax.index = len(CMAPS) - 1
326 328 elif ax.index == len(CMAPS):
327 329 ax.index = 0
328 330 cmap = CMAPS[ax.index]
329 331 ax.cbar.set_cmap(cmap)
330 332 ax.cbar.draw_all()
331 333 ax.plt.set_cmap(cmap)
332 334 ax.cbar.patch.figure.canvas.draw()
333 335 self.colormap = cmap.name
334 336
335 337 def OnBtnScroll(self, event):
336 338 '''
337 339 Event for scrolling, scale figure
338 340 '''
339 341 cb_ax = event.inaxes
340 342 if cb_ax in [ax.cbar.ax for ax in self.axes if ax.cbar]:
341 343 ax = [ax for ax in self.axes if cb_ax == ax.cbar.ax][0]
342 344 pt = ax.cbar.ax.bbox.get_points()[:, 1]
343 345 nrm = ax.cbar.norm
344 346 vmin, vmax, p0, p1, pS = (
345 347 nrm.vmin, nrm.vmax, pt[0], pt[1], event.y)
346 348 scale = 2 if event.step == 1 else 0.5
347 349 point = vmin + (vmax - vmin) / (p1 - p0) * (pS - p0)
348 350 ax.cbar.norm.vmin = point - scale * (point - vmin)
349 351 ax.cbar.norm.vmax = point - scale * (point - vmax)
350 352 ax.plt.set_norm(ax.cbar.norm)
351 353 ax.cbar.draw_all()
352 354 ax.cbar.patch.figure.canvas.draw()
353 355
354 356 def onBtnPress(self, event):
355 357 '''
356 358 Event for mouse button press
357 359 '''
358 360 cb_ax = event.inaxes
359 361 if cb_ax is None:
360 362 return
361 363
362 364 if cb_ax in [ax.cbar.ax for ax in self.axes if ax.cbar]:
363 365 cb_ax.press = event.x, event.y
364 366 else:
365 367 cb_ax.press = None
366 368
367 369 def onMotion(self, event):
368 370 '''
369 371 Event for move inside colorbar
370 372 '''
371 373 cb_ax = event.inaxes
372 374 if cb_ax is None:
373 375 return
374 376 if cb_ax not in [ax.cbar.ax for ax in self.axes if ax.cbar]:
375 377 return
376 378 if cb_ax.press is None:
377 379 return
378 380
379 381 ax = [ax for ax in self.axes if cb_ax == ax.cbar.ax][0]
380 382 xprev, yprev = cb_ax.press
381 383 dx = event.x - xprev
382 384 dy = event.y - yprev
383 385 cb_ax.press = event.x, event.y
384 386 scale = ax.cbar.norm.vmax - ax.cbar.norm.vmin
385 387 perc = 0.03
386 388
387 389 if event.button == 1:
388 390 ax.cbar.norm.vmin -= (perc * scale) * numpy.sign(dy)
389 391 ax.cbar.norm.vmax -= (perc * scale) * numpy.sign(dy)
390 392 elif event.button == 3:
391 393 ax.cbar.norm.vmin -= (perc * scale) * numpy.sign(dy)
392 394 ax.cbar.norm.vmax += (perc * scale) * numpy.sign(dy)
393 395
394 396 ax.cbar.draw_all()
395 397 ax.plt.set_norm(ax.cbar.norm)
396 398 ax.cbar.patch.figure.canvas.draw()
397 399
398 400 def onBtnRelease(self, event):
399 401 '''
400 402 Event for mouse button release
401 403 '''
402 404 cb_ax = event.inaxes
403 405 if cb_ax is not None:
404 406 cb_ax.press = None
405 407
406 408 def __add_axes(self, ax, size='30%', pad='8%'):
407 409 '''
408 410 Add new axes to the given figure
409 411 '''
410 412 divider = make_axes_locatable(ax)
411 413 nax = divider.new_horizontal(size=size, pad=pad)
412 414 ax.figure.add_axes(nax)
413 415 return nax
414 416
415 417 def fill_gaps(self, x_buffer, y_buffer, z_buffer):
416 418 '''
417 419 Create a masked array for missing data
418 420 '''
419 421 if x_buffer.shape[0] < 2:
420 422 return x_buffer, y_buffer, z_buffer
421 423
422 424 deltas = x_buffer[1:] - x_buffer[0:-1]
423 425 x_median = numpy.median(deltas)
424 426
425 427 index = numpy.where(deltas > 5 * x_median)
426 428
427 429 if len(index[0]) != 0:
428 430 z_buffer[::, index[0], ::] = self.__missing
429 431 z_buffer = numpy.ma.masked_inside(z_buffer,
430 432 0.99 * self.__missing,
431 433 1.01 * self.__missing)
432 434
433 435 return x_buffer, y_buffer, z_buffer
434 436
435 437 def decimate(self):
436 438
437 439 # dx = int(len(self.x)/self.__MAXNUMX) + 1
438 440 dy = int(len(self.y) / self.decimation) + 1
439 441
440 442 # x = self.x[::dx]
441 443 x = self.x
442 444 y = self.y[::dy]
443 445 z = self.z[::, ::, ::dy]
444 446
445 447 return x, y, z
446 448
447 449 def format(self):
448 450 '''
449 451 Set min and max values, labels, ticks and titles
450 452 '''
451 453
452 454 if self.xmin is None:
453 455 xmin = self.data.min_time
454 456 else:
455 457 if self.xaxis is 'time':
456 458 dt = self.getDateTime(self.data.min_time)
457 459 xmin = (dt.replace(hour=int(self.xmin), minute=0, second=0) -
458 460 datetime.datetime(1970, 1, 1)).total_seconds()
459 461 if self.data.localtime:
460 462 xmin += time.timezone
461 463 else:
462 464 xmin = self.xmin
463 465
464 466 if self.xmax is None:
465 467 xmax = xmin + self.xrange * 60 * 60
466 468 else:
467 469 if self.xaxis is 'time':
468 470 dt = self.getDateTime(self.data.max_time)
469 471 xmax = (dt.replace(hour=int(self.xmax), minute=59, second=59) -
470 472 datetime.datetime(1970, 1, 1) + datetime.timedelta(seconds=1)).total_seconds()
471 473 if self.data.localtime:
472 474 xmax += time.timezone
473 475 else:
474 476 xmax = self.xmax
475 477
476 478 ymin = self.ymin if self.ymin else numpy.nanmin(self.y)
477 479 ymax = self.ymax if self.ymax else numpy.nanmax(self.y)
478 480 #Y = numpy.array([1, 2, 5, 10, 20, 50, 100, 200, 500, 1000, 2000, 5000, 10000, 20000, 50000])
479 481
480 482 #i = 1 if numpy.where(
481 483 # abs(ymax-ymin) <= Y)[0][0] < 0 else numpy.where(abs(ymax-ymin) <= Y)[0][0]
482 484 #ystep = Y[i] / 10.
483 485 dig = int(numpy.log10(ymax))
484 486 if dig == 0:
485 487 digD = len(str(ymax)) - 2
486 488 ydec = ymax*(10**digD)
487 489
488 490 dig = int(numpy.log10(ydec))
489 491 ystep = ((ydec + (10**(dig)))//10**(dig))*(10**(dig))
490 492 ystep = ystep/5
491 493 ystep = ystep/(10**digD)
492 494
493 495 else:
494 496 ystep = ((ymax + (10**(dig)))//10**(dig))*(10**(dig))
495 497 ystep = ystep/5
496 498
497 499 if self.xaxis is not 'time':
498 500
499 501 dig = int(numpy.log10(xmax))
500 502
501 503 if dig <= 0:
502 504 digD = len(str(xmax)) - 2
503 505 xdec = xmax*(10**digD)
504 506
505 507 dig = int(numpy.log10(xdec))
506 508 xstep = ((xdec + (10**(dig)))//10**(dig))*(10**(dig))
507 509 xstep = xstep*0.5
508 510 xstep = xstep/(10**digD)
509 511
510 512 else:
511 513 xstep = ((xmax + (10**(dig)))//10**(dig))*(10**(dig))
512 514 xstep = xstep/5
513 515
514 516 for n, ax in enumerate(self.axes):
515 517 if ax.firsttime:
516 518 ax.set_facecolor(self.bgcolor)
517 519 ax.yaxis.set_major_locator(MultipleLocator(ystep))
518 520 if self.xscale:
519 521 ax.xaxis.set_major_formatter(FuncFormatter(
520 522 lambda x, pos: '{0:g}'.format(x*self.xscale)))
521 523 if self.xscale:
522 524 ax.yaxis.set_major_formatter(FuncFormatter(
523 525 lambda x, pos: '{0:g}'.format(x*self.yscale)))
524 526 if self.xaxis is 'time':
525 527 ax.xaxis.set_major_formatter(FuncFormatter(self.__fmtTime))
526 528 ax.xaxis.set_major_locator(LinearLocator(9))
527 529 else:
528 530 ax.xaxis.set_major_locator(MultipleLocator(xstep))
529 531 if self.xlabel is not None:
530 532 ax.set_xlabel(self.xlabel)
531 533 ax.set_ylabel(self.ylabel)
532 534 ax.firsttime = False
533 535 if self.showprofile:
534 536 self.pf_axes[n].set_ylim(ymin, ymax)
535 537 self.pf_axes[n].set_xlim(self.zmin, self.zmax)
536 538 self.pf_axes[n].set_xlabel('dB')
537 539 self.pf_axes[n].grid(b=True, axis='x')
538 540 [tick.set_visible(False)
539 541 for tick in self.pf_axes[n].get_yticklabels()]
540 542 if self.colorbar:
541 543 ax.cbar = plt.colorbar(
542 544 ax.plt, ax=ax, fraction=0.05, pad=0.02, aspect=10)
543 545 ax.cbar.ax.tick_params(labelsize=8)
544 546 ax.cbar.ax.press = None
545 547 if self.cb_label:
546 548 ax.cbar.set_label(self.cb_label, size=8)
547 549 elif self.cb_labels:
548 550 ax.cbar.set_label(self.cb_labels[n], size=8)
549 551 else:
550 552 ax.cbar = None
551 553 if self.grid:
552 554 ax.grid(True)
553 555
554 556 if not self.polar:
555 557 ax.set_xlim(xmin, xmax)
556 558 ax.set_ylim(ymin, ymax)
557 559 ax.set_title('{} {} {}'.format(
558 560 self.titles[n],
559 561 self.getDateTime(self.data.max_time).strftime(
560 562 '%Y-%m-%d %H:%M:%S'),
561 563 self.time_label),
562 564 size=8)
563 565 else:
564 566 ax.set_title('{}'.format(self.titles[n]), size=8)
565 567 ax.set_ylim(0, 90)
566 568 ax.set_yticks(numpy.arange(0, 90, 20))
567 569 ax.yaxis.labelpad = 40
568 570
569 571 def clear_figures(self):
570 572 '''
571 573 Reset axes for redraw plots
572 574 '''
573 575
574 576 for ax in self.axes:
575 577 ax.clear()
576 578 ax.firsttime = True
577 579 if ax.cbar:
578 580 ax.cbar.remove()
579 581
580 582 def __plot(self):
581 583 '''
582 584 Main function to plot, format and save figures
583 585 '''
584 586
585 587 try:
586 588 self.plot()
587 589 self.format()
588 590 except Exception as e:
589 591 log.warning('{} Plot could not be updated... check data'.format(
590 592 self.CODE), self.name)
591 593 log.error(str(e), '')
592 594 return
593 595
594 596 for n, fig in enumerate(self.figures):
595 597 if self.nrows == 0 or self.nplots == 0:
596 598 log.warning('No data', self.name)
597 599 fig.text(0.5, 0.5, 'No Data', fontsize='large', ha='center')
598 600 fig.canvas.manager.set_window_title(self.CODE)
599 601 continue
600 602
601 603 fig.tight_layout()
602 604 fig.canvas.manager.set_window_title('{} - {}'.format(self.title,
603 605 self.getDateTime(self.data.max_time).strftime('%Y/%m/%d')))
604 606 fig.canvas.draw()
605 607 if self.show:
606 608 fig.show()
607 609 figpause(0.1)
608 610
609 611 if self.save:
610 612 self.save_figure(n)
611 613
612 614 if self.plot_server:
613 615 self.send_to_server()
614 616 # t = Thread(target=self.send_to_server)
615 617 # t.start()
616 618
617 619 def save_figure(self, n):
618 620 '''
619 621 '''
620 622
621 623 if self.save_counter < self.save_period:
622 624 self.save_counter += 1
623 625 return
624 626
625 627 self.save_counter = 1
626 628
627 629 fig = self.figures[n]
628 630
629 631 if self.save_labels:
630 632 labels = self.save_labels
631 633 else:
632 634 labels = list(range(self.nrows))
633 635
634 636 if self.oneFigure:
635 637 label = ''
636 638 else:
637 639 label = '-{}'.format(labels[n])
638 640 figname = os.path.join(
639 641 self.save,
640 642 self.CODE,
641 643 '{}{}_{}.png'.format(
642 644 self.CODE,
643 645 label,
644 646 self.getDateTime(self.data.max_time).strftime(
645 647 '%Y%m%d_%H%M%S'
646 648 ),
647 649 )
648 650 )
649 651 log.log('Saving figure: {}'.format(figname), self.name)
650 652 if not os.path.isdir(os.path.dirname(figname)):
651 653 os.makedirs(os.path.dirname(figname))
652 654 fig.savefig(figname)
653 655
654 656 if self.realtime:
655 657 figname = os.path.join(
656 658 self.save,
657 659 '{}{}_{}.png'.format(
658 660 self.CODE,
659 661 label,
660 662 self.getDateTime(self.data.min_time).strftime(
661 663 '%Y%m%d'
662 664 ),
663 665 )
664 666 )
665 667 fig.savefig(figname)
666 668
667 669 def send_to_server(self):
668 670 '''
669 671 '''
670 672
671 673 if self.sender_counter < self.sender_period:
672 674 self.sender_counter += 1
673 675
674 676 self.sender_counter = 1
675
677 self.data.meta['titles'] = self.titles
676 678 retries = 2
677 679 while True:
678 self.socket.send_string(self.data.jsonify())
680 self.socket.send_string(self.data.jsonify(self.plot_name, self.plot_type))
679 681 socks = dict(self.poll.poll(5000))
680 682 if socks.get(self.socket) == zmq.POLLIN:
681 683 reply = self.socket.recv_string()
682 684 if reply == 'ok':
683 685 log.log("Response from server ok", self.name)
684 686 break
685 687 else:
686 688 log.warning(
687 689 "Malformed reply from server: {}".format(reply), self.name)
688 690
689 691 else:
690 692 log.warning(
691 693 "No response from server, retrying...", self.name)
692 694 self.socket.setsockopt(zmq.LINGER, 0)
693 695 self.socket.close()
694 696 self.poll.unregister(self.socket)
695 697 retries -= 1
696 698 if retries == 0:
697 699 log.error(
698 700 "Server seems to be offline, abandoning", self.name)
699 701 self.socket = self.context.socket(zmq.REQ)
700 702 self.socket.connect(self.plot_server)
701 703 self.poll.register(self.socket, zmq.POLLIN)
702 704 time.sleep(1)
703 705 break
704 706 self.socket = self.context.socket(zmq.REQ)
705 707 self.socket.connect(self.plot_server)
706 708 self.poll.register(self.socket, zmq.POLLIN)
707 709 time.sleep(0.5)
708 710
709 711 def setup(self):
710 712 '''
711 713 This method should be implemented in the child class, the following
712 714 attributes should be set:
713 715
714 716 self.nrows: number of rows
715 717 self.ncols: number of cols
716 718 self.nplots: number of plots (channels or pairs)
717 719 self.ylabel: label for Y axes
718 720 self.titles: list of axes title
719 721
720 722 '''
721 723 raise NotImplementedError
722 724
723 725 def plot(self):
724 726 '''
725 727 Must be defined in the child class
726 728 '''
727 729 raise NotImplementedError
728 730
729 731 def run(self, dataOut, **kwargs):
730 732 '''
731 733 Main plotting routine
732 734 '''
733 735
734 736 if self.isConfig is False:
735 737 self.__setup(**kwargs)
736 738 if dataOut.type == 'Parameters':
737 739 t = dataOut.utctimeInit
738 740 else:
739 741 t = dataOut.utctime
740 742
741 743 if dataOut.useLocalTime:
742 744 self.getDateTime = datetime.datetime.fromtimestamp
743 745 if not self.localtime:
744 746 t += time.timezone
745 747 else:
746 748 self.getDateTime = datetime.datetime.utcfromtimestamp
747 749 if self.localtime:
748 750 t -= time.timezone
749 751
750 752 if self.xmin is None:
751 753 self.tmin = t
752 754 else:
753 755 self.tmin = (
754 756 self.getDateTime(t).replace(
755 757 hour=self.xmin,
756 758 minute=0,
757 759 second=0) - self.getDateTime(0)).total_seconds()
758 760
759 761 self.data.setup()
760 762 self.isConfig = True
761 763 if self.plot_server:
762 764 self.context = zmq.Context()
763 765 self.socket = self.context.socket(zmq.REQ)
764 766 self.socket.connect(self.plot_server)
765 767 self.poll = zmq.Poller()
766 768 self.poll.register(self.socket, zmq.POLLIN)
767 769
768 770 if dataOut.type == 'Parameters':
769 771 tm = dataOut.utctimeInit
770 772 else:
771 773 tm = dataOut.utctime
772 774
773 775 if not dataOut.useLocalTime and self.localtime:
774 776 tm -= time.timezone
775 777 if dataOut.useLocalTime and not self.localtime:
776 778 tm += time.timezone
777 779
778 780 if self.xaxis is 'time' and self.data and (tm - self.tmin) >= self.xrange*60*60:
779 781 self.save_counter = self.save_period
780 782 self.__plot()
781 783 self.xmin += self.xrange
782 784 if self.xmin >= 24:
783 785 self.xmin -= 24
784 786 self.tmin += self.xrange*60*60
785 787 self.data.setup()
786 788 self.clear_figures()
787 789
788 790 self.data.update(dataOut, tm)
789 791
790 792 if self.isPlotConfig is False:
791 793 self.__setup_plot()
792 794 self.isPlotConfig = True
793 795
794 796 if self.realtime:
795 797 self.__plot()
796 798 else:
797 799 self.__throttle_plot(self.__plot)#, coerce=coerce)
798 800
799 801 def close(self):
800 802
801 803 if self.data:
802 804 self.save_counter = self.save_period
803 805 self.__plot()
804 806 if self.data and self.pause:
805 807 figpause(10)
806 808
Show file before | Show file after | Hide comments
@@ -1,748 +1,759
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 plot_name = 'Spectra'
46 plot_type = 'pcolor'
45 47
46 48 def setup(self):
47 49 self.nplots = len(self.data.channels)
48 50 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
49 51 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
50 52 self.width = 3.4 * self.ncols
51 53 self.height = 3 * self.nrows
52 54 self.cb_label = 'dB'
53 55 if self.showprofile:
54 56 self.width += 0.8 * self.ncols
55 57
56 58 self.ylabel = 'Range [km]'
57 59
58 60 def plot(self):
59 61 if self.xaxis == "frequency":
60 62 x = self.data.xrange[0]
61 63 self.xlabel = "Frequency (kHz)"
62 64 elif self.xaxis == "time":
63 65 x = self.data.xrange[1]
64 66 self.xlabel = "Time (ms)"
65 67 else:
66 68 x = self.data.xrange[2]
67 69 self.xlabel = "Velocity (m/s)"
68 70
69 71 if self.CODE == 'spc_moments':
70 72 x = self.data.xrange[2]
71 73 self.xlabel = "Velocity (m/s)"
72 74
73 75 self.titles = []
74 76
75 77 y = self.data.heights
76 78 self.y = y
77 79 z = self.data['spc']
78 80
79 81 for n, ax in enumerate(self.axes):
80 82 noise = self.data['noise'][n][-1]
81 83 if self.CODE == 'spc_moments':
82 84 mean = self.data['moments'][n, :, 1, :][-1]
83 85 if ax.firsttime:
84 86 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
85 87 self.xmin = self.xmin if self.xmin else -self.xmax
86 88 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
87 89 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
88 90 ax.plt = ax.pcolormesh(x, y, z[n].T,
89 91 vmin=self.zmin,
90 92 vmax=self.zmax,
91 93 cmap=plt.get_cmap(self.colormap)
92 94 )
93 95
94 96 if self.showprofile:
95 97 ax.plt_profile = self.pf_axes[n].plot(
96 98 self.data['rti'][n][-1], y)[0]
97 99 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
98 100 color="k", linestyle="dashed", lw=1)[0]
99 101 if self.CODE == 'spc_moments':
100 102 ax.plt_mean = ax.plot(mean, y, color='k')[0]
101 103 else:
102 104 ax.plt.set_array(z[n].T.ravel())
103 105 if self.showprofile:
104 106 ax.plt_profile.set_data(self.data['rti'][n][-1], y)
105 107 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
106 108 if self.CODE == 'spc_moments':
107 109 ax.plt_mean.set_data(mean, y)
108 110 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
109 111
110 112
111 113 class CrossSpectraPlot(Plot):
112 114
113 115 CODE = 'cspc'
114 116 colormap = 'jet'
117 plot_name = 'CrossSpectra'
118 plot_type = 'pcolor'
115 119 zmin_coh = None
116 120 zmax_coh = None
117 121 zmin_phase = None
118 122 zmax_phase = None
119 123
120 124 def setup(self):
121 125
122 126 self.ncols = 4
123 127 self.nrows = len(self.data.pairs)
124 128 self.nplots = self.nrows * 4
125 129 self.width = 3.4 * self.ncols
126 130 self.height = 3 * self.nrows
127 131 self.ylabel = 'Range [km]'
128 132 self.showprofile = False
129 133
130 134 def plot(self):
131 135
132 136 if self.xaxis == "frequency":
133 137 x = self.data.xrange[0]
134 138 self.xlabel = "Frequency (kHz)"
135 139 elif self.xaxis == "time":
136 140 x = self.data.xrange[1]
137 141 self.xlabel = "Time (ms)"
138 142 else:
139 143 x = self.data.xrange[2]
140 144 self.xlabel = "Velocity (m/s)"
141 145
142 146 self.titles = []
143 147
144 148 y = self.data.heights
145 149 self.y = y
146 150 spc = self.data['spc']
147 151 cspc = self.data['cspc']
148 152
149 153 for n in range(self.nrows):
150 154 noise = self.data['noise'][n][-1]
151 155 pair = self.data.pairs[n]
152 156 ax = self.axes[4 * n]
153 157 spc0 = 10.*numpy.log10(spc[pair[0]]/self.data.factor)
154 158 if ax.firsttime:
155 159 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
156 160 self.xmin = self.xmin if self.xmin else -self.xmax
157 161 self.zmin = self.zmin if self.zmin else numpy.nanmin(spc)
158 162 self.zmax = self.zmax if self.zmax else numpy.nanmax(spc)
159 163 ax.plt = ax.pcolormesh(x , y , spc0.T,
160 164 vmin=self.zmin,
161 165 vmax=self.zmax,
162 166 cmap=plt.get_cmap(self.colormap)
163 167 )
164 168 else:
165 169 ax.plt.set_array(spc0.T.ravel())
166 170 self.titles.append('CH {}: {:3.2f}dB'.format(pair[0], noise))
167 171
168 172 ax = self.axes[4 * n + 1]
169 173 spc1 = 10.*numpy.log10(spc[pair[1]]/self.data.factor)
170 174 if ax.firsttime:
171 175 ax.plt = ax.pcolormesh(x , y, spc1.T,
172 176 vmin=self.zmin,
173 177 vmax=self.zmax,
174 178 cmap=plt.get_cmap(self.colormap)
175 179 )
176 180 else:
177 181 ax.plt.set_array(spc1.T.ravel())
178 182 self.titles.append('CH {}: {:3.2f}dB'.format(pair[1], noise))
179 183
180 184 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
181 185 coh = numpy.abs(out)
182 186 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
183 187
184 188 ax = self.axes[4 * n + 2]
185 189 if ax.firsttime:
186 190 ax.plt = ax.pcolormesh(x, y, coh.T,
187 191 vmin=0,
188 192 vmax=1,
189 193 cmap=plt.get_cmap(self.colormap_coh)
190 194 )
191 195 else:
192 196 ax.plt.set_array(coh.T.ravel())
193 197 self.titles.append(
194 198 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
195 199
196 200 ax = self.axes[4 * n + 3]
197 201 if ax.firsttime:
198 202 ax.plt = ax.pcolormesh(x, y, phase.T,
199 203 vmin=-180,
200 204 vmax=180,
201 205 cmap=plt.get_cmap(self.colormap_phase)
202 206 )
203 207 else:
204 208 ax.plt.set_array(phase.T.ravel())
205 209 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
206 210
207 211
208 212 class SpectralMomentsPlot(SpectraPlot):
209 213 '''
210 214 Plot for Spectral Moments
211 215 '''
212 216 CODE = 'spc_moments'
213 217 colormap = 'jro'
218 plot_name = 'SpectralMoments'
219 plot_type = 'pcolor'
214 220
215 221
216 222 class RTIPlot(Plot):
217 223 '''
218 224 Plot for RTI data
219 225 '''
220 226
221 227 CODE = 'rti'
222 228 colormap = 'jro'
229 plot_name = 'RTI'
230 plot_type = 'pcolorbuffer'
223 231
224 232 def setup(self):
225 233 self.xaxis = 'time'
226 234 self.ncols = 1
227 235 self.nrows = len(self.data.channels)
228 236 self.nplots = len(self.data.channels)
229 237 self.ylabel = 'Range [km]'
230 238 self.cb_label = 'dB'
231 239 self.titles = ['{} Channel {}'.format(
232 240 self.CODE.upper(), x) for x in range(self.nrows)]
233 241
234 242 def plot(self):
235 243 self.x = self.data.times
236 244 self.y = self.data.heights
237 245 self.z = self.data[self.CODE]
238 246 self.z = numpy.ma.masked_invalid(self.z)
239 247
240 248 if self.decimation is None:
241 249 x, y, z = self.fill_gaps(self.x, self.y, self.z)
242 250 else:
243 251 x, y, z = self.fill_gaps(*self.decimate())
244 252
245 253 for n, ax in enumerate(self.axes):
246 254 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
247 255 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
248 256 if ax.firsttime:
249 257 ax.plt = ax.pcolormesh(x, y, z[n].T,
250 258 vmin=self.zmin,
251 259 vmax=self.zmax,
252 260 cmap=plt.get_cmap(self.colormap)
253 261 )
254 262 if self.showprofile:
255 263 ax.plot_profile = self.pf_axes[n].plot(
256 264 self.data['rti'][n][-1], self.y)[0]
257 265 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(self.data['noise'][n][-1], len(self.y)), self.y,
258 266 color="k", linestyle="dashed", lw=1)[0]
259 267 else:
260 268 ax.collections.remove(ax.collections[0])
261 269 ax.plt = ax.pcolormesh(x, y, z[n].T,
262 270 vmin=self.zmin,
263 271 vmax=self.zmax,
264 272 cmap=plt.get_cmap(self.colormap)
265 273 )
266 274 if self.showprofile:
267 275 ax.plot_profile.set_data(self.data['rti'][n][-1], self.y)
268 276 ax.plot_noise.set_data(numpy.repeat(
269 277 self.data['noise'][n][-1], len(self.y)), self.y)
270 278
271 279
272 280 class CoherencePlot(RTIPlot):
273 281 '''
274 282 Plot for Coherence data
275 283 '''
276 284
277 285 CODE = 'coh'
286 plot_name = 'Coherence'
278 287
279 288 def setup(self):
280 289 self.xaxis = 'time'
281 290 self.ncols = 1
282 291 self.nrows = len(self.data.pairs)
283 292 self.nplots = len(self.data.pairs)
284 293 self.ylabel = 'Range [km]'
285 294 if self.CODE == 'coh':
286 295 self.cb_label = ''
287 296 self.titles = [
288 297 'Coherence Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
289 298 else:
290 299 self.cb_label = 'Degrees'
291 300 self.titles = [
292 301 'Phase Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
293 302
294 303
295 304 class PhasePlot(CoherencePlot):
296 305 '''
297 306 Plot for Phase map data
298 307 '''
299 308
300 309 CODE = 'phase'
301 310 colormap = 'seismic'
311 plot_name = 'Phase'
302 312
303 313
304 314 class NoisePlot(Plot):
305 315 '''
306 316 Plot for noise
307 317 '''
308 318
309 319 CODE = 'noise'
320 plot_name = 'Noise'
321 plot_type = 'scatterbuffer'
322
310 323
311 324 def setup(self):
312 325 self.xaxis = 'time'
313 326 self.ncols = 1
314 327 self.nrows = 1
315 328 self.nplots = 1
316 329 self.ylabel = 'Intensity [dB]'
317 330 self.titles = ['Noise']
318 331 self.colorbar = False
319 332
320 333 def plot(self):
321 334
322 335 x = self.data.times
323 336 xmin = self.data.min_time
324 337 xmax = xmin + self.xrange * 60 * 60
325 338 Y = self.data[self.CODE]
326 339
327 340 if self.axes[0].firsttime:
328 341 for ch in self.data.channels:
329 342 y = Y[ch]
330 343 self.axes[0].plot(x, y, lw=1, label='Ch{}'.format(ch))
331 344 plt.legend()
332 345 else:
333 346 for ch in self.data.channels:
334 347 y = Y[ch]
335 348 self.axes[0].lines[ch].set_data(x, y)
336 349
337 350 self.ymin = numpy.nanmin(Y) - 5
338 351 self.ymax = numpy.nanmax(Y) + 5
339 352
340 353
341 354 class SnrPlot(RTIPlot):
342 355 '''
343 356 Plot for SNR Data
344 357 '''
345 358
346 359 CODE = 'snr'
347 360 colormap = 'jet'
361 plot_name = 'SNR'
348 362
349 363
350 364 class DopplerPlot(RTIPlot):
351 365 '''
352 366 Plot for DOPPLER Data
353 367 '''
354 368
355 369 CODE = 'dop'
356 370 colormap = 'jet'
371 plot_name = 'Doppler'
357 372
358 373
359 374 class SkyMapPlot(Plot):
360 375 '''
361 376 Plot for meteors detection data
362 377 '''
363 378
364 379 CODE = 'param'
365 380
366 381 def setup(self):
367 382
368 383 self.ncols = 1
369 384 self.nrows = 1
370 385 self.width = 7.2
371 386 self.height = 7.2
372 387 self.nplots = 1
373 388 self.xlabel = 'Zonal Zenith Angle (deg)'
374 389 self.ylabel = 'Meridional Zenith Angle (deg)'
375 390 self.polar = True
376 391 self.ymin = -180
377 392 self.ymax = 180
378 393 self.colorbar = False
379 394
380 395 def plot(self):
381 396
382 397 arrayParameters = numpy.concatenate(self.data['param'])
383 398 error = arrayParameters[:, -1]
384 399 indValid = numpy.where(error == 0)[0]
385 400 finalMeteor = arrayParameters[indValid, :]
386 401 finalAzimuth = finalMeteor[:, 3]
387 402 finalZenith = finalMeteor[:, 4]
388 403
389 404 x = finalAzimuth * numpy.pi / 180
390 405 y = finalZenith
391 406
392 407 ax = self.axes[0]
393 408
394 409 if ax.firsttime:
395 410 ax.plot = ax.plot(x, y, 'bo', markersize=5)[0]
396 411 else:
397 412 ax.plot.set_data(x, y)
398 413
399 414 dt1 = self.getDateTime(self.data.min_time).strftime('%y/%m/%d %H:%M:%S')
400 415 dt2 = self.getDateTime(self.data.max_time).strftime('%y/%m/%d %H:%M:%S')
401 416 title = 'Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n' % (dt1,
402 417 dt2,
403 418 len(x))
404 419 self.titles[0] = title
405 420
406 421
407 422 class ParametersPlot(RTIPlot):
408 423 '''
409 424 Plot for data_param object
410 425 '''
411 426
412 427 CODE = 'param'
413 428 colormap = 'seismic'
429 plot_name = 'Parameters'
414 430
415 431 def setup(self):
416 432 self.xaxis = 'time'
417 433 self.ncols = 1
418 434 self.nrows = self.data.shape(self.CODE)[0]
419 435 self.nplots = self.nrows
420 436 if self.showSNR:
421 437 self.nrows += 1
422 438 self.nplots += 1
423 439
424 440 self.ylabel = 'Height [km]'
425 441 if not self.titles:
426 442 self.titles = self.data.parameters \
427 443 if self.data.parameters else ['Param {}'.format(x) for x in range(self.nrows)]
428 444 if self.showSNR:
429 445 self.titles.append('SNR')
430 446
431 447 def plot(self):
432 448 self.data.normalize_heights()
433 449 self.x = self.data.times
434 450 self.y = self.data.heights
435 451 if self.showSNR:
436 452 self.z = numpy.concatenate(
437 453 (self.data[self.CODE], self.data['snr'])
438 454 )
439 455 else:
440 456 self.z = self.data[self.CODE]
441 457
442 458 self.z = numpy.ma.masked_invalid(self.z)
443 459
444 460 if self.decimation is None:
445 461 x, y, z = self.fill_gaps(self.x, self.y, self.z)
446 462 else:
447 463 x, y, z = self.fill_gaps(*self.decimate())
448 464
449 465 for n, ax in enumerate(self.axes):
450 466
451 467 self.zmax = self.zmax if self.zmax is not None else numpy.max(
452 468 self.z[n])
453 469 self.zmin = self.zmin if self.zmin is not None else numpy.min(
454 470 self.z[n])
455 471
456 472 if ax.firsttime:
457 473 if self.zlimits is not None:
458 474 self.zmin, self.zmax = self.zlimits[n]
459 475
460 476 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
461 477 vmin=self.zmin,
462 478 vmax=self.zmax,
463 479 cmap=self.cmaps[n]
464 480 )
465 481 else:
466 482 if self.zlimits is not None:
467 483 self.zmin, self.zmax = self.zlimits[n]
468 484 ax.collections.remove(ax.collections[0])
469 485 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
470 486 vmin=self.zmin,
471 487 vmax=self.zmax,
472 488 cmap=self.cmaps[n]
473 489 )
474 490
475 491
476 492 class OutputPlot(ParametersPlot):
477 493 '''
478 494 Plot data_output object
479 495 '''
480 496
481 497 CODE = 'output'
482 498 colormap = 'seismic'
499 plot_name = 'Output'
483 500
484 501
485 502 class PolarMapPlot(Plot):
486 503 '''
487 504 Plot for weather radar
488 505 '''
489 506
490 507 CODE = 'param'
491 508 colormap = 'seismic'
492 509
493 510 def setup(self):
494 511 self.ncols = 1
495 512 self.nrows = 1
496 513 self.width = 9
497 514 self.height = 8
498 515 self.mode = self.data.meta['mode']
499 516 if self.channels is not None:
500 517 self.nplots = len(self.channels)
501 518 self.nrows = len(self.channels)
502 519 else:
503 520 self.nplots = self.data.shape(self.CODE)[0]
504 521 self.nrows = self.nplots
505 522 self.channels = list(range(self.nplots))
506 523 if self.mode == 'E':
507 524 self.xlabel = 'Longitude'
508 525 self.ylabel = 'Latitude'
509 526 else:
510 527 self.xlabel = 'Range (km)'
511 528 self.ylabel = 'Height (km)'
512 529 self.bgcolor = 'white'
513 530 self.cb_labels = self.data.meta['units']
514 531 self.lat = self.data.meta['latitude']
515 532 self.lon = self.data.meta['longitude']
516 533 self.xmin, self.xmax = float(
517 534 km2deg(self.xmin) + self.lon), float(km2deg(self.xmax) + self.lon)
518 535 self.ymin, self.ymax = float(
519 536 km2deg(self.ymin) + self.lat), float(km2deg(self.ymax) + self.lat)
520 537 # self.polar = True
521 538
522 539 def plot(self):
523 540
524 541 for n, ax in enumerate(self.axes):
525 542 data = self.data['param'][self.channels[n]]
526 543
527 544 zeniths = numpy.linspace(
528 545 0, self.data.meta['max_range'], data.shape[1])
529 546 if self.mode == 'E':
530 547 azimuths = -numpy.radians(self.data.heights)+numpy.pi/2
531 548 r, theta = numpy.meshgrid(zeniths, azimuths)
532 549 x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])), r*numpy.sin(
533 550 theta)*numpy.cos(numpy.radians(self.data.meta['elevation']))
534 551 x = km2deg(x) + self.lon
535 552 y = km2deg(y) + self.lat
536 553 else:
537 554 azimuths = numpy.radians(self.data.heights)
538 555 r, theta = numpy.meshgrid(zeniths, azimuths)
539 556 x, y = r*numpy.cos(theta), r*numpy.sin(theta)
540 557 self.y = zeniths
541 558
542 559 if ax.firsttime:
543 560 if self.zlimits is not None:
544 561 self.zmin, self.zmax = self.zlimits[n]
545 562 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
546 563 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
547 564 vmin=self.zmin,
548 565 vmax=self.zmax,
549 566 cmap=self.cmaps[n])
550 567 else:
551 568 if self.zlimits is not None:
552 569 self.zmin, self.zmax = self.zlimits[n]
553 570 ax.collections.remove(ax.collections[0])
554 571 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
555 572 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
556 573 vmin=self.zmin,
557 574 vmax=self.zmax,
558 575 cmap=self.cmaps[n])
559 576
560 577 if self.mode == 'A':
561 578 continue
562 579
563 580 # plot district names
564 581 f = open('/data/workspace/schain_scripts/distrito.csv')
565 582 for line in f:
566 583 label, lon, lat = [s.strip() for s in line.split(',') if s]
567 584 lat = float(lat)
568 585 lon = float(lon)
569 586 # ax.plot(lon, lat, '.b', ms=2)
570 587 ax.text(lon, lat, label.decode('utf8'), ha='center',
571 588 va='bottom', size='8', color='black')
572 589
573 590 # plot limites
574 591 limites = []
575 592 tmp = []
576 593 for line in open('/data/workspace/schain_scripts/lima.csv'):
577 594 if '#' in line:
578 595 if tmp:
579 596 limites.append(tmp)
580 597 tmp = []
581 598 continue
582 599 values = line.strip().split(',')
583 600 tmp.append((float(values[0]), float(values[1])))
584 601 for points in limites:
585 602 ax.add_patch(
586 603 Polygon(points, ec='k', fc='none', ls='--', lw=0.5))
587 604
588 605 # plot Cuencas
589 606 for cuenca in ('rimac', 'lurin', 'mala', 'chillon', 'chilca', 'chancay-huaral'):
590 607 f = open('/data/workspace/schain_scripts/{}.csv'.format(cuenca))
591 608 values = [line.strip().split(',') for line in f]
592 609 points = [(float(s[0]), float(s[1])) for s in values]
593 610 ax.add_patch(Polygon(points, ec='b', fc='none'))
594 611
595 612 # plot grid
596 613 for r in (15, 30, 45, 60):
597 614 ax.add_artist(plt.Circle((self.lon, self.lat),
598 615 km2deg(r), color='0.6', fill=False, lw=0.2))
599 616 ax.text(
600 617 self.lon + (km2deg(r))*numpy.cos(60*numpy.pi/180),
601 618 self.lat + (km2deg(r))*numpy.sin(60*numpy.pi/180),
602 619 '{}km'.format(r),
603 620 ha='center', va='bottom', size='8', color='0.6', weight='heavy')
604 621
605 622 if self.mode == 'E':
606 623 title = 'El={}$^\circ$'.format(self.data.meta['elevation'])
607 624 label = 'E{:02d}'.format(int(self.data.meta['elevation']))
608 625 else:
609 626 title = 'Az={}$^\circ$'.format(self.data.meta['azimuth'])
610 627 label = 'A{:02d}'.format(int(self.data.meta['azimuth']))
611 628
612 629 self.save_labels = ['{}-{}'.format(lbl, label) for lbl in self.labels]
613 630 self.titles = ['{} {}'.format(
614 631 self.data.parameters[x], title) for x in self.channels]
615 632
616 633
617 634 class ScopePlot(Plot):
618 635
619 636 '''
620 637 Plot for Scope
621 638 '''
622 639
623 640 CODE = 'scope'
641 plot_name = 'Scope'
642 plot_type = 'scatter'
624 643
625 644 def setup(self):
626 645
627 646 self.xaxis = 'Range (Km)'
628 647 self.ncols = 1
629 648 self.nrows = 1
630 649 self.nplots = 1
631 650 self.ylabel = 'Intensity [dB]'
632 651 self.titles = ['Scope']
633 652 self.colorbar = False
634 653 colspan = 3
635 654 rowspan = 1
636 655
637 656 def plot_iq(self, x, y, channelIndexList, thisDatetime, wintitle):
638 657
639 658 yreal = y[channelIndexList,:].real
640 659 yimag = y[channelIndexList,:].imag
641 660 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y"))
642 661 self.xlabel = "Range (Km)"
643 662 self.ylabel = "Intensity - IQ"
644 663
645 664 self.y = yreal
646 665 self.x = x
647 666 self.xmin = min(x)
648 667 self.xmax = max(x)
649 668
650 669
651 670 self.titles[0] = title
652 671
653 672 for i,ax in enumerate(self.axes):
654 673 title = "Channel %d" %(i)
655 674 if ax.firsttime:
656 675 ax.plt_r = ax.plot(x, yreal[i,:], color='b')[0]
657 676 ax.plt_i = ax.plot(x, yimag[i,:], color='r')[0]
658 677 else:
659 678 #pass
660 679 ax.plt_r.set_data(x, yreal[i,:])
661 680 ax.plt_i.set_data(x, yimag[i,:])
662 681
663 682 def plot_power(self, x, y, channelIndexList, thisDatetime, wintitle):
664 683 y = y[channelIndexList,:] * numpy.conjugate(y[channelIndexList,:])
665 684 yreal = y.real
666 685 self.y = yreal
667 686 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y"))
668 687 self.xlabel = "Range (Km)"
669 688 self.ylabel = "Intensity"
670 689 self.xmin = min(x)
671 690 self.xmax = max(x)
672 691
673 692
674 693 self.titles[0] = title
675 694
676 695 for i,ax in enumerate(self.axes):
677 696 title = "Channel %d" %(i)
678 697
679 698 ychannel = yreal[i,:]
680 699
681 700 if ax.firsttime:
682 701 ax.plt_r = ax.plot(x, ychannel)[0]
683 702 else:
684 703 #pass
685 704 ax.plt_r.set_data(x, ychannel)
686 705
687 706
688 707 def plot(self):
689 708
690 709 if self.channels:
691 710 channels = self.channels
692 711 else:
693 712 channels = self.data.channels
694 713
695 714
696 715
697 716 thisDatetime = datetime.datetime.utcfromtimestamp(self.data.times[-1])
698 717
699 718 scope = self.data['scope']
700 719
701 720
702 721 if self.data.flagDataAsBlock:
703 722
704 723 for i in range(self.data.nProfiles):
705 724
706 725 wintitle1 = " [Profile = %d] " %i
707 726
708 727 if self.type == "power":
709 728 self.plot_power(self.data.heights,
710 729 scope[:,i,:],
711 730 channels,
712 731 thisDatetime,
713 732 wintitle1
714 733 )
715 734
716 735 if self.type == "iq":
717 736 self.plot_iq(self.data.heights,
718 737 scope[:,i,:],
719 738 channels,
720 739 thisDatetime,
721 740 wintitle1
722 741 )
723
724
725
726
727
728 742 else:
729 743 wintitle = " [Profile = %d] " %self.data.profileIndex
730 744
731 745 if self.type == "power":
732 746 self.plot_power(self.data.heights,
733 747 scope,
734 748 channels,
735 749 thisDatetime,
736 750 wintitle
737 751 )
738 752
739 753 if self.type == "iq":
740 754 self.plot_iq(self.data.heights,
741 755 scope,
742 756 channels,
743 757 thisDatetime,
744 758 wintitle
745 759 )
746
747
748 No newline at end of file
General Comments 0
You need to be logged in to leave comments. Login now