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