##// END OF EJS Templates
schain-jc
RSS
Fix publish and plots operations issue #929
Juan C. Espinoza -
r1062:8048843f4edf
parent child
Show More
Side by Side Unified
Context file:
r1062:8048843f4edf Loading diff...:

The requested changes are too big and content was truncated. Show full diff

Show file before | Show file after | Hide comments
@@ -1,1220 +1,1220
1 1 '''
2 2
3 3 $Author: murco $
4 4 $Id: JROData.py 173 2012-11-20 15:06:21Z murco $
5 5 '''
6 6
7 7 import copy
8 8 import numpy
9 9 import datetime
10 10
11 11 from jroheaderIO import SystemHeader, RadarControllerHeader
12 12 from schainpy import cSchain
13 13
14 14
15 15 def getNumpyDtype(dataTypeCode):
16 16
17 17 if dataTypeCode == 0:
18 18 numpyDtype = numpy.dtype([('real','<i1'),('imag','<i1')])
19 19 elif dataTypeCode == 1:
20 20 numpyDtype = numpy.dtype([('real','<i2'),('imag','<i2')])
21 21 elif dataTypeCode == 2:
22 22 numpyDtype = numpy.dtype([('real','<i4'),('imag','<i4')])
23 23 elif dataTypeCode == 3:
24 24 numpyDtype = numpy.dtype([('real','<i8'),('imag','<i8')])
25 25 elif dataTypeCode == 4:
26 26 numpyDtype = numpy.dtype([('real','<f4'),('imag','<f4')])
27 27 elif dataTypeCode == 5:
28 28 numpyDtype = numpy.dtype([('real','<f8'),('imag','<f8')])
29 29 else:
30 30 raise ValueError, 'dataTypeCode was not defined'
31 31
32 32 return numpyDtype
33 33
34 34 def getDataTypeCode(numpyDtype):
35 35
36 36 if numpyDtype == numpy.dtype([('real','<i1'),('imag','<i1')]):
37 37 datatype = 0
38 38 elif numpyDtype == numpy.dtype([('real','<i2'),('imag','<i2')]):
39 39 datatype = 1
40 40 elif numpyDtype == numpy.dtype([('real','<i4'),('imag','<i4')]):
41 41 datatype = 2
42 42 elif numpyDtype == numpy.dtype([('real','<i8'),('imag','<i8')]):
43 43 datatype = 3
44 44 elif numpyDtype == numpy.dtype([('real','<f4'),('imag','<f4')]):
45 45 datatype = 4
46 46 elif numpyDtype == numpy.dtype([('real','<f8'),('imag','<f8')]):
47 47 datatype = 5
48 48 else:
49 49 datatype = None
50 50
51 51 return datatype
52 52
53 53 def hildebrand_sekhon(data, navg):
54 54 """
55 55 This method is for the objective determination of the noise level in Doppler spectra. This
56 56 implementation technique is based on the fact that the standard deviation of the spectral
57 57 densities is equal to the mean spectral density for white Gaussian noise
58 58
59 59 Inputs:
60 60 Data : heights
61 61 navg : numbers of averages
62 62
63 63 Return:
64 64 -1 : any error
65 65 anoise : noise's level
66 66 """
67 67
68 68 sortdata = numpy.sort(data, axis=None)
69 69 # lenOfData = len(sortdata)
70 70 # nums_min = lenOfData*0.2
71 71 #
72 72 # if nums_min <= 5:
73 73 # nums_min = 5
74 74 #
75 75 # sump = 0.
76 76 #
77 77 # sumq = 0.
78 78 #
79 79 # j = 0
80 80 #
81 81 # cont = 1
82 82 #
83 83 # while((cont==1)and(j<lenOfData)):
84 84 #
85 85 # sump += sortdata[j]
86 86 #
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 return cSchain.hildebrand_sekhon(sortdata, navg)
104 104
105 105
106 106 class Beam:
107 107
108 108 def __init__(self):
109 109 self.codeList = []
110 110 self.azimuthList = []
111 111 self.zenithList = []
112 112
113 113 class GenericData(object):
114 114
115 115 flagNoData = True
116 116
117 117 def copy(self, inputObj=None):
118 118
119 119 if inputObj == None:
120 120 return copy.deepcopy(self)
121 121
122 122 for key in inputObj.__dict__.keys():
123 123
124 124 attribute = inputObj.__dict__[key]
125 125
126 126 #If this attribute is a tuple or list
127 127 if type(inputObj.__dict__[key]) in (tuple, list):
128 128 self.__dict__[key] = attribute[:]
129 129 continue
130 130
131 131 #If this attribute is another object or instance
132 132 if hasattr(attribute, '__dict__'):
133 133 self.__dict__[key] = attribute.copy()
134 134 continue
135 135
136 136 self.__dict__[key] = inputObj.__dict__[key]
137 137
138 138 def deepcopy(self):
139 139
140 140 return copy.deepcopy(self)
141 141
142 142 def isEmpty(self):
143 143
144 144 return self.flagNoData
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
153 153 radarControllerHeaderObj = RadarControllerHeader()
154 154
155 155 # data = None
156 156
157 157 type = None
158 158
159 159 datatype = None #dtype but in string
160 160
161 161 # dtype = None
162 162
163 163 # nChannels = None
164 164
165 165 # nHeights = None
166 166
167 167 nProfiles = None
168 168
169 169 heightList = None
170 170
171 171 channelList = None
172 172
173 173 flagDiscontinuousBlock = False
174 174
175 175 useLocalTime = False
176 176
177 177 utctime = None
178 178
179 179 timeZone = None
180 180
181 181 dstFlag = None
182 182
183 183 errorCount = None
184 184
185 185 blocksize = None
186 186
187 187 # nCode = None
188 188 #
189 189 # nBaud = None
190 190 #
191 191 # code = None
192 192
193 193 flagDecodeData = False #asumo q la data no esta decodificada
194 194
195 195 flagDeflipData = False #asumo q la data no esta sin flip
196 196
197 197 flagShiftFFT = False
198 198
199 199 # ippSeconds = None
200 200
201 201 # timeInterval = None
202 202
203 203 nCohInt = None
204 204
205 205 # noise = None
206 206
207 207 windowOfFilter = 1
208 208
209 209 #Speed of ligth
210 210 C = 3e8
211 211
212 212 frequency = 49.92e6
213 213
214 214 realtime = False
215 215
216 216 beacon_heiIndexList = None
217 217
218 218 last_block = None
219 219
220 220 blocknow = None
221 221
222 222 azimuth = None
223 223
224 224 zenith = None
225 225
226 226 beam = Beam()
227 227
228 228 profileIndex = None
229 229
230 230 def getNoise(self):
231 231
232 232 raise NotImplementedError
233 233
234 234 def getNChannels(self):
235 235
236 236 return len(self.channelList)
237 237
238 238 def getChannelIndexList(self):
239 239
240 240 return range(self.nChannels)
241 241
242 242 def getNHeights(self):
243 243
244 244 return len(self.heightList)
245 245
246 246 def getHeiRange(self, extrapoints=0):
247 247
248 248 heis = self.heightList
249 249 # deltah = self.heightList[1] - self.heightList[0]
250 250 #
251 251 # heis.append(self.heightList[-1])
252 252
253 253 return heis
254 254
255 255 def getDeltaH(self):
256 256
257 257 delta = self.heightList[1] - self.heightList[0]
258 258
259 259 return delta
260 260
261 261 def getltctime(self):
262 262
263 263 if self.useLocalTime:
264 264 return self.utctime - self.timeZone*60
265 265
266 266 return self.utctime
267 267
268 268 def getDatatime(self):
269 269
270 270 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
271 271 return datatimeValue
272 272
273 273 def getTimeRange(self):
274 274
275 275 datatime = []
276 276
277 277 datatime.append(self.ltctime)
278 278 datatime.append(self.ltctime + self.timeInterval+1)
279 279
280 280 datatime = numpy.array(datatime)
281 281
282 282 return datatime
283 283
284 284 def getFmaxTimeResponse(self):
285 285
286 286 period = (10**-6)*self.getDeltaH()/(0.15)
287 287
288 288 PRF = 1./(period * self.nCohInt)
289 289
290 290 fmax = PRF
291 291
292 292 return fmax
293 293
294 294 def getFmax(self):
295 295
296 296 PRF = 1./(self.ippSeconds * self.nCohInt)
297 297
298 298 fmax = PRF
299 299
300 300 return fmax
301 301
302 302 def getVmax(self):
303 303
304 304 _lambda = self.C/self.frequency
305 305
306 306 vmax = self.getFmax() * _lambda/2
307 307
308 308 return vmax
309 309
310 310 def get_ippSeconds(self):
311 311 '''
312 312 '''
313 313 return self.radarControllerHeaderObj.ippSeconds
314 314
315 315 def set_ippSeconds(self, ippSeconds):
316 316 '''
317 317 '''
318 318
319 319 self.radarControllerHeaderObj.ippSeconds = ippSeconds
320 320
321 321 return
322 322
323 323 def get_dtype(self):
324 324 '''
325 325 '''
326 326 return getNumpyDtype(self.datatype)
327 327
328 328 def set_dtype(self, numpyDtype):
329 329 '''
330 330 '''
331 331
332 332 self.datatype = getDataTypeCode(numpyDtype)
333 333
334 334 def get_code(self):
335 335 '''
336 336 '''
337 337 return self.radarControllerHeaderObj.code
338 338
339 339 def set_code(self, code):
340 340 '''
341 341 '''
342 342 self.radarControllerHeaderObj.code = code
343 343
344 344 return
345 345
346 346 def get_ncode(self):
347 347 '''
348 348 '''
349 349 return self.radarControllerHeaderObj.nCode
350 350
351 351 def set_ncode(self, nCode):
352 352 '''
353 353 '''
354 354 self.radarControllerHeaderObj.nCode = nCode
355 355
356 356 return
357 357
358 358 def get_nbaud(self):
359 359 '''
360 360 '''
361 361 return self.radarControllerHeaderObj.nBaud
362 362
363 363 def set_nbaud(self, nBaud):
364 364 '''
365 365 '''
366 366 self.radarControllerHeaderObj.nBaud = nBaud
367 367
368 368 return
369 369
370 370 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
371 371 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
372 372 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
373 373 #noise = property(getNoise, "I'm the 'nHeights' property.")
374 374 datatime = property(getDatatime, "I'm the 'datatime' property")
375 375 ltctime = property(getltctime, "I'm the 'ltctime' property")
376 376 ippSeconds = property(get_ippSeconds, set_ippSeconds)
377 377 dtype = property(get_dtype, set_dtype)
378 378 # timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
379 379 code = property(get_code, set_code)
380 380 nCode = property(get_ncode, set_ncode)
381 381 nBaud = property(get_nbaud, set_nbaud)
382 382
383 383 class Voltage(JROData):
384 384
385 385 #data es un numpy array de 2 dmensiones (canales, alturas)
386 386 data = None
387 387
388 388 def __init__(self):
389 389 '''
390 390 Constructor
391 391 '''
392 392
393 393 self.useLocalTime = True
394 394
395 395 self.radarControllerHeaderObj = RadarControllerHeader()
396 396
397 397 self.systemHeaderObj = SystemHeader()
398 398
399 399 self.type = "Voltage"
400 400
401 401 self.data = None
402 402
403 403 # self.dtype = None
404 404
405 405 # self.nChannels = 0
406 406
407 407 # self.nHeights = 0
408 408
409 409 self.nProfiles = None
410 410
411 411 self.heightList = None
412 412
413 413 self.channelList = None
414 414
415 415 # self.channelIndexList = None
416 416
417 417 self.flagNoData = True
418 418
419 419 self.flagDiscontinuousBlock = False
420 420
421 421 self.utctime = None
422 422
423 423 self.timeZone = None
424 424
425 425 self.dstFlag = None
426 426
427 427 self.errorCount = None
428 428
429 429 self.nCohInt = None
430 430
431 431 self.blocksize = None
432 432
433 433 self.flagDecodeData = False #asumo q la data no esta decodificada
434 434
435 435 self.flagDeflipData = False #asumo q la data no esta sin flip
436 436
437 437 self.flagShiftFFT = False
438 438
439 439 self.flagDataAsBlock = False #Asumo que la data es leida perfil a perfil
440 440
441 441 self.profileIndex = 0
442 442
443 443 def getNoisebyHildebrand(self, channel = None):
444 444 """
445 445 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
446 446
447 447 Return:
448 448 noiselevel
449 449 """
450 450
451 451 if channel != None:
452 452 data = self.data[channel]
453 453 nChannels = 1
454 454 else:
455 455 data = self.data
456 456 nChannels = self.nChannels
457 457
458 458 noise = numpy.zeros(nChannels)
459 459 power = data * numpy.conjugate(data)
460 460
461 461 for thisChannel in range(nChannels):
462 462 if nChannels == 1:
463 463 daux = power[:].real
464 464 else:
465 465 daux = power[thisChannel,:].real
466 466 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
467 467
468 468 return noise
469 469
470 470 def getNoise(self, type = 1, channel = None):
471 471
472 472 if type == 1:
473 473 noise = self.getNoisebyHildebrand(channel)
474 474
475 475 return noise
476 476
477 477 def getPower(self, channel = None):
478 478
479 479 if channel != None:
480 480 data = self.data[channel]
481 481 else:
482 482 data = self.data
483 483
484 484 power = data * numpy.conjugate(data)
485 485 powerdB = 10*numpy.log10(power.real)
486 486 powerdB = numpy.squeeze(powerdB)
487 487
488 488 return powerdB
489 489
490 490 def getTimeInterval(self):
491 491
492 492 timeInterval = self.ippSeconds * self.nCohInt
493 493
494 494 return timeInterval
495 495
496 496 noise = property(getNoise, "I'm the 'nHeights' property.")
497 497 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
498 498
499 499 class Spectra(JROData):
500 500
501 501 #data spc es un numpy array de 2 dmensiones (canales, perfiles, alturas)
502 502 data_spc = None
503 503
504 504 #data cspc es un numpy array de 2 dmensiones (canales, pares, alturas)
505 505 data_cspc = None
506 506
507 507 #data dc es un numpy array de 2 dmensiones (canales, alturas)
508 508 data_dc = None
509 509
510 510 #data power
511 511 data_pwr = None
512 512
513 513 nFFTPoints = None
514 514
515 515 # nPairs = None
516 516
517 517 pairsList = None
518 518
519 519 nIncohInt = None
520 520
521 521 wavelength = None #Necesario para cacular el rango de velocidad desde la frecuencia
522 522
523 523 nCohInt = None #se requiere para determinar el valor de timeInterval
524 524
525 525 ippFactor = None
526 526
527 527 profileIndex = 0
528 528
529 529 plotting = "spectra"
530 530
531 531 def __init__(self):
532 532 '''
533 533 Constructor
534 534 '''
535 535
536 536 self.useLocalTime = True
537 537
538 538 self.radarControllerHeaderObj = RadarControllerHeader()
539 539
540 540 self.systemHeaderObj = SystemHeader()
541 541
542 542 self.type = "Spectra"
543 543
544 544 # self.data = None
545 545
546 546 # self.dtype = None
547 547
548 548 # self.nChannels = 0
549 549
550 550 # self.nHeights = 0
551 551
552 552 self.nProfiles = None
553 553
554 554 self.heightList = None
555 555
556 556 self.channelList = None
557 557
558 558 # self.channelIndexList = None
559 559
560 560 self.pairsList = None
561 561
562 562 self.flagNoData = True
563 563
564 564 self.flagDiscontinuousBlock = False
565 565
566 566 self.utctime = None
567 567
568 568 self.nCohInt = None
569 569
570 570 self.nIncohInt = None
571 571
572 572 self.blocksize = None
573 573
574 574 self.nFFTPoints = None
575 575
576 576 self.wavelength = None
577 577
578 578 self.flagDecodeData = False #asumo q la data no esta decodificada
579 579
580 580 self.flagDeflipData = False #asumo q la data no esta sin flip
581 581
582 582 self.flagShiftFFT = False
583 583
584 584 self.ippFactor = 1
585 585
586 586 #self.noise = None
587 587
588 588 self.beacon_heiIndexList = []
589 589
590 590 self.noise_estimation = None
591 591
592 592
593 593 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
594 594 """
595 595 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
596 596
597 597 Return:
598 598 noiselevel
599 599 """
600 600
601 601 noise = numpy.zeros(self.nChannels)
602 602
603 603 for channel in range(self.nChannels):
604 604 daux = self.data_spc[channel,xmin_index:xmax_index,ymin_index:ymax_index]
605 605 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
606 606
607 607 return noise
608 608
609 609 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
610 610
611 611 if self.noise_estimation is not None:
612 612 return self.noise_estimation #this was estimated by getNoise Operation defined in jroproc_spectra.py
613 613 else:
614 614 noise = self.getNoisebyHildebrand(xmin_index, xmax_index, ymin_index, ymax_index)
615 615 return noise
616 616
617 617 def getFreqRangeTimeResponse(self, extrapoints=0):
618 618
619 619 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints*self.ippFactor)
620 620 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
621 621
622 622 return freqrange
623 623
624 624 def getAcfRange(self, extrapoints=0):
625 625
626 626 deltafreq = 10./(self.getFmax() / (self.nFFTPoints*self.ippFactor))
627 627 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
628 628
629 629 return freqrange
630 630
631 631 def getFreqRange(self, extrapoints=0):
632 632
633 633 deltafreq = self.getFmax() / (self.nFFTPoints*self.ippFactor)
634 634 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
635 635
636 636 return freqrange
637 637
638 638 def getVelRange(self, extrapoints=0):
639 639
640 640 deltav = self.getVmax() / (self.nFFTPoints*self.ippFactor)
641 641 velrange = deltav*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) #- deltav/2
642 642
643 643 return velrange
644 644
645 645 def getNPairs(self):
646 646
647 647 return len(self.pairsList)
648 648
649 649 def getPairsIndexList(self):
650 650
651 651 return range(self.nPairs)
652 652
653 653 def getNormFactor(self):
654 654
655 655 pwcode = 1
656 656
657 657 if self.flagDecodeData:
658 658 pwcode = numpy.sum(self.code[0]**2)
659 659 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
660 660 normFactor = self.nProfiles*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
661 661
662 662 return normFactor
663 663
664 664 def getFlagCspc(self):
665 665
666 666 if self.data_cspc is None:
667 667 return True
668 668
669 669 return False
670 670
671 671 def getFlagDc(self):
672 672
673 673 if self.data_dc is None:
674 674 return True
675 675
676 676 return False
677 677
678 678 def getTimeInterval(self):
679 679
680 680 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles
681 681
682 682 return timeInterval
683 683
684 684 def getPower(self):
685 685
686 686 factor = self.normFactor
687 687 z = self.data_spc/factor
688 688 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
689 689 avg = numpy.average(z, axis=1)
690 690
691 691 return 10*numpy.log10(avg)
692 692
693 693 def getCoherence(self, pairsList=None, phase=False):
694 694
695 695 z = []
696 696 if pairsList is None:
697 697 pairsIndexList = self.pairsIndexList
698 698 else:
699 699 pairsIndexList = []
700 700 for pair in pairsList:
701 701 if pair not in self.pairsList:
702 702 raise ValueError, "Pair %s is not in dataOut.pairsList" %(pair)
703 pairsIndexList.append(self.pairsList.index(pair))
703 pairsIndexList.append(self.pairsList.index(pair))
704 704 for i in range(len(pairsIndexList)):
705 705 pair = self.pairsList[pairsIndexList[i]]
706 706 ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
707 707 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
708 708 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
709 709 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
710 710 if phase:
711 711 data = numpy.arctan2(avgcoherenceComplex.imag,
712 712 avgcoherenceComplex.real)*180/numpy.pi
713 713 else:
714 714 data = numpy.abs(avgcoherenceComplex)
715 715
716 716 z.append(data)
717 717
718 718 return numpy.array(z)
719 719
720 720 def setValue(self, value):
721 721
722 722 print "This property should not be initialized"
723 723
724 724 return
725 725
726 726 nPairs = property(getNPairs, setValue, "I'm the 'nPairs' property.")
727 727 pairsIndexList = property(getPairsIndexList, setValue, "I'm the 'pairsIndexList' property.")
728 728 normFactor = property(getNormFactor, setValue, "I'm the 'getNormFactor' property.")
729 729 flag_cspc = property(getFlagCspc, setValue)
730 730 flag_dc = property(getFlagDc, setValue)
731 731 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
732 732 timeInterval = property(getTimeInterval, setValue, "I'm the 'timeInterval' property")
733 733
734 734 class SpectraHeis(Spectra):
735 735
736 736 data_spc = None
737 737
738 738 data_cspc = None
739 739
740 740 data_dc = None
741 741
742 742 nFFTPoints = None
743 743
744 744 # nPairs = None
745 745
746 746 pairsList = None
747 747
748 748 nCohInt = None
749 749
750 750 nIncohInt = None
751 751
752 752 def __init__(self):
753 753
754 754 self.radarControllerHeaderObj = RadarControllerHeader()
755 755
756 756 self.systemHeaderObj = SystemHeader()
757 757
758 758 self.type = "SpectraHeis"
759 759
760 760 # self.dtype = None
761 761
762 762 # self.nChannels = 0
763 763
764 764 # self.nHeights = 0
765 765
766 766 self.nProfiles = None
767 767
768 768 self.heightList = None
769 769
770 770 self.channelList = None
771 771
772 772 # self.channelIndexList = None
773 773
774 774 self.flagNoData = True
775 775
776 776 self.flagDiscontinuousBlock = False
777 777
778 778 # self.nPairs = 0
779 779
780 780 self.utctime = None
781 781
782 782 self.blocksize = None
783 783
784 784 self.profileIndex = 0
785 785
786 786 self.nCohInt = 1
787 787
788 788 self.nIncohInt = 1
789 789
790 790 def getNormFactor(self):
791 791 pwcode = 1
792 792 if self.flagDecodeData:
793 793 pwcode = numpy.sum(self.code[0]**2)
794 794
795 795 normFactor = self.nIncohInt*self.nCohInt*pwcode
796 796
797 797 return normFactor
798 798
799 799 def getTimeInterval(self):
800 800
801 801 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
802 802
803 803 return timeInterval
804 804
805 805 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
806 806 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
807 807
808 808 class Fits(JROData):
809 809
810 810 heightList = None
811 811
812 812 channelList = None
813 813
814 814 flagNoData = True
815 815
816 816 flagDiscontinuousBlock = False
817 817
818 818 useLocalTime = False
819 819
820 820 utctime = None
821 821
822 822 timeZone = None
823 823
824 824 # ippSeconds = None
825 825
826 826 # timeInterval = None
827 827
828 828 nCohInt = None
829 829
830 830 nIncohInt = None
831 831
832 832 noise = None
833 833
834 834 windowOfFilter = 1
835 835
836 836 #Speed of ligth
837 837 C = 3e8
838 838
839 839 frequency = 49.92e6
840 840
841 841 realtime = False
842 842
843 843
844 844 def __init__(self):
845 845
846 846 self.type = "Fits"
847 847
848 848 self.nProfiles = None
849 849
850 850 self.heightList = None
851 851
852 852 self.channelList = None
853 853
854 854 # self.channelIndexList = None
855 855
856 856 self.flagNoData = True
857 857
858 858 self.utctime = None
859 859
860 860 self.nCohInt = 1
861 861
862 862 self.nIncohInt = 1
863 863
864 864 self.useLocalTime = True
865 865
866 866 self.profileIndex = 0
867 867
868 868 # self.utctime = None
869 869 # self.timeZone = None
870 870 # self.ltctime = None
871 871 # self.timeInterval = None
872 872 # self.header = None
873 873 # self.data_header = None
874 874 # self.data = None
875 875 # self.datatime = None
876 876 # self.flagNoData = False
877 877 # self.expName = ''
878 878 # self.nChannels = None
879 879 # self.nSamples = None
880 880 # self.dataBlocksPerFile = None
881 881 # self.comments = ''
882 882 #
883 883
884 884
885 885 def getltctime(self):
886 886
887 887 if self.useLocalTime:
888 888 return self.utctime - self.timeZone*60
889 889
890 890 return self.utctime
891 891
892 892 def getDatatime(self):
893 893
894 894 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
895 895 return datatime
896 896
897 897 def getTimeRange(self):
898 898
899 899 datatime = []
900 900
901 901 datatime.append(self.ltctime)
902 902 datatime.append(self.ltctime + self.timeInterval)
903 903
904 904 datatime = numpy.array(datatime)
905 905
906 906 return datatime
907 907
908 908 def getHeiRange(self):
909 909
910 910 heis = self.heightList
911 911
912 912 return heis
913 913
914 914 def getNHeights(self):
915 915
916 916 return len(self.heightList)
917 917
918 918 def getNChannels(self):
919 919
920 920 return len(self.channelList)
921 921
922 922 def getChannelIndexList(self):
923 923
924 924 return range(self.nChannels)
925 925
926 926 def getNoise(self, type = 1):
927 927
928 928 #noise = numpy.zeros(self.nChannels)
929 929
930 930 if type == 1:
931 931 noise = self.getNoisebyHildebrand()
932 932
933 933 if type == 2:
934 934 noise = self.getNoisebySort()
935 935
936 936 if type == 3:
937 937 noise = self.getNoisebyWindow()
938 938
939 939 return noise
940 940
941 941 def getTimeInterval(self):
942 942
943 943 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
944 944
945 945 return timeInterval
946 946
947 947 datatime = property(getDatatime, "I'm the 'datatime' property")
948 948 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
949 949 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
950 950 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
951 951 noise = property(getNoise, "I'm the 'nHeights' property.")
952 952
953 953 ltctime = property(getltctime, "I'm the 'ltctime' property")
954 954 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
955 955
956 956
957 957 class Correlation(JROData):
958 958
959 959 noise = None
960 960
961 961 SNR = None
962 962
963 963 #--------------------------------------------------
964 964
965 965 mode = None
966 966
967 967 split = False
968 968
969 969 data_cf = None
970 970
971 971 lags = None
972 972
973 973 lagRange = None
974 974
975 975 pairsList = None
976 976
977 977 normFactor = None
978 978
979 979 #--------------------------------------------------
980 980
981 981 # calculateVelocity = None
982 982
983 983 nLags = None
984 984
985 985 nPairs = None
986 986
987 987 nAvg = None
988 988
989 989
990 990 def __init__(self):
991 991 '''
992 992 Constructor
993 993 '''
994 994 self.radarControllerHeaderObj = RadarControllerHeader()
995 995
996 996 self.systemHeaderObj = SystemHeader()
997 997
998 998 self.type = "Correlation"
999 999
1000 1000 self.data = None
1001 1001
1002 1002 self.dtype = None
1003 1003
1004 1004 self.nProfiles = None
1005 1005
1006 1006 self.heightList = None
1007 1007
1008 1008 self.channelList = None
1009 1009
1010 1010 self.flagNoData = True
1011 1011
1012 1012 self.flagDiscontinuousBlock = False
1013 1013
1014 1014 self.utctime = None
1015 1015
1016 1016 self.timeZone = None
1017 1017
1018 1018 self.dstFlag = None
1019 1019
1020 1020 self.errorCount = None
1021 1021
1022 1022 self.blocksize = None
1023 1023
1024 1024 self.flagDecodeData = False #asumo q la data no esta decodificada
1025 1025
1026 1026 self.flagDeflipData = False #asumo q la data no esta sin flip
1027 1027
1028 1028 self.pairsList = None
1029 1029
1030 1030 self.nPoints = None
1031 1031
1032 1032 def getPairsList(self):
1033 1033
1034 1034 return self.pairsList
1035 1035
1036 1036 def getNoise(self, mode = 2):
1037 1037
1038 1038 indR = numpy.where(self.lagR == 0)[0][0]
1039 1039 indT = numpy.where(self.lagT == 0)[0][0]
1040 1040
1041 1041 jspectra0 = self.data_corr[:,:,indR,:]
1042 1042 jspectra = copy.copy(jspectra0)
1043 1043
1044 1044 num_chan = jspectra.shape[0]
1045 1045 num_hei = jspectra.shape[2]
1046 1046
1047 1047 freq_dc = jspectra.shape[1]/2
1048 1048 ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
1049 1049
1050 1050 if ind_vel[0]<0:
1051 1051 ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
1052 1052
1053 1053 if mode == 1:
1054 1054 jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
1055 1055
1056 1056 if mode == 2:
1057 1057
1058 1058 vel = numpy.array([-2,-1,1,2])
1059 1059 xx = numpy.zeros([4,4])
1060 1060
1061 1061 for fil in range(4):
1062 1062 xx[fil,:] = vel[fil]**numpy.asarray(range(4))
1063 1063
1064 1064 xx_inv = numpy.linalg.inv(xx)
1065 1065 xx_aux = xx_inv[0,:]
1066 1066
1067 1067 for ich in range(num_chan):
1068 1068 yy = jspectra[ich,ind_vel,:]
1069 1069 jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy)
1070 1070
1071 1071 junkid = jspectra[ich,freq_dc,:]<=0
1072 1072 cjunkid = sum(junkid)
1073 1073
1074 1074 if cjunkid.any():
1075 1075 jspectra[ich,freq_dc,junkid.nonzero()] = (jspectra[ich,ind_vel[1],junkid] + jspectra[ich,ind_vel[2],junkid])/2
1076 1076
1077 1077 noise = jspectra0[:,freq_dc,:] - jspectra[:,freq_dc,:]
1078 1078
1079 1079 return noise
1080 1080
1081 1081 def getTimeInterval(self):
1082 1082
1083 1083 timeInterval = self.ippSeconds * self.nCohInt * self.nProfiles
1084 1084
1085 1085 return timeInterval
1086 1086
1087 1087 def splitFunctions(self):
1088 1088
1089 1089 pairsList = self.pairsList
1090 1090 ccf_pairs = []
1091 1091 acf_pairs = []
1092 1092 ccf_ind = []
1093 1093 acf_ind = []
1094 1094 for l in range(len(pairsList)):
1095 1095 chan0 = pairsList[l][0]
1096 1096 chan1 = pairsList[l][1]
1097 1097
1098 1098 #Obteniendo pares de Autocorrelacion
1099 1099 if chan0 == chan1:
1100 1100 acf_pairs.append(chan0)
1101 1101 acf_ind.append(l)
1102 1102 else:
1103 1103 ccf_pairs.append(pairsList[l])
1104 1104 ccf_ind.append(l)
1105 1105
1106 1106 data_acf = self.data_cf[acf_ind]
1107 1107 data_ccf = self.data_cf[ccf_ind]
1108 1108
1109 1109 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
1110 1110
1111 1111 def getNormFactor(self):
1112 1112 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
1113 1113 acf_pairs = numpy.array(acf_pairs)
1114 1114 normFactor = numpy.zeros((self.nPairs,self.nHeights))
1115 1115
1116 1116 for p in range(self.nPairs):
1117 1117 pair = self.pairsList[p]
1118 1118
1119 1119 ch0 = pair[0]
1120 1120 ch1 = pair[1]
1121 1121
1122 1122 ch0_max = numpy.max(data_acf[acf_pairs==ch0,:,:], axis=1)
1123 1123 ch1_max = numpy.max(data_acf[acf_pairs==ch1,:,:], axis=1)
1124 1124 normFactor[p,:] = numpy.sqrt(ch0_max*ch1_max)
1125 1125
1126 1126 return normFactor
1127 1127
1128 1128 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
1129 1129 normFactor = property(getNormFactor, "I'm the 'normFactor property'")
1130 1130
1131 1131 class Parameters(Spectra):
1132 1132
1133 1133 experimentInfo = None #Information about the experiment
1134 1134
1135 1135 #Information from previous data
1136 1136
1137 1137 inputUnit = None #Type of data to be processed
1138 1138
1139 1139 operation = None #Type of operation to parametrize
1140 1140
1141 1141 #normFactor = None #Normalization Factor
1142 1142
1143 1143 groupList = None #List of Pairs, Groups, etc
1144 1144
1145 1145 #Parameters
1146 1146
1147 1147 data_param = None #Parameters obtained
1148 1148
1149 1149 data_pre = None #Data Pre Parametrization
1150 1150
1151 1151 data_SNR = None #Signal to Noise Ratio
1152 1152
1153 1153 # heightRange = None #Heights
1154 1154
1155 1155 abscissaList = None #Abscissa, can be velocities, lags or time
1156 1156
1157 1157 # noise = None #Noise Potency
1158 1158
1159 1159 utctimeInit = None #Initial UTC time
1160 1160
1161 1161 paramInterval = None #Time interval to calculate Parameters in seconds
1162 1162
1163 1163 useLocalTime = True
1164 1164
1165 1165 #Fitting
1166 1166
1167 1167 data_error = None #Error of the estimation
1168 1168
1169 1169 constants = None
1170 1170
1171 1171 library = None
1172 1172
1173 1173 #Output signal
1174 1174
1175 1175 outputInterval = None #Time interval to calculate output signal in seconds
1176 1176
1177 1177 data_output = None #Out signal
1178 1178
1179 1179 nAvg = None
1180 1180
1181 1181 noise_estimation = None
1182 1182
1183 1183
1184 1184 def __init__(self):
1185 1185 '''
1186 1186 Constructor
1187 1187 '''
1188 1188 self.radarControllerHeaderObj = RadarControllerHeader()
1189 1189
1190 1190 self.systemHeaderObj = SystemHeader()
1191 1191
1192 1192 self.type = "Parameters"
1193 1193
1194 1194 def getTimeRange1(self, interval):
1195 1195
1196 1196 datatime = []
1197 1197
1198 1198 if self.useLocalTime:
1199 1199 time1 = self.utctimeInit - self.timeZone*60
1200 1200 else:
1201 1201 time1 = self.utctimeInit
1202 1202
1203 1203 datatime.append(time1)
1204 1204 datatime.append(time1 + interval)
1205 1205 datatime = numpy.array(datatime)
1206 1206
1207 1207 return datatime
1208 1208
1209 1209 def getTimeInterval(self):
1210 1210
1211 1211 if hasattr(self, 'timeInterval1'):
1212 1212 return self.timeInterval1
1213 1213 else:
1214 1214 return self.paramInterval
1215 1215
1216 1216 def getNoise(self):
1217 1217
1218 1218 return self.spc_noise
1219 1219
1220 1220 timeInterval = property(getTimeInterval)
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (1208 lines changed) Show them Hide them
@@ -1,964 +1,782
1 1
2 2 import os
3 import zmq
4 3 import time
5 import numpy
4 import glob
6 5 import datetime
7 import numpy as np
6 from multiprocessing import Process
7
8 import zmq
9 import numpy
8 10 import matplotlib
9 import glob
10 matplotlib.use('TkAgg')
11 11 import matplotlib.pyplot as plt
12 12 from mpl_toolkits.axes_grid1 import make_axes_locatable
13 from matplotlib.ticker import FuncFormatter, LinearLocator
14 from multiprocessing import Process
13 from matplotlib.ticker import FuncFormatter, LinearLocator, MultipleLocator
15 14
16 15 from schainpy.model.proc.jroproc_base import Operation
17
18 plt.ion()
16 from schainpy.utils import log
19 17
20 18 func = lambda x, pos: ('%s') %(datetime.datetime.fromtimestamp(x).strftime('%H:%M'))
21 fromtimestamp = lambda x, mintime : (datetime.datetime.utcfromtimestamp(mintime).replace(hour=(x + 5), minute=0) - d1970).total_seconds()
22 19
20 d1970 = datetime.datetime(1970, 1, 1)
23 21
24 d1970 = datetime.datetime(1970,1,1)
25 22
26 23 class PlotData(Operation, Process):
24 '''
25 Base class for Schain plotting operations
26 '''
27 27
28 28 CODE = 'Figure'
29 29 colormap = 'jro'
30 bgcolor = 'white'
30 31 CONFLATE = False
31 32 __MAXNUMX = 80
32 33 __missing = 1E30
33 34
34 35 def __init__(self, **kwargs):
35 36
36 37 Operation.__init__(self, plot=True, **kwargs)
37 38 Process.__init__(self)
38 39 self.kwargs['code'] = self.CODE
39 40 self.mp = False
40 self.dataOut = None
41 self.isConfig = False
42 self.figure = None
41 self.data = None
42 self.isConfig = False
43 self.figures = []
43 44 self.axes = []
45 self.cb_axes = []
44 46 self.localtime = kwargs.pop('localtime', True)
45 47 self.show = kwargs.get('show', True)
46 48 self.save = kwargs.get('save', False)
47 49 self.colormap = kwargs.get('colormap', self.colormap)
48 50 self.colormap_coh = kwargs.get('colormap_coh', 'jet')
49 51 self.colormap_phase = kwargs.get('colormap_phase', 'RdBu_r')
50 self.showprofile = kwargs.get('showprofile', True)
51 self.title = kwargs.get('wintitle', '')
52 self.colormaps = kwargs.get('colormaps', None)
53 self.bgcolor = kwargs.get('bgcolor', self.bgcolor)
54 self.showprofile = kwargs.get('showprofile', False)
55 self.title = kwargs.get('wintitle', self.CODE.upper())
56 self.cb_label = kwargs.get('cb_label', None)
57 self.cb_labels = kwargs.get('cb_labels', None)
52 58 self.xaxis = kwargs.get('xaxis', 'frequency')
53 59 self.zmin = kwargs.get('zmin', None)
54 60 self.zmax = kwargs.get('zmax', None)
61 self.zlimits = kwargs.get('zlimits', None)
55 62 self.xmin = kwargs.get('xmin', None)
63 if self.xmin is not None:
64 self.xmin += 5
56 65 self.xmax = kwargs.get('xmax', None)
57 66 self.xrange = kwargs.get('xrange', 24)
58 67 self.ymin = kwargs.get('ymin', None)
59 68 self.ymax = kwargs.get('ymax', None)
60 self.__MAXNUMY = kwargs.get('decimation', 5000)
61 self.throttle_value = 5
62 self.times = []
63 #self.interactive = self.kwargs['parent']
69 self.xlabel = kwargs.get('xlabel', None)
70 self.__MAXNUMY = kwargs.get('decimation', 100)
71 self.showSNR = kwargs.get('showSNR', False)
72 self.oneFigure = kwargs.get('oneFigure', True)
73 self.width = kwargs.get('width', None)
74 self.height = kwargs.get('height', None)
75 self.colorbar = kwargs.get('colorbar', True)
76 self.factors = kwargs.get('factors', [1, 1, 1, 1, 1, 1, 1, 1])
77 self.titles = ['' for __ in range(16)]
78
79 def __setup(self):
80 '''
81 Common setup for all figures, here figures and axes are created
82 '''
83
84 self.setup()
85
86 if self.width is None:
87 self.width = 8
64 88
89 self.figures = []
90 self.axes = []
91 self.cb_axes = []
92 self.pf_axes = []
93 self.cmaps = []
94
95 size = '15%' if self.ncols==1 else '30%'
96 pad = '4%' if self.ncols==1 else '8%'
97
98 if self.oneFigure:
99 if self.height is None:
100 self.height = 1.4*self.nrows + 1
101 fig = plt.figure(figsize=(self.width, self.height),
102 edgecolor='k',
103 facecolor='w')
104 self.figures.append(fig)
105 for n in range(self.nplots):
106 ax = fig.add_subplot(self.nrows, self.ncols, n+1)
107 ax.tick_params(labelsize=8)
108 ax.firsttime = True
109 self.axes.append(ax)
110 if self.showprofile:
111 cax = self.__add_axes(ax, size=size, pad=pad)
112 cax.tick_params(labelsize=8)
113 self.pf_axes.append(cax)
114 else:
115 if self.height is None:
116 self.height = 3
117 for n in range(self.nplots):
118 fig = plt.figure(figsize=(self.width, self.height),
119 edgecolor='k',
120 facecolor='w')
121 ax = fig.add_subplot(1, 1, 1)
122 ax.tick_params(labelsize=8)
123 ax.firsttime = True
124 self.figures.append(fig)
125 self.axes.append(ax)
126 if self.showprofile:
127 cax = self.__add_axes(ax, size=size, pad=pad)
128 cax.tick_params(labelsize=8)
129 self.pf_axes.append(cax)
130
131 for n in range(self.nrows):
132 if self.colormaps is not None:
133 cmap = plt.get_cmap(self.colormaps[n])
134 else:
135 cmap = plt.get_cmap(self.colormap)
136 cmap.set_bad(self.bgcolor, 1.)
137 self.cmaps.append(cmap)
138
139 def __add_axes(self, ax, size='30%', pad='8%'):
65 140 '''
66 this new parameter is created to plot data from varius channels at different figures
67 1. crear una lista de figuras donde se puedan plotear las figuras,
68 2. dar las opciones de configuracion a cada figura, estas opciones son iguales para ambas figuras
69 3. probar?
141 Add new axes to the given figure
70 142 '''
71 self.ind_plt_ch = kwargs.get('ind_plt_ch', False)
72 self.figurelist = None
143 divider = make_axes_locatable(ax)
144 nax = divider.new_horizontal(size=size, pad=pad)
145 ax.figure.add_axes(nax)
146 return nax
73 147
74 148
75 def fill_gaps(self, x_buffer, y_buffer, z_buffer):
149 def setup(self):
150 '''
151 This method should be implemented in the child class, the following
152 attributes should be set:
153
154 self.nrows: number of rows
155 self.ncols: number of cols
156 self.nplots: number of plots (channels or pairs)
157 self.ylabel: label for Y axes
158 self.titles: list of axes title
159
160 '''
161 raise(NotImplementedError, 'Implement this method in child class')
76 162
163 def fill_gaps(self, x_buffer, y_buffer, z_buffer):
164 '''
165 Create a masked array for missing data
166 '''
77 167 if x_buffer.shape[0] < 2:
78 168 return x_buffer, y_buffer, z_buffer
79 169
80 170 deltas = x_buffer[1:] - x_buffer[0:-1]
81 x_median = np.median(deltas)
171 x_median = numpy.median(deltas)
82 172
83 index = np.where(deltas > 5*x_median)
173 index = numpy.where(deltas > 5*x_median)
84 174
85 175 if len(index[0]) != 0:
86 176 z_buffer[::, index[0], ::] = self.__missing
87 z_buffer = np.ma.masked_inside(z_buffer,
177 z_buffer = numpy.ma.masked_inside(z_buffer,
88 178 0.99*self.__missing,
89 179 1.01*self.__missing)
90 180
91 181 return x_buffer, y_buffer, z_buffer
92 182
93 183 def decimate(self):
94 184
95 185 # dx = int(len(self.x)/self.__MAXNUMX) + 1
96 186 dy = int(len(self.y)/self.__MAXNUMY) + 1
97 187
98 188 # x = self.x[::dx]
99 189 x = self.x
100 190 y = self.y[::dy]
101 191 z = self.z[::, ::, ::dy]
102
192
103 193 return x, y, z
104 194
105 '''
106 JM:
107 elimana las otras imagenes generadas debido a que lso workers no llegan en orden y le pueden
108 poner otro tiempo a la figura q no necesariamente es el ultimo.
109 Solo se realiza cuando termina la imagen.
110 Problemas:
195 def format(self):
196 '''
197 Set min and max values, labels, ticks and titles
198 '''
111 199
112 File "/home/ci-81/workspace/schainv2.3/schainpy/model/graphics/jroplot_data.py", line 145, in __plot
113 for n, eachfigure in enumerate(self.figurelist):
114 TypeError: 'NoneType' object is not iterable
200 if self.xmin is None:
201 xmin = self.min_time
202 else:
203 if self.xaxis is 'time':
204 dt = datetime.datetime.fromtimestamp(self.min_time)
205 xmin = (datetime.datetime.combine(dt.date(),
206 datetime.time(int(self.xmin), 0, 0))-d1970).total_seconds()
207 else:
208 xmin = self.xmin
115 209
116 '''
117 def deleteanotherfiles(self):
118 figurenames=[]
119 if self.figurelist != None:
120 for n, eachfigure in enumerate(self.figurelist):
121 #add specific name for each channel in channelList
122 ghostfigname = os.path.join(self.save, '{}_{}_{}'.format(self.titles[n].replace(' ',''),self.CODE,
123 datetime.datetime.fromtimestamp(self.saveTime).strftime('%y%m%d')))
124 figname = os.path.join(self.save, '{}_{}_{}.png'.format(self.titles[n].replace(' ',''),self.CODE,
125 datetime.datetime.fromtimestamp(self.saveTime).strftime('%y%m%d_%H%M%S')))
126
127 for ghostfigure in glob.glob(ghostfigname+'*'): #ghostfigure will adopt all posible names of figures
128 if ghostfigure != figname:
129 os.remove(ghostfigure)
130 print 'Removing GhostFigures:' , figname
131 else :
132 '''Erasing ghost images for just on******************'''
133 ghostfigname = os.path.join(self.save, '{}_{}'.format(self.CODE,datetime.datetime.fromtimestamp(self.saveTime).strftime('%y%m%d')))
134 figname = os.path.join(self.save, '{}_{}.png'.format(self.CODE,datetime.datetime.fromtimestamp(self.saveTime).strftime('%y%m%d_%H%M%S')))
135 for ghostfigure in glob.glob(ghostfigname+'*'): #ghostfigure will adopt all posible names of figures
136 if ghostfigure != figname:
137 os.remove(ghostfigure)
138 print 'Removing GhostFigures:' , figname
210 if self.xmax is None:
211 xmax = xmin+self.xrange*60*60
212 else:
213 if self.xaxis is 'time':
214 dt = datetime.datetime.fromtimestamp(self.min_time)
215 xmax = (datetime.datetime.combine(dt.date(),
216 datetime.time(int(self.xmax), 0, 0))-d1970).total_seconds()
217 else:
218 xmax = self.xmax
219
220 ymin = self.ymin if self.ymin else numpy.nanmin(self.y)
221 ymax = self.ymax if self.ymax else numpy.nanmax(self.y)
222
223 ystep = 200 if ymax>= 800 else 100 if ymax>=400 else 50 if ymax>=200 else 20
224
225 for n, ax in enumerate(self.axes):
226 if ax.firsttime:
227 ax.set_facecolor(self.bgcolor)
228 ax.yaxis.set_major_locator(MultipleLocator(ystep))
229 if self.xaxis is 'time':
230 ax.xaxis.set_major_formatter(FuncFormatter(func))
231 ax.xaxis.set_major_locator(LinearLocator(9))
232 if self.xlabel is not None:
233 ax.set_xlabel(self.xlabel)
234 ax.set_ylabel(self.ylabel)
235 ax.firsttime = False
236 if self.showprofile:
237 self.pf_axes[n].set_ylim(ymin, ymax)
238 self.pf_axes[n].set_xlim(self.zmin, self.zmax)
239 self.pf_axes[n].set_xlabel('dB')
240 self.pf_axes[n].grid(b=True, axis='x')
241 [tick.set_visible(False) for tick in self.pf_axes[n].get_yticklabels()]
242 if self.colorbar:
243 cb = plt.colorbar(ax.plt, ax=ax, pad=0.02)
244 cb.ax.tick_params(labelsize=8)
245 if self.cb_label:
246 cb.set_label(self.cb_label, size=8)
247 elif self.cb_labels:
248 cb.set_label(self.cb_labels[n], size=8)
249
250 ax.set_title('{} - {} UTC'.format(
251 self.titles[n],
252 datetime.datetime.fromtimestamp(self.max_time).strftime('%H:%M:%S')),
253 size=8)
254 ax.set_xlim(xmin, xmax)
255 ax.set_ylim(ymin, ymax)
256
139 257
140 258 def __plot(self):
141
142 print 'plotting...{}'.format(self.CODE)
143 if self.ind_plt_ch is False : #standard
259 '''
260 '''
261 log.success('Plotting', self.name)
262
263 self.plot()
264 self.format()
265
266 for n, fig in enumerate(self.figures):
267 if self.nrows == 0 or self.nplots == 0:
268 log.warning('No data', self.name)
269 continue
144 270 if self.show:
145 self.figure.show()
146 self.plot()
147 plt.tight_layout()
148 self.figure.canvas.manager.set_window_title('{} {} - {}'.format(self.title, self.CODE.upper(),
149 datetime.datetime.fromtimestamp(self.max_time).strftime('%Y/%m/%d')))
150 else :
151 print 'len(self.figurelist): ',len(self.figurelist)
152 for n, eachfigure in enumerate(self.figurelist):
153 if self.show:
154 eachfigure.show()
155
156 self.plot()
157 eachfigure.tight_layout() # ajuste de cada subplot
158 eachfigure.canvas.manager.set_window_title('{} {} - {}'.format(self.title[n], self.CODE.upper(),
159 datetime.datetime.fromtimestamp(self.max_time).strftime('%Y/%m/%d')))
160
161 # if self.save:
162 # if self.ind_plt_ch is False : #standard
163 # figname = os.path.join(self.save, '{}_{}.png'.format(self.CODE,
164 # datetime.datetime.fromtimestamp(self.saveTime).strftime('%y%m%d_%H%M%S')))
165 # print 'Saving figure: {}'.format(figname)
166 # self.figure.savefig(figname)
167 # else :
168 # for n, eachfigure in enumerate(self.figurelist):
169 # #add specific name for each channel in channelList
170 # figname = os.path.join(self.save, '{}_{}_{}.png'.format(self.titles[n],self.CODE,
171 # datetime.datetime.fromtimestamp(self.saveTime).strftime('%y%m%d_%H%M%S')))
172 #
173 # print 'Saving figure: {}'.format(figname)
174 # eachfigure.savefig(figname)
175
176 if self.ind_plt_ch is False :
177 self.figure.canvas.draw()
178 else :
179 for eachfigure in self.figurelist:
180 eachfigure.canvas.draw()
181
182 if self.save:
183 if self.ind_plt_ch is False : #standard
184 figname = os.path.join(self.save, '{}_{}.png'.format(self.CODE,
185 datetime.datetime.fromtimestamp(self.saveTime).strftime('%y%m%d_%H%M%S')))
271 fig.show()
272
273 fig.tight_layout()
274 fig.canvas.manager.set_window_title('{} - {}'.format(self.title,
275 datetime.datetime.fromtimestamp(self.max_time).strftime('%Y/%m/%d')))
276 # fig.canvas.draw()
277
278 if self.save and self.data.ended:
279 channels = range(self.nrows)
280 if self.oneFigure:
281 label = ''
282 else:
283 label = '_{}'.format(channels[n])
284 figname = os.path.join(
285 self.save,
286 '{}{}_{}.png'.format(
287 self.CODE,
288 label,
289 datetime.datetime.fromtimestamp(self.saveTime).strftime('%y%m%d_%H%M%S')
290 )
291 )
186 292 print 'Saving figure: {}'.format(figname)
187 self.figure.savefig(figname)
188 else :
189 for n, eachfigure in enumerate(self.figurelist):
190 #add specific name for each channel in channelList
191 figname = os.path.join(self.save, '{}_{}_{}.png'.format(self.titles[n].replace(' ',''),self.CODE,
192 datetime.datetime.fromtimestamp(self.saveTime).strftime('%y%m%d_%H%M%S')))
193
194 print 'Saving figure: {}'.format(figname)
195 eachfigure.savefig(figname)
196
293 fig.savefig(figname)
197 294
198 295 def plot(self):
199
200 print 'plotting...{}'.format(self.CODE.upper())
201 return
296 '''
297 '''
298 raise(NotImplementedError, 'Implement this method in child class')
202 299
203 300 def run(self):
204 301
205 print '[Starting] {}'.format(self.name)
302 log.success('Starting', self.name)
206 303
207 304 context = zmq.Context()
208 305 receiver = context.socket(zmq.SUB)
209 306 receiver.setsockopt(zmq.SUBSCRIBE, '')
210 307 receiver.setsockopt(zmq.CONFLATE, self.CONFLATE)
211 308
212 309 if 'server' in self.kwargs['parent']:
213 310 receiver.connect('ipc:///tmp/{}.plots'.format(self.kwargs['parent']['server']))
214 311 else:
215 receiver.connect("ipc:///tmp/zmq.plots")
216
217 seconds_passed = 0
312 receiver.connect("ipc:///tmp/zmq.plots")
218 313
219 314 while True:
220 315 try:
221 self.data = receiver.recv_pyobj(flags=zmq.NOBLOCK)#flags=zmq.NOBLOCK
222 self.started = self.data['STARTED']
223 self.dataOut = self.data['dataOut']
224
225 if (len(self.times) < len(self.data['times']) and not self.started and self.data['ENDED']):
226 continue
227
228 self.times = self.data['times']
229 self.times.sort()
230 self.throttle_value = self.data['throttle']
231 self.min_time = self.times[0]
232 self.max_time = self.times[-1]
316 self.data = receiver.recv_pyobj(flags=zmq.NOBLOCK)
317
318 self.min_time = self.data.times[0]
319 self.max_time = self.data.times[-1]
233 320
234 321 if self.isConfig is False:
235 print 'setting up'
236 self.setup()
322 self.__setup()
237 323 self.isConfig = True
238 self.__plot()
239
240 if self.data['ENDED'] is True:
241 print '********GRAPHIC ENDED********'
242 self.ended = True
243 self.isConfig = False
244 self.__plot()
245 self.deleteanotherfiles() #CLPDG
246 elif seconds_passed >= self.data['throttle']:
247 print 'passed', seconds_passed
248 self.__plot()
249 seconds_passed = 0
324
325 self.__plot()
250 326
251 327 except zmq.Again as e:
252 print 'Waiting for data...'
253 plt.pause(2)
254 seconds_passed += 2
328 log.log('Waiting for data...')
329 if self.data:
330 plt.pause(self.data.throttle)
331 else:
332 time.sleep(2)
255 333
256 334 def close(self):
257 if self.dataOut:
335 if self.data:
258 336 self.__plot()
259 337
260 338
261 339 class PlotSpectraData(PlotData):
340 '''
341 Plot for Spectra data
342 '''
262 343
263 344 CODE = 'spc'
264 colormap = 'jro'
265 CONFLATE = False
345 colormap = 'jro'
266 346
267 347 def setup(self):
268
269 ncolspan = 1
270 colspan = 1
271 self.ncols = int(numpy.sqrt(self.dataOut.nChannels)+0.9)
272 self.nrows = int(self.dataOut.nChannels*1./self.ncols + 0.9)
273 self.width = 3.6*self.ncols
274 self.height = 3.2*self.nrows
275 if self.showprofile:
276 ncolspan = 3
277 colspan = 2
278 self.width += 1.2*self.ncols
348 self.nplots = len(self.data.channels)
349 self.ncols = int(numpy.sqrt(self.nplots)+ 0.9)
350 self.nrows = int((1.0*self.nplots/self.ncols) + 0.9)
351 self.width = 3.4*self.ncols
352 self.height = 3*self.nrows
353 self.cb_label = 'dB'
354 if self.showprofile:
355 self.width += 0.8*self.ncols
279 356
280 357 self.ylabel = 'Range [Km]'
281 self.titles = ['Channel {}'.format(x) for x in self.dataOut.channelList]
282
283 if self.figure is None:
284 self.figure = plt.figure(figsize=(self.width, self.height),
285 edgecolor='k',
286 facecolor='w')
287 else:
288 self.figure.clf()
289
290 n = 0
291 for y in range(self.nrows):
292 for x in range(self.ncols):
293 if n >= self.dataOut.nChannels:
294 break
295 ax = plt.subplot2grid((self.nrows, self.ncols*ncolspan), (y, x*ncolspan), 1, colspan)
296 if self.showprofile:
297 ax.ax_profile = plt.subplot2grid((self.nrows, self.ncols*ncolspan), (y, x*ncolspan+colspan), 1, 1)
298
299 ax.firsttime = True
300 self.axes.append(ax)
301 n += 1
302 358
303 359 def plot(self):
304
305 360 if self.xaxis == "frequency":
306 x = self.dataOut.getFreqRange(1)/1000.
307 xlabel = "Frequency (kHz)"
361 x = self.data.xrange[0]
362 self.xlabel = "Frequency (kHz)"
308 363 elif self.xaxis == "time":
309 x = self.dataOut.getAcfRange(1)
310 xlabel = "Time (ms)"
364 x = self.data.xrange[1]
365 self.xlabel = "Time (ms)"
311 366 else:
312 x = self.dataOut.getVelRange(1)
313 xlabel = "Velocity (m/s)"
367 x = self.data.xrange[2]
368 self.xlabel = "Velocity (m/s)"
369
370 if self.CODE == 'spc_mean':
371 x = self.data.xrange[2]
372 self.xlabel = "Velocity (m/s)"
314 373
315 y = self.dataOut.getHeiRange()
316 z = self.data[self.CODE]
374 self.titles = []
317 375
376 y = self.data.heights
377 self.y = y
378 z = self.data['spc']
379
318 380 for n, ax in enumerate(self.axes):
381 noise = self.data['noise'][n][-1]
382 if self.CODE == 'spc_mean':
383 mean = self.data['mean'][n][-1]
319 384 if ax.firsttime:
320 self.xmax = self.xmax if self.xmax else np.nanmax(x)
385 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
321 386 self.xmin = self.xmin if self.xmin else -self.xmax
322 self.ymin = self.ymin if self.ymin else np.nanmin(y)
323 self.ymax = self.ymax if self.ymax else np.nanmax(y)
324 self.zmin = self.zmin if self.zmin else np.nanmin(z)
325 self.zmax = self.zmax if self.zmax else np.nanmax(z)
326 ax.plot = ax.pcolormesh(x, y, z[n].T,
327 vmin=self.zmin,
328 vmax=self.zmax,
329 cmap=plt.get_cmap(self.colormap)
330 )
331 divider = make_axes_locatable(ax)
332 cax = divider.new_horizontal(size='3%', pad=0.05)
333 self.figure.add_axes(cax)
334 plt.colorbar(ax.plot, cax)
335
336 ax.set_xlim(self.xmin, self.xmax)
337 ax.set_ylim(self.ymin, self.ymax)
338
339 ax.set_ylabel(self.ylabel)
340 ax.set_xlabel(xlabel)
341
342 ax.firsttime = False
387 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
388 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
389 ax.plt = ax.pcolormesh(x, y, z[n].T,
390 vmin=self.zmin,
391 vmax=self.zmax,
392 cmap=plt.get_cmap(self.colormap)
393 )
343 394
344 395 if self.showprofile:
345 ax.plot_profile= ax.ax_profile.plot(self.data['rti'][self.max_time][n], y)[0]
346 ax.ax_profile.set_xlim(self.zmin, self.zmax)
347 ax.ax_profile.set_ylim(self.ymin, self.ymax)
348 ax.ax_profile.set_xlabel('dB')
349 ax.ax_profile.grid(b=True, axis='x')
350 ax.plot_noise = ax.ax_profile.plot(numpy.repeat(self.data['noise'][self.max_time][n], len(y)), y,
351 color="k", linestyle="dashed", lw=2)[0]
352 [tick.set_visible(False) for tick in ax.ax_profile.get_yticklabels()]
396 ax.plt_profile= self.pf_axes[n].plot(self.data['rti'][n][-1], y)[0]
397 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
398 color="k", linestyle="dashed", lw=1)[0]
399 if self.CODE == 'spc_mean':
400 ax.plt_mean = ax.plot(mean, y, color='k')[0]
353 401 else:
354 ax.plot.set_array(z[n].T.ravel())
402 ax.plt.set_array(z[n].T.ravel())
355 403 if self.showprofile:
356 ax.plot_profile.set_data(self.data['rti'][self.max_time][n], y)
357 ax.plot_noise.set_data(numpy.repeat(self.data['noise'][self.max_time][n], len(y)), y)
404 ax.plt_profile.set_data(self.data['rti'][n][-1], y)
405 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
406 if self.CODE == 'spc_mean':
407 ax.plt_mean.set_data(mean, y)
358 408
359 ax.set_title('{} - Noise: {:.2f} dB'.format(self.titles[n], self.data['noise'][self.max_time][n]),
360 size=8)
409 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
361 410 self.saveTime = self.max_time
362 411
363 412
364 413 class PlotCrossSpectraData(PlotData):
365 414
366 415 CODE = 'cspc'
367 416 zmin_coh = None
368 417 zmax_coh = None
369 418 zmin_phase = None
370 zmax_phase = None
371 CONFLATE = False
419 zmax_phase = None
372 420
373 421 def setup(self):
374 422
375 ncolspan = 1
376 colspan = 1
377 self.ncols = 2
378 self.nrows = self.dataOut.nPairs
379 self.width = 3.6*self.ncols
380 self.height = 3.2*self.nrows
381
423 self.ncols = 4
424 self.nrows = len(self.data.pairs)
425 self.nplots = self.nrows*4
426 self.width = 3.4*self.ncols
427 self.height = 3*self.nrows
382 428 self.ylabel = 'Range [Km]'
383 self.titles = ['Channel {}'.format(x) for x in self.dataOut.channelList]
384
385 if self.figure is None:
386 self.figure = plt.figure(figsize=(self.width, self.height),
387 edgecolor='k',
388 facecolor='w')
389 else:
390 self.figure.clf()
391
392 for y in range(self.nrows):
393 for x in range(self.ncols):
394 ax = plt.subplot2grid((self.nrows, self.ncols), (y, x), 1, 1)
395 ax.firsttime = True
396 self.axes.append(ax)
429 self.showprofile = False
397 430
398 431 def plot(self):
399 432
400 433 if self.xaxis == "frequency":
401 x = self.dataOut.getFreqRange(1)/1000.
402 xlabel = "Frequency (kHz)"
434 x = self.data.xrange[0]
435 self.xlabel = "Frequency (kHz)"
403 436 elif self.xaxis == "time":
404 x = self.dataOut.getAcfRange(1)
405 xlabel = "Time (ms)"
437 x = self.data.xrange[1]
438 self.xlabel = "Time (ms)"
406 439 else:
407 x = self.dataOut.getVelRange(1)
408 xlabel = "Velocity (m/s)"
440 x = self.data.xrange[2]
441 self.xlabel = "Velocity (m/s)"
442
443 self.titles = []
409 444
410 y = self.dataOut.getHeiRange()
411 z_coh = self.data['cspc_coh']
412 z_phase = self.data['cspc_phase']
445 y = self.data.heights
446 self.y = y
447 spc = self.data['spc']
448 cspc = self.data['cspc']
413 449
414 450 for n in range(self.nrows):
415 ax = self.axes[2*n]
416 ax1 = self.axes[2*n+1]
451 noise = self.data['noise'][n][-1]
452 pair = self.data.pairs[n]
453 ax = self.axes[4*n]
454 ax3 = self.axes[4*n+3]
417 455 if ax.firsttime:
418 self.xmax = self.xmax if self.xmax else np.nanmax(x)
456 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
419 457 self.xmin = self.xmin if self.xmin else -self.xmax
420 self.ymin = self.ymin if self.ymin else np.nanmin(y)
421 self.ymax = self.ymax if self.ymax else np.nanmax(y)
422 self.zmin_coh = self.zmin_coh if self.zmin_coh else 0.0
423 self.zmax_coh = self.zmax_coh if self.zmax_coh else 1.0
424 self.zmin_phase = self.zmin_phase if self.zmin_phase else -180
425 self.zmax_phase = self.zmax_phase if self.zmax_phase else 180
426
427 ax.plot = ax.pcolormesh(x, y, z_coh[n].T,
428 vmin=self.zmin_coh,
429 vmax=self.zmax_coh,
430 cmap=plt.get_cmap(self.colormap_coh)
431 )
432 divider = make_axes_locatable(ax)
433 cax = divider.new_horizontal(size='3%', pad=0.05)
434 self.figure.add_axes(cax)
435 plt.colorbar(ax.plot, cax)
436
437 ax.set_xlim(self.xmin, self.xmax)
438 ax.set_ylim(self.ymin, self.ymax)
439
440 ax.set_ylabel(self.ylabel)
441 ax.set_xlabel(xlabel)
442 ax.firsttime = False
443
444 ax1.plot = ax1.pcolormesh(x, y, z_phase[n].T,
445 vmin=self.zmin_phase,
446 vmax=self.zmax_phase,
447 cmap=plt.get_cmap(self.colormap_phase)
448 )
449 divider = make_axes_locatable(ax1)
450 cax = divider.new_horizontal(size='3%', pad=0.05)
451 self.figure.add_axes(cax)
452 plt.colorbar(ax1.plot, cax)
453
454 ax1.set_xlim(self.xmin, self.xmax)
455 ax1.set_ylim(self.ymin, self.ymax)
456
457 ax1.set_ylabel(self.ylabel)
458 ax1.set_xlabel(xlabel)
459 ax1.firsttime = False
458 self.zmin = self.zmin if self.zmin else numpy.nanmin(spc)
459 self.zmax = self.zmax if self.zmax else numpy.nanmax(spc)
460 ax.plt = ax.pcolormesh(x, y, spc[pair[0]].T,
461 vmin=self.zmin,
462 vmax=self.zmax,
463 cmap=plt.get_cmap(self.colormap)
464 )
460 465 else:
461 ax.plot.set_array(z_coh[n].T.ravel())
462 ax1.plot.set_array(z_phase[n].T.ravel())
463
464 ax.set_title('Coherence Ch{} * Ch{}'.format(self.dataOut.pairsList[n][0], self.dataOut.pairsList[n][1]), size=8)
465 ax1.set_title('Phase Ch{} * Ch{}'.format(self.dataOut.pairsList[n][0], self.dataOut.pairsList[n][1]), size=8)
466 self.saveTime = self.max_time
467
466 ax.plt.set_array(spc[pair[0]].T.ravel())
467 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
468 468
469 class PlotSpectraMeanData(PlotSpectraData):
470
471 CODE = 'spc_mean'
472 colormap = 'jet'
473
474 def plot(self):
475
476 if self.xaxis == "frequency":
477 x = self.dataOut.getFreqRange(1)/1000.
478 xlabel = "Frequency (kHz)"
479 elif self.xaxis == "time":
480 x = self.dataOut.getAcfRange(1)
481 xlabel = "Time (ms)"
482 else:
483 x = self.dataOut.getVelRange(1)
484 xlabel = "Velocity (m/s)"
485
486 y = self.dataOut.getHeiRange()
487 z = self.data['spc']
488 mean = self.data['mean'][self.max_time]
489
490 for n, ax in enumerate(self.axes):
491
492 if ax.firsttime:
493 self.xmax = self.xmax if self.xmax else np.nanmax(x)
494 self.xmin = self.xmin if self.xmin else -self.xmax
495 self.ymin = self.ymin if self.ymin else np.nanmin(y)
496 self.ymax = self.ymax if self.ymax else np.nanmax(y)
497 self.zmin = self.zmin if self.zmin else np.nanmin(z)
498 self.zmax = self.zmax if self.zmax else np.nanmax(z)
499 ax.plt = ax.pcolormesh(x, y, z[n].T,
469 ax = self.axes[4*n+1]
470 if ax.firsttime:
471 ax.plt = ax.pcolormesh(x, y, spc[pair[1]].T,
500 472 vmin=self.zmin,
501 473 vmax=self.zmax,
502 474 cmap=plt.get_cmap(self.colormap)
503 475 )
504 ax.plt_dop = ax.plot(mean[n], y,
505 color='k')[0]
506
507 divider = make_axes_locatable(ax)
508 cax = divider.new_horizontal(size='3%', pad=0.05)
509 self.figure.add_axes(cax)
510 plt.colorbar(ax.plt, cax)
511
512 ax.set_xlim(self.xmin, self.xmax)
513 ax.set_ylim(self.ymin, self.ymax)
514
515 ax.set_ylabel(self.ylabel)
516 ax.set_xlabel(xlabel)
517
518 ax.firsttime = False
519
520 if self.showprofile:
521 ax.plt_profile= ax.ax_profile.plot(self.data['rti'][self.max_time][n], y)[0]
522 ax.ax_profile.set_xlim(self.zmin, self.zmax)
523 ax.ax_profile.set_ylim(self.ymin, self.ymax)
524 ax.ax_profile.set_xlabel('dB')
525 ax.ax_profile.grid(b=True, axis='x')
526 ax.plt_noise = ax.ax_profile.plot(numpy.repeat(self.data['noise'][self.max_time][n], len(y)), y,
527 color="k", linestyle="dashed", lw=2)[0]
528 [tick.set_visible(False) for tick in ax.ax_profile.get_yticklabels()]
529 476 else:
530 ax.plt.set_array(z[n].T.ravel())
531 ax.plt_dop.set_data(mean[n], y)
532 if self.showprofile:
533 ax.plt_profile.set_data(self.data['rti'][self.max_time][n], y)
534 ax.plt_noise.set_data(numpy.repeat(self.data['noise'][self.max_time][n], len(y)), y)
477 ax.plt.set_array(spc[pair[1]].T.ravel())
478 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
479
480 out = cspc[n]/numpy.sqrt(spc[pair[0]]*spc[pair[1]])
481 coh = numpy.abs(out)
482 phase = numpy.arctan2(out.imag, out.real)*180/numpy.pi
483
484 ax = self.axes[4*n+2]
485 if ax.firsttime:
486 ax.plt = ax.pcolormesh(x, y, coh.T,
487 vmin=0,
488 vmax=1,
489 cmap=plt.get_cmap(self.colormap_coh)
490 )
491 else:
492 ax.plt.set_array(coh.T.ravel())
493 self.titles.append('Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
535 494
536 ax.set_title('{} - Noise: {:.2f} dB'.format(self.titles[n], self.data['noise'][self.max_time][n]),
537 size=8)
495 ax = self.axes[4*n+3]
496 if ax.firsttime:
497 ax.plt = ax.pcolormesh(x, y, phase.T,
498 vmin=-180,
499 vmax=180,
500 cmap=plt.get_cmap(self.colormap_phase)
501 )
502 else:
503 ax.plt.set_array(phase.T.ravel())
504 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
505
538 506 self.saveTime = self.max_time
539 507
540 508
509 class PlotSpectraMeanData(PlotSpectraData):
510 '''
511 Plot for Spectra and Mean
512 '''
513 CODE = 'spc_mean'
514 colormap = 'jro'
515
516
541 517 class PlotRTIData(PlotData):
518 '''
519 Plot for RTI data
520 '''
542 521
543 522 CODE = 'rti'
544 523 colormap = 'jro'
545 524
546 525 def setup(self):
547 self.ncols = 1
548 self.nrows = self.dataOut.nChannels
549 self.width = 10
550 #TODO : arreglar la altura de la figura, esta hardcodeada.
551 #Se arreglo, testear!
552 if self.ind_plt_ch:
553 self.height = 3.2#*self.nrows if self.nrows<6 else 12
554 else:
555 self.height = 2.2*self.nrows if self.nrows<6 else 12
556
557 '''
558 if self.nrows==1:
559 self.height += 1
560 '''
526 self.xaxis = 'time'
527 self.ncols = 1
528 self.nrows = len(self.data.channels)
529 self.nplots = len(self.data.channels)
561 530 self.ylabel = 'Range [Km]'
562 self.titles = ['Channel {}'.format(x) for x in self.dataOut.channelList]
563
564 '''
565 Logica:
566 1) Si la variable ind_plt_ch es True, va a crear mas de 1 figura
567 2) guardamos "Figures" en una lista y "axes" en otra, quizas se deberia guardar el
568 axis dentro de "Figures" como un diccionario.
569 '''
570 if self.ind_plt_ch is False: #standard mode
571
572 if self.figure is None: #solo para la priemra vez
573 self.figure = plt.figure(figsize=(self.width, self.height),
574 edgecolor='k',
575 facecolor='w')
576 else:
577 self.figure.clf()
578 self.axes = []
579
580
581 for n in range(self.nrows):
582 ax = self.figure.add_subplot(self.nrows, self.ncols, n+1)
583 #ax = self.figure(n+1)
584 ax.firsttime = True
585 self.axes.append(ax)
586
587 else : #append one figure foreach channel in channelList
588 if self.figurelist == None:
589 self.figurelist = []
590 for n in range(self.nrows):
591 self.figure = plt.figure(figsize=(self.width, self.height),
592 edgecolor='k',
593 facecolor='w')
594 #add always one subplot
595 self.figurelist.append(self.figure)
596
597 else : # cada dia nuevo limpia el axes, pero mantiene el figure
598 for eachfigure in self.figurelist:
599 eachfigure.clf() # eliminaria todas las figuras de la lista?
600 self.axes = []
601
602 for eachfigure in self.figurelist:
603 ax = eachfigure.add_subplot(1,1,1) #solo 1 axis por figura
604 #ax = self.figure(n+1)
605 ax.firsttime = True
606 #Cada figura tiene un distinto puntero
607 self.axes.append(ax)
608 #plt.close(eachfigure)
609
531 self.cb_label = 'dB'
532 self.titles = ['{} Channel {}'.format(self.CODE.upper(), x) for x in range(self.nrows)]
610 533
611 534 def plot(self):
535 self.x = self.data.times
536 self.y = self.data.heights
537 self.z = self.data[self.CODE]
538 self.z = numpy.ma.masked_invalid(self.z)
612 539
613 if self.ind_plt_ch is False: #standard mode
614 self.x = np.array(self.times)
615 self.y = self.dataOut.getHeiRange()
616 self.z = []
617
618 for ch in range(self.nrows):
619 self.z.append([self.data[self.CODE][t][ch] for t in self.times])
620
621 self.z = np.array(self.z)
622 for n, ax in enumerate(self.axes):
623 x, y, z = self.fill_gaps(*self.decimate())
624 if self.xmin is None:
625 xmin = self.min_time
626 else:
627 xmin = fromtimestamp(int(self.xmin), self.min_time)
628 if self.xmax is None:
629 xmax = xmin + self.xrange*60*60
630 else:
631 xmax = xmin + (self.xmax - self.xmin) * 60 * 60
632 self.zmin = self.zmin if self.zmin else np.min(self.z)
633 self.zmax = self.zmax if self.zmax else np.max(self.z)
634 if ax.firsttime:
635 self.ymin = self.ymin if self.ymin else np.nanmin(self.y)
636 self.ymax = self.ymax if self.ymax else np.nanmax(self.y)
637 plot = ax.pcolormesh(x, y, z[n].T,
638 vmin=self.zmin,
639 vmax=self.zmax,
640 cmap=plt.get_cmap(self.colormap)
641 )
642 divider = make_axes_locatable(ax)
643 cax = divider.new_horizontal(size='2%', pad=0.05)
644 self.figure.add_axes(cax)
645 plt.colorbar(plot, cax)
646 ax.set_ylim(self.ymin, self.ymax)
647 ax.xaxis.set_major_formatter(FuncFormatter(func))
648 ax.xaxis.set_major_locator(LinearLocator(6))
649 ax.set_ylabel(self.ylabel)
650 # if self.xmin is None:
651 # xmin = self.min_time
652 # else:
653 # xmin = (datetime.datetime.combine(self.dataOut.datatime.date(),
654 # datetime.time(self.xmin, 0, 0))-d1970).total_seconds()
655
656 ax.set_xlim(xmin, xmax)
657 ax.firsttime = False
658 else:
659 ax.collections.remove(ax.collections[0])
660 ax.set_xlim(xmin, xmax)
661 plot = ax.pcolormesh(x, y, z[n].T,
662 vmin=self.zmin,
663 vmax=self.zmax,
664 cmap=plt.get_cmap(self.colormap)
665 )
666 ax.set_title('{} {}'.format(self.titles[n],
667 datetime.datetime.fromtimestamp(self.max_time).strftime('%y/%m/%d %H:%M:%S')),
668 size=8)
669
670 self.saveTime = self.min_time
671 else :
672 self.x = np.array(self.times)
673 self.y = self.dataOut.getHeiRange()
674 self.z = []
675
676 for ch in range(self.nrows):
677 self.z.append([self.data[self.CODE][t][ch] for t in self.times])
678
679 self.z = np.array(self.z)
680 for n, eachfigure in enumerate(self.figurelist): #estaba ax in axes
681
682 x, y, z = self.fill_gaps(*self.decimate())
683 xmin = self.min_time
684 xmax = xmin+self.xrange*60*60
685 self.zmin = self.zmin if self.zmin else np.min(self.z)
686 self.zmax = self.zmax if self.zmax else np.max(self.z)
687 if self.axes[n].firsttime:
688 self.ymin = self.ymin if self.ymin else np.nanmin(self.y)
689 self.ymax = self.ymax if self.ymax else np.nanmax(self.y)
690 plot = self.axes[n].pcolormesh(x, y, z[n].T,
691 vmin=self.zmin,
692 vmax=self.zmax,
693 cmap=plt.get_cmap(self.colormap)
694 )
695 divider = make_axes_locatable(self.axes[n])
696 cax = divider.new_horizontal(size='2%', pad=0.05)
697 eachfigure.add_axes(cax)
698 #self.figure2.add_axes(cax)
699 plt.colorbar(plot, cax)
700 self.axes[n].set_ylim(self.ymin, self.ymax)
701
702 self.axes[n].xaxis.set_major_formatter(FuncFormatter(func))
703 self.axes[n].xaxis.set_major_locator(LinearLocator(6))
704
705 self.axes[n].set_ylabel(self.ylabel)
706
707 if self.xmin is None:
708 xmin = self.min_time
709 else:
710 xmin = (datetime.datetime.combine(self.dataOut.datatime.date(),
711 datetime.time(self.xmin, 0, 0))-d1970).total_seconds()
712
713 self.axes[n].set_xlim(xmin, xmax)
714 self.axes[n].firsttime = False
715 else:
716 self.axes[n].collections.remove(self.axes[n].collections[0])
717 self.axes[n].set_xlim(xmin, xmax)
718 plot = self.axes[n].pcolormesh(x, y, z[n].T,
719 vmin=self.zmin,
720 vmax=self.zmax,
721 cmap=plt.get_cmap(self.colormap)
722 )
723 self.axes[n].set_title('{} {}'.format(self.titles[n],
724 datetime.datetime.fromtimestamp(self.max_time).strftime('%y/%m/%d %H:%M:%S')),
725 size=8)
540 for n, ax in enumerate(self.axes):
541 x, y, z = self.fill_gaps(*self.decimate())
542 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
543 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
544 if ax.firsttime:
545 ax.plt = ax.pcolormesh(x, y, z[n].T,
546 vmin=self.zmin,
547 vmax=self.zmax,
548 cmap=plt.get_cmap(self.colormap)
549 )
550 if self.showprofile:
551 ax.plot_profile= self.pf_axes[n].plot(self.data['rti'][n][-1], self.y)[0]
552 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(self.data['noise'][n][-1], len(self.y)), self.y,
553 color="k", linestyle="dashed", lw=1)[0]
554 else:
555 ax.collections.remove(ax.collections[0])
556 ax.plt = ax.pcolormesh(x, y, z[n].T,
557 vmin=self.zmin,
558 vmax=self.zmax,
559 cmap=plt.get_cmap(self.colormap)
560 )
561 if self.showprofile:
562 ax.plot_profile.set_data(self.data['rti'][n][-1], self.y)
563 ax.plot_noise.set_data(numpy.repeat(self.data['noise'][n][-1], len(self.y)), self.y)
726 564
727 self.saveTime = self.min_time
565 self.saveTime = self.min_time
728 566
729 567
730 568 class PlotCOHData(PlotRTIData):
569 '''
570 Plot for Coherence data
571 '''
731 572
732 573 CODE = 'coh'
733 574
734 575 def setup(self):
735
576 self.xaxis = 'time'
736 577 self.ncols = 1
737 self.nrows = self.dataOut.nPairs
738 self.width = 10
739 self.height = 2.2*self.nrows if self.nrows<6 else 12
740 self.ind_plt_ch = False #just for coherence and phase
741 if self.nrows==1:
742 self.height += 1
743 self.ylabel = 'Range [Km]'
744 self.titles = ['{} Ch{} * Ch{}'.format(self.CODE.upper(), x[0], x[1]) for x in self.dataOut.pairsList]
745
746 if self.figure is None:
747 self.figure = plt.figure(figsize=(self.width, self.height),
748 edgecolor='k',
749 facecolor='w')
578 self.nrows = len(self.data.pairs)
579 self.nplots = len(self.data.pairs)
580 self.ylabel = 'Range [Km]'
581 if self.CODE == 'coh':
582 self.cb_label = ''
583 self.titles = ['Coherence Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
750 584 else:
751 self.figure.clf()
752 self.axes = []
585 self.cb_label = 'Degrees'
586 self.titles = ['Phase Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
753 587
754 for n in range(self.nrows):
755 ax = self.figure.add_subplot(self.nrows, self.ncols, n+1)
756 ax.firsttime = True
757 self.axes.append(ax)
588
589 class PlotPHASEData(PlotCOHData):
590 '''
591 Plot for Phase map data
592 '''
593
594 CODE = 'phase'
595 colormap = 'seismic'
758 596
759 597
760 598 class PlotNoiseData(PlotData):
599 '''
600 Plot for noise
601 '''
602
761 603 CODE = 'noise'
762 604
763 605 def setup(self):
764
606 self.xaxis = 'time'
765 607 self.ncols = 1
766 608 self.nrows = 1
767 self.width = 10
768 self.height = 3.2
609 self.nplots = 1
769 610 self.ylabel = 'Intensity [dB]'
770 611 self.titles = ['Noise']
771
772 if self.figure is None:
773 self.figure = plt.figure(figsize=(self.width, self.height),
774 edgecolor='k',
775 facecolor='w')
776 else:
777 self.figure.clf()
778 self.axes = []
779
780 self.ax = self.figure.add_subplot(self.nrows, self.ncols, 1)
781 self.ax.firsttime = True
612 self.colorbar = False
782 613
783 614 def plot(self):
784 615
785 x = self.times
616 x = self.data.times
786 617 xmin = self.min_time
787 618 xmax = xmin+self.xrange*60*60
788 if self.ax.firsttime:
789 for ch in self.dataOut.channelList:
790 y = [self.data[self.CODE][t][ch] for t in self.times]
791 self.ax.plot(x, y, lw=1, label='Ch{}'.format(ch))
792 self.ax.firsttime = False
793 self.ax.xaxis.set_major_formatter(FuncFormatter(func))
794 self.ax.xaxis.set_major_locator(LinearLocator(6))
795 self.ax.set_ylabel(self.ylabel)
619 Y = self.data[self.CODE]
620
621 if self.axes[0].firsttime:
622 for ch in self.data.channels:
623 y = Y[ch]
624 self.axes[0].plot(x, y, lw=1, label='Ch{}'.format(ch))
796 625 plt.legend()
797 626 else:
798 for ch in self.dataOut.channelList:
799 y = [self.data[self.CODE][t][ch] for t in self.times]
800 self.ax.lines[ch].set_data(x, y)
801
802 self.ax.set_xlim(xmin, xmax)
803 self.ax.set_ylim(min(y)-5, max(y)+5)
627 for ch in self.data.channels:
628 y = Y[ch]
629 self.axes[0].lines[ch].set_data(x, y)
630
631 self.ymin = numpy.nanmin(Y) - 5
632 self.ymax = numpy.nanmax(Y) + 5
804 633 self.saveTime = self.min_time
805 634
806 635
807 class PlotWindProfilerData(PlotRTIData):
808
809 CODE = 'wind'
810 colormap = 'seismic'
811
812 def setup(self):
813 self.ncols = 1
814 self.nrows = self.dataOut.data_output.shape[0]
815 self.width = 10
816 self.height = 2.2*self.nrows
817 self.ylabel = 'Height [Km]'
818 self.titles = ['Zonal Wind' ,'Meridional Wind', 'Vertical Wind']
819 self.clabels = ['Velocity (m/s)','Velocity (m/s)','Velocity (cm/s)']
820 self.windFactor = [1, 1, 100]
821
822 if self.figure is None:
823 self.figure = plt.figure(figsize=(self.width, self.height),
824 edgecolor='k',
825 facecolor='w')
826 else:
827 self.figure.clf()
828 self.axes = []
829
830 for n in range(self.nrows):
831 ax = self.figure.add_subplot(self.nrows, self.ncols, n+1)
832 ax.firsttime = True
833 self.axes.append(ax)
834
835 def plot(self):
836
837 self.x = np.array(self.times)
838 self.y = self.dataOut.heightList
839 self.z = []
840
841 for ch in range(self.nrows):
842 self.z.append([self.data['output'][t][ch] for t in self.times])
843
844 self.z = np.array(self.z)
845 self.z = numpy.ma.masked_invalid(self.z)
846
847 cmap=plt.get_cmap(self.colormap)
848 cmap.set_bad('black', 1.)
849
850 for n, ax in enumerate(self.axes):
851 x, y, z = self.fill_gaps(*self.decimate())
852 xmin = self.min_time
853 xmax = xmin+self.xrange*60*60
854 if ax.firsttime:
855 self.ymin = self.ymin if self.ymin else np.nanmin(self.y)
856 self.ymax = self.ymax if self.ymax else np.nanmax(self.y)
857 self.zmax = self.zmax if self.zmax else numpy.nanmax(abs(self.z[:-1, :]))
858 self.zmin = self.zmin if self.zmin else -self.zmax
859
860 plot = ax.pcolormesh(x, y, z[n].T*self.windFactor[n],
861 vmin=self.zmin,
862 vmax=self.zmax,
863 cmap=cmap
864 )
865 divider = make_axes_locatable(ax)
866 cax = divider.new_horizontal(size='2%', pad=0.05)
867 self.figure.add_axes(cax)
868 cb = plt.colorbar(plot, cax)
869 cb.set_label(self.clabels[n])
870 ax.set_ylim(self.ymin, self.ymax)
871
872 ax.xaxis.set_major_formatter(FuncFormatter(func))
873 ax.xaxis.set_major_locator(LinearLocator(6))
874
875 ax.set_ylabel(self.ylabel)
876
877 ax.set_xlim(xmin, xmax)
878 ax.firsttime = False
879 else:
880 ax.collections.remove(ax.collections[0])
881 ax.set_xlim(xmin, xmax)
882 plot = ax.pcolormesh(x, y, z[n].T*self.windFactor[n],
883 vmin=self.zmin,
884 vmax=self.zmax,
885 cmap=plt.get_cmap(self.colormap)
886 )
887 ax.set_title('{} {}'.format(self.titles[n],
888 datetime.datetime.fromtimestamp(self.max_time).strftime('%y/%m/%d %H:%M:%S')),
889 size=8)
890
891 self.saveTime = self.min_time
892
893
894 636 class PlotSNRData(PlotRTIData):
637 '''
638 Plot for SNR Data
639 '''
640
895 641 CODE = 'snr'
896 642 colormap = 'jet'
897 643
644
898 645 class PlotDOPData(PlotRTIData):
646 '''
647 Plot for DOPPLER Data
648 '''
649
899 650 CODE = 'dop'
900 651 colormap = 'jet'
901 652
902 653
903 class PlotPHASEData(PlotCOHData):
904 CODE = 'phase'
905 colormap = 'seismic'
906
907
908 654 class PlotSkyMapData(PlotData):
655 '''
656 Plot for meteors detection data
657 '''
909 658
910 659 CODE = 'met'
911 660
912 661 def setup(self):
913 662
914 663 self.ncols = 1
915 664 self.nrows = 1
916 665 self.width = 7.2
917 666 self.height = 7.2
918 667
919 668 self.xlabel = 'Zonal Zenith Angle (deg)'
920 669 self.ylabel = 'Meridional Zenith Angle (deg)'
921 670
922 671 if self.figure is None:
923 672 self.figure = plt.figure(figsize=(self.width, self.height),
924 673 edgecolor='k',
925 674 facecolor='w')
926 675 else:
927 676 self.figure.clf()
928 677
929 678 self.ax = plt.subplot2grid((self.nrows, self.ncols), (0, 0), 1, 1, polar=True)
930 679 self.ax.firsttime = True
931 680
932 681
933 682 def plot(self):
934 683
935 arrayParameters = np.concatenate([self.data['param'][t] for t in self.times])
684 arrayParameters = numpy.concatenate([self.data['param'][t] for t in self.data.times])
936 685 error = arrayParameters[:,-1]
937 686 indValid = numpy.where(error == 0)[0]
938 687 finalMeteor = arrayParameters[indValid,:]
939 688 finalAzimuth = finalMeteor[:,3]
940 689 finalZenith = finalMeteor[:,4]
941 690
942 691 x = finalAzimuth*numpy.pi/180
943 692 y = finalZenith
944 693
945 694 if self.ax.firsttime:
946 695 self.ax.plot = self.ax.plot(x, y, 'bo', markersize=5)[0]
947 696 self.ax.set_ylim(0,90)
948 697 self.ax.set_yticks(numpy.arange(0,90,20))
949 698 self.ax.set_xlabel(self.xlabel)
950 699 self.ax.set_ylabel(self.ylabel)
951 700 self.ax.yaxis.labelpad = 40
952 701 self.ax.firsttime = False
953 702 else:
954 703 self.ax.plot.set_data(x, y)
955 704
956 705
957 706 dt1 = datetime.datetime.fromtimestamp(self.min_time).strftime('%y/%m/%d %H:%M:%S')
958 707 dt2 = datetime.datetime.fromtimestamp(self.max_time).strftime('%y/%m/%d %H:%M:%S')
959 708 title = 'Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n' % (dt1,
960 709 dt2,
961 710 len(x))
962 711 self.ax.set_title(title, size=8)
963 712
964 713 self.saveTime = self.max_time
714
715 class PlotParamData(PlotRTIData):
716 '''
717 Plot for data_param object
718 '''
719
720 CODE = 'param'
721 colormap = 'seismic'
722
723 def setup(self):
724 self.xaxis = 'time'
725 self.ncols = 1
726 self.nrows = self.data.shape(self.CODE)[0]
727 self.nplots = self.nrows
728 if self.showSNR:
729 self.nrows += 1
730
731 self.ylabel = 'Height [Km]'
732 self.titles = self.data.parameters \
733 if self.data.parameters else ['Param {}'.format(x) for x in xrange(self.nrows)]
734 if self.showSNR:
735 self.titles.append('SNR')
736
737 def plot(self):
738 self.data.normalize_heights()
739 self.x = self.data.times
740 self.y = self.data.heights
741 if self.showSNR:
742 self.z = numpy.concatenate(
743 (self.data[self.CODE], self.data['snr'])
744 )
745 else:
746 self.z = self.data[self.CODE]
747
748 self.z = numpy.ma.masked_invalid(self.z)
749
750 for n, ax in enumerate(self.axes):
751
752 x, y, z = self.fill_gaps(*self.decimate())
753
754 if ax.firsttime:
755 if self.zlimits is not None:
756 self.zmin, self.zmax = self.zlimits[n]
757 self.zmax = self.zmax if self.zmax is not None else numpy.nanmax(abs(self.z[:-1, :]))
758 self.zmin = self.zmin if self.zmin is not None else -self.zmax
759 ax.plt = ax.pcolormesh(x, y, z[n, :, :].T*self.factors[n],
760 vmin=self.zmin,
761 vmax=self.zmax,
762 cmap=self.cmaps[n]
763 )
764 else:
765 if self.zlimits is not None:
766 self.zmin, self.zmax = self.zlimits[n]
767 ax.collections.remove(ax.collections[0])
768 ax.plt = ax.pcolormesh(x, y, z[n, :, :].T*self.factors[n],
769 vmin=self.zmin,
770 vmax=self.zmax,
771 cmap=self.cmaps[n]
772 )
773
774 self.saveTime = self.min_time
775
776 class PlotOuputData(PlotParamData):
777 '''
778 Plot data_output object
779 '''
780
781 CODE = 'output'
782 colormap = 'seismic' No newline at end of file
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
The requested commit or file is too big and content was truncated. Show full diff
Show file before | Show file after | Hide comments
@@ -1,904 +1,903
1 import itertools
2
1 3 import numpy
2 4
3 5 from jroproc_base import ProcessingUnit, Operation
4 6 from schainpy.model.data.jrodata import Spectra
5 7 from schainpy.model.data.jrodata import hildebrand_sekhon
6 8
7 9 class SpectraProc(ProcessingUnit):
8 10
9 11 def __init__(self, **kwargs):
10 12
11 13 ProcessingUnit.__init__(self, **kwargs)
12 14
13 15 self.buffer = None
14 16 self.firstdatatime = None
15 17 self.profIndex = 0
16 18 self.dataOut = Spectra()
17 19 self.id_min = None
18 20 self.id_max = None
19 21
20 22 def __updateSpecFromVoltage(self):
21 23
22 24 self.dataOut.timeZone = self.dataIn.timeZone
23 25 self.dataOut.dstFlag = self.dataIn.dstFlag
24 26 self.dataOut.errorCount = self.dataIn.errorCount
25 27 self.dataOut.useLocalTime = self.dataIn.useLocalTime
26 28
27 29 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
28 30 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
29 31 self.dataOut.channelList = self.dataIn.channelList
30 32 self.dataOut.heightList = self.dataIn.heightList
31 33 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
32 34
33 35 self.dataOut.nBaud = self.dataIn.nBaud
34 36 self.dataOut.nCode = self.dataIn.nCode
35 37 self.dataOut.code = self.dataIn.code
36 38 self.dataOut.nProfiles = self.dataOut.nFFTPoints
37 39
38 40 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
39 41 self.dataOut.utctime = self.firstdatatime
40 42 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
41 43 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
42 44 self.dataOut.flagShiftFFT = False
43 45
44 46 self.dataOut.nCohInt = self.dataIn.nCohInt
45 47 self.dataOut.nIncohInt = 1
46 48
47 49 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
48 50
49 51 self.dataOut.frequency = self.dataIn.frequency
50 52 self.dataOut.realtime = self.dataIn.realtime
51 53
52 54 self.dataOut.azimuth = self.dataIn.azimuth
53 55 self.dataOut.zenith = self.dataIn.zenith
54 56
55 57 self.dataOut.beam.codeList = self.dataIn.beam.codeList
56 58 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
57 59 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
58 60
59 61 def __getFft(self):
60 62 """
61 63 Convierte valores de Voltaje a Spectra
62 64
63 65 Affected:
64 66 self.dataOut.data_spc
65 67 self.dataOut.data_cspc
66 68 self.dataOut.data_dc
67 69 self.dataOut.heightList
68 70 self.profIndex
69 71 self.buffer
70 72 self.dataOut.flagNoData
71 73 """
72 74 fft_volt = numpy.fft.fft(self.buffer,n=self.dataOut.nFFTPoints,axis=1)
73 75 fft_volt = fft_volt.astype(numpy.dtype('complex'))
74 76 dc = fft_volt[:,0,:]
75 77
76 78 #calculo de self-spectra
77 79 fft_volt = numpy.fft.fftshift(fft_volt,axes=(1,))
78 80 spc = fft_volt * numpy.conjugate(fft_volt)
79 81 spc = spc.real
80 82
81 83 blocksize = 0
82 84 blocksize += dc.size
83 85 blocksize += spc.size
84 86
85 87 cspc = None
86 88 pairIndex = 0
87 89 if self.dataOut.pairsList != None:
88 90 #calculo de cross-spectra
89 91 cspc = numpy.zeros((self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
90 92 for pair in self.dataOut.pairsList:
91 93 if pair[0] not in self.dataOut.channelList:
92 94 raise ValueError, "Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" %(str(pair), str(self.dataOut.channelList))
93 95 if pair[1] not in self.dataOut.channelList:
94 96 raise ValueError, "Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" %(str(pair), str(self.dataOut.channelList))
95 97
96 98 cspc[pairIndex,:,:] = fft_volt[pair[0],:,:] * numpy.conjugate(fft_volt[pair[1],:,:])
97 99 pairIndex += 1
98 100 blocksize += cspc.size
99 101
100 102 self.dataOut.data_spc = spc
101 103 self.dataOut.data_cspc = cspc
102 104 self.dataOut.data_dc = dc
103 105 self.dataOut.blockSize = blocksize
104 106 self.dataOut.flagShiftFFT = True
105 107
106 108 def run(self, nProfiles=None, nFFTPoints=None, pairsList=[], ippFactor=None):
107 109
108 110 self.dataOut.flagNoData = True
109 111
110 112 if self.dataIn.type == "Spectra":
111 113 self.dataOut.copy(self.dataIn)
112 # self.__selectPairs(pairsList)
114 if not pairsList:
115 pairsList = itertools.combinations(self.dataOut.channelList, 2)
116 if self.dataOut.data_cspc is not None:
117 self.__selectPairs(pairsList)
113 118 return True
114 119
115 120 if self.dataIn.type == "Voltage":
116 121
117 122 if nFFTPoints == None:
118 123 raise ValueError, "This SpectraProc.run() need nFFTPoints input variable"
119 124
120 125 if nProfiles == None:
121 126 nProfiles = nFFTPoints
122 127
123 128 if ippFactor == None:
124 129 ippFactor = 1
125 130
126 131 self.dataOut.ippFactor = ippFactor
127 132
128 133 self.dataOut.nFFTPoints = nFFTPoints
129 134 self.dataOut.pairsList = pairsList
130 135
131 136 if self.buffer is None:
132 137 self.buffer = numpy.zeros( (self.dataIn.nChannels,
133 138 nProfiles,
134 139 self.dataIn.nHeights),
135 140 dtype='complex')
136 141
137 142 if self.dataIn.flagDataAsBlock:
138 143 #data dimension: [nChannels, nProfiles, nSamples]
139 144 nVoltProfiles = self.dataIn.data.shape[1]
140 145 # nVoltProfiles = self.dataIn.nProfiles
141 146
142 147 if nVoltProfiles == nProfiles:
143 148 self.buffer = self.dataIn.data.copy()
144 149 self.profIndex = nVoltProfiles
145 150
146 151 elif nVoltProfiles < nProfiles:
147 152
148 153 if self.profIndex == 0:
149 154 self.id_min = 0
150 155 self.id_max = nVoltProfiles
151 156
152 157 self.buffer[:,self.id_min:self.id_max,:] = self.dataIn.data
153 158 self.profIndex += nVoltProfiles
154 159 self.id_min += nVoltProfiles
155 160 self.id_max += nVoltProfiles
156 161 else:
157 162 raise ValueError, "The type object %s has %d profiles, it should just has %d profiles"%(self.dataIn.type,self.dataIn.data.shape[1],nProfiles)
158 163 self.dataOut.flagNoData = True
159 164 return 0
160 165 else:
161 166 self.buffer[:,self.profIndex,:] = self.dataIn.data.copy()
162 167 self.profIndex += 1
163 168
164 169 if self.firstdatatime == None:
165 170 self.firstdatatime = self.dataIn.utctime
166 171
167 172 if self.profIndex == nProfiles:
168 173 self.__updateSpecFromVoltage()
169 174 self.__getFft()
170 175
171 176 self.dataOut.flagNoData = False
172 177 self.firstdatatime = None
173 178 self.profIndex = 0
174 179
175 180 return True
176 181
177 182 raise ValueError, "The type of input object '%s' is not valid"%(self.dataIn.type)
178 183
179 184 def __selectPairs(self, pairsList):
180 185
181 if channelList == None:
186 if not pairsList:
182 187 return
183 188
184 pairsIndexListSelected = []
185
186 for thisPair in pairsList:
189 pairs = []
190 pairsIndex = []
187 191
188 if thisPair not in self.dataOut.pairsList:
192 for pair in pairsList:
193 if pair[0] not in self.dataOut.channelList or pair[1] not in self.dataOut.channelList:
189 194 continue
190
191 pairIndex = self.dataOut.pairsList.index(thisPair)
192
193 pairsIndexListSelected.append(pairIndex)
194
195 if not pairsIndexListSelected:
196 self.dataOut.data_cspc = None
197 self.dataOut.pairsList = []
198 return
199
200 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
201 self.dataOut.pairsList = [self.dataOut.pairsList[i] for i in pairsIndexListSelected]
195 pairs.append(pair)
196 pairsIndex.append(pairs.index(pair))
197
198 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndex]
199 self.dataOut.pairsList = pairs
200 self.dataOut.pairsIndexList = pairsIndex
202 201
203 202 return
204 203
205 204 def __selectPairsByChannel(self, channelList=None):
206 205
207 206 if channelList == None:
208 207 return
209 208
210 209 pairsIndexListSelected = []
211 210 for pairIndex in self.dataOut.pairsIndexList:
212 211 #First pair
213 212 if self.dataOut.pairsList[pairIndex][0] not in channelList:
214 213 continue
215 214 #Second pair
216 215 if self.dataOut.pairsList[pairIndex][1] not in channelList:
217 216 continue
218 217
219 218 pairsIndexListSelected.append(pairIndex)
220 219
221 220 if not pairsIndexListSelected:
222 221 self.dataOut.data_cspc = None
223 222 self.dataOut.pairsList = []
224 223 return
225 224
226 225 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
227 226 self.dataOut.pairsList = [self.dataOut.pairsList[i] for i in pairsIndexListSelected]
228 227
229 228 return
230 229
231 230 def selectChannels(self, channelList):
232 231
233 232 channelIndexList = []
234 233
235 234 for channel in channelList:
236 235 if channel not in self.dataOut.channelList:
237 236 raise ValueError, "Error selecting channels, Channel %d is not valid.\nAvailable channels = %s" %(channel, str(self.dataOut.channelList))
238 237
239 238 index = self.dataOut.channelList.index(channel)
240 239 channelIndexList.append(index)
241 240
242 241 self.selectChannelsByIndex(channelIndexList)
243 242
244 243 def selectChannelsByIndex(self, channelIndexList):
245 244 """
246 245 Selecciona un bloque de datos en base a canales segun el channelIndexList
247 246
248 247 Input:
249 248 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
250 249
251 250 Affected:
252 251 self.dataOut.data_spc
253 252 self.dataOut.channelIndexList
254 253 self.dataOut.nChannels
255 254
256 255 Return:
257 256 None
258 257 """
259 258
260 259 for channelIndex in channelIndexList:
261 260 if channelIndex not in self.dataOut.channelIndexList:
262 261 raise ValueError, "Error selecting channels: The value %d in channelIndexList is not valid.\nAvailable channel indexes = " %(channelIndex, self.dataOut.channelIndexList)
263 262
264 263 # nChannels = len(channelIndexList)
265 264
266 265 data_spc = self.dataOut.data_spc[channelIndexList,:]
267 266 data_dc = self.dataOut.data_dc[channelIndexList,:]
268 267
269 268 self.dataOut.data_spc = data_spc
270 269 self.dataOut.data_dc = data_dc
271 270
272 271 self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
273 272 # self.dataOut.nChannels = nChannels
274 273
275 274 self.__selectPairsByChannel(self.dataOut.channelList)
276 275
277 276 return 1
278 277
279 278 def selectHeights(self, minHei, maxHei):
280 279 """
281 280 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
282 281 minHei <= height <= maxHei
283 282
284 283 Input:
285 284 minHei : valor minimo de altura a considerar
286 285 maxHei : valor maximo de altura a considerar
287 286
288 287 Affected:
289 288 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
290 289
291 290 Return:
292 291 1 si el metodo se ejecuto con exito caso contrario devuelve 0
293 292 """
294 293
295 294 if (minHei > maxHei):
296 295 raise ValueError, "Error selecting heights: Height range (%d,%d) is not valid" % (minHei, maxHei)
297 296
298 297 if (minHei < self.dataOut.heightList[0]):
299 298 minHei = self.dataOut.heightList[0]
300 299
301 300 if (maxHei > self.dataOut.heightList[-1]):
302 301 maxHei = self.dataOut.heightList[-1]
303 302
304 303 minIndex = 0
305 304 maxIndex = 0
306 305 heights = self.dataOut.heightList
307 306
308 307 inda = numpy.where(heights >= minHei)
309 308 indb = numpy.where(heights <= maxHei)
310 309
311 310 try:
312 311 minIndex = inda[0][0]
313 312 except:
314 313 minIndex = 0
315 314
316 315 try:
317 316 maxIndex = indb[0][-1]
318 317 except:
319 318 maxIndex = len(heights)
320 319
321 320 self.selectHeightsByIndex(minIndex, maxIndex)
322 321
323 322 return 1
324 323
325 324 def getBeaconSignal(self, tauindex = 0, channelindex = 0, hei_ref=None):
326 325 newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
327 326
328 327 if hei_ref != None:
329 328 newheis = numpy.where(self.dataOut.heightList>hei_ref)
330 329
331 330 minIndex = min(newheis[0])
332 331 maxIndex = max(newheis[0])
333 332 data_spc = self.dataOut.data_spc[:,:,minIndex:maxIndex+1]
334 333 heightList = self.dataOut.heightList[minIndex:maxIndex+1]
335 334
336 335 # determina indices
337 336 nheis = int(self.dataOut.radarControllerHeaderObj.txB/(self.dataOut.heightList[1]-self.dataOut.heightList[0]))
338 337 avg_dB = 10*numpy.log10(numpy.sum(data_spc[channelindex,:,:],axis=0))
339 338 beacon_dB = numpy.sort(avg_dB)[-nheis:]
340 339 beacon_heiIndexList = []
341 340 for val in avg_dB.tolist():
342 341 if val >= beacon_dB[0]:
343 342 beacon_heiIndexList.append(avg_dB.tolist().index(val))
344 343
345 344 #data_spc = data_spc[:,:,beacon_heiIndexList]
346 345 data_cspc = None
347 346 if self.dataOut.data_cspc is not None:
348 347 data_cspc = self.dataOut.data_cspc[:,:,minIndex:maxIndex+1]
349 348 #data_cspc = data_cspc[:,:,beacon_heiIndexList]
350 349
351 350 data_dc = None
352 351 if self.dataOut.data_dc is not None:
353 352 data_dc = self.dataOut.data_dc[:,minIndex:maxIndex+1]
354 353 #data_dc = data_dc[:,beacon_heiIndexList]
355 354
356 355 self.dataOut.data_spc = data_spc
357 356 self.dataOut.data_cspc = data_cspc
358 357 self.dataOut.data_dc = data_dc
359 358 self.dataOut.heightList = heightList
360 359 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
361 360
362 361 return 1
363 362
364 363
365 364 def selectHeightsByIndex(self, minIndex, maxIndex):
366 365 """
367 366 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
368 367 minIndex <= index <= maxIndex
369 368
370 369 Input:
371 370 minIndex : valor de indice minimo de altura a considerar
372 371 maxIndex : valor de indice maximo de altura a considerar
373 372
374 373 Affected:
375 374 self.dataOut.data_spc
376 375 self.dataOut.data_cspc
377 376 self.dataOut.data_dc
378 377 self.dataOut.heightList
379 378
380 379 Return:
381 380 1 si el metodo se ejecuto con exito caso contrario devuelve 0
382 381 """
383 382
384 383 if (minIndex < 0) or (minIndex > maxIndex):
385 384 raise ValueError, "Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex)
386 385
387 386 if (maxIndex >= self.dataOut.nHeights):
388 387 maxIndex = self.dataOut.nHeights-1
389 388
390 389 #Spectra
391 390 data_spc = self.dataOut.data_spc[:,:,minIndex:maxIndex+1]
392 391
393 392 data_cspc = None
394 393 if self.dataOut.data_cspc is not None:
395 394 data_cspc = self.dataOut.data_cspc[:,:,minIndex:maxIndex+1]
396 395
397 396 data_dc = None
398 397 if self.dataOut.data_dc is not None:
399 398 data_dc = self.dataOut.data_dc[:,minIndex:maxIndex+1]
400 399
401 400 self.dataOut.data_spc = data_spc
402 401 self.dataOut.data_cspc = data_cspc
403 402 self.dataOut.data_dc = data_dc
404 403
405 404 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex+1]
406 405
407 406 return 1
408 407
409 408 def removeDC(self, mode = 2):
410 409 jspectra = self.dataOut.data_spc
411 410 jcspectra = self.dataOut.data_cspc
412 411
413 412
414 413 num_chan = jspectra.shape[0]
415 414 num_hei = jspectra.shape[2]
416 415
417 416 if jcspectra is not None:
418 417 jcspectraExist = True
419 418 num_pairs = jcspectra.shape[0]
420 419 else: jcspectraExist = False
421 420
422 421 freq_dc = jspectra.shape[1]/2
423 422 ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
424 423
425 424 if ind_vel[0]<0:
426 425 ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
427 426
428 427 if mode == 1:
429 428 jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
430 429
431 430 if jcspectraExist:
432 431 jcspectra[:,freq_dc,:] = (jcspectra[:,ind_vel[1],:] + jcspectra[:,ind_vel[2],:])/2
433 432
434 433 if mode == 2:
435 434
436 435 vel = numpy.array([-2,-1,1,2])
437 436 xx = numpy.zeros([4,4])
438 437
439 438 for fil in range(4):
440 439 xx[fil,:] = vel[fil]**numpy.asarray(range(4))
441 440
442 441 xx_inv = numpy.linalg.inv(xx)
443 442 xx_aux = xx_inv[0,:]
444 443
445 444 for ich in range(num_chan):
446 445 yy = jspectra[ich,ind_vel,:]
447 446 jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy)
448 447
449 448 junkid = jspectra[ich,freq_dc,:]<=0
450 449 cjunkid = sum(junkid)
451 450
452 451 if cjunkid.any():
453 452 jspectra[ich,freq_dc,junkid.nonzero()] = (jspectra[ich,ind_vel[1],junkid] + jspectra[ich,ind_vel[2],junkid])/2
454 453
455 454 if jcspectraExist:
456 455 for ip in range(num_pairs):
457 456 yy = jcspectra[ip,ind_vel,:]
458 457 jcspectra[ip,freq_dc,:] = numpy.dot(xx_aux,yy)
459 458
460 459
461 460 self.dataOut.data_spc = jspectra
462 461 self.dataOut.data_cspc = jcspectra
463 462
464 463 return 1
465 464
466 465 def removeInterference(self, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None):
467 466
468 467 jspectra = self.dataOut.data_spc
469 468 jcspectra = self.dataOut.data_cspc
470 469 jnoise = self.dataOut.getNoise()
471 470 num_incoh = self.dataOut.nIncohInt
472 471
473 472 num_channel = jspectra.shape[0]
474 473 num_prof = jspectra.shape[1]
475 474 num_hei = jspectra.shape[2]
476 475
477 476 #hei_interf
478 477 if hei_interf is None:
479 478 count_hei = num_hei/2 #Como es entero no importa
480 479 hei_interf = numpy.asmatrix(range(count_hei)) + num_hei - count_hei
481 480 hei_interf = numpy.asarray(hei_interf)[0]
482 481 #nhei_interf
483 482 if (nhei_interf == None):
484 483 nhei_interf = 5
485 484 if (nhei_interf < 1):
486 485 nhei_interf = 1
487 486 if (nhei_interf > count_hei):
488 487 nhei_interf = count_hei
489 488 if (offhei_interf == None):
490 489 offhei_interf = 0
491 490
492 491 ind_hei = range(num_hei)
493 492 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
494 493 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
495 494 mask_prof = numpy.asarray(range(num_prof))
496 495 num_mask_prof = mask_prof.size
497 496 comp_mask_prof = [0, num_prof/2]
498 497
499 498
500 499 #noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
501 500 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
502 501 jnoise = numpy.nan
503 502 noise_exist = jnoise[0] < numpy.Inf
504 503
505 504 #Subrutina de Remocion de la Interferencia
506 505 for ich in range(num_channel):
507 506 #Se ordena los espectros segun su potencia (menor a mayor)
508 507 power = jspectra[ich,mask_prof,:]
509 508 power = power[:,hei_interf]
510 509 power = power.sum(axis = 0)
511 510 psort = power.ravel().argsort()
512 511
513 512 #Se estima la interferencia promedio en los Espectros de Potencia empleando
514 513 junkspc_interf = jspectra[ich,:,hei_interf[psort[range(offhei_interf, nhei_interf + offhei_interf)]]]
515 514
516 515 if noise_exist:
517 516 # tmp_noise = jnoise[ich] / num_prof
518 517 tmp_noise = jnoise[ich]
519 518 junkspc_interf = junkspc_interf - tmp_noise
520 519 #junkspc_interf[:,comp_mask_prof] = 0
521 520
522 521 jspc_interf = junkspc_interf.sum(axis = 0) / nhei_interf
523 522 jspc_interf = jspc_interf.transpose()
524 523 #Calculando el espectro de interferencia promedio
525 524 noiseid = numpy.where(jspc_interf <= tmp_noise/ numpy.sqrt(num_incoh))
526 525 noiseid = noiseid[0]
527 526 cnoiseid = noiseid.size
528 527 interfid = numpy.where(jspc_interf > tmp_noise/ numpy.sqrt(num_incoh))
529 528 interfid = interfid[0]
530 529 cinterfid = interfid.size
531 530
532 531 if (cnoiseid > 0): jspc_interf[noiseid] = 0
533 532
534 533 #Expandiendo los perfiles a limpiar
535 534 if (cinterfid > 0):
536 535 new_interfid = (numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof)%num_prof
537 536 new_interfid = numpy.asarray(new_interfid)
538 537 new_interfid = {x for x in new_interfid}
539 538 new_interfid = numpy.array(list(new_interfid))
540 539 new_cinterfid = new_interfid.size
541 540 else: new_cinterfid = 0
542 541
543 542 for ip in range(new_cinterfid):
544 543 ind = junkspc_interf[:,new_interfid[ip]].ravel().argsort()
545 544 jspc_interf[new_interfid[ip]] = junkspc_interf[ind[nhei_interf/2],new_interfid[ip]]
546 545
547 546
548 547 jspectra[ich,:,ind_hei] = jspectra[ich,:,ind_hei] - jspc_interf #Corregir indices
549 548
550 549 #Removiendo la interferencia del punto de mayor interferencia
551 550 ListAux = jspc_interf[mask_prof].tolist()
552 551 maxid = ListAux.index(max(ListAux))
553 552
554 553
555 554 if cinterfid > 0:
556 555 for ip in range(cinterfid*(interf == 2) - 1):
557 556 ind = (jspectra[ich,interfid[ip],:] < tmp_noise*(1 + 1/numpy.sqrt(num_incoh))).nonzero()
558 557 cind = len(ind)
559 558
560 559 if (cind > 0):
561 560 jspectra[ich,interfid[ip],ind] = tmp_noise*(1 + (numpy.random.uniform(cind) - 0.5)/numpy.sqrt(num_incoh))
562 561
563 562 ind = numpy.array([-2,-1,1,2])
564 563 xx = numpy.zeros([4,4])
565 564
566 565 for id1 in range(4):
567 566 xx[:,id1] = ind[id1]**numpy.asarray(range(4))
568 567
569 568 xx_inv = numpy.linalg.inv(xx)
570 569 xx = xx_inv[:,0]
571 570 ind = (ind + maxid + num_mask_prof)%num_mask_prof
572 571 yy = jspectra[ich,mask_prof[ind],:]
573 572 jspectra[ich,mask_prof[maxid],:] = numpy.dot(yy.transpose(),xx)
574 573
575 574
576 575 indAux = (jspectra[ich,:,:] < tmp_noise*(1-1/numpy.sqrt(num_incoh))).nonzero()
577 576 jspectra[ich,indAux[0],indAux[1]] = tmp_noise * (1 - 1/numpy.sqrt(num_incoh))
578 577
579 578 #Remocion de Interferencia en el Cross Spectra
580 579 if jcspectra is None: return jspectra, jcspectra
581 580 num_pairs = jcspectra.size/(num_prof*num_hei)
582 581 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
583 582
584 583 for ip in range(num_pairs):
585 584
586 585 #-------------------------------------------
587 586
588 587 cspower = numpy.abs(jcspectra[ip,mask_prof,:])
589 588 cspower = cspower[:,hei_interf]
590 589 cspower = cspower.sum(axis = 0)
591 590
592 591 cspsort = cspower.ravel().argsort()
593 592 junkcspc_interf = jcspectra[ip,:,hei_interf[cspsort[range(offhei_interf, nhei_interf + offhei_interf)]]]
594 593 junkcspc_interf = junkcspc_interf.transpose()
595 594 jcspc_interf = junkcspc_interf.sum(axis = 1)/nhei_interf
596 595
597 596 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
598 597
599 598 median_real = numpy.median(numpy.real(junkcspc_interf[mask_prof[ind[range(3*num_prof/4)]],:]))
600 599 median_imag = numpy.median(numpy.imag(junkcspc_interf[mask_prof[ind[range(3*num_prof/4)]],:]))
601 600 junkcspc_interf[comp_mask_prof,:] = numpy.complex(median_real, median_imag)
602 601
603 602 for iprof in range(num_prof):
604 603 ind = numpy.abs(junkcspc_interf[iprof,:]).ravel().argsort()
605 604 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf/2]]
606 605
607 606 #Removiendo la Interferencia
608 607 jcspectra[ip,:,ind_hei] = jcspectra[ip,:,ind_hei] - jcspc_interf
609 608
610 609 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
611 610 maxid = ListAux.index(max(ListAux))
612 611
613 612 ind = numpy.array([-2,-1,1,2])
614 613 xx = numpy.zeros([4,4])
615 614
616 615 for id1 in range(4):
617 616 xx[:,id1] = ind[id1]**numpy.asarray(range(4))
618 617
619 618 xx_inv = numpy.linalg.inv(xx)
620 619 xx = xx_inv[:,0]
621 620
622 621 ind = (ind + maxid + num_mask_prof)%num_mask_prof
623 622 yy = jcspectra[ip,mask_prof[ind],:]
624 623 jcspectra[ip,mask_prof[maxid],:] = numpy.dot(yy.transpose(),xx)
625 624
626 625 #Guardar Resultados
627 626 self.dataOut.data_spc = jspectra
628 627 self.dataOut.data_cspc = jcspectra
629 628
630 629 return 1
631 630
632 631 def setRadarFrequency(self, frequency=None):
633 632
634 633 if frequency != None:
635 634 self.dataOut.frequency = frequency
636 635
637 636 return 1
638 637
639 638 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
640 639 #validacion de rango
641 640 if minHei == None:
642 641 minHei = self.dataOut.heightList[0]
643 642
644 643 if maxHei == None:
645 644 maxHei = self.dataOut.heightList[-1]
646 645
647 646 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
648 647 print 'minHei: %.2f is out of the heights range'%(minHei)
649 648 print 'minHei is setting to %.2f'%(self.dataOut.heightList[0])
650 649 minHei = self.dataOut.heightList[0]
651 650
652 651 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
653 652 print 'maxHei: %.2f is out of the heights range'%(maxHei)
654 653 print 'maxHei is setting to %.2f'%(self.dataOut.heightList[-1])
655 654 maxHei = self.dataOut.heightList[-1]
656 655
657 656 # validacion de velocidades
658 657 velrange = self.dataOut.getVelRange(1)
659 658
660 659 if minVel == None:
661 660 minVel = velrange[0]
662 661
663 662 if maxVel == None:
664 663 maxVel = velrange[-1]
665 664
666 665 if (minVel < velrange[0]) or (minVel > maxVel):
667 666 print 'minVel: %.2f is out of the velocity range'%(minVel)
668 667 print 'minVel is setting to %.2f'%(velrange[0])
669 668 minVel = velrange[0]
670 669
671 670 if (maxVel > velrange[-1]) or (maxVel < minVel):
672 671 print 'maxVel: %.2f is out of the velocity range'%(maxVel)
673 672 print 'maxVel is setting to %.2f'%(velrange[-1])
674 673 maxVel = velrange[-1]
675 674
676 675 # seleccion de indices para rango
677 676 minIndex = 0
678 677 maxIndex = 0
679 678 heights = self.dataOut.heightList
680 679
681 680 inda = numpy.where(heights >= minHei)
682 681 indb = numpy.where(heights <= maxHei)
683 682
684 683 try:
685 684 minIndex = inda[0][0]
686 685 except:
687 686 minIndex = 0
688 687
689 688 try:
690 689 maxIndex = indb[0][-1]
691 690 except:
692 691 maxIndex = len(heights)
693 692
694 693 if (minIndex < 0) or (minIndex > maxIndex):
695 694 raise ValueError, "some value in (%d,%d) is not valid" % (minIndex, maxIndex)
696 695
697 696 if (maxIndex >= self.dataOut.nHeights):
698 697 maxIndex = self.dataOut.nHeights-1
699 698
700 699 # seleccion de indices para velocidades
701 700 indminvel = numpy.where(velrange >= minVel)
702 701 indmaxvel = numpy.where(velrange <= maxVel)
703 702 try:
704 703 minIndexVel = indminvel[0][0]
705 704 except:
706 705 minIndexVel = 0
707 706
708 707 try:
709 708 maxIndexVel = indmaxvel[0][-1]
710 709 except:
711 710 maxIndexVel = len(velrange)
712 711
713 712 #seleccion del espectro
714 713 data_spc = self.dataOut.data_spc[:,minIndexVel:maxIndexVel+1,minIndex:maxIndex+1]
715 714 #estimacion de ruido
716 715 noise = numpy.zeros(self.dataOut.nChannels)
717 716
718 717 for channel in range(self.dataOut.nChannels):
719 718 daux = data_spc[channel,:,:]
720 719 noise[channel] = hildebrand_sekhon(daux, self.dataOut.nIncohInt)
721 720
722 721 self.dataOut.noise_estimation = noise.copy()
723 722
724 723 return 1
725 724
726 725 class IncohInt(Operation):
727 726
728 727
729 728 __profIndex = 0
730 729 __withOverapping = False
731 730
732 731 __byTime = False
733 732 __initime = None
734 733 __lastdatatime = None
735 734 __integrationtime = None
736 735
737 736 __buffer_spc = None
738 737 __buffer_cspc = None
739 738 __buffer_dc = None
740 739
741 740 __dataReady = False
742 741
743 742 __timeInterval = None
744 743
745 744 n = None
746 745
747 746
748 747
749 748 def __init__(self, **kwargs):
750 749
751 750 Operation.__init__(self, **kwargs)
752 751 # self.isConfig = False
753 752
754 753 def setup(self, n=None, timeInterval=None, overlapping=False):
755 754 """
756 755 Set the parameters of the integration class.
757 756
758 757 Inputs:
759 758
760 759 n : Number of coherent integrations
761 760 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
762 761 overlapping :
763 762
764 763 """
765 764
766 765 self.__initime = None
767 766 self.__lastdatatime = 0
768 767
769 768 self.__buffer_spc = 0
770 769 self.__buffer_cspc = 0
771 770 self.__buffer_dc = 0
772 771
773 772 self.__profIndex = 0
774 773 self.__dataReady = False
775 774 self.__byTime = False
776 775
777 776 if n is None and timeInterval is None:
778 777 raise ValueError, "n or timeInterval should be specified ..."
779 778
780 779 if n is not None:
781 780 self.n = int(n)
782 781 else:
783 782 self.__integrationtime = int(timeInterval) #if (type(timeInterval)!=integer) -> change this line
784 783 self.n = None
785 784 self.__byTime = True
786 785
787 786 def putData(self, data_spc, data_cspc, data_dc):
788 787
789 788 """
790 789 Add a profile to the __buffer_spc and increase in one the __profileIndex
791 790
792 791 """
793 792
794 793 self.__buffer_spc += data_spc
795 794
796 795 if data_cspc is None:
797 796 self.__buffer_cspc = None
798 797 else:
799 798 self.__buffer_cspc += data_cspc
800 799
801 800 if data_dc is None:
802 801 self.__buffer_dc = None
803 802 else:
804 803 self.__buffer_dc += data_dc
805 804
806 805 self.__profIndex += 1
807 806
808 807 return
809 808
810 809 def pushData(self):
811 810 """
812 811 Return the sum of the last profiles and the profiles used in the sum.
813 812
814 813 Affected:
815 814
816 815 self.__profileIndex
817 816
818 817 """
819 818
820 819 data_spc = self.__buffer_spc
821 820 data_cspc = self.__buffer_cspc
822 821 data_dc = self.__buffer_dc
823 822 n = self.__profIndex
824 823
825 824 self.__buffer_spc = 0
826 825 self.__buffer_cspc = 0
827 826 self.__buffer_dc = 0
828 827 self.__profIndex = 0
829 828
830 829 return data_spc, data_cspc, data_dc, n
831 830
832 831 def byProfiles(self, *args):
833 832
834 833 self.__dataReady = False
835 834 avgdata_spc = None
836 835 avgdata_cspc = None
837 836 avgdata_dc = None
838 837
839 838 self.putData(*args)
840 839
841 840 if self.__profIndex == self.n:
842 841
843 842 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
844 843 self.n = n
845 844 self.__dataReady = True
846 845
847 846 return avgdata_spc, avgdata_cspc, avgdata_dc
848 847
849 848 def byTime(self, datatime, *args):
850 849
851 850 self.__dataReady = False
852 851 avgdata_spc = None
853 852 avgdata_cspc = None
854 853 avgdata_dc = None
855 854
856 855 self.putData(*args)
857 856
858 857 if (datatime - self.__initime) >= self.__integrationtime:
859 858 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
860 859 self.n = n
861 860 self.__dataReady = True
862 861
863 862 return avgdata_spc, avgdata_cspc, avgdata_dc
864 863
865 864 def integrate(self, datatime, *args):
866 865
867 866 if self.__profIndex == 0:
868 867 self.__initime = datatime
869 868
870 869 if self.__byTime:
871 870 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(datatime, *args)
872 871 else:
873 872 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
874 873
875 874 if not self.__dataReady:
876 875 return None, None, None, None
877 876
878 877 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
879 878
880 879 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
881 880
882 881 if n==1:
883 882 return
884 883
885 884 dataOut.flagNoData = True
886 885
887 886 if not self.isConfig:
888 887 self.setup(n, timeInterval, overlapping)
889 888 self.isConfig = True
890 889
891 890 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
892 891 dataOut.data_spc,
893 892 dataOut.data_cspc,
894 893 dataOut.data_dc)
895 894
896 895 if self.__dataReady:
897 896
898 897 dataOut.data_spc = avgdata_spc
899 898 dataOut.data_cspc = avgdata_cspc
900 899 dataOut.data_dc = avgdata_dc
901 900
902 901 dataOut.nIncohInt *= self.n
903 902 dataOut.utctime = avgdatatime
904 903 dataOut.flagNoData = False
Show file before | Show file after | Hide comments
@@ -1,501 +1,604
1 1 '''
2 2 @author: Juan C. Espinoza
3 3 '''
4 4
5 5 import time
6 6 import json
7 7 import numpy
8 8 import paho.mqtt.client as mqtt
9 9 import zmq
10 10 import datetime
11 11 from zmq.utils.monitor import recv_monitor_message
12 12 from functools import wraps
13 13 from threading import Thread
14 14 from multiprocessing import Process
15 15
16 16 from schainpy.model.proc.jroproc_base import Operation, ProcessingUnit
17 17 from schainpy.model.data.jrodata import JROData
18 from schainpy.utils import log
18 19
19 20 MAXNUMX = 100
20 21 MAXNUMY = 100
21 22
22 23 class PrettyFloat(float):
23 24 def __repr__(self):
24 25 return '%.2f' % self
25 26
26 27 def roundFloats(obj):
27 28 if isinstance(obj, list):
28 29 return map(roundFloats, obj)
29 30 elif isinstance(obj, float):
30 31 return round(obj, 2)
31 32
32 33 def decimate(z, MAXNUMY):
33 # dx = int(len(self.x)/self.__MAXNUMX) + 1
34
35 34 dy = int(len(z[0])/MAXNUMY) + 1
36 35
37 36 return z[::, ::dy]
38 37
39 38 class throttle(object):
40 """Decorator that prevents a function from being called more than once every
39 '''
40 Decorator that prevents a function from being called more than once every
41 41 time period.
42 42 To create a function that cannot be called more than once a minute, but
43 43 will sleep until it can be called:
44 44 @throttle(minutes=1)
45 45 def foo():
46 46 pass
47 47
48 48 for i in range(10):
49 49 foo()
50 50 print "This function has run %s times." % i
51 """
51 '''
52 52
53 53 def __init__(self, seconds=0, minutes=0, hours=0):
54 54 self.throttle_period = datetime.timedelta(
55 55 seconds=seconds, minutes=minutes, hours=hours
56 56 )
57 57
58 58 self.time_of_last_call = datetime.datetime.min
59 59
60 60 def __call__(self, fn):
61 61 @wraps(fn)
62 62 def wrapper(*args, **kwargs):
63 63 now = datetime.datetime.now()
64 64 time_since_last_call = now - self.time_of_last_call
65 65 time_left = self.throttle_period - time_since_last_call
66 66
67 67 if time_left > datetime.timedelta(seconds=0):
68 68 return
69 69
70 70 self.time_of_last_call = datetime.datetime.now()
71 71 return fn(*args, **kwargs)
72 72
73 73 return wrapper
74 74
75 class Data(object):
76 '''
77 Object to hold data to be plotted
78 '''
79
80 def __init__(self, plottypes, throttle_value):
81 self.plottypes = plottypes
82 self.throttle = throttle_value
83 self.ended = False
84 self.__times = []
85
86 def __str__(self):
87 dum = ['{}{}'.format(key, self.shape(key)) for key in self.data]
88 return 'Data[{}][{}]'.format(';'.join(dum), len(self.__times))
89
90 def __len__(self):
91 return len(self.__times)
92
93 def __getitem__(self, key):
94 if key not in self.data:
95 raise KeyError(log.error('Missing key: {}'.format(key)))
96
97 if 'spc' in key:
98 ret = self.data[key]
99 else:
100 ret = numpy.array([self.data[key][x] for x in self.times])
101 if ret.ndim > 1:
102 ret = numpy.swapaxes(ret, 0, 1)
103 return ret
104
105 def setup(self):
106 '''
107 Configure object
108 '''
109
110 self.ended = False
111 self.data = {}
112 self.__times = []
113 self.__heights = []
114 self.__all_heights = set()
115 for plot in self.plottypes:
116 self.data[plot] = {}
117
118 def shape(self, key):
119 '''
120 Get the shape of the one-element data for the given key
121 '''
122
123 if len(self.data[key]):
124 if 'spc' in key:
125 return self.data[key].shape
126 return self.data[key][self.__times[0]].shape
127 return (0,)
128
129 def update(self, dataOut):
130 '''
131 Update data object with new dataOut
132 '''
133
134 tm = dataOut.utctime
135 if tm in self.__times:
136 return
137
138 self.parameters = getattr(dataOut, 'parameters', [])
139 self.pairs = dataOut.pairsList
140 self.channels = dataOut.channelList
141 self.xrange = (dataOut.getFreqRange(1)/1000. , dataOut.getAcfRange(1) , dataOut.getVelRange(1))
142 self.interval = dataOut.getTimeInterval()
143 self.__heights.append(dataOut.heightList)
144 self.__all_heights.update(dataOut.heightList)
145 self.__times.append(tm)
146
147 for plot in self.plottypes:
148 if plot == 'spc':
149 z = dataOut.data_spc/dataOut.normFactor
150 self.data[plot] = 10*numpy.log10(z)
151 if plot == 'cspc':
152 self.data[plot] = dataOut.data_cspc
153 if plot == 'noise':
154 self.data[plot][tm] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
155 if plot == 'rti':
156 self.data[plot][tm] = dataOut.getPower()
157 if plot == 'snr_db':
158 self.data['snr'][tm] = dataOut.data_SNR
159 if plot == 'snr':
160 self.data[plot][tm] = 10*numpy.log10(dataOut.data_SNR)
161 if plot == 'dop':
162 self.data[plot][tm] = 10*numpy.log10(dataOut.data_DOP)
163 if plot == 'mean':
164 self.data[plot][tm] = dataOut.data_MEAN
165 if plot == 'std':
166 self.data[plot][tm] = dataOut.data_STD
167 if plot == 'coh':
168 self.data[plot][tm] = dataOut.getCoherence()
169 if plot == 'phase':
170 self.data[plot][tm] = dataOut.getCoherence(phase=True)
171 if plot == 'output':
172 self.data[plot][tm] = dataOut.data_output
173 if plot == 'param':
174 self.data[plot][tm] = dataOut.data_param
175
176 def normalize_heights(self):
177 '''
178 Ensure same-dimension of the data for different heighList
179 '''
180
181 H = numpy.array(list(self.__all_heights))
182 H.sort()
183 for key in self.data:
184 shape = self.shape(key)[:-1] + H.shape
185 for tm, obj in self.data[key].items():
186 h = self.__heights[self.__times.index(tm)]
187 if H.size == h.size:
188 continue
189 index = numpy.where(numpy.in1d(H, h))[0]
190 dummy = numpy.zeros(shape) + numpy.nan
191 if len(shape) == 2:
192 dummy[:, index] = obj
193 else:
194 dummy[index] = obj
195 self.data[key][tm] = dummy
196
197 self.__heights = [H for tm in self.__times]
198
199 def jsonify(self, decimate=False):
200 '''
201 Convert data to json
202 '''
203
204 ret = {}
205 tm = self.times[-1]
206
207 for key, value in self.data:
208 if key in ('spc', 'cspc'):
209 ret[key] = roundFloats(self.data[key].to_list())
210 else:
211 ret[key] = roundFloats(self.data[key][tm].to_list())
212
213 ret['timestamp'] = tm
214 ret['interval'] = self.interval
215
216 @property
217 def times(self):
218 '''
219 Return the list of times of the current data
220 '''
221
222 ret = numpy.array(self.__times)
223 ret.sort()
224 return ret
225
226 @property
227 def heights(self):
228 '''
229 Return the list of heights of the current data
230 '''
231
232 return numpy.array(self.__heights[-1])
75 233
76 234 class PublishData(Operation):
77 """Clase publish."""
235 '''
236 Operation to send data over zmq.
237 '''
78 238
79 239 def __init__(self, **kwargs):
80 240 """Inicio."""
81 241 Operation.__init__(self, **kwargs)
82 242 self.isConfig = False
83 243 self.client = None
84 244 self.zeromq = None
85 245 self.mqtt = None
86 246
87 247 def on_disconnect(self, client, userdata, rc):
88 248 if rc != 0:
89 print("Unexpected disconnection.")
249 log.warning('Unexpected disconnection.')
90 250 self.connect()
91 251
92 252 def connect(self):
93 print 'trying to connect'
253 log.warning('trying to connect')
94 254 try:
95 255 self.client.connect(
96 256 host=self.host,
97 257 port=self.port,
98 258 keepalive=60*10,
99 259 bind_address='')
100 260 self.client.loop_start()
101 261 # self.client.publish(
102 262 # self.topic + 'SETUP',
103 263 # json.dumps(setup),
104 264 # retain=True
105 265 # )
106 266 except:
107 print "MQTT Conection error."
267 log.error('MQTT Conection error.')
108 268 self.client = False
109 269
110 270 def setup(self, port=1883, username=None, password=None, clientId="user", zeromq=1, verbose=True, **kwargs):
111 271 self.counter = 0
112 272 self.topic = kwargs.get('topic', 'schain')
113 273 self.delay = kwargs.get('delay', 0)
114 274 self.plottype = kwargs.get('plottype', 'spectra')
115 275 self.host = kwargs.get('host', "10.10.10.82")
116 276 self.port = kwargs.get('port', 3000)
117 277 self.clientId = clientId
118 278 self.cnt = 0
119 279 self.zeromq = zeromq
120 280 self.mqtt = kwargs.get('plottype', 0)
121 281 self.client = None
122 self.verbose = verbose
123 self.dataOut.firstdata = True
282 self.verbose = verbose
124 283 setup = []
125 284 if mqtt is 1:
126 285 self.client = mqtt.Client(
127 286 client_id=self.clientId + self.topic + 'SCHAIN',
128 287 clean_session=True)
129 288 self.client.on_disconnect = self.on_disconnect
130 289 self.connect()
131 290 for plot in self.plottype:
132 291 setup.append({
133 292 'plot': plot,
134 293 'topic': self.topic + plot,
135 294 'title': getattr(self, plot + '_' + 'title', False),
136 295 'xlabel': getattr(self, plot + '_' + 'xlabel', False),
137 296 'ylabel': getattr(self, plot + '_' + 'ylabel', False),
138 297 'xrange': getattr(self, plot + '_' + 'xrange', False),
139 298 'yrange': getattr(self, plot + '_' + 'yrange', False),
140 299 'zrange': getattr(self, plot + '_' + 'zrange', False),
141 300 })
142 301 if zeromq is 1:
143 302 context = zmq.Context()
144 303 self.zmq_socket = context.socket(zmq.PUSH)
145 304 server = kwargs.get('server', 'zmq.pipe')
146 305
147 306 if 'tcp://' in server:
148 307 address = server
149 308 else:
150 309 address = 'ipc:///tmp/%s' % server
151 310
152 311 self.zmq_socket.connect(address)
153 312 time.sleep(1)
154 313
155 314
156 315 def publish_data(self):
157 316 self.dataOut.finished = False
158 317 if self.mqtt is 1:
159 318 yData = self.dataOut.heightList[:2].tolist()
160 319 if self.plottype == 'spectra':
161 320 data = getattr(self.dataOut, 'data_spc')
162 321 z = data/self.dataOut.normFactor
163 322 zdB = 10*numpy.log10(z)
164 323 xlen, ylen = zdB[0].shape
165 324 dx = int(xlen/MAXNUMX) + 1
166 325 dy = int(ylen/MAXNUMY) + 1
167 326 Z = [0 for i in self.dataOut.channelList]
168 327 for i in self.dataOut.channelList:
169 328 Z[i] = zdB[i][::dx, ::dy].tolist()
170 329 payload = {
171 330 'timestamp': self.dataOut.utctime,
172 331 'data': roundFloats(Z),
173 332 'channels': ['Ch %s' % ch for ch in self.dataOut.channelList],
174 333 'interval': self.dataOut.getTimeInterval(),
175 334 'type': self.plottype,
176 335 'yData': yData
177 336 }
178 # print payload
179 337
180 338 elif self.plottype in ('rti', 'power'):
181 339 data = getattr(self.dataOut, 'data_spc')
182 340 z = data/self.dataOut.normFactor
183 341 avg = numpy.average(z, axis=1)
184 342 avgdB = 10*numpy.log10(avg)
185 343 xlen, ylen = z[0].shape
186 344 dy = numpy.floor(ylen/self.__MAXNUMY) + 1
187 345 AVG = [0 for i in self.dataOut.channelList]
188 346 for i in self.dataOut.channelList:
189 347 AVG[i] = avgdB[i][::dy].tolist()
190 348 payload = {
191 349 'timestamp': self.dataOut.utctime,
192 350 'data': roundFloats(AVG),
193 351 'channels': ['Ch %s' % ch for ch in self.dataOut.channelList],
194 352 'interval': self.dataOut.getTimeInterval(),
195 353 'type': self.plottype,
196 354 'yData': yData
197 355 }
198 356 elif self.plottype == 'noise':
199 357 noise = self.dataOut.getNoise()/self.dataOut.normFactor
200 358 noisedB = 10*numpy.log10(noise)
201 359 payload = {
202 360 'timestamp': self.dataOut.utctime,
203 361 'data': roundFloats(noisedB.reshape(-1, 1).tolist()),
204 362 'channels': ['Ch %s' % ch for ch in self.dataOut.channelList],
205 363 'interval': self.dataOut.getTimeInterval(),
206 364 'type': self.plottype,
207 365 'yData': yData
208 366 }
209 367 elif self.plottype == 'snr':
210 368 data = getattr(self.dataOut, 'data_SNR')
211 369 avgdB = 10*numpy.log10(data)
212 370
213 371 ylen = data[0].size
214 372 dy = numpy.floor(ylen/self.__MAXNUMY) + 1
215 373 AVG = [0 for i in self.dataOut.channelList]
216 374 for i in self.dataOut.channelList:
217 375 AVG[i] = avgdB[i][::dy].tolist()
218 376 payload = {
219 377 'timestamp': self.dataOut.utctime,
220 378 'data': roundFloats(AVG),
221 379 'channels': ['Ch %s' % ch for ch in self.dataOut.channelList],
222 380 'type': self.plottype,
223 381 'yData': yData
224 382 }
225 383 else:
226 384 print "Tipo de grafico invalido"
227 385 payload = {
228 386 'data': 'None',
229 387 'timestamp': 'None',
230 388 'type': None
231 389 }
232 # print 'Publishing data to {}'.format(self.host)
390
233 391 self.client.publish(self.topic + self.plottype, json.dumps(payload), qos=0)
234 392
235 393 if self.zeromq is 1:
236 394 if self.verbose:
237 print '[Sending] {} - {}'.format(self.dataOut.type, self.dataOut.datatime)
395 log.log(
396 '{} - {}'.format(self.dataOut.type, self.dataOut.datatime),
397 'Sending'
398 )
238 399 self.zmq_socket.send_pyobj(self.dataOut)
239 self.dataOut.firstdata = False
240
241 400
242 401 def run(self, dataOut, **kwargs):
243 402 self.dataOut = dataOut
244 403 if not self.isConfig:
245 404 self.setup(**kwargs)
246 405 self.isConfig = True
247 406
248 407 self.publish_data()
249 408 time.sleep(self.delay)
250 409
251 410 def close(self):
252 411 if self.zeromq is 1:
253 412 self.dataOut.finished = True
254 413 self.zmq_socket.send_pyobj(self.dataOut)
414 time.sleep(0.1)
255 415 self.zmq_socket.close()
256 416 if self.client:
257 417 self.client.loop_stop()
258 418 self.client.disconnect()
259 419
260 420
261 421 class ReceiverData(ProcessingUnit):
262 422
263 423 def __init__(self, **kwargs):
264 424
265 425 ProcessingUnit.__init__(self, **kwargs)
266 426
267 427 self.isConfig = False
268 428 server = kwargs.get('server', 'zmq.pipe')
269 429 if 'tcp://' in server:
270 430 address = server
271 431 else:
272 432 address = 'ipc:///tmp/%s' % server
273 433
274 434 self.address = address
275 435 self.dataOut = JROData()
276 436
277 437 def setup(self):
278 438
279 439 self.context = zmq.Context()
280 440 self.receiver = self.context.socket(zmq.PULL)
281 441 self.receiver.bind(self.address)
282 442 time.sleep(0.5)
283 print '[Starting] ReceiverData from {}'.format(self.address)
443 log.success('ReceiverData from {}'.format(self.address))
284 444
285 445
286 446 def run(self):
287 447
288 448 if not self.isConfig:
289 449 self.setup()
290 450 self.isConfig = True
291 451
292 452 self.dataOut = self.receiver.recv_pyobj()
293 print '[Receiving] {} - {}'.format(self.dataOut.type,
294 self.dataOut.datatime.ctime())
453 log.log('{} - {}'.format(self.dataOut.type,
454 self.dataOut.datatime.ctime(),),
455 'Receiving')
295 456
296 457
297 458 class PlotterReceiver(ProcessingUnit, Process):
298 459
299 460 throttle_value = 5
300 461
301 462 def __init__(self, **kwargs):
302 463
303 464 ProcessingUnit.__init__(self, **kwargs)
304 465 Process.__init__(self)
305 466 self.mp = False
306 467 self.isConfig = False
307 468 self.isWebConfig = False
308 self.plottypes = []
309 469 self.connections = 0
310 470 server = kwargs.get('server', 'zmq.pipe')
311 471 plot_server = kwargs.get('plot_server', 'zmq.web')
312 472 if 'tcp://' in server:
313 473 address = server
314 474 else:
315 475 address = 'ipc:///tmp/%s' % server
316 476
317 477 if 'tcp://' in plot_server:
318 478 plot_address = plot_server
319 479 else:
320 480 plot_address = 'ipc:///tmp/%s' % plot_server
321 481
322 482 self.address = address
323 483 self.plot_address = plot_address
324 484 self.plottypes = [s.strip() for s in kwargs.get('plottypes', 'rti').split(',')]
325 485 self.realtime = kwargs.get('realtime', False)
326 486 self.throttle_value = kwargs.get('throttle', 5)
327 487 self.sendData = self.initThrottle(self.throttle_value)
488 self.dates = []
328 489 self.setup()
329 490
330 491 def setup(self):
331 492
332 self.data = {}
333 self.data['times'] = []
334 for plottype in self.plottypes:
335 self.data[plottype] = {}
336 self.data['noise'] = {}
337 self.data['throttle'] = self.throttle_value
338 self.data['ENDED'] = False
339 self.isConfig = True
340 self.data_web = {}
493 self.data = Data(self.plottypes, self.throttle_value)
494 self.isConfig = True
341 495
342 496 def event_monitor(self, monitor):
343 497
344 498 events = {}
345 499
346 500 for name in dir(zmq):
347 501 if name.startswith('EVENT_'):
348 502 value = getattr(zmq, name)
349 503 events[value] = name
350 504
351 505 while monitor.poll():
352 506 evt = recv_monitor_message(monitor)
353 507 if evt['event'] == 32:
354 508 self.connections += 1
355 509 if evt['event'] == 512:
356 510 pass
357 if self.connections == 0 and self.started is True:
358 self.ended = True
359 511
360 512 evt.update({'description': events[evt['event']]})
361 513
362 514 if evt['event'] == zmq.EVENT_MONITOR_STOPPED:
363 515 break
364 516 monitor.close()
365 print("event monitor thread done!")
517 print('event monitor thread done!')
366 518
367 519 def initThrottle(self, throttle_value):
368 520
369 521 @throttle(seconds=throttle_value)
370 522 def sendDataThrottled(fn_sender, data):
371 523 fn_sender(data)
372 524
373 525 return sendDataThrottled
374 526
375
376 527 def send(self, data):
377 # print '[sending] data=%s size=%s' % (data.keys(), len(data['times']))
528 log.success('Sending {}'.format(data), self.name)
378 529 self.sender.send_pyobj(data)
379 530
380
381 def update(self):
382 t = self.dataOut.utctime
383
384 if t in self.data['times']:
385 return
386
387 self.data['times'].append(t)
388 self.data['dataOut'] = self.dataOut
389
390 for plottype in self.plottypes:
391 if plottype == 'spc':
392 z = self.dataOut.data_spc/self.dataOut.normFactor
393 self.data[plottype] = 10*numpy.log10(z)
394 self.data['noise'][t] = 10*numpy.log10(self.dataOut.getNoise()/self.dataOut.normFactor)
395 if plottype == 'cspc':
396 jcoherence = self.dataOut.data_cspc/numpy.sqrt(self.dataOut.data_spc*self.dataOut.data_spc)
397 self.data['cspc_coh'] = numpy.abs(jcoherence)
398 self.data['cspc_phase'] = numpy.arctan2(jcoherence.imag, jcoherence.real)*180/numpy.pi
399 if plottype == 'rti':
400 self.data[plottype][t] = self.dataOut.getPower()
401 if plottype == 'snr':
402 self.data[plottype][t] = 10*numpy.log10(self.dataOut.data_SNR)
403 if plottype == 'dop':
404 self.data[plottype][t] = 10*numpy.log10(self.dataOut.data_DOP)
405 if plottype == 'mean':
406 self.data[plottype][t] = self.dataOut.data_MEAN
407 if plottype == 'std':
408 self.data[plottype][t] = self.dataOut.data_STD
409 if plottype == 'coh':
410 self.data[plottype][t] = self.dataOut.getCoherence()
411 if plottype == 'phase':
412 self.data[plottype][t] = self.dataOut.getCoherence(phase=True)
413 if plottype == 'output':
414 self.data[plottype][t] = self.dataOut.data_output
415 if plottype == 'param':
416 self.data[plottype][t] = self.dataOut.data_param
417 if self.realtime:
418 self.data_web['timestamp'] = t
419 if plottype == 'spc':
420 self.data_web[plottype] = roundFloats(decimate(self.data[plottype]).tolist())
421 elif plottype == 'cspc':
422 self.data_web['cspc_coh'] = roundFloats(decimate(self.data['cspc_coh']).tolist())
423 self.data_web['cspc_phase'] = roundFloats(decimate(self.data['cspc_phase']).tolist())
424 elif plottype == 'noise':
425 self.data_web['noise'] = roundFloats(self.data['noise'][t].tolist())
426 else:
427 self.data_web[plottype] = roundFloats(decimate(self.data[plottype][t]).tolist())
428 self.data_web['interval'] = self.dataOut.getTimeInterval()
429 self.data_web['type'] = plottype
430
431 531 def run(self):
432 532
433 print '[Starting] {} from {}'.format(self.name, self.address)
533 log.success(
534 'Starting from {}'.format(self.address),
535 self.name
536 )
434 537
435 538 self.context = zmq.Context()
436 539 self.receiver = self.context.socket(zmq.PULL)
437 540 self.receiver.bind(self.address)
438 541 monitor = self.receiver.get_monitor_socket()
439 542 self.sender = self.context.socket(zmq.PUB)
440 543 if self.realtime:
441 544 self.sender_web = self.context.socket(zmq.PUB)
442 545 self.sender_web.connect(self.plot_address)
443 546 time.sleep(1)
444 547
445 548 if 'server' in self.kwargs:
446 549 self.sender.bind("ipc:///tmp/{}.plots".format(self.kwargs['server']))
447 550 else:
448 551 self.sender.bind("ipc:///tmp/zmq.plots")
449 552
450 time.sleep(3)
553 time.sleep(2)
451 554
452 555 t = Thread(target=self.event_monitor, args=(monitor,))
453 556 t.start()
454 557
455 558 while True:
456 self.dataOut = self.receiver.recv_pyobj()
457 # print '[Receiving] {} - {}'.format(self.dataOut.type,
458 # self.dataOut.datatime.ctime())
459
460 self.update()
559 dataOut = self.receiver.recv_pyobj()
560 dt = datetime.datetime.fromtimestamp(dataOut.utctime).date()
561 sended = False
562 if dt not in self.dates:
563 if self.data:
564 self.data.ended = True
565 self.send(self.data)
566 sended = True
567 self.data.setup()
568 self.dates.append(dt)
461 569
462 if self.dataOut.firstdata is True:
463 self.data['STARTED'] = True
570 self.data.update(dataOut)
464 571
465 if self.dataOut.finished is True:
466 self.send(self.data)
572 if dataOut.finished is True:
467 573 self.connections -= 1
468 if self.connections == 0 and self.started:
469 self.ended = True
470 self.data['ENDED'] = True
574 if self.connections == 0 and dt in self.dates:
575 self.data.ended = True
471 576 self.send(self.data)
472 self.setup()
473 self.started = False
577 self.data.setup()
474 578 else:
475 579 if self.realtime:
476 580 self.send(self.data)
477 self.sender_web.send_string(json.dumps(self.data_web))
581 # self.sender_web.send_string(self.data.jsonify())
478 582 else:
479 self.sendData(self.send, self.data)
480 self.started = True
583 if not sended:
584 self.sendData(self.send, self.data)
481 585
482 self.data['STARTED'] = False
483 586 return
484 587
485 588 def sendToWeb(self):
486 589
487 590 if not self.isWebConfig:
488 591 context = zmq.Context()
489 592 sender_web_config = context.socket(zmq.PUB)
490 593 if 'tcp://' in self.plot_address:
491 594 dum, address, port = self.plot_address.split(':')
492 595 conf_address = '{}:{}:{}'.format(dum, address, int(port)+1)
493 596 else:
494 597 conf_address = self.plot_address + '.config'
495 598 sender_web_config.bind(conf_address)
496 599 time.sleep(1)
497 600 for kwargs in self.operationKwargs.values():
498 601 if 'plot' in kwargs:
499 print '[Sending] Config data to web for {}'.format(kwargs['code'].upper())
602 log.success('[Sending] Config data to web for {}'.format(kwargs['code'].upper()))
500 603 sender_web_config.send_string(json.dumps(kwargs))
501 self.isWebConfig = True
604 self.isWebConfig = True No newline at end of file
General Comments 0
You need to be logged in to leave comments. Login now