##// END OF EJS Templates
schain
RSS
Add SkyMapPlotData, operation can access parent kwargs, fix server plot for multiple ReceiverData
Juan C. Espinoza -
r937:6cdcc42f2694
parent child
Show More
Show file before | Show file after | Hide comments
@@ -1,1225 +1,1228
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 __init__(self):
118 118
119 119 raise NotImplementedError
120 120
121 121 def copy(self, inputObj=None):
122 122
123 123 if inputObj == None:
124 124 return copy.deepcopy(self)
125 125
126 126 for key in inputObj.__dict__.keys():
127 127
128 128 attribute = inputObj.__dict__[key]
129 129
130 130 #If this attribute is a tuple or list
131 131 if type(inputObj.__dict__[key]) in (tuple, list):
132 132 self.__dict__[key] = attribute[:]
133 133 continue
134 134
135 135 #If this attribute is another object or instance
136 136 if hasattr(attribute, '__dict__'):
137 137 self.__dict__[key] = attribute.copy()
138 138 continue
139 139
140 140 self.__dict__[key] = inputObj.__dict__[key]
141 141
142 142 def deepcopy(self):
143 143
144 144 return copy.deepcopy(self)
145 145
146 146 def isEmpty(self):
147 147
148 148 return self.flagNoData
149 149
150 150 class JROData(GenericData):
151 151
152 152 # m_BasicHeader = BasicHeader()
153 153 # m_ProcessingHeader = ProcessingHeader()
154 154
155 155 systemHeaderObj = SystemHeader()
156 156
157 157 radarControllerHeaderObj = RadarControllerHeader()
158 158
159 159 # data = None
160 160
161 161 type = None
162 162
163 163 datatype = None #dtype but in string
164 164
165 165 # dtype = None
166 166
167 167 # nChannels = None
168 168
169 169 # nHeights = None
170 170
171 171 nProfiles = None
172 172
173 173 heightList = None
174 174
175 175 channelList = None
176 176
177 177 flagDiscontinuousBlock = False
178 178
179 179 useLocalTime = False
180 180
181 181 utctime = None
182 182
183 183 timeZone = None
184 184
185 185 dstFlag = None
186 186
187 187 errorCount = None
188 188
189 189 blocksize = None
190 190
191 191 # nCode = None
192 192 #
193 193 # nBaud = None
194 194 #
195 195 # code = None
196 196
197 197 flagDecodeData = False #asumo q la data no esta decodificada
198 198
199 199 flagDeflipData = False #asumo q la data no esta sin flip
200 200
201 201 flagShiftFFT = False
202 202
203 203 # ippSeconds = None
204 204
205 205 # timeInterval = None
206 206
207 207 nCohInt = None
208 208
209 209 # noise = None
210 210
211 211 windowOfFilter = 1
212 212
213 213 #Speed of ligth
214 214 C = 3e8
215 215
216 216 frequency = 49.92e6
217 217
218 218 realtime = False
219 219
220 220 beacon_heiIndexList = None
221 221
222 222 last_block = None
223 223
224 224 blocknow = None
225 225
226 226 azimuth = None
227 227
228 228 zenith = None
229 229
230 230 beam = Beam()
231 231
232 232 profileIndex = None
233 233
234 234 def __init__(self):
235 235
236 236 raise NotImplementedError
237 237
238 238 def getNoise(self):
239 239
240 240 raise NotImplementedError
241 241
242 242 def getNChannels(self):
243 243
244 244 return len(self.channelList)
245 245
246 246 def getChannelIndexList(self):
247 247
248 248 return range(self.nChannels)
249 249
250 250 def getNHeights(self):
251 251
252 252 return len(self.heightList)
253 253
254 254 def getHeiRange(self, extrapoints=0):
255 255
256 256 heis = self.heightList
257 257 # deltah = self.heightList[1] - self.heightList[0]
258 258 #
259 259 # heis.append(self.heightList[-1])
260 260
261 261 return heis
262 262
263 263 def getDeltaH(self):
264 264
265 265 delta = self.heightList[1] - self.heightList[0]
266 266
267 267 return delta
268 268
269 269 def getltctime(self):
270 270
271 271 if self.useLocalTime:
272 272 return self.utctime - self.timeZone*60
273 273
274 274 return self.utctime
275 275
276 276 def getDatatime(self):
277 277
278 278 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
279 279 return datatimeValue
280 280
281 281 def getTimeRange(self):
282 282
283 283 datatime = []
284 284
285 285 datatime.append(self.ltctime)
286 286 datatime.append(self.ltctime + self.timeInterval+1)
287 287
288 288 datatime = numpy.array(datatime)
289 289
290 290 return datatime
291 291
292 292 def getFmaxTimeResponse(self):
293 293
294 294 period = (10**-6)*self.getDeltaH()/(0.15)
295 295
296 296 PRF = 1./(period * self.nCohInt)
297 297
298 298 fmax = PRF
299 299
300 300 return fmax
301 301
302 302 def getFmax(self):
303 303
304 304 PRF = 1./(self.ippSeconds * self.nCohInt)
305 305
306 306 fmax = PRF
307 307
308 308 return fmax
309 309
310 310 def getVmax(self):
311 311
312 312 _lambda = self.C/self.frequency
313 313
314 314 vmax = self.getFmax() * _lambda/2
315 315
316 316 return vmax
317 317
318 318 def get_ippSeconds(self):
319 319 '''
320 320 '''
321 321 return self.radarControllerHeaderObj.ippSeconds
322 322
323 323 def set_ippSeconds(self, ippSeconds):
324 324 '''
325 325 '''
326 326
327 327 self.radarControllerHeaderObj.ippSeconds = ippSeconds
328 328
329 329 return
330 330
331 331 def get_dtype(self):
332 332 '''
333 333 '''
334 334 return getNumpyDtype(self.datatype)
335 335
336 336 def set_dtype(self, numpyDtype):
337 337 '''
338 338 '''
339 339
340 340 self.datatype = getDataTypeCode(numpyDtype)
341 341
342 342 def get_code(self):
343 343 '''
344 344 '''
345 345 return self.radarControllerHeaderObj.code
346 346
347 347 def set_code(self, code):
348 348 '''
349 349 '''
350 350 self.radarControllerHeaderObj.code = code
351 351
352 352 return
353 353
354 354 def get_ncode(self):
355 355 '''
356 356 '''
357 357 return self.radarControllerHeaderObj.nCode
358 358
359 359 def set_ncode(self, nCode):
360 360 '''
361 361 '''
362 362 self.radarControllerHeaderObj.nCode = nCode
363 363
364 364 return
365 365
366 366 def get_nbaud(self):
367 367 '''
368 368 '''
369 369 return self.radarControllerHeaderObj.nBaud
370 370
371 371 def set_nbaud(self, nBaud):
372 372 '''
373 373 '''
374 374 self.radarControllerHeaderObj.nBaud = nBaud
375 375
376 376 return
377 377
378 378 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
379 379 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
380 380 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
381 381 #noise = property(getNoise, "I'm the 'nHeights' property.")
382 382 datatime = property(getDatatime, "I'm the 'datatime' property")
383 383 ltctime = property(getltctime, "I'm the 'ltctime' property")
384 384 ippSeconds = property(get_ippSeconds, set_ippSeconds)
385 385 dtype = property(get_dtype, set_dtype)
386 386 # timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
387 387 code = property(get_code, set_code)
388 388 nCode = property(get_ncode, set_ncode)
389 389 nBaud = property(get_nbaud, set_nbaud)
390 390
391 391 class Voltage(JROData):
392 392
393 393 #data es un numpy array de 2 dmensiones (canales, alturas)
394 394 data = None
395 395
396 396 def __init__(self):
397 397 '''
398 398 Constructor
399 399 '''
400 400
401 401 self.useLocalTime = True
402 402
403 403 self.radarControllerHeaderObj = RadarControllerHeader()
404 404
405 405 self.systemHeaderObj = SystemHeader()
406 406
407 407 self.type = "Voltage"
408 408
409 409 self.data = None
410 410
411 411 # self.dtype = None
412 412
413 413 # self.nChannels = 0
414 414
415 415 # self.nHeights = 0
416 416
417 417 self.nProfiles = None
418 418
419 419 self.heightList = None
420 420
421 421 self.channelList = None
422 422
423 423 # self.channelIndexList = None
424 424
425 425 self.flagNoData = True
426 426
427 427 self.flagDiscontinuousBlock = False
428 428
429 429 self.utctime = None
430 430
431 431 self.timeZone = None
432 432
433 433 self.dstFlag = None
434 434
435 435 self.errorCount = None
436 436
437 437 self.nCohInt = None
438 438
439 439 self.blocksize = None
440 440
441 441 self.flagDecodeData = False #asumo q la data no esta decodificada
442 442
443 443 self.flagDeflipData = False #asumo q la data no esta sin flip
444 444
445 445 self.flagShiftFFT = False
446 446
447 447 self.flagDataAsBlock = False #Asumo que la data es leida perfil a perfil
448 448
449 449 self.profileIndex = 0
450 450
451 451 def getNoisebyHildebrand(self, channel = None):
452 452 """
453 453 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
454 454
455 455 Return:
456 456 noiselevel
457 457 """
458 458
459 459 if channel != None:
460 460 data = self.data[channel]
461 461 nChannels = 1
462 462 else:
463 463 data = self.data
464 464 nChannels = self.nChannels
465 465
466 466 noise = numpy.zeros(nChannels)
467 467 power = data * numpy.conjugate(data)
468 468
469 469 for thisChannel in range(nChannels):
470 470 if nChannels == 1:
471 471 daux = power[:].real
472 472 else:
473 473 daux = power[thisChannel,:].real
474 474 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
475 475
476 476 return noise
477 477
478 478 def getNoise(self, type = 1, channel = None):
479 479
480 480 if type == 1:
481 481 noise = self.getNoisebyHildebrand(channel)
482 482
483 483 return noise
484 484
485 485 def getPower(self, channel = None):
486 486
487 487 if channel != None:
488 488 data = self.data[channel]
489 489 else:
490 490 data = self.data
491 491
492 492 power = data * numpy.conjugate(data)
493 493 powerdB = 10*numpy.log10(power.real)
494 494 powerdB = numpy.squeeze(powerdB)
495 495
496 496 return powerdB
497 497
498 498 def getTimeInterval(self):
499 499
500 500 timeInterval = self.ippSeconds * self.nCohInt
501 501
502 502 return timeInterval
503 503
504 504 noise = property(getNoise, "I'm the 'nHeights' property.")
505 505 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
506 506
507 507 class Spectra(JROData):
508 508
509 509 #data spc es un numpy array de 2 dmensiones (canales, perfiles, alturas)
510 510 data_spc = None
511 511
512 512 #data cspc es un numpy array de 2 dmensiones (canales, pares, alturas)
513 513 data_cspc = None
514 514
515 515 #data dc es un numpy array de 2 dmensiones (canales, alturas)
516 516 data_dc = None
517 517
518 518 #data power
519 519 data_pwr = None
520 520
521 521 nFFTPoints = None
522 522
523 523 # nPairs = None
524 524
525 525 pairsList = None
526 526
527 527 nIncohInt = None
528 528
529 529 wavelength = None #Necesario para cacular el rango de velocidad desde la frecuencia
530 530
531 531 nCohInt = None #se requiere para determinar el valor de timeInterval
532 532
533 533 ippFactor = None
534 534
535 535 profileIndex = 0
536 536
537 537 plotting = "spectra"
538 538
539 539 def __init__(self):
540 540 '''
541 541 Constructor
542 542 '''
543 543
544 544 self.useLocalTime = True
545 545
546 546 self.radarControllerHeaderObj = RadarControllerHeader()
547 547
548 548 self.systemHeaderObj = SystemHeader()
549 549
550 550 self.type = "Spectra"
551 551
552 552 # self.data = None
553 553
554 554 # self.dtype = None
555 555
556 556 # self.nChannels = 0
557 557
558 558 # self.nHeights = 0
559 559
560 560 self.nProfiles = None
561 561
562 562 self.heightList = None
563 563
564 564 self.channelList = None
565 565
566 566 # self.channelIndexList = None
567 567
568 568 self.pairsList = None
569 569
570 570 self.flagNoData = True
571 571
572 572 self.flagDiscontinuousBlock = False
573 573
574 574 self.utctime = None
575 575
576 576 self.nCohInt = None
577 577
578 578 self.nIncohInt = None
579 579
580 580 self.blocksize = None
581 581
582 582 self.nFFTPoints = None
583 583
584 584 self.wavelength = None
585 585
586 586 self.flagDecodeData = False #asumo q la data no esta decodificada
587 587
588 588 self.flagDeflipData = False #asumo q la data no esta sin flip
589 589
590 590 self.flagShiftFFT = False
591 591
592 592 self.ippFactor = 1
593 593
594 594 #self.noise = None
595 595
596 596 self.beacon_heiIndexList = []
597 597
598 598 self.noise_estimation = None
599 599
600 600
601 601 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
602 602 """
603 603 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
604 604
605 605 Return:
606 606 noiselevel
607 607 """
608 608
609 609 noise = numpy.zeros(self.nChannels)
610 610
611 611 for channel in range(self.nChannels):
612 612 daux = self.data_spc[channel,xmin_index:xmax_index,ymin_index:ymax_index]
613 613 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
614 614
615 615 return noise
616 616
617 617 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
618 618
619 619 if self.noise_estimation is not None:
620 620 return self.noise_estimation #this was estimated by getNoise Operation defined in jroproc_spectra.py
621 621 else:
622 622 noise = self.getNoisebyHildebrand(xmin_index, xmax_index, ymin_index, ymax_index)
623 623 return noise
624 624
625 625 def getFreqRangeTimeResponse(self, extrapoints=0):
626 626
627 627 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints*self.ippFactor)
628 628 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
629 629
630 630 return freqrange
631 631
632 632 def getAcfRange(self, extrapoints=0):
633 633
634 634 deltafreq = 10./(self.getFmax() / (self.nFFTPoints*self.ippFactor))
635 635 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
636 636
637 637 return freqrange
638 638
639 639 def getFreqRange(self, extrapoints=0):
640 640
641 641 deltafreq = self.getFmax() / (self.nFFTPoints*self.ippFactor)
642 642 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
643 643
644 644 return freqrange
645 645
646 646 def getVelRange(self, extrapoints=0):
647 647
648 648 deltav = self.getVmax() / (self.nFFTPoints*self.ippFactor)
649 649 velrange = deltav*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) #- deltav/2
650 650
651 651 return velrange
652 652
653 653 def getNPairs(self):
654 654
655 655 return len(self.pairsList)
656 656
657 657 def getPairsIndexList(self):
658 658
659 659 return range(self.nPairs)
660 660
661 661 def getNormFactor(self):
662 662
663 663 pwcode = 1
664 664
665 665 if self.flagDecodeData:
666 666 pwcode = numpy.sum(self.code[0]**2)
667 667 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
668 668 normFactor = self.nProfiles*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
669 669
670 670 return normFactor
671 671
672 672 def getFlagCspc(self):
673 673
674 674 if self.data_cspc is None:
675 675 return True
676 676
677 677 return False
678 678
679 679 def getFlagDc(self):
680 680
681 681 if self.data_dc is None:
682 682 return True
683 683
684 684 return False
685 685
686 686 def getTimeInterval(self):
687 687
688 688 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles
689 689
690 690 return timeInterval
691 691
692 692 def getPower(self):
693 693
694 694 factor = self.normFactor
695 695 z = self.data_spc/factor
696 696 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
697 697 avg = numpy.average(z, axis=1)
698 698
699 699 return 10*numpy.log10(avg)
700 700
701 701 def getCoherence(self, pairsList=None, phase=False):
702 702
703 703 z = []
704 704 if pairsList is None:
705 705 pairsIndexList = self.pairsIndexList
706 706 else:
707 707 pairsIndexList = []
708 708 for pair in pairsList:
709 709 if pair not in self.pairsList:
710 710 raise ValueError, "Pair %s is not in dataOut.pairsList" %(pair)
711 711 pairsIndexList.append(self.pairsList.index(pair))
712 712 for i in range(len(pairsIndexList)):
713 713 pair = self.pairsList[pairsIndexList[i]]
714 714 ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
715 715 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
716 716 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
717 717 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
718 718 if phase:
719 719 data = numpy.arctan2(avgcoherenceComplex.imag,
720 720 avgcoherenceComplex.real)*180/numpy.pi
721 721 else:
722 722 data = numpy.abs(avgcoherenceComplex)
723 723
724 724 z.append(data)
725 725
726 726 return numpy.array(z)
727 727
728 728 def setValue(self, value):
729 729
730 730 print "This property should not be initialized"
731 731
732 732 return
733 733
734 734 nPairs = property(getNPairs, setValue, "I'm the 'nPairs' property.")
735 735 pairsIndexList = property(getPairsIndexList, setValue, "I'm the 'pairsIndexList' property.")
736 736 normFactor = property(getNormFactor, setValue, "I'm the 'getNormFactor' property.")
737 737 flag_cspc = property(getFlagCspc, setValue)
738 738 flag_dc = property(getFlagDc, setValue)
739 739 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
740 740 timeInterval = property(getTimeInterval, setValue, "I'm the 'timeInterval' property")
741 741
742 742 class SpectraHeis(Spectra):
743 743
744 744 data_spc = None
745 745
746 746 data_cspc = None
747 747
748 748 data_dc = None
749 749
750 750 nFFTPoints = None
751 751
752 752 # nPairs = None
753 753
754 754 pairsList = None
755 755
756 756 nCohInt = None
757 757
758 758 nIncohInt = None
759 759
760 760 def __init__(self):
761 761
762 762 self.radarControllerHeaderObj = RadarControllerHeader()
763 763
764 764 self.systemHeaderObj = SystemHeader()
765 765
766 766 self.type = "SpectraHeis"
767 767
768 768 # self.dtype = None
769 769
770 770 # self.nChannels = 0
771 771
772 772 # self.nHeights = 0
773 773
774 774 self.nProfiles = None
775 775
776 776 self.heightList = None
777 777
778 778 self.channelList = None
779 779
780 780 # self.channelIndexList = None
781 781
782 782 self.flagNoData = True
783 783
784 784 self.flagDiscontinuousBlock = False
785 785
786 786 # self.nPairs = 0
787 787
788 788 self.utctime = None
789 789
790 790 self.blocksize = None
791 791
792 792 self.profileIndex = 0
793 793
794 794 self.nCohInt = 1
795 795
796 796 self.nIncohInt = 1
797 797
798 798 def getNormFactor(self):
799 799 pwcode = 1
800 800 if self.flagDecodeData:
801 801 pwcode = numpy.sum(self.code[0]**2)
802 802
803 803 normFactor = self.nIncohInt*self.nCohInt*pwcode
804 804
805 805 return normFactor
806 806
807 807 def getTimeInterval(self):
808 808
809 809 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
810 810
811 811 return timeInterval
812 812
813 813 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
814 814 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
815 815
816 816 class Fits(JROData):
817 817
818 818 heightList = None
819 819
820 820 channelList = None
821 821
822 822 flagNoData = True
823 823
824 824 flagDiscontinuousBlock = False
825 825
826 826 useLocalTime = False
827 827
828 828 utctime = None
829 829
830 830 timeZone = None
831 831
832 832 # ippSeconds = None
833 833
834 834 # timeInterval = None
835 835
836 836 nCohInt = None
837 837
838 838 nIncohInt = None
839 839
840 840 noise = None
841 841
842 842 windowOfFilter = 1
843 843
844 844 #Speed of ligth
845 845 C = 3e8
846 846
847 847 frequency = 49.92e6
848 848
849 849 realtime = False
850 850
851 851
852 852 def __init__(self):
853 853
854 854 self.type = "Fits"
855 855
856 856 self.nProfiles = None
857 857
858 858 self.heightList = None
859 859
860 860 self.channelList = None
861 861
862 862 # self.channelIndexList = None
863 863
864 864 self.flagNoData = True
865 865
866 866 self.utctime = None
867 867
868 868 self.nCohInt = 1
869 869
870 870 self.nIncohInt = 1
871 871
872 872 self.useLocalTime = True
873 873
874 874 self.profileIndex = 0
875 875
876 876 # self.utctime = None
877 877 # self.timeZone = None
878 878 # self.ltctime = None
879 879 # self.timeInterval = None
880 880 # self.header = None
881 881 # self.data_header = None
882 882 # self.data = None
883 883 # self.datatime = None
884 884 # self.flagNoData = False
885 885 # self.expName = ''
886 886 # self.nChannels = None
887 887 # self.nSamples = None
888 888 # self.dataBlocksPerFile = None
889 889 # self.comments = ''
890 890 #
891 891
892 892
893 893 def getltctime(self):
894 894
895 895 if self.useLocalTime:
896 896 return self.utctime - self.timeZone*60
897 897
898 898 return self.utctime
899 899
900 900 def getDatatime(self):
901 901
902 902 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
903 903 return datatime
904 904
905 905 def getTimeRange(self):
906 906
907 907 datatime = []
908 908
909 909 datatime.append(self.ltctime)
910 910 datatime.append(self.ltctime + self.timeInterval)
911 911
912 912 datatime = numpy.array(datatime)
913 913
914 914 return datatime
915 915
916 916 def getHeiRange(self):
917 917
918 918 heis = self.heightList
919 919
920 920 return heis
921 921
922 922 def getNHeights(self):
923 923
924 924 return len(self.heightList)
925 925
926 926 def getNChannels(self):
927 927
928 928 return len(self.channelList)
929 929
930 930 def getChannelIndexList(self):
931 931
932 932 return range(self.nChannels)
933 933
934 934 def getNoise(self, type = 1):
935 935
936 936 #noise = numpy.zeros(self.nChannels)
937 937
938 938 if type == 1:
939 939 noise = self.getNoisebyHildebrand()
940 940
941 941 if type == 2:
942 942 noise = self.getNoisebySort()
943 943
944 944 if type == 3:
945 945 noise = self.getNoisebyWindow()
946 946
947 947 return noise
948 948
949 949 def getTimeInterval(self):
950 950
951 951 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
952 952
953 953 return timeInterval
954 954
955 955 datatime = property(getDatatime, "I'm the 'datatime' property")
956 956 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
957 957 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
958 958 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
959 959 noise = property(getNoise, "I'm the 'nHeights' property.")
960 960
961 961 ltctime = property(getltctime, "I'm the 'ltctime' property")
962 962 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
963 963
964 964
965 965 class Correlation(JROData):
966 966
967 967 noise = None
968 968
969 969 SNR = None
970 970
971 971 #--------------------------------------------------
972 972
973 973 mode = None
974 974
975 975 split = False
976 976
977 977 data_cf = None
978 978
979 979 lags = None
980 980
981 981 lagRange = None
982 982
983 983 pairsList = None
984 984
985 985 normFactor = None
986 986
987 987 #--------------------------------------------------
988 988
989 989 # calculateVelocity = None
990 990
991 991 nLags = None
992 992
993 993 nPairs = None
994 994
995 995 nAvg = None
996 996
997 997
998 998 def __init__(self):
999 999 '''
1000 1000 Constructor
1001 1001 '''
1002 1002 self.radarControllerHeaderObj = RadarControllerHeader()
1003 1003
1004 1004 self.systemHeaderObj = SystemHeader()
1005 1005
1006 1006 self.type = "Correlation"
1007 1007
1008 1008 self.data = None
1009 1009
1010 1010 self.dtype = None
1011 1011
1012 1012 self.nProfiles = None
1013 1013
1014 1014 self.heightList = None
1015 1015
1016 1016 self.channelList = None
1017 1017
1018 1018 self.flagNoData = True
1019 1019
1020 1020 self.flagDiscontinuousBlock = False
1021 1021
1022 1022 self.utctime = None
1023 1023
1024 1024 self.timeZone = None
1025 1025
1026 1026 self.dstFlag = None
1027 1027
1028 1028 self.errorCount = None
1029 1029
1030 1030 self.blocksize = None
1031 1031
1032 1032 self.flagDecodeData = False #asumo q la data no esta decodificada
1033 1033
1034 1034 self.flagDeflipData = False #asumo q la data no esta sin flip
1035 1035
1036 1036 self.pairsList = None
1037 1037
1038 1038 self.nPoints = None
1039 1039
1040 1040 def getPairsList(self):
1041 1041
1042 1042 return self.pairsList
1043 1043
1044 1044 def getNoise(self, mode = 2):
1045 1045
1046 1046 indR = numpy.where(self.lagR == 0)[0][0]
1047 1047 indT = numpy.where(self.lagT == 0)[0][0]
1048 1048
1049 1049 jspectra0 = self.data_corr[:,:,indR,:]
1050 1050 jspectra = copy.copy(jspectra0)
1051 1051
1052 1052 num_chan = jspectra.shape[0]
1053 1053 num_hei = jspectra.shape[2]
1054 1054
1055 1055 freq_dc = jspectra.shape[1]/2
1056 1056 ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
1057 1057
1058 1058 if ind_vel[0]<0:
1059 1059 ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
1060 1060
1061 1061 if mode == 1:
1062 1062 jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
1063 1063
1064 1064 if mode == 2:
1065 1065
1066 1066 vel = numpy.array([-2,-1,1,2])
1067 1067 xx = numpy.zeros([4,4])
1068 1068
1069 1069 for fil in range(4):
1070 1070 xx[fil,:] = vel[fil]**numpy.asarray(range(4))
1071 1071
1072 1072 xx_inv = numpy.linalg.inv(xx)
1073 1073 xx_aux = xx_inv[0,:]
1074 1074
1075 1075 for ich in range(num_chan):
1076 1076 yy = jspectra[ich,ind_vel,:]
1077 1077 jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy)
1078 1078
1079 1079 junkid = jspectra[ich,freq_dc,:]<=0
1080 1080 cjunkid = sum(junkid)
1081 1081
1082 1082 if cjunkid.any():
1083 1083 jspectra[ich,freq_dc,junkid.nonzero()] = (jspectra[ich,ind_vel[1],junkid] + jspectra[ich,ind_vel[2],junkid])/2
1084 1084
1085 1085 noise = jspectra0[:,freq_dc,:] - jspectra[:,freq_dc,:]
1086 1086
1087 1087 return noise
1088 1088
1089 1089 def getTimeInterval(self):
1090 1090
1091 1091 timeInterval = self.ippSeconds * self.nCohInt * self.nProfiles
1092 1092
1093 1093 return timeInterval
1094 1094
1095 1095 def splitFunctions(self):
1096 1096
1097 1097 pairsList = self.pairsList
1098 1098 ccf_pairs = []
1099 1099 acf_pairs = []
1100 1100 ccf_ind = []
1101 1101 acf_ind = []
1102 1102 for l in range(len(pairsList)):
1103 1103 chan0 = pairsList[l][0]
1104 1104 chan1 = pairsList[l][1]
1105 1105
1106 1106 #Obteniendo pares de Autocorrelacion
1107 1107 if chan0 == chan1:
1108 1108 acf_pairs.append(chan0)
1109 1109 acf_ind.append(l)
1110 1110 else:
1111 1111 ccf_pairs.append(pairsList[l])
1112 1112 ccf_ind.append(l)
1113 1113
1114 1114 data_acf = self.data_cf[acf_ind]
1115 1115 data_ccf = self.data_cf[ccf_ind]
1116 1116
1117 1117 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
1118 1118
1119 1119 def getNormFactor(self):
1120 1120 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
1121 1121 acf_pairs = numpy.array(acf_pairs)
1122 1122 normFactor = numpy.zeros((self.nPairs,self.nHeights))
1123 1123
1124 1124 for p in range(self.nPairs):
1125 1125 pair = self.pairsList[p]
1126 1126
1127 1127 ch0 = pair[0]
1128 1128 ch1 = pair[1]
1129 1129
1130 1130 ch0_max = numpy.max(data_acf[acf_pairs==ch0,:,:], axis=1)
1131 1131 ch1_max = numpy.max(data_acf[acf_pairs==ch1,:,:], axis=1)
1132 1132 normFactor[p,:] = numpy.sqrt(ch0_max*ch1_max)
1133 1133
1134 1134 return normFactor
1135 1135
1136 1136 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
1137 1137 normFactor = property(getNormFactor, "I'm the 'normFactor property'")
1138 1138
1139 1139 class Parameters(Spectra):
1140 1140
1141 1141 experimentInfo = None #Information about the experiment
1142 1142
1143 1143 #Information from previous data
1144 1144
1145 1145 inputUnit = None #Type of data to be processed
1146 1146
1147 1147 operation = None #Type of operation to parametrize
1148 1148
1149 1149 #normFactor = None #Normalization Factor
1150 1150
1151 1151 groupList = None #List of Pairs, Groups, etc
1152 1152
1153 1153 #Parameters
1154 1154
1155 1155 data_param = None #Parameters obtained
1156 1156
1157 1157 data_pre = None #Data Pre Parametrization
1158 1158
1159 1159 data_SNR = None #Signal to Noise Ratio
1160 1160
1161 1161 # heightRange = None #Heights
1162 1162
1163 1163 abscissaList = None #Abscissa, can be velocities, lags or time
1164 1164
1165 1165 # noise = None #Noise Potency
1166 1166
1167 1167 utctimeInit = None #Initial UTC time
1168 1168
1169 1169 paramInterval = None #Time interval to calculate Parameters in seconds
1170 1170
1171 1171 useLocalTime = True
1172 1172
1173 1173 #Fitting
1174 1174
1175 1175 data_error = None #Error of the estimation
1176 1176
1177 1177 constants = None
1178 1178
1179 1179 library = None
1180 1180
1181 1181 #Output signal
1182 1182
1183 1183 outputInterval = None #Time interval to calculate output signal in seconds
1184 1184
1185 1185 data_output = None #Out signal
1186 1186
1187 1187 nAvg = None
1188 1188
1189 1189 noise_estimation = None
1190 1190
1191 1191
1192 1192 def __init__(self):
1193 1193 '''
1194 1194 Constructor
1195 1195 '''
1196 1196 self.radarControllerHeaderObj = RadarControllerHeader()
1197 1197
1198 1198 self.systemHeaderObj = SystemHeader()
1199 1199
1200 1200 self.type = "Parameters"
1201 1201
1202 1202 def getTimeRange1(self, interval):
1203 1203
1204 1204 datatime = []
1205 1205
1206 1206 if self.useLocalTime:
1207 1207 time1 = self.utctimeInit - self.timeZone*60
1208 1208 else:
1209 1209 time1 = self.utctimeInit
1210 1210
1211 1211 datatime.append(time1)
1212 1212 datatime.append(time1 + interval)
1213 1213 datatime = numpy.array(datatime)
1214 1214
1215 1215 return datatime
1216 1216
1217 1217 def getTimeInterval(self):
1218 1218
1219 return self.timeInterval1
1219 if hasattr(self, 'timeInterval1'):
1220 return self.timeInterval1
1221 else:
1222 return self.paramInterval
1220 1223
1221 1224 def getNoise(self):
1222 1225
1223 1226 return self.spc_noise
1224 1227
1225 1228 timeInterval = property(getTimeInterval)
Show file before | Show file after | Hide comments
@@ -1,693 +1,756
1 1
2 2 import os
3 3 import zmq
4 4 import time
5 5 import numpy
6 6 import datetime
7 7 import numpy as np
8 8 import matplotlib
9 9 matplotlib.use('TkAgg')
10 10 import matplotlib.pyplot as plt
11 11 from mpl_toolkits.axes_grid1 import make_axes_locatable
12 12 from matplotlib.ticker import FuncFormatter, LinearLocator
13 13 from multiprocessing import Process
14 14
15 15 from schainpy.model.proc.jroproc_base import Operation
16 16
17 17 plt.ioff()
18 18
19 19 func = lambda x, pos: ('%s') %(datetime.datetime.fromtimestamp(x).strftime('%H:%M'))
20 20
21 21 d1970 = datetime.datetime(1970,1,1)
22 22
23 23 class PlotData(Operation, Process):
24 24
25 25 CODE = 'Figure'
26 26 colormap = 'jro'
27 27 CONFLATE = True
28 28 __MAXNUMX = 80
29 29 __MAXNUMY = 80
30 30 __missing = 1E30
31 31
32 32 def __init__(self, **kwargs):
33 33
34 34 Operation.__init__(self, plot=True, **kwargs)
35 35 Process.__init__(self)
36 36 self.kwargs['code'] = self.CODE
37 37 self.mp = False
38 38 self.dataOut = None
39 39 self.isConfig = False
40 40 self.figure = None
41 41 self.axes = []
42 42 self.localtime = kwargs.pop('localtime', True)
43 43 self.show = kwargs.get('show', True)
44 44 self.save = kwargs.get('save', False)
45 45 self.colormap = kwargs.get('colormap', self.colormap)
46 46 self.colormap_coh = kwargs.get('colormap_coh', 'jet')
47 47 self.colormap_phase = kwargs.get('colormap_phase', 'RdBu_r')
48 48 self.showprofile = kwargs.get('showprofile', True)
49 49 self.title = kwargs.get('wintitle', '')
50 50 self.xaxis = kwargs.get('xaxis', 'frequency')
51 51 self.zmin = kwargs.get('zmin', None)
52 52 self.zmax = kwargs.get('zmax', None)
53 53 self.xmin = kwargs.get('xmin', None)
54 54 self.xmax = kwargs.get('xmax', None)
55 55 self.xrange = kwargs.get('xrange', 24)
56 56 self.ymin = kwargs.get('ymin', None)
57 57 self.ymax = kwargs.get('ymax', None)
58 58 self.throttle_value = 5
59 59
60 60 def fill_gaps(self, x_buffer, y_buffer, z_buffer):
61 61
62 62 if x_buffer.shape[0] < 2:
63 63 return x_buffer, y_buffer, z_buffer
64 64
65 65 deltas = x_buffer[1:] - x_buffer[0:-1]
66 66 x_median = np.median(deltas)
67 67
68 68 index = np.where(deltas > 5*x_median)
69 69
70 70 if len(index[0]) != 0:
71 71 z_buffer[::, index[0], ::] = self.__missing
72 72 z_buffer = np.ma.masked_inside(z_buffer,
73 73 0.99*self.__missing,
74 74 1.01*self.__missing)
75 75
76 76 return x_buffer, y_buffer, z_buffer
77 77
78 78 def decimate(self):
79 79
80 80 # dx = int(len(self.x)/self.__MAXNUMX) + 1
81 81 dy = int(len(self.y)/self.__MAXNUMY) + 1
82 82
83 83 # x = self.x[::dx]
84 84 x = self.x
85 85 y = self.y[::dy]
86 86 z = self.z[::, ::, ::dy]
87 87
88 88 return x, y, z
89 89
90 90 def __plot(self):
91 91
92 92 print 'plotting...{}'.format(self.CODE)
93 93
94 94 if self.show:
95 print 'showing'
96 95 self.figure.show()
97 96
98 97 self.plot()
99 98 plt.tight_layout()
100 self.figure.canvas.manager.set_window_title('{} {} - Date:{}'.format(self.title, self.CODE.upper(),
101 datetime.datetime.fromtimestamp(self.max_time).strftime('%y/%m/%d %H:%M:%S')))
99 self.figure.canvas.manager.set_window_title('{} {} - {}'.format(self.title, self.CODE.upper(),
100 datetime.datetime.fromtimestamp(self.max_time).strftime('%Y/%m/%d')))
102 101
103 102 if self.save:
104 103 figname = os.path.join(self.save, '{}_{}.png'.format(self.CODE,
105 104 datetime.datetime.fromtimestamp(self.saveTime).strftime('%y%m%d_%H%M%S')))
106 105 print 'Saving figure: {}'.format(figname)
107 106 self.figure.savefig(figname)
108 107
109 108 self.figure.canvas.draw()
110 109
111 110 def plot(self):
112 111
113 112 print 'plotting...{}'.format(self.CODE.upper())
114 113 return
115 114
116 115 def run(self):
117 116
118 117 print '[Starting] {}'.format(self.name)
118
119 119 context = zmq.Context()
120 120 receiver = context.socket(zmq.SUB)
121 121 receiver.setsockopt(zmq.SUBSCRIBE, '')
122 122 receiver.setsockopt(zmq.CONFLATE, self.CONFLATE)
123 receiver.connect("ipc:///tmp/zmq.plots")
123 if 'server' in self.kwargs['parent']:
124 receiver.connect('ipc:///tmp/{}.plots'.format(self.kwargs['parent']['server']))
125 else:
126 receiver.connect("ipc:///tmp/zmq.plots")
124 127
125 128 while True:
126 129 try:
127 130 self.data = receiver.recv_pyobj(flags=zmq.NOBLOCK)
128 131 self.dataOut = self.data['dataOut']
129 132 self.times = self.data['times']
130 133 self.times.sort()
131 134 self.throttle_value = self.data['throttle']
132 135 self.min_time = self.times[0]
133 136 self.max_time = self.times[-1]
134 137
135 138 if self.isConfig is False:
136 139 self.setup()
137 140 self.isConfig = True
138 141 self.__plot()
139 142
140 143 if self.data['ENDED'] is True:
141 144 self.isConfig = False
142 145
143 146 except zmq.Again as e:
144 147 print 'Waiting for data...'
145 148 plt.pause(self.throttle_value)
146 149
147 150 def close(self):
148 151 if self.dataOut:
149 152 self.__plot()
150 153
151 154
152 155 class PlotSpectraData(PlotData):
153 156
154 157 CODE = 'spc'
155 158 colormap = 'jro'
156 159 CONFLATE = False
157 160
158 161 def setup(self):
159 162
160 163 ncolspan = 1
161 164 colspan = 1
162 165 self.ncols = int(numpy.sqrt(self.dataOut.nChannels)+0.9)
163 166 self.nrows = int(self.dataOut.nChannels*1./self.ncols + 0.9)
164 167 self.width = 3.6*self.ncols
165 168 self.height = 3.2*self.nrows
166 169 if self.showprofile:
167 170 ncolspan = 3
168 171 colspan = 2
169 172 self.width += 1.2*self.ncols
170 173
171 174 self.ylabel = 'Range [Km]'
172 175 self.titles = ['Channel {}'.format(x) for x in self.dataOut.channelList]
173 176
174 177 if self.figure is None:
175 178 self.figure = plt.figure(figsize=(self.width, self.height),
176 179 edgecolor='k',
177 180 facecolor='w')
178 181 else:
179 182 self.figure.clf()
180 183
181 184 n = 0
182 185 for y in range(self.nrows):
183 186 for x in range(self.ncols):
184 187 if n >= self.dataOut.nChannels:
185 188 break
186 189 ax = plt.subplot2grid((self.nrows, self.ncols*ncolspan), (y, x*ncolspan), 1, colspan)
187 190 if self.showprofile:
188 191 ax.ax_profile = plt.subplot2grid((self.nrows, self.ncols*ncolspan), (y, x*ncolspan+colspan), 1, 1)
189 192
190 193 ax.firsttime = True
191 194 self.axes.append(ax)
192 195 n += 1
193 196
194 197 def plot(self):
195 198
196 199 if self.xaxis == "frequency":
197 200 x = self.dataOut.getFreqRange(1)/1000.
198 201 xlabel = "Frequency (kHz)"
199 202 elif self.xaxis == "time":
200 203 x = self.dataOut.getAcfRange(1)
201 204 xlabel = "Time (ms)"
202 205 else:
203 206 x = self.dataOut.getVelRange(1)
204 207 xlabel = "Velocity (m/s)"
205 208
206 209 y = self.dataOut.getHeiRange()
207 210 z = self.data[self.CODE]
208 211
209 212 for n, ax in enumerate(self.axes):
210 213
211 214 if ax.firsttime:
212 215 self.xmax = self.xmax if self.xmax else np.nanmax(x)
213 216 self.xmin = self.xmin if self.xmin else -self.xmax
214 217 self.ymin = self.ymin if self.ymin else np.nanmin(y)
215 218 self.ymax = self.ymax if self.ymax else np.nanmax(y)
216 219 self.zmin = self.zmin if self.zmin else np.nanmin(z)
217 220 self.zmax = self.zmax if self.zmax else np.nanmax(z)
218 221 ax.plot = ax.pcolormesh(x, y, z[n].T,
219 222 vmin=self.zmin,
220 223 vmax=self.zmax,
221 224 cmap=plt.get_cmap(self.colormap)
222 225 )
223 226 divider = make_axes_locatable(ax)
224 227 cax = divider.new_horizontal(size='3%', pad=0.05)
225 228 self.figure.add_axes(cax)
226 229 plt.colorbar(ax.plot, cax)
227 230
228 231 ax.set_xlim(self.xmin, self.xmax)
229 232 ax.set_ylim(self.ymin, self.ymax)
230 233
231 234 ax.set_ylabel(self.ylabel)
232 235 ax.set_xlabel(xlabel)
233 236
234 237 ax.firsttime = False
235 238
236 239 if self.showprofile:
237 240 ax.plot_profile= ax.ax_profile.plot(self.data['rti'][self.max_time][n], y)[0]
238 241 ax.ax_profile.set_xlim(self.zmin, self.zmax)
239 242 ax.ax_profile.set_ylim(self.ymin, self.ymax)
240 243 ax.ax_profile.set_xlabel('dB')
241 244 ax.ax_profile.grid(b=True, axis='x')
242 245 ax.plot_noise = ax.ax_profile.plot(numpy.repeat(self.data['noise'][self.max_time][n], len(y)), y,
243 246 color="k", linestyle="dashed", lw=2)[0]
244 247 [tick.set_visible(False) for tick in ax.ax_profile.get_yticklabels()]
245 248 else:
246 249 ax.plot.set_array(z[n].T.ravel())
247 250 if self.showprofile:
248 251 ax.plot_profile.set_data(self.data['rti'][self.max_time][n], y)
249 252 ax.plot_noise.set_data(numpy.repeat(self.data['noise'][self.max_time][n], len(y)), y)
250 253
251 254 ax.set_title('{} - Noise: {:.2f} dB'.format(self.titles[n], self.data['noise'][self.max_time][n]),
252 255 size=8)
253 256 self.saveTime = self.max_time
254 257
255 258
256 259 class PlotCrossSpectraData(PlotData):
257 260
258 261 CODE = 'cspc'
259 262 zmin_coh = None
260 263 zmax_coh = None
261 264 zmin_phase = None
262 265 zmax_phase = None
263 266 CONFLATE = False
264 267
265 268 def setup(self):
266 269
267 270 ncolspan = 1
268 271 colspan = 1
269 272 self.ncols = 2
270 273 self.nrows = self.dataOut.nPairs
271 274 self.width = 3.6*self.ncols
272 275 self.height = 3.2*self.nrows
273 276
274 277 self.ylabel = 'Range [Km]'
275 278 self.titles = ['Channel {}'.format(x) for x in self.dataOut.channelList]
276 279
277 280 if self.figure is None:
278 281 self.figure = plt.figure(figsize=(self.width, self.height),
279 282 edgecolor='k',
280 283 facecolor='w')
281 284 else:
282 285 self.figure.clf()
283 286
284 287 for y in range(self.nrows):
285 288 for x in range(self.ncols):
286 289 ax = plt.subplot2grid((self.nrows, self.ncols), (y, x), 1, 1)
287 290 ax.firsttime = True
288 291 self.axes.append(ax)
289 292
290 293 def plot(self):
291 294
292 295 if self.xaxis == "frequency":
293 296 x = self.dataOut.getFreqRange(1)/1000.
294 297 xlabel = "Frequency (kHz)"
295 298 elif self.xaxis == "time":
296 299 x = self.dataOut.getAcfRange(1)
297 300 xlabel = "Time (ms)"
298 301 else:
299 302 x = self.dataOut.getVelRange(1)
300 303 xlabel = "Velocity (m/s)"
301 304
302 305 y = self.dataOut.getHeiRange()
303 306 z_coh = self.data['cspc_coh']
304 307 z_phase = self.data['cspc_phase']
305 308
306 309 for n in range(self.nrows):
307 310 ax = self.axes[2*n]
308 311 ax1 = self.axes[2*n+1]
309 312 if ax.firsttime:
310 313 self.xmax = self.xmax if self.xmax else np.nanmax(x)
311 314 self.xmin = self.xmin if self.xmin else -self.xmax
312 315 self.ymin = self.ymin if self.ymin else np.nanmin(y)
313 316 self.ymax = self.ymax if self.ymax else np.nanmax(y)
314 317 self.zmin_coh = self.zmin_coh if self.zmin_coh else 0.0
315 318 self.zmax_coh = self.zmax_coh if self.zmax_coh else 1.0
316 319 self.zmin_phase = self.zmin_phase if self.zmin_phase else -180
317 320 self.zmax_phase = self.zmax_phase if self.zmax_phase else 180
318 321
319 322 ax.plot = ax.pcolormesh(x, y, z_coh[n].T,
320 323 vmin=self.zmin_coh,
321 324 vmax=self.zmax_coh,
322 325 cmap=plt.get_cmap(self.colormap_coh)
323 326 )
324 327 divider = make_axes_locatable(ax)
325 328 cax = divider.new_horizontal(size='3%', pad=0.05)
326 329 self.figure.add_axes(cax)
327 330 plt.colorbar(ax.plot, cax)
328 331
329 332 ax.set_xlim(self.xmin, self.xmax)
330 333 ax.set_ylim(self.ymin, self.ymax)
331 334
332 335 ax.set_ylabel(self.ylabel)
333 336 ax.set_xlabel(xlabel)
334 337 ax.firsttime = False
335 338
336 339 ax1.plot = ax1.pcolormesh(x, y, z_phase[n].T,
337 340 vmin=self.zmin_phase,
338 341 vmax=self.zmax_phase,
339 342 cmap=plt.get_cmap(self.colormap_phase)
340 343 )
341 344 divider = make_axes_locatable(ax1)
342 345 cax = divider.new_horizontal(size='3%', pad=0.05)
343 346 self.figure.add_axes(cax)
344 347 plt.colorbar(ax1.plot, cax)
345 348
346 349 ax1.set_xlim(self.xmin, self.xmax)
347 350 ax1.set_ylim(self.ymin, self.ymax)
348 351
349 352 ax1.set_ylabel(self.ylabel)
350 353 ax1.set_xlabel(xlabel)
351 354 ax1.firsttime = False
352 355 else:
353 356 ax.plot.set_array(z_coh[n].T.ravel())
354 357 ax1.plot.set_array(z_phase[n].T.ravel())
355 358
356 359 ax.set_title('Coherence Ch{} * Ch{}'.format(self.dataOut.pairsList[n][0], self.dataOut.pairsList[n][1]), size=8)
357 360 ax1.set_title('Phase Ch{} * Ch{}'.format(self.dataOut.pairsList[n][0], self.dataOut.pairsList[n][1]), size=8)
358 361 self.saveTime = self.max_time
359 362
360 363
361 364 class PlotSpectraMeanData(PlotSpectraData):
362 365
363 366 CODE = 'spc_mean'
364 367 colormap = 'jet'
365 368
366 369 def plot(self):
367 370
368 371 if self.xaxis == "frequency":
369 372 x = self.dataOut.getFreqRange(1)/1000.
370 373 xlabel = "Frequency (kHz)"
371 374 elif self.xaxis == "time":
372 375 x = self.dataOut.getAcfRange(1)
373 376 xlabel = "Time (ms)"
374 377 else:
375 378 x = self.dataOut.getVelRange(1)
376 379 xlabel = "Velocity (m/s)"
377 380
378 381 y = self.dataOut.getHeiRange()
379 382 z = self.data['spc']
380 383 mean = self.data['mean'][self.max_time]
381 384
382 385 for n, ax in enumerate(self.axes):
383 386
384 387 if ax.firsttime:
385 388 self.xmax = self.xmax if self.xmax else np.nanmax(x)
386 389 self.xmin = self.xmin if self.xmin else -self.xmax
387 390 self.ymin = self.ymin if self.ymin else np.nanmin(y)
388 391 self.ymax = self.ymax if self.ymax else np.nanmax(y)
389 392 self.zmin = self.zmin if self.zmin else np.nanmin(z)
390 393 self.zmax = self.zmax if self.zmax else np.nanmax(z)
391 394 ax.plt = ax.pcolormesh(x, y, z[n].T,
392 395 vmin=self.zmin,
393 396 vmax=self.zmax,
394 397 cmap=plt.get_cmap(self.colormap)
395 398 )
396 399 ax.plt_dop = ax.plot(mean[n], y,
397 400 color='k')[0]
398 401
399 402 divider = make_axes_locatable(ax)
400 403 cax = divider.new_horizontal(size='3%', pad=0.05)
401 404 self.figure.add_axes(cax)
402 405 plt.colorbar(ax.plt, cax)
403 406
404 407 ax.set_xlim(self.xmin, self.xmax)
405 408 ax.set_ylim(self.ymin, self.ymax)
406 409
407 410 ax.set_ylabel(self.ylabel)
408 411 ax.set_xlabel(xlabel)
409 412
410 413 ax.firsttime = False
411 414
412 415 if self.showprofile:
413 416 ax.plt_profile= ax.ax_profile.plot(self.data['rti'][self.max_time][n], y)[0]
414 417 ax.ax_profile.set_xlim(self.zmin, self.zmax)
415 418 ax.ax_profile.set_ylim(self.ymin, self.ymax)
416 419 ax.ax_profile.set_xlabel('dB')
417 420 ax.ax_profile.grid(b=True, axis='x')
418 421 ax.plt_noise = ax.ax_profile.plot(numpy.repeat(self.data['noise'][self.max_time][n], len(y)), y,
419 422 color="k", linestyle="dashed", lw=2)[0]
420 423 [tick.set_visible(False) for tick in ax.ax_profile.get_yticklabels()]
421 424 else:
422 425 ax.plt.set_array(z[n].T.ravel())
423 426 ax.plt_dop.set_data(mean[n], y)
424 427 if self.showprofile:
425 428 ax.plt_profile.set_data(self.data['rti'][self.max_time][n], y)
426 429 ax.plt_noise.set_data(numpy.repeat(self.data['noise'][self.max_time][n], len(y)), y)
427 430
428 431 ax.set_title('{} - Noise: {:.2f} dB'.format(self.titles[n], self.data['noise'][self.max_time][n]),
429 432 size=8)
430 433 self.saveTime = self.max_time
431 434
432 435
433 436 class PlotRTIData(PlotData):
434 437
435 438 CODE = 'rti'
436 439 colormap = 'jro'
437 440
438 441 def setup(self):
439 442 self.ncols = 1
440 443 self.nrows = self.dataOut.nChannels
441 444 self.width = 10
442 445 self.height = 2.2*self.nrows if self.nrows<6 else 12
443 446 if self.nrows==1:
444 447 self.height += 1
445 448 self.ylabel = 'Range [Km]'
446 449 self.titles = ['Channel {}'.format(x) for x in self.dataOut.channelList]
447 450
448 451 if self.figure is None:
449 452 self.figure = plt.figure(figsize=(self.width, self.height),
450 453 edgecolor='k',
451 454 facecolor='w')
452 455 else:
453 456 self.figure.clf()
454 457 self.axes = []
455 458
456 459 for n in range(self.nrows):
457 460 ax = self.figure.add_subplot(self.nrows, self.ncols, n+1)
458 461 ax.firsttime = True
459 462 self.axes.append(ax)
460 463
461 464 def plot(self):
462 465
463 466 self.x = np.array(self.times)
464 467 self.y = self.dataOut.getHeiRange()
465 468 self.z = []
466 469
467 470 for ch in range(self.nrows):
468 471 self.z.append([self.data[self.CODE][t][ch] for t in self.times])
469 472
470 473 self.z = np.array(self.z)
471 474 for n, ax in enumerate(self.axes):
472 475
473 476 x, y, z = self.fill_gaps(*self.decimate())
474 477 xmin = self.min_time
475 478 xmax = xmin+self.xrange*60*60
476 479 if ax.firsttime:
477 480 self.ymin = self.ymin if self.ymin else np.nanmin(self.y)
478 481 self.ymax = self.ymax if self.ymax else np.nanmax(self.y)
479 482 self.zmin = self.zmin if self.zmin else np.nanmin(self.z)
480 483 self.zmax = self.zmax if self.zmax else np.nanmax(self.z)
481 484 plot = ax.pcolormesh(x, y, z[n].T,
482 485 vmin=self.zmin,
483 486 vmax=self.zmax,
484 487 cmap=plt.get_cmap(self.colormap)
485 488 )
486 489 divider = make_axes_locatable(ax)
487 490 cax = divider.new_horizontal(size='2%', pad=0.05)
488 491 self.figure.add_axes(cax)
489 492 plt.colorbar(plot, cax)
490 493 ax.set_ylim(self.ymin, self.ymax)
491 494
492 495 ax.xaxis.set_major_formatter(FuncFormatter(func))
493 496 ax.xaxis.set_major_locator(LinearLocator(6))
494 497
495 498 ax.set_ylabel(self.ylabel)
496 499
497 500 # if self.xmin is None:
498 501 # xmin = self.min_time
499 502 # else:
500 503 # xmin = (datetime.datetime.combine(self.dataOut.datatime.date(),
501 504 # datetime.time(self.xmin, 0, 0))-d1970).total_seconds()
502 505
503 506 ax.set_xlim(xmin, xmax)
504 507 ax.firsttime = False
505 508 else:
506 509 ax.collections.remove(ax.collections[0])
507 510 ax.set_xlim(xmin, xmax)
508 511 plot = ax.pcolormesh(x, y, z[n].T,
509 512 vmin=self.zmin,
510 513 vmax=self.zmax,
511 514 cmap=plt.get_cmap(self.colormap)
512 515 )
513 516 ax.set_title('{} {}'.format(self.titles[n],
514 517 datetime.datetime.fromtimestamp(self.max_time).strftime('%y/%m/%d %H:%M:%S')),
515 518 size=8)
516 519
517 520 self.saveTime = self.min_time
518 521
519 522
520 523 class PlotCOHData(PlotRTIData):
521 524
522 525 CODE = 'coh'
523 526
524 527 def setup(self):
525 528
526 529 self.ncols = 1
527 530 self.nrows = self.dataOut.nPairs
528 531 self.width = 10
529 532 self.height = 2.2*self.nrows if self.nrows<6 else 12
530 533 if self.nrows==1:
531 534 self.height += 1
532 535 self.ylabel = 'Range [Km]'
533 536 self.titles = ['{} Ch{} * Ch{}'.format(self.CODE.upper(), x[0], x[1]) for x in self.dataOut.pairsList]
534 537
535 538 if self.figure is None:
536 539 self.figure = plt.figure(figsize=(self.width, self.height),
537 540 edgecolor='k',
538 541 facecolor='w')
539 542 else:
540 543 self.figure.clf()
541 544 self.axes = []
542 545
543 546 for n in range(self.nrows):
544 547 ax = self.figure.add_subplot(self.nrows, self.ncols, n+1)
545 548 ax.firsttime = True
546 549 self.axes.append(ax)
547 550
548 551
549 552 class PlotNoiseData(PlotData):
550 553 CODE = 'noise'
551 554
552 555 def setup(self):
553 556
554 557 self.ncols = 1
555 558 self.nrows = 1
556 559 self.width = 10
557 560 self.height = 3.2
558 561 self.ylabel = 'Intensity [dB]'
559 562 self.titles = ['Noise']
560 563
561 564 if self.figure is None:
562 565 self.figure = plt.figure(figsize=(self.width, self.height),
563 566 edgecolor='k',
564 567 facecolor='w')
565 568 else:
566 569 self.figure.clf()
567 570 self.axes = []
568 571
569 572 self.ax = self.figure.add_subplot(self.nrows, self.ncols, 1)
570 573 self.ax.firsttime = True
571 574
572 575 def plot(self):
573 576
574 577 x = self.times
575 578 xmin = self.min_time
576 579 xmax = xmin+self.xrange*60*60
577 580 if self.ax.firsttime:
578 581 for ch in self.dataOut.channelList:
579 582 y = [self.data[self.CODE][t][ch] for t in self.times]
580 583 self.ax.plot(x, y, lw=1, label='Ch{}'.format(ch))
581 584 self.ax.firsttime = False
582 585 self.ax.xaxis.set_major_formatter(FuncFormatter(func))
583 586 self.ax.xaxis.set_major_locator(LinearLocator(6))
584 587 self.ax.set_ylabel(self.ylabel)
585 588 plt.legend()
586 589 else:
587 590 for ch in self.dataOut.channelList:
588 591 y = [self.data[self.CODE][t][ch] for t in self.times]
589 592 self.ax.lines[ch].set_data(x, y)
590 593
591 594 self.ax.set_xlim(xmin, xmax)
592 595 self.ax.set_ylim(min(y)-5, max(y)+5)
593 596 self.saveTime = self.min_time
594 597
595 598
596 599 class PlotWindProfilerData(PlotRTIData):
600
597 601 CODE = 'wind'
598 602 colormap = 'seismic'
599 603
600 604 def setup(self):
601 605 self.ncols = 1
602 606 self.nrows = self.dataOut.data_output.shape[0]
603 607 self.width = 10
604 608 self.height = 2.2*self.nrows
605 609 self.ylabel = 'Height [Km]'
606 self.titles = ['Zonal' ,'Meridional', 'Vertical']
610 self.titles = ['Zonal Wind' ,'Meridional Wind', 'Vertical Wind']
607 611 self.clabels = ['Velocity (m/s)','Velocity (m/s)','Velocity (cm/s)']
608 612 self.windFactor = [1, 1, 100]
609 613
610 614 if self.figure is None:
611 615 self.figure = plt.figure(figsize=(self.width, self.height),
612 616 edgecolor='k',
613 617 facecolor='w')
614 618 else:
615 619 self.figure.clf()
616 620 self.axes = []
617 621
618 622 for n in range(self.nrows):
619 623 ax = self.figure.add_subplot(self.nrows, self.ncols, n+1)
620 624 ax.firsttime = True
621 625 self.axes.append(ax)
622 626
623 627 def plot(self):
624 628
625 629 self.x = np.array(self.times)
626 630 self.y = self.dataOut.heightList
627 631 self.z = []
628 632
629 633 for ch in range(self.nrows):
630 self.z.append([self.data[self.CODE][t][ch] for t in self.times])
634 self.z.append([self.data['output'][t][ch] for t in self.times])
631 635
632 636 self.z = np.array(self.z)
633 637 self.z = numpy.ma.masked_invalid(self.z)
634 638
635 639 cmap=plt.get_cmap(self.colormap)
636 cmap.set_bad('white', 1.)
640 cmap.set_bad('black', 1.)
637 641
638 642 for n, ax in enumerate(self.axes):
639 643 x, y, z = self.fill_gaps(*self.decimate())
640 644 xmin = self.min_time
641 645 xmax = xmin+self.xrange*60*60
642 646 if ax.firsttime:
643 647 self.ymin = self.ymin if self.ymin else np.nanmin(self.y)
644 648 self.ymax = self.ymax if self.ymax else np.nanmax(self.y)
645 649 self.zmax = self.zmax if self.zmax else numpy.nanmax(abs(self.z[:-1, :]))
646 650 self.zmin = self.zmin if self.zmin else -self.zmax
647 651
648 652 plot = ax.pcolormesh(x, y, z[n].T*self.windFactor[n],
649 653 vmin=self.zmin,
650 654 vmax=self.zmax,
651 655 cmap=cmap
652 656 )
653 657 divider = make_axes_locatable(ax)
654 658 cax = divider.new_horizontal(size='2%', pad=0.05)
655 cax.set_ylabel(self.clabels[n])
656 659 self.figure.add_axes(cax)
657 plt.colorbar(plot, cax)
660 cb = plt.colorbar(plot, cax)
661 cb.set_label(self.clabels[n])
658 662 ax.set_ylim(self.ymin, self.ymax)
659 663
660 664 ax.xaxis.set_major_formatter(FuncFormatter(func))
661 665 ax.xaxis.set_major_locator(LinearLocator(6))
662 666
663 667 ax.set_ylabel(self.ylabel)
664 668
665 669 ax.set_xlim(xmin, xmax)
666 670 ax.firsttime = False
667 671 else:
668 672 ax.collections.remove(ax.collections[0])
669 673 ax.set_xlim(xmin, xmax)
670 674 plot = ax.pcolormesh(x, y, z[n].T*self.windFactor[n],
671 675 vmin=self.zmin,
672 676 vmax=self.zmax,
673 677 cmap=plt.get_cmap(self.colormap)
674 678 )
675 679 ax.set_title('{} {}'.format(self.titles[n],
676 680 datetime.datetime.fromtimestamp(self.max_time).strftime('%y/%m/%d %H:%M:%S')),
677 681 size=8)
678 682
679 683 self.saveTime = self.min_time
680 684
681 685
682 686 class PlotSNRData(PlotRTIData):
683 687 CODE = 'snr'
684 688 colormap = 'jet'
685 689
686 690 class PlotDOPData(PlotRTIData):
687 691 CODE = 'dop'
688 692 colormap = 'jet'
689 693
690 694
691 695 class PlotPHASEData(PlotCOHData):
692 696 CODE = 'phase'
693 697 colormap = 'seismic'
698
699
700 class PlotSkyMapData(PlotData):
701
702 CODE = 'met'
703
704 def setup(self):
705
706 self.ncols = 1
707 self.nrows = 1
708 self.width = 7.2
709 self.height = 7.2
710
711 self.xlabel = 'Zonal Zenith Angle (deg)'
712 self.ylabel = 'Meridional Zenith Angle (deg)'
713
714 if self.figure is None:
715 self.figure = plt.figure(figsize=(self.width, self.height),
716 edgecolor='k',
717 facecolor='w')
718 else:
719 self.figure.clf()
720
721 self.ax = plt.subplot2grid((self.nrows, self.ncols), (0, 0), 1, 1, polar=True)
722 self.ax.firsttime = True
723
724
725 def plot(self):
726
727 arrayParameters = np.concatenate([self.data['param'][t] for t in self.times])
728 error = arrayParameters[:,-1]
729 indValid = numpy.where(error == 0)[0]
730 finalMeteor = arrayParameters[indValid,:]
731 finalAzimuth = finalMeteor[:,3]
732 finalZenith = finalMeteor[:,4]
733
734 x = finalAzimuth*numpy.pi/180
735 y = finalZenith
736
737 if self.ax.firsttime:
738 self.ax.plot = self.ax.plot(x, y, 'bo', markersize=5)[0]
739 self.ax.set_ylim(0,90)
740 self.ax.set_yticks(numpy.arange(0,90,20))
741 self.ax.set_xlabel(self.xlabel)
742 self.ax.set_ylabel(self.ylabel)
743 self.ax.yaxis.labelpad = 40
744 self.ax.firsttime = False
745 else:
746 self.ax.plot.set_data(x, y)
747
748
749 dt1 = datetime.datetime.fromtimestamp(self.min_time).strftime('%y/%m/%d %H:%M:%S')
750 dt2 = datetime.datetime.fromtimestamp(self.max_time).strftime('%y/%m/%d %H:%M:%S')
751 title = 'Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n' % (dt1,
752 dt2,
753 len(x))
754 self.ax.set_title(title, size=8)
755
756 self.saveTime = self.max_time
Show file before | Show file after | Hide comments
@@ -1,315 +1,318
1 1 '''
2 2
3 3 $Author: murco $
4 4 $Id: jroproc_base.py 1 2012-11-12 18:56:07Z murco $
5 5 '''
6 6
7 7 class ProcessingUnit(object):
8 8
9 9 """
10 10 Esta es la clase base para el procesamiento de datos.
11 11
12 12 Contiene el metodo "call" para llamar operaciones. Las operaciones pueden ser:
13 13 - Metodos internos (callMethod)
14 14 - Objetos del tipo Operation (callObject). Antes de ser llamados, estos objetos
15 15 tienen que ser agreagados con el metodo "add".
16 16
17 17 """
18 18 # objeto de datos de entrada (Voltage, Spectra o Correlation)
19 19 dataIn = None
20 20 dataInList = []
21 21
22 22 # objeto de datos de entrada (Voltage, Spectra o Correlation)
23 23 dataOut = None
24 24
25 25 operations2RunDict = None
26 26
27 27 isConfig = False
28 28
29 29
30 30 def __init__(self, *args, **kwargs):
31 31
32 32 self.dataIn = None
33 33 self.dataInList = []
34 34
35 35 self.dataOut = None
36 36
37 37 self.operations2RunDict = {}
38 38 self.operationKwargs = {}
39 39
40 40 self.isConfig = False
41 41
42 42 self.args = args
43 43 self.kwargs = kwargs
44 44
45 45 def addOperationKwargs(self, objId, **kwargs):
46 46 '''
47 47 '''
48 48
49 49 self.operationKwargs[objId] = kwargs
50 50
51 51
52 52 def addOperation(self, opObj, objId):
53 53
54 54 """
55 55 Agrega un objeto del tipo "Operation" (opObj) a la lista de objetos "self.objectList" y retorna el
56 56 identificador asociado a este objeto.
57 57
58 58 Input:
59 59
60 60 object : objeto de la clase "Operation"
61 61
62 62 Return:
63 63
64 64 objId : identificador del objeto, necesario para ejecutar la operacion
65 65 """
66 66
67 67 self.operations2RunDict[objId] = opObj
68 68
69 69 return objId
70 70
71 71 def getOperationObj(self, objId):
72 72
73 73 if objId not in self.operations2RunDict.keys():
74 74 return None
75 75
76 76 return self.operations2RunDict[objId]
77 77
78 78 def operation(self, **kwargs):
79 79
80 80 """
81 81 Operacion directa sobre la data (dataOut.data). Es necesario actualizar los valores de los
82 82 atributos del objeto dataOut
83 83
84 84 Input:
85 85
86 86 **kwargs : Diccionario de argumentos de la funcion a ejecutar
87 87 """
88 88
89 89 raise NotImplementedError
90 90
91 91 def callMethod(self, name, opId):
92 92
93 93 """
94 94 Ejecuta el metodo con el nombre "name" y con argumentos **kwargs de la propia clase.
95 95
96 96 Input:
97 97 name : nombre del metodo a ejecutar
98 98
99 99 **kwargs : diccionario con los nombres y valores de la funcion a ejecutar.
100 100
101 101 """
102 102
103 103 #Checking the inputs
104 104 if name == 'run':
105 105
106 106 if not self.checkInputs():
107 107 self.dataOut.flagNoData = True
108 108 return False
109 109 else:
110 110 #Si no es un metodo RUN la entrada es la misma dataOut (interna)
111 111 if self.dataOut is not None and self.dataOut.isEmpty():
112 112 return False
113 113
114 114 #Getting the pointer to method
115 115 methodToCall = getattr(self, name)
116 116
117 117 #Executing the self method
118 118
119 119 if hasattr(self, 'mp'):
120 120 if name=='run':
121 121 if self.mp is False:
122 122 self.mp = True
123 123 self.start()
124 else:
124 else:
125 self.operationKwargs[opId]['parent'] = self.kwargs
125 126 methodToCall(**self.operationKwargs[opId])
126 127 else:
127 128 if name=='run':
128 129 methodToCall(**self.kwargs)
129 else:
130 else:
130 131 methodToCall(**self.operationKwargs[opId])
131 132
132 133 if self.dataOut is None:
133 134 return False
134 135
135 136 if self.dataOut.isEmpty():
136 137 return False
137 138
138 139 return True
139 140
140 141 def callObject(self, objId):
141 142
142 143 """
143 144 Ejecuta la operacion asociada al identificador del objeto "objId"
144 145
145 146 Input:
146 147
147 148 objId : identificador del objeto a ejecutar
148 149
149 150 **kwargs : diccionario con los nombres y valores de la funcion a ejecutar.
150 151
151 152 Return:
152 153
153 154 None
154 155 """
155 156
156 157 if self.dataOut is not None and self.dataOut.isEmpty():
157 158 return False
158 159
159 160 externalProcObj = self.operations2RunDict[objId]
160 161
161 162 if hasattr(externalProcObj, 'mp'):
162 163 if externalProcObj.mp is False:
164 externalProcObj.kwargs['parent'] = self.kwargs
163 165 self.operationKwargs[objId] = externalProcObj.kwargs
164 166 externalProcObj.mp = True
165 167 externalProcObj.start()
166 168 else:
167 169 externalProcObj.run(self.dataOut, **externalProcObj.kwargs)
168 170 self.operationKwargs[objId] = externalProcObj.kwargs
171
169 172
170 173 return True
171 174
172 175 def call(self, opType, opName=None, opId=None):
173 176
174 177 """
175 178 Return True si ejecuta la operacion interna nombrada "opName" o la operacion externa
176 179 identificada con el id "opId"; con los argumentos "**kwargs".
177 180
178 181 False si la operacion no se ha ejecutado.
179 182
180 183 Input:
181 184
182 185 opType : Puede ser "self" o "external"
183 186
184 187 Depende del tipo de operacion para llamar a:callMethod or callObject:
185 188
186 189 1. If opType = "self": Llama a un metodo propio de esta clase:
187 190
188 191 name_method = getattr(self, name)
189 192 name_method(**kwargs)
190 193
191 194
192 195 2. If opType = "other" o"external": Llama al metodo "run()" de una instancia de la
193 196 clase "Operation" o de un derivado de ella:
194 197
195 198 instanceName = self.operationList[opId]
196 199 instanceName.run(**kwargs)
197 200
198 201 opName : Si la operacion es interna (opType = 'self'), entonces el "opName" sera
199 202 usada para llamar a un metodo interno de la clase Processing
200 203
201 204 opId : Si la operacion es externa (opType = 'other' o 'external), entonces el
202 205 "opId" sera usada para llamar al metodo "run" de la clase Operation
203 206 registrada anteriormente con ese Id
204 207
205 208 Exception:
206 209 Este objeto de tipo Operation debe de haber sido agregado antes con el metodo:
207 210 "addOperation" e identificado con el valor "opId" = el id de la operacion.
208 211 De lo contrario retornara un error del tipo ValueError
209 212
210 213 """
211 214
212 215 if opType == 'self':
213 216
214 217 if not opName:
215 218 raise ValueError, "opName parameter should be defined"
216 219
217 220 sts = self.callMethod(opName, opId)
218 221
219 222 elif opType == 'other' or opType == 'external' or opType == 'plotter':
220 223
221 224 if not opId:
222 225 raise ValueError, "opId parameter should be defined"
223 226
224 227 if opId not in self.operations2RunDict.keys():
225 228 raise ValueError, "Any operation with id=%s has been added" %str(opId)
226 229
227 230 sts = self.callObject(opId)
228 231
229 232 else:
230 233 raise ValueError, "opType should be 'self', 'external' or 'plotter'; and not '%s'" %opType
231 234
232 235 return sts
233 236
234 237 def setInput(self, dataIn):
235 238
236 239 self.dataIn = dataIn
237 240 self.dataInList.append(dataIn)
238 241
239 242 def getOutputObj(self):
240 243
241 244 return self.dataOut
242 245
243 246 def checkInputs(self):
244 247
245 248 for thisDataIn in self.dataInList:
246 249
247 250 if thisDataIn.isEmpty():
248 251 return False
249 252
250 253 return True
251 254
252 255 def setup(self):
253 256
254 257 raise NotImplementedError
255 258
256 259 def run(self):
257 260
258 261 raise NotImplementedError
259 262
260 263 def close(self):
261 264 #Close every thread, queue or any other object here is it is neccesary.
262 265 return
263 266
264 267 class Operation(object):
265 268
266 269 """
267 270 Clase base para definir las operaciones adicionales que se pueden agregar a la clase ProcessingUnit
268 271 y necesiten acumular informacion previa de los datos a procesar. De preferencia usar un buffer de
269 272 acumulacion dentro de esta clase
270 273
271 274 Ejemplo: Integraciones coherentes, necesita la informacion previa de los n perfiles anteriores (bufffer)
272 275
273 276 """
274 277
275 278 __buffer = None
276 279 isConfig = False
277 280
278 281 def __init__(self, **kwargs):
279 282
280 283 self.__buffer = None
281 284 self.isConfig = False
282 285 self.kwargs = kwargs
283 286
284 287 def setup(self):
285 288
286 289 self.isConfig = True
287 290
288 291 raise NotImplementedError
289 292
290 293 def run(self, dataIn, **kwargs):
291 294
292 295 """
293 296 Realiza las operaciones necesarias sobre la dataIn.data y actualiza los
294 297 atributos del objeto dataIn.
295 298
296 299 Input:
297 300
298 301 dataIn : objeto del tipo JROData
299 302
300 303 Return:
301 304
302 305 None
303 306
304 307 Affected:
305 308 __buffer : buffer de recepcion de datos.
306 309
307 310 """
308 311 if not self.isConfig:
309 312 self.setup(**kwargs)
310 313
311 314 raise NotImplementedError
312 315
313 316 def close(self):
314 317
315 318 pass
Show file before | Show file after | Hide comments
@@ -1,445 +1,450
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 cPickle as pickle
11 11 import datetime
12 12 from zmq.utils.monitor import recv_monitor_message
13 13 from functools import wraps
14 14 from threading import Thread
15 15 from multiprocessing import Process
16 16
17 17 from schainpy.model.proc.jroproc_base import Operation, ProcessingUnit
18 18
19 19 MAXNUMX = 100
20 20 MAXNUMY = 100
21 21
22 22 class PrettyFloat(float):
23 23 def __repr__(self):
24 24 return '%.2f' % self
25 25
26 26 def roundFloats(obj):
27 27 if isinstance(obj, list):
28 28 return map(roundFloats, obj)
29 29 elif isinstance(obj, float):
30 30 return round(obj, 2)
31 31
32 32 def decimate(z):
33 33 # dx = int(len(self.x)/self.__MAXNUMX) + 1
34 34
35 35 dy = int(len(z[0])/MAXNUMY) + 1
36 36
37 37 return z[::, ::dy]
38 38
39 39 class throttle(object):
40 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 75
76 76 class PublishData(Operation):
77 77 """Clase publish."""
78 78
79 79 def __init__(self, **kwargs):
80 80 """Inicio."""
81 81 Operation.__init__(self, **kwargs)
82 82 self.isConfig = False
83 83 self.client = None
84 84 self.zeromq = None
85 85 self.mqtt = None
86 86
87 87 def on_disconnect(self, client, userdata, rc):
88 88 if rc != 0:
89 89 print("Unexpected disconnection.")
90 90 self.connect()
91 91
92 92 def connect(self):
93 93 print 'trying to connect'
94 94 try:
95 95 self.client.connect(
96 96 host=self.host,
97 97 port=self.port,
98 98 keepalive=60*10,
99 99 bind_address='')
100 100 self.client.loop_start()
101 101 # self.client.publish(
102 102 # self.topic + 'SETUP',
103 103 # json.dumps(setup),
104 104 # retain=True
105 105 # )
106 106 except:
107 107 print "MQTT Conection error."
108 108 self.client = False
109 109
110 110 def setup(self, port=1883, username=None, password=None, clientId="user", zeromq=1, **kwargs):
111 111 self.counter = 0
112 112 self.topic = kwargs.get('topic', 'schain')
113 113 self.delay = kwargs.get('delay', 0)
114 114 self.plottype = kwargs.get('plottype', 'spectra')
115 115 self.host = kwargs.get('host', "10.10.10.82")
116 116 self.port = kwargs.get('port', 3000)
117 117 self.clientId = clientId
118 118 self.cnt = 0
119 119 self.zeromq = zeromq
120 120 self.mqtt = kwargs.get('plottype', 0)
121 121 self.client = None
122 122 setup = []
123 123 if mqtt is 1:
124 124 self.client = mqtt.Client(
125 125 client_id=self.clientId + self.topic + 'SCHAIN',
126 126 clean_session=True)
127 127 self.client.on_disconnect = self.on_disconnect
128 128 self.connect()
129 129 for plot in self.plottype:
130 130 setup.append({
131 131 'plot': plot,
132 132 'topic': self.topic + plot,
133 133 'title': getattr(self, plot + '_' + 'title', False),
134 134 'xlabel': getattr(self, plot + '_' + 'xlabel', False),
135 135 'ylabel': getattr(self, plot + '_' + 'ylabel', False),
136 136 'xrange': getattr(self, plot + '_' + 'xrange', False),
137 137 'yrange': getattr(self, plot + '_' + 'yrange', False),
138 138 'zrange': getattr(self, plot + '_' + 'zrange', False),
139 139 })
140 140 if zeromq is 1:
141 141 context = zmq.Context()
142 142 self.zmq_socket = context.socket(zmq.PUSH)
143 143 server = kwargs.get('server', 'zmq.pipe')
144 144
145 145 if 'tcp://' in server:
146 146 address = server
147 147 else:
148 148 address = 'ipc:///tmp/%s' % server
149 149
150 150 self.zmq_socket.connect(address)
151 151 time.sleep(1)
152 152
153 153 def publish_data(self):
154 154 self.dataOut.finished = False
155 155 if self.mqtt is 1:
156 156 yData = self.dataOut.heightList[:2].tolist()
157 157 if self.plottype == 'spectra':
158 158 data = getattr(self.dataOut, 'data_spc')
159 159 z = data/self.dataOut.normFactor
160 160 zdB = 10*numpy.log10(z)
161 161 xlen, ylen = zdB[0].shape
162 162 dx = int(xlen/MAXNUMX) + 1
163 163 dy = int(ylen/MAXNUMY) + 1
164 164 Z = [0 for i in self.dataOut.channelList]
165 165 for i in self.dataOut.channelList:
166 166 Z[i] = zdB[i][::dx, ::dy].tolist()
167 167 payload = {
168 168 'timestamp': self.dataOut.utctime,
169 169 'data': roundFloats(Z),
170 170 'channels': ['Ch %s' % ch for ch in self.dataOut.channelList],
171 171 'interval': self.dataOut.getTimeInterval(),
172 172 'type': self.plottype,
173 173 'yData': yData
174 174 }
175 175 # print payload
176 176
177 177 elif self.plottype in ('rti', 'power'):
178 178 data = getattr(self.dataOut, 'data_spc')
179 179 z = data/self.dataOut.normFactor
180 180 avg = numpy.average(z, axis=1)
181 181 avgdB = 10*numpy.log10(avg)
182 182 xlen, ylen = z[0].shape
183 183 dy = numpy.floor(ylen/self.__MAXNUMY) + 1
184 184 AVG = [0 for i in self.dataOut.channelList]
185 185 for i in self.dataOut.channelList:
186 186 AVG[i] = avgdB[i][::dy].tolist()
187 187 payload = {
188 188 'timestamp': self.dataOut.utctime,
189 189 'data': roundFloats(AVG),
190 190 'channels': ['Ch %s' % ch for ch in self.dataOut.channelList],
191 191 'interval': self.dataOut.getTimeInterval(),
192 192 'type': self.plottype,
193 193 'yData': yData
194 194 }
195 195 elif self.plottype == 'noise':
196 196 noise = self.dataOut.getNoise()/self.dataOut.normFactor
197 197 noisedB = 10*numpy.log10(noise)
198 198 payload = {
199 199 'timestamp': self.dataOut.utctime,
200 200 'data': roundFloats(noisedB.reshape(-1, 1).tolist()),
201 201 'channels': ['Ch %s' % ch for ch in self.dataOut.channelList],
202 202 'interval': self.dataOut.getTimeInterval(),
203 203 'type': self.plottype,
204 204 'yData': yData
205 205 }
206 206 elif self.plottype == 'snr':
207 207 data = getattr(self.dataOut, 'data_SNR')
208 208 avgdB = 10*numpy.log10(data)
209 209
210 210 ylen = data[0].size
211 211 dy = numpy.floor(ylen/self.__MAXNUMY) + 1
212 212 AVG = [0 for i in self.dataOut.channelList]
213 213 for i in self.dataOut.channelList:
214 214 AVG[i] = avgdB[i][::dy].tolist()
215 215 payload = {
216 216 'timestamp': self.dataOut.utctime,
217 217 'data': roundFloats(AVG),
218 218 'channels': ['Ch %s' % ch for ch in self.dataOut.channelList],
219 219 'type': self.plottype,
220 220 'yData': yData
221 221 }
222 222 else:
223 223 print "Tipo de grafico invalido"
224 224 payload = {
225 225 'data': 'None',
226 226 'timestamp': 'None',
227 227 'type': None
228 228 }
229 229 # print 'Publishing data to {}'.format(self.host)
230 230 self.client.publish(self.topic + self.plottype, json.dumps(payload), qos=0)
231 231
232 232 if self.zeromq is 1:
233 233 print '[Sending] {} - {}'.format(self.dataOut.type, self.dataOut.datatime)
234 234 self.zmq_socket.send_pyobj(self.dataOut)
235 235
236 236 def run(self, dataOut, **kwargs):
237 237 self.dataOut = dataOut
238 238 if not self.isConfig:
239 239 self.setup(**kwargs)
240 240 self.isConfig = True
241 241
242 242 self.publish_data()
243 243 time.sleep(self.delay)
244 244
245 245 def close(self):
246 246 if self.zeromq is 1:
247 247 self.dataOut.finished = True
248 248 self.zmq_socket.send_pyobj(self.dataOut)
249 249
250 250 if self.client:
251 251 self.client.loop_stop()
252 252 self.client.disconnect()
253 253
254 254
255 255 class ReceiverData(ProcessingUnit, Process):
256 256
257 257 throttle_value = 5
258 258
259 259 def __init__(self, **kwargs):
260 260
261 261 ProcessingUnit.__init__(self, **kwargs)
262 262 Process.__init__(self)
263 263 self.mp = False
264 264 self.isConfig = False
265 265 self.isWebConfig = False
266 266 self.plottypes =[]
267 267 self.connections = 0
268 268 server = kwargs.get('server', 'zmq.pipe')
269 269 plot_server = kwargs.get('plot_server', 'zmq.web')
270 270 if 'tcp://' in server:
271 271 address = server
272 272 else:
273 273 address = 'ipc:///tmp/%s' % server
274 274
275 275 if 'tcp://' in plot_server:
276 276 plot_address = plot_server
277 277 else:
278 278 plot_address = 'ipc:///tmp/%s' % plot_server
279 279
280 280 self.address = address
281 281 self.plot_address = plot_address
282 282 self.plottypes = [s.strip() for s in kwargs.get('plottypes', 'rti').split(',')]
283 283 self.realtime = kwargs.get('realtime', False)
284 284 self.throttle_value = kwargs.get('throttle', 5)
285 285 self.sendData = self.initThrottle(self.throttle_value)
286 286 self.setup()
287 287
288 288 def setup(self):
289 289
290 290 self.data = {}
291 291 self.data['times'] = []
292 292 for plottype in self.plottypes:
293 293 self.data[plottype] = {}
294 294 self.data['noise'] = {}
295 295 self.data['throttle'] = self.throttle_value
296 296 self.data['ENDED'] = False
297 297 self.isConfig = True
298 298 self.data_web = {}
299 299
300 300 def event_monitor(self, monitor):
301 301
302 302 events = {}
303 303
304 304 for name in dir(zmq):
305 305 if name.startswith('EVENT_'):
306 306 value = getattr(zmq, name)
307 307 events[value] = name
308 308
309 309 while monitor.poll():
310 310 evt = recv_monitor_message(monitor)
311 311 if evt['event'] == 32:
312 312 self.connections += 1
313 313 if evt['event'] == 512:
314 314 pass
315 315 if self.connections == 0 and self.started is True:
316 316 self.ended = True
317 317
318 318 evt.update({'description': events[evt['event']]})
319 319
320 320 if evt['event'] == zmq.EVENT_MONITOR_STOPPED:
321 321 break
322 322 monitor.close()
323 323 print("event monitor thread done!")
324 324
325 325 def initThrottle(self, throttle_value):
326 326
327 327 @throttle(seconds=throttle_value)
328 328 def sendDataThrottled(fn_sender, data):
329 329 fn_sender(data)
330 330
331 331 return sendDataThrottled
332 332
333 333 def send(self, data):
334 334 # print '[sending] data=%s size=%s' % (data.keys(), len(data['times']))
335 335 self.sender.send_pyobj(data)
336 336
337 337 def update(self):
338 338
339 339 t = self.dataOut.utctime
340 340
341 341 if t in self.data['times']:
342 342 return
343 343
344 344 self.data['times'].append(t)
345 345 self.data['dataOut'] = self.dataOut
346 346
347 347 for plottype in self.plottypes:
348 348 if plottype == 'spc':
349 349 z = self.dataOut.data_spc/self.dataOut.normFactor
350 350 self.data[plottype] = 10*numpy.log10(z)
351 351 self.data['noise'][t] = 10*numpy.log10(self.dataOut.getNoise()/self.dataOut.normFactor)
352 352 if plottype == 'cspc':
353 353 jcoherence = self.dataOut.data_cspc/numpy.sqrt(self.dataOut.data_spc*self.dataOut.data_spc)
354 354 self.data['cspc_coh'] = numpy.abs(jcoherence)
355 355 self.data['cspc_phase'] = numpy.arctan2(jcoherence.imag, jcoherence.real)*180/numpy.pi
356 356 if plottype == 'rti':
357 357 self.data[plottype][t] = self.dataOut.getPower()
358 358 if plottype == 'snr':
359 359 self.data[plottype][t] = 10*numpy.log10(self.dataOut.data_SNR)
360 360 if plottype == 'dop':
361 361 self.data[plottype][t] = 10*numpy.log10(self.dataOut.data_DOP)
362 362 if plottype == 'mean':
363 363 self.data[plottype][t] = self.dataOut.data_MEAN
364 364 if plottype == 'std':
365 365 self.data[plottype][t] = self.dataOut.data_STD
366 366 if plottype == 'coh':
367 367 self.data[plottype][t] = self.dataOut.getCoherence()
368 368 if plottype == 'phase':
369 369 self.data[plottype][t] = self.dataOut.getCoherence(phase=True)
370 if plottype == 'wind':
370 if plottype == 'output':
371 371 self.data[plottype][t] = self.dataOut.data_output
372 if plottype == 'param':
373 self.data[plottype][t] = self.dataOut.data_param
372 374 if self.realtime:
373 375 self.data_web['timestamp'] = t
374 376 if plottype == 'spc':
375 377 self.data_web[plottype] = roundFloats(decimate(self.data[plottype]).tolist())
376 378 elif plottype == 'cspc':
377 379 self.data_web['cspc_coh'] = roundFloats(decimate(self.data['cspc_coh']).tolist())
378 380 self.data_web['cspc_phase'] = roundFloats(decimate(self.data['cspc_phase']).tolist())
379 381 elif plottype == 'noise':
380 382 self.data_web['noise'] = roundFloats(self.data['noise'][t].tolist())
381 383 else:
382 384 self.data_web[plottype] = roundFloats(decimate(self.data[plottype][t]).tolist())
383 385 self.data_web['interval'] = self.dataOut.getTimeInterval()
384 386 self.data_web['type'] = plottype
385 387
386 388 def run(self):
387 389
388 390 print '[Starting] {} from {}'.format(self.name, self.address)
389 391
390 392 self.context = zmq.Context()
391 393 self.receiver = self.context.socket(zmq.PULL)
392 394 self.receiver.bind(self.address)
393 395 monitor = self.receiver.get_monitor_socket()
394 396 self.sender = self.context.socket(zmq.PUB)
395 397 if self.realtime:
396 398 self.sender_web = self.context.socket(zmq.PUB)
397 399 self.sender_web.connect(self.plot_address)
398 400 time.sleep(1)
399 self.sender.bind("ipc:///tmp/zmq.plots")
401 if 'server' in self.kwargs:
402 self.sender.bind("ipc:///tmp/{}.plots".format(self.kwargs['server']))
403 else:
404 self.sender.bind("ipc:///tmp/zmq.plots")
400 405
401 406 t = Thread(target=self.event_monitor, args=(monitor,))
402 407 t.start()
403 408
404 409 while True:
405 410 self.dataOut = self.receiver.recv_pyobj()
406 411 # print '[Receiving] {} - {}'.format(self.dataOut.type,
407 412 # self.dataOut.datatime.ctime())
408 413
409 414 self.update()
410 415
411 416 if self.dataOut.finished is True:
412 417 self.send(self.data)
413 418 self.connections -= 1
414 419 if self.connections == 0 and self.started:
415 420 self.ended = True
416 421 self.data['ENDED'] = True
417 422 self.send(self.data)
418 423 self.setup()
419 424 else:
420 425 if self.realtime:
421 426 self.send(self.data)
422 427 self.sender_web.send_string(json.dumps(self.data_web))
423 428 else:
424 429 self.sendData(self.send, self.data)
425 430 self.started = True
426 431
427 432 return
428 433
429 434 def sendToWeb(self):
430 435
431 436 if not self.isWebConfig:
432 437 context = zmq.Context()
433 438 sender_web_config = context.socket(zmq.PUB)
434 439 if 'tcp://' in self.plot_address:
435 440 dum, address, port = self.plot_address.split(':')
436 441 conf_address = '{}:{}:{}'.format(dum, address, int(port)+1)
437 442 else:
438 443 conf_address = self.plot_address + '.config'
439 444 sender_web_config.bind(conf_address)
440 445 time.sleep(1)
441 446 for kwargs in self.operationKwargs.values():
442 447 if 'plot' in kwargs:
443 448 print '[Sending] Config data to web for {}'.format(kwargs['code'].upper())
444 449 sender_web_config.send_string(json.dumps(kwargs))
445 450 self.isWebConfig = True
General Comments 0
You need to be logged in to leave comments. Login now